Multilayer Perceptrons and Backpropagation › Formula

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.