| 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)