38. Machine Learning Principles and Algorithms

Machine Learning Principles and Algorithms¶

Overview¶

Most predictive brain models draw on machine learning, a discipline at the intersection of computer science and statistics. Machine learning is a large umbrella of techniques. Supervised learning algorithms optimize prediction of known outcomes from labeled data — regression when the outcome is continuous, classification when it is categorical. Unsupervised algorithms instead group unlabeled data into clusters or coherent components (clustering, PCA, ICA; see Chapter 31). In neuroimaging, the observations are often participants or trials, the features ($x$) are voxels, connections, or graph metrics, and the outcome ($y$) might be patient status, stimulus category, or a continuous symptom score.

One useful way to understand machine learning is as a set of variations on classical statistical models. The goal is to learn a function $f(x)$ that accurately predicts $y$; multiple linear regression is the special case $f(x) = w^T x + b$, where $w$ is a vector of feature weights and $b$ a scalar intercept (the bias), with weights estimated by minimizing a loss function. Machine learning is flexible about how that loss is defined. A general and enormously influential form adds a regularization (penalty) term to the measure of fit:

Other algorithms swap the fit term itself — support vector machines use a hinge loss, and absolute-error losses resist influential data points better than squared error.

The deepest difference from classical statistics is the explicit focus on generalizability. In classical regression, error is evaluated on the same data used for fitting, so adding parameters always improves apparent fit — and if the number of features $p$ reaches the number of observations $n$ ($p \geq n$), a linear model can fit any dataset perfectly, with zero error, without capturing anything that generalizes. This is overfitting: the model absorbs noise idiosyncratic to the sample. The machine learning remedy is to partition data into independent training and test sets. Even a wildly overfit model can then be fairly evaluated on held-out data, and models of different complexity — including non-nested models with different structures — can be compared on equal footing. Cross-validation (Chapter 39) systematizes this within a single dataset, and is typically followed by validation on fully independent data.

Complexity itself trades off two sources of error. Simple models underfit: they are systematically wrong in the same way in every sample (high bias). Complex models fit the training data beautifully but change drastically from sample to sample (high variance). Expected prediction error on new data decomposes, approximately, into $\text{bias}^2 + \text{variance} + \text{irreducible noise}$, so test error follows a characteristic U-shaped curve as complexity grows — falling while added flexibility captures real structure, then rising as the model begins to memorize noise. Training error, by contrast, only decreases. Regularization, ensemble averaging (bagging, boosting, random forests), and dimension reduction are all ways of managing this tradeoff — accepting a little bias to buy a large reduction in variance.

Classification algorithms use an object’s features to decide which class it belongs to. A binary linear classifier computes a score $z = w^T x + b$ and predicts class +1 when $z > 0$ and -1 otherwise. Geometrically, $w$ and the bias $b$ define a hyperplane that splits feature space in two; algorithms differ in how they place it. Linear discriminant analysis maximizes between-class relative to within-class variance under normality assumptions; logistic regression models the log odds of class membership, yielding calibrated probabilities. The support vector machine (SVM) chooses the hyperplane with the largest margin — the greatest distance to the nearest training points, which are called the support vectors because they alone determine the solution. When classes overlap, slack variables $\xi_i$ measure each point’s degree of misclassification, and the objective becomes

SVMs are accurate, fast, and stable with very large numbers of features — one reason for their popularity in fMRI, where kernels map the data into an $n \times n$ similarity space so computation scales with observations rather than voxels. Simple nonlinear alternatives include $k$-nearest neighbors (majority vote among the $k$ most similar training examples) and naive Bayes (combining per-feature class probabilities under an independence assumption — often unrealistic, yet remarkably effective when features are many and examples few).

Regression algorithms predict continuous outcomes, and here regularization does the heavy lifting in high dimensions. The two most common penalties are L1 and L2:

L2 regularization defines ridge regression: all coefficients shrink smoothly toward zero, stabilizing estimates when predictors are highly multicollinear — but none become exactly zero, so every feature stays in the model. L1 regularization defines the LASSO: as $\lambda$ grows, coefficients of the least important features hit exactly zero, producing a sparse model. Sparsity is attractive for simple measurement models, but applied directly to brain voxels it tends to pick one arbitrary representative from each set of correlated voxels, scattering “speckles” of predictive features across the brain. Structured variants (group LASSO, fused LASSO, GraphNet) constrain selection to spatially coherent sets; the elastic net combines both penalties and behaves better than LASSO with correlated features; and LASSO-PCR applies sparsity to principal-component scores rather than voxels, selecting distributed components instead of isolated speckles. Shrinkage deliberately introduces bias — the training-data fit is no longer optimal — in exchange for lower variance: more stable, more reproducible, and often more accurate weights in new data.

Beyond these workhorses, decision trees segment feature space into regions with distinct predicted values (nonparametric, nonlinear, fast), and random forests average many trees trained on random subsamples. Ensemble strategies — bagging, boosting, stochastic gradient descent, dropout — build robustness by averaging over stochastically perturbed models. Extensions such as partial least squares (PLS) and canonical correlation analysis (CCA) decode multiple outcomes simultaneously by finding maximally covarying latent components of brain and behavior; encoding–decoding models predict each voxel’s activity from rich stimulus feature spaces and then invert the mapping; and deep neural networks learn hierarchies of intermediate features directly from data (Chapter 42). Whatever the algorithm, the principles above are constant: define a loss, control complexity, and evaluate on data the model has never seen.

Hands-on tutorial¶

These exercises make the chapter’s two central lessons concrete with simulated data: (1) test error is U-shaped in model complexity, and (2) any data-dependent choice made outside cross-validation poisons the estimate of accuracy. The full labs extend both, adding lasso vs. ridge coefficient paths and ROC curves.

Step 1 — Overfitting and the U-shaped test error curve. We fit polynomials of increasing degree to noisy data from a smooth true function. Training error always falls; test error falls, bottoms out near the true complexity, then rises.

% Requires only base MATLAB
rng(1);                                     % seed for reproducibility
f = @(x) sin(2*x);                          % true function
xtr = 4*rand(30,1); ytr = f(xtr) + .4*randn(30,1);   % training set: 30 points, noise SD = 0.4
xte = 4*rand(200,1); yte = f(xte) + .4*randn(200,1); % test set: 200 independent points

degrees = 1:12;                             % complexity levels to sweep
for d = degrees
    p = polyfit(xtr, ytr, d);               % fit on training data only
    train_mse(d) = mean((ytr - polyval(p, xtr)).^2);
    test_mse(d)  = mean((yte - polyval(p, xte)).^2);
end

figure; plot(degrees, train_mse, 'b.-', degrees, test_mse, 'r.-');
set(gca, 'YScale', 'log');                  % log scale keeps the U-shape visible
legend({'Training error' 'Test error'});
xlabel('Polynomial degree (model complexity)'); ylabel('MSE (log scale)');

import numpy as np, matplotlib.pyplot as plt

rng = np.random.default_rng(1)                       # seed for reproducibility
f = lambda x: np.sin(2 * x)                          # true function
xtr = 4 * rng.random(30); ytr = f(xtr) + .4 * rng.standard_normal(30)    # training set: 30 points, noise SD = 0.4
xte = 4 * rng.random(200); yte = f(xte) + .4 * rng.standard_normal(200)  # test set: 200 independent points

degrees = np.arange(1, 13)                           # complexity levels to sweep
train_mse, test_mse = [], []
for d in degrees:
    p = np.polyfit(xtr, ytr, d)                      # fit on training data only
    train_mse.append(np.mean((ytr - np.polyval(p, xtr)) ** 2))
    test_mse.append(np.mean((yte - np.polyval(p, xte)) ** 2))

plt.plot(degrees, train_mse, "b.-", label="Training error")
plt.plot(degrees, test_mse, "r.-", label="Test error")
plt.yscale("log")                                    # log scale keeps the U-shape visible
plt.xlabel("Polynomial degree (model complexity)"); plt.ylabel("MSE (log scale)")
plt.legend()

Example output:

Step 2 — Cross-validation done wrong vs. right. The features here are pure noise — there is nothing to find. Selecting the most outcome-correlated features on the full dataset before cross-validating leaks test information into training and yields impressively wrong accuracy; selecting inside each fold gives the honest answer (~50%).

rng(7);                     % seed for reproducibility
n = 50; p = 2000; k = 20;   % n = participants, p = noise features, k = features to select
X = randn(n, p); y = [ones(n/2,1); -ones(n/2,1)];  % features are pure noise

% WRONG: select features using ALL data, then cross-validate
r = corr(X, y); [~, sel] = maxk(abs(r), k);
acc_wrong = 1 - crossval('mcr', X(:, sel), y, ...
    'Predfun', @(xt,yt,xv) sign(xv * (xt \ yt)), 'KFold', 10);

% RIGHT: feature selection re-done inside every training fold
cv = cvpartition(n, 'KFold', 10); correct = 0;
for i = 1:cv.NumTestSets
    tr = training(cv, i); te = test(cv, i);
    r = corr(X(tr,:), y(tr)); [~, sel] = maxk(abs(r), k);
    yhat = sign(X(te, sel) * (X(tr, sel) \ y(tr)));
    correct = correct + sum(yhat == y(te));
end
acc_right = correct / n;
fprintf('Wrong CV: %.0f%%  Right CV: %.0f%%  (chance = 50%%)\n', ...
    100*acc_wrong, 100*acc_right);

import numpy as np
from sklearn.model_selection import StratifiedKFold

rng = np.random.default_rng(7)                     # seed for reproducibility
n, p, k = 50, 2000, 20                             # n = participants, p = noise features, k = features to select
X = rng.standard_normal((n, p))                    # features are pure noise
y = np.repeat([1, -1], n // 2)

def top_k(Xa, ya, k):                              # |correlation| screening
    r = (Xa - Xa.mean(0)).T @ (ya - ya.mean())
    r /= Xa.std(0) * ya.std() * len(ya)
    return np.argsort(np.abs(r))[-k:]

def cv_accuracy(select_first):
    sel_all = top_k(X, y, k) if select_first else None
    correct = 0
    for tr, te in StratifiedKFold(10, shuffle=True, random_state=0).split(X, y):
        sel = sel_all if select_first else top_k(X[tr], y[tr], k)
        w = np.linalg.lstsq(X[tr][:, sel], y[tr], rcond=None)[0]
        correct += np.sum(np.sign(X[te][:, sel] @ w) == y[te])
    return correct / n

print(f"Wrong CV: {cv_accuracy(True):.0%}   Right CV: {cv_accuracy(False):.0%}"
      "   (chance = 50%)")

Example output:

Wrong CV: 88%   Right CV: 54%   (chance = 50%)

The gap between the two numbers is the optimism purchased by leakage — often 25–40 percentage points on pure noise. The same mistake with real fMRI data inflates accuracy just as silently. The full labs go on to trace ridge and lasso coefficient paths across the regularization strength $\lambda$, and to build ROC curves and AUC for informative vs. uninformative classifiers.

Thought questions¶

1. In a typical fMRI decoding study, $p \approx 100{,}000$ voxels and $n \approx 30$ participants. Explain, using the loss-plus-penalty framework $\mathcal{L}(w) = E(w) + \lambda R(w)$, why an unregularized linear model is not merely suboptimal here but mathematically ill-posed — and what different choices of $R(w)$ implicitly assume about how predictive signal is distributed across the brain.

2. Ridge regression deliberately produces biased coefficient estimates, yet often predicts new data better than unbiased OLS. Reconcile this with the classical statistical emphasis on unbiasedness, using the bias–variance decomposition. When would you expect the advantage of shrinkage to be largest, and when negligible?

3. A colleague reports 85% cross-validated accuracy discriminating patients from controls. List at least three distinct ways information could have leaked from test folds into training decisions in a realistic fMRI pipeline (consider preprocessing, feature selection, hyperparameter tuning, and participant-level structure), and rank them by how much optimism each would likely produce.

4. LASSO applied to voxels tends to select scattered “speckles” — one arbitrary voxel per correlated cluster — while elastic net and LASSO-PCR spread weight over correlated sets. If two models achieve identical prediction accuracy, what scientific and practical considerations would lead you to prefer one weight map over the other?

5. The SVM hyperparameter $C$ and the regularization coefficient $\lambda$ both govern model complexity, and both are often tuned by cross-validation. If the tuning itself uses the same folds that produce the final accuracy estimate, what goes wrong, and how does the logic parallel the feature-selection leakage demonstrated in the tutorial?

Quiz yourself¶

⌂ Back to home · Table of contents · How to use this site