A Neural Network from Scratch, With Every Matrix Explained

What We Are Building

In this notebook we will build a small neural network in pure NumPy and explain every major tensor along the way.

The goal is not just to get a model that works. The goal is to make the matrix math feel obvious.

By the end, you will have:

• a vectorized 2 -> 4 -> 1 neural network
• a clear mental model for what X, W1, b1, Z1, A1, W2, and A2 actually mean
• an explanation for why we transpose the data matrix
• interactive visualizations for the data, training loss, and decision boundary over time

We will use a two-class make_moons dataset because it is simple to plot but nonlinear enough that a hidden layer actually matters.

Our Matrix Convention

I will use the neural-network convention below for the full post:

• X has shape (n_features, m)
• Y has shape (1, m)
• each column is one training example
• each row is one feature tracked across all examples

That means if we have 2 features and 300 training examples, X.shape == (2, 300).

Why transpose into this layout? Because it makes the forward pass and the gradients line up cleanly:

\[ Z_1 = W_1 X + b_1 \]

If W1 has shape (n_h, n_x) and X has shape (n_x, m), then one matrix multiplication computes the hidden-layer pre-activation for all m examples at once.

import numpy as np
import pandas as pd
import plotly.io as pio
pio.renderers.default = "notebook"

import plotly.express as px
import plotly.graph_objects as go

from IPython.display import display
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

np.random.seed(42)

X_raw, y_raw = make_moons(n_samples=400, noise=0.20, random_state=42)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)

X_train_rows, X_test_rows, y_train, y_test = train_test_split(
    X_scaled,
    y_raw,
    test_size=0.25,
    random_state=42,
    stratify=y_raw,
)

X_train = X_train_rows.T
X_test = X_test_rows.T
Y_train = y_train.reshape(1, -1)
Y_test = y_test.reshape(1, -1)

n_x = X_train.shape[0]
m_train = X_train.shape[1]

print(f"X_train_rows shape (scikit-learn default): {X_train_rows.shape}")
print(f"X_train shape (vectorized neural-network layout): {X_train.shape}")
print(f"Y_train shape: {Y_train.shape}")
print(f"Each column in X_train is one example with {n_x} features.")

X_train_rows shape (scikit-learn default): (300, 2)
X_train shape (vectorized neural-network layout): (2, 300)
Y_train shape: (1, 300)
Each column in X_train is one example with 2 features.

shape_guide = pd.DataFrame(
    [
        ("X_train_rows", X_train_rows.shape, "rows = examples, columns = features", "default shape returned by scikit-learn"),
        ("X_train", X_train.shape, "rows = features, columns = examples", "best layout for vectorized neural-network equations"),
        ("Y_train", Y_train.shape, "one label row across all examples", "binary targets for the full batch"),
    ],
    columns=["object", "shape", "what_it_represents", "why_it_exists"],
)
display(shape_guide)

example_0 = pd.DataFrame(
    {
        "feature": ["x1", "x2"],
        "value_for_example_0": np.round(X_train[:, 0], 3),
        "meaning": [
            "standardized horizontal coordinate",
            "standardized vertical coordinate",
        ],
    }
)
display(example_0)

plot_df = pd.DataFrame(X_train_rows, columns=["x1", "x2"])
plot_df["label"] = y_train.astype(str)

fig = px.scatter(
    plot_df,
    x="x1",
    y="x2",
    color="label",
    title="Training Data in Feature Space",
    opacity=0.85,
    template="plotly_white",
)
fig.update_traces(marker=dict(size=9, line=dict(width=1, color="white")))
fig.update_layout(height=500)
fig.show()

	object	shape	what_it_represents	why_it_exists
0	X_train_rows	(300, 2)	rows = examples, columns = features	default shape returned by scikit-learn
1	X_train	(2, 300)	rows = features, columns = examples	best layout for vectorized neural-network equa...
2	Y_train	(1, 300)	one label row across all examples	binary targets for the full batch

	feature	value_for_example_0	meaning
0	x1	-0.914	standardized horizontal coordinate
1	x2	0.647	standardized vertical coordinate

Forward Pass: What Every Value Means

We will use a 2 -> 4 -> 1 network:

• 2 input features: x1, x2
• 4 hidden neurons
• 1 output neuron for the probability of class 1

The forward pass is:

\[ Z_1 = W_1 X + b_1 \] \[ A_1 = tanh(Z_1) \] \[ Z_2 = W_2 A_1 + b_2 \] \[ A_2 = \sigma(Z_2) \]

Interpretation:

• W1[i, j] tells us how strongly feature j contributes to hidden neuron i
• b1[i, 0] shifts hidden neuron i for every example in the batch
• Z1[i, k] is the raw score of hidden neuron i on example k
• A1[i, k] is the activated hidden representation after tanh
• Z2[0, k] is the final raw score for example k
• A2[0, k] is the predicted probability for example k

Notice the bias shapes: b1 is (4, 1) and b2 is (1, 1). They are not transposed because we want NumPy broadcasting to copy each bias across all m columns.

def sigmoid(z):
    return 1 / (1 + np.exp(-z))


def initialize_parameters(n_x, n_h, n_y, seed=7):
    rng = np.random.default_rng(seed)
    return {
        "W1": rng.normal(0, 0.6, size=(n_h, n_x)),
        "b1": np.zeros((n_h, 1)),
        "W2": rng.normal(0, 0.6, size=(n_y, n_h)),
        "b2": np.zeros((n_y, 1)),
    }


def forward_propagation(X, parameters):
    W1, b1 = parameters["W1"], parameters["b1"]
    W2, b2 = parameters["W2"], parameters["b2"]

    Z1 = W1 @ X + b1
    A1 = np.tanh(Z1)
    Z2 = W2 @ A1 + b2
    A2 = sigmoid(Z2)

    cache = {"Z1": Z1, "A1": A1, "Z2": Z2, "A2": A2}
    return A2, cache


def compute_loss(A2, Y):
    eps = 1e-9
    return -np.mean(Y * np.log(A2 + eps) + (1 - Y) * np.log(1 - A2 + eps))


parameters = initialize_parameters(n_x=n_x, n_h=4, n_y=1, seed=7)
X_demo = X_train[:, :5]
Y_demo = Y_train[:, :5]
A2_demo, cache_demo = forward_propagation(X_demo, parameters)

forward_shapes = pd.DataFrame(
    [
        ("W1", parameters["W1"].shape, "weights from inputs to hidden layer"),
        ("b1", parameters["b1"].shape, "hidden-layer bias"),
        ("X_demo", X_demo.shape, "5 examples stacked as columns"),
        ("Z1", cache_demo["Z1"].shape, "pre-activation hidden scores"),
        ("A1", cache_demo["A1"].shape, "hidden activations after tanh"),
        ("W2", parameters["W2"].shape, "weights from hidden layer to output"),
        ("Z2", cache_demo["Z2"].shape, "output pre-activation"),
        ("A2", cache_demo["A2"].shape, "predicted probabilities"),
    ],
    columns=["object", "shape", "meaning"],
)
display(forward_shapes)

display(
    pd.DataFrame(
        np.round(cache_demo["Z1"], 3),
        index=[f"hidden_{i}" for i in range(1, 5)],
        columns=[f"example_{i}" for i in range(5)],
    )
)

display(
    pd.DataFrame(
        np.round(A2_demo, 3),
        index=["predicted_prob"],
        columns=[f"example_{i}" for i in range(5)],
    )
)

	object	shape	meaning
0	W1	(4, 2)	weights from inputs to hidden layer
1	b1	(4, 1)	hidden-layer bias
2	X_demo	(2, 5)	5 examples stacked as columns
3	Z1	(4, 5)	pre-activation hidden scores
4	A1	(4, 5)	hidden activations after tanh
5	W2	(1, 4)	weights from hidden layer to output
6	Z2	(1, 5)	output pre-activation
7	A2	(1, 5)	predicted probabilities

	example_0	example_1	example_2	example_3	example_4
hidden_1	0.115	0.229	0.137	-0.191	0.289
hidden_2	-0.196	-0.575	-0.141	0.452	-0.637
hidden_3	-0.136	-0.581	-0.009	0.437	-0.584
hidden_4	0.487	1.007	0.560	-0.835	1.253

	example_0	example_1	example_2	example_3	example_4
predicted_prob	0.524	0.534	0.529	0.468	0.538

Why the Transposes Show Up

There are two transposes that data scientists usually ask about first:

1. X_train = X_train_rows.T
2. gradient expressions like dW1 = dZ1 @ X.T / m

The first transpose changes the data layout so that every example becomes a column. That lets one matrix multiplication process the full batch.

The second transpose exists because gradients for the weights are outer products accumulated across examples.

• dZ1 has shape (n_h, m)
• X has shape (n_x, m)
• X.T has shape (m, n_x)
• so dZ1 @ X.T becomes (n_h, n_x), which matches W1

That is the real test for whether your transpose is correct: the resulting gradient must have the exact same shape as the parameter you want to update.

sample_idx = 0
hidden_idx = 0

x_sample = X_demo[:, sample_idx]
w_hidden = parameters["W1"][hidden_idx, :]
b_hidden = parameters["b1"][hidden_idx, 0]

contribution_df = pd.DataFrame(
    {
        "feature": ["x1", "x2"],
        "x_value": np.round(x_sample, 3),
        "weight_into_hidden_1": np.round(w_hidden, 3),
        "weight_times_value": np.round(w_hidden * x_sample, 3),
    }
)
display(contribution_df)

z1_manual = float(w_hidden @ x_sample + b_hidden)
a1_manual = float(np.tanh(z1_manual))
z2_manual = float(parameters["W2"][0, :] @ cache_demo["A1"][:, sample_idx] + parameters["b2"][0, 0])
a2_manual = float(sigmoid(z2_manual))

trace_df = pd.DataFrame(
    [
        ("z1 for hidden neuron 1", round(z1_manual, 4), "weighted sum of x1 and x2 plus bias"),
        ("a1 for hidden neuron 1", round(a1_manual, 4), "tanh activation for that hidden neuron"),
        ("z2 for output neuron", round(z2_manual, 4), "weighted sum across all hidden activations"),
        ("a2 final probability", round(a2_manual, 4), "sigmoid output interpreted as P(class = 1)"),
    ],
    columns=["quantity", "value", "interpretation"],
)
display(trace_df)

	feature	x_value	weight_into_hidden_1	weight_times_value
0	x1	-0.914	0.001	-0.001
1	x2	0.647	0.179	0.116

	quantity	value	interpretation
0	z1 for hidden neuron 1	0.1153	weighted sum of x1 and x2 plus bias
1	a1 for hidden neuron 1	0.1148	tanh activation for that hidden neuron
2	z2 for output neuron	0.0951	weighted sum across all hidden activations
3	a2 final probability	0.5238	sigmoid output interpreted as P(class = 1)

Backpropagation

For binary classification with sigmoid output and log loss, the vectorized gradients are:

\[ dZ_2 = A_2 - Y \] \[ dW_2 = \frac{1}{m} dZ_2 A_1^T \] \[ db_2 = \frac{1}{m} \sum dZ_2 \] \[ dZ_1 = (W_2^T dZ_2) \odot (1 - A_1^2) \] \[ dW_1 = \frac{1}{m} dZ_1 X^T \] \[ db_1 = \frac{1}{m} \sum dZ_1 \]

The transpose in W2.T @ dZ2 moves the output-layer signal back into hidden-layer space. The transpose in A1.T or X.T turns examples into rows so each example contributes one outer-product slice to the weight gradient.

def backward_propagation(parameters, cache, X, Y):
    m = X.shape[1]
    W2 = parameters["W2"]
    A1, A2 = cache["A1"], cache["A2"]

    dZ2 = A2 - Y
    dW2 = (dZ2 @ A1.T) / m
    db2 = np.sum(dZ2, axis=1, keepdims=True) / m

    dZ1 = (W2.T @ dZ2) * (1 - A1 ** 2)
    dW1 = (dZ1 @ X.T) / m
    db1 = np.sum(dZ1, axis=1, keepdims=True) / m

    return {"dW1": dW1, "db1": db1, "dW2": dW2, "db2": db2, "dZ2": dZ2, "dZ1": dZ1}


def update_parameters(parameters, grads, learning_rate):
    return {
        "W1": parameters["W1"] - learning_rate * grads["dW1"],
        "b1": parameters["b1"] - learning_rate * grads["db1"],
        "W2": parameters["W2"] - learning_rate * grads["dW2"],
        "b2": parameters["b2"] - learning_rate * grads["db2"],
    }


grads_demo = backward_propagation(parameters, cache_demo, X_demo, Y_demo)

gradient_shapes = pd.DataFrame(
    [
        ("dZ2", grads_demo["dZ2"].shape, "prediction error at the output layer"),
        ("dW2", grads_demo["dW2"].shape, "gradient for W2; same shape as W2"),
        ("db2", grads_demo["db2"].shape, "gradient for b2; same shape as b2"),
        ("dZ1", grads_demo["dZ1"].shape, "hidden-layer error after backprop"),
        ("dW1", grads_demo["dW1"].shape, "gradient for W1; same shape as W1"),
        ("db1", grads_demo["db1"].shape, "gradient for b1; same shape as b1"),
    ],
    columns=["object", "shape", "meaning"],
)
display(gradient_shapes)

display(pd.DataFrame(np.round(grads_demo["dW1"], 4), columns=["x1", "x2"], index=[f"hidden_{i}" for i in range(1, 5)]))
display(pd.DataFrame(np.round(grads_demo["dW2"], 4), columns=[f"hidden_{i}" for i in range(1, 5)], index=["output"]))

	object	shape	meaning
0	dZ2	(1, 5)	prediction error at the output layer
1	dW2	(1, 4)	gradient for W2; same shape as W2
2	db2	(1, 1)	gradient for b2; same shape as b2
3	dZ1	(4, 5)	hidden-layer error after backprop
4	dW1	(4, 2)	gradient for W1; same shape as W1
5	db1	(4, 1)	gradient for b1; same shape as b1

	x1	x2
hidden_1	0.1082	-0.1240
hidden_2	0.1138	-0.1241
hidden_3	-0.0923	0.1004
hidden_4	-0.0371	0.0381

	hidden_1	hidden_2	hidden_3	hidden_4
output	0.078	-0.1575	-0.1444	0.2565

Training the Network

Now we put the forward pass, loss, backpropagation, and parameter updates into one training loop.

The most important detail is that every step remains vectorized. We never loop over individual examples during learning. Each iteration processes the entire batch with matrix multiplications.

def predict_probabilities(X, parameters):
    probs, _ = forward_propagation(X, parameters)
    return probs


def predict_classes(X, parameters, threshold=0.5):
    return (predict_probabilities(X, parameters) >= threshold).astype(int)


def train_network(X, Y, n_h=4, learning_rate=1.0, iterations=4000, snapshot_every=200, seed=7):
    parameters = initialize_parameters(n_x=X.shape[0], n_h=n_h, n_y=1, seed=seed)
    history = []
    snapshots = []

    for iteration in range(iterations + 1):
        A2, cache = forward_propagation(X, parameters)
        loss = compute_loss(A2, Y)
        grads = backward_propagation(parameters, cache, X, Y)
        parameters = update_parameters(parameters, grads, learning_rate)

        if iteration % snapshot_every == 0:
            train_preds = (A2 >= 0.5).astype(int)
            train_accuracy = np.mean(train_preds == Y)
            history.append(
                {
                    "iteration": iteration,
                    "loss": float(loss),
                    "train_accuracy": float(train_accuracy),
                }
            )
            snapshots.append(
                {
                    "iteration": iteration,
                    "parameters": {key: value.copy() for key, value in parameters.items()},
                }
            )

    return parameters, pd.DataFrame(history), snapshots


trained_parameters, history_df, snapshots = train_network(
    X_train,
    Y_train,
    n_h=4,
    learning_rate=1.0,
    iterations=4000,
    snapshot_every=200,
    seed=7,
)

train_accuracy = np.mean(predict_classes(X_train, trained_parameters) == Y_train)
test_accuracy = np.mean(predict_classes(X_test, trained_parameters) == Y_test)

summary_df = pd.DataFrame(
    [
        ("train_accuracy", round(float(train_accuracy), 4)),
        ("test_accuracy", round(float(test_accuracy), 4)),
        ("final_loss", round(float(history_df["loss"].iloc[-1]), 4)),
    ],
    columns=["metric", "value"],
)
display(summary_df)

metrics_long = history_df.melt(id_vars="iteration", var_name="metric", value_name="value")
loss_fig = px.line(
    metrics_long,
    x="iteration",
    y="value",
    color="metric",
    markers=True,
    template="plotly_white",
    title="Training Curves",
)
loss_fig.update_layout(height=500)
loss_fig.show()

	metric	value
0	train_accuracy	0.9933
1	test_accuracy	0.9200
2	final_loss	0.0432

Interactive Decision Boundary

This chart shows the classifier surface as training progresses.

• Blue and red shading represent the model probability across the feature space.
• The slider lets you inspect how the hidden layer gradually bends the boundary.
• Because the data was standardized, both axes are in standardized feature units rather than the original raw coordinates.

def boundary_surface(parameters, x_min, x_max, y_min, y_max, steps=120):
    xx, yy = np.meshgrid(
        np.linspace(x_min, x_max, steps),
        np.linspace(y_min, y_max, steps),
    )
    grid = np.c_[xx.ravel(), yy.ravel()].T
    probs, _ = forward_propagation(grid, parameters)
    zz = probs.reshape(xx.shape)
    return xx, yy, zz


x_min, x_max = X_train_rows[:, 0].min() - 1.0, X_train_rows[:, 0].max() + 1.0
y_min, y_max = X_train_rows[:, 1].min() - 1.0, X_train_rows[:, 1].max() + 1.0

xx, yy, zz0 = boundary_surface(snapshots[0]["parameters"], x_min, x_max, y_min, y_max)
x_axis = xx[0]
y_axis = yy[:, 0]

contour = go.Contour(
    x=x_axis,
    y=y_axis,
    z=zz0,
    colorscale="RdBu",
    opacity=0.5,
    showscale=True,
    zmin=0,
    zmax=1,
    colorbar=dict(title="P(class = 1)"),
)

scatter = go.Scatter(
    x=X_train_rows[:, 0],
    y=X_train_rows[:, 1],
    mode="markers",
    marker=dict(
        color=y_train,
        colorscale="Viridis",
        size=8,
        line=dict(color="white", width=1),
    ),
    text=[f"class={label}" for label in y_train],
    hovertemplate="x1=%{x:.2f}<br>x2=%{y:.2f}<br>%{text}<extra></extra>",
    showlegend=False,
)

frames = []
for snapshot in snapshots:
    _, _, zz = boundary_surface(snapshot["parameters"], x_min, x_max, y_min, y_max)
    frames.append(
        go.Frame(
            name=str(snapshot["iteration"]),
            data=[
                go.Contour(
                    x=x_axis,
                    y=y_axis,
                    z=zz,
                    colorscale="RdBu",
                    opacity=0.5,
                    showscale=True,
                    zmin=0,
                    zmax=1,
                )
            ],
        )
    )

slider_steps = [
    {
        "args": [[str(snapshot["iteration"])], {"frame": {"duration": 0, "redraw": True}, "mode": "immediate"}],
        "label": str(snapshot["iteration"]),
        "method": "animate",
    }
    for snapshot in snapshots
]

boundary_fig = go.Figure(data=[contour, scatter], frames=frames)
boundary_fig.update_layout(
    title="How the Decision Boundary Changes During Training",
    template="plotly_white",
    height=650,
    xaxis_title="x1 (standardized)",
    yaxis_title="x2 (standardized)",
    updatemenus=[
        {
            "type": "buttons",
            "buttons": [
                {
                    "label": "Play",
                    "method": "animate",
                    "args": [None, {"frame": {"duration": 300, "redraw": True}, "fromcurrent": True}],
                },
                {
                    "label": "Pause",
                    "method": "animate",
                    "args": [[None], {"frame": {"duration": 0, "redraw": False}, "mode": "immediate"}],
                },
            ],
            "x": 0.02,
            "y": 1.08,
        }
    ],
    sliders=[
        {
            "currentvalue": {"prefix": "iteration="},
            "pad": {"t": 50},
            "steps": slider_steps,
        }
    ],
)
boundary_fig.show()

Final Takeaways

A few patterns are worth remembering:

• X is transposed so that one matrix multiply can process every example in the batch.
• every row in a weight matrix belongs to one neuron in the next layer.
• every column in an activation matrix belongs to one example.
• every gradient must have the exact same shape as the parameter it updates.

If those four ideas stay clear, most neural-network notation becomes much easier to read.

From here, the next natural step is to add more hidden layers, swap tanh for ReLU, or introduce mini-batch training. The shape logic stays the same. Only the architecture gets deeper.