From Linear Models to Multi-Layer Perceptrons

"Give me one matrix and I will draw you a straight line through anything. Give me two matrices with nothing between them and I will draw you the exact same line, twice the work, none the wiser."

A Linear Classifier With Delusions of Depth

Part II ended with classical recognition pipelines that fed hand-crafted features into a linear classifier, the support vector machine and its relatives in Chapter 16. The classifier was a single matrix; all the cleverness lived upstream in the features a human designed. Deep learning moves the cleverness inside the model, but the journey starts from that same linear classifier. We will build it, ask what we gain by stacking copies of it, discover the answer is "nothing", and then make one small change that turns it into a function approximator powerful enough to motivate the entire rest of this book. The illustration below captures that one move in a single picture.

1. The Linear Classifier, Revisited Beginner

A linear classifier scores an input vector $\mathbf{x} \in \mathbb{R}^d$ against $C$ classes with a single affine map. We stack the per-class weight vectors into a matrix $W \in \mathbb{R}^{C \times d}$ and add a bias $\mathbf{b} \in \mathbb{R}^C$:

$$ \mathbf{z} = W\mathbf{x} + \mathbf{b}, \qquad \mathbf{z} \in \mathbb{R}^C. $$

Each component $z_c$ is the unnormalized score (the logit) for class $c$. To turn logits into a probability distribution we apply the softmax function, which exponentiates and normalizes:

$$ p_c = \frac{e^{z_c}}{\sum_{k=1}^{C} e^{z_k}}. $$

For a vision input this $\mathbf{x}$ is, for now, a flattened image: a $28 \times 28$ grayscale picture becomes a length-784 vector by reading its pixels in row-major order. That flattening discards every spatial relationship between neighboring pixels, which is exactly the sin Chapter 19 will repair. For this section it is a feature, not a bug: it keeps the model a plain matrix multiply so we can reason about it cleanly. The code below builds the classifier with nothing but NumPy, to show there is no framework magic in a forward pass.

# Linear image classifier in pure NumPy: one affine map (W x + b) turns a
# flattened image into logits, then a numerically stable softmax converts
# those logits into a per-image probability distribution over the C classes.
import numpy as np

def softmax(z):
    # subtract the row max for numerical stability before exponentiating
    z = z - z.max(axis=1, keepdims=True)
    e = np.exp(z)
    return e / e.sum(axis=1, keepdims=True)

rng = np.random.default_rng(0)
batch, d, C = 4, 784, 10            # 4 flattened 28x28 images, 10 classes
X = rng.standard_normal((batch, d))
W = rng.standard_normal((C, d)) * 0.01
b = np.zeros(C)

Z = X @ W.T + b                     # logits, shape (4, 10)
P = softmax(Z)                      # class probabilities, shape (4, 10)
print(P.shape, P.sum(axis=1))       # (4, 10) [1. 1. 1. 1.]

The output confirms the shape and that each row of P is a valid distribution. This single linear layer is the entire model. The only question deep learning answers differently from Chapter 16 is where $W$ comes from: there it was fit by a convex solver on hand-crafted features; here it will be discovered by gradient descent (Section 18.2) directly on pixels.

2. Why Stacking Linear Layers Gains Nothing Beginner

The obvious idea for more power is depth: send $\mathbf{x}$ through one linear layer, then another. Let the first layer be $W_1, \mathbf{b}_1$ producing a hidden vector $\mathbf{h} = W_1\mathbf{x} + \mathbf{b}_1$, and the second be $W_2, \mathbf{b}_2$ producing $\mathbf{z} = W_2\mathbf{h} + \mathbf{b}_2$. Substituting one into the other:

$$ \mathbf{z} = W_2(W_1\mathbf{x} + \mathbf{b}_1) + \mathbf{b}_2 = \underbrace{(W_2 W_1)}_{W'}\,\mathbf{x} + \underbrace{(W_2\mathbf{b}_1 + \mathbf{b}_2)}_{\mathbf{b}'}. $$

The two-layer network is identically a one-layer network with weights $W' = W_2 W_1$ and bias $\mathbf{b}' = W_2\mathbf{b}_1 + \mathbf{b}_2$. No matter how many linear layers you stack, the composition is one matrix; the family of functions you can represent never grows beyond "straight cuts through input space". This is not an approximation or a training difficulty: it is an algebraic identity. The experiment below verifies it numerically.

# Numerical check of the collapse identity: run a batch through two stacked
# linear layers with no activation, then through the single equivalent layer
# whose weight is W2 @ W1, and confirm the two outputs agree exactly.
import numpy as np
rng = np.random.default_rng(1)
x = rng.standard_normal((5, 6))           # batch of 5, input dim 6
W1, b1 = rng.standard_normal((8, 6)), rng.standard_normal(8)
W2, b2 = rng.standard_normal((3, 8)), rng.standard_normal(3)

# two linear layers, no nonlinearity between them
h = x @ W1.T + b1
z_deep = h @ W2.T + b2

# the equivalent single linear layer
W_eq = W2 @ W1                            # (3, 6)
b_eq = W2 @ b1 + b2                       # (3,)
z_flat = x @ W_eq.T + b_eq

print(np.allclose(z_deep, z_flat))        # True

3. The Nonlinearity That Breaks the Collapse Beginner

Insert a nonlinear function $\sigma$ applied element-wise after the first linear layer, and the collapse argument fails because $W_2 \sigma(W_1 \mathbf{x})$ cannot be folded into a single matrix. The simplest and now-dominant choice is the rectified linear unit, $\mathrm{ReLU}(u) = \max(0, u)$, which zeroes negative values and passes positive ones unchanged. Its appeal is brutal simplicity: cheap to compute, a gradient of exactly $1$ for positive inputs (so signal does not shrink as it backpropagates through depth), and a hard zero elsewhere that yields sparse activations.

ReLU was not the first choice historically. Two older alternatives, the logistic sigmoid $\sigma(u) = 1/(1+e^{-u})$ and $\tanh$, squash inputs into a bounded range but saturate: for large magnitudes their gradient is near zero, which starves the deep layers of learning signal, the vanishing-gradient problem ReLU was introduced to sidestep.

Geometrically, each ReLU unit folds input space along a hyperplane, and a layer of them tiles space into linear regions. A linear region here just means a patch of input space inside which the same set of ReLU units is active, so the network behaves like one fixed linear map; cross into the next patch and a different unit switches on, giving a different linear map and thus a bend. Stacking layers multiplies the number of regions, which is how a network builds a curved decision boundary out of many flat pieces. Figure 18.1.1 contrasts the three activation shapes that the prose just described.

With a nonlinearity in place, a two-layer network $\mathbf{z} = W_2\,\sigma(W_1\mathbf{x} + \mathbf{b}_1) + \mathbf{b}_2$ is genuinely more expressive than a single layer. This is the multi-layer perceptron (MLP), and the layer producing $\mathbf{h}$ is the hidden layer because its outputs are never observed directly, only consumed by the next layer.

4. The Multi-Layer Perceptron Intermediate

An MLP is a chain of (linear layer, activation) pairs ending in a final linear layer that produces logits. With one hidden layer of width $H$, the forward pass on a flattened image is two matrix multiplies with a ReLU between them. The code makes a complete forward pass explicit, still in NumPy, so the only new ingredient beyond Section 1 is the single maximum(0, ...) call.

# One-hidden-layer MLP forward pass in NumPy: two weight matrices with an
# elementwise ReLU between them, plus the He-scaled initialization that keeps
# activations well-conditioned. This is the linear classifier plus one max(0, .).
import numpy as np

def softmax(z):
    z = z - z.max(axis=1, keepdims=True)
    e = np.exp(z)
    return e / e.sum(axis=1, keepdims=True)

rng = np.random.default_rng(2)
batch, d, H, C = 8, 784, 128, 10           # hidden layer of width 128
X  = rng.standard_normal((batch, d))
W1 = rng.standard_normal((H, d)) * np.sqrt(2.0 / d)   # He init for ReLU
b1 = np.zeros(H)
W2 = rng.standard_normal((C, H)) * np.sqrt(2.0 / H)
b2 = np.zeros(C)

h = np.maximum(0.0, X @ W1.T + b1)         # hidden activations (ReLU)
Z = h @ W2.T + b2                          # output logits
P = softmax(Z)
print(h.shape, P.shape)                    # (8, 128) (8, 10)

Notice the initialization scale $\sqrt{2/d}$, the He initialization matched to ReLU. Section 18.2 will explain why the scale of initial weights matters for gradient flow; for now treat it as the recipe that keeps activations from exploding or vanishing as they pass through the layer. The remarkable fact about this humble architecture is the universal approximation theorem: a single hidden layer of sufficient width can approximate any continuous function on a bounded domain to arbitrary accuracy. The theorem says nothing about how wide (it can be impractically large) or how to find the weights, which is exactly why depth and good optimization, not just width, drive practice.

5. Decision Boundaries: What the Hidden Layer Buys You Intermediate

The clearest way to feel the difference is the exclusive-or problem: four points labeled so that no straight line separates the two classes. A linear classifier provably cannot solve it (this is the historical result that froze neural-network research for a decade). A two-layer MLP solves it easily, because the hidden ReLU units carve input space into regions the output layer can then combine. The snippet evaluates an MLP whose weights are set by hand to the known XOR solution, to show the mechanism without yet invoking training.

# Solve XOR with a hand-built MLP: two hidden ReLU units detect the "exactly
# one input on" cases, and a single output combines them. This shows the
# mechanism (depth plus nonlinearity) without any training.
import numpy as np

# XOR: class 1 when exactly one input is 1
X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]], dtype=float)
y = np.array([0, 1, 1, 0])

# Hand-set weights: two hidden ReLU units detect the two "exactly one on" cases
W1 = np.array([[1.0, 1.0], [1.0, 1.0]])
b1 = np.array([0.0, -1.0])
W2 = np.array([[1.0, -2.0]])      # combine the two hidden features
b2 = np.array([0.0])

h = np.maximum(0.0, X @ W1.T + b1)
score = (h @ W2.T + b2).ravel()
pred = (score > 0.5).astype(int)
print(pred, "matches labels:", np.array_equal(pred, y))   # [0 1 1 0] matches labels: True

To see the mechanism on one input, trace $[1, 1]$: both hidden units receive $1+1=2$ before their biases, so the first unit outputs $\max(0, 2)=2$ and the second outputs $\max(0, 2-1)=1$; the output then computes $1\cdot 2 - 2\cdot 1 = 0$, which is below the $0.5$ threshold, giving class 0, exactly as XOR requires. The output [0 1 1 0] matches the labels exactly. The single hidden layer transformed the inputs into a new space where the classes are linearly separable, and the output layer drew the now-trivial line. Figure 18.1.2 makes this transformation literal: it plots all four XOR points in the original input plane, where no line separates the two classes, beside the same four points in the hidden layer's $(h_1, h_2)$ coordinates, where a single straight line does.

That is the entire principle behind every network in Part III: learn a representation in which the final decision is easy. In Section 18.2 we stop setting weights by hand and let gradient descent discover them; in Chapter 25 we will see networks learn representations so good that the final linear layer barely has to do anything.