Deep Neural Network Backpropagation

The fully connected network developed here has several hidden layers followed by a binary output. Each layer transforms the preceding activations. We extend the one-hidden-layer implementation into loops over layers, then check that the backward loop differentiates the same loss as the forward loop. Run the body code blocks in order.

Describing an L-layer network

Let \(L\) count parameterized layers, including the output but excluding the input. Let \(n_l\) be the number of units in layer \(l\), with \(n_0=n_x\) input features. For \(m\) examples stored in columns, \(A^{[0]}=X\in\mathbb{R}^{n_0\times m}\). Each layer owns \(W^{[l]}\in\mathbb{R}^{n_l\times n_{l-1}}\) and \(b^{[l]}\in\mathbb{R}^{n_l\times1}\). Thus sizes \([2,3,2,1]\) specify two inputs, two hidden layers of widths 3 and 2, and one output: \(L=3\).

The total parameter count is \(\sum_{l=1}^{L}(n_ln_{l-1}+n_l)\). For \([2,3,2,1]\), this is \((3\times2+3)+(2\times3+2)+(1\times2+1)=20\). It depends on layer widths, not on batch size.

A complete one-hidden-layer implementation of this is in Building a Shallow Neural Network with NumPy.

Forward propagation

For each fully connected layer, compute \(Z^{[l]}=W^{[l]}A^{[l-1]}+b^{[l]}\), followed by \(A^{[l]}=g^{[l]}(Z^{[l]})\). Both arrays have shape \((n_l,m)\). This implementation uses elementwise ReLU in hidden layers and sigmoid at its single binary output. Softmax or regression outputs require different output and loss rules.

import numpy as np

def linear_activation_forward(A_prev, W, b, activation):
    assert A_prev.ndim == W.ndim == 2
    assert W.shape[1] == A_prev.shape[0]
    assert b.shape == (W.shape[0], 1)
    Z = W @ A_prev + b
    if activation == "relu":
        A = np.maximum(0, Z)
    elif activation == "sigmoid":
        A = np.empty_like(Z, dtype=float)
        positive = Z >= 0
        A[positive] = 1 / (1 + np.exp(-Z[positive]))
        exp_z = np.exp(Z[~positive])
        A[~positive] = exp_z / (1 + exp_z)
    else:
        raise ValueError("unsupported activation")
    return A, (A_prev, W, b, Z)

Why caches matter

This implementation stores the previous activation, current weights, bias, and pre-activation in each layer’s cache. The backward rules use the first two values for matrix products and the pre-activation for the ReLU mask; the bias is retained for checking its shape. These are references to arrays, not frozen copies. Finish the backward pass before updating weights or overwriting activations, or the cache may no longer describe the forward calculation. Other implementations can save less data or recompute it.

Backward propagation layer by layer

For binary labels \(Y_i\in\{0,1\}\), the objective is the mean cross-entropy \(J=-\frac1m\sum_i[Y_i\log A^{[L]}_i+(1-Y_i)\log(1-A^{[L]}_i)]\). The code below evaluates it from logits for numerical stability.

Write \(dZ^{[l]}=\partial J/\partial Z^{[l]}\), where \(J\) is the mean loss over the batch. Then the affine backward equations are \(dW^{[l]}=dZ^{[l]}(A^{[l-1]})^T\), \(db^{[l]}=\sum_{i=1}^m dZ^{[l]}_{:,i}\), and \(dA^{[l-1]}=(W^{[l]})^T dZ^{[l]}\). The batch averaging will enter once at the output. It is already included in every gradient passed toward earlier layers.

def linear_backward(dZ, cache):
    A_prev, W, b, Z = cache
    assert dZ.shape == Z.shape
    dW = dZ @ A_prev.T
    db = np.sum(dZ, axis=1, keepdims=True)
    dA_prev = W.T @ dZ
    assert dW.shape == W.shape and db.shape == b.shape
    assert dA_prev.shape == A_prev.shape
    return dA_prev, dW, db

For sigmoid and mean binary cross-entropy, begin with \(dZ^{[L]}=(A^{[L]}-Y)/m\). For a hidden ReLU, use \(dZ^{[l]}=dA^{[l]}\odot\mathbf{1}(Z^{[l]}>0)\): multiply by 1 at positive pre-activations and 0 at negative ones. At zero, ReLU has no ordinary derivative; this code chooses 0 as its backward convention. Process layers from \(L\) down to 1 so each receives the gradient already computed by the layer after it.

Getting the matrix dimensions right

For these dense array gradients, each derivative array has the shape of its parameter or activation. This catches dimension mistakes but cannot detect a wrong sign, missing averaging factor, or incorrect values. For layer \(l\) with \(n_l\) units:

Tensor	Shape	Depends on \(m\)?
\(W^{[l]}\)	\((n_l,\,n_{l-1})\)	no
\(b^{[l]}\)	\((n_l,\,1)\)	no
\(Z^{[l]},A^{[l]}\)	\((n_l,\,m)\)	yes
\(dW^{[l]}\)	same as \(W^{[l]}\)	no
\(db^{[l]}\)	same as \(b^{[l]}\)	no

Check parameter shapes before the forward calculation and gradient shapes before the update. Matching shapes are necessary for this implementation; numerical gradient checks test the values as well.

import numpy as np

def check_shapes(params, grads, layer_dims):
    L = len(layer_dims) - 1
    for l in range(1, L + 1):
        n_out, n_in = layer_dims[l], layer_dims[l - 1]
        assert params[f"W{l}"].shape == (n_out, n_in), f"W{l}"
        assert params[f"b{l}"].shape == (n_out, 1), f"b{l}"
        assert grads[f"dW{l}"].shape == params[f"W{l}"].shape, f"dW{l}"
        assert grads[f"db{l}"].shape == params[f"b{l}"].shape, f"db{l}"
    return "shapes ok"

def parameter_count(layer_dims):
    return sum(layer_dims[l] * layer_dims[l - 1] + layer_dims[l]
               for l in range(1, len(layer_dims)))

print(parameter_count([784, 256, 128, 64, 10]))
# 242762

A bias of shape \((n_l,1)\) adds one value per unit across examples. A one-dimensional bias \((n_l,)\) aligns with the last dimension of \((n_l,m)\). With \((n_l,m)=(4,10)\), it raises an error. With \((4,4)\), it runs but adds per-example offsets. With \((4,1)\), it can expand the result to \((4,4)\). The outcome depends on the dimensions, so validate the intended bias shape explicitly.

Summing dZ over axis=1 collects example contributions for each bias. keepdims=True preserves shape \((n_l,1)\); an explicit reshape could also preserve that shape. No additional division by \(m\) is needed with the gradient convention used here.

Connecting the layer functions

The parameter dictionary contains one weight matrix and bias per layer. The cache list follows forward order: Python index 0 is mathematical layer 1. After differentiating layer 3, the loop applies layer 2’s ReLU mask before differentiating layer 2, then repeats for layer 1. The initial scale is variance 2/fan_in for ReLU weights and 1/fan_in at the output; these are starting choices, not guarantees that activations remain well scaled.

def initialize(layer_dims, seed=3):
    rng = np.random.default_rng(seed)
    params = {}
    L = len(layer_dims) - 1
    for l in range(1, L + 1):
        fan_in, fan_out = layer_dims[l-1], layer_dims[l]
        variance = (2.0 if l < L else 1.0) / fan_in
        params[f"W{l}"] = rng.normal(size=(fan_out, fan_in)) * np.sqrt(variance)
        params[f"b{l}"] = np.zeros((fan_out, 1))
    return params

def model_forward(X, params):
    L = len(params) // 2
    A, caches = X, []
    for l in range(1, L + 1):
        activation = "sigmoid" if l == L else "relu"
        A, cache = linear_activation_forward(
            A, params[f"W{l}"], params[f"b{l}"], activation
        )
        caches.append(cache)
    return A, caches

def mean_loss(caches, Y):
    Z = caches[-1][3]
    assert Z.shape == Y.shape and Z.shape[0] == 1
    assert np.all((Y == 0) | (Y == 1))
    return float(np.mean(np.maximum(Z, 0) - Z * Y
                         + np.log1p(np.exp(-np.abs(Z)))))

def model_backward(AL, Y, caches):
    assert AL.shape == Y.shape
    dZ = (AL - Y) / Y.shape[1]
    grads = {}
    for index in reversed(range(len(caches))):
        dA_prev, dW, db = linear_backward(dZ, caches[index])
        l = index + 1
        grads[f"dW{l}"], grads[f"db{l}"] = dW, db
        if index > 0:
            previous_Z = caches[index-1][3]
            dZ = dA_prev * (previous_Z > 0)
    return grads

layer_dims = [2, 3, 2, 1]
X = np.array([[0.2, -0.4, 0.7], [0.5, 0.3, -0.2]])
Y = np.array([[0., 1., 1.]])
params = initialize(layer_dims)
for l in range(1, len(layer_dims) - 1):
    params[f"b{l}"].fill(0.1)
AL, caches = model_forward(X, params)
grads = model_backward(AL, Y, caches)
print("layer shapes", [cache[3].shape for cache in caches])
print(check_shapes(params, grads, layer_dims))
print("parameters", parameter_count(layer_dims))
print("loss", round(mean_loss(caches, Y), 6))
# layer shapes [(3, 3), (2, 3), (1, 3)]
# shapes ok
# parameters 20
# loss 1.377523

The mean loss is evaluated from logits with a stable binary cross-entropy expression. For the later finite-difference check, we set the hidden biases to 0.1 to keep this example away from ReLU kinks. This first call runs one forward and backward pass. Changing the hidden widths or adding entries before the final 1 in layer_dims changes how many times the loops run; the single-output sigmoid assumption stays the same.

Why depth can help

Depth enables hierarchical composition. In image tasks, earlier layers may respond to edges, middle layers to textures or parts, and later layers to task-specific combinations. In language tasks, layers can progressively combine token context. These interpretations are useful, but they are not guaranteed labels assigned to individual units.

Common implementation failures

• An output activation that does not match the loss.
• Incorrect transposes in \(dW\) or \(dA_{\mathrm{prev}}\).
• Omitting division by \(m\) or applying it twice.
• Updating parameters before all gradients are computed.
• Reusing a cache from the wrong layer.
• In-place operations that overwrite values needed by backward propagation.
• Large or tiny activation scales caused by poor initialization.

The gradient check that verifies this, and an ordered debugging procedure, are in Gradient Checking and a Systematic Debugging Method for Neural Networks.

Verification strategy

Start with a deterministic tiny network. Assert every activation, parameter, and gradient shape. Confirm that the cost and parameter norms remain finite. Compare selected analytical gradients with central-difference estimates. Then overfit a small batch before attempting a full dataset. Together these checks help narrow down whether a failure comes from the calculation, the data, or the optimization settings.

Here we check all 20 parameters of the three-layer example. The perturbations must keep every hidden ReLU on the same side of zero, which the mask assertions verify. Each parameter is restored before testing the next. Once all gradients have been calculated at the original parameter values, one gradient-descent update changes them together.

h = 1e-5
passed, checked = True, 0
for name, values in params.items():
    numerical = np.zeros_like(values)
    for index in np.ndindex(values.shape):
        original = values[index]
        values[index] = original + h
        _, plus = model_forward(X, params)
        loss_plus = mean_loss(plus, Y)
        values[index] = original - h
        _, minus = model_forward(X, params)
        loss_minus = mean_loss(minus, Y)
        values[index] = original
        for base, up, down in zip(caches[:-1], plus[:-1], minus[:-1]):
            assert np.array_equal(base[3] > 0, up[3] > 0)
            assert np.array_equal(base[3] > 0, down[3] > 0)
        numerical[index] = (loss_plus - loss_minus) / (2 * h)
        checked += 1
    passed = passed and np.allclose(numerical, grads["d" + name],
                                    rtol=1e-6, atol=1e-8)
print("gradient check", passed, "parameters", checked)
loss_before = mean_loss(caches, Y)
learning_rate = 0.01
for name in params:
    params[name] -= learning_rate * grads["d" + name]
_, updated_caches = model_forward(X, params)
print("one update", round(loss_before, 6), round(mean_loss(updated_caches, Y), 6))
# gradient check True parameters 20
# one update 1.377523 1.341426

The check agrees within the stated tolerances at these inputs, and the single update lowers this batch’s loss. Neither result establishes convergence or predictive quality. The layer loops reuse one set of local rules; extending the architecture still requires checking activation scales, numerical behavior, and performance on the task.

The multiclass output layer and its fused gradient are covered in Softmax and Multiclass Classification Explained.