Dev Vaibhav

Posted on Oct 19 Originally published at Medium on Oct 19

Building a Neural Network From Scratch: The Mathematics Behind Learning XOR

#machinelearning #artificialintelligen #deeplearning #python

Understanding forward propagation, backpropagation, and gradient descent through a neural network built using Python, NumPy, and mathematics.

Introduction

Making a neural network produce predictions is relatively straightforward when using modern deep-learning frameworks. Understanding how the network actually learns is a different challenge.

How does it calculate its output? How does it measure an incorrect prediction? And how do derivatives determine which weights and biases need to change?

I wanted to understand these operations without hiding them behind automatic differentiation.

So I built a fully trainable neural network from scratch using Python and NumPy, implementing the mathematical operations myself without relying on PyTorch or TensorFlow.

The network follows a 2 → 3 → 1 architecture and learns the XOR function.

Despite its small size, the implementation brings together several fundamental concepts: matrix multiplication, activation functions, loss functions, the chain rule, backpropagation, and gradient descent.

In this article, I’ll break down these concepts, explain the mathematics behind them, and connect the equations to their actual Python implementation.

AI-assisted writing and visuals.

The goal isn’t just to build a working neural network. It’s to understand how it learns.

1. Why XOR?

XOR (exclusive OR) is a logical operation that returns 1 when its two binary inputs are different and 0 when they are the same.

Despite having only four possible input combinations, XOR presents an interesting challenge.

Figure 1 — XOR truth table and input-space visualization. AI-generated illustration.

A single linear classifier cannot correctly separate the two classes using one straight decision boundary.

Consider a neuron computing:

z = w₁x₁ + w₂x₂ + b

Here, x₁ and x₂ represent the input features, w₁ and w₂ are their corresponding weights, and b is the bias.

Applying a threshold to this weighted sum produces a linear decision boundary. However, the positive XOR examples occupy opposite corners of the input space, as do the negative examples.

No single straight line can correctly separate both classes.

XOR is not linearly separable.

To solve this problem, we need a network capable of learning nonlinear relationships. A hidden layer combined with nonlinear activation functions makes this possible.

2. Network Architecture

To solve XOR, I built a small fully connected feedforward neural network with a 2 → 3 → 1 architecture.

The network consists of three layers:

Input layer: Two input features (x₁ and x₂), representing the binary XOR inputs.
Hidden layer: Three neurons using the Rectified Linear Unit (ReLU) activation function.
Output layer: One neuron using sigmoid activation to produce a value between 0 and 1.

Every input connects to each hidden neuron, and all three hidden neurons connect to the output neuron.

Figure 2 — The 2 → 3 → 1 neural network architecture, showing two inputs, three hidden neurons, and one output. AI-generated illustration.

Each connection has a trainable weight, while every hidden neuron and the output neuron have their own bias.

In total, the network contains 13 trainable parameters :

6 hidden-layer weights
3 hidden-layer biases
3 output-layer weights
1 output bias

These parameters determine how the network transforms inputs into predictions. During training, their values are gradually adjusted to minimize prediction error.

The hidden layer introduces nonlinearity through ReLU, allowing the network to learn patterns that a single linear classifier cannot represent.

The output layer uses sigmoid because we’re solving a binary classification problem. Its output is interpreted as a prediction score, with 0.5 used as the classification threshold.

Despite its small architecture, this network is sufficient to learn the four XOR combinations.

3. Preparing the Dataset and Initializing Parameters

Before implementing forward propagation, we need to prepare our training data and initialize the network’s weights and biases.

Preparing the XOR Dataset

Our dataset contains all four possible XOR input combinations.

import numpy as np

X = np.array([
    [0.0, 0.0],
    [0.0, 1.0],
    [1.0, 0.0],
    [1.0, 1.0]
])

y = np.array([
    [0.0],
    [1.0],
    [1.0],
    [0.0]
])

Here:

X represents the input features, with shape (4, 2).
y represents the expected outputs, with shape (4, 1).

The four rows represent training examples, while the two columns of X represent input features.

Instead of processing each example separately, we’ll use NumPy vectorization to perform calculations across the entire dataset.

Initializing Weights and Biases

Our network contains two sets of weights and biases: one for the hidden layer and another for the output layer.

rng = np.random.default_rng(42)

# Hidden layer
W_h = rng.normal(0, 0.5, size=(2, 3))
b_h = np.zeros(3)

# Output layer
W_o = rng.normal(0, 0.5, size=(3, 1))
b_o = np.zeros(1)

We initialize the weights using a normal distribution with a mean of 0 and a standard deviation of 0.5.

The random seed 42 makes our initialization reproducible, allowing us to obtain the same starting weights across runs.

The biases are initially set to zero.

Why random weights?

If all hidden neurons started with identical weights and biases, they would receive identical gradients and could continue learning the same features.

Random initialization helps break this symmetry, allowing different neurons to learn different representations.

With our dataset and parameters ready, we can begin the first forward pass and observe how the network produces its initial predictions.

4. Forward Propagation: How the Network Makes Predictions

Forward propagation is the process of passing input data through a neural network to calculate its predictions.

At each trainable layer, two operations occur:

Calculate a weighted sum and add a bias.
Apply an activation function.

Let’s break down how these operations work in our 2 → 3 → 1 network.

Figure 3 — Forward propagation through the hidden and output layers. AI-generated illustration.

The diagram uses single-example column-vector notation, while the NumPy implementation processes row-wise batches.

Step 1: Hidden Layer - Weighted Sum

The first operation is:

Zₕ = XWₕ + bₕ

Where:

X: Input matrix with shape (4, 2).
Wₕ: Hidden-layer weights with shape (2, 3).
bₕ: Hidden-layer biases with shape (3,).
Zₕ: Resulting weighted sums with shape (4, 3).

In NumPy:

z_h = X @ W_h + b_h

The @ operator performs matrix multiplication.

The dimensions are:

(4, 2) @ (2, 3) → (4, 3)

NumPy then broadcasts the three biases across all four training examples.

Each row now contains three weighted sums, one for each hidden neuron.

Step 2: Applying ReLU

Next, we introduce nonlinearity using the Rectified Linear Unit (ReLU).

ReLU(z) = max(0, z)

ReLU replaces negative values with zero while keeping positive values unchanged.

h = np.maximum(0, z_h)

For example:

ReLU([-2, 3, -1]) → [0, 3, 0]

Why is this important?

Without nonlinear activation functions, stacking multiple affine layers would still produce an affine transformation. The network would remain unable to learn XOR’s nonlinear decision boundary.

The resulting hidden activation matrix has shape (4, 3).

Step 3: Output Layer - Weighted Sum

The hidden activations now become inputs to the output neuron.

Zₒ = HWₒ + bₒ

Where:

H: Hidden activations with shape (4, 3).
Wₒ: Output weights with shape (3, 1).
bₒ: Output bias with shape (1,).
Zₒ: Output weighted sums with shape (4, 1).

In NumPy:

z_o = h @ W_o + b_o

Matrix multiplication:

(4, 3) @ (3, 1) → (4, 1)

We now have one output value for each training example, before activation.

Step 4: Applying Sigmoid

Finally, we apply the sigmoid activation function.

σ(z) = 1 / (1 + e⁻ᶻ)

Sigmoid maps values into the interval between 0 and 1, making it useful for producing scores in binary classification.

a_o = 1 / (1 + np.exp(-z_o))

We use 0.5 as our classification threshold:

Output ≥ 0.5 → Class 1
Output below 0.5 → Class 0

At this point, our network has successfully completed a forward pass.

However, making predictions is not the same as learning.

The initial outputs depend entirely on the randomly initialized parameters. To improve those predictions, the network must first measure how far they are from the expected values.

Before discussing the loss, let’s follow one actual input through our initialized network.

4.1 A Forward Pass With Actual Numbers

Let’s examine what happens when our randomly initialized network receives the input:

X = [0, 1]

Using our initial weights (seed 42), the hidden-layer weighted sums are approximately:

Zₕ = [0.470282, −0.975518, −0.651090]

Applying ReLU:

H = [0.470282, 0, 0]

Only the first hidden neuron remains active. The other two produce zero because their weighted sums are negative.

Next, we calculate the output weighted sum:

zₒ = 0.470282 × 0.063920 ≈ 0.030061

Applying sigmoid:

ŷ = 1 / (1 + e⁻⁰·⁰³⁰⁰⁶¹) ≈ 0.507515

The correct XOR target for (0, 1) is 1.

Our untrained network produces approximately 0.5075. Although this is already above the classification threshold of 0.5, the prediction still has considerable error.

This is why classification accuracy alone doesn’t tell us how well the network has learned.

We need a loss function.

5. Measuring Error With Mean Squared Error

To measure the difference between predictions and their actual targets, we’ll use Mean Squared Error (MSE).

The formula is:

MSE = (1/N) Σ(ŷᵢ − yᵢ)²

Where:

N represents the number of training examples.
ŷᵢ represents the predicted output.
yᵢ represents the actual target.
Σ represents summation across all examples.

For our network, N = 4 because we have four XOR combinations.

Consider our previous example:

Target = 1

Prediction ≈ 0.507515

The squared error becomes:

(0.507515 − 1)² ≈ 0.242542

We repeat this calculation for all four examples and compute their average.

Figure 4 — Initial XOR predictions, squared errors, and MSE before training. Values reproduced from the actual NumPy implementation.

In NumPy:

error = a_o - y
squared_error = error ** 2
loss = np.mean(squared_error)

At initialization, the network’s actual loss is approximately:

MSE = 0.25023636

Our objective is to reduce this value by adjusting the weights and biases.

However, calculating the loss only tells us how wrong the predictions are. It doesn’t tell us which parameters to change or by how much.

To answer that question, we need derivatives and the chain rule.

That’s where backpropagation begins.

6. Backpropagation: Understanding the Chain Rule

So far, our network can produce predictions and calculate their errors using Mean Squared Error.

But one important question remains:

How does the network know which weights and biases to change to reduce the loss?

This is where backpropagation comes in.

Backpropagation calculates the gradients of the loss with respect to each trainable parameter using the chain rule of differentiation.

A gradient tells us how sensitive the loss is to a small change in a particular parameter.

Before calculating these gradients, we need to understand how the operations depend on one another.

Figure 5 — Gradient flow from the loss through the output and hidden layers. AI-generated illustration.

6.1 Understanding the Chain Rule

Recall the equations from forward propagation:

z = wh + b

a = σ(z)

L = (a − y)²

Here:

h represents an activation entering the neuron.
w represents its weight.
b represents its bias.
z represents the weighted sum.
a represents the activated output.
y represents the expected target.
L represents the squared error for one example.

Notice how these operations depend on one another.

The weight influences the weighted sum, the weighted sum influences the activation, and the activation influences the loss.

w → z → a → L

To calculate how changing the weight affects the loss, we apply the chain rule:

∂L/∂w = (∂L/∂a) × (∂a/∂z) × (∂z/∂w)

This expression contains three derivatives:

∂L/∂a: How the loss changes when the activation changes.
∂a/∂z: How the activation changes when the weighted sum changes.
∂z/∂w: How the weighted sum changes when the weight changes.

Multiplying them gives us the gradient of the loss with respect to the weight.

Our actual network uses four training examples, so we’ll apply this same principle to the batch MSE.

6.2 Differentiating Mean Squared Error

Recall our loss function:

L = (1/N) Σ(aᵢ − yᵢ)²

We want to calculate its derivative with respect to each prediction, aᵢ.

For one example, the squared error is:

(aᵢ − yᵢ)²

Using the power rule and chain rule:

d/daᵢ [(aᵢ − yᵢ)²] = 2(aᵢ − yᵢ)

Since MSE averages the squared errors across N examples, we include the factor 1/N:

∂L/∂aᵢ = (2/N)(aᵢ − yᵢ)

For our complete batch, the NumPy implementation is:

dl_da_o = 2 / len(X) * (a_o - y)

Here, len(X) equals 4, and the resulting gradient array has shape (4, 1).

We have now calculated the first derivative needed for backpropagation.

Next, we’ll differentiate the sigmoid activation and combine both derivatives using the chain rule.

6.3 Differentiating the Sigmoid Activation

From forward propagation, recall the sigmoid function:

a = σ(z) = 1 / (1 + e⁻ᶻ)

Here, z represents the output neuron’s weighted sum, and a represents its activated output.

To understand how changes in z affect a, we differentiate sigmoid with respect to z.

The derivative is:

∂a/∂z = e⁻ᶻ / (1 + e⁻ᶻ)²

This expression can be simplified using the original sigmoid function:

∂a/∂z = σ(z)(1 − σ(z))

Since we already calculated the sigmoid output during forward propagation, we can reuse it directly:

∂a/∂z = a(1 − a)

In NumPy:

da_od_zo = a_o * (1 - a_o)

This calculates the sigmoid derivative for all four training examples simultaneously.

The resulting array has shape (4, 1).

6.4 Applying the Chain Rule

We now have two derivatives.

From Section 6.2, the MSE derivative:

∂L/∂a = (2/N)(a − y)

And from Section 6.3, the sigmoid derivative:

∂a/∂z = a(1 − a)

We want to calculate how the loss changes with respect to the output weighted sum, z.

Using the chain rule:

∂L/∂z = (∂L/∂a) × (∂a/∂z)

Substituting both expressions:

∂L/∂z = (2/N)(a − y) × a(1 − a)

For the entire batch, multiplication is performed elementwise.

In NumPy:

dl_dz_o = dl_da_o * da_od_zo

The resulting gradient has shape (4, 1).

We have now propagated the loss gradient backward through the sigmoid activation.

The next step is to determine how the output weights and bias contribute to this gradient.

6.5 Calculating Output Weight and Bias Gradients

Recall our output-layer equation:

Zₒ = HWₒ + bₒ

Here:

H: Hidden-layer activations.
Wₒ: Output-layer weights.
bₒ: Output-layer bias.
Zₒ: Output pre-activation values.

For a single training example, the output weighted sum is:

zₒ = h₁w₁ + h₂w₂ + h₃w₃ + bₒ

Let’s differentiate it with respect to the first output weight, w₁.

Since the other terms do not depend on w₁:

∂zₒ/∂w₁ = h₁

Similarly:

∂zₒ/∂w₂ = h₂

∂zₒ/∂w₃ = h₃

Now apply the chain rule:

∂L/∂w₁ = (∂L/∂zₒ) × h₁

The same principle applies to the remaining output weights.

Because we’re processing four examples, we sum their contributions using matrix multiplication.

∂L/∂Wₒ = Hᵀ @ (∂L/∂Zₒ)

In NumPy:

dW_o = h.T @ dl_dz_o

Let’s verify the matrix dimensions:

(3, 4) @ (4, 1) → (3, 1)

The result contains three gradients, one for each output weight.

Calculating the Output Bias Gradient

The output weighted sum also depends on its bias:

zₒ = h₁w₁ + h₂w₂ + h₃w₃ + bₒ

Differentiating with respect to bₒ:

∂zₒ/∂bₒ = 1

Applying the chain rule:

∂L/∂bₒ = (∂L/∂zₒ) × 1

For the complete batch, we sum the contributions from all four examples:

db_o = np.sum(dl_dz_o, axis=0)

The result has shape (1,), since our output layer contains only one bias.

Notice that we don’t divide by N again. The averaging factor was already included when differentiating MSE.

At this point, we’ve calculated the gradients for all three output weights and the output bias.

However, our network also contains six hidden-layer weights and three hidden biases.

To calculate their gradients, we need to propagate the error backward through the hidden layer and its ReLU activation.

6.6 Backpropagating Through the Hidden Layer

So far, we’ve calculated the gradients for the output layer: three weights and one bias.

But our network also contains six hidden-layer weights and three hidden biases.

To calculate their gradients, we need to propagate the loss backward through the hidden layer.

Recall the equations from forward propagation:

Zₕ = XWₕ + bₕ

H = ReLU(Zₕ)

Zₒ = HWₒ + bₒ

The hidden activations influence the output weighted sum, which influences the prediction and ultimately the loss.

We already calculated ∂L/∂Zₒ in Section 6.4.

Now we want to determine how the loss changes with respect to the hidden activations, H.

For a single example, the output weighted sum is:

zₒ = h₁w₁ + h₂w₂ + h₃w₃ + bₒ

Differentiating with respect to each hidden activation:

∂zₒ/∂h₁ = w₁

∂zₒ/∂h₂ = w₂

∂zₒ/∂h₃ = w₃

Applying the chain rule:

∂L/∂hⱼ = (∂L/∂zₒ) × wⱼ

For the entire batch, we can calculate these gradients using matrix multiplication:

∂L/∂H = (∂L/∂Zₒ) @ Wₒᵀ

In NumPy:

dl_dh = dl_dz_o @ W_o.T

Let’s verify the dimensions:

(4, 1) @ (1, 3) → (4, 3)

The resulting matrix contains the loss gradients for all three hidden activations across all four training examples.

Notice that we use Wₒᵀ , the transpose of the output weight matrix, not the output weight gradients.

6.7 Differentiating ReLU

Our hidden layer uses ReLU:

h = ReLU(z) = max(0, z)

Its derivative is:

ReLU’(z) = 1 when z > 0, otherwise 0

At exactly zero, ReLU is not differentiable. Our implementation uses a derivative value of 0.

This means gradients pass through positive hidden pre-activations, while inactive ReLU units receive zero gradient through that activation.

We can calculate the derivative using a Boolean mask:

dh_dhz = (z_h > 0).astype(float)

The expression (z_h > 0) returns True for positive values and False otherwise.

Using .astype(float) converts these Boolean values into 1.0 and 0.0.

Now we apply the chain rule:

∂L/∂Zₕ = (∂L/∂H) × ReLU’(Zₕ)

In NumPy:

dl_dzh = dl_dh * dh_dhz

This operation uses elementwise multiplication.

Both arrays have shape (4, 3), so the result also has shape (4, 3).

We have now propagated the loss gradient backward through the hidden activation function.

6.8 Calculating Hidden Weight Gradients

Recall the hidden-layer weighted sum:

Zₕ = XWₕ + bₕ

For a single hidden neuron, its weighted sum is:

zₕ = x₁w₁ + x₂w₂ + bₕ

Differentiating with respect to the first weight:

∂zₕ/∂w₁ = x₁

Similarly:

∂zₕ/∂w₂ = x₂

Applying the chain rule:

∂L/∂w₁ = (∂L/∂zₕ) × x₁

∂L/∂w₂ = (∂L/∂zₕ) × x₂

We apply this principle to all three hidden neurons and sum contributions across the four training examples.

The matrix equation becomes:

∂L/∂Wₕ = Xᵀ @ (∂L/∂Zₕ)

In NumPy:

dW_h = X.T @ dl_dzh

Checking the dimensions:

(2, 4) @ (4, 3) → (2, 3)

The resulting gradient matrix contains six values, corresponding to our six hidden-layer weights.

6.9 Calculating Hidden Bias Gradients

Finally, we calculate the hidden bias gradients.

For a hidden neuron:

zₕ = x₁w₁ + x₂w₂ + bₕ

Differentiating with respect to its bias:

∂zₕ/∂bₕ = 1

Applying the chain rule:

∂L/∂bₕ = (∂L/∂zₕ) × 1

For all four training examples, we sum the gradients along the batch dimension:

db_h = np.sum(dl_dzh, axis=0)

The result has shape (3,), giving us one gradient for each hidden neuron.

6.10 Completing Backpropagation

We’ve now calculated all the gradients required to train our network.

Gradient Shape Parameters
----------------------------------
dW_h (2, 3) 6
db_h (3,) 3
dW_o (3, 1) 3
db_o (1,) 1
----------------------------------
Total 13

The complete hidden-layer backpropagation implementation is:

# Gradient through output weights
dl_dh = dl_dz_o @ W_o.T

# Gradient through ReLU
dh_dhz = (z_h > 0).astype(float)
dl_dzh = dl_dh * dh_dhz

# Hidden weight gradients
dW_h = X.T @ dl_dzh

# Hidden bias gradients
db_h = np.sum(dl_dzh, axis=0)

All these gradients are calculated using the current weights and the values saved during forward propagation.

An important distinction: backpropagation calculates gradients; it does not update the network’s parameters.

We now know how each of our 13 parameters influences the loss.

The next step is to use these gradients to adjust the weights and biases through gradient descent.

7. Gradient Descent: Updating the Parameters

We’ve now calculated the gradients for all 13 trainable parameters through backpropagation.

But calculating gradients alone doesn’t improve our network. We need to use them to update the weights and biases.

This is where gradient descent comes in.

Gradient descent is an optimization algorithm that adjusts parameters in the direction opposite to their gradients, aiming to minimize the loss function.

7.1 Understanding the Update Rule

The general gradient descent equation is:

θ_new = θ_old − α(∂L/∂θ)

Where:

θ represents a trainable parameter, such as a weight or bias.
α represents the learning rate.
∂L/∂θ represents the gradient of the loss with respect to that parameter.

The learning rate controls the size of each update.

A very large learning rate may cause unstable training, while a very small learning rate can make convergence slow.

For our XOR network, we use:

α = 0.05

Figure 6 — Gradient descent updates weights and biases using the calculated gradients. AI-generated illustration.

In this diagram, η represents the learning rate, denoted by α in our equations. Both refer to the same parameter.

7.2 Updating Weights and Biases

We apply the same update rule to all four parameter groups.

In NumPy:

# Hidden layer updates
W_h = W_h - A * dW_h
b_h = b_h - A * db_h

# Output layer updates
W_o = W_o - A * dW_o
b_o = b_o - A * db_o

Here, A stores our learning rate of 0.05.

Each parameter is adjusted using its corresponding gradient.

An important detail is that all gradients must be calculated before updating any weights or biases.

This ensures that every gradient corresponds to the same forward pass and parameter state.

7.3 A Real Weight Update

Let’s examine one actual update from our network’s initial training step.

The first output-layer weight begins at approximately:

w = 0.06392020

After backpropagation, its gradient is:

∂L/∂w ≈ 0.00123833

Using the gradient descent equation:

w_new = 0.06392020 − (0.05 × 0.00123833)

Therefore:

w_new ≈ 0.06385829

The weight changes only slightly, but this is how learning progresses: through repeated updates across all trainable parameters.

After updating all 13 parameters in our first training step, the MSE decreases from approximately 0.25023636 to 0.25020869.

This is a small improvement, but it demonstrates that the calculated gradients are being used to adjust the network.

7.4 Repeating the Learning Process

One update is not enough to learn the complete XOR mapping.

We repeat the training cycle:

Forward Propagation → Loss Calculation → Backpropagation → Gradient Descent

After each update, the new weights and biases are used to calculate the next predictions.

Over many iterations, the network gradually adjusts its parameters to reduce prediction errors.

In our implementation, we repeat this process for 5,000 epochs , using all four XOR examples for each update.

This method is called full-batch gradient descent.

Now that we’ve implemented forward propagation, loss calculation, backpropagation, and gradient descent, we can combine everything into one complete training loop.

8. Bringing Everything Together: The Training Loop

We’ve implemented all four fundamental operations required to train our neural network:

Forward propagation — Calculate predictions using the current weights and biases.
Loss calculation — Measure prediction errors using Mean Squared Error.
Backpropagation — Calculate gradients using the chain rule.
Gradient descent — Update weights and biases to reduce the loss.

Now we can combine these operations into a single training loop.

8.1 Implementing the Complete Training Loop

We use a learning rate of 0.05 and train the network for 5,000 epochs.

A = 0.05
epochs = 5000

for epoch in range(epochs):

    # Forward propagation
    z_h = X @ W_h + b_h
    h = np.maximum(0, z_h)

    z_o = h @ W_o + b_o
    a_o = 1 / (1 + np.exp(-z_o))

    # Mean Squared Error
    error = a_o - y
    squared_error = error ** 2
    loss = np.mean(squared_error)

    # Output layer backpropagation
    dl_da_o = 2 / len(X) * (a_o - y)
    da_od_zo = a_o * (1 - a_o)
    dl_dz_o = dl_da_o * da_od_zo

    dW_o = h.T @ dl_dz_o
    db_o = np.sum(dl_dz_o, axis=0)

    # Hidden layer backpropagation
    dl_dh = dl_dz_o @ W_o.T

    dh_dhz = (z_h > 0).astype(float)
    dl_dzh = dl_dh * dh_dhz

    dW_h = X.T @ dl_dzh
    db_h = np.sum(dl_dzh, axis=0)

    # Gradient descent
    W_h = W_h - A * dW_h
    b_h = b_h - A * db_h

    W_o = W_o - A * dW_o
    b_o = b_o - A * db_o

    # Display training progress
    if epoch % 500 == 0:
        print(f"Epoch {epoch} | Loss: {loss:.6f}")

This is the complete training loop from our NumPy implementation.

Notice that every operation we previously derived appears here. There is no automatic differentiation or built-in neural network optimizer.

8.2 Understanding an Epoch

An epoch represents one complete pass through the training dataset.

Since our dataset contains only four XOR combinations, each epoch processes all four examples together.

During each epoch, the network:

Calculates the predictions for all four examples.
Computes the average squared error.
Calculates gradients for all 13 parameters.
Updates the weights and biases once.

Because the complete dataset contributes to every gradient calculation, this approach is called full-batch gradient descent.

8.3 Tracking the Training Loss

We use the following condition to print the training progress:

if epoch % 500 == 0:
    print(f"Epoch {epoch} | Loss: {loss:.6f}")

This prints the loss every 500 epochs, allowing us to observe how the network changes during training.

The training loop still performs an update during every epoch. The condition only controls when the loss is displayed.

An important detail is that the displayed loss is calculated before the parameter update for that epoch.

After training finishes, we need to perform another forward pass using the final updated weights and biases to obtain the actual final predictions and loss.

8.4 Why the Order Matters

The sequence of operations is important.

All gradients must be calculated before updating the parameters because backpropagation depends on the weights and activations from the same forward pass.

Updating the output weights before calculating the hidden-layer gradients would introduce inconsistencies.

By following the correct sequence, every gradient corresponds to the same parameter state.

After 5,000 epochs, our network has repeatedly adjusted its parameters to minimize prediction errors.

The next step is to evaluate the trained network and examine whether it has successfully learned the XOR mapping.

9. Evaluating the Trained Neural Network

After 5,000 epochs of forward propagation, backpropagation, and gradient descent, it’s time to evaluate whether our network has successfully learned the XOR function.

We’ll examine how the loss changed during training, inspect the final predictions, and calculate classification accuracy.

9.1 Evaluating the Final Parameters

Once training finishes, we perform another forward pass using the updated weights and biases.

z_h_new = X @ W_h + b_h
h_new = np.maximum(0, z_h_new)

z_o_new = h_new @ W_o + b_o
a_o_new = 1 / (1 + np.exp(-z_o_new))

# Final Mean Squared Error
new_loss = np.mean((a_o_new - y) ** 2)

# Convert probabilities into binary classes
predictions = (a_o_new >= 0.5).astype(int)

Why perform another forward pass?

The predictions calculated during the last training epoch were produced before the final parameter update.

To evaluate the fully trained network, we need predictions calculated using its final weights and biases.

9.2 Visualizing the Training Loss

At the beginning of training, our network’s Mean Squared Error was approximately:

Initial MSE = 0.25023636

As training progressed, the loss gradually decreased.

Actual loss across 5,000 training updates, reproduced from our NumPy experiment.

Here are some recorded loss values:

Epoch Training MSE
-------------------------
0 0.25023636
500 0.22134667
1000 0.18205232
1500 0.16114958
2000 0.08436386
2500 0.02993619
3000 0.01605965
3500 0.01045661
4000 0.00756720
4500 0.00584314

These values represent the loss before the parameter update at each listed epoch.

The loss decreases gradually at first, followed by a more noticeable improvement around epochs 1,500–2,500.

This indicates that the network is progressively finding parameters that fit the XOR targets more closely.

9.3 Final XOR Predictions

After completing all 5,000 updates, the network produces the following results:

Input Target Output Class
-------------------------------------
(0,0) 0 0.08389788 0
(0,1) 1 0.95051648 1
(1,0) 1 0.95173893 1
(1,1) 0 0.08393200 0

Final sigmoid outputs compared with the expected XOR targets.

The model outputs values close to zero for the negative XOR examples and values close to one for the positive examples.

Using a classification threshold of 0.5, all four predictions match their expected targets.

9.4 Final Loss and Accuracy

Our final Mean Squared Error is:

MSE = 0.00471529585312808

And our classification accuracy is:

Accuracy = (4 / 4) × 100 = 100%

The network has successfully learned to classify all four XOR combinations in its training dataset.

However, an important distinction needs to be made.

100% training accuracy does not necessarily imply strong generalization.

Our evaluation uses the same four combinations that were used during training. This experiment demonstrates that our small network can fit the XOR mapping, but it does not evaluate performance on an independent dataset.

Additionally, training behaviour can vary with random initialization, learning rate, and network architecture.

For this particular configuration, the results demonstrate that our manually implemented forward propagation, backpropagation, and gradient descent work together to learn a nonlinear mapping without relying on a deep-learning framework.

10. What I Learned From Building This Network

Building this neural network from scratch helped me understand how the mathematical concepts behind neural networks connect to one another.

Instead of relying on a framework to handle the calculations, I had to work through every step myself.

A few things became much clearer during the process:

Forward propagation determines how inputs are transformed into predictions through weights, biases, and activation functions.
Loss functions measure how far predictions are from their expected targets.
Backpropagation applies the chain rule to calculate how each trainable parameter influences the loss.
Gradient descent uses those gradients to update parameters and improve predictions.
Matrix dimensions and vectorization are essential for implementing these operations efficiently using NumPy.

One particularly important realization was that backpropagation and gradient descent are not the same operation.

Backpropagation calculates gradients. Gradient descent uses them to update the parameters.

Understanding this distinction made the training process much easier to reason about.

10.1 Challenges and Limitations

The most challenging part was understanding how gradients propagate backward through multiple layers.

Calculating derivatives for a single neuron was manageable, but extending those calculations across an entire batch required careful attention to the chain rule, matrix multiplication, transposes, and array shapes.

The project also has a few limitations:

Small dataset: The network was trained and evaluated on the same four XOR combinations.
Initialization sensitivity: Different random weights can lead to different training outcomes, particularly with a small ReLU network.
Limited evaluation: Achieving 100% training accuracy does not demonstrate generalization to unseen data.
Gradient verification: Although the training results were reproduced, a complete numerical gradient check of all 13 parameters remains future work.

These limitations don’t invalidate the experiment. They define what the results demonstrate and what still needs further investigation.

11. Conclusion

This project started with a simple question:

How does a neural network actually learn?

Building a small XOR network helped me understand the answer through mathematics and implementation rather than relying entirely on abstractions.

Using Python and NumPy, I implemented a 2 → 3 → 1 neural network with 13 trainable parameters, ReLU and sigmoid activations, Mean Squared Error, manually derived backpropagation, and full-batch gradient descent.

After 5,000 epochs, the network successfully classified all four XOR combinations, achieving a final training MSE of approximately 0.004715.

The architecture is small, but the fundamental principles involved extend to much larger neural networks.

There’s still plenty to explore, from numerical gradient checking and alternative activation functions to deeper architectures and more challenging datasets.

But this experiment established an important foundation for understanding those systems.

Frameworks make neural networks easier to build. Understanding the mathematics makes them easier to reason about.

Source Code and Further Reading

The complete implementation, mathematical derivations, training results, and diagrams are available in my GitHub repository:

Neural-Network — GitHub

If you find the project helpful, consider starring the repository. Feedback and suggestions are always welcome.