| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 6: Black-Box Methods — Neural Networks, Support Vector
Machines, and Ensembles
Chapter 6: Black-Box Methods — Neural Networks, Support Vector Machines, and Ensembles
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
August 26, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly) # interactive figures and ALL 3-D graphics
library(caret)
library(kernlab)
library(e1071)
library(ranger)How this chapter uses graphics
Every two-dimensional figure is drawn with
ggplot2and rendered statically. Immediately after each one, the equivalentplot_ly()code appears in a chunk markedeval=FALSE, echo=TRUE.Every three-dimensional figure is drawn with
plot_ly()and evaluated. This chapter leans on that heavily and for a specific reason: its central objects are surfaces. A decision boundary is a level set of a posterior surface; a loss landscape is a surface over weight space; a hyperparameter grid is a surface over \((C,\gamma)\). Each is unreadable from a single fixed viewpoint.
1 Learning objectives
After completing this chapter you will be able to:
- Write down the forward pass of a feed-forward network in matrix form and count its floating-point operations.
- Derive backpropagation from the chain rule, implement it, and verify the gradients numerically.
- Explain the vanishing-gradient problem quantitatively, and justify why ReLU-family activations replaced sigmoid and tanh.
- State the universal approximation theorem and explain what it does and does not promise.
- Evaluate a regression model with metrics that detect bias and scale error, and explain why correlation cannot.
- Design an honest validation protocol for time-series data using forward chaining.
- Derive the SVM maximum-margin primal, its dual, and the KKT condition that makes the solution sparse in support vectors.
- State Mercer’s condition, explain the kernel trick, and describe the feature space induced by the RBF kernel.
- Analyze SVM training complexity and explain why SVMs stop scaling where tree ensembles do not.
- Distinguish bagging, AdaBoost, and gradient boosting by their loss functions and what each reduces, and derive AdaBoost’s weight-update coefficient.
Estimated time: 12–15 hours including exercises. Prerequisites: Chapter 3 (matrix computing, conditioning), Chapter 4 (kernels via kernel PCA), and Chapter 5 — especially the evaluation apparatus (§5.3), the leakage taxonomy (§5.4.1), and the ensemble variance identity (§5.22.1), all of which this chapter uses without re-deriving.
2 PART I — NEURAL NETWORKS
3 From biological to artificial neurons
An artificial neural network is a directed graph of simple computational units. The biological analogy — dendrites collect signals, the soma aggregates them, the axon fires if a threshold is crossed — motivated the architecture but does not constrain it; modern networks are best understood as compositions of parameterized affine maps and nonlinearities.
A single unit takes inputs \(x_1,\dots,x_n\), weights them, adds a bias, and passes the result through an activation function \(f\):
\[\boxed{\;y(\mathbf{x})=f\!\left(b+\sum_{i=1}^{n}w_i x_i\right)=f\big(\mathbf{w}^\top\mathbf{x}+b\big)\;}\]
The weights \(w_i\) set the relative influence of each input; the bias \(b\) shifts the activation left or right, which is what lets the unit fire at a threshold other than zero. Without a bias, every decision boundary would be forced through the origin.
Three components define a network:
- Activation function — the nonlinearity applied at each unit.
- Topology — how many layers, how many units per layer, how connected.
- Training algorithm — how the weights are estimated.
Common misconception: “a neural network is just a brain simulation.” The analogy ends quickly. Biological neurons are stochastic, spike in time, and learn by mechanisms nothing like gradient descent. What makes an artificial network work is that it is a differentiable function approximator — a composition \(f_L\circ\cdots\circ f_1\) whose parameters can be fitted by following a gradient. Understanding it as calculus rather than neuroscience is what lets you reason about why it succeeds and fails.
4 Activation functions
Without a nonlinearity, a network collapses. If every \(f\) were the identity, then \(W_2(W_1\mathbf{x}+\mathbf{b}_1)+\mathbf{b}_2=(W_2W_1)\mathbf{x}+(W_2\mathbf{b}_1+\mathbf{b}_2)\) — a single affine map, no matter how many layers. Depth buys nothing without nonlinearity.
4.1 The classical family, and its derivatives
\[ \begin{aligned} \textbf{Threshold: }\ & f(x)=\mathbb{1}\{x\ge 0\} & f'(x)&=0\ \ (x\ne 0)\\[1mm] \textbf{Sigmoid: }\ & \sigma(x)=\frac{1}{1+e^{-x}} & \sigma'(x)&=\sigma(x)\big(1-\sigma(x)\big)\ \le\ \tfrac14\\[1mm] \textbf{Tanh: }\ & \tanh(x)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}} & \tanh'(x)&=1-\tanh^2(x)\ \le\ 1\\[1mm] \textbf{Gaussian (RBF): }\ & f(x)=e^{-x^2/2} & f'(x)&=-x\,e^{-x^2/2}\\[1mm] \textbf{Linear: }\ & f(x)=x & f'(x)&=1 \end{aligned} \]
The threshold function is instructive and unusable: its derivative is zero everywhere it is defined, so gradient-based training has nothing to follow. That single fact is why smooth activations exist.
4.2 The saturation problem
Look at \(\sigma'\le\frac14\). Backpropagation multiplies one derivative per layer (§6.6), so a gradient reaching layer \(\ell\) from a network of depth \(L\) carries a factor
\[\prod_{k=\ell}^{L}f'(z_k)\ \le\ \left(\tfrac14\right)^{L-\ell+1}.\]
After ten sigmoid layers the gradient is attenuated by at most \(4^{-10}\approx 10^{-6}\) — and that is the best case, at \(x=0\). In the saturated regime (\(|x|>4\)), \(\sigma'<0.018\) and the attenuation is catastrophic. This is the vanishing gradient problem, and it is why deep networks were effectively untrainable before 2010.
4.3 The modern family
\[ \begin{aligned} \textbf{ReLU: }\ & f(x)=\max(0,x) & f'(x)&=\mathbb{1}\{x>0\}\in\{0,1\}\\[1mm] \textbf{Leaky ReLU: }\ & f(x)=\max(\alpha x,\ x),\ \alpha\approx0.01 & f'(x)&\in\{\alpha,1\}\\[1mm] \textbf{ELU: }\ & f(x)=\begin{cases}x & x>0\\ \alpha(e^x-1) & x\le0\end{cases} & f'(x)&=\begin{cases}1 & x>0\\ f(x)+\alpha & x\le 0\end{cases}\\[1mm] \textbf{GELU: }\ & f(x)=x\,\Phi(x) & f'(x)&=\Phi(x)+x\,\phi(x) \end{aligned} \]
ReLU’s derivative is exactly 1 on the positive half-line, so the product above does not shrink along active paths — the gradient flows undiminished through arbitrarily many layers. That single property, not any biological motivation, is why ReLU became the default in 2011 and remains so.
Its cost is the dying ReLU: a unit whose pre-activation is negative for every training example has zero gradient forever and never recovers. Leaky ReLU, ELU, and GELU all address this by keeping the derivative non-zero on the negative side.
For a classification output layer, the softmax converts scores to a probability vector:
\[\mathrm{softmax}(\mathbf{z})_k=\frac{e^{z_k}}{\sum_{j=1}^{K}e^{z_j}}, \qquad \frac{\partial\,\mathrm{softmax}_k}{\partial z_j}=\mathrm{softmax}_k\big(\delta_{kj}-\mathrm{softmax}_j\big).\]
In practice compute it as \(z_k-\log\sum_j e^{z_j}\) with the maximum subtracted (log-sum-exp, Chapter 5, §5.14) — the direct form overflows for \(z_k>709\).
act <- list(
Sigmoid = function(x) 1 / (1 + exp(-x)),
Tanh = function(x) tanh(x),
Gaussian = function(x) exp(-x^2 / 2),
ReLU = function(x) pmax(0, x),
`Leaky ReLU` = function(x) pmax(0.05 * x, x),
ELU = function(x) ifelse(x > 0, x, exp(x) - 1),
GELU = function(x) x * pnorm(x))
deriv <- list(
Sigmoid = function(x) { s <- 1/(1+exp(-x)); s * (1 - s) },
Tanh = function(x) 1 - tanh(x)^2,
Gaussian = function(x) -x * exp(-x^2 / 2),
ReLU = function(x) as.numeric(x > 0),
`Leaky ReLU` = function(x) ifelse(x > 0, 1, 0.05),
ELU = function(x) ifelse(x > 0, 1, exp(x)),
GELU = function(x) pnorm(x) + x * dnorm(x))
xg <- seq(-6, 6, length.out = 600)
act_df <- bind_rows(lapply(names(act), \(n)
data.frame(x = xg, y = act[[n]](xg), fn = n, panel = "Activation f(x)")))
der_df <- bind_rows(lapply(names(deriv), \(n)
data.frame(x = xg, y = deriv[[n]](xg), fn = n, panel = "Derivative f'(x)")))
bind_rows(act_df, der_df) |>
mutate(fn = factor(fn, levels = names(act))) |>
ggplot(aes(x, y, color = fn)) +
geom_hline(yintercept = 0, color = "grey75") +
geom_line(linewidth = 0.9) +
facet_wrap(~ panel, scales = "free_y") +
scale_color_brewer(palette = "Dark2") +
coord_cartesian(ylim = c(-1.5, 3)) +
labs(title = "Activation functions and their derivatives",
subtitle = "The right panel is what backpropagation multiplies. Sigmoid never exceeds 1/4; ReLU is exactly 1 where active",
x = "Pre-activation z", y = NULL, color = NULL) +
theme_dspa(11)
# --- Interactive equivalent ------------------------------------------------
p <- plot_ly()
for (n in names(act))
p <- add_trace(p, x = xg, y = act[[n]](xg), type = "scatter",
mode = "lines", name = n)
p |> layout(title = "Various activation functions",
xaxis = list(title = "Input signal"),
yaxis = list(title = "Output signal", range = c(-1.5, 2)),
legend = list(orientation = "h"))The attenuation is a two-variable phenomenon — depth \(\times\) activation — so it is worth seeing as a surface:
depths <- 1:20
acts3 <- c("Sigmoid", "Tanh", "ReLU")
# Best-case attenuation: the LARGEST derivative each activation can contribute
best <- c(Sigmoid = 0.25, Tanh = 1.0, ReLU = 1.0)
# Typical-case: the mean derivative over a standard-normal pre-activation
set.seed(11); zsamp <- rnorm(20000)
typ <- vapply(acts3, \(n) mean(deriv[[n]](zsamp)), numeric(1))
Zatt <- outer(seq_along(acts3), depths, \(i, d) typ[i]^d)
plot_ly(x = depths, y = acts3, z = log10(Zatt), type = "surface",
colorscale = "Viridis",
colorbar = list(title = "log10(gradient\nattenuation)")) |>
layout(title = "Gradient attenuation by depth and activation (log scale)",
scene = list(xaxis = list(title = "Layers traversed"),
yaxis = list(title = "Activation"),
zaxis = list(title = "log10 attenuation factor")))data.frame(activation = acts3,
mean_derivative = round(typ, 4),
attenuation_10_layers = signif(typ^10, 3),
attenuation_20_layers = signif(typ^20, 3))At 20 layers a sigmoid network attenuates the gradient by roughly ten orders of magnitude; a ReLU network by a factor of two or three. That difference is the whole reason deep learning became possible.
5 Network topology
Inputs are fixed by the data and outputs by the task. Everything between — depth and width — is a modelling choice.
In a feed-forward network information flows one way, from input to output. A recurrent network adds connections backwards in time, giving the model a short-term memory over sequences — the architecture behind language and time-series models (Chapter 12, Chapter 14).
Bias nodes are constant-valued units (typically 1) connected forward into every unit of the next layer. They are what supply the \(b\) in \(f(\mathbf{w}^\top\mathbf{x}+b)\), and without them every activation would be pinned to fire at zero.
5.1 The universal approximation theorem
Universal approximation (Cybenko, 1989; Hornik, 1991). Let \(f\) be a continuous function on a compact set \(K\subset\mathbb{R}^d\) and let \(\sigma\) be any non-polynomial continuous activation. Then for every \(\varepsilon>0\) there exist \(N\), weights \(\mathbf{w}_j\), biases \(b_j\), and coefficients \(\alpha_j\) such that \[\sup_{\mathbf{x}\in K}\left|f(\mathbf{x})-\sum_{j=1}^{N}\alpha_j\,\sigma\big(\mathbf{w}_j^\top\mathbf{x}+b_j\big)\right|<\varepsilon .\]
A single hidden layer of sufficient width can approximate any continuous function arbitrarily well. Note carefully what this does not say:
- It gives no bound on \(N\). The required width can be exponential in \(d\).
- It says nothing about finding the weights — the theorem is about existence, and the optimization problem is non-convex.
- It says nothing about generalization from finite data.
- Depth is not covered. Deep networks are used because they represent many functions with exponentially fewer units than a shallow one would need, not because shallow ones cannot in principle.
Common misconception: “universal approximation means neural networks can learn anything.” Existence is not attainability. The theorem guarantees a good approximator exists in the hypothesis class; it says nothing about whether gradient descent on finite noisy data will find it, or whether the one it finds will generalize. Every practical difficulty in training — local minima, initialization, learning rates, overfitting — lives in the gap between existence and attainability.
library(MASS)
# A deliberately wiggly target on [0, 1]
f_target <- function(x) sin(6 * pi * x) * exp(-2 * x) + 0.4 * x
set.seed(21)
n_uat <- 400
x_uat <- sort(runif(n_uat))
y_uat <- f_target(x_uat) # noiseless: this is an APPROXIMATION test
## OLD:
# fit_width <- function(N, seed = 31) {
# set.seed(seed)
# # Random hidden layer (fixed), least squares on the output weights.
# # This is a "random features" network: it isolates APPROXIMATION capacity
# # from the difficulty of optimizing the hidden layer.
# W <- runif(N, -25, 25); b <- runif(N, -12, 12)
# H <- tanh(outer(x_uat, W) + matrix(b, n_uat, N, byrow = TRUE))
# H <- cbind(1, H)
# # beta <- qr.solve(H, y_uat) # QR, not normal equations (Ch. 3)
# # The design matrix H (random features + intercept) can be rank‑deficient or
# # numerically singular, especially for small N (like 3 or 10) or when the random
# # features are highly correlated.
# # The qr.solve function is more stable than solve, but it can still fail
# # A robust fix is to use the Moore‑Penrose generalized inverse (MASS::ginv) or
# # add a small ridge penalty to ensure the matrix is invertible:
# # add a tiny ridge penalty
# # lambda <- 1e-8 # small positive value
# # beta <- solve(crossprod(H) + lambda * diag(ncol(H)), crossprod(H, y_uat))
# # This makes the matrix H'H + λI positive definite and thus non‑singular.
# # This gives a least‑squares solution even when H is rank‑deficient.
# beta <- ginv(H) %*% y_uat
# list(pred = as.vector(H %*% beta),
# rmse = sqrt(mean((as.vector(H %*% beta) - y_uat)^2)))
# }
fit_width <- function(N, seed = 31) {
set.seed(seed)
W <- runif(N, -25, 25)
bias <- runif(N, -12, 12) # avoid name clash with function 'b'
H <- tanh(outer(x_uat, W) + matrix(bias, n_uat, N, byrow = TRUE))
H <- cbind(1, H)
beta <- MASS::ginv(H) %*% y_uat # robust to rank deficiency
pred <- as.vector(H %*% beta)
rmse <- sqrt(mean((pred - y_uat)^2))
list(pred = pred, rmse = rmse)
}
widths <- c(3, 10, 40, 150)
uat <- bind_rows(lapply(widths, \(N) {
f <- fit_width(N)
data.frame(x = x_uat, y = f$pred,
panel = sprintf("N = %d (RMSE %.4f)", N, f$rmse))
}))
ggplot(uat, aes(x, y)) +
geom_line(data = data.frame(x = x_uat, y = y_uat), aes(x, y),
color = "grey35", linewidth = 1.1, inherit.aes = FALSE) +
geom_line(color = "firebrick", linewidth = 0.9) +
facet_wrap(~ panel, nrow = 1) +
labs(title = "Universal approximation: error falls as width grows",
subtitle = "Grey: target function. Red: single hidden layer of N tanh units",
x = "x", y = "f(x)") +
theme_dspa(10)
rates <- data.frame(N = c(2, 4, 8, 16, 32, 64, 128, 256))
rates$rmse <- vapply(rates$N, \(N) fit_width(N)$rmse, numeric(1))
ratesggplot(rates, aes(N, rmse)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
scale_x_log10() + scale_y_log10() +
labs(title = "Approximation error against hidden-layer width",
subtitle = "Log-log: a straight line indicates a power-law rate",
x = "Number of hidden units N (log scale)", y = "RMSE (log scale)") +
theme_dspa()
The error falls steadily with width, as the theorem promises — and note how many units the approximation requires even with the hidden layer chosen at random and the output solved exactly. Real training must find both layers by gradient descent, which is strictly harder.
6 The forward pass, in matrix form
For a network with layers \(\ell=1,\dots,L\), let \(\mathbf{a}^{(0)}=\mathbf{x}\) be the input. Then
\[\boxed{\;\mathbf{z}^{(\ell)}=W^{(\ell)}\mathbf{a}^{(\ell-1)}+\mathbf{b}^{(\ell)}, \qquad \mathbf{a}^{(\ell)}=f^{(\ell)}\!\big(\mathbf{z}^{(\ell)}\big)\;}\]
with \(W^{(\ell)}\in\mathbb{R}^{n_\ell\times n_{\ell-1}}\) and \(\mathbf{b}^{(\ell)}\in\mathbb{R}^{n_\ell}\). The prediction is \(\hat{\mathbf{y}}=\mathbf{a}^{(L)}\).
Complexity. Layer \(\ell\) costs one matrix–vector product, \(2n_\ell n_{\ell-1}\) flops, plus \(n_\ell\) activation evaluations. For one example,
\[\text{forward cost}=O\!\left(\sum_{\ell=1}^{L}n_\ell n_{\ell-1}\right),\]
and for a minibatch of \(m\) examples the matrix–vector products become matrix–matrix products of the same total order times \(m\) — which is why minibatching is fast: it converts memory-bound BLAS level-2 operations into compute-bound level-3 ones (Chapter 3, §3.2.2).
Parameter count is \(\sum_\ell (n_\ell n_{\ell-1}+n_\ell)\).
net_cost <- function(sizes) {
L <- length(sizes) - 1
flops <- sum(2 * sizes[-1] * sizes[-length(sizes)])
params <- sum(sizes[-1] * sizes[-length(sizes)] + sizes[-1])
c(layers = L, params = params, flops_per_example = flops)
}
rbind(`7 -> 1` = net_cost(c(7, 1)),
`7 -> 4 -> 1` = net_cost(c(7, 4, 1)),
`7 -> 4 -> 3 -> 3 -> 1` = net_cost(c(7, 4, 3, 3, 1)),
`784 -> 256 -> 128 -> 10` = net_cost(c(784, 256, 128, 10)))#> layers params flops_per_example
#> 7 -> 1 1 8 14
#> 7 -> 4 -> 1 2 37 64
#> 7 -> 4 -> 3 -> 3 -> 1 4 63 104
#> 784 -> 256 -> 128 -> 10 3 235146 4695047 Backpropagation
Training means choosing \(\{W^{(\ell)},\mathbf{b}^{(\ell)}\}\) to minimize a loss \(\mathcal{L}\). Gradient descent needs \(\partial\mathcal{L}/\partial W^{(\ell)}\) for every layer, and backpropagation is the chain rule organized so that all of them cost about as much as one forward pass.
7.1 Derivation
Define the error signal at layer \(\ell\) as the gradient of the loss with respect to that layer’s pre-activation:
\[\boldsymbol\delta^{(\ell)}\ :=\ \frac{\partial\mathcal{L}}{\partial\mathbf{z}^{(\ell)}}\ \in\mathbb{R}^{n_\ell}.\]
Output layer. For squared loss \(\mathcal{L}=\tfrac12\lVert\mathbf{a}^{(L)}-\mathbf{y}\rVert^2\), the chain rule through the activation gives
\[\boldsymbol\delta^{(L)}=\big(\mathbf{a}^{(L)}-\mathbf{y}\big)\odot f'\!\big(\mathbf{z}^{(L)}\big),\]
where \(\odot\) is the elementwise product.
Recursion. Since \(\mathbf{z}^{(\ell+1)}=W^{(\ell+1)}f(\mathbf{z}^{(\ell)})+\mathbf{b}^{(\ell+1)}\),
\[\boxed{\;\boldsymbol\delta^{(\ell)}=\Big(\big(W^{(\ell+1)}\big)^{\!\top}\boldsymbol\delta^{(\ell+1)}\Big)\odot f'\!\big(\mathbf{z}^{(\ell)}\big)\;}\]
This is where the activation derivative enters, and where it multiplies once per layer — the mechanism behind §6.3.2.
Parameter gradients. Because \(\partial z_j^{(\ell)}/\partial W_{jk}^{(\ell)}=a_k^{(\ell-1)}\),
\[\frac{\partial\mathcal{L}}{\partial W^{(\ell)}}=\boldsymbol\delta^{(\ell)}\big(\mathbf{a}^{(\ell-1)}\big)^{\!\top}, \qquad \frac{\partial\mathcal{L}}{\partial\mathbf{b}^{(\ell)}}=\boldsymbol\delta^{(\ell)} .\]
A convenient simplification. With a softmax output and cross-entropy loss, the two Jacobians cancel and
\[\boldsymbol\delta^{(L)}=\mathbf{a}^{(L)}-\mathbf{y}\]
exactly — no activation derivative at all. The same holds for a sigmoid output with binary cross-entropy, and for a linear output with squared loss. These pairings are chosen for precisely this reason: they keep the output-layer gradient from saturating.
Complexity. The backward pass costs one transposed matrix–vector product per layer, the same \(O(\sum_\ell n_\ell n_{\ell-1})\) as the forward pass. All gradients cost roughly twice a forward pass, regardless of the number of parameters — which is what makes networks with \(10^9\) weights trainable at all. The naive alternative, finite differences, would need one forward pass per parameter.
7.2 Implementation from scratch
# A minimal feed-forward network: fully vectorized over a minibatch.
# Columns of A are examples, so A^(l) is (n_l x m).
nn_init <- function(sizes, seed = 1) {
set.seed(seed)
L <- length(sizes) - 1
# He initialization for ReLU: Var(w) = 2 / fan_in keeps activation scale stable
W <- lapply(seq_len(L), \(l) matrix(rnorm(sizes[l+1] * sizes[l],
sd = sqrt(2 / sizes[l])),
sizes[l+1], sizes[l]))
b <- lapply(seq_len(L), \(l) matrix(0, sizes[l+1], 1))
list(W = W, b = b, sizes = sizes, L = L)
}
relu <- function(z) pmax(0, z)
d_relu <- function(z) (z > 0) * 1
# nn_forward <- function(net, X) {
# A <- list(X); Z <- list()
# print(dim(net$W[[l]]))
# print(dim(A[[l]]))
# for (l in seq_len(net$L)) {
# Z[[l]] <- net$W[[l]] %*% A[[l]] + matrix(net$b[[l]], nrow(net$b[[l]]), ncol(X))
# # Hidden layers: ReLU. Output layer: linear (regression).
# A[[l + 1]] <- if (l < net$L) relu(Z[[l]]) else Z[[l]]
# }
# list(A = A, Z = Z)
# }
nn_forward <- function(net, X) {
X <- as.matrix(X)
L <- net$L
A <- list()
Z <- list()
A[[1]] <- X
for (l in seq_len(L)) {
W <- as.matrix(net$W[[l]])
b <- as.vector(net$b[[l]])
# Linear transformation
Z[[l]] <- W %*% A[[l]] + b # b is recycled column-wise, but ensure length matches nrow(W)
# Activation: tanh for hidden layers, linear for output
if (l < L) {
A[[l + 1]] <- tanh(Z[[l]])
} else {
A[[l + 1]] <- Z[[l]]
}
}
list(A = A, Z = Z)
}
## OLD:
# nn_backward <- function(net, fp, Y) {
# m <- ncol(Y); L <- net$L
# dW <- vector("list", L); db <- vector("list", L)
#
# # Output layer, linear activation + squared loss: delta^(L) = a^(L) - y
# delta <- (fp$A[[L + 1]] - Y) / m
#
# for (l in L:1) {
# dW[[l]] <- delta %*% t(fp$A[[l]])
# db[[l]] <- matrix(rowSums(delta), ncol = 1)
# if (l > 1) delta <- (t(net$W[[l]]) %*% delta) * d_relu(fp$Z[[l - 1]])
# }
# list(dW = dW, db = db)
# }
nn_backward <- function(net, fp, Y) {
L <- net$L
m <- ncol(Y)
# Ensure Y is a matrix
Y <- as.matrix(Y)
dW <- vector("list", L)
db <- vector("list", L)
# Output layer error (squared loss + linear output)
delta <- (fp$A[[L + 1]] - Y) / m # shape: n_L x m
for (l in L:1) {
# Ensure matrices for multiplication
A_l <- as.matrix(fp$A[[l]]) # n_l x m
delta_mat <- as.matrix(delta) # n_{l+1} x m
# Gradient w.r.t. weights
dW[[l]] <- delta_mat %*% t(A_l) # n_{l+1} x n_l
# Gradient w.r.t. bias (sum over samples)
db[[l]] <- matrix(rowSums(delta_mat), ncol = 1)
# Backprop to previous layer if not first
if (l > 1) {
W_l <- as.matrix(net$W[[l]]) # n_{l+1} x n_l
Z_prev <- as.matrix(fp$Z[[l - 1]]) # n_l x m
# Derivative of tanh
dtanh <- 1 - tanh(Z_prev)^2
delta <- (t(W_l) %*% delta_mat) * dtanh # n_l x m
}
}
list(dW = dW, db = db)
}
nn_loss <- function(net, X, Y) {
0.5 * mean(colSums((nn_forward(net, X)$A[[net$L + 1]] - Y)^2))
}Never trust a hand-written gradient without checking it. The finite-difference approximation \(\frac{\mathcal{L}(w+\epsilon)-\mathcal{L}(w-\epsilon)}{2\epsilon}\) is \(O(\epsilon^2)\) accurate and costs two forward passes per parameter — far too slow for training, exactly right for verification.
set.seed(41)
Xc <- matrix(rnorm(4 * 25), 4, 25)
Yc <- matrix(rnorm(2 * 25), 2, 25)
netc <- nn_init(c(4, 6, 5, 2), seed = 7)
grads <- nn_backward(netc, nn_forward(netc, Xc), Yc)
num_grad <- function(net, X, Y, l, i, j, eps = 1e-6) {
up <- net; up$W[[l]][i, j] <- up$W[[l]][i, j] + eps
dn <- net; dn$W[[l]][i, j] <- dn$W[[l]][i, j] - eps
(nn_loss(up, X, Y) - nn_loss(dn, X, Y)) / (2 * eps)
}
set.seed(43)
checks <- do.call(rbind, lapply(1:8, function(k) {
l <- sample(netc$L, 1)
i <- sample(nrow(netc$W[[l]]), 1); j <- sample(ncol(netc$W[[l]]), 1)
a <- grads$dW[[l]][i, j]; n <- num_grad(netc, Xc, Yc, l, i, j)
data.frame(layer = l, i = i, j = j,
analytic = a, numeric = n,
rel_error = abs(a - n) / max(1e-12, abs(a) + abs(n)))
}))
checks |> mutate(across(where(is.numeric), \(z) signif(z, 6)))#> max_relative_error
#> 8.38e-10Relative errors near \(10^{-8}\) confirm the derivation. Anything above \(10^{-5}\) signals a bug — and gradient checking is the first thing to run whenever a network refuses to learn.
8 Optimization
Given gradients, the update rule matters as much as the gradient itself.
\[ \begin{aligned} \textbf{SGD: }\quad & \theta_{t+1}=\theta_t-\eta\,g_t\\[2mm] \textbf{Momentum: }\quad & v_{t+1}=\beta v_t+g_t, & \theta_{t+1}&=\theta_t-\eta\,v_{t+1}\\[2mm] \textbf{RMSProp: }\quad & s_{t+1}=\rho s_t+(1-\rho)g_t^2, & \theta_{t+1}&=\theta_t-\frac{\eta}{\sqrt{s_{t+1}}+\epsilon}\,g_t\\[2mm] \textbf{Adam: }\quad & \hat m_{t+1}=\frac{\beta_1 m_t+(1-\beta_1)g_t}{1-\beta_1^{t+1}},\ \ \hat s_{t+1}=\frac{\beta_2 s_t+(1-\beta_2)g_t^2}{1-\beta_2^{t+1}}, & \theta_{t+1}&=\theta_t-\frac{\eta\,\hat m_{t+1}}{\sqrt{\hat s_{t+1}}+\epsilon} \end{aligned} \]
Momentum accumulates a velocity, which damps oscillation across a narrow valley and accelerates along its floor — precisely the ill-conditioned geometry of Chapter 3, §3.7, now appearing in weight space. RMSProp divides by a running root-mean-square of the gradient, giving each parameter its own effective step size. Adam combines both, with bias correction for the initialization at zero; defaults \(\beta_1=0.9\), \(\beta_2=0.999\), \(\eta=10^{-3}\) work across a remarkable range of problems.
The learning rate \(\eta\) is the single most consequential hyperparameter. Too small and training crawls; too large and the iterates diverge. The threshold is set by the curvature: for a quadratic with Hessian \(H\), gradient descent converges iff \(\eta<2/\lambda_{\max}(H)\).
# An ill-conditioned quadratic bowl: kappa = 20
Hq <- diag(c(20, 1))
f_q <- function(th) 0.5 * drop(t(th) %*% Hq %*% th)
g_q <- function(th) drop(Hq %*% th)
run_opt <- function(method, eta, steps = 60, th0 = c(-4, 4)) {
th <- th0; v <- c(0, 0); s <- c(0, 0); m <- c(0, 0)
path <- matrix(NA_real_, steps + 1, 2); path[1, ] <- th
for (t in seq_len(steps)) {
g <- g_q(th)
th <- switch(method,
sgd = th - eta * g,
momentum = { v <<- 0.9 * v + g; th - eta * v },
adam = { m <<- 0.9 * m + 0.1 * g; s <<- 0.999 * s + 0.001 * g^2
mh <- m / (1 - 0.9^t); sh <- s / (1 - 0.999^t)
th - eta * mh / (sqrt(sh) + 1e-8) })
path[t + 1, ] <- th
}
data.frame(step = 0:steps, w1 = path[, 1], w2 = path[, 2],
loss = apply(path, 1, f_q), method = method)
}
paths <- bind_rows(run_opt("sgd", 0.05), run_opt("momentum", 0.01),
run_opt("adam", 0.30))
c(max_stable_eta_for_SGD = 2 / max(diag(Hq)))#> max_stable_eta_for_SGD
#> 0.1The loss landscape and the paths across it are a surface plus curves on it — so this figure is interactive:
w1g <- seq(-5, 5, length.out = 90); w2g <- seq(-5, 5, length.out = 90)
Lz <- outer(w1g, w2g, \(a, b) 0.5 * (20 * a^2 + b^2))
p <- plot_ly() |>
add_surface(x = w2g, y = w1g, z = Lz, opacity = 0.7, showscale = FALSE,
colorscale = "Viridis")
for (m in unique(paths$method)) {
d <- filter(paths, method == m)
p <- add_trace(p, x = d$w2, y = d$w1, z = d$loss + 1,
type = "scatter3d", mode = "lines+markers", name = m,
line = list(width = 6), marker = list(size = 2))
}
p |> layout(title = "Gradient descent on an ill-conditioned quadratic (condition number 20)",
scene = list(xaxis = list(title = "w2"), yaxis = list(title = "w1"),
zaxis = list(title = "Loss")))ggplot(paths, aes(step, loss, color = method)) +
geom_line(linewidth = 1) +
scale_y_log10() +
scale_color_manual(values = c(sgd = "#3B7DD8", momentum = "#7FB069",
adam = "#D8433B")) +
labs(title = "Convergence on an ill-conditioned quadratic",
subtitle = "Plain SGD zig-zags across the narrow direction; momentum and Adam do not",
x = "Iteration", y = "Loss (log scale)", color = NULL) +
theme_dspa()
9 Regularization
A network with more parameters than observations can interpolate the training data exactly. Four standard controls:
Weight decay (\(L_2\)). Add \(\frac{\lambda}{2}\lVert\theta\rVert^2\) to the loss; the gradient gains \(\lambda\theta\), shrinking every weight toward zero each step. This is ridge regression’s penalty (Chapter 3, Problem 8) applied to network weights.
Dropout. During training, zero each hidden unit independently with probability \(p\) and rescale the survivors by \(1/(1-p)\). At test time use the full network. Dropout approximates averaging over an exponential family of thinned sub-networks — an implicit ensemble (Chapter 5, §5.22).
Early stopping. Monitor validation loss and stop when it starts rising. Under gradient descent on a quadratic this is provably similar to \(L_2\) regularization, with the number of steps playing the role of \(1/\lambda\).
Batch normalization. Standardize each layer’s pre-activations within the minibatch. It stabilizes the scale of gradients across layers and permits substantially larger learning rates.
set.seed(53)
n_r <- 60
x_r <- sort(runif(n_r, 0, 1))
y_r <- f_target(x_r) + rnorm(n_r, sd = 0.12) # NOISY this time
ridge_fit <- function(lambda, N = 120, seed = 61) {
set.seed(seed)
W <- runif(N, -25, 25); b <- runif(N, -12, 12)
H <- cbind(1, tanh(outer(x_r, W) + matrix(b, n_r, N, byrow = TRUE)))
# Ridge solution: (H'H + lambda I)^{-1} H'y -- see Chapter 3, Section 3.6.4
beta <- solve(crossprod(H) + lambda * diag(ncol(H)), crossprod(H, y_r))
xs <- seq(0, 1, length.out = 400)
Hs <- cbind(1, tanh(outer(xs, W) + matrix(b, length(xs), N, byrow = TRUE)))
list(x = xs, y = as.vector(Hs %*% beta),
train_rmse = sqrt(mean((as.vector(H %*% beta) - y_r)^2)),
true_rmse = sqrt(mean((as.vector(Hs %*% beta) - f_target(xs))^2)))
}
lams <- c(1e-8, 1e-4, 1e-1, 10)
reg_df <- bind_rows(lapply(lams, \(l) {
f <- ridge_fit(l)
data.frame(x = f$x, y = f$y,
panel = sprintf("lambda = %g\ntrain %.3f | true %.3f",
l, f$train_rmse, f$true_rmse))
}))
ggplot(reg_df, aes(x, y)) +
geom_line(data = data.frame(x = seq(0, 1, length.out = 400),
y = f_target(seq(0, 1, length.out = 400))),
aes(x, y), color = "grey35", linewidth = 1, inherit.aes = FALSE) +
geom_point(data = data.frame(x = x_r, y = y_r), aes(x, y),
color = "grey20", size = 0.9, alpha = 0.6, inherit.aes = FALSE) +
geom_line(color = "firebrick", linewidth = 0.9) +
facet_wrap(~ panel, nrow = 1) +
coord_cartesian(ylim = c(-1, 1.2)) +
labs(title = "Weight decay controls the bias-variance tradeoff",
subtitle = "Grey line: truth. Points: noisy training data. Red: fitted network",
x = "x", y = "f(x)") +
theme_dspa(9)
Note the two RMSE figures in each panel. At \(\lambda=10^{-8}\) the training error is smallest and the error against the truth is largest — the signature of overfitting, and the reason training loss can never be used to select a regularization strength.
10 Case study: Google Trends and the stock market
10.1 The data
The Google Trends and stock market dataset records daily search-volume indices alongside the Dow Jones Industrial Average, 2008–2009. Metadata and CSV are on the Case Studies site.
Variables. Index, Date;
search indices for Unemployment, Rental,
RealEstate, Mortgage, Jobs,
Investing; DJI_Index; and StdDJI
\(=3+(\mathrm{DJI}-11091)/1501\), the
DJI standardized by its 2005–2011 mean and standard deviation. Each of
the eight also appears as a 30-day and a 180-day moving average.
google_raw <- dspa_read(
"https://umich.instructure.com/files/416274/download?download_frd=1",
"GoogleTrends_Data.csv")
dim(google_raw); names(google_raw)[1:12]#> [1] 731 26#> [1] "Index" "Date" "Unemployment"
#> [4] "Rental" "RealEstate" "Mortgage"
#> [7] "Jobs" "Investing" "DJI_Index"
#> [10] "StdDJI" "Unemployment_30MA" "Rental_30MA"#> 'data.frame': 731 obs. of 8 variables:
#> $ Unemployment: num 1.54 1.56 1.59 1.62 1.64 1.64 1.71 1.85 1.82 1.78 ...
#> $ Rental : num 0.88 0.9 0.92 0.92 0.94 0.96 0.99 1.02 1.02 1.01 ...
#> $ RealEstate : num 0.79 0.81 0.82 0.82 0.83 0.84 0.86 0.89 0.89 0.89 ...
#> $ Mortgage : num 1 1.05 1.07 1.08 1.1 1.11 1.15 1.22 1.23 1.24 ...
#> $ Jobs : num 0.99 1.05 1.1 1.14 1.17 1.2 1.3 1.41 1.43 1.44 ...
#> $ Investing : num 0.92 0.94 0.96 0.98 0.99 0.99 1.02 1.09 1.1 1.1 ...
#> $ DJI_Index : num 13044 13044 13057 12800 12827 ...
#> $ StdDJI : num 4.3 4.3 4.31 4.14 4.16 4.16 4.16 4 4.1 4.17 ...10.2 These are daily observations of a time series
Common misconception: “randomly split the rows into 75% train and 25% test.” For a time series that is temporal leakage. Consecutive days of a search index are nearly identical, so a random split places day \(t\) in training and day \(t+1\) in test — the model is interpolating between points it has already seen, not forecasting. Reported performance can approach perfection while genuine predictive skill is close to zero.
The honest protocol is forward chaining: train on \([1,\tau]\), test on \((\tau,T]\), never shuffling. For repeated estimates, roll the origin forward.
ac <- acf(google$RealEstate, lag.max = 30, plot = FALSE)
data.frame(lag = as.numeric(ac$lag), acf = as.numeric(ac$acf)) |>
ggplot(aes(lag, acf)) +
geom_col(fill = "steelblue", width = 0.5) +
geom_hline(yintercept = c(-1, 1) * 1.96 / sqrt(nrow(google)),
linetype = "dashed", color = "firebrick") +
labs(title = "Autocorrelation of the Real Estate search index",
subtitle = sprintf("Lag-1 autocorrelation = %.3f. Adjacent days carry almost the same information",
ac$acf[2]),
x = "Lag (days)", y = "Autocorrelation") +
theme_dspa()
With a lag-1 autocorrelation this high, a randomly held-out day is nearly a duplicate of one in training.
n_g <- nrow(google)
cut_g <- floor(0.75 * n_g)
train_idx_time <- 1:cut_g # chronological
test_idx_time <- (cut_g + 1):n_g
set.seed(1234)
train_idx_rand <- sample(n_g, cut_g) # for the comparison only
test_idx_rand <- setdiff(seq_len(n_g), train_idx_rand)
c(n = n_g, train = length(train_idx_time), test = length(test_idx_time))#> n train test
#> 731 548 183Scaling constants come from the training window only — the same discipline as Chapter 5, §5.4.1, and more consequential here because min–max scaling lets a single extreme test value define the scale.
make_scaler <- function(train_df) {
mn <- sapply(train_df, min); mx <- sapply(train_df, max)
rng <- mx - mn; rng[rng == 0] <- 1
function(df) as.data.frame(Map(\(v, a, b) (v - a) / b, df, mn[names(df)],
rng[names(df)]))
}
scale_time <- make_scaler(google[train_idx_time, ])
g_tr <- scale_time(google[train_idx_time, ])
g_te <- scale_time(google[test_idx_time, ])
# Test values may fall outside [0, 1] -- and that is correct: the future can
# exceed anything seen in the training window.
round(range(g_te$RealEstate), 3)#> [1] -0.156 0.66710.3 Fitting
library(neuralnet)
preds_g <- c("Unemployment", "Rental", "Mortgage", "Jobs", "Investing",
"DJI_Index", "StdDJI")
form_g <- as.formula(paste("RealEstate ~", paste(preds_g, collapse = " + ")))
set.seed(1234)
nn1 <- neuralnet(form_g, data = g_tr, hidden = 1, linear.output = TRUE)
set.seed(1234)
nn4 <- neuralnet(form_g, data = g_tr, hidden = 4, linear.output = TRUE)
set.seed(1234)
nn433 <- neuralnet(form_g, data = g_tr, hidden = c(4, 3, 3), linear.output = TRUE)
sapply(list(`hidden = 1` = nn1, `hidden = 4` = nn4, `hidden = c(4,3,3)` = nn433),
\(m) c(SSE = m$result.matrix["error", 1],
steps = m$result.matrix["steps", 1]))#> hidden = 1 hidden = 4 hidden = c(4,3,3)
#> SSE.error 0.792759 0.486876 0.407794
#> steps.steps 1840.000000 1976.000000 1635.000000
The figure is rendered from the fitted object, so it
reflects the model actually being discussed. Error is the
training sum of squared errors and Steps the iterations to
convergence. The blue singletons are bias nodes —
constant-valued units supplying the \(b\) of \(f(\mathbf{w}^\top\mathbf{x}+b)\) (§6.4).
10.4 Evaluating a regression model
Common misconception: “a correlation above 0.9 between predicted and observed means the model is good.” Correlation is invariant to affine transformation: \(\operatorname{cor}(a+b\hat y,\ y)=\operatorname{cor}(\hat y,y)\) for any \(b>0\). A model returning \(\hat y=100+0.01\,y\) has correlation exactly 1 and is useless — it gets every value wrong, by a lot, in a perfectly ordered way.
Report RMSE and MAE (in the response’s units), \(R^2\) on held-out data (which can be negative — informatively so), and the calibration slope and intercept from regressing \(y\) on \(\hat y\): agreement means intercept 0 and slope 1, which correlation cannot detect.
set.seed(67)
y_demo <- rnorm(200, 50, 10)
yhat_bad <- 100 + 0.01 * y_demo # perfectly correlated, wildly wrong
c(correlation = cor(yhat_bad, y_demo),
RMSE = sqrt(mean((yhat_bad - y_demo)^2)),
R2 = 1 - sum((yhat_bad - y_demo)^2) / sum((y_demo - mean(y_demo))^2),
calibration_slope = unname(coef(lm(y_demo ~ yhat_bad))[2]))#> correlation RMSE R2 calibration_slope
#> 1.0000 51.5496 -25.3661 100.0000Correlation 1.0, RMSE 50, and an \(R^2\) of \(-24\). The three honest metrics all detect what correlation cannot.
reg_metrics <- function(pred, obs) {
cal <- coef(lm(obs ~ pred))
ccc <- 2 * cov(pred, obs) /
(var(pred) + var(obs) + (mean(pred) - mean(obs))^2) # Lin's CCC
c(RMSE = sqrt(mean((pred - obs)^2)),
MAE = mean(abs(pred - obs)),
R2 = 1 - sum((pred - obs)^2) / sum((obs - mean(obs))^2),
correlation = cor(pred, obs),
calib_intercept = unname(cal[1]),
calib_slope = unname(cal[2]),
CCC = ccc)
}# predict() matches by NAME; compute() matches by column POSITION, which fails
# silently if the data frame's column order differs from the formula's.
pred_nn <- function(m, newdata) as.vector(predict(m, newdata[, preds_g, drop = FALSE]))
res_g <- rbind(
`hidden = 1` = reg_metrics(pred_nn(nn1, g_te), g_te$RealEstate),
`hidden = 4` = reg_metrics(pred_nn(nn4, g_te), g_te$RealEstate),
`hidden = c(4,3,3)` = reg_metrics(pred_nn(nn433, g_te), g_te$RealEstate),
`linear model` = reg_metrics(
as.vector(predict(lm(form_g, data = g_tr), g_te)), g_te$RealEstate),
`mean baseline` = reg_metrics(
rep(mean(g_tr$RealEstate), nrow(g_te)), g_te$RealEstate))
round(res_g, 4)#> RMSE MAE R2 correlation calib_intercept calib_slope
#> hidden = 1 0.1426 0.1224 0.6068 0.9587 0.1579 0.6911
#> hidden = 4 0.1793 0.1475 0.3777 0.9710 0.1954 0.6548
#> hidden = c(4,3,3) 0.0677 0.0541 0.9113 0.9557 0.0012 0.9726
#> linear model 0.0756 0.0642 0.8895 0.9677 0.0730 0.8327
#> mean baseline 0.3958 0.3241 -2.0307 NA 0.3287 NA
#> CCC
#> hidden = 1 0.8714
#> hidden = 4 0.8227
#> hidden = c(4,3,3) 0.9549
#> linear model 0.9532
#> mean baseline 0.0000Read the calibration slope alongside the correlation. A model can rank the test days almost perfectly (high correlation) while being systematically biased, which the slope and intercept expose immediately.
10.5 How much did the random split flatter the model?
scale_rand <- make_scaler(google[train_idx_rand, ])
r_tr <- scale_rand(google[train_idx_rand, ])
r_te <- scale_rand(google[test_idx_rand, ])
set.seed(1234)
nn_rand <- neuralnet(form_g, data = r_tr, hidden = 4, linear.output = TRUE)
comparison <- rbind(
`Random split (leaky)` = reg_metrics(pred_nn(nn_rand, r_te), r_te$RealEstate),
`Chronological split (honest)` = reg_metrics(pred_nn(nn4, g_te), g_te$RealEstate))
round(comparison[, c("RMSE", "R2", "correlation", "calib_slope")], 4)#> RMSE R2 correlation calib_slope
#> Random split (leaky) 0.0414 0.9677 0.9838 0.9973
#> Chronological split (honest) 0.1793 0.3777 0.9710 0.6548The random split reports a far better model. It is the same architecture on the same data — the difference is entirely the validation protocol. On an autocorrelated series the random split measures interpolation; the chronological split measures forecasting, which is what anyone would actually want.
10.6 Rolling-origin evaluation
A single chronological cut gives one estimate. Rolling the origin forward gives several, and their spread is the uncertainty in the estimate.
rolling_origin <- function(data, n_splits = 6, min_train = 0.4, horizon = 40) {
n <- nrow(data)
starts <- round(seq(min_train * n, n - horizon, length.out = n_splits))
do.call(rbind, lapply(starts, function(tau) {
tr <- data[1:tau, ]; te <- data[(tau + 1):(tau + horizon), ]
sc <- make_scaler(tr); trs <- sc(tr); tes <- sc(te)
set.seed(1234)
m <- neuralnet(form_g, data = trs, hidden = 4, linear.output = TRUE)
p <- as.vector(predict(m, tes[, preds_g, drop = FALSE]))
data.frame(train_end = tau,
RMSE = sqrt(mean((p - tes$RealEstate)^2)),
R2 = 1 - sum((p - tes$RealEstate)^2) /
sum((tes$RealEstate - mean(tes$RealEstate))^2))
}))
}
roll <- rolling_origin(google)
roll |> mutate(across(where(is.numeric), \(z) round(z, 4)))#> mean_RMSE sd_RMSE
#> 0.0827 0.0288roll |> pivot_longer(c(RMSE, R2), names_to = "metric", values_to = "value") |>
ggplot(aes(train_end, value)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
facet_wrap(~ metric, scales = "free_y") +
labs(title = "Rolling-origin evaluation",
subtitle = "Each point trains on all days up to the origin and forecasts the next 40",
x = "Last training day", y = NULL) +
theme_dspa(10)
sc_df <- data.frame(pred = pred_nn(nn4, g_te), obs = g_te$RealEstate)
m_g <- reg_metrics(sc_df$pred, sc_df$obs)
ggplot(sc_df, aes(pred, obs)) +
geom_abline(slope = 1, intercept = 0, color = "firebrick", linewidth = 1) +
geom_abline(slope = m_g[["calib_slope"]], intercept = m_g[["calib_intercept"]],
color = "steelblue", linetype = "dashed", linewidth = 0.9) +
geom_point(alpha = 0.5, size = 1.6, color = "grey25") +
coord_fixed() +
labs(title = "Predicted vs. observed on the held-out future",
subtitle = sprintf("Red: ideal agreement. Blue dashed: calibration line (slope %.2f, intercept %.2f). RMSE %.3f, R2 %.3f",
m_g[["calib_slope"]], m_g[["calib_intercept"]],
m_g[["RMSE"]], m_g[["R2"]]),
x = "Neural network prediction", y = "Observed (scaled) Real Estate index") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_markers(x = sc_df$pred, y = sc_df$obs, name = "Data scatter") |>
add_trace(x = c(0, 1), y = c(0, 1), type = "scatter", mode = "lines",
line = list(width = 4), name = "Ideal agreement") |>
layout(title = sprintf("Observed vs. predicted Real Estate (RMSE = %.3f, R2 = %.3f)",
m_g[["RMSE"]], m_g[["R2"]]),
xaxis = list(title = "NN prediction"),
yaxis = list(title = "Observed (scaled)"),
legend = list(orientation = "h"))11 Learning \(\sqrt{\cdot}\): what a network can and cannot do
Networks appear at the interface of experimental, theoretical, and computational science — including applications to string theory. A clean test of the machinery: learn \(\sqrt{\cdot}:\mathbb{R}^{+}\to\mathbb{R}^{+}\) from samples.
set.seed(1234) # seed BEFORE generating the data
rand_data <- runif(1000, 0, 100) # already non-negative
sqrt_df <- data.frame(rand_data, sqrt_data = sqrt(rand_data))
ggplot(data.frame(x = seq(0, 100, length.out = 500)), aes(x, sqrt(x))) +
geom_line(linewidth = 1, color = "steelblue") +
geom_rug(data = data.frame(x = rand_data[1:200]), aes(x, y = 0),
sides = "b", alpha = 0.3, inherit.aes = FALSE) +
labs(title = "The target function and the training range",
subtitle = "Ticks show training inputs: all in [0, 100]",
x = "x", y = expression(sqrt(x))) +
theme_dspa()
set.seed(1234)
net_sqrt <- neuralnet(sqrt_data ~ rand_data, data = sqrt_df,
hidden = 10, threshold = 0.01, stepmax = 1e6)
c(training_SSE = net_sqrt$result.matrix["error", 1],
steps = net_sqrt$result.matrix["steps", 1])#> training_SSE.error steps.steps
#> 1.08303e-03 5.96020e+04threshold is the stopping criterion on the
partial derivatives of the error: training halts when
every one falls below it. A loose value stops early with a visibly worse
fit.
test_x <- seq(0, 200, by = 0.5) # deliberately EXCEEDS [0, 100]
test_df <- data.frame(rand_data = test_x)
pred_sq <- as.vector(predict(net_sqrt, test_df))
sq <- data.frame(x = test_x, truth = sqrt(test_x), pred = pred_sq,
region = ifelse(test_x <= 100, "Interpolation [0, 100]",
"Extrapolation (100, 200]"))
ggplot(sq, aes(x)) +
geom_line(aes(y = truth, color = "True sqrt(x)"), linewidth = 1) +
geom_line(aes(y = pred, color = "Network prediction"), linewidth = 1) +
geom_vline(xintercept = 100, linetype = "dashed", color = "grey40") +
ggplot2::annotate("text", x = 103, y = 3, hjust = 0, size = 3.2, color = "grey30",
label = "training range ends") +
scale_color_manual(values = c("True sqrt(x)" = "grey30",
"Network prediction" = "firebrick")) +
labs(title = "The network learns the function, then stops learning at the edge of its data",
x = "x", y = "y", color = NULL) +
theme_dspa()
sq |> summarise(RMSE = sqrt(mean((pred - truth)^2)),
max_abs_error = max(abs(pred - truth)), .by = region) |>
mutate(across(where(is.numeric), \(z) round(z, 4)))# --- Interactive equivalents ----------------------------------------------
plot_ly(x = ~pred_sq, y = ~sqrt(test_x), type = "scatter", mode = "markers",
name = "Predicted vs. actual") |>
add_trace(x = c(0, 15), y = c(0, 15), mode = "lines",
line = list(width = 4), name = "Ideal agreement") |>
layout(title = "Predicted vs. actual square root",
xaxis = list(title = "NN predicted", scaleanchor = "y"),
yaxis = list(title = "Actual sqrt(x)"),
legend = list(orientation = "h"))
plot_ly(x = ~test_x, y = ~sqrt(test_x), type = "scatter", mode = "lines",
name = "sqrt(x)") |>
add_trace(y = ~pred_sq, mode = "lines", name = "NN prediction") |>
layout(title = "Predicted vs. actual square root",
xaxis = list(title = "Input x"), yaxis = list(title = "Output"),
legend = list(orientation = "h"))Common misconception: “the network learned the square-root function.” It learned an approximation on \([0,100]\), where it saw data. Beyond that boundary the prediction flattens and the error grows without limit. Universal approximation (§6.4.1) is stated on a compact set for exactly this reason.
Neural networks — like the tree ensembles of Chapter 5, §5.18 — do not extrapolate. They interpolate within the convex hull of their training inputs. Any deployment where new data can fall outside the training range needs either a parametric form that encodes the right functional shape, or an explicit out-of-distribution check that refuses to predict.
12 Case study: Google Trends classification
The same data, recast as a three-class problem: is the Real Estate index in its top quartile, bottom quartile, or the middle?
# cut() with quantile breaks: the label names come from the factor itself,
# so codes and names cannot desynchronize.
q_re <- quantile(google$RealEstate, c(0.25, 0.75))
re_class <- cut(google$RealEstate,
breaks = c(-Inf, q_re[1], q_re[2], Inf),
labels = c("Low", "Middle", "High"))
table(re_class)#> re_class
#> Low Middle High
#> 190 362 179#> 25% 75%
#> 0.72 0.90# Chronological split again -- the outcome is still a time series
cls_tr <- g_tr; cls_tr$y <- re_class[train_idx_time]
cls_te <- g_te; cls_te$y <- re_class[test_idx_time]
round(rbind(train = prop.table(table(cls_tr$y)),
test = prop.table(table(cls_te$y))), 3)#> Low Middle High
#> train 0.142 0.531 0.327
#> test 0.612 0.388 0.000#> no_information_rate
#> 0.612The class balance differs sharply between the two windows — the market moved. That is a real property of the problem and the reason a random split would have hidden it.
# One-hot encoding, with column names taken FROM the factor levels
Y_ind <- model.matrix(~ y - 1, data = cls_tr)
colnames(Y_ind) <- levels(cls_tr$y)
train_nn <- cbind(cls_tr[, preds_g], Y_ind)
form_cls <- as.formula(paste(paste(levels(cls_tr$y), collapse = " + "), "~",
paste(preds_g, collapse = " + ")))
set.seed(2017)
nn_cls <- neuralnet(form_cls, data = train_nn, hidden = 4,
linear.output = FALSE, stepmax = 1e6)
predict_class <- function(m, newdata, lev) {
scores <- predict(m, newdata[, preds_g, drop = FALSE])
factor(lev[max.col(scores)], levels = lev)
}
pred_cls <- predict_class(nn_cls, cls_te, levels(cls_tr$y))
cm_g <- confusionMatrix(data = pred_cls, reference = cls_te$y)
cm_g$overall[1:6] |> round(4)#> Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
#> 0.9672 0.9306 0.9300 0.9879 0.6120
#> AccuracyPValue
#> 0.0000#> Reference
#> Prediction Low Middle High
#> Low 110 4 0
#> Middle 2 67 0
#> High 0 0 0set.seed(2017)
nn_cls2 <- neuralnet(form_cls, data = train_nn, hidden = c(4, 5),
linear.output = FALSE, stepmax = 1e6)
pred_cls2 <- predict_class(nn_cls2, cls_te, levels(cls_tr$y))
data.frame(
model = c("Majority class", "hidden = 4", "hidden = c(4, 5)"),
accuracy = round(c(NIR_g, mean(pred_cls == cls_te$y),
mean(pred_cls2 == cls_te$y)), 4),
kappa = round(c(0,
cm_g$overall[["Kappa"]],
confusionMatrix(pred_cls2, cls_te$y)$overall[["Kappa"]]), 4))A deeper network does not automatically help. Capacity without more information buys variance, not accuracy — and with a distribution shift between the training and test windows, extra flexibility fits the past more tightly rather than the future more accurately.
12.1 Beyond
neuralnet
neuralnet is well suited to teaching: small, plottable,
pure R. It fits by resilient backpropagation on the full batch, which
limits it to modest problems.
Production work uses torch or
keras3, which supply minibatch SGD and
Adam, ReLU-family activations, dropout and batch normalization,
automatic differentiation, and GPU execution. The mathematics of §6.5–§6.8
is unchanged — those libraries implement exactly the recursion derived
above. Restricted Boltzmann machines, a stochastic
alternative used for dimensionality reduction, collaborative filtering,
and generative modelling, are developed in DSPA
Appendix 12. Chapter 14 covers
convolutional and recurrent architectures.
13 PART II — SUPPORT VECTOR MACHINES
Neural networks approximate the posterior by composing nonlinearities. Support vector machines take a different route: find the separating surface with the largest margin, and reach nonlinearity through an inner-product substitution rather than through depth.
14 Separating hyperplanes
In \(\mathbb{R}^n\) a hyperplane is the solution set of a single linear equation
\[\mathbf{w}^\top\mathbf{x}+b=0,\]
where \(\mathbf{w}\) is the normal vector and \(b\) fixes the offset. In \(\mathbb{R}^3\) this is the familiar \(ax+by+cz+d=0\) with \(d=-(ax_0+by_0+cz_0)\) for any point \((x_0,y_0,z_0)\) on the plane. Such a surface is \((n-1)\)-dimensional and cuts the space in two, which makes it a binary classifier: predict by the sign of \(\mathbf{w}^\top\mathbf{x}+b\).
When two classes are linearly separable, infinitely many hyperplanes separate them. Which should we choose?
A <- c(1, 4, 3, 2, 4, 8, 6, 10, 9)
B <- c(1, 5, 3, 2, 3, 8, 8, 7, 10)
grp <- factor(c(rep("Class -1", 5), rep("Class +1", 4)))
hp <- data.frame(x = A, y = B, grp)
ggplot(hp, aes(x, y, shape = grp, color = grp)) +
geom_point(size = 3.5) +
geom_vline(xintercept = 5, color = "#3B7DD8", linetype = "dashed", linewidth = 0.9) +
geom_abline(slope = -1, intercept = 12, color = "#D8433B",
linetype = "dashed", linewidth = 0.9) +
ggplot2::annotate("text", x = 5.3, y = 9.6, label="A", color="#3B7DD8", size=5) +
ggplot2::annotate("text", x = 8.3, y = 4.2, label="B", color="#D8433B", size=5) +
scale_color_manual(values = c("Class -1" = "grey25", "Class +1" = "grey25")) +
coord_fixed() +
labs(title = "Many hyperplanes separate the same data",
subtitle = "Which one should a classifier prefer, and why?",
x = "X", y = "Y", shape = NULL, color = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = A, y = B, type = "scatter", mode = "markers", name = "Data") |>
add_lines(x = c(5, 5), y = c(1, 10), name = "Line A") |>
add_lines(x = c(10, 2), y = c(2, 10), name = "Line B") |>
layout(title = "Hyperplane (line) separation of 2-D data",
xaxis = list(title = "X", scaleanchor = "y"),
yaxis = list(title = "Y"), legend = list(orientation = "h"))15 The maximum-margin principle
15.1 The geometric margin
Encode the classes as \(y_i\in\{-1,+1\}\). The signed distance from a point \(\mathbf{x}_i\) to the hyperplane is
\[\frac{\mathbf{w}^\top\mathbf{x}_i+b}{\lVert\mathbf{w}\rVert},\]
so the distance to the correct side is \(y_i(\mathbf{w}^\top\mathbf{x}_i+b)/\lVert\mathbf{w}\rVert\). The margin is the smallest such distance over all points:
\[\gamma=\min_{i}\ \frac{y_i\big(\mathbf{w}^\top\mathbf{x}_i+b\big)}{\lVert\mathbf{w}\rVert}.\]
The pair \((\mathbf{w},b)\) is determined only up to scale — multiplying both by \(c>0\) leaves the hyperplane unchanged. Fix the scale by requiring the closest points to satisfy \(y_i(\mathbf{w}^\top\mathbf{x}_i+b)=1\). This is the canonical form, and under it
\[\gamma=\frac{1}{\lVert\mathbf{w}\rVert}, \qquad\text{so the two margin boundaries are separated by }\ \frac{2}{\lVert\mathbf{w}\rVert}.\]
Maximizing \(2/\lVert\mathbf{w}\rVert\) is minimizing \(\lVert\mathbf{w}\rVert\).
15.2 The hard-margin primal
\[\boxed{\;\min_{\mathbf{w},\,b}\ \ \frac{1}{2}\lVert\mathbf{w}\rVert^{2} \qquad\text{subject to}\qquad y_i\big(\mathbf{w}^\top\mathbf{x}_i+b\big)\ \ge\ 1,\quad i=1,\dots,n\;}\]
Why the square, and why the half? Minimizing \(\lVert\mathbf{w}\rVert\) and minimizing \(\tfrac12\lVert\mathbf{w}\rVert^2\) have the same argmin, but only the second is a quadratic program: \(\lVert\mathbf{w}\rVert=\sqrt{\mathbf{w}^\top\mathbf{w}}\) is not differentiable at the origin, whereas \(\tfrac12\mathbf{w}^\top\mathbf{w}\) is smooth and strictly convex, so the problem has a unique global minimum reachable by standard QP solvers. The \(\tfrac12\) exists purely so that \(\nabla_{\mathbf{w}}=\mathbf{w}\) rather than \(2\mathbf{w}\). This is the objective used throughout — including in the Lagrangian below.
The constraint set is the intersection of \(n\) half-spaces (convex), and the objective is strictly convex, so the problem is a convex QP with a unique solution whenever it is feasible.
Geometrically, the maximum-margin hyperplane is perpendicular to the shortest segment joining the two classes’ convex hulls, and passes through its midpoint.
g1 <- hp[hp$grp == "Class -1", c("x", "y")]
g2 <- hp[hp$grp == "Class +1", c("x", "y")]
h1 <- g1[chull(g1), ]; h2 <- g2[chull(g2), ]
ggplot() +
geom_polygon(data = h1, aes(x, y), fill = "#3B7DD8", alpha = 0.15,
color = "#3B7DD8", linewidth = 0.7) +
geom_polygon(data = h2, aes(x, y), fill = "#D8433B", alpha = 0.15,
color = "#D8433B", linewidth = 0.7) +
geom_segment(aes(x = 4, y = 5, xend = 6, yend = 8), linetype = "dashed",
color = "darkgreen", linewidth = 1) +
geom_abline(slope = -2/3, intercept = 9.833, color = "black", linewidth = 1.1) +
geom_point(data = hp, aes(x, y, shape = grp), size = 3, color = "grey20") +
geom_point(data = data.frame(x = c(4, 6), y = c(5, 8)), aes(x, y),
size = 5, shape = 21, fill = "gold", stroke = 1.2) +
coord_fixed() +
labs(title = "The maximum-margin hyperplane",
subtitle = "Shaded: convex hulls. Green dashed: shortest segment between them. Black: MMH, its perpendicular bisector. Gold: support vectors",
x = "X", y = "Y", shape = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly() |>
add_lines(x = c(6, 4), y = c(8, 5), name = "Shortest line between hulls",
line = list(dash = "dash")) |>
add_lines(x = c(10, 2), y = c(3, 8.7), name = "Maximum-margin hyperplane") |>
add_polygons(x = h1$x, y = h1$y, name = "Hull, class -1", opacity = 0.2) |>
add_polygons(x = h2$x, y = h2$y, name = "Hull, class +1", opacity = 0.2) |>
add_markers(x = A, y = B, name = "Data", marker = list(color = "black")) |>
add_markers(x = c(4, 6), y = c(5, 8), name = "Support vectors",
marker = list(size = 18, color = "gold",
line = list(color = "black", width = 2))) |>
layout(title = "Hyperplane separation of 2-D data",
xaxis = list(title = "X", scaleanchor = "y"),
yaxis = list(title = "Y"), legend = list(orientation = "h"))16 The dual problem and support vectors
16.1 Lagrangian and stationarity
Introduce multipliers \(\alpha_i\ge0\) for the \(n\) inequality constraints:
\[L_P(\mathbf{w},b,\boldsymbol\alpha)=\frac12\lVert\mathbf{w}\rVert^2-\sum_{i=1}^{n}\alpha_i\Big[y_i\big(\mathbf{w}^\top\mathbf{x}_i+b\big)-1\Big].\]
Setting the partial derivatives to zero:
\[\frac{\partial L_P}{\partial\mathbf{w}}=\mathbf{w}-\sum_{i=1}^{n}\alpha_iy_i\mathbf{x}_i=\mathbf{0} \ \Longrightarrow\ \boxed{\ \mathbf{w}=\sum_{i=1}^{n}\alpha_iy_i\mathbf{x}_i\ }\]
\[\frac{\partial L_P}{\partial b}=-\sum_{i=1}^{n}\alpha_iy_i=0 \ \Longrightarrow\ \sum_{i=1}^{n}\alpha_iy_i=0 .\]
The first identity is the representer result: the optimal normal vector is a linear combination of the training points, weighted by \(\alpha_iy_i\).
16.2 The dual
Substituting both back into \(L_P\) eliminates \(\mathbf{w}\) and \(b\) entirely:
\[\boxed{\;\max_{\boldsymbol\alpha}\ \ L_D(\boldsymbol\alpha)=\sum_{i=1}^{n}\alpha_i-\frac12\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_i\alpha_j\,y_iy_j\,\mathbf{x}_i^\top\mathbf{x}_j \quad\text{s.t.}\quad \alpha_i\ge0,\ \ \sum_i\alpha_iy_i=0\;}\]
Two features of this expression drive everything that follows.
The data enters only through inner products \(\mathbf{x}_i^\top\mathbf{x}_j\). Never through the coordinates themselves. That is what makes the kernel substitution of §6.16 possible.
The dual has \(n\) variables, one per observation, regardless of the number of features \(d\). SVMs therefore handle \(d\gg n\) gracefully and \(n\) large poorly — the reverse of most methods.
16.3 KKT complementary slackness and sparsity
The Karush–Kuhn–Tucker conditions require, at the optimum,
\[\boxed{\;\alpha_i\Big[y_i\big(\mathbf{w}^\top\mathbf{x}_i+b\big)-1\Big]=0\quad\text{for every }i\;}\]
The product is zero, so for each \(i\) at least one factor vanishes:
- If \(y_i(\mathbf{w}^\top\mathbf{x}_i+b)>1\) — the point sits strictly beyond the margin — then the bracket is non-zero, forcing \(\alpha_i=0\).
- If \(\alpha_i>0\) then the bracket must vanish, so \(y_i(\mathbf{w}^\top\mathbf{x}_i+b)=1\) and the point lies exactly on the margin boundary.
Combining with \(\mathbf{w}=\sum_i\alpha_iy_i\mathbf{x}_i\): only the points on the margin contribute to \(\mathbf{w}\). Those are the support vectors, and they give the method its name.
The practical consequence is decisive. Delete every non-support-vector from the training set and refit: the solution is identical. The decision boundary is determined by a handful of points, typically a small fraction of \(n\).
set.seed(71)
n_sp <- 300
Xsp <- matrix(rnorm(n_sp * 2), n_sp, 2)
ysp <- factor(ifelse(Xsp[, 1] + Xsp[, 2] > 0.9, "pos", "neg"))
m_full <- ksvm(Xsp, ysp, kernel = "vanilladot", C = 10, scaled = FALSE)#> Setting default kernel parameterssv_idx <- SVindex(m_full)
c(n = n_sp, support_vectors = length(sv_idx),
fraction = round(length(sv_idx) / n_sp, 4))#> n support_vectors fraction
#> 300.0000 23.0000 0.0767# Refit using ONLY the support vectors
m_sv <- ksvm(Xsp[sv_idx, , drop = FALSE], ysp[sv_idx],
kernel = "vanilladot", C = 10, scaled = FALSE)#> Setting default kernel parametersgrid_sp <- as.matrix(expand.grid(x1 = seq(-3, 3, length.out = 60),
x2 = seq(-3, 3, length.out = 60)))
c(predictions_identical =
mean(predict(m_full, grid_sp) == predict(m_sv, grid_sp)))#> predictions_identical
#> 1Roughly a tenth of the data determines the boundary, and a model fitted on that tenth alone reproduces every prediction.
16.4 The decision function
Substituting the representer identity into \(f(\mathbf{x})=\mathbf{w}^\top\mathbf{x}+b\):
\[\boxed{\;\hat y(\mathbf{x}_\star)=\operatorname{sign}\!\left(\sum_{i\in\mathcal{SV}}\alpha_iy_i\,\mathbf{x}_i^\top\mathbf{x}_\star\ +\ b\right)\;}\]
Note the indices carefully: the sum runs over training support vectors \(\mathbf{x}_i\), each paired with the single test point \(\mathbf{x}_\star\) through an inner product. The bias is recovered from any support vector via \(b=y_k-\sum_i\alpha_iy_i\mathbf{x}_i^\top\mathbf{x}_k\).
17 Soft margins
Real data is rarely separable, and even when it is, the maximum-margin solution can be driven by a single mislabelled point. Slack variables \(\xi_i\ge0\) permit violations at a price:
\[\boxed{\;\min_{\mathbf{w},b,\boldsymbol\xi}\ \ \frac12\lVert\mathbf{w}\rVert^2+C\sum_{i=1}^{n}\xi_i \quad\text{s.t.}\quad y_i\big(\mathbf{w}^\top\mathbf{x}_i+b\big)\ \ge\ 1-\xi_i,\ \ \xi_i\ge0\;}\]
The slack must appear in the constraint, not only in the objective. It is the term \(1-\xi_i\) that lets a point sit inside the margin (\(0<\xi_i\le1\)) or on the wrong side of the hyperplane (\(\xi_i>1\)). Without it, the minimizer would simply set every \(\xi_i=0\) and the formulation would collapse back to the hard margin — infeasible on non-separable data.
Hinge-loss form. At the optimum \(\xi_i=\max\big(0,\,1-y_if(\mathbf{x}_i)\big)\), so the problem is equivalent to unconstrained regularized empirical risk minimization:
\[\min_{\mathbf{w},b}\ \ \underbrace{\sum_{i=1}^{n}\max\big(0,\,1-y_if(\mathbf{x}_i)\big)}_{\text{hinge loss}}\ +\ \underbrace{\frac{1}{2C}\lVert\mathbf{w}\rVert^2}_{L_2\text{ penalty}} .\]
The SVM is ridge-penalized hinge-loss classification. Its dual differs from the hard-margin dual in exactly one respect — a box constraint:
\[\max_{\boldsymbol\alpha}\ \sum_i\alpha_i-\frac12\sum_{i,j}\alpha_i\alpha_jy_iy_j\mathbf{x}_i^\top\mathbf{x}_j \quad\text{s.t.}\quad \boxed{0\le\alpha_i\le C},\ \ \sum_i\alpha_iy_i=0 .\]
The ceiling \(C\) caps how much influence any single point can exert — which is what makes soft-margin SVMs robust to outliers.
What \(C\) controls. Large \(C\) penalizes violations heavily, giving a narrow margin that fits the training data tightly (low bias, high variance). Small \(C\) tolerates violations, giving a wide margin and a smoother boundary (high bias, low variance). \(C\) is the regularization dial, and \(1/(2C)\) is the ridge parameter.
set.seed(73)
n_sm <- 200
Xsm <- matrix(rnorm(n_sm * 2), n_sm, 2)
ysm <- factor(ifelse(Xsm[, 1] + Xsm[, 2] + rnorm(n_sm, sd = 0.9) > 0, "pos", "neg"))
sm_df <- data.frame(x1 = Xsm[, 1], x2 = Xsm[, 2], y = ysm)
gx6 <- seq(-3, 3, length.out = 120); gy6 <- seq(-3, 3, length.out = 120)
gd6 <- expand.grid(x1 = gx6, x2 = gy6)
Cs <- c(0.01, 1, 100)
sm_grid <- bind_rows(lapply(Cs, function(cc) {
m <- ksvm(y ~ ., data = sm_df, kernel = "vanilladot", C = cc, scaled = FALSE)
data.frame(gd6, pred = predict(m, gd6),
panel = sprintf("C = %g (%d support vectors)", cc, nSV(m)))
}))#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parametersggplot(sm_grid, aes(x1, x2)) +
geom_raster(aes(fill = pred), alpha = 0.28) +
geom_point(data = sm_df, aes(x1, x2, shape = y), size = 1.3, color = "grey15") +
scale_fill_manual(values = c(neg = "#3B7DD8", pos = "#D8433B"), guide = "none") +
facet_wrap(~ panel, nrow = 1) + coord_fixed() +
labs(title = "The cost parameter C trades margin width against training fit",
subtitle = "Small C: wide margin, many support vectors, smooth boundary. Large C: narrow margin, tight fit",
x = expression(x[1]), y = expression(x[2]), shape = NULL) +
theme_dspa(10)
18 Kernels
18.1 The trick
The dual and the decision function touch the data only through \(\mathbf{x}_i^\top\mathbf{x}_j\). Suppose we first map into a higher-dimensional space, \(\phi:\mathbb{R}^d\to\mathcal{H}\), and run a linear SVM there. Every occurrence of \(\mathbf{x}_i^\top\mathbf{x}_j\) becomes \(\phi(\mathbf{x}_i)^\top\phi(\mathbf{x}_j)\) — and if that inner product can be evaluated without ever computing \(\phi\), the lift is free.
\[K(\mathbf{x}_i,\mathbf{x}_j):=\big\langle\phi(\mathbf{x}_i),\phi(\mathbf{x}_j)\big\rangle .\]
\[\hat y(\mathbf{x}_\star)=\operatorname{sign}\!\left(\sum_{i\in\mathcal{SV}}\alpha_iy_i\,K(\mathbf{x}_i,\mathbf{x}_\star)+b\right)\]
18.2 Mercer’s condition
Mercer’s theorem. A symmetric continuous function \(K\) corresponds to an inner product in some feature space if and only if it is positive semi-definite: for every finite set \(\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\), the Gram matrix \(\mathbf{K}\) with \(K_{ij}=K(\mathbf{x}_i,\mathbf{x}_j)\) satisfies \(\mathbf{c}^\top\mathbf{K}\mathbf{c}\ge0\) for all \(\mathbf{c}\in\mathbb{R}^n\).
This matters for optimization, not aesthetics: a PSD Gram matrix makes the dual a concave maximization with a unique optimum. An indefinite kernel — the sigmoid kernel for most parameter choices — destroys that guarantee, and solvers may return whatever local point they reach.
18.3 The standard kernels
\[ \begin{aligned} \textbf{Linear: }\quad & K(\mathbf{x},\mathbf{z})=\mathbf{x}^\top\mathbf{z}\\[1mm] \textbf{Polynomial: }\quad & K(\mathbf{x},\mathbf{z})=\big(\mathbf{x}^\top\mathbf{z}+c\big)^{p}\\[1mm] \textbf{Gaussian (RBF): }\quad & K(\mathbf{x},\mathbf{z})=\exp\!\left(-\frac{\lVert\mathbf{x}-\mathbf{z}\rVert^2}{2\sigma^2}\right)=\exp\!\big(-\gamma\lVert\mathbf{x}-\mathbf{z}\rVert^2\big)\\[1mm] \textbf{Laplacian: }\quad & K(\mathbf{x},\mathbf{z})=\exp\!\big(-\gamma\lVert\mathbf{x}-\mathbf{z}\rVert_1\big)\\[1mm] \textbf{Sigmoid: }\quad & K(\mathbf{x},\mathbf{z})=\tanh\!\big(\kappa\,\mathbf{x}^\top\mathbf{z}-\delta\big)\quad\text{(not PSD in general)} \end{aligned} \]
Note the reparameterization \(\gamma=1/(2\sigma^2)\) used by
e1071 and kernlab. Large \(\gamma\) means a narrow
kernel and a wiggly boundary; small \(\gamma\) means a wide kernel and a nearly
linear one.
The RBF kernel’s feature space is infinite-dimensional. Expanding the exponential for \(d=1\), \(\gamma=\tfrac12\):
\[e^{-\frac{(x-z)^2}{2}}=e^{-\frac{x^2}{2}}e^{-\frac{z^2}{2}}\sum_{k=0}^{\infty}\frac{(xz)^k}{k!} =\sum_{k=0}^{\infty}\underbrace{\left(e^{-x^2/2}\frac{x^k}{\sqrt{k!}}\right)}_{\phi_k(x)}\underbrace{\left(e^{-z^2/2}\frac{z^k}{\sqrt{k!}}\right)}_{\phi_k(z)},\]
an inner product between two infinite sequences. The kernel trick computes it in \(O(d)\).
set.seed(79)
Xk <- matrix(rnorm(40), 20, 2)
gram <- function(f) outer(1:20, 1:20, Vectorize(\(i, j) f(Xk[i, ], Xk[j, ])))
min_eig <- function(K) min(eigen((K + t(K)) / 2, symmetric = TRUE, only.values = TRUE)$values)
data.frame(
kernel = c("Linear", "Polynomial (p = 3)", "RBF (gamma = 0.5)", "Sigmoid (kappa = 1)"),
min_eigenvalue = signif(c(
min_eig(gram(\(a, b) sum(a * b))),
min_eig(gram(\(a, b) (sum(a * b) + 1)^3)),
min_eig(gram(\(a, b) exp(-0.5 * sum((a - b)^2)))),
min_eig(gram(\(a, b) tanh(sum(a * b) - 1)))), 4),
valid_Mercer_kernel = c(TRUE, TRUE, TRUE, FALSE))The first three Gram matrices are PSD; the sigmoid’s smallest eigenvalue is negative, so it is not a valid Mercer kernel at these parameters and the dual is no longer concave.
18.4 Lifting, done honestly
The standard illustration takes 2-D data that no line separates and shows it becoming separable after a lift. The lift must be an explicit map, and the separating plane must be fitted in the lifted space — otherwise the figure demonstrates nothing.
set.seed(83)
n_ring <- 260
r_in <- sqrt(runif(n_ring / 2, 0, 1)); th_in <- runif(n_ring / 2, 0, 2 * pi)
r_out <- sqrt(runif(n_ring / 2, 2.2, 4)); th_out <- runif(n_ring / 2, 0, 2 * pi)
ring <- data.frame(
x1 = c(r_in * cos(th_in), r_out * cos(th_out)),
x2 = c(r_in * sin(th_in), r_out * sin(th_out)),
y = factor(rep(c("inner", "outer"), each = n_ring / 2)))
ggplot(ring, aes(x1, x2, color = y, shape = y)) +
geom_point(size = 1.9, alpha = 0.85) +
scale_color_manual(values = c(inner = "#D8433B", outer = "#3B7DD8")) +
coord_fixed() +
labs(title = "No line separates these two classes",
subtitle = "Concentric structure: a linear SVM in 2-D cannot do better than chance",
x = expression(x[1]), y = expression(x[2]), color = NULL, shape = NULL) +
theme_dspa()
#> Setting default kernel parameters#> linear_SVM_in_2D_training_error majority_class
#> 0.3923 0.5000Apply the explicit quadratic map \(\phi(x_1,x_2)=(x_1,\ x_2,\ x_1^2+x_2^2)\) and fit a linear SVM in the lifted space:
ring3 <- ring |> mutate(x3 = x1^2 + x2^2)
lin3d <- ksvm(y ~ x1 + x2 + x3, data = ring3, kernel = "vanilladot",
C = 10, scaled = FALSE)#> Setting default kernel parameters#> linear_SVM_in_3D_training_error support_vectors
#> 0 4# Extract the FITTED plane: w = sum_i alpha_i y_i x_i, and b = -b(model)
Wv <- colSums(coef(lin3d)[[1]] * xmatrix(lin3d)[[1]])
bv <- -b(lin3d)
round(c(w1 = Wv[1], w2 = Wv[2], w3 = Wv[3], b = bv), 4)#> w1.x1 w2.x2 w3.x3 b
#> -0.0271 -0.0531 1.6392 -2.6885# The plane w1*x1 + w2*x2 + w3*x3 + b = 0, solved for x3
gx3 <- seq(min(ring3$x1), max(ring3$x1), length.out = 40)
gy3 <- seq(min(ring3$x2), max(ring3$x2), length.out = 40)
Zpl <- outer(gy3, gx3, \(b2, a1) -(bv + Wv[1] * a1 + Wv[2] * b2) / Wv[3])
plot_ly() |>
add_surface(x = gx3, y = gy3, z = Zpl, opacity = 0.55, showscale = FALSE,
colorscale = list(c(0, "#BBBBBB"), c(1, "#BBBBBB")),
name = "Fitted separating plane") |>
add_trace(data = ring3, x = ~x1, y = ~x2, z = ~x3, type = "scatter3d",
mode = "markers", color = ~y,
colors = c(inner = "#D8433B", outer = "#3B7DD8"),
marker = list(size = 3, opacity = 0.85)) |>
layout(title = sprintf("Lifted by phi(x) = (x1, x2, x1^2 + x2^2): now linearly separable (training error %.3f)",
error(lin3d)),
scene = list(xaxis = list(title = "x1"), yaxis = list(title = "x2"),
zaxis = list(title = "x3 = x1^2 + x2^2",
range = c(0, max(ring3$x3)))))The plane is fitted, not placed — its coefficients come from the SVM’s own \(\alpha_i\) via \(\mathbf{w}=\sum_i\alpha_iy_i\mathbf{x}_i\). In the original 2-D space that plane is a circle, which is exactly the boundary an RBF kernel would find without ever constructing \(\phi\).
mt <- mtcars |>
mutate(engine = factor(ifelse(vs == 0, "V-shaped", "Straight")))
ggplot(mt, aes(wt, hp, color = engine, shape = engine)) +
geom_point(size = 3) +
scale_color_manual(values = c(`V-shaped` = "#D8433B", Straight = "#3B7DD8")) +
labs(title = "Automobile weight vs. horsepower, colored by engine type",
subtitle = "Not cleanly separable by a line in these two features",
x = "Weight (1000 lbs)", y = "Gross horsepower",
color = NULL, shape = NULL) +
theme_dspa()
### OLD: works, but clumsy
# mt_svm <- ksvm(engine ~ wt + hp + qsec, data = mt, kernel = "vanilladot",
# C = 10, scaled = TRUE)
# Wm <- colSums(coef(mt_svm)[[1]] * xmatrix(mt_svm)[[1]])
# bm <- -b(mt_svm)
# sc <- scaling(mt_svm)$x.scale
# sc_c <- sc[["scaled:center"]]
# sc_s <- sc[["scaled:scale"]]
#
# gw <- seq(min(mt$wt), max(mt$wt), length.out = 40)
# gh <- seq(min(mt$hp), max(mt$hp), length.out = 40)
# # Solve the plane for qsec, undoing ksvm's internal scaling
# Zq <- outer(gh, gw, function(h, w) {
# ws <- (w - sc_c["wt"]) / sc_s["wt"]; hs <- (h - sc_c["hp"]) / sc_s["hp"]
# qs <- -(bm + Wm[1] * ws + Wm[2] * hs) / Wm[3]
# qs * sc_s["qsec"] + sc_c["qsec"]
# })
### New approach: compute directly in original space
## using a linear kernel, we can compute the weight vector and intercept directly
## in the original (unscaled) space using the support vectors returned by xmatrix().
## This avoids manual scaling entirely and is less error‑prone.
library(kernlab)
# Fit the model
mt_svm <- ksvm(engine ~ wt + hp + qsec, data = mt,
kernel = "vanilladot", C = 10, scaled = TRUE)#> Setting default kernel parameters# Extract coefficients and support vectors (already in original units)
coefs <- coef(mt_svm)[[1]]
SVs <- xmatrix(mt_svm)[[1]]
# Weight vector for linear kernel in original space
Wm <- colSums(coefs * SVs) # length 3
b0 <- -b(mt_svm) # intercept (decision boundary: Wm'x + b0 = 0)
# Verify
print(Wm)#> wt hp qsec
#> 1.772554 0.736698 -4.592907#> [1] 0.0109947# Create grid for wt and hp
gw <- seq(min(mt$wt), max(mt$wt), length.out = 40)
gh <- seq(min(mt$hp), max(mt$hp), length.out = 40)
# Compute qsec on the separating hyperplane
Zq <- outer(gh, gw, function(h, w) {
-(b0 + Wm[1]*w + Wm[2]*h) / Wm[3]
})
plot_ly() |>
add_surface(x = gw, y = gh, z = Zq, opacity = 0.5, showscale = FALSE,
colorscale = list(c(0, "#BBBBBB"), c(1, "#BBBBBB"))) |>
add_trace(data = mt, x = ~wt, y = ~hp, z = ~qsec, type = "scatter3d",
mode = "markers", color = ~engine,
colors = c(`V-shaped` = "#D8433B", Straight = "#3B7DD8"),
marker = list(size = 5)) |>
layout(title = sprintf("Adding quarter-mile time separates the engine types (training error %.3f)",
error(mt_svm)),
scene = list(xaxis = list(title = "Weight"),
yaxis = list(title = "Horsepower"),
zaxis = list(title = "1/4 mile time",
range = range(mt$qsec))))19 The \(C\times\gamma\) interaction
An RBF SVM has two hyperparameters, and they interact. \(\gamma\) sets the kernel’s reach — how far a single training point’s influence extends — while \(C\) sets how hard the fit is pushed. Large values of either produce overfitting, and their effects partly compensate, so they must be tuned jointly. A sequential search over one and then the other will land in the wrong place.
set.seed(89)
gammas <- c(0.05, 0.5, 10)
gd_rbf <- expand.grid(x1 = seq(-3, 3, length.out = 70),
x2 = seq(-3, 3, length.out = 70))
p <- plot_ly()
for (i in seq_along(gammas)) {
m <- ksvm(y ~ ., data = sm_df, kernel = "rbfdot",
kpar = list(sigma = gammas[i]), C = 1, scaled = FALSE)
dv <- matrix(predict(m, gd_rbf, type = "decision"), 70)
p <- add_surface(p, x = seq(-3, 3, length.out = 70),
y = seq(-3, 3, length.out = 70),
z = pmin(pmax(dv, -3), 3) + (i - 1) * 7,
showscale = FALSE, opacity = 0.95, colorscale = "RdBu",
reversescale = TRUE, name = paste("gamma =", gammas[i]))
}
p |> layout(title = "RBF decision surfaces: gamma = 0.05 (bottom), 0.5 (middle), 10 (top)",
scene = list(xaxis = list(title = "x1"), yaxis = list(title = "x2"),
zaxis = list(title = "Decision value + offset")))At \(\gamma=0.05\) the surface is a near-planar ramp — the kernel is so wide that the model is effectively linear. At \(\gamma=10\) it is a field of isolated bumps, one around each training point: the model has memorized the data and will predict the majority class everywhere else.
The cross-validated error over the two-dimensional grid is itself a surface:
set.seed(97)
C_grid <- 2^seq(-4, 8, by = 2)
g_grid <- 2^seq(-8, 4, by = 2)
cv_err <- outer(seq_along(C_grid), seq_along(g_grid), Vectorize(function(i, j) {
m <- ksvm(y ~ ., data = sm_df, kernel = "rbfdot",
kpar = list(sigma = g_grid[j]), C = C_grid[i],
cross = 5, scaled = FALSE)
cross(m)
}))
best <- which(cv_err == min(cv_err), arr.ind = TRUE)[1, ]
c(best_C = C_grid[best[1]], best_gamma = g_grid[best[2]],
cv_error = round(min(cv_err), 4))#> best_C best_gamma cv_error
#> 16.0000 0.0625 0.2300plot_ly(x = log2(g_grid), y = log2(C_grid), z = cv_err, type = "surface",
colorscale = "Viridis", colorbar = list(title = "5-fold CV error")) |>
add_trace(x = log2(g_grid[best[2]]), y = log2(C_grid[best[1]]),
z = min(cv_err), type = "scatter3d", mode = "markers",
name = "CV optimum",
marker = list(size = 7, color = "red")) |>
layout(title = "Cross-validated error over the (C, gamma) grid",
scene = list(xaxis = list(title = "log2(gamma)"),
yaxis = list(title = "log2(C)"),
zaxis = list(title = "CV error")))The valley runs diagonally, which is the geometric statement of the interaction: an increase in \(\gamma\) can be partly offset by a decrease in \(C\). A one-at-a-time search follows an axis and misses the diagonal.
20 Multi-class SVM
SVMs are intrinsically binary — the margin is defined between two classes. Multi-class problems are handled by reduction:
| Strategy | Machines fitted | Prediction | Cost |
|---|---|---|---|
| One-vs-one | \(\binom{K}{2}\) | Majority vote | Each machine sees only 2 classes, so each is small |
| One-vs-rest | \(K\) | Largest decision value | Each machine sees all \(n\), and classes are imbalanced |
| Crammer–Singer | 1 joint problem | Direct argmax | Single QP, larger and slower |
kernlab::ksvm and e1071::svm both default
to one-vs-one. For the 26-class OCR problem in §6.20 that is \(\binom{26}{2}=325\) binary machines — a
hidden cost worth knowing before you wait for it.
21 Computational complexity
| Quantity | Cost |
|---|---|
| Kernel (Gram) matrix | \(O(n^2 d)\) time, \(\mathbf{O(n^2)}\) memory |
| Training (SMO, typical) | \(O(n^2 d)\) |
| Training (worst case) | \(O(n^3)\) |
| Prediction, one point | \(O(n_{SV}\,d)\) |
| Linear SVM (LIBLINEAR) | \(O(nd)\) — linear in \(n\) |
This is why SVMs stopped scaling. The Gram matrix alone is \(n^2\): at \(n=10^5\) that is 80 GB in double precision, before any optimization. Training cost between \(O(n^2)\) and \(O(n^3)\) means a tenfold increase in data costs 100–1000× the time. Tree ensembles are \(O(Bmn\log n)\) — essentially linear — which is the practical reason they displaced SVMs on large datasets after 2010.
Three escapes, in order of preference. Use a linear SVM when \(d\) is large and the data is nearly separable — LIBLINEAR is \(O(nd)\) and needs no Gram matrix. Approximate the kernel with Nyström or random Fourier features, reducing to a linear problem in \(m\ll n\) dimensions. Or subsample, since the solution depends only on support vectors anyway.
set.seed(101)
time_svm <- function(n, d = 10) {
X <- matrix(rnorm(n * d), n, d)
y <- factor(ifelse(rowSums(X[, 1:3]) + rnorm(n) > 0, "a", "b"))
c(n = n,
gram_memory_MB = round(8 * n^2 / 1e6, 1),
seconds = round(system.time(
ksvm(X, y, kernel = "rbfdot", C = 1, scaled = FALSE))[["elapsed"]], 3))
}
as.data.frame(do.call(rbind, lapply(c(500, 1000, 2000, 4000), time_svm)))Time grows faster than linearly and the Gram-matrix footprint grows quadratically. Extrapolate the last row to \(n=10^5\) and the memory requirement alone rules the method out.
22 Case study: optical character recognition
Chapter 4 embedded handwritten digits; here we classify handwritten letters.
Protocol. Divide a scanned page into a grid, one glyph per cell; match each glyph to a character; reassemble the characters into words.
The UCI letter-recognition data supplies 20,000 pre-gridded examples of the 26 English capitals, rendered in 20 randomly distorted fonts and summarized by 16 numeric shape features.
letters_df <- dspa_read(
"https://umich.instructure.com/files/2837863/download?download_frd=1",
"HandwrittenLetters.csv", header = TRUE)
letters_df$letter <- factor(letters_df$letter)
dim(letters_df); str(letters_df[, 1:6])#> [1] 20000 17#> 'data.frame': 20000 obs. of 6 variables:
#> $ letter: Factor w/ 26 levels "A","B","C","D",..: 20 9 4 14 7 19 2 1 10 13 ...
#> $ xbox : int 2 5 4 7 2 4 4 1 2 11 ...
#> $ ybox : int 8 12 11 11 1 11 2 1 2 15 ...
#> $ width : int 3 3 6 6 3 5 5 3 4 13 ...
#> $ height: int 5 7 8 6 1 8 4 2 4 9 ...
#> $ onpix : int 1 2 6 3 1 3 4 1 2 7 ...c(classes = nlevels(letters_df$letter),
one_vs_one_machines = choose(nlevels(letters_df$letter), 2))#> classes one_vs_one_machines
#> 26 325set.seed(123)
# Stratified RANDOM split: a positional cut risks a systematic train/test
# difference if the file carries any ordering.
split_ocr <- rsample::initial_split(letters_df, prop = 0.75, strata = letter)
ocr_train <- rsample::training(split_ocr)
ocr_test <- rsample::testing(split_ocr)
c(train = nrow(ocr_train), test = nrow(ocr_test),
no_information_rate = round(max(prop.table(table(ocr_test$letter))), 4))#> train test no_information_rate
#> 1.5e+04 5.0e+03 4.6e-02ksvm centres and scales internally with
scaled = TRUE, learning the constants from the
training data and applying them to any new data — the
one place in a typical pipeline where the discipline of Chapter 5, §5.4.1 is handled for
you.
set.seed(123)
t_lin <- system.time(
ocr_linear <- ksvm(letter ~ ., data = ocr_train, kernel = "vanilladot"))#> Setting default kernel parameters#> Support Vector Machine object of class "ksvm"
#>
#> SV type: C-svc (classification)
#> parameter : cost C = 1
#>
#> Linear (vanilla) kernel function.
#>
#> Number of Support Vectors : 6678
#>
#> Objective Function Value : -12.0757 -22.0282 -26.0057 -6.906 -6.7466 -36.1311 -49.578 -18.4749 -52.6397 -33.831 -18.8742 -31.9876 -30.3967 -50.4276 -4.9888 -37.9956 -28.4969 -17.4952 -14.8849 -37.1153 -27.0704 -7.1844 -11.194 -29.8494 -12.6595 -8.5352 -142.457 -46.8674 -54.3687 -117.056 -144.304 -59.496 -45.4833 -63.4898 -22.2419 -24.6064 -19.6607 -33.8424 -34.7173 -122.151 -184.208 -199.233 -21.2684 -10.3033 -54.617 -10.7174 -48.1005 -9.2705 -18.6785 -11.0539 -109.554 -29.2836 -212.733 -67.4051 -7.1255 -4.6424 -131.291 -82.8987 -21.1632 -15.7521 -75.2044 -11.3993 -28.4599 -18.9696 -18.3726 -25.2114 -49.7602 -10.0309 -4.6463 -12.6668 -4.8506 -3.4229 -7.6366 -33.8567 -53.2168 -160.27 -45.1946 -50.335 -44.0366 -16.9098 -16.6457 -84.1637 -111.178 -34.6761 -34.1307 -108.474 -31.6224 -26.008 -29.5989 -16.1122 -3.9382 -39.6862 -9.4329 -18.2096 -41.4656 -154.896 -43.3248 -36.2807 -32.2206 -68.7536 -122.723 -10.1966 -5.5764 -13.0949 -21.4152 -123.464 -48.9412 -161.236 -93.1671 -9.6963 -17.0088 -3.8622 -67.5234 -7.587 -87.7238 -43.0385 -74.3929 -72.9996 -66.3438 -16.1875 -9.9742 -7.154 -25.6414 -12.2905 -206.052 -27.871 -19.4275 -119.507 -116.111 -8.2955 -31.8314 -6.0869 -44.8874 -67.5061 -22.782 -186.238 -33.3417 -17.4044 -121.957 -166.745 -46.1363 -22.1868 -149.378 -67.0377 -327.417 -147.632 -145.593 -29.1456 -34.7619 -52.8868 -27.7469 -39.5101 -7.5607 -10.4002 -33.074 -58.3279 -171.347 -54.2702 -89.269 -137.173 -563.328 -95.5255 -122.559 -287.761 -34.7784 -59.4846 -144.421 -103.476 -30.4757 -57.2788 -48.43 -7.2501 -211.407 -14.4372 -39.5369 -2.2074 -6.5635 -16.6399 -23.0932 -49.5148 -21.4599 -178.968 -19.866 -4.5403 -4.5757 -0.7949 -115.518 -7.299 -60.4205 -20.3351 -14.7527 -4.441 -13.4378 -29.2422 -19.463 -84.8545 -25.1506 -100.317 -14.9194 -10.4036 -7.0877 -1.2784 -74.8735 -7.4496 -98.4933 -100.85 -40.1475 -24.2442 -52.5388 -19.292 -48.8228 -236.591 -39.1974 -34.8078 -29.5919 -16.4216 -11.7809 -108.322 -6.2962 -5.3582 -9.1017 -15.091 -26.9628 -17.8214 -130.216 -32.7434 -97.4543 -34.315 -14.5081 -9.5127 -3.1255 -86.1746 -7.5975 -14.9142 -63.8149 -103.359 -10.5868 -13.6621 -53.2336 -2.8738 -7.8404 -77.6159 -37.9807 -105.471 -3.1303 -6.8417 -1.2321 -92.0831 -18.2633 -9.5768 -46.6064 -3.4108 -19.4482 -66.3247 -42.3009 -48.206 -4.964 -17.5083 -2.3064 -68.4819 -116.3 -107.984 -26.6367 -21.2674 -52.5151 -38.445 -63.8797 -18.8368 -6.3928 -4.7918 -47.9023 -26.4668 -48.0111 -22.8189 -5.2968 -46.3513 -12.3037 -16.411 -57.1823 -3.079 -59.3924 -235.702 -14.5906 -11.8576 -17.474 -9.2204 -54.4917 -15.8317 -38.5418 -50.5338 -25.6608 -15.1738 -43.5911 -14.5403 -57.9278 -5.9907 -5.8003 -79.5683 -3.2532 -6.5919 -1.0693 -124.544 -25.2817 -353.582 -28.8275 -30.8154 -4.8432 -71.3641 -128.012 -72.7254 -28.1138 -36.5386 -10.3768 -21.5268 -1.9431 -61.6687 -6.9603 -134.924 -1.7064 -1.303 -8.9234 -0.5131 -25.9642 -30.5132 -5.8562
#> Training error : 0.1334#> seconds
#> 2.1pred_lin <- predict(ocr_linear, ocr_test)
acc_lin <- mean(pred_lin == ocr_test$letter)
c(accuracy = round(acc_lin, 4),
kappa = round(confusionMatrix(pred_lin, ocr_test$letter)$overall[["Kappa"]], 4))#> accuracy kappa
#> 0.8570 0.8513set.seed(123)
t_rbf <- system.time(
ocr_rbf <- ksvm(letter ~ ., data = ocr_train, kernel = "rbfdot"))
pred_rbf <- predict(ocr_rbf, ocr_test)
data.frame(
kernel = c("Majority class", "Linear (vanilladot)", "RBF (rbfdot)"),
accuracy = round(c(max(prop.table(table(ocr_test$letter))),
acc_lin, mean(pred_rbf == ocr_test$letter)), 4),
support_vectors = c(NA, nSV(ocr_linear), nSV(ocr_rbf)),
fit_seconds = round(c(NA, t_lin[["elapsed"]], t_rbf[["elapsed"]]), 1))The RBF kernel gains a large margin over the linear one, which says the class boundaries in these 16 shape features are genuinely curved. Note the support vector counts — a substantial fraction of the training set — and that the RBF fit takes noticeably longer, both consequences of §6.19.
cm_ocr <- table(Predicted = pred_rbf, Actual = ocr_test$letter)
as.data.frame(cm_ocr) |>
ggplot(aes(Actual, Predicted, fill = log1p(Freq))) +
geom_tile() +
scale_fill_viridis_c(option = "magma", direction = -1, name = "log(1+count)") +
coord_fixed() +
labs(title = "OCR confusion matrix (RBF kernel, held-out data)",
subtitle = "Off-diagonal mass concentrates on visually similar letters",
x = "True letter", y = "Predicted letter") +
theme_dspa(8)
# --- Interactive equivalent ------------------------------------------------
# Note the axes: rownames are PREDICTIONS, colnames are TRUTH.
plot_ly(x = colnames(cm_ocr), y = rownames(cm_ocr),
z = matrix(as.numeric(cm_ocr), nrow(cm_ocr)), type = "heatmap") |>
layout(title = "OCR confusion matrix",
xaxis = list(title = "True letter"),
yaxis = list(title = "Predicted letter"))# Which letters are hardest, and what are they confused with?
off <- as.data.frame(cm_ocr) |> filter(Predicted != Actual) |> arrange(desc(Freq))
head(off, 8)per_letter <- as.data.frame(cm_ocr) |>
summarise(recall = Freq[Predicted == Actual] / sum(Freq), .by = Actual) |>
arrange(recall)
head(per_letter, 6) |> mutate(recall = round(recall, 4))The confusions are the ones a human would make — visually similar glyph pairs. That is reassuring: the model has learned shape, not an artifact.
23 Case study: iris flowers
#>
#> setosa versicolor virginica
#> 50 50 50set.seed(1234)
split_ir <- rsample::initial_split(iris, prop = 0.75, strata = Species)
iris_train <- rsample::training(split_ir)
iris_test <- rsample::testing(split_ir)# Fit on TRAINING data; draw the boundary over a grid; overlay training points.
gp <- expand.grid(
Petal.Length = seq(min(iris$Petal.Length), max(iris$Petal.Length), length.out = 220),
Petal.Width = seq(min(iris$Petal.Width), max(iris$Petal.Width), length.out = 220))
bd <- bind_rows(lapply(c("linear", "radial"), function(k) {
m <- svm(Species ~ Petal.Length + Petal.Width, data = iris_train,
kernel = k, cost = 1)
data.frame(gp, pred = predict(m, gp), kernel = paste(k, "kernel"))
}))
ggplot(bd, aes(Petal.Length, Petal.Width)) +
geom_raster(aes(fill = pred), alpha = 0.30) +
geom_point(data = iris_train, aes(Petal.Length, Petal.Width, shape = Species),
size = 1.9, color = "grey15") +
scale_fill_brewer(palette = "Set2", name = "Region") +
facet_wrap(~ kernel) + coord_fixed() +
labs(title = "Linear and RBF decision regions on two iris features",
subtitle = "Fitted on the training split; regions drawn over the full feature range",
x = "Petal length", y = "Petal width", shape = NULL) +
theme_dspa(10)
#> Setting default kernel parametersiris_rbf <- ksvm(Species ~ ., data = iris_train, kernel = "rbfdot")
data.frame(
model = c("Majority class", "Linear kernel", "RBF kernel"),
test_accuracy = round(c(
max(prop.table(table(iris_test$Species))),
mean(predict(iris_lin, iris_test) == iris_test$Species),
mean(predict(iris_rbf, iris_test) == iris_test$Species)), 4),
support_vectors = c(NA, nSV(iris_lin), nSV(iris_rbf)))#> Reference
#> Prediction setosa versicolor virginica
#> setosa 13 0 0
#> versicolor 0 13 2
#> virginica 0 0 11The RBF kernel does not improve on the linear one here, because the iris classes are close to linearly separable in these four features. A more flexible model is not automatically a better one — flexibility costs variance and must be paid for by structure the data actually contains.
23.1 Tuning, without transcribing
set.seed(2020)
costs <- 2^(-4:8)
gammas <- 2^(-8:2)
tuned <- tune.svm(Species ~ ., data = iris_train, kernel = "radial",
cost = costs, gamma = gammas,
tunecontrol = tune.control(sampling = "cross", cross = 10))
tuned$best.parameters#> [1] 0.0273# Read the tuned values from the object -- never transcribe them by hand,
# or they go stale the moment the data or the seed changes.
iris_tuned <- svm(Species ~ ., data = iris_train, kernel = "radial",
cost = tuned$best.parameters$cost,
gamma = tuned$best.parameters$gamma)
c(test_accuracy = round(mean(predict(iris_tuned, iris_test) == iris_test$Species), 4))#> test_accuracy
#> 0.9487# Cross-validation INSIDE the training data only; the test set is never seen.
set.seed(2017)
folds_ir <- rsample::vfold_cv(iris_train, v = 5, strata = Species)
cost_curve <- do.call(rbind, lapply(costs, function(cc) {
cv <- mean(vapply(folds_ir$splits, function(s) {
a <- rsample::analysis(s); b <- rsample::assessment(s)
m <- svm(Species ~ ., data = a, kernel = "radial", cost = cc)
mean(predict(m, b) != b$Species)
}, numeric(1)))
m_full <- svm(Species ~ ., data = iris_train, kernel = "radial", cost = cc)
data.frame(cost = cc,
Train = mean(m_full$fitted != iris_train$Species),
CV = cv,
Test = mean(predict(m_full, iris_test) != iris_test$Species))
}))
cost_curve |>
pivot_longer(-cost, names_to = "set", values_to = "error") |>
ggplot(aes(cost, error, color = set)) +
geom_line(linewidth = 1) + geom_point(size = 2) +
scale_x_log10() +
scale_color_manual(values = c(Train = "#7FB069", CV = "#D8433B",
Test = "#3B7DD8")) +
labs(title = "SVM error against the cost parameter",
subtitle = "CV computed inside the training split; Test is held out and used only for the final report",
x = "Cost C (log scale)", y = "Classification error", color = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
dl <- cost_curve |> pivot_longer(-cost, names_to = "variable", values_to = "value")
plot_ly(dl, x = ~log(cost), y = ~value, color = ~variable,
type = "scatter", mode = "lines+markers") |>
layout(title = "SVM CV plot of model performance (iris data)",
xaxis = list(title = "log(Cost)"),
yaxis = list(title = "Classifier error"),
legend = list(orientation = "h"))Training error falls monotonically toward zero as \(C\) grows — it always will, and it is therefore useless for selection. The CV curve has an interior minimum and is the one to read.
24 PART III — ENSEMBLES
Ensembles combine many models into one. The two classical families attack different terms of the error decomposition, and confusing them is the most common conceptual error in this area.
25 Why ensembling works
For squared loss, the expected error of a predictor decomposes as
\[\mathbb{E}\big[(y-\hat f(x))^2\big]=\underbrace{\big(\mathbb{E}\hat f(x)-f(x)\big)^2}_{\text{bias}^2}+\underbrace{\operatorname{Var}\big(\hat f(x)\big)}_{\text{variance}}+\underbrace{\sigma^2}_{\text{irreducible}} .\]
Bagging attacks the variance term. Boosting attacks the bias term.
Averaging \(B\) predictors each with variance \(\sigma^2\) and pairwise correlation \(\rho\) gives (Chapter 5, §5.22.1)
\[\boxed{\;\operatorname{Var}\!\left(\frac1B\sum_b\hat f_b\right)=\rho\sigma^2+\frac{1-\rho}{B}\,\sigma^2\;\xrightarrow[B\to\infty]{}\;\rho\sigma^2\;}\]
Note that this identity follows from bilinearity of variance, not from the central limit theorem — the CLT concerns the asymptotic normality of a standardized mean, which is a different statement.
The second term vanishes with \(B\); the first does not. The floor is set by how correlated the members are, which is why decorrelating them matters more than adding more of them.
Bagged predictors are identically distributed, so \(\mathbb{E}[\bar f]=\mathbb{E}[\hat f_b]\) — bagging does not change the bias. Only variance falls. Boosting, by fitting each new member to what the current ensemble gets wrong, reduces bias instead.
26 Bagging
Bootstrap aggregating: draw \(B\) bootstrap resamples of size \(n\), fit a model to each, and average (regression) or vote (classification).
Each bootstrap sample omits a fraction \(\left(1-\frac1n\right)^n\to e^{-1}\approx0.368\) of the data. Those out-of-bag observations give a free validation estimate: predict each case using only the trees that did not see it.
Bagging helps most for unstable learners — those whose fit changes a lot under resampling. Deep decision trees are the archetype; linear regression is nearly unaffected, because it is already stable.
set.seed(107)
n_b <- 150
x_b <- sort(runif(n_b, 0, 1))
y_b <- f_target(x_b) + rnorm(n_b, sd = 0.15)
xs_b <- seq(0, 1, length.out = 300)
resample_fits <- function(fitter, reps = 25) {
sapply(seq_len(reps), function(r) {
idx <- sample(n_b, replace = TRUE)
fitter(x_b[idx], y_b[idx], xs_b)
})
}
tree_fit <- function(x, y, xs)
predict(rpart::rpart(y ~ x, data.frame(x, y),
control = rpart::rpart.control(cp = 0, minsplit = 5)),
data.frame(x = xs))
lm_fit <- function(x, y, xs)
predict(lm(y ~ poly(x, 3)), data.frame(x = xs))
set.seed(109); Ptree <- resample_fits(tree_fit)
set.seed(109); Plm <- resample_fits(lm_fit)
bind_rows(
data.frame(x = rep(xs_b, 25), y = as.vector(Ptree),
rep = rep(1:25, each = 300), learner = "Deep tree (unstable)"),
data.frame(x = rep(xs_b, 25), y = as.vector(Plm),
rep = rep(1:25, each = 300), learner = "Cubic polynomial (stable)")) |>
ggplot(aes(x, y, group = rep)) +
geom_line(alpha = 0.25, color = "steelblue", linewidth = 0.4) +
geom_line(data = data.frame(x = xs_b, y = f_target(xs_b)), aes(x, y),
inherit.aes = FALSE, color = "firebrick", linewidth = 1) +
facet_wrap(~ learner) +
coord_cartesian(ylim = c(-1, 1.2)) +
labs(title = "Bagging helps unstable learners",
subtitle = "25 bootstrap fits each. Red: the truth. Wide spread means high variance -- and room for averaging to help",
x = "x", y = "f(x)") +
theme_dspa(10)
c(tree_variance = round(mean(apply(Ptree, 1, var)), 5),
lm_variance = round(mean(apply(Plm, 1, var)), 5))#> tree_variance lm_variance
#> 0.01403 0.0021427 Random forests
Common misconception: “random forests are a boosting method.” Random forests are bagging, with one addition. Boosting fits members sequentially, each correcting its predecessors’ errors, and reduces bias. Bagging fits members independently and in parallel, and reduces variance. A random forest’s trees never see each other’s residuals.
The addition is feature subsampling: at every split, only a random subset of \(m\) predictors is considered. Its purpose is to lower the correlation \(\rho\) in the variance identity, since without it a single dominant predictor would occupy the top split of nearly every tree and \(\rho\) would stay high. Typical defaults are \(m=\lfloor\sqrt p\rfloor\) for classification and \(p/3\) for regression.
Algorithm. For \(b=1,\dots,B\): draw a bootstrap sample \(Z^*_b\) of size \(n\); grow a tree \(T_b\) by recursively selecting \(m\ll p\) features at random, choosing the best split among them, and splitting, until the minimum node size is reached. Then
\[\hat f^{RF}_B(\mathbf{x})=\frac1B\sum_{b=1}^{B}T_b(\mathbf{x})\ \ \text{(regression)}, \qquad \hat C^{RF}_B(\mathbf{x})=\text{majority vote}\big\{\hat C_b(\mathbf{x})\big\}\ \ \text{(classification)}.\]
Because the trees are independent, the whole procedure is embarrassingly parallel — a property boosting cannot have.
28 Boosting
Boosting builds an additive model sequentially:
\[F_M(\mathbf{x})=\sum_{m=1}^{M}\nu\,h_m(\mathbf{x}),\]
where each \(h_m\) is a weak learner (commonly a shallow tree) and \(\nu\in(0,1]\) is a shrinkage or learning rate. The two classical algorithms differ in what each new learner is fitted to, and that difference traces back to the loss.
28.1 AdaBoost, from exponential loss
AdaBoost minimizes exponential loss \(\mathcal{L}(y,F)=e^{-yF(\mathbf{x})}\) over \(y\in\{-1,+1\}\), by forward stagewise additive modelling. At stage \(m\) with the current ensemble \(F_{m-1}\), solve
\[(\beta_m,h_m)=\arg\min_{\beta,h}\ \sum_{i=1}^{n}\exp\!\big(-y_i\big[F_{m-1}(\mathbf{x}_i)+\beta h(\mathbf{x}_i)\big]\big) =\arg\min_{\beta,h}\ \sum_{i=1}^{n}w_i^{(m)}e^{-\beta y_i h(\mathbf{x}_i)},\]
where \(w_i^{(m)}=e^{-y_iF_{m-1}(\mathbf{x}_i)}\) is the observation weight — it appears automatically, not by fiat. Cases the ensemble currently misclassifies have \(y_iF_{m-1}<0\) and therefore large weight.
For \(h\in\{-1,+1\}\), split the sum by whether \(h\) is right or wrong:
\[\sum_i w_i^{(m)}e^{-\beta y_ih(\mathbf{x}_i)}=e^{-\beta}\!\!\sum_{y_i=h(\mathbf{x}_i)}\!\!w_i^{(m)}+e^{\beta}\!\!\sum_{y_i\ne h(\mathbf{x}_i)}\!\!w_i^{(m)} =\big(e^{\beta}-e^{-\beta}\big)\sum_i w_i^{(m)}\mathbb{1}\{y_i\ne h_i\}+e^{-\beta}\sum_i w_i^{(m)} .\]
The first factor is positive for \(\beta>0\), so the optimal \(h_m\) is the one minimizing the weighted error rate \(\epsilon_m=\sum_iw_i^{(m)}\mathbb{1}\{y_i\ne h_i\}/\sum_iw_i^{(m)}\). Differentiating with respect to \(\beta\) and setting to zero:
\[\boxed{\;\alpha_m=\frac12\ln\!\left(\frac{1-\epsilon_m}{\epsilon_m}\right)\;}\]
This is the exact line search for exponential loss, not a heuristic. Weights then update multiplicatively, \(w_i^{(m+1)}=w_i^{(m)}e^{-\alpha_my_ih_m(\mathbf{x}_i)}\), and are renormalized.
The variants differ in how \(\alpha\) is defined for multi-class problems:
| Variant | Coefficient | Requires |
|---|---|---|
| Breiman (AdaBoost.M1, default) | \(\alpha=\tfrac12\ln\frac{1-\epsilon}{\epsilon}\) | \(\epsilon<\tfrac12\) |
| Freund | \(\alpha=\ln\frac{1-\epsilon}{\epsilon}\) | \(\epsilon<\tfrac12\) |
| Zhu (SAMME) | \(\alpha=\ln\frac{1-\epsilon}{\epsilon}+\ln(K-1)\) | \(\epsilon<1-\tfrac1K\) |
The \(\ln(K-1)\) term in SAMME is what makes multi-class boosting work: with \(K\) classes, a learner beating chance (\(\epsilon<1-1/K\)) should still contribute, whereas AdaBoost.M1 demands the far stricter \(\epsilon<1/2\) and stalls when no weak learner meets it.
eps <- seq(0.01, 0.99, length.out = 400)
data.frame(
epsilon = rep(eps, 3),
alpha = c(0.5 * log((1 - eps) / eps),
log((1 - eps) / eps),
log((1 - eps) / eps) + log(3 - 1)),
variant = rep(c("Breiman (M1)", "Freund", "Zhu (SAMME, K = 3)"), each = length(eps))) |>
ggplot(aes(epsilon, alpha, color = variant)) +
geom_hline(yintercept = 0, color = "grey60") +
geom_vline(xintercept = c(0.5, 2/3), linetype = "dashed", color = "grey45") +
geom_line(linewidth = 1) +
coord_cartesian(ylim = c(-3, 4)) +
scale_color_brewer(palette = "Set1") +
labs(title = "The boosting coefficient as a function of weighted error",
subtitle = "Dashed lines: eps = 1/2 and eps = 1 - 1/K. Alpha crosses zero where a learner stops helping",
x = expression(epsilon[m]), y = expression(alpha[m]), color = NULL) +
theme_dspa()
Read where each curve crosses zero. Breiman and Freund stop rewarding a learner at \(\epsilon=1/2\); SAMME keeps rewarding it until \(\epsilon=1-1/K=2/3\) for three classes, which is exactly the range where a multi-class weak learner lives.
28.2 Gradient boosting, as functional gradient descent
Gradient boosting generalizes to any differentiable loss. View the ensemble \(F\) as a point in function space and take a gradient step: fit the next learner to the negative gradient of the loss with respect to the current predictions,
\[r_i^{(m)}=-\left.\frac{\partial\mathcal{L}\big(y_i,F(\mathbf{x}_i)\big)}{\partial F(\mathbf{x}_i)}\right|_{F=F_{m-1}}, \qquad h_m=\arg\min_h\sum_i\big(r_i^{(m)}-h(\mathbf{x}_i)\big)^2,\]
then update \(F_m=F_{m-1}+\nu\,\gamma_m h_m\) with \(\gamma_m\) from a line search.
| Loss | Negative gradient \(r_i\) | Yields |
|---|---|---|
| Squared, \(\tfrac12(y-F)^2\) | \(y_i-F(\mathbf{x}_i)\) — the residual | \(L_2\) boosting |
| Absolute, \(\lvert y-F\rvert\) | \(\operatorname{sign}(y_i-F(\mathbf{x}_i))\) | Robust regression |
| Exponential, \(e^{-yF}\) | \(y_ie^{-y_iF(\mathbf{x}_i)}\) | AdaBoost |
| Logistic, \(\ln(1+e^{-yF})\) | \(y_i/(1+e^{y_iF(\mathbf{x}_i)})\) | LogitBoost |
AdaBoost and gradient boosting are not the same algorithm. AdaBoost reweights observations and is specific to exponential loss; gradient boosting fits each learner to a pseudo-residual and works for any differentiable loss. They coincide only in the sense that AdaBoost is gradient boosting under exponential loss.
The residual recursion \(r=y-F_{m-1}\) that people often write for “boosting” is the squared-loss case specifically — and note the offset is the cumulative ensemble \(F_{m-1}=\sum_{j<m}\nu\gamma_jh_j\), not just the previous learner alone.
Common misconception: “boosting is guaranteed to beat the best individual learner.” There is no such guarantee, and boosting demonstrably overfits given enough rounds. The Schapire–Freund result is that weak learnability implies strong learnability: a learner beating chance by a fixed margin can be boosted to arbitrarily low training error, with generalization controlled by margin-based bounds. That is a PAC-learnability statement, not a promise about held-out performance. Boosting needs \(M\) tuned by validation exactly as any other capacity parameter does.
Modern implementations. xgboost,
lightgbm, and catboost add second-order
(Newton) steps using the Hessian, explicit \(L_1\)/\(L_2\) penalties on leaf weights, column and
row subsampling, histogram-based split finding, and native missing-value
handling. They are the current standard for tabular prediction, and
their cost is \(O(M\cdot m\cdot n\log
n)\) — sequential in \(M\),
which is why forests parallelize and boosters do not.
set.seed(113)
n_bo <- 400
Xbo <- matrix(rnorm(n_bo * 5), n_bo, 5)
ybo <- as.numeric(Xbo[, 1] + Xbo[, 2]^2 + rnorm(n_bo, sd = 1.2))
bo_df <- data.frame(Xbo, y = ybo)
tr_bo <- 1:250; te_bo <- 251:n_bo
library(gbm)
set.seed(113)
gb <- gbm(y ~ ., data = bo_df[tr_bo, ], distribution = "gaussian",
n.trees = 1200, interaction.depth = 3, shrinkage = 0.05,
bag.fraction = 1, train.fraction = 1, verbose = FALSE)
Ms <- seq(10, 1200, by = 25)
curve_bo <- do.call(rbind, lapply(Ms, function(M) {
data.frame(M = M,
Train = sqrt(mean((predict(gb, bo_df[tr_bo, ], n.trees = M) - ybo[tr_bo])^2)),
Test = sqrt(mean((predict(gb, bo_df[te_bo, ], n.trees = M) - ybo[te_bo])^2)))
}))
curve_bo |> pivot_longer(-M, names_to = "set", values_to = "rmse") |>
ggplot(aes(M, rmse, color = set)) +
geom_line(linewidth = 1) +
geom_vline(xintercept = curve_bo$M[which.min(curve_bo$Test)],
linetype = "dashed", color = "grey40") +
scale_color_manual(values = c(Train = "#7FB069", Test = "#3B7DD8")) +
labs(title = "Boosting overfits: training error falls forever, test error turns around",
subtitle = sprintf("Test minimum at M = %d. The number of rounds is a capacity parameter and must be tuned",
curve_bo$M[which.min(curve_bo$Test)]),
x = "Boosting rounds M", y = "RMSE", color = NULL) +
theme_dspa()
29 Stacking
Bagging and boosting combine models of the same type. Stacking combines models of different types by training a meta-learner on their predictions.
The essential detail is that the meta-learner must be trained on out-of-fold predictions. Using in-sample predictions lets a base learner that memorized the training data dominate the meta-learner — a leakage that looks like excellent stacking performance and collapses on new data.
set.seed(127)
stack_data <- data.frame(sm_df)
folds_st <- rsample::vfold_cv(stack_data, v = 5, strata = y)
base_learners <- list(
svm_rbf = function(tr, te) predict(ksvm(y ~ ., data = tr, kernel = "rbfdot",
C = 1, prob.model = TRUE), te,
type = "probabilities")[, "pos"],
rf = function(tr, te) predict(ranger(y ~ ., data = tr, num.trees = 300,
probability = TRUE, num.threads = 1),
te)$predictions[, "pos"],
logistic = function(tr, te) predict(glm(y ~ ., data = tr, family = binomial()),
te, type = "response"))
# Out-of-fold predictions: each row is predicted by a model that never saw it
oof <- matrix(NA_real_, nrow(stack_data), length(base_learners),
dimnames = list(NULL, names(base_learners)))
for (s in folds_st$splits) {
a <- rsample::analysis(s); idx <- setdiff(seq_len(nrow(stack_data)),
as.integer(rownames(a)))
b <- stack_data[idx, , drop = FALSE]
for (nm in names(base_learners)) oof[idx, nm] <- base_learners[[nm]](a, b)
}
meta_df <- data.frame(oof, y = stack_data$y)
meta <- glm(y ~ ., data = meta_df, family = binomial())
round(coef(meta), 4)#> (Intercept) svm_rbf rf logistic
#> -7.9025 -17.4413 25.3910 9.2609set.seed(131)
sp_st <- rsample::initial_split(stack_data, prop = 0.7, strata = y)
tr_st <- rsample::training(sp_st); te_st <- rsample::testing(sp_st)
base_te <- sapply(base_learners, \(f) f(tr_st, te_st))
stack_p <- predict(meta, as.data.frame(base_te), type = "response")
acc <- function(p) mean(factor(ifelse(p > 0.5, "pos", "neg"),
levels = levels(te_st$y)) == te_st$y)
data.frame(
model = c(names(base_learners), "Stacked ensemble"),
test_accuracy = round(c(apply(base_te, 2, acc), acc(stack_p)), 4))The meta-learner’s coefficients show how much weight each base model earns. Stacking helps most when the base learners make different mistakes — which is the same decorrelation principle as \(\rho\) in the variance identity, applied across model families rather than within one.
30 Case study: quality of life and chronic disease
qol_raw <- dspa_read("https://umich.instructure.com/files/481332/download?download_frd=1",
"Case06_QoL_Symptom_ChronicIllness.csv")
dim(qol_raw)#> [1] 2356 41c(missing_charlson = sum(qol_raw$CHARLSONSCORE == -9),
missing_cds = sum(qol_raw$CHRONICDISEASESCORE == -9))#> missing_charlson missing_cds
#> 28 142# Filter ONCE, select columns by NAME. Repeating a positional drop removes
# real features on the second pass.
qol <- qol_raw |>
filter(CHARLSONSCORE != -9, CHRONICDISEASESCORE != -9) |>
dplyr::select(-any_of(c("ID", "INTERVIEWDATE")))
table(qol$CHARLSONSCORE)#>
#> 0 1 2 3 4 5 6 7 8 9 10
#> 799 942 263 123 24 7 5 2 15 8 2The 11 Charlson levels are severely unbalanced — several have fewer than ten cases, which makes stratified resampling fragile and per-class metrics unestimable. Collapsing the sparse tail is a modelling decision, so it is made explicitly.
qol$charlson <- if (CHARLSON_COLLAPSE == "collapse") {
factor(pmin(qol$CHARLSONSCORE, 4),
levels = 0:4, labels = c("0", "1", "2", "3", "4+"))
} else {
factor(qol$CHARLSONSCORE)
}
table(qol$charlson)#>
#> 0 1 2 3 4+
#> 799 942 263 123 63c(smallest_class = min(table(qol$charlson)),
no_information_rate = round(max(prop.table(table(qol$charlson))), 4))#> smallest_class no_information_rate
#> 63.0000 0.4301# Threshold computed from the data (median gives balanced classes, so accuracy
# is interpretable against a ~50% baseline)
cut_cds <- median(qol$CHRONICDISEASESCORE)
qol$cd <- factor(qol$CHRONICDISEASESCORE > cut_cds,
levels = c(FALSE, TRUE),
labels = c("minor_disease", "severe_disease"))
c(cut_point = round(cut_cds, 4))#> cut_point
#> 1.39#>
#> minor_disease severe_disease
#> 0.5009 0.4991model_df <- qol |> dplyr::select(-any_of(c("CHARLSONSCORE", "CHRONICDISEASESCORE", "charlson")))
set.seed(1234)
sp_qol <- rsample::initial_split(model_df, prop = 0.8, strata = cd)
qol_train <- rsample::training(sp_qol)
qol_test <- rsample::testing(sp_qol)
NIR_q <- max(prop.table(table(qol_test$cd)))
c(train = nrow(qol_train), test = nrow(qol_test),
no_information_rate = round(NIR_q, 4))#> train test no_information_rate
#> 1751.0000 439.0000 0.501130.1 Bagging, forests, and boosting
library(ipred)
set.seed(123)
bag_fit <- ipred::bagging(cd ~ ., data = qol_train, nbagg = 50, coob = TRUE)
c(resubstitution_accuracy = round(mean(predict(bag_fit, qol_train) == qol_train$cd), 4),
out_of_bag_error = round(bag_fit$err, 4),
test_accuracy = round(mean(predict(bag_fit, qol_test) == qol_test$cd), 4))#> resubstitution_accuracy out_of_bag_error test_accuracy
#> 1.0000 0.4637 0.5490The three numbers say different things. Resubstitution accuracy is near 1 and means nothing — a 50-tree bag of unpruned trees memorizes its training data by construction. The out-of-bag error is a genuine internal estimate. The test accuracy is the one to report.
set.seed(123)
rf_fit <- ranger(cd ~ ., data = qol_train, num.trees = 500,
importance = "permutation", num.threads = 1)
rf_fit#> Ranger result
#>
#> Call:
#> ranger(cd ~ ., data = qol_train, num.trees = 500, importance = "permutation", num.threads = 1)
#>
#> Type: Classification
#> Number of trees: 500
#> Sample size: 1751
#> Number of independent variables: 37
#> Mtry: 6
#> Target node size: 1
#> Variable importance mode: permutation
#> Splitrule: gini
#> OOB prediction error: 45.00 %rf_pred <- predict(rf_fit, qol_test)$predictions
cm_rf <- confusionMatrix(rf_pred, qol_test$cd, positive = "severe_disease")
cm_rf$overall[1:6] |> round(4)#> Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
#> 0.6059 0.2120 0.5585 0.6519 0.5011
#> AccuracyPValue
#> 0.0000mtry is the only real tuning parameter, and it is
exactly the \(\rho\)-lowering knob of
§6.21.
p_q <- ncol(qol_train) - 1
mtries <- unique(round(c(2, sqrt(p_q), p_q / 3, p_q / 2, p_q)))
oob_by_mtry <- vapply(mtries, function(m) {
set.seed(123)
ranger(cd ~ ., data = qol_train, num.trees = 400, mtry = m,
num.threads = 1)$prediction.error
}, numeric(1))
data.frame(mtry = mtries, oob_error = round(oob_by_mtry, 4))ggplot(data.frame(mtry = mtries, oob = oob_by_mtry), aes(mtry, oob)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
geom_point(data = data.frame(mtry = mtries[which.min(oob_by_mtry)],
oob = min(oob_by_mtry)),
color = "firebrick", size = 4) +
labs(title = "Out-of-bag error against mtry",
subtitle = sprintf("p = %d features. mtry = p is bagging; smaller mtry decorrelates the trees", p_q),
x = "mtry (features considered per split)", y = "OOB error") +
theme_dspa()
imp_q <- sort(ranger::importance(rf_fit), decreasing = TRUE)[1:15]
ggplot(data.frame(v = factor(names(imp_q), levels = rev(names(imp_q))),
i = as.numeric(imp_q)), aes(i, v)) +
geom_col(fill = "steelblue") +
labs(title = "Permutation variable importance (top 15)",
subtitle = "OOB accuracy loss when each feature is shuffled",
x = "Mean decrease in accuracy", y = NULL) +
theme_dspa(10)
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = ~as.numeric(imp_q),
y = ~reorder(names(imp_q), as.numeric(imp_q)),
type = "bar", name = "Importance") |>
layout(title = "Random forest variable importance (permutation)",
xaxis = list(title = "Mean decrease in accuracy"),
yaxis = list(title = "Variable"))library(adabag)
# All three variants fitted on TRAINING data and evaluated on the untouched
# test set. The `$class` slot holds FITTED classes -- resubstitution, not a
# performance estimate.
boost_variants <- c("Breiman", "Freund", "Zhu")
boost_res <- do.call(rbind, lapply(boost_variants, function(cl) {
set.seed(123)
m <- boosting(cd ~ ., data = qol_train, mfinal = 60, coeflearn = cl)
p <- predict(m, qol_test)
data.frame(coeflearn = cl,
train_accuracy = mean(m$class == qol_train$cd),
test_accuracy = 1 - p$error,
test_kappa = confusionMatrix(
factor(p$class, levels = levels(qol_test$cd)),
qol_test$cd)$overall[["Kappa"]])
}))
boost_res |> mutate(across(where(is.numeric), \(z) round(z, 4)))Read the two accuracy columns together. Training accuracy is near 1 for every variant — that is what 60 boosting rounds do, regardless of signal. Test accuracy is far lower and is the only number that describes the model. The three variants differ in their \(\alpha\) formula (§6.25.1); on a binary outcome the difference is modest, and SAMME’s \(\ln(K-1)\) term contributes nothing when \(K=2\).
31 An honest comparison
set.seed(1234)
ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3,
classProbs = TRUE, summaryFunction = twoClassSummary,
savePredictions = "final")
# Every model on TRAINING data, with the same resamples. `preProcess` recomputes
# centring and scaling inside each fold.
fit_svm <- train(cd ~ ., data = qol_train, method = "svmRadial", metric = "ROC",
trControl = ctrl, tuneLength = 4,
preProcess = c("center", "scale"))
fit_nn <- train(cd ~ ., data = qol_train, method = "nnet", metric = "ROC",
trControl = ctrl, tuneLength = 3, trace = FALSE, maxit = 300,
preProcess = c("center", "scale"))
fit_rf2 <- train(cd ~ ., data = qol_train, method = "ranger", metric = "ROC",
trControl = ctrl, tuneLength = 3, num.trees = 300,
num.threads = 1)
fit_gbm <- train(cd ~ ., data = qol_train, method = "gbm", metric = "ROC",
trControl = ctrl, tuneLength = 3, verbose = FALSE)res_q <- resamples(list(SVM = fit_svm, NeuralNet = fit_nn,
RandomForest = fit_rf2, GBM = fit_gbm))
summary(res_q)$statistics$ROC |> round(4)#> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
#> SVM 0.4936 0.5480 0.5845 0.5842 0.6165 0.6921 0
#> NeuralNet 0.4882 0.5444 0.5705 0.5723 0.6023 0.6639 0
#> RandomForest 0.4978 0.5484 0.5675 0.5752 0.6031 0.6790 0
#> GBM 0.5188 0.5709 0.6107 0.6000 0.6256 0.6735 0res_q$values |>
dplyr::select(Resample, ends_with("~ROC")) |>
pivot_longer(-Resample, names_to = "model", values_to = "ROC") |>
mutate(model = sub("~ROC$", "", model)) |>
ggplot(aes(reorder(model, ROC, median), ROC)) +
geom_boxplot(fill = "#9EC5E8", width = 0.55, outlier.alpha = 0.4) +
geom_jitter(width = 0.08, alpha = 0.25, size = 0.8) +
coord_flip() +
labs(title = "Cross-validated AUC across 30 common resamples",
subtitle = "Identical folds for every model, so the comparison is paired",
x = NULL, y = "AUC") +
theme_dspa()
#>
#> Call:
#> summary.diff.resamples(object = diff(res_q))
#>
#> p-value adjustment: bonferroni
#> Upper diagonal: estimates of the difference
#> Lower diagonal: p-value for H0: difference = 0
#>
#> ROC
#> SVM NeuralNet RandomForest GBM
#> SVM 0.01191 0.00894 -0.01584
#> NeuralNet 1.0000 -0.00298 -0.02775
#> RandomForest 1.0000 1.0000 -0.02478
#> GBM 0.8109 0.0543 0.2306
#>
#> Sens
#> SVM NeuralNet RandomForest GBM
#> SVM -0.027016 -0.006069 -0.005273
#> NeuralNet 0.485 0.020947 0.021743
#> RandomForest 1.000 1.000 0.000797
#> GBM 1.000 0.799 1.000
#>
#> Spec
#> SVM NeuralNet RandomForest GBM
#> SVM 0.03082 0.02511 -0.01190
#> NeuralNet 0.3353 -0.00571 -0.04272
#> RandomForest 0.9051 1.0000 -0.03701
#> GBM 1.0000 0.0911 0.1479models_q <- list(SVM = fit_svm, NeuralNet = fit_nn,
RandomForest = fit_rf2, GBM = fit_gbm)
holdout_q <- data.frame(
model = c("Majority class", names(models_q)),
accuracy = round(c(NIR_q, vapply(models_q,
\(m) mean(predict(m, qol_test) == qol_test$cd), numeric(1))), 4),
kappa = round(c(0, vapply(models_q,
\(m) confusionMatrix(predict(m, qol_test), qol_test$cd)$overall[["Kappa"]],
numeric(1))), 4),
AUC = round(c(0.5, vapply(models_q, \(m) as.numeric(pROC::auc(
pROC::roc(qol_test$cd, predict(m, qol_test, type = "prob")[, "severe_disease"],
quiet = TRUE))), numeric(1))), 4))
holdout_q$acc_CI <- c(NA, vapply(models_q, function(m) {
ci <- binom.test(sum(predict(m, qol_test) == qol_test$cd),
nrow(qol_test))$conf.int
sprintf("[%.3f, %.3f]", ci[1], ci[2])
}, character(1)))
holdout_qRead the confidence intervals before the ranking. On a test set of this size the half-width is several percentage points, so models whose intervals overlap substantially are not distinguishable and ordering them is over-reading.
32 Computational complexity summary
\(n\) = observations, \(d\) = features, \(B\) = trees, \(M\) = boosting rounds, \(n_\ell\) = units in layer \(\ell\), \(E\) = epochs, \(n_{SV}\) = support vectors.
| Method | Training | Prediction (one case) | Memory | Parallel? |
|---|---|---|---|---|
| NN forward pass | — | \(O\!\left(\sum_\ell n_\ell n_{\ell-1}\right)\) | \(O\!\left(\sum_\ell n_\ell n_{\ell-1}\right)\) | Over the batch |
| NN backward pass | \(\approx 2\times\) forward | — | Same | Over the batch |
| NN full training | \(O\!\left(E\,n\sum_\ell n_\ell n_{\ell-1}\right)\) | as above | \(O(\text{params})\) | Data-parallel |
| Gram matrix | \(O(n^2 d)\) | — | \(\mathbf{O(n^2)}\) | Yes |
| SVM (SMO, typical) | \(O(n^2 d)\) | \(O(n_{SV}d)\) | \(O(n^2)\) | Poorly |
| SVM (worst case) | \(O(n^3)\) | \(O(n_{SV}d)\) | \(O(n^2)\) | Poorly |
| Linear SVM (LIBLINEAR) | \(\mathbf{O(nd)}\) | \(O(d)\) | \(O(nd)\) | Yes |
| Multi-class SVM (1-vs-1) | \(\binom{K}{2}\times\) binary cost | \(\binom{K}{2}\) evaluations | — | Over pairs |
| Nyström / random features | \(O(nm^2)\), \(m\ll n\) | \(O(md)\) | \(O(nm)\) | Yes |
| Decision tree (CART) | \(O(dn\log n)\) | \(O(\text{depth})\) | \(O(L)\) | Over features |
| Bagging / random forest | \(O(B\,m\,n\log n)\) | \(O(B\cdot\text{depth})\) | \(O(BL)\) | Embarrassingly |
| AdaBoost / gradient boosting | \(O(M\,d\,n\log n)\) | \(O(M\cdot\text{depth})\) | \(O(ML)\) | No — sequential |
| XGBoost (histogram) | \(O(M\,d\,n)\) | \(O(M\cdot\text{depth})\) | \(O(ML+nd)\) | Within each round |
| Stacking | \(\sum_j c_j + c_{\text{meta}}\) | \(\sum_j\) base costs | — | Over base learners |
Four practical consequences.
Backpropagation costs about twice a forward pass regardless of parameter count. That single fact is why networks with \(10^9\) weights are trainable; finite differences would need one forward pass per parameter.
The SVM’s \(O(n^2)\) memory is a harder wall than its time. At \(n=10^5\) the Gram matrix alone is 80 GB. Time can be bought; memory of that order usually cannot.
Forests parallelize; boosters do not. Trees in a forest are independent, so \(B\) divides across cores. Each boosting round depends on the previous one, so \(M\) is inherently sequential — which is why forests often win on wall-clock time even when boosting wins on accuracy.
Kernel approximation converts \(O(n^2)\) into \(O(n)\). Nyström and random Fourier features build an explicit \(m\)-dimensional feature map and then run a linear model. This is the standard escape when an RBF SVM is the right model and \(n\) is too large for it.
33 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | Reporting regression fit as a correlation | Cannot detect bias or scale error; \(\hat y=100+0.01y\) scores 1.0 | RMSE, MAE, held-out \(R^2\), calibration slope |
| 2 | Randomly splitting a time series | Temporal leakage; near-perfect apparent skill | Chronological or rolling-origin split |
| 3 | Scaling before splitting | Preprocessing leakage, acute with min–max | Learn constants on training; recompute per fold |
| 4 | Reporting $class or $fitted as
accuracy |
Resubstitution; near 1 for any ensemble | Predict on a held-out set |
| 5 | Cross-validating a model fitted on the full data | The “test” rows were training rows | Fit inside each fold, on that fold’s training rows only |
| 6 | Using sigmoid or tanh in a deep network | Gradients attenuate by \(\le(1/4)^L\) | ReLU family, or residual connections |
| 7 | Reading universal approximation as a guarantee | Existence is not attainability | The theorem bounds nothing about \(N\), optimization, or generalization |
| 8 | Expecting a network to extrapolate | Predictions flatten outside the training range | Parametric form, or an out-of-distribution check |
| 9 | Selecting weight decay by training loss | Always picks the least regularization | Validation loss |
| 10 | Writing the SVM objective as \(\min\lVert w\rVert/2\) | Not the QP that is solved | \(\min\frac12\lVert w\rVert^2\) |
| 11 | Omitting \(\xi_i\) from the soft-margin constraint | Optimum forces \(\xi_i=0\); collapses to hard margin | \(y_i(w^\top x_i+b)\ge 1-\xi_i\) |
| 12 | Reading \(\alpha_i\) as an importance weight | KKT forces \(\alpha_i=0\) off the margin | \(\alpha_i>0\) iff the point is a support vector |
| 13 | Tuning \(C\) and \(\gamma\) one at a time | The CV valley runs diagonally | Joint grid or random search |
| 14 | Transcribing tuned hyperparameters by hand | Silently stale after any change | Read from tune$best.parameters |
| 15 | Using the sigmoid kernel casually | Not PSD; the dual stops being concave | Check the Gram matrix’s smallest eigenvalue |
| 16 | Fitting an RBF SVM on \(n>10^5\) | \(O(n^2)\) memory before any computation | Linear SVM, Nyström, or subsample |
| 17 | Ignoring the one-vs-one cost for many classes | \(\binom{K}{2}\) machines fitted silently | 325 binary SVMs at \(K=26\) |
| 18 | Calling random forests boosting | Reverses which error term is reduced | Bagging is parallel and cuts variance |
| 19 | Attributing \(\sigma^2/B\) to the CLT | Misidentifies the result | Bilinearity of variance; the CLT is about normality |
| 20 | Conflating AdaBoost with gradient boosting | Different losses, different updates | AdaBoost reweights cases; GBM fits pseudo-residuals |
| 21 | Believing boosting cannot overfit | It does, reliably, given enough rounds | Tune \(M\) on validation data |
| 22 | Stacking on in-sample base predictions | Memorizing base learners dominate the meta-model | Out-of-fold predictions only |
| 23 | Reading impurity importance with correlated features | Importance split arbitrarily within groups | Permutation importance; group-aware methods |
| 24 | Ranking models by point accuracy | Differences below the CI width are noise | Paired tests on common resamples; report CIs |
34 Practice problems
34.1 Problem 1 — Verify backpropagation on a deeper network
Extend the gradient check of §6.6.2 to biases and to a five-layer network, and report the largest relative error.
Solution
set.seed(151)
Xp <- matrix(rnorm(6 * 20), 6, 20); Yp <- matrix(rnorm(3 * 20), 3, 20)
netp <- nn_init(c(6, 8, 7, 5, 3), seed = 17)
gp <- nn_backward(netp, nn_forward(netp, Xp), Yp)
num_b <- function(net, X, Y, l, i, eps = 1e-6) {
up <- net; up$b[[l]][i, 1] <- up$b[[l]][i, 1] + eps
dn <- net; dn$b[[l]][i, 1] <- dn$b[[l]][i, 1] - eps
(nn_loss(up, X, Y) - nn_loss(dn, X, Y)) / (2 * eps)
}
set.seed(153)
errs <- replicate(15, {
l <- sample(netp$L, 1)
if (runif(1) < 0.5) {
i <- sample(nrow(netp$W[[l]]), 1); j <- sample(ncol(netp$W[[l]]), 1)
a <- gp$dW[[l]][i, j]; n <- num_grad(netp, Xp, Yp, l, i, j)
} else {
i <- sample(nrow(netp$b[[l]]), 1)
a <- gp$db[[l]][i, 1]; n <- num_b(netp, Xp, Yp, l, i)
}
abs(a - n) / max(1e-12, abs(a) + abs(n))
})
c(max_relative_error = signif(max(errs), 3),
median_relative_error = signif(median(errs), 3))#> max_relative_error median_relative_error
#> 8.98e-10 1.44e-1034.2 Problem 2 — Measure the vanishing gradient directly
Build networks of increasing depth with sigmoid and with ReLU activations, and measure the gradient magnitude reaching the first layer.
Solution
## OLD:
# grad_at_layer1 <- function(depth, act = c("relu", "sigmoid"), width = 12, seed = 23) {
# act <- match.arg(act)
# set.seed(seed)
# sizes <- c(8, rep(width, depth), 1)
# L <- length(sizes) - 1
# W <- lapply(seq_len(L), \(l) matrix(rnorm(sizes[l+1]*sizes[l],
# sd = sqrt(2/sizes[l])),
# sizes[l+1], sizes[l]))
# X <- matrix(rnorm(8 * 32), 8, 32); Y <- matrix(rnorm(32), 1, 32)
# f <- if (act == "relu") \(z) pmax(0, z) else \(z) 1/(1+exp(-z))
# df <- if (act == "relu") \(z) (z > 0)*1 else \(z) { s <- 1/(1+exp(-z)); s*(1-s) }
#
# A <- list(X); Z <- list()
# for (l in seq_len(L)) {
# Z[[l]] <- W[[l]] %*% A[[l]]
# A[[l+1]] <- if (l < L) f(Z[[l]]) else Z[[l]]
# }
# delta <- (A[[L+1]] - Y) / 32
# for (l in L:2) delta <- (t(W[[l]]) %*% delta) * df(Z[[l-1]])
# sqrt(mean((delta %*% t(A[[1]]))^2)) # RMS gradient at layer 1
# }
grad_at_layer1 <- function(depth, act = c("relu", "sigmoid"), width = 12, seed = 23) {
act <- match.arg(act)
set.seed(seed)
sizes <- c(8, rep(width, depth), 1)
L <- length(sizes) - 1
W <- vector("list", L)
for (l in seq_len(L)) {
W[[l]] <- matrix(rnorm(sizes[l + 1] * sizes[l],
sd = sqrt(2 / sizes[l])),
nrow = sizes[l + 1], ncol = sizes[l])
}
X <- matrix(rnorm(8 * 32), nrow = 8, ncol = 32)
Y <- matrix(rnorm(32), nrow = 1, ncol = 32)
if (act == "relu") {
f <- function(z) pmax(0, z)
df <- function(z) ifelse(z > 0, 1, 0)
} else {
f <- function(z) 1 / (1 + exp(-z))
df <- function(z) {
s <- 1 / (1 + exp(-z))
s * (1 - s)
}
}
A <- list(X)
Z <- list()
for (l in seq_len(L)) {
W_l <- as.matrix(W[[l]])
A_l <- as.matrix(A[[l]])
if (ncol(W_l) != nrow(A_l)) {
stop(sprintf("Dimension mismatch at layer %d: W is %dx%d, A is %dx%d",
l, nrow(W_l), ncol(W_l), nrow(A_l), ncol(A_l)))
}
Z[[l]] <- W_l %*% A_l
# Force the activation to be a matrix with the same dimensions as Z[[l]]
if (l < L) {
A[[l + 1]] <- matrix(f(Z[[l]]), nrow = nrow(Z[[l]]), ncol = ncol(Z[[l]]))
} else {
A[[l + 1]] <- matrix(Z[[l]], nrow = nrow(Z[[l]]), ncol = ncol(Z[[l]]))
}
}
delta <- (A[[L + 1]] - Y) / 32
for (l in L:2) {
W_l <- as.matrix(W[[l]])
delta_prev <- as.matrix(delta)
Z_prev <- as.matrix(Z[[l - 1]])
delta <- (t(W_l) %*% delta_prev) * df(Z_prev)
}
sqrt(mean((delta %*% t(A[[1]]))^2))
}
depths <- c(1, 2, 4, 8, 16, 24)
gv <-
data.frame(depth = rep(depths, 2),
activation = rep(c("relu", "sigmoid"), each = length(depths)),
grad = c(vapply(depths, grad_at_layer1, numeric(1), act = "relu"),
vapply(depths, grad_at_layer1, numeric(1), act = "sigmoid")))
gv |> mutate(grad = signif(grad, 3))ggplot(gv, aes(depth, grad, color = activation)) +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_y_log10() +
scale_color_manual(values = c(relu = "#3B7DD8", sigmoid = "#D8433B")) +
labs(title = "Gradient magnitude reaching the first layer",
subtitle = "Log scale. The sigmoid network's first layer stops receiving a usable signal",
x = "Hidden layers", y = "RMS gradient at layer 1 (log scale)", color = NULL) +
theme_dspa()
34.3 Problem 3 — Quantify temporal leakage
Simulate an autocorrelated series with no predictive relationship, and show that a random split reports skill where none exists.
Solution
leak_ts <- function(phi, n = 600, seed = 157) {
set.seed(seed)
ar1 <- function() as.numeric(arima.sim(list(ar = phi), n))
X <- cbind(ar1(), ar1(), ar1())
y <- ar1() # INDEPENDENT of X by construction
d <- data.frame(X, y = y)
cut <- floor(0.75 * n)
m_t <- lm(y ~ ., d[1:cut, ])
r2_time <- 1 - sum((predict(m_t, d[(cut+1):n, ]) - y[(cut+1):n])^2) /
sum((y[(cut+1):n] - mean(y[(cut+1):n]))^2)
set.seed(seed + 1); tr <- sample(n, cut)
m_r <- lm(y ~ ., d[tr, ])
r2_rand <- 1 - sum((predict(m_r, d[-tr, ]) - y[-tr])^2) /
sum((y[-tr] - mean(y[-tr]))^2)
c(phi = phi, R2_random_split = r2_rand, R2_temporal_split = r2_time)
}
as.data.frame(do.call(rbind, lapply(c(0, 0.5, 0.9, 0.99), leak_ts))) |> round(4)X and y —
both are independent AR(1) series. As \(\phi\) rises, the random split reports an
increasingly positive \(R^2\), while
the temporal split stays near zero or negative, which is correct.
The leakage grows with the autocorrelation, and daily
search or market series sit at the far end of that range.
34.4 Problem 4 — Confirm SVM sparsity
Show that the SVM solution depends only on support vectors, and that the support-vector fraction falls as \(C\) rises.
Solution
set.seed(163)
n4 <- 400
X4 <- matrix(rnorm(n4 * 2), n4, 2)
y4 <- factor(ifelse(X4[, 1] - X4[, 2] + rnorm(n4, sd = 0.8) > 0, "a", "b"))
d4 <- data.frame(x1 = X4[, 1], x2 = X4[, 2], y = y4)
g4 <- expand.grid(x1 = seq(-3, 3, length.out = 50), x2 = seq(-3, 3, length.out = 50))
Cs4 <- c(0.01, 0.1, 1, 10, 100)
tab4 <- do.call(rbind, lapply(Cs4, function(cc) {
m <- ksvm(y ~ ., data = d4, kernel = "vanilladot", C = cc, scaled = FALSE)
sv <- SVindex(m)
m2 <- ksvm(y ~ ., data = d4[sv, ], kernel = "vanilladot", C = cc, scaled = FALSE)
data.frame(C = cc, n_SV = length(sv), sv_fraction = length(sv) / n4,
refit_identical = mean(predict(m, g4) == predict(m2, g4)))
}))#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parametersggplot(tab4, aes(C, sv_fraction)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
scale_x_log10() + scale_y_continuous(labels = scales::percent) +
labs(title = "Support-vector fraction against C",
subtitle = "Small C tolerates violations, so more points sit on or inside the margin",
x = "C (log scale)", y = "Fraction of training points that are support vectors") +
theme_dspa()
34.5 Problem 5 — Check Mercer’s condition
Test candidate kernels for positive semi-definiteness and show what an indefinite kernel does to the dual.
Solution
set.seed(167)
X5 <- matrix(rnorm(60), 30, 2)
gm <- function(f) { K <- outer(1:30, 1:30, Vectorize(\(i,j) f(X5[i,], X5[j,])))
(K + t(K)) / 2 }
me <- function(K) min(eigen(K, symmetric = TRUE, only.values = TRUE)$values)
cands <- list(
`Linear` = \(a,b) sum(a*b),
`Polynomial p=2` = \(a,b) (sum(a*b) + 1)^2,
`RBF gamma=1` = \(a,b) exp(-sum((a-b)^2)),
`Laplacian` = \(a,b) exp(-sum(abs(a-b))),
`Sigmoid k=1 d=0` = \(a,b) tanh(sum(a*b)),
`Sigmoid k=2 d=1` = \(a,b) tanh(2*sum(a*b) - 1),
`"Distance" (invalid)` = \(a,b) -sqrt(sum((a-b)^2)))
data.frame(kernel = names(cands),
min_eigenvalue = signif(vapply(cands, \(f) me(gm(f)), numeric(1)), 4),
PSD = vapply(cands, \(f) me(gm(f)) > -1e-8, logical(1)))kernlab warns about
tanhdot.
34.6 Problem 6 — Bagging helps unstable learners only
Quantify the improvement bagging gives to a deep tree, a shallow tree, and a linear model.
Solution
set.seed(173)
n6 <- 300
X6 <- matrix(rnorm(n6 * 4), n6, 4)
y6 <- factor(ifelse(X6[,1]^2 + X6[,2] + rnorm(n6, sd = 0.8) > 1, "a", "b"))
d6 <- data.frame(X6, y = y6); tr6 <- 1:200; te6 <- 201:n6
bag_it <- function(fitter, B = 50) {
votes <- replicate(B, {
idx <- sample(tr6, replace = TRUE)
as.character(fitter(d6[idx, ], d6[te6, ]))
})
apply(votes, 1, \(v) names(which.max(table(v))))
}
fitters <- list(
`Deep tree` = \(tr, te) predict(rpart::rpart(y ~ ., tr,
control = rpart::rpart.control(cp = 0, minsplit = 2)),
te, type = "class"),
`Shallow stump` = \(tr, te) predict(rpart::rpart(y ~ ., tr,
control = rpart::rpart.control(maxdepth = 1)),
te, type = "class"),
`Logistic` = \(tr, te) factor(ifelse(predict(glm(y ~ ., tr,
family = binomial()), te, type = "response") > 0.5, "b", "a"),
levels = levels(y6)))
set.seed(179)
do.call(rbind, lapply(names(fitters), function(nm) {
single <- mean(fitters[[nm]](d6[tr6, ], d6[te6, ]) == d6$y[te6])
bagged <- mean(bag_it(fitters[[nm]]) == as.character(d6$y[te6]))
data.frame(learner = nm, single = round(single, 4), bagged = round(bagged, 4),
gain = round(bagged - single, 4))
}))34.7 Problem 7 — Derive and verify AdaBoost’s \(\alpha\)
Implement AdaBoost with decision stumps and confirm that the closed-form \(\alpha_m\) matches a numerical line search on exponential loss.
Solution
set.seed(181)
n7 <- 300
X7 <- matrix(rnorm(n7 * 2), n7, 2)
y7 <- ifelse(X7[, 1] + X7[, 2] > 0.4, 1, -1)
stump <- function(X, y, w) {
best <- list(err = Inf)
for (j in 1:2) for (thr in quantile(X[, j], seq(0.05, 0.95, by = 0.05)))
for (s in c(1, -1)) {
pred <- s * ifelse(X[, j] > thr, 1, -1)
e <- sum(w * (pred != y)) / sum(w)
if (e < best$err) best <- list(err = e, j = j, thr = thr, s = s, pred = pred)
}
best
}
w <- rep(1 / n7, n7); F7 <- numeric(n7); trace7 <- NULL
for (m in 1:25) {
h <- stump(X7, y7, w)
alpha_closed <- 0.5 * log((1 - h$err) / h$err)
# Numerical line search on the SAME exponential objective
alpha_num <- optimize(\(a) sum(w * exp(-a * y7 * h$pred)),
interval = c(-5, 5))$minimum
trace7 <- rbind(trace7, data.frame(m, err = h$err,
alpha_closed, alpha_num,
diff = abs(alpha_closed - alpha_num)))
F7 <- F7 + alpha_closed * h$pred
w <- w * exp(-alpha_closed * y7 * h$pred); w <- w / sum(w)
}
head(trace7, 5) |> mutate(across(where(is.numeric), \(z) round(z, 5)))#> max_abs_difference final_training_error
#> 3.55e-06 2.00e-02ggplot(trace7, aes(m)) +
geom_line(aes(y = alpha_closed, color = "Closed form"), linewidth = 1) +
geom_point(aes(y = alpha_num, color = "Numerical line search"), size = 2) +
scale_color_manual(values = c("Closed form" = "#3B7DD8",
"Numerical line search" = "#D8433B")) +
labs(title = "AdaBoost's alpha is the exact line search for exponential loss",
x = "Boosting round m", y = expression(alpha[m]), color = NULL) +
theme_dspa()
34.8 Problem 8 — When does the kernel trick pay?
Compare an explicit polynomial feature expansion against the equivalent kernel, in both accuracy and cost, as the degree grows.
Solution
set.seed(191)
n8 <- 400; d8 <- 6
X8 <- matrix(rnorm(n8 * d8), n8, d8)
y8 <- factor(ifelse(rowSums(X8[, 1:3]^2) - 3 + rnorm(n8, sd = 0.6) > 0, "a", "b"))
expand_poly <- function(X, p) {
# Explicit monomial map up to total degree p
out <- X
if (p >= 2) for (i in 1:ncol(X)) for (j in i:ncol(X)) out <- cbind(out, X[,i]*X[,j])
if (p >= 3) for (i in 1:ncol(X)) for (j in i:ncol(X)) for (k in j:ncol(X))
out <- cbind(out, X[,i]*X[,j]*X[,k])
out
}
do.call(rbind, lapply(1:3, function(p) {
Xe <- expand_poly(X8, p)
t_exp <- system.time(m_e <- ksvm(Xe, y8, kernel = "vanilladot", C = 1,
cross = 5, scaled = TRUE))[["elapsed"]]
t_ker <- system.time(m_k <- ksvm(X8, y8, kernel = "polydot",
kpar = list(degree = p, scale = 1, offset = 1),
C = 1, cross = 5, scaled = TRUE))[["elapsed"]]
data.frame(degree = p, explicit_features = ncol(Xe),
explicit_cv_error = round(cross(m_e), 4),
kernel_cv_error = round(cross(m_k), 4),
explicit_sec = round(t_exp, 3), kernel_sec = round(t_ker, 3))
}))#> Setting default kernel parameters
#> Setting default kernel parameters
#> Setting default kernel parameters35 Checkpoint
- A colleague reports a correlation of 0.97 between predicted and observed values on held-out data. What two things do you ask for?
- Why does replacing sigmoid with ReLU make deep networks trainable, in one sentence with a number in it?
- An SVM has 40 support vectors out of 2,000 training points. What happens to the decision boundary if you delete the other 1,960 and refit?
- Your RBF SVM is perfect on training data and at chance on test data. Which hyperparameter is the likely culprit and in which direction?
- You add 500 more trees to a random forest and nothing improves. What does the variance identity say to change instead?
- Adding boosting rounds keeps lowering training error but test error started rising 200 rounds ago. Is the model broken?
Answers
- The calibration slope and intercept, and the RMSE. Correlation is invariant to affine transformation, so \(\hat y = a + b\hat y\) scores identically for any \(b>0\) — a model that is systematically biased or on the wrong scale can correlate at 0.97 and be useless. Ask also whether the data is a time series and how it was split; on an autocorrelated series a random split reports interpolation skill, not forecasting skill.
- Backpropagation multiplies one activation derivative per layer, and \(\sigma'\le\frac14\), so a 20-layer sigmoid network attenuates the gradient by at most \(4^{-20}\approx 10^{-12}\) — while ReLU’s derivative is exactly 1 on active paths, so the product does not shrink.
- Nothing. KKT complementary slackness forces \(\alpha_i=0\) for every point strictly beyond the margin, and \(\mathbf{w}=\sum_i\alpha_iy_i\mathbf{x}_i\), so non-support-vectors contribute exactly zero. Refitting on the 40 support vectors reproduces the identical boundary. This is what “the solution is sparse in the data” means.
- \(\gamma\) is too large. A wide \(\gamma\) makes the kernel narrow, so each training point influences only its immediate neighbourhood; the model places a bump on each training point and predicts the majority class everywhere else. Decrease \(\gamma\) — and tune it jointly with \(C\), since the CV valley runs diagonally and a large \(C\) compounds the same failure.
- The variance of the ensemble mean is \(\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\). If
more trees do not help, the second term is already negligible and you
are sitting on the floor \(\rho\sigma^2\). Improvement must come from
lowering \(\rho\)
(reduce
mtry, or use extremely randomized trees) or lowering \(\sigma^2\) (better features, better individual trees). If the floor is irreducible noise, no ensemble change helps. - Not broken — overfitting, which is exactly what
boosting does given enough rounds. Training error falls monotonically
because each round fits the current errors; the number of rounds \(M\) is a capacity parameter and must be
selected on validation data. Stop at the test-error minimum (or use
gbm.perf, or a validation-based early stop). There is no guarantee that more boosting is better, and no guarantee that a boosted ensemble beats its best individual member on held-out data.
36 Summary
Neural networks
- A network is a composition of affine maps and nonlinearities. Without the nonlinearity, depth collapses to a single linear map.
- Backpropagation is the chain rule organized so that all gradients cost about twice a forward pass, independent of the parameter count.
- Sigmoid and tanh saturate: \(\sigma'\le\frac14\), so gradients vanish with depth. ReLU’s derivative is 1 on active paths, which is why it replaced them.
- Universal approximation guarantees existence, not attainability, and says nothing about width, optimization, or generalization. Networks do not extrapolate beyond their training range.
- Evaluate regression with RMSE, MAE, held-out \(R^2\), and calibration — never with correlation alone.
- A time series must be split chronologically. A random split measures interpolation between adjacent, nearly identical days.
Support vector machines
- Maximizing the margin \(2/\lVert w\rVert\) means minimizing \(\frac12\lVert w\rVert^2\) — a convex QP with a unique solution.
- The dual depends on the data only through inner products, which is what makes the kernel substitution possible, and has \(n\) variables regardless of \(d\).
- KKT complementary slackness forces \(\alpha_i=0\) off the margin, so the solution is sparse in support vectors.
- The soft margin needs \(\xi_i\) in the constraint; it is equivalent to ridge-penalized hinge loss with \(\lambda=1/(2C)\).
- Kernels must satisfy Mercer’s condition or the dual stops being concave. The RBF’s feature space is infinite-dimensional and costs \(O(d)\) to use.
- \(C\) and \(\gamma\) interact and must be tuned jointly.
- Training is \(O(n^2d)\)–\(O(n^3)\) with \(O(n^2)\) memory. That memory wall, not accuracy, is why tree ensembles displaced SVMs on large data.
Ensembles
- \(\operatorname{Var}(\bar f)=\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\): more members remove the second term, decorrelation lowers the first.
- Bagging reduces variance and cannot reduce bias; it helps unstable learners and does little for stable ones.
- Random forests are bagging plus feature subsampling, fitted in parallel.
- AdaBoost minimizes exponential loss by reweighting observations; \(\alpha_m=\frac12\ln\frac{1-\epsilon_m}{\epsilon_m}\) is the exact line search.
- Gradient boosting fits each learner to the negative gradient of any differentiable loss; AdaBoost is its exponential-loss special case.
- Boosting overfits. \(M\) is a capacity parameter, and there is no guarantee the ensemble beats its best member out of sample.
- Stacking requires out-of-fold base predictions, or the meta-learner rewards memorization.
Where these threads continue
| Thread | Continues in |
|---|---|
| Association rules and text mining | NLP and rule learning |
| Clustering without labels | Unsupervised clustering |
| Nested resampling, calibration, cost-sensitive tuning | Model assessment |
| Regularization paths, LASSO, stability selection | Feature selection |
| Sequence models and recurrent architectures | Longitudinal analysis |
| Gradient descent, duality, constrained optimization | Function optimization |
| Convolutional networks, embeddings, transfer learning | Deep learning |
37 Chapter roadmap
- Chapter 1 —
Foundations. R toolchain, reproducibility conventions,
dspa_read(), simulation. - Chapter 2 — Data quality and exploratory visual analytics. Missingness, contingency tables, class imbalance.
- Chapter 3 — Linear algebra, matrix computing, and regression. Projection, SVD, conditioning, trees, BART.
- Chapter 4 — Dimensionality reduction. PCA, ICA, kernel PCA, distance concentration, t-SNE, UMAP.
- Chapter 5 — Supervised classification. Bayes error, evaluation metrics, leakage, kNN, naive Bayes, trees, forests.
- Text mining and association rules. NLP pipelines, Apriori, qualitative learning.
- Unsupervised clustering. k-means, hierarchical, spectral, Gaussian mixtures.
- Model assessment, validation, improvement. Nested resampling, calibration, hyperparameter optimization.
- Variable importance and feature selection. Ridge, LASSO, elastic net, stability selection, FDR control.
- Longitudinal and time-series analysis. Mixed models, GEE, ARIMA, recurrent architectures.
- Function optimization. Gradient descent, Newton and quasi-Newton methods, duality, Bayesian optimization.
- Deep learning. Convolutional and recurrent networks, autoencoders, embeddings, transfer learning.
38 Session information
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] parallel stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] adabag_5.0 doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2
#> [5] rpart_4.1.23 ipred_0.9-14 gbm_2.3.1 neuralnet_1.44.2
#> [9] MASS_7.3-60.0.1 ranger_0.16.0 e1071_1.7-14 kernlab_0.9-32
#> [13] caret_6.0-94 lattice_0.22-6 plotly_4.12.0 patchwork_1.3.0
#> [17] tidyr_1.3.1 dplyr_1.1.4 ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 viridisLite_0.4.2 timeDate_4032.109
#> [4] farver_2.1.2 S7_0.2.1 fastmap_1.2.0
#> [7] lazyeval_0.2.2 pROC_1.18.5 digest_0.6.37
#> [10] timechange_0.3.0 lifecycle_1.0.5 survival_3.7-0
#> [13] magrittr_2.0.3 compiler_4.3.3 rlang_1.1.5
#> [16] sass_0.4.9 tools_4.3.3 yaml_2.3.10
#> [19] data.table_1.16.4 knitr_1.51 labeling_0.4.3
#> [22] htmlwidgets_1.6.4 plyr_1.8.9 RColorBrewer_1.1-3
#> [25] withr_3.0.2 purrr_1.0.2 rgl_1.3.17
#> [28] nnet_7.3-19 grid_4.3.3 stats4_4.3.3
#> [31] future_1.33.2 gtools_3.9.5 globals_0.16.3
#> [34] scales_1.4.0 cli_3.6.3 rmarkdown_2.31
#> [37] generics_0.1.3 rlist_0.4.6.2 otel_0.2.0
#> [40] rstudioapi_0.18.0 future.apply_1.11.2 httr_1.4.7
#> [43] reshape2_1.4.4 cachem_1.1.0 proxy_0.4-27
#> [46] stringr_1.5.1 splines_4.3.3 base64enc_0.1-3
#> [49] vctrs_0.6.5 hardhat_1.4.3 Matrix_1.6-5
#> [52] jsonlite_1.8.9 listenv_0.9.1 crosstalk_1.2.1
#> [55] gower_1.0.1 jquerylib_0.1.4 recipes_1.4.0
#> [58] glue_1.8.0 parallelly_1.37.1 ConsRank_2.1.4
#> [61] codetools_0.2-20 rsample_1.2.1 lubridate_1.9.3
#> [64] stringi_1.8.4 gtable_0.3.6 tibble_3.2.1
#> [67] furrr_0.3.1 pillar_1.10.1 htmltools_0.5.8.1
#> [70] lava_1.8.0 R6_2.6.1 evaluate_1.0.3
#> [73] bslib_0.9.0 class_7.3-22 Rcpp_1.0.14
#> [76] nlme_3.1-165 prodlim_2024.06.25 xfun_0.52
#> [79] pkgconfig_2.0.3 ModelMetrics_1.2.2.2