Feature Selection & Dimensionality Reduction

Assumes you know from ML-101

This lesson builds on ML-101 Lesson 2 (Data & Features), ML-101 Lesson 10 (Overfitting, Bias & Variance), and ML-101 Lesson 13 (Unsupervised Learning). The reader is assumed to know that features are columns of a design matrix $X \in \mathbb{R}^{n \times p}$, that the variance of a column is the average squared deviation from its mean, and that the covariance between two columns is the average product of their deviations. The reader is also expected to know what unsupervised learning means — that the algorithm receives $X$ with no labels $y$ — and that k-means clusters points by Euclidean distance to a centroid.

What ML-101 did not do is derive which linear combination of features captures the most variance, why an L1 penalty produces exact zeros rather than small coefficients, or how to attribute a single prediction to its inputs in a way that is mathematically fair. This lesson derives each from first principles, names the assumption each one rests on, and shows what fails when the assumption is violated.

Learning Objectives

1. Derive PCA as the orthogonal projection that maximises retained variance, and show the algebraic identity between maximum variance and minimum reconstruction error.
2. Compute the explained-variance ratio and decide how many components to keep, naming the rule whose violation leads to either over-compression or noise retention.
3. Derive the L1-regularised least-squares solution (Lasso) and prove that the L1 constraint set has corners, which is why it yields exact zeros while L2 does not.
4. Implement permutation importance and state the assumption whose violation makes the reported importances a story rather than a measurement.
5. Define the Shapley value as the unique attribution that satisfies efficiency, symmetry, null player and additivity, and show why a 201-level practitioner treats it as the default interpretation tool.

PCA as a variance-maximising rotation

Take a centred design matrix $X \in \mathbb{R}^{n \times p}$ with columns of mean zero. The sample covariance matrix is

$\Sigma$ is symmetric positive-semidefinite, so it admits an eigendecomposition $\Sigma = W \Lambda W^{\top}$ where $W$ is orthogonal ($W^{\top} W = I$) and $\Lambda = \mathrm{diag}(\lambda_1, \dots, \lambda_p)$ with $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p \ge 0$. The columns of $W$ are the principal axes.

For a unit direction $v \in \mathbb{R}^p$, the variance of the data projected onto $v$ is

Maximise $v^{\top}\Sigma v$ subject to $v^{\top}v = 1$. The Lagrangian is $\mathcal{L}(v, \mu) = v^{\top}\Sigma v - \mu(v^{\top}v - 1)$. The first-order condition is $2\Sigma v - 2\mu v = 0$, i.e.

The optimum $v$ must be an eigenvector of $\Sigma$, and the optimum value $\mu$ is the corresponding eigenvalue. Picking the largest eigenvalue gives the direction of greatest variance.

This is the Rayleigh quotient argument, and it is the same one that underlies Fisher's linear discriminant, the power iteration, and the spectral clustering normalised-cut relaxation. The constraint $v^{\top}v=1$ is what makes the problem well-posed: scale is not informative for a direction.

import numpy as np

def pca_first_component(X: np.ndarray) -> np.ndarray:
    """
    Return the first principal axis of a centred design matrix.

    The power iteration converges geometrically; here we use the closed form
    because `numpy.linalg.eigh` is numerically stable for symmetric matrices.
    """
    Xc = X - X.mean(axis=0, keepdims=True)
    cov = (Xc.T @ Xc) / Xc.shape[0]
    eigvals, eigvecs = np.linalg.eigh(cov)
    # eigh returns ascending order; the largest is the last
    return eigvecs[:, -1]


# 3-D data stretched along one direction: PC1 should align with that axis.
rng = np.random.default_rng(0)
X = rng.normal(size=(2000, 3))
X[:, 0] *= 5.0
w1 = pca_first_component(X)
print("PC1 (should be ~[1, 0, 0]):", w1.round(3))

The eigendecomposition is the cleanest path. For streaming or huge $p$, the same fixed point can be reached by power iteration: multiply by $\Sigma$ repeatedly and re-normalise, since $\Sigma^k v$ converges to the dominant eigenvector for any non-orthogonal starting $v$. The convergence rate is $|\lambda_2/\lambda_1|^k$, which is geometric but slow when eigenvalues are close.

Explained variance and the reconstruction view

Projecting $X$ onto the top-$k$ principal axes $W_k \in \mathbb{R}^{p \times k}$ produces scores $Z = X W_k \in \mathbb{R}^{n \times k}$. The reconstruction is $\hat{X} = Z W_k^{\top}$.

The reconstruction error is the squared Frobenius norm of the residual:

The explained variance ratio of the top-$k$ components is

PCA is both the projection that maximises variance and the projection that minimises reconstruction error. The two are not separate statements; they are algebraically the same:

and these two quantities sum to $\mathrm{tr}(\Sigma) = \sum_{i=1}^{p}\lambda_i$, the total variance. Maximising one is exactly minimising the other.

The decision of how many components to keep is a model-selection problem with a clear elbow rule. Plot $\mathrm{EVR}_k$ against $k$ and look for the point where adding components stops paying for itself. A common quantitative choice is to keep enough components to explain 90% or 95% of the variance. The threshold is arbitrary; the choice is not — picking too few loses signal, picking too many keeps noise.

The Kaiser rule deserves a sharp warning. If features are on different scales, the covariance is dominated by the largest-scale feature and the eigenvalues inherit that scale. The fix is to standardise each column to unit variance before computing PCA — but doing so means the components are no longer directions in the original feature space, and "PC1" loses its interpretability as a linear combination of the raw measurements.

L1 selection: why the geometry gives zeros

PCA keeps combinations of features. Lasso drops whole features. The Lasso estimator solves

The penalty is the L1 norm $\lVert \beta \rVert_1 = \sum_j \lvert \beta_j \rvert$. Compare with Ridge, which uses $\lVert \beta \rVert_2^2$.

The reason L1 produces zeros and L2 does not is geometric. Both penalties carve a constraint set out of $\mathbb{R}^p$ — the L1 ball is a diamond with corners on the axes; the L2 ball is a sphere. The Lasso solution is the first point where the diamond touches a level set of the squared-error loss. Because the diamond has corners where one coordinate is zero, the solution is much more likely to land at a corner than the sphere is, and landing at a corner means a zero coefficient.

Mathematically, the subgradient optimality condition at the solution is

where $s \in \partial \lVert \hat{\beta} \rVert_1$ is a subgradient with $s_j = \mathrm{sign}(\hat{\beta}_j)$ when $\hat{\beta}_j \ne 0$ and $s_j \in [-1, 1]$ when $\hat{\beta}_j = 0$. If the data-driven term is small enough in magnitude that even its largest plausible $s_j$ cannot balance it, the only solution is $\hat{\beta}_j = 0$.

import numpy as np

def soft_threshold(z: np.ndarray, alpha: float) -> np.ndarray:
    """
    The proximal operator of alpha * ||.||_1.

    Apply coordinate-wise soft thresholding: shrink |z| by alpha and set to
    zero if the result is non-positive. The Lasso coordinate update reduces
    to this operation in the orthogonal-design case.
    """
    return np.sign(z) * np.maximum(np.abs(z) - alpha, 0.0)


# Demonstrate sparsity: 1000 features, only 5 have signal.
rng = np.random.default_rng(1)
n, p = 200, 1000
X = rng.normal(size=(n, p))
beta_true = np.zeros(p)
beta_true[:5] = [3.0, -2.0, 1.5, 0.0, 4.0]
y = X @ beta_true + 0.1 * rng.normal(size=n)

# Coordinate-descent Lasso (illustrative; uses sklearn under the hood)
from sklearn.linear_model import Lasso
model = Lasso(alpha=0.1, max_iter=20000).fit(X, y)
nonzero = np.flatnonzero(model.coef_)
print(f"selected {len(nonzero)} features; expected ~5")

The output is a handful of non-zero coefficients. The Lasso has identified the support of $\beta_{\text{true}}$ in the limit of large $n$ and small $\alpha$, but the guarantee is conditional on the design matrix satisfying a restricted eigenvalue condition — a quantitative statement that the columns of $X$ are not too correlated. When features are highly correlated, Lasso's selection is unstable: it picks one and ignores the rest, arbitrarily.

The intuition worth holding onto: Lasso's sparsity is a geometric accident. L2 is the natural penalty under a Gaussian prior (the maximum-a-posteriori estimator is Ridge), but L2 does not yield zeros. The L1 norm is not the natural penalty for any Gaussian likelihood — it corresponds to a Laplace prior — and the corners of the L1 ball are what deliver sparsity as a side effect of a non-differentiable geometry, not as a deeper statistical truth.

Permutation importance

Model-agnostic feature importance starts from a fitted model $f$ and a held-out evaluation set $(X_{\text{test}}, y_{\text{test}})$. The baseline score is $s_0 = \mathrm{score}(f, X_{\text{test}}, y_{\text{test}})$. To assess feature $j$:

1. Shuffle column $j$ of $X_{\text{test}}$ to produce $X_{\text{test}}^{(j)}$.
2. Re-score: $s_j = \mathrm{score}(f, X_{\text{test}}^{(j)}, y_{\text{test}})$.
3. Report $\Delta_j = s_0 - s_j$.

Large $\Delta_j$ means the model depended on column $j$ for its score; small $\Delta_j$ means it did not. The procedure is repeated $K$ times and the mean and standard deviation are reported, giving a confidence interval without parametric assumptions.

import numpy as np
from sklearn.inspection import permutation_importance

def permutation_importance_report(model, X_test, y_test, n_repeats=30):
    """
    Compute permutation importances with confidence intervals.

    The `n_repeats` argument controls the Monte Carlo error on the estimate;
    30 is enough to see whether the lower bound of the CI excludes zero.
    """
    r = permutation_importance(model, X_test, y_test,
                               n_repeats=n_repeats, random_state=0)
    out = sorted(zip(r.importances_mean, r.importances_std,
                     X_test.columns), reverse=True)
    for mean, std, name in out:
        lo, hi = mean - 2 * std, mean + 2 * std
        print(f"{name:>16s}  {mean:+.4f}  CI [{lo:+.4f}, {hi:+.4f}]")

The assumption is feature independence: shuffling column $j$ while leaving the others intact breaks only the dependence between $X_j$ and $y$, not any dependence between $X_j$ and $X_{j'}$ that the model might have exploited. When features are correlated, the test set has $X_j$ and $X_{j'}$ in their natural co-occurrence pattern. Shuffling $X_j$ produces a row where $X_j$ is no longer consistent with $X_{j'}$, which the model has never seen, and the resulting drop in score is partly a measurement of the model's lack of robustness rather than $X_j$'s true contribution.

The 201-level judgement is that permutation importance on correlated inputs does not average out — it is structurally biased toward under-counting important-but-correlated features and over-counting unique-but-noisy ones. The fix is conditional permutation: build a model $\hat{X}_j = g(X_{\setminus j})$, then replace $X_j$ with samples from $g(\cdot \mid X_{\setminus j})$ instead of unconditional shuffling. The feature's marginal contribution is now measured against the joint distribution the model actually sees.

Shapley values and SHAP

Shapley values come from cooperative game theory and answer the question: in a coalition of $p$ players that produces a total payoff $v(S)$, how much of that payoff is attributable to each individual player? The answer, uniquely defined by four axioms (efficiency, symmetry, null player, additivity), is

For an ML model with features $X_1, \dots, X_p$ and a prediction $f(x)$, the value function $v(S)$ is the expected prediction when only the features in $S$ are known and the others are integrated out over a background distribution. Shapley values are then the unique attribution that satisfies:

• Efficiency: $\sum_j \phi_j = f(x) - \mathbb{E}[f(X)]$, so the attributions sum to the deviation from the average prediction.
• Symmetry: two features that contribute identically across all coalitions receive identical attributions.
• Null player: a feature that never changes the prediction gets $\phi_j = 0$.
• Additivity: for an ensemble of models $f = f_a + f_b$, the Shapley values add.

The additive axiom is why SHAP is the default interpretation tool at the 201 level: when you ensemble models, you do not have to re-derive the attributions; they compose. This is not a property any other commonly-used importance method has.

The exact Shapley computation is exponential in $p$, so practical implementations approximate. Kernel SHAP estimates a linear surrogate around $f$ that is consistent with the Shapley axioms. Tree SHAP computes exact Shapley values for tree ensembles in polynomial time by exploiting the additive structure of the splits. Deep SHAP propagates Shapley values through a neural network by linearising each layer. Each approximation trades off accuracy, compute, and the assumptions it can afford to make.

import shap

def explain_with_shap(model, X_background, X_explain):
    """
    Compute SHAP values for a fitted model on a small set of examples.

    `X_background` is a sample from the training distribution used to
    estimate E[f(X)]; `X_explain` is the set of points to attribute.
    """
    explainer = shap.KernelExplainer(model.predict, X_background)
    shap_values = explainer.shap_values(X_explain, nsamples=200)
    # shap_values has shape (n_explain, n_features)
    return shap_values

The output is a matrix whose row $i$ sums (within numerical error) to $f(x_i) - \mathbb{E}[f(X)]$. The columns tell you how each feature pushed that specific prediction above or below the average — a local, signed, model-faithful attribution. It is not a global importance ranking; the global importance of feature $j$ is the mean absolute Shapley value $\frac{1}{n}\sum_i \lvert \phi_j(x_i) \rvert$, which is a defensible summary but loses the per-instance sign.

The assumption whose violation breaks SHAP is feature independence in the background distribution. The value function $v(S) = \mathbb{E}[f(X) \mid X_S = x_S]$ integrates over $X_{\setminus S}$ using whatever distribution you provided as a reference. If the reference is uncorrelated draws from marginals and the true features are correlated, $v(S)$ is no longer the model's behaviour on realistic inputs and the attributions drift. The fix is to use a background sample drawn from the joint training distribution and to be explicit about which subset you chose, because every choice is a counterfactual.

Choosing between methods

The four tools answer different questions and are not interchangeable.

The 201-level failure modes:

• Using PCA on raw unscaled features and reading the components as "the directions the data lives along". The directions are dominated by whichever feature has the largest variance.
• Using Lasso on highly correlated features and trusting the support. The support is one of many possible equivalent ones.
• Reporting permutation importance on correlated features as if it were a clean measurement. It is not; the procedure produces a number but the number is not what its name suggests.
• Reporting mean absolute SHAP as if it were the only correct global ranking. It is the correct summary under the Shapley axioms, but the axioms depend on a chosen background distribution, and changing the background changes the ranking.

The unifying frame: each method is optimal under a stated assumption (orthogonal Gaussian data, restricted eigenvalue condition, feature independence, well-specified background distribution), and each method degrades gracefully up to the point where the assumption fails and then degrades catastrophically. The practitioner's job is to know which assumption is load-bearing for the decision they are about to make.

Mutual information and feature ranking

PCA and Lasso are both linear. The mutual information between a feature and the target is a model-free measure of dependence:

$I = 0$ iff $X_j$ and $Y$ are independent under the joint distribution. Unlike correlation, MI captures non-monotone dependence: a feature that is informative through a step function or a sign flip has $I > 0$ even though its Pearson correlation is zero. The estimator in scikit-learn (mutual_info_classif) bins continuous features and applies a correction for the bias that small bins have toward large $I$.

import numpy as np
from sklearn.feature_selection import mutual_info_classif

def mi_ranking(X: np.ndarray, y: np.ndarray) -> np.ndarray:
    """
    Rank features by estimated mutual information with the label.

    `mutual_info_classif` bins continuous features internally; the
    `random_state` argument controls the bin assignments for reproducibility.
    """
    mi = mutual_info_classif(X, y, random_state=0, n_neighbors=3)
    return np.argsort(-mi)

The failure mode of MI on a small sample is variance. The plug-in estimator $\hat{I} = \sum_{x,y} \hat{p}(x,y)\log[\hat{p}(x,y)/(\hat{p}(x)\hat{p}(y))]$ is biased upward for sparse bins; the corrected estimator subtracts a per-bin term and is consistent as $n \to \infty$, but for $n < 1000$ the ranking is unstable across random subsamples. The 201-level practice is to bootstrap the ranking: sample rows with replacement, recompute MI, and report the features whose rank stays in the top-$K$ across resamples. Anything that survives bootstrap MI is robustly informative.

Recursive feature elimination and stability

RFE starts from a fitted model $m$, computes an importance ranking (the absolute value of the coefficient for linear models, the impurity decrease for trees), drops the bottom feature, refits, and repeats. The number of features to keep is itself a hyperparameter and is set by nested cross-validation.

The deeper issue is stability: if the ranking of the bottom-10 features changes when the training set is resampled, RFE's choice is noise. The stability score

averages the Jaccard similarity between the chosen support $S$ and the support $S_m$ found on the $m$-th bootstrap. A score near 1 means the selection is data-invariant; a score near 0 means RFE is selecting whichever feature happened to be slightly less noisy this run.

The 201-level rule: do not trust a single RFE run on a single train/test split. Run it on $M = 30$ bootstraps, report the support that appears in at least 80% of them, and acknowledge that the "selected features" are a consensus, not a fact. The opposite mistake — selecting a feature because RFE picked it once — is silent overfitting at the meta level: the selection procedure is being fit to a specific noise realisation in the training set.

Tree impurity importance and its known flaw

Random forests expose feature_importances_, computed as the total decrease in node impurity (Gini or entropy) weighted by the probability of reaching the node. The number is fast and global, but it has a documented structural bias: it over-weights high-cardinality features and features with many possible split points. A feature that can take many distinct values gets more split opportunities, and each split that improves the impurity by a little contributes to the total — even when the feature is less informative than a low-cardinality competitor.

The fix is permutation importance on top of the impurity ranking: do the impurity ranking for an initial screen, then use permutation importance to validate that the top features remain important after the corruption. A feature that is high-impurity but low-permutation is one the tree splits on but the test score does not depend on — it is fitting noise. The two rankings together are what the practitioner reports.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance

def impurity_and_perm(X_tr, y_tr, X_te, y_te) -> tuple[np.ndarray, np.ndarray]:
    rf = RandomForestClassifier(n_estimators=300, random_state=0,
                                n_jobs=-1).fit(X_tr, y_tr)
    imp = rf.feature_importances_
    perm = permutation_importance(rf, X_te, y_te, n_repeats=20,
                                  random_state=0, n_jobs=-1)
    return imp, perm.importances_mean

The function returns two arrays. The 201-level judgement is the intersection: features that are high in both arrays are robustly important; features high in one but not the other need a hypothesis about why (cardinality bias for high-impurity-only; leakage or surrogate split for high-permutation-only). Reporting only impurity importances is a tell that the practitioner is using a default ranking they have not audited.

Worked example: 200 features, 5 informative, choosing the support

The following walks the full pipeline: build a sparse signal, run Lasso, run MI ranking, run RFE, and check whether the supports agree.

import numpy as np
from sklearn.linear_model import Lasso
from sklearn.feature_selection import mutual_info_classif, RFE
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(0)
n, p = 300, 200
X = rng.normal(size=(n, p))
beta = np.zeros(p)
beta[:5] = [3.0, -2.5, 1.8, -1.2, 2.2]
y = (X @ beta + 0.5 * rng.normal(size=n) > 0).astype(int)

# Lasso support
lasso = Lasso(alpha=0.05, max_iter=20000).fit(X, y)
lasso_support = np.flatnonzero(lasso.coef_)

# MI ranking, top 5
mi = mutual_info_classif(X, y, random_state=0)
mi_support = np.argsort(-mi)[:5]

# RFE with logistic regression
rfe = RFE(LogisticRegression(max_iter=2000),
          n_features_to_select=5, step=0.1).fit(X, y)
rfe_support = np.flatnonzero(rfe.support_)

print(f"true support  : {np.flatnonzero(beta).tolist()}")
print(f"lasso support : {lasso_support.tolist()}")
print(f"mi top-5      : {mi_support.tolist()}")
print(f"rfe support   : {rfe_support.tolist()}")

The three methods converge on the same support when $n$ is large and features are independent. The point of the exercise is what happens when features are correlated: Lasso picks one and drops the rest, MI splits its weight across the correlated cluster, and RFE is whichever the underlying model happened to find first. The 201-level practitioner reports all three supports and labels the features that survive every method as "robust". The features that survive one but not the others are flagged for manual review.

The unifying judgement across the lesson: feature selection and dimensionality reduction are not the same problem. PCA compresses; Lasso selects; MI ranks; RFE bootstraps; Shapley attributes. The choice depends on whether the downstream model wants compact inputs, a sparse representation, a defensible ranking, or a faithful per-prediction attribution. Conflating them is the 101-level mistake that produces a "selected features" list that no one can defend at a model review.

Key Takeaways

• PCA is the unique orthogonal rotation whose top-$k$ subspace simultaneously maximises retained variance and minimises reconstruction error; the two are the same Rayleigh quotient.
• Lasso's sparsity comes from the corners of the L1 ball, not from any likelihood: L2 (Ridge) corresponds to a Gaussian prior and gives small coefficients, never zeros.
• Permutation importance is a fast screen but assumes feature independence; on correlated inputs it under-counts the features that mattered.
• Shapley values are the unique attribution under four game-theoretic axioms; SHAP is the practical estimator and composes under ensembling, which is what makes it the default at the 201 level.
• Choose the method by the question you are asking, not by the tool's reputation: PCA compresses, Lasso selects, permutation ranks, Shapley attributes. The 201-level mistake is to swap them.
• Mutual information captures non-monotone dependence but is high-variance on small samples; tree impurity importance over-counts high-cardinality features; recursive feature elimination's stability must itself be bootstrapped to be trusted.