Building a Shallow Neural Network with NumPy

The shallow network built here has one hidden layer between its input and output. That small change lets the model learn nonlinear decision boundaries, provided that the hidden layer uses a nonlinear activation. This guide builds the complete model with NumPy and keeps every matrix shape visible.

Architecture and dimensions

Each column of \(X\in\mathbb{R}^{n_x\times m}\) is one example: \(n_x\) counts input features and \(m\) counts examples. We use \(n_h\) hidden units and one sigmoid output for binary classification. The parameters are \(W^{[1]}\in\mathbb{R}^{n_h\times n_x}\), \(b^{[1]}\in\mathbb{R}^{n_h\times1}\), \(W^{[2]}\in\mathbb{R}^{1\times n_h}\), and \(b^{[2]}\in\mathbb{R}^{1\times1}\). Superscripts identify layers. Run the body code blocks in order; later blocks reuse the functions defined earlier.

The network contains \(n_hn_x+n_h+n_h+1\) trainable scalar parameters. The batch size changes activation shapes, but it does not change the parameter count.

The chain-rule machinery behind every gradient here is derived in Derivatives and Computation Graphs for Neural Network Learning.

Why nonlinear activation is necessary

Without a nonlinear hidden activation, the two affine transformations collapse into one affine transformation, even if a sigmoid follows at the output. Here we use tanh in the hidden layer because it provides both positive and negative outputs, and sigmoid to produce a binary probability estimate. Tanh is symmetric about zero, but its observed activation mean depends on its inputs and need not be zero. The other activations in the table are references for later models; this implementation uses only tanh and sigmoid.

The table uses \(a\) for the activation output and assumes \(0<\alpha<1\) for Leaky ReLU. ReLU and Leaky ReLU have no ordinary derivative at zero; the stated values there are backward-pass conventions used below. For GELU, \(\Phi\) and \(\phi\) are the standard normal cumulative distribution function and density.

How activations affect gradients

For an elementwise activation, backpropagation multiplies the incoming gradient by the local derivative. The variable dA below is that incoming gradient, with the same shape as the activation. Small or zero local derivatives can limit what passes through this operation; weights and sums elsewhere in the network also affect the full gradient.

Saturation. The sigmoid derivative \(a(1-a)\) reaches its maximum of 0.25 at \(z=0\) and falls below 0.01 once \(|z|>4.6\). Tanh peaks at 1 but also approaches zero in the tails. Products of small activation derivatives contribute to vanishing gradients. These local properties alone do not establish which activation will train a model faster; input scales, initialization, and the optimizer also matter.

Inactive ReLU units. For \(z<0\), ReLU blocks the incoming gradient through that activation. If a unit is inactive for every example in a batch, that batch supplies no data-loss gradient through the unit to its incoming weights. This does not prove permanent inactivity: other inputs or changes in earlier layers can change its pre-activation. Leaky ReLU uses a small positive slope in the negative region, which permits a gradient to pass when the incoming gradient is nonzero.

This implementation caches activations so their derivatives can reuse forward-pass results. Caching avoids recomputation here; it is an implementation choice, not a mathematical requirement.

import numpy as np

def tanh_backward(dA, A):
    return dA * (1 - A ** 2)          # uses the cached activation

def relu_backward(dA, Z):
    return dA * (Z > 0)               # uses the cached pre-activation

def leaky_relu_backward(dA, Z, alpha=0.01):
    return dA * np.where(Z > 0, 1.0, alpha)

Z = np.array([[-2.0, 0.0, 3.0]])
A = np.tanh(Z)
print(tanh_backward(np.ones_like(A), A))        # [[0.07065082 1.  0.00986604]]
print(relu_backward(np.ones_like(Z), Z))        # [[0. 0. 1.]]
print(leaky_relu_backward(np.ones_like(Z), Z))  # [[0.01 0.01 1.  ]]

At \(z=3\), an incoming gradient of 1 leaves the tanh operation as approximately 0.0099. Three such activation factors have product about \(9.6\times10^{-7}\). This is only their contribution: intervening weight matrices and sums over paths determine the full network gradient.

Initialization and symmetry

With all weights and biases initialized to zero in this tanh network, hidden units remain identical under ordinary gradient descent. More specifically, the hidden activations are zero, so the outgoing weight gradients are zero; the zero outgoing weights also block gradients to the hidden parameters. Only the output bias can change, giving an input-independent prediction. On the balanced XOR labels used below, even that bias has zero gradient at initialization. Setting only W1 to zero while keeping distinct random entries in W2 does not generally give identical hidden gradients, because backpropagation multiplies by W2.T. The initialization below randomizes both weight matrices and scales each by \(1/\sqrt{\text{fan-in}}\), where fan-in is the number of inputs to a unit. Biases start at zero.

def initialize(n_x, n_h, seed=7):
    rng = np.random.default_rng(seed)
    W1 = rng.standard_normal((n_h, n_x)) * np.sqrt(1 / n_x)
    b1 = np.zeros((n_h, 1))
    W2 = rng.standard_normal((1, n_h)) * np.sqrt(1 / n_h)
    b2 = np.zeros((1, 1))
    return {"W1": W1, "b1": b1, "W2": W2, "b2": b2}

Forward propagation

The hidden layer computes \(Z^{[1]}=W^{[1]}X+b^{[1]}\) and \(A^{[1]}=\tanh(Z^{[1]})\), both shaped (n_h, m). The output computes \(Z^{[2]}=W^{[2]}A^{[1]}+b^{[2]}\) and \(A^{[2]}=\sigma(Z^{[2]})\), shaped (1, m). In code, p is the parameter dictionary; the returned dictionary is the forward cache, later called c. The two-branch sigmoid avoids overflow for large negative logits.

def sigmoid(Z):
    A = np.empty_like(Z, dtype=float)
    positive = Z >= 0
    A[positive] = 1 / (1 + np.exp(-Z[positive]))
    e = np.exp(Z[~positive])
    A[~positive] = e / (1 + e)
    return A

def forward(X, p):
    assert X.ndim == 2
    assert p["W1"].ndim == 2 and p["W1"].shape[1] == X.shape[0]
    n_h = p["W1"].shape[0]
    assert p["b1"].shape == (n_h, 1)
    assert p["W2"].shape == (1, n_h) and p["b2"].shape == (1, 1)
    Z1 = p["W1"] @ X + p["b1"]
    A1 = np.tanh(Z1)
    Z2 = p["W2"] @ A1 + p["b2"]
    A2 = sigmoid(Z2)
    return A2, {"Z1": Z1, "A1": A1, "Z2": Z2, "A2": A2}

The loss being differentiated

Labels have shape (1, m) and values 0 or 1. For probabilities \(a_i\), the mean binary cross-entropy is \(J=-\frac1m\sum_i[y_i\log a_i+(1-y_i)\log(1-a_i)]\). We evaluate the same objective from logits to avoid taking logs of rounded probabilities. For each logit \(z\), the stable expression is \(\max(z,0)-zy+\log(1+e^{-|z|})\), whose exponential never receives a positive argument. This protects the calculation for large finite logits; it does not repair nonfinite model values.

def cost_from_logits(Z, Y):
    assert Z.shape == Y.shape
    assert np.all((Y == 0) | (Y == 1))
    return float(np.mean(np.maximum(Z, 0) - Z * Y
                         + np.log1p(np.exp(-np.abs(Z)))))

The shape checks in these functions guard the examples-as-columns convention before arithmetic can silently broadcast the wrong labels or biases. Use a fresh forward cache from the same inputs and parameter values when calling backward.

Backpropagation

For the mean binary cross-entropy defined above, dZ2 = (A2 - Y) / m is the loss derivative with respect to the output logits. The division by \(m\) accounts for averaging over examples. Then \(dW^{[2]}=dZ^{[2]}(A^{[1]})^T\), and summing dZ2 over example columns gives the bias gradient with shape (1, 1).

To reach the hidden layer, W2.T @ dZ2 distributes the output gradient to hidden units, producing shape (n_h, m). Multiplying elementwise by 1 - A1**2 applies each unit’s tanh derivative. Finally, dZ1 @ X.T combines those contributions with the input features into dW1, shaped (n_h, n_x). Summing dZ1 across examples, with keepdims=True, gives db1 shaped (n_h, 1). These gradients already contain the averaging factor inherited from dZ2; no second division is needed.

def backward(X, Y, p, c):
    m = X.shape[1]
    assert Y.shape == c["A2"].shape == (1, m)
    dZ2 = (c["A2"] - Y) / m
    dW2 = dZ2 @ c["A1"].T
    db2 = np.sum(dZ2, axis=1, keepdims=True)
    dZ1 = (p["W2"].T @ dZ2) * (1 - c["A1"] ** 2)
    dW1 = dZ1 @ X.T
    db1 = np.sum(dZ1, axis=1, keepdims=True)
    return {"dW1": dW1, "db1": db1, "dW2": dW2, "db2": db2}

Checking the implemented gradients

Before training, compare the actual backward result with central differences of the actual loss. This three-example batch gives hidden activations shape (4, 3), so the unit count and batch size differ. The check perturbs each of the 17 parameters in turn while keeping all the others fixed.

X_check = np.array([[0.2, -0.4, 0.7], [0.5, 0.3, -0.2]])
Y_check = np.array([[0., 1., 1.]])
p_check = initialize(n_x=2, n_h=4, seed=7)
_, cache_check = forward(X_check, p_check)
grads_check = backward(X_check, Y_check, p_check, cache_check)
h = 1e-5
checked = 0
for name, values in p_check.items():
    numerical = np.zeros_like(values)
    for index in np.ndindex(values.shape):
        original = values[index]
        values[index] = original + h
        _, plus = forward(X_check, p_check)
        values[index] = original - h
        _, minus = forward(X_check, p_check)
        values[index] = original
        numerical[index] = (
            cost_from_logits(plus["Z2"], Y_check)
            - cost_from_logits(minus["Z2"], Y_check)
        ) / (2 * h)
        checked += 1
    assert numerical.shape == grads_check["d" + name].shape
    print(name, np.allclose(numerical, grads_check["d" + name],
                            rtol=1e-6, atol=1e-8))
print("parameters checked", checked)
# W1 True
# b1 True
# W2 True
# b2 True
# parameters checked 17

Each parameter is restored before checking the next one. The four True results mean the returned gradients agree with central differences within the stated tolerances on this batch. This is a value check of the implemented forward, loss, and backward functions; it does not establish correctness at every input. The XOR training run initializes its own parameter dictionary; it does not reuse the check parameters.

Training and diagnosis

Forward propagation computes predictions with the current parameters. Backpropagation calculates how the loss changes with each weight and bias, and gradient descent subtracts a learning-rate-scaled gradient. Below, XOR labels are 1 when the two binary inputs differ and 0 when they match. No single straight decision boundary separates these four points. The model has two input features, four hidden units, and 17 parameters.

def update(p, gradients, learning_rate):
    for name in p:
        p[name] -= learning_rate * gradients["d" + name]

X = np.array([[0., 0., 1., 1.], [0., 1., 0., 1.]])
Y = np.array([[0., 1., 1., 0.]])
p = initialize(n_x=2, n_h=4, seed=7)
learning_rate = 0.1
for step in range(10001):
    A2, c = forward(X, p)
    cost = cost_from_logits(c["Z2"], Y)
    if step in (0, 1000, 10000):
        print("step", step, "cost", round(cost, 6))
    if step == 10000:
        break
    gradients = backward(X, Y, p, c)
    if step == 0:
        assert X.shape == (2, 4) and Y.shape == A2.shape == (1, 4)
        assert c["A1"].shape == (4, 4)
        for name in p:
            assert gradients["d" + name].shape == p[name].shape
    assert np.isfinite(cost)
    assert all(np.isfinite(v).all() for v in gradients.values())
    update(p, gradients, learning_rate)
    assert all(np.isfinite(v).all() for v in p.values())
print("output shape", A2.shape)
print("probabilities", np.round(A2, 4))
print("predictions", (A2 >= 0.5).astype(int))
print("training accuracy", float(np.mean((A2 >= 0.5) == Y)))
# step 0 cost 0.694948
# step 1000 cost 0.173361
# step 10000 cost 0.002566
# output shape (1, 4)
# probabilities [[0.0026 0.9963 0.9995 0.0034]]
# predictions [[0 1 1 0]]
# training accuracy 1.0

The final output is evaluated after 10,000 updates. At threshold 0.5 it matches all four XOR labels. The recorded loss falls during this run; that observation does not guarantee convergence for other seeds, learning rates, or datasets.

Track the cost at intervals rather than printing every iteration. A decreasing training cost confirms that optimization is making progress, but it does not show that the model will work on new examples. Evaluation on suitable data not used for fitting is needed to assess predictive generalization, and a small training subset is useful when checking whether the implementation can learn at all.

• Assert every parameter and gradient shape.
• Confirm that hidden units do not start with identical weights.
• Check that all values remain finite.
• Compare analytical gradients with numerical gradients on a tiny model.
• Try to overfit a very small dataset before scaling up.
• Keep preprocessing identical for training and evaluation data.

The L-layer forward and backward passes are written out in Deep Neural Networks: Architecture and Backpropagation.

Checking the training loop

Separate functions make it easier to locate an error. The forward function can be checked against a hand calculation, and the backward function against finite differences of the loss. An incorrect sign in the update code can reverse the intended direction. A positive learning rate that is too large keeps the negative-gradient direction but can still increase the loss by overshooting. Check the gradient, update sign, and step size separately.

The XOR run above checks whether these functions can fit a small nonlinearly separable dataset. All four possible binary input pairs are used for fitting, so this run exhausts the Boolean XOR domain. It provides no evaluation of an unseen region or noisy observations. To investigate such behavior, first define a larger data distribution, reserve evaluation data from it, and compare with a logistic-regression baseline.

When the model fails, change one assumption at a time. First confirm the data layout and labels. Next check parameter shapes and a single forward result. Then verify the cost and analytical gradients. Only after these checks should learning rate, hidden-unit count, or iteration count be tuned. Changing everything simultaneously can hide the original defect.

A one-hidden-layer network is small enough to inspect completely, yet it contains the essential mechanics of deeper models: parameter initialization, nonlinear representation, cached forward values, chain-rule gradients, and iterative updates.

How initialization affects signal and gradient scales is explained in Weight Initialization and Gradient Flow in Deep Networks.

Exercises

1. Which activation saturates first. For \(z\) in \([-6,6]\), compute the derivative of sigmoid and of tanh, and find the value of \(|z|\) at which each drops below \(0.01\). Which derivative stays above this chosen cutoff over a wider range?

You should get: two thresholds, with tanh’s roughly 1.6 units smaller, plus a four-fold difference in peak value.

import numpy as np
z = np.linspace(-6, 6, 12001)
sig = 1/(1+np.exp(-z)); dsig = sig*(1-sig)
dtanh = 1 - np.tanh(z)**2
for name, d in (("sigmoid", dsig), ("tanh", dtanh)):
    thr = z[(z >= 0) & (d < 0.01)][0]
    print(name, "peak", round(d.max(), 4), "| <0.01 beyond |z| =", round(thr, 3))
# sigmoid peak 0.25 | <0.01 beyond |z| = 4.585
# tanh peak 1.0 | <0.01 beyond |z| = 2.994

2. Batch inactivity under ReLU. Build a small ReLU layer, push its bias strongly negative, and measure the fraction of units that output zero for every input in a batch. Then swap in Leaky ReLU with \(\alpha=0.01\) and measure the fraction of activation entries whose local derivative is zero.

You should get: all units inactive on this batch, with zero local derivatives for ReLU and nonzero local derivatives for Leaky ReLU.

import numpy as np
rng = np.random.default_rng(0)
X = rng.normal(size=(64, 100))
W = rng.normal(size=(32, 64)) * 0.1
b = np.full((32, 1), -5.0)                  # negative bias for this batch
Z = W @ X + b
print("ReLU units always zero:", float(np.mean((Z <= 0).all(axis=1))))
print("ReLU zero local-derivative entries:", float(np.mean(Z <= 0)))
print("LeakyReLU zero local-derivative entries:", float(np.mean(np.where(Z > 0, 1.0, 0.01) == 0)))

3. Gradient signs for one example and a batch. Consider two weights with a single-example gradient \(\delta a\), where both components of \(a\) are positive and \(\delta\ne0\). What signs can the components have? Does the same restriction hold after averaging gradients from several examples?

You should get: matching signs for one example, but potentially mixed signs for a batch.

References

• Nair and Hinton (2010). Rectified Linear Units Improve Restricted Boltzmann Machines. ICML.
• Maas, Hannun, and Ng (2013). Rectifier Nonlinearities Improve Neural Network Acoustic Models.
• Hendrycks and Gimpel (2016). Gaussian Error Linear Units (GELUs). arXiv preprint.