Multilayer Perceptrons and Backpropagation

A logistic regression model using the original input features draws a linear decision boundary. In the noisy two-moons example below, a straight line cuts across the crescent-shaped groups. A multilayer perceptron (MLP) can represent a curved boundary by stacking affine transformations with nonlinear activation functions. The first layer learns a new representation of the inputs, and the final layer makes a prediction from that representation.

The central difficulty is training. The hidden-layer network has more parameters than the linear baseline, and its loss is generally nonconvex in the combined layer parameters. We fit it with iterative gradient updates. Backpropagation computes the loss gradient with respect to every weight and bias by applying the chain rule systematically to the computation graph.

In this chapter we build a one-hidden-layer network from primitive NumPy operations. We derive the forward pass, the stable binary cross-entropy loss, and every gradient used in training. We check the gradients with finite differences, train the network with mini-batches, and compare it with logistic regression on a noisy nonlinear classification task. The chapter assumes familiarity with gradient descent and numerical optimization and with train/test discipline from validation strategies and data leakage.

The main advantage of an MLP is representational flexibility. A hidden layer applies several learned affine projections followed by a nonlinearity, so the network can construct features that are not present in the raw input coordinates. In a two-dimensional classification problem, this lets the decision boundary bend. In higher-dimensional tabular problems, it lets the model represent interactions such as "high income matters only when the account is new." Logistic regression on the original features has additive log-odds; to represent an income-by-account-age interaction in the log-odds, it needs additional interaction or nonlinear features.

The same hidden-layer machinery extends to several prediction tasks. A sigmoid output with binary cross-entropy gives binary classification; a softmax output with cross-entropy gives multiclass classification; a linear output with squared error gives regression. The output dimension is chosen for the target representation: one logit for binary classification, one logit per class for a softmax output, and one value per regression target. The output activation and loss must match that task. The hidden-layer chain-rule calculation follows the same principle.

The forward and backward passes can be expressed with matrix multiplications and elementwise operations, with examples arranged in batches. Mini-batch training uses fewer examples for each update than full-batch gradient descent and averages more per-example gradients than one-example stochastic gradient descent. Regularization can be added directly to the loss, and validation data can monitor the resulting training loop.

MLPs are nonconvex. The loss surface can contain many local minima, flat regions, saddle points, and steep cliffs. Gradient descent does not come with a global convergence guarantee. In practice, careful initialization, input scaling, sensible learning rates, and mini-batch noise make training work well on many problems, but the optimization is still more fragile than fitting a linear model or a convex regularized model.

MLPs are sensitive to preprocessing and hyperparameters. Features measured on very different scales can produce gradients with very different magnitudes, making a single learning rate awkward to choose. Training can stall or diverge. Hidden width, activation, learning rate, batch size, regularization strength, and training duration all affect the fitted model. A small network can underfit nonlinear structure; a large network can overfit, especially with few labeled examples.

Finally, MLPs are less interpretable than linear models. The weights of the first layer are not directly the effects of the original features, because they are mixed by the hidden units and the activation function. Interpretation usually requires post hoc methods, which we treat later in the book. For small tabular datasets where a linear model or a tree ensemble already performs well, the extra complexity of an MLP may not be justified.

With the motivation in place, the formula section answers a practical question: what exactly does the network compute, and how do we adjust every parameter so that its predictions improve? We will proceed in the same order as the implementation. First we define the forward pass and keep track of shapes. Then we define a numerically stable loss. Finally we differentiate that loss layer by layer. The order matters: backpropagation is not a separate learning algorithm, it is the chain rule applied to the forward computation we have just written down.

We begin with notation. Suppose we have $n$ training examples, each with $d$ input features. Let $x_i \in \mathbb{R}^d$ be the $i$-th input and $y_i \in \{0,1\}$ its binary label. We use a network with one hidden layer containing $h$ hidden units. The hidden layer has weight matrix $\mathbf{W} \in \mathbb{R}^{h \times d}$ and bias vector $\mathbf{b} \in \mathbb{R}^h$. The output layer has weight vector $\mathbf{v} \in \mathbb{R}^h$ and bias $c \in \mathbb{R}$.

where:

• $n$: number of training examples, with one example per row when stored in the matrix $\mathbf{X}$
• $d$: number of input features per example
• $x_i$: the $i$-th input vector, an element of $\mathbb{R}^d$
• $y_i$: the binary label of the $i$-th example, either $0$ or $1$
• $h$: number of hidden units in the single hidden layer
• $\mathbf{W}$: hidden-layer weight matrix, whose $j$-th row $\mathbf{W}_j \in \mathbb{R}^d$ holds the weights feeding hidden unit $j$
• $\mathbf{b}$: hidden-layer bias vector, whose $j$-th entry $b_j$ shifts hidden unit $j$
• $\mathbf{v}$: output-layer weight vector, whose $j$-th entry $v_j$ scales hidden activation $a_j$ in the output logit
• $c$: output-layer bias, which shifts the logit before the sigmoid

The symbol $\mathbb{R}^{h \times d}$ means the set of real matrices with $h$ rows and $d$ columns. Each row of $\mathbf{W}$ is the weight vector for one hidden unit. The bias $b_j$ shifts the $j$-th hidden unit, just as the intercept shifts a regression line. One way to read the notation is as a contract: the input has $d$ numbers, the hidden layer is allowed to create $h$ new numbers, and the output layer must combine those $h$ numbers into one prediction. The rest of the derivation is about making that contract precise and keeping the dimensions consistent.

We start with the simplest part: an affine transformation. For a single input $x$, the hidden preactivation is

$$z = \mathbf{W}x + \mathbf{b},$$

where:

• $x$: a single input vector, so that $\mathbf{W}x$ applies the $h \times d$ hidden weight matrix
• $z$: the hidden preactivation vector, the affine output before the nonlinearity

The vector $z \in \mathbb{R}^h$ contains one number per hidden unit. The $j$-th component is

$$z_j = \sum_{k=1}^{d} W_{jk}x_k + b_j.$$

where:

• $z_j$: the preactivation of hidden unit $j$
• $W_{jk}$: the weight connecting input feature $k$ to hidden unit $j$, the entry in row $j$ and column $k$ of $\mathbf{W}$
• $x_k$: the $k$-th feature of the input
• $b_j$: the bias of hidden unit $j$
• $d$: the number of input features

We read $\sum_{k=1}^{d}$ as “sum over the input features $k=1,\dots,d$.” The entry $W_{jk}$ is the weight connecting input feature $k$ to hidden unit $j$. This is exactly the same weighted sum that appears in linear regression, repeated $h$ times with its own weights for each hidden unit. The hidden layer is a bank of learned affine projections, which can respond to different directions in the input space. The bias $b_j$ shifts the point at which the $j$-th projection is centered.

Without a nonlinearity, stacking two affine layers would collapse into one affine layer:

$$\mathbf{v}^\top(\mathbf{W}x+\mathbf{b})+c = (\mathbf{v}^\top\mathbf{W})x + (\mathbf{v}^\top\mathbf{b}+c),$$

which is still affine (linear plus a constant) in $x$. A nonlinear activation allows the composition to represent functions that a single affine layer cannot. A common choice is the hyperbolic tangent,

$$a_j = \tanh(z_j) = \frac{e^{z_j}-e^{-z_j}}{e^{z_j}+e^{-z_j}},$$

where:

• $a_j$: the activation of hidden unit $j$
• $z_j$: the preactivation of hidden unit $j$ entering the activation

The derivative of tanh is

$$\frac{d}{dz_j}\tanh(z_j)=1-\tanh^2(z_j)=1-a_j^2.$$

Another common choice is the rectified linear unit,

$$a_j = \operatorname{ReLU}(z_j)=\max(0,z_j),$$

where:

• $a_j$: the activation of hidden unit $j$
• $z_j$: the preactivation of hidden unit $j$ entering the activation

The ReLU derivative is $1$ when $z_j>0$ and $0$ when $z_j<0$. At exactly zero the derivative is not defined, and implementations choose a convention, usually $0$. The sigmoid activation

$$\sigma(z)=\frac{1}{1+e^{-z}}$$

where:

• $\sigma$: the logistic (sigmoid) function mapping any real input to $(0,1)$
• $z$: the input to the sigmoid, for example an output logit

The sigmoid is often used at the output for binary classification, and sometimes in hidden layers. Its derivative is $\sigma(z)(1-\sigma(z))$. These derivatives tell us how a small change in a preactivation changes the activation. The backward pass multiplies them by the upstream loss gradient to obtain the gradient with respect to the preactivation. In a saturated region, tanh or sigmoid has a small derivative and attenuates a fixed upstream gradient. A ReLU has derivative zero on the negative side, so its local data-loss gradient is zero for that input.

The figure below makes the local exchange rates visible. Solid curves are the activations, and dashed curves are their derivatives. Where a derivative is near zero, the backward pass transmits very little gradient. The figure also shows that sigmoid and tanh flatten at large magnitudes, while ReLU stays linear on the positive side and becomes exactly zero on the negative side.

The hidden setup below makes the shared plotting module available and seeds the book plots once. Each plotting cell imports the style helpers it uses. The cell after it only prepares the data that the activation figure will display.

For a mini-batch of $B$ examples, we store the inputs in a matrix $\mathbf{X} \in \mathbb{R}^{B \times d}$, with one example per row, and use the same hidden-layer mapping as in the single-example derivation. The hidden preactivation matrix is

$$\mathbf{Z} = \mathbf{X}\mathbf{W}^\top + \mathbf{1}\mathbf{b}^\top,$$

where:

• $B$: the number of examples in the mini-batch
• $\mathbf{X}$: the $B \times d$ input matrix, with one example per row
• $\mathbf{W}^\top$: the $d \times h$ transpose of the hidden weight matrix
• $\mathbf{1}$: a column vector of ones with length $B$
• $\mathbf{b}^\top$: the $1 \times h$ transpose of the hidden bias, broadcast across rows
• $\mathbf{Z}$: the $B \times h$ hidden preactivation matrix

The transpose appears because $\mathbf{X}$ has examples as rows while $\mathbf{W}$ has hidden units as rows. The product $\mathbf{X}\mathbf{W}^\top$ has shape $B \times h$: one row per example, one column per hidden unit. Broadcasting adds the same bias vector $\mathbf{b}^\top$ to every row. The activation matrix is $\mathbf{A}=g(\mathbf{Z})$, where $g$ is applied elementwise.

The output logit is

$$s = \mathbf{A}\mathbf{v}+c,$$

where:

• $\mathbf{A}$: the $B \times h$ hidden activation matrix
• $\mathbf{v}$: the $h \times 1$ output weight vector
• $c$: the output bias
• $s$: the output logit vector of length $B$

For binary classification we convert the logit to a probability with the sigmoid:

$$p_i = \sigma(s_i)=\frac{1}{1+e^{-s_i}}.$$

where:

• $s_i$: the output logit of example $i$
• $p_i$: the predicted probability that example $i$ belongs to class $1$

The shape chain is worth checking: $\mathbf{X}$ is $B \times d$, $\mathbf{W}^\top$ is $d \times h$, and $\mathbf{Z}$ and $\mathbf{A}$ are $B \times h$. In the column-vector notation used here, $\mathbf{v}$ is $h \times 1$ and $s$ is $B \times 1$. The NumPy implementation stores these vectors as one-dimensional arrays with shapes (h,) and (B,); the scalar formulas are the same. Writing each shape beside the corresponding operation helps find dimension mismatches before coding.

The forward pass is the sequence of operations that turns inputs into predictions. For one mini-batch we compute $\mathbf{Z}$, then $\mathbf{A}=g(\mathbf{Z})$, then $s=\mathbf{A}\mathbf{v}+c$, then $p=\sigma(s)$. The forward pass stores $\mathbf{X}$, $\mathbf{Z}$, and $\mathbf{A}$ in a cache. This tanh backward pass reuses $\mathbf{A}$ to compute the activation derivatives elementwise as $1-a_{ij}^2$; it receives $\mathbf{X}$ separately and does not read the cached $\mathbf{Z}$. It is useful to think of the forward pass as building a cache: every quantity that will be needed later is stored now, so the backward pass can reuse it instead of recomputing it.

For binary labels, the negative log-likelihood of one example is

$$\ell_i = -\left[y_i\log p_i + (1-y_i)\log(1-p_i)\right].$$

where:

• $\ell_i$: the binary cross-entropy loss for example $i$
• $y_i$: the binary label, either $0$ or $1$
• $p_i$: the predicted probability of class $1$ for example $i$

This is the binary cross-entropy loss. If $y_i=1$, only $-\log p_i$ remains; if $y_i=0$, only $-\log(1-p_i)$ remains. The loss is small when the predicted probability agrees with the label and large when it disagrees. In that sense the loss is a precise way of saying “how surprised should the model be by the observed label?”

For a very negative logit, a naive sigmoid calculation can overflow in $e^{-s_i}$ and round the predicted probability to zero. Taking its logarithm then produces negative infinity even though the mathematical log probability is finite. We compute the loss from the logit instead. Since $p=\sigma(s)=1/(1+e^{-s})$, we have

$$-\log p = \log(1+e^{-s})$$

and

$$-\log(1-p) = \log(1+e^{s}).$$

where:

• $p$: the predicted probability for the example, equal to $\sigma(s)$
• $s$: the output logit for the example

Combining the two cases gives

$$\ell_i = \log(1+e^{s_i}) - y_i s_i.$$

where:

• $\ell_i$: the binary cross-entropy loss for example $i$, computed from the logit
• $s_i$: the output logit of example $i$
• $y_i$: the binary label, either $0$ or $1$

The expression $\log(1+e^{s_i})$ is exactly $\log(e^0+e^{s_i})$, which is a log-sum-exp written in the form $\log(e^0 + e^{s_i})$. Numerical libraries provide a stable function often called logaddexp, which computes $\log(e^a+e^b)$ without overflowing. Thus we compute

$$\ell_i = \operatorname{logaddexp}(0,s_i) - y_i s_i.$$

For label 1 the loss decreases as the logit grows; for label 0 it increases. The logaddexp calculation avoids forming an overflowing exponential directly. The plotted losses remain finite over the displayed logit range.

where:

• $\operatorname{logaddexp}(a,b)$: a numerically stable computation of $\log(e^a+e^b)$
• $s_i$: the output logit of example $i$
• $y_i$: the binary label, either $0$ or $1$

For a batch we average:

$$L=\frac{1}{B}\sum_{i=1}^{B}\ell_i.$$

where:

• $L$: the mean binary cross-entropy over the $B$ examples
• $B$: the number of examples in the batch
• $\ell_i$: the binary cross-entropy loss for example $i$

The same idea extends to multiclass classification. If the network produces logits $s_{i1},\dots,s_{iK}$ for $K$ classes, the softmax probability of class $k$ is

$$\hat p_{ik} = \frac{e^{s_{ik}}}{\sum_{m=1}^{K}e^{s_{im}}},$$

where:

• $\hat p_{ik}$: the predicted probability that example $i$ belongs to class $k$
• $s_{ik}$: the logit for class $k$ on example $i$
• $K$: the number of classes

For this multiclass formula, let $y_i\in\{1,\dots,K\}$ denote the true class index, replacing the earlier binary encoding. The multiclass cross-entropy is

$$\ell_i = -s_{i,y_i} + \log\sum_{m=1}^{K}e^{s_{im}}.$$

where:

• $\ell_i$: the cross-entropy loss for example $i$
• $K$: the number of classes
• $s_{im}$: the logit for class $m$ on example $i$
• $s_{i,y_i}$: the logit for the true class $y_i$ on example $i$

The second term is computed with a stable log-sum-exp operation. In this chapter we focus on the binary case, but the backward pass for softmax cross-entropy has the same clean form: the gradient with respect to the logits is predicted probability minus one-hot label.

Backpropagation is the chain rule applied from the output backward. We first compute the derivative of the loss with respect to the output logit $s_i$. A useful discipline is to derive the gradient for one scalar example first, then average over the batch; the only change is the factor $1/B$, which appears in the delta defined below. For one example,

$$\frac{\partial \ell}{\partial s} = \frac{\partial \ell}{\partial p}\frac{\partial p}{\partial s}.$$

where:

• $\ell$: the binary cross-entropy loss for one example
• $s$: the output logit for that example
• $p=\sigma(s)$: the predicted probability for that example

We read $\partial \ell/\partial s$ as the partial derivative of the loss with respect to the logit $s$. The derivative of the binary cross-entropy with respect to $p$ is

$$\frac{\partial \ell}{\partial p} = -\frac{y}{p}+\frac{1-y}{1-p}.$$

where:

• $\ell$: the binary cross-entropy loss for one example
• $p$: the predicted probability for that example
• $y$: the binary label of that example

The derivative of the sigmoid is $\partial p/\partial s=p(1-p)$. Multiplying gives

$$\frac{\partial \ell}{\partial s} = \left(-\frac{y}{p}+\frac{1-y}{1-p}\right)p(1-p) = p-y.$$

where:

• $\partial \ell/\partial s$: the derivative of the loss with respect to the logit
• $p$: the predicted probability
• $y$: the binary label

This is a simple result: the gradient at the logit is the predicted probability minus the label. If the network predicts $p=0.9$ and the label is $1$, the gradient is $-0.1$, so gradient descent increases the logit slightly. If the label is $0$, the gradient is $0.9$, a larger correction. The same pattern will reappear in multiclass classification, where the gradient is the predicted probability vector minus the one-hot label.

For a batch with mean loss, we divide by $B$. Define

$$\delta_i = \frac{p_i-y_i}{B}.$$

where:

• $\delta_i$: the scaled output delta for example $i$
• $p_i$: the predicted probability for example $i$
• $y_i$: the binary label for example $i$
• $B$: the number of examples in the batch, so division by $B$ reflects the mean loss

The gradient with respect to the output weights is

$$\frac{\partial L}{\partial v_j} = \sum_{i=1}^{B}\delta_i a_{ij},$$

where:

• $L$: the mean loss over the batch
• $v_j$: the $j$-th output weight
• $a_{ij}$: the activation of hidden unit $j$ on example $i$

This holds because $s_i=\sum_j v_j a_{ij}+c$, so $\partial s_i/\partial v_j=a_{ij}$. In matrix form,

$$\nabla_{\mathbf{v}}L = \mathbf{A}^\top \delta,$$

where:

• $\nabla_{\mathbf{v}}L$: the gradient of the mean loss with respect to the output weight vector
• $\mathbf{A}$: the $B \times h$ hidden activation matrix
• $\delta$: the $B \times 1$ vector of $\delta_i$ values

The transpose is needed because $\mathbf{A}$ is $B \times h$ and $\delta$ is $B \times 1$, so $\mathbf{A}^\top\delta$ is $h \times 1$, matching $\mathbf{v}$. The gradient with respect to the output bias is

$$\frac{\partial L}{\partial c}=\sum_{i=1}^{B}\delta_i,$$

where:

• $c$: the output bias
• $\delta_i$: the scaled output delta for example $i$

Now we propagate backward into the hidden layer. For hidden unit $j$ and example $i$, the logit depends on $a_{ij}$ through $s_i=\sum_j v_j a_{ij}+c$, so

$$\frac{\partial L}{\partial a_{ij}}=\delta_i v_j.$$

where:

• $\partial L/\partial a_{ij}$: the derivative of the mean loss with respect to the activation of hidden unit $j$ on example $i$
• $\delta_i$: the scaled output delta for example $i$
• $v_j$: the output weight for hidden unit $j$

The activation $a_{ij}=g(z_{ij})$ contributes its local derivative:

$$\frac{\partial L}{\partial z_{ij}} = \delta_i v_j\,g'(z_{ij}).$$

where:

• $\partial L/\partial z_{ij}$: the derivative of the mean loss with respect to the preactivation of hidden unit $j$ on example $i$
• $\delta_i$: the scaled output delta for example $i$
• $v_j$: the output weight for hidden unit $j$
• $g'(z_{ij})$: the derivative of the hidden activation evaluated at the preactivation

For $\tanh$, $g'(z_{ij})=1-a_{ij}^2$. Define the hidden delta matrix

$$\Delta_{ij} = \delta_i v_j\,(1-a_{ij}^2),$$

where:

• $\Delta_{ij}$: the hidden delta for hidden unit $j$ on example $i$
• $\delta_i$: the scaled output delta for example $i$
• $v_j$: the output weight for hidden unit $j$
• $a_{ij}$: the activation of hidden unit $j$ on example $i$

In matrix notation,

$$\Delta = (\delta\mathbf{v}^\top)\odot(1-\mathbf{A}\odot\mathbf{A}),$$

where:

• $\Delta$: the $B \times h$ hidden delta matrix
• $\delta$: the $B \times 1$ vector of scaled output deltas
• $\mathbf{v}$: the $h \times 1$ output weight vector, so $\delta\mathbf{v}^\top$ is a $B \times h$ outer product
• $\mathbf{A}$: the $B \times h$ hidden activation matrix
• $\odot$: elementwise (Hadamard) multiplication

The outer product $\delta\mathbf{v}^\top$ has shape $B \times h$, and $1-\mathbf{A}\odot\mathbf{A}$ has the same shape. The hidden preactivation is $z_{ij}=\sum_k W_{jk}x_{ik}+b_j$, so

$$\frac{\partial L}{\partial W_{jk}} = \sum_{i=1}^{B}\Delta_{ij}x_{ik}, \qquad \frac{\partial L}{\partial b_j} = \sum_{i=1}^{B}\Delta_{ij}.$$

where:

• $\partial L/\partial W_{jk}$: the derivative of the mean loss with respect to the hidden weight connecting input feature $k$ to hidden unit $j$
• $\partial L/\partial b_j$: the derivative of the mean loss with respect to the hidden bias of unit $j$
• $\Delta_{ij}$: the hidden delta for hidden unit $j$ on example $i$
• $x_{ik}$: the $k$-th feature of example $i$

In matrix form,

$$\nabla_{\mathbf{W}}L = \Delta^\top \mathbf{X}, \qquad \nabla_{\mathbf{b}}L = \Delta^\top \mathbf{1}.$$

where:

• $\nabla_{\mathbf{W}}L$: the gradient of the mean loss with respect to the hidden weight matrix
• $\nabla_{\mathbf{b}}L$: the gradient of the mean loss with respect to the hidden bias vector
• $\Delta^\top$: the $h \times B$ transpose of the hidden delta matrix
• $\mathbf{X}$: the $B \times d$ input matrix
• $\mathbf{1}$: the $B \times 1$ column vector of ones

Again the shapes match: $\Delta^\top$ is $h \times B$, $\mathbf{X}$ is $B \times d$, so $\Delta^\top\mathbf{X}$ is $h \times d$, the same shape as $\mathbf{W}$. The gradient with respect to $\mathbf{b}$ is the column sum of $\Delta$. Notice how each matrix expression encodes a sum: the transpose aligns examples with the quantity being summed, and the matrix product performs the sum over examples in one operation.

If we add L2 regularization, the loss becomes

$$L_{\text{reg}} = L + \frac{\lambda}{2}\left(\sum_{j,k}W_{jk}^2+\sum_j v_j^2\right).$$

where:

• $L_{\text{reg}}$: the regularized mean loss
• $L$: the unregularized mean binary cross-entropy
• $\lambda$: the L2 regularization strength, with $\lambda \geq 0$
• $W_{jk}$: the hidden weight connecting input feature $k$ to hidden unit $j$
• $v_j$: the output weight for hidden unit $j$

The factor $1/2$ makes the derivative of the penalty with respect to $W_{jk}$ equal to $\lambda W_{jk}$, and similarly for $v_j$. We do not regularize the biases by default, because they do not control the scale of the learned features in the same way.

Gradient checking verifies the backward pass by comparing analytic gradients with finite differences. It is the unit test of a hand-written backward pass. For a parameter $\theta_i$,

$$\frac{\partial L}{\partial \theta_i} \approx \frac{L(\theta_i+\epsilon)-L(\theta_i-\epsilon)}{2\epsilon}.$$

where:

• $\theta_i$: the $i$-th scalar parameter being checked
• $L(\theta_i+\epsilon)$: the loss with parameter $\theta_i$ increased by the perturbation size $\epsilon$
• $L(\theta_i-\epsilon)$: the loss with parameter $\theta_i$ decreased by $\epsilon$
• $\epsilon$: the finite-difference step size

Perturb one parameter by $\epsilon$ in each direction, compute the two losses, and divide their difference by $2\epsilon$. For a sufficiently smooth loss, the centered formula has truncation error of order $\epsilon^2$. Floating-point cancellation grows as the step becomes too small, so this is not a guarantee of greater accuracy for every step size. A relative error below about $10^{-4}$ is a useful diagnostic for a smooth float64 example at a suitable step size, but near-zero gradients and activation kinks need separate inspection.

With a sufficiently wide hidden layer, a suitable continuous sigmoidal activation, and a linear output, a network can approximate continuous functions on a compact input domain to arbitrary accuracy. This is an existence statement about a family of functions, not a guarantee that gradient descent will find the required parameters from finite data. Our binary classifier instead uses a sigmoid output to keep predictions between zero and one. Approximation capacity does not remove the need for regularization and validation.

The joint loss is generally nonconvex in $\mathbf{W}$, $\mathbf{b}$, $\mathbf{v}$, and $c$. It also has a permutation symmetry: permuting hidden units and applying the same permutation to rows of $\mathbf{W}$, entries of $\mathbf{b}$, and entries of $\mathbf{v}$ leaves the network output unchanged. Distinct parameter vectors can therefore represent the same function. Random initialization helps give hidden units different starting points. If every weight and bias in this tanh network starts at zero, hidden activations and output weights are zero; the hidden features cannot begin learning under the stated gradient updates.

The gradient with respect to $W_{jk}$ contains the input factor $x_k$. Holding the other gradient factors fixed, rescaling that feature rescales this gradient. Standardization brings marginal feature scales closer together, but it does not remove correlations or guarantee a well-conditioned loss. Saturating tanh and sigmoid activations have small local derivatives at large magnitudes. These multipliers can attenuate gradient flow; ReLU avoids positive-side saturation but has zero derivative on its negative side.

To make the formula concrete, we compute it by hand. The following example is deliberately tiny: two input features, two hidden units, and one binary output. Its purpose is not to be realistic but to expose every number that the code will later compute. We will carry the forward pass, the backward pass, and one gradient step by hand, so that each matrix in the implementation has a concrete scalar counterpart. Calculations retain full precision; displayed results are rounded to four decimal places.

Suppose we have two input features, two hidden units, and one binary output. Let the input be

$$x = \begin{bmatrix}1.0\\-1.0\end{bmatrix}, \qquad y=1.$$

Let the hidden weights and biases be

$$\mathbf{W} = \begin{bmatrix} 0.5 & -0.5\\ 1.0 & 0.5 \end{bmatrix}, \qquad \mathbf{b}= \begin{bmatrix} 0.1\\-0.2 \end{bmatrix},$$

and let the output weights and bias be

$$\mathbf{v}= \begin{bmatrix} 1.2\\-0.8 \end{bmatrix}, \qquad c=0.05.$$

We use $\tanh$ in the hidden layer and sigmoid at the output.

We compute the hidden preactivations one unit at a time. For the first hidden unit,

$$z_1 = 0.5(1.0)+(-0.5)(-1.0)+0.1 = 0.5+0.5+0.1=1.1.$$

Then

$$a_1=\tanh(1.1)\approx 0.8005.$$

For the second hidden unit,

$$z_2 = 1.0(1.0)+0.5(-1.0)-0.2 = 1.0-0.5-0.2=0.3,$$

so

$$a_2=\tanh(0.3)\approx 0.2913.$$

The output logit is

$$s = 1.2a_1-0.8a_2+0.05\approx 0.7775.$$

The predicted probability is

$$p=\sigma(s)=\frac{1}{1+e^{-s}}\approx 0.6852.$$

where:

• $\sigma$: the sigmoid function mapping the logit to a probability
• $p$: the predicted probability that the example belongs to class $1$

The loss is

$$\ell = -\log(p)\approx 0.3781.$$

where:

• $\ell$: the binary cross-entropy loss for this single example
• $\log$: the natural logarithm

Equivalently, using the stable logit form,

$$\ell=\operatorname{logaddexp}(0,s)-s\approx 0.3781.$$

We can read this forward pass as a pipeline: the input is projected into two hidden numbers, each hidden number is squashed by tanh, and the two squashed numbers are combined into a logit. The predicted probability is above $0.5$, so the network already leans toward class $1$, and the observed label is $1$. Since the assigned probability $p\approx0.6852$ is below $1$, its binary cross-entropy $-\log p\approx0.3781$ is positive.

We now reverse the forward pass. The output gradient is

$$\delta = p-y \approx -0.3148,$$

The output weight gradients are

$$\frac{\partial \ell}{\partial v_1}=\delta a_1\approx -0.2520,$$ $$\frac{\partial \ell}{\partial v_2}=\delta a_2\approx -0.0917,$$

where:

• $\delta = p - y$: the output delta for this single example, without the $1/B$ scaling used for a batch mean
• $a_1, a_2$: the hidden activations computed in the forward pass
• $v_1, v_2$: the output weights

and the output bias gradient is

$$\frac{\partial \ell}{\partial c}=\delta\approx -0.3148.$$

For the hidden units, we first compute

$$\frac{\partial \ell}{\partial a_1}=\delta v_1\approx -0.3778,$$ $$\frac{\partial \ell}{\partial a_2}=\delta v_2\approx 0.2519.$$

The tanh derivatives are

$$1-a_1^2\approx 0.3592, \qquad 1-a_2^2\approx 0.9151,$$

Therefore the hidden preactivation gradients are

$$\Delta_1=\delta v_1(1-a_1^2)\approx -0.1357,$$ $$\Delta_2=\delta v_2(1-a_2^2)\approx 0.2305.$$

Now we can compute the weight gradients. Since $z_1=W_{11}x_1+W_{12}x_2+b_1$,

$$\frac{\partial \ell}{\partial W_{11}}=\Delta_1 x_1\approx -0.1357,$$ $$\frac{\partial \ell}{\partial W_{12}}=\Delta_1 x_2\approx 0.1357.$$

Similarly,

$$\frac{\partial \ell}{\partial W_{21}}=\Delta_2 x_1\approx 0.2305,$$ $$\frac{\partial \ell}{\partial W_{22}}=\Delta_2 x_2\approx -0.2305.$$

The hidden bias gradients are simply $\Delta_1$ and $\Delta_2$:

$$\frac{\partial \ell}{\partial b_1}\approx -0.1357, \qquad \frac{\partial \ell}{\partial b_2}\approx 0.2305.$$

At this point every gradient has been computed by following the chain rule one link at a time. The output delta is the derivative of the loss with respect to the logit; its sign tells us which local change in the logit would reduce the loss. Multiplying by the hidden activations tells us how to change the output weights. Multiplying the output delta by the output weights and tanh derivative gives the loss gradient with respect to each hidden preactivation. Finally, multiplying by the input tells us how to change the hidden weights.

Gradient descent subtracts a multiple of the gradient. With learning rate $\eta=0.5$, the new first hidden weight is

$$W_{11}^{\text{new}}=W_{11}-\eta\frac{\partial\ell}{\partial W_{11}}\approx 0.5679.$$

where:

• $W_{11}^{\text{new}}$: the updated hidden weight after the gradient step
• $0.5$: the initial value of $W_{11}$
• $0.5$: the learning rate $\eta$
• $-0.1357$: the gradient $\partial \ell/\partial W_{11}$

The new first output weight is

$$v_1^{\text{new}}=v_1-\eta\frac{\partial\ell}{\partial v_1}\approx 1.3260.$$

where:

• $v_1^{\text{new}}$: the updated first output weight after the gradient step
• $1.2$: the initial value of $v_1$
• $0.5$: the learning rate $\eta$
• $-0.2520$: the rounded gradient $\partial \ell/\partial v_1$

Applying the same rule at full precision and rounding the displayed values gives

$$\begin{aligned} W_{11}&\approx0.5679, & W_{12}&\approx-0.5679,\\ W_{21}&\approx0.8847, & W_{22}&\approx0.6153,\\ b_1&\approx0.1679, & b_2&\approx-0.3153,\\ v_1&\approx1.3260, & v_2&\approx-0.7541,\\ c&\approx0.2074. \end{aligned}$$

A new forward pass with these updated parameters gives $z_1\approx 1.3036$, $a_1\approx 0.8626$, $z_2\approx -0.0458$, $a_2\approx -0.0457$, $s\approx 1.3858$, $p\approx 0.7999$, and loss $\approx 0.2232$. The loss decreased from about $0.378$ to about $0.223$, which is the expected effect of a step in the negative gradient direction. It is a small confirmation that the signs and scales in the backward pass are consistent with the goal of reducing the loss.

We choose NumPy for the network because the goal is to expose every operation in the forward and backward passes. We use scikit-learn for data generation, train/test splitting, standardization, evaluation metrics, and the logistic regression baseline. The network itself is implemented with matrix operations and explicit gradients, so we can compare its predictions and check its gradients against finite differences. The point is not to hide the network behind a high-level API; it is to make the path from a mathematical expression to an executable update visible. All code cells below declare explicit Quarto directives so that computation is shown, while any printed results are displayed from calculated variables in a separate output cell.

We generate a noisy two-moons classification problem with $600$ examples, split it into training and test sets, and fit the scaler on the training set only. Fitting the scaler on the full dataset would leak test-set information into training, a point we emphasized in validation strategies and data leakage.