How Neural Networks Work: Math, Intuition, and Code

Introduction

In the last post of our Deep Learning Series ↗, we built a Perceptron, a single neuron that draws a straight line to separate two classes. It has one hard limitation: the moment data is not linearly separable, it completely fails.

Real-world data almost never is. So what if you stacked many neurons together, where each layer builds on the one before it?

That is the neural network, the foundation of nearly everything exciting in modern machine learning. Let's build one from scratch.

What is a Neural Network?

A neural network is a parametric function that maps inputs to outputs through a chain of layers. Each layer performs the same two-step computation: a linear transformation followed by a non-linear activation. Stacking many such layers lets the network learn hierarchical representations — early layers detect simple patterns, later layers combine those into complex concepts.

Every fully-connected neural network has three kinds of layers:

• Input layer — receives raw features. No computation happens here. For MNIST ↗ it is 784 neurons, one per pixel of a 28×28 image.
• Hidden layers — the core of the network. Each neuron takes a weighted sum of every neuron in the previous layer, adds a bias, and passes the result through a non-linear activation. You can stack as many as you want.
• Output layer — produces the final prediction. For multi-class classification we use a softmax activation to convert raw scores into a probability distribution over classes.

The architecture we use in this post:

Input (784)  →  Hidden (128, ReLU)  →  Hidden (64, ReLU)  →  Output (10, Softmax)

In our implementation this is expressed as a single list passed at initialization:

nn = NeuralNetwork(
    layer_sizes   = [784, 128, 64, 10],
    learning_rate = 0.1,
    iterations    = 500,
)

Why Do Layers Work?

Consider image recognition. The first hidden layer sees raw pixel intensities. Its neurons learn to detect local patterns like edges at various orientations, color transitions, gradients. The second layer sees those edge detectors and learns to combine them into shapes: corners, curves, circles. A third layer combines shapes into parts. The output layer says "this combination of parts corresponds to class 7."

This is the power of hierarchical composition. A single linear layer can only draw hyperplanes. Stack non-linear transformations and you can carve the input space into arbitrarily complex regions.

The Universal Approximation Theorem formalizes this: a neural network with at least one hidden layer and a non-linear activation can approximate any continuous function on a compact domain to any desired precision, given enough neurons. Depth is not strictly required by the theorem, but in practice deeper networks need exponentially fewer neurons to represent the same functions.

Mathematical Intuition

Forward Pass

The forward pass is the network making a prediction. Data goes in at the left, passes through every layer one by one, and probabilities come out at the right.

Each layer does exactly two things:

Step 1: Linear combination. Every neuron takes a weighted sum of all the neurons from the previous layer, plus a bias term. For each layer $i$ from $1$ to $L$:

Think of it like this: $W$ holds the weights (how much each input matters), $A^{[i-1]}$ is the output of the previous layer, and $b$ is a bias that shifts the result. The dot product $W \cdot A$ does all the multiplications and additions in one go, and $b$ (shape (n_out, 1)) is added to every column via broadcasting — the same bias vector is applied to every example in the batch.

Step 2: Activation. Pass the result through a non-linear function:

Hidden layers use ReLU, and the output layer uses Softmax.

Why do we need the activation function at all?

If we skip it and just use $f(z) = z$, the two steps above collapse into a single linear map. Ignoring biases for a moment to see the core issue:

where $X$ is the input matrix (shape $n_0 \times m$, one column per example) and $A^{[L]}$ is the final output.

All those weight matrices multiplied together collapse into one big matrix — call it $W_{\text{eff}} = W^{[L]} \cdots W^{[1]}$. No matter how many layers you add, the whole network behaves like a single linear layer: $A^{[L]} = W_{\text{eff}} \cdot X$ (plus a combined bias term that also collapses in the same way). It can only learn straight-line decision boundaries, just like the Perceptron. The activation function is what breaks this collapse and gives the network the ability to learn curves and complex shapes.

ReLU (Rectified Linear Unit)

ReLU is the simplest possible non-linearity: if the input is positive, keep it; if it is negative, set it to zero.

Why is this useful? A neuron with a large positive $z$ is "firing" — it detected something in the input. A neuron with a negative $z$ is "off." ReLU makes this binary distinction crisp.

Why ReLU specifically? Compared to sigmoid or tanh:

• It does not saturate for positive inputs, so gradients do not vanish in that region.
• It is computationally trivial — just a threshold.
• Empirically it trains faster and works well in deep networks.

The one risk is the dying ReLU problem: if a neuron's pre-activation $z$ is always negative, its output is always zero, its gradient is always zero, and it never updates. Small initial weights and a modest learning rate prevent this in practice.

Softmax (Output Layer)

After the last hidden layer, the network produces a vector of raw real-valued scores, one per class — the logits:

where $K$ is the number of classes (10 for MNIST digits 0–9) and $z_k$ is the raw score assigned to class $k$. A logit can be any real number — positive, negative, large, small — and has no direct probabilistic meaning on its own.

For the output to be interpretable as a probability distribution, we need values that are all in $[0, 1]$ and sum to exactly $1$ across all classes. Softmax does exactly this, mapping $z$ to a new vector $\hat y$ where the $k$-th entry is:

Why exponentiate? Two reasons. First, exponentiating maps any real number to a positive number, which is required before you can interpret values as probabilities. Second, it makes the function differentiable everywhere, which backpropagation requires. A hard max (1 for the top class, 0 for everything else) is not differentiable.

Applied column-wise over a batch. In our implementation, $Z^{[L]}$ is a matrix of shape $(K, m)$ — each of the $m$ columns is the logit vector for one training example. Softmax is applied independently to each column, turning the entire matrix into a probability matrix in one vectorized operation:

Z[L] shape:  (K, m)   — raw scores, K classes, m examples
A[L] shape:  (K, m)   — probabilities, each column sums to 1

Numerical stability. There is a hidden danger in the formula as written. If any logit $z_k$ is large (say, 1000), $e^{1000}$ overflows to inf in floating point. The fix: subtract the maximum value in the column before exponentiating. This does not change the output mathematically, because the shift cancels in numerator and denominator:

In code:

def softmax(Z):
    shifted = Z - np.max(Z, axis=0, keepdims=True)  # subtract max per example
    exp_Z = np.exp(shifted)
    return exp_Z / np.sum(exp_Z, axis=0, keepdims=True)

The Loss Function: Cross-Entropy

After the forward pass, the network outputs a probability for each class. We need a single number that says how wrong those probabilities are — that is the loss.

We use cross-entropy loss. The idea is simple: look only at the probability the network assigned to the correct class, and penalize it if that probability is low.

For a single example with true class $c$:

where $A^{[L]}_c$ is the predicted probability assigned to the correct class. The negative log is a clever choice: it is $0$ when the network is perfectly confident and correct, and it shoots toward infinity as the predicted probability for the right class drops toward zero.

Over the full batch of $m$ examples, the average cross-entropy is:

where $c_j$ is the true class of the $j$-th example.

Why cross-entropy and not mean squared error? MSE works poorly with softmax because the gradient becomes very small when the network is confidently wrong (due to the exponential saturating). Cross-entropy avoids this — as we'll see, its gradient with softmax simplifies to something clean and strong.

Backpropagation

The forward pass gives us a prediction. The loss tells us how wrong we are. Now we need: which weights and biases caused the error, and by how much?

That is what backpropagation does. It applies the chain rule layer by layer, moving from the output toward the input. For every parameter $\theta$ we want $\partial \mathcal{L} / \partial \theta$ — how much a small change in $\theta$ would change the loss. Gradient descent then moves $\theta$ in the opposite direction.

Step 1 — Output Layer Gradient

Let $\hat y_k = A^{[L]}_k$ be the predicted probability for class $k$, and let the true class be $c$ (one-hot: $y_k = 1$ if $k=c$, else $0$). The loss is $\mathcal{L} = -\log \hat y_c$. We need $\partial \mathcal{L} / \partial z_k$ for each $k$.

Case 1: $k = c$ (the correct class).

The softmax self-derivative is $\frac{\partial \hat y_k}{\partial z_k} = \hat y_k(1-\hat y_k)$, so:

Case 2: $k \neq c$ (wrong class).

The softmax cross-derivative is $\frac{\partial \hat y_j}{\partial z_k} = -\hat y_j \hat y_k$ (for $k \neq j$), so:

Putting both cases together, $\frac{\partial \mathcal{L}}{\partial z_k} = \hat y_k - y_k$ for every $k$. In matrix form over the full batch:

This is the famous clean result: the gradient of softmax + cross-entropy is simply the difference between predicted probabilities and true one-hot labels. The two messy derivatives cancel each other out.

Step 2 — Gradients for the Output Layer Parameters

Since $Z^{[L]} = W^{[L]} \cdot A^{[L-1]} + b^{[L]}$, the chain rule gives:

The $\frac{1}{m}$ averages the gradient over the batch so the learning rate doesn't depend on batch size.

Step 3 — Propagating the Error Backward Through Hidden Layers

For each hidden layer $i$, two things happen:

Route the error back through the weights. Transposing $W$ reverses the forward direction of information flow:

Gate through ReLU. ReLU zeroed out any neuron with a negative pre-activation; those neurons contributed nothing in the forward pass, so they receive no gradient in the backward pass either:

Then $dW^{[i]}$ and $db^{[i]}$ follow the exact same formula as the output layer. We repeat this for every hidden layer, moving right to left.

Gradient Descent: The Update Rule

Once we have all gradients, the parameter update is:

Applied to every layer:

$\alpha$ is the learning rate — the step size. The gradient points toward steepest increase in loss, so subtracting it moves parameters toward a minimum.

Why Initialize Weights Randomly

Before training starts, all weights must be set to some initial values. Why not all zeros?

If $W^{[i]} = 0$ for all $i$, then every neuron in a layer computes exactly the same output, receives exactly the same gradient, and updates identically. No amount of training can break this symmetry — all neurons stay identical forever.

The fix is random initialization: sample weights from a distribution with small variance, e.g. $W^{[i]} \sim \mathcal{N}(0, 0.01)$. In practice this means sampling from the standard normal distribution (np.random.randn, mean 0, variance 1) and scaling by $0.1$ — since scaling a variance-1 distribution by $0.1$ gives it variance $0.1^2 = 0.01$, exactly matching the target. This breaks symmetry while keeping initial values small enough to avoid saturating ReLU on the very first forward pass.

Implementation

Now let's implement the neural network from scratch in Python. This version is a corrected, complete pass over the earlier draft — every method referenced below is actually defined, and the naming is consistent throughout (all math helpers live on a Utils class and are called via self.utils.<fn>, not as bare globals). Full code: GitHub ↗.

Step 1 — Utility Functions

All static, pure math with no object state:

class Utils:
 
    def ReLU(Z):
        return np.maximum(0, Z)
 
    def derivative_ReLU(Z):
        # 1 where Z > 0, 0 elsewhere — gates the backward gradient
        return Z > 0
 
    def softmax(Z):
        shifted = Z - np.max(Z, axis=0, keepdims=True)
        exp_Z   = np.exp(shifted)
        return exp_Z / np.sum(exp_Z, axis=0, keepdims=True)
 
    def one_hot(Y):
        # label 3 with 10 classes -> [0,0,0,1,0,0,0,0,0,0]
        encoded_Y = np.zeros((Y.size, Y.max() + 1))
        encoded_Y[np.arange(Y.size), Y] = 1
        return encoded_Y.T  # -> (num_classes, m)
 
    def accuracy(y_pred, y):
        return np.sum(y_pred == y) / y.size

Step 2 — Initialize the Network Class

class NeuralNetwork:
 
    def __init__(self, layer_sizes, learning_rate=0.1, iterations=500):
        if len(layer_sizes) < 2:
            raise ValueError("layer_sizes must have at least 2 entries (input + output).")
 
        self.layer_sizes = layer_sizes
        self.alpha       = learning_rate
        self.n_iters     = iterations
        self.n_layers    = len(layer_sizes) - 1   # number of weight matrices
        self.utils       = Utils
 
        self.weights = []   # weights[i] shape: (layer_sizes[i+1], layer_sizes[i])
        self.biases  = []   # biases[i]  shape: (layer_sizes[i+1], 1)
 
        self.Z = []   # pre-activation values per layer
        self.A = []   # post-activation values per layer (A[0] == X)

layer_sizes is the one knob that defines the whole architecture: [784, 10] gives a network with zero hidden layers (a plain linear classifier), [784, 128, 64, 10] gives two hidden layers, and so on.

Step 3 — Initialize Parameters

def init_params(self):
    self.weights = []
    self.biases  = []
 
    for i in range(self.n_layers):
        n_in, n_out = self.layer_sizes[i], self.layer_sizes[i + 1]
        W = np.random.randn(n_out, n_in) * 0.1   # small random weights
        b = np.zeros((n_out, 1))                 # zero biases
        self.weights.append(W)
        self.biases.append(b)

Rows of $W$ correspond to neurons in the next layer; columns correspond to neurons in the current layer feeding into it.

Step 4 — Forward Propagation

def forward(self, X):
    self.Z = [None] * self.n_layers
    self.A = [None] * (self.n_layers + 1)
    self.A[0] = X
 
    for i in range(self.n_layers):
        self.Z[i] = self.weights[i].dot(self.A[i]) + self.biases[i]
 
        if i == self.n_layers - 1:
            self.A[i + 1] = self.utils.softmax(self.Z[i])   # output layer
        else:
            self.A[i + 1] = self.utils.ReLU(self.Z[i])      # hidden layers

Step 5 — Backward Propagation

def backward(self, X, Y):
    m  = Y.size
    dZ = [None] * self.n_layers
    self.dW = [None] * self.n_layers
    self.db = [None] * self.n_layers
 
    one_hot_Y = self.utils.one_hot(Y)
    dZ[-1] = self.A[-1] - one_hot_Y   # softmax + cross-entropy gradient
 
    for i in reversed(range(self.n_layers)):
        self.dW[i] = 1 / m * dZ[i].dot(self.A[i].T)
        self.db[i] = 1 / m * np.sum(dZ[i], axis=1, keepdims=True)
 
        if i > 0:
            dZ[i - 1] = (
                self.weights[i].T.dot(dZ[i])
                * self.utils.derivative_ReLU(self.Z[i - 1])
            )

Step 6 — Update Parameters

def update_params(self):
    for i in range(self.n_layers):
        self.weights[i] -= self.alpha * self.dW[i]
        self.biases[i]  -= self.alpha * self.db[i]

Step 7 — Inference Helpers

These were missing from the earlier draft, but they're needed for train() to report accuracy and for the final evaluation to work at all:

def predict(self, A_out):
    # argmax along axis=0 picks the highest-probability class per example
    return np.argmax(A_out, axis=0)
 
def make_predictions(self, X):
    self.forward(X)
    return self.predict(self.A[-1])
 
def test_prediction(self, index, X, Y):
    image = X[:, index, None]   # (784, 1)
    pred  = self.make_predictions(image)
    label = Y[index]
    print(f"Prediction : {pred[0]}")
    print(f"True label : {label}")
    plt.gray()
    plt.imshow(image.reshape(28, 28), interpolation="nearest")
    plt.title(f"Predicted: {pred[0]}  |  Label: {label}")
    plt.axis("off")
    plt.show()

Step 8 — Training Loop

def train(self, X, Y):
    self.init_params()
 
    for i in range(self.n_iters):
        self.forward(X)
        self.backward(X, Y)
        self.update_params()
 
        if i % 50 == 0:
            preds = self.predict(self.A[-1])
            acc   = self.utils.accuracy(preds, Y)
            print(f"Iteration {i:>4}  |  Training accuracy: {acc:.4f}")

Note this is full-batch gradient descent — every iteration computes the forward and backward pass over all 60,000 training examples at once, rather than mini-batches. That keeps the implementation simple and is fine at MNIST's scale on CPU, but it's worth naming: it's a deliberate simplification, not the industry-standard approach (which is mini-batch SGD or Adam, typically with batch sizes of 32–256).

Step 9 — Architecture Summary (optional, for inspection)

def summary(self):
    print("=" * 44)
    print(f"  Neural Network  ,  {self.n_layers} layer(s)")
    print("=" * 44)
    print(f"  {'Layer':<12} {'Shape (out x in)':<20} Params")
    print("-" * 44)
 
    labels = (
        ["Input -> H1"]
        + [f"H{i} -> H{i+1}" for i in range(1, self.n_layers - 1)]
        + ["Hidden -> Out"]
    ) if self.n_layers > 1 else ["Input -> Out"]
 
    total = 0
    for label, W, b in zip(labels, self.weights, self.biases):
        params = W.size + b.size
        total += params
        print(f"  {label:<12} {str(W.shape):<20} {params:,}")
 
    print("-" * 44)
    print(f"  {'Total':<12} {'':<20} {total:,}")
    print("=" * 44)

Step 10 — Run It

if __name__ == "__main__":
    from sklearn.datasets import fetch_openml
 
    mnist = fetch_openml('mnist_784', version=1, as_frame=False)
    X = mnist.data / 255.0          # shape (70000, 784)
    y = mnist.target.astype(int)    # shape (70000,)
 
    indices = np.random.permutation(len(X))
    X, y = X[indices], y[indices]
 
    X_test  = X[:1000].T            # (784, 1000)
    Y_test  = y[:1000]              # (1000,)
    X_train = X[1000:].T            # (784, 60000)
    Y_train = y[1000:]              # (60000,)
 
    nn = NeuralNetwork(
        layer_sizes   = [784, 128, 64, 10],
        learning_rate = 0.1,
        iterations    = 500,
    )
 
    nn.summary()
    nn.train(X_train, Y_train)
    nn.summary()
 
    test_preds    = nn.make_predictions(X_test)
    test_accuracy = Utils.accuracy(test_preds, Y_test)
    print(f"\nTest accuracy: {test_accuracy:.4f}")
 
    nn.test_prediction(index=2, X=X_test, Y=Y_test)

Every method this block calls — summary, train, make_predictions, test_prediction, Utils.accuracy — is now actually defined above, so this runs end to end without the AttributeErrors the earlier draft would have hit.

Output

Common Misconceptions

• "More layers always means better performance." Depth only helps when you have enough data and regularization to support the extra capacity. A 10-layer network trained on 500 examples will almost certainly memorize the training set and fail on anything new.
• "Gradient descent finds the global minimum." It finds a local minimum, often a saddle point. For large overparameterized networks most local minima happen to have similar loss values, so this matters less than the theory suggests — but the global optimum is not guaranteed.
• "Neural networks are complete black boxes." Partially true, but overstated. Activation maximization, saliency maps, SHAP values, and probing classifiers all give genuine insight into what networks learn. Mechanistic interpretability is an active and growing field.
• "You need a GPU to train neural networks." For the MNIST model in this post, a CPU is perfectly fine. Neural networks scale, but they start small.

When to Apply Neural Networks

Use neural networks when:

• Input is high-dimensional and unstructured — images, audio, raw text, sensor signals.
• Data is plentiful — neural networks are data-hungry; thousands to millions of labelled examples is the comfortable range.
• The underlying function is highly complex — too complex for manual feature engineering.
• You can afford hyperparameter tuning — learning rate, architecture, regularization, and batch size all interact.

Prefer simpler models when:

• Data is small or tabular. Logistic regression, random forests, and XGBoost routinely beat neural networks here.
• Interpretability is non-negotiable (medicine, law, finance).
• Inference speed is critical and the deployment environment is constrained.
• You want statistical guarantees neural networks can't provide.

Advantages of Neural Networks

• Universal approximation — given enough capacity, a neural network can approximate any continuous function to any precision.
• Automatic feature learning — no manual feature engineering.
• Scales with data — performance tends to keep improving with more training examples.
• Transfer learning — a network trained on one task can be fine-tuned for a related task with far less data.
• State of the art across domains — vision, NLP, protein folding, game playing.
• Hardware parallelism — the core operation is matrix multiplication, which GPUs execute extremely efficiently.

Disadvantages & Limitations

• Data hunger — strong generalization typically requires large labelled datasets.
• Computational cost — training large networks requires significant hardware and time.
• Hyperparameter sensitivity — learning rate, architecture, initialization, batch size, regularization all interact in non-obvious ways.
• Overfitting — high-capacity networks can memorize training data; requires dropout, weight decay, data augmentation, early stopping.
• Limited interpretability — understanding why a specific prediction was made is hard, and a serious problem in safety-critical applications.
• Gradient pathologies — vanishing and exploding gradients are genuine problems in deep networks. Batch normalization, residual connections, and careful initialization help, but add complexity.
• Not data-efficient — a child generalizes from two examples; neural networks need thousands. This gap remains an open research problem.

Conclusion

In this post we built a neural network from scratch, derived the math end to end, implemented backpropagation, and made the architecture fully configurable in pure NumPy.

But these networks are memoryless. They process each input independently, with no sense of sequence or time. Language and speech don't work that way — the meaning of a word depends on everything before it.

That's where Recurrent Neural Networks come in. See you in the next one.