A beginner-friendly walk through how a neural network actually works under the hood, built in plain NumPy with no deep learning library. Every equation is followed by the code that implements it.
Full code: https://github.com/siva-rgb/Neural_Network_Scratch
Why I wrote this
When most of us start with machine learning, the path looks the same. You learn how an algorithm works in theory, you implement it with a library like scikit-learn, TensorFlow, or PyTorch, you pick up a few optimization and regularization tricks, and you move on to the next algorithm. I learned it that way too.
But neural networks always nagged at me. The library made them feel like a black box. Input goes into a hidden layer, weights and biases do something, that repeats until the final layer, and then some backpropagation step updates all the weights using gradient descent. It worked, but I could not have told you what was really happening between model.fit() and the accuracy going up.
So I sat down and rebuilt a neural network from scratch, using nothing but NumPy for the array math. A lecture series by Dr. Raj Dandekar helped a lot along the way. This article is the result: the theory and the math behind a neural network, translated line by line into Python. By the end, the black box should feel like a glass box.
First, the easy way
Here is a small neural network in TensorFlow, so we have something to compare against.
import numpy as np
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
# 1. Generate dummy data (1000 samples, 20 features)
X_train = np.random.random((1000, 20))
y_train = np.random.randint(2, size=(1000, 1)) # binary labels (0 or 1)
# 2. Build the network
model = Sequential([
Dense(64, activation='relu', input_shape=(20,)), # first hidden layer
Dense(32, activation='relu'), # second hidden layer
Dense(1, activation='sigmoid') # output layer
])
# 3. Compile
model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])
# 4. Train
model.fit(X_train, y_train, epochs=5, batch_size=32)
This is easy to read. Sequential stacks the layers, each Dense sets the number of neurons and the activation, compile picks the optimizer and the loss, and fit trains for 5 epochs in batches of 32. Four steps and you have a working network.
But that is exactly the problem. Every hard question, how the forward pass is computed, what the loss really measures, how backpropagation finds better weights, is hidden inside those four calls. So let us throw the library away and build all of it ourselves.
What a neural network actually is
Once I actually dug into it, a neural network turned out to be just two operations, repeated over and over:
- A linear step: multiply the inputs by some weights, add a bias.
- A non-linear step: pass the result through an activation function.
One round of those two is a layer. Stack a handful of them and you have got a deep network. And training? That is just nudging all those weights and biases, again and again, until the output starts matching the answer you actually wanted.
That is really the whole idea. Everything from here is just making those two steps precise and then writing them in code.
We will use a small, concrete example throughout: a fully connected network with 4 input features and a hidden layer of 4 neurons.
The forward pass: the math
Our single input example is a vector of four features:
Every one of the 4 neurons has its own set of 4 weights, one per input, so the weights form a matrix:
Each neuron also has one bias:
The output of the layer, before any activation, is:
Read that in plain words: for each neuron, multiply every input by its matching weight, add them up, then add the neuron's bias. That single number is the neuron's pre-activation output. Do it for all four neurons and you get a vector of four values.
The forward pass: the code
Here is the exact same computation in NumPy.
import numpy as np
x_input = [2, 3, 4, 5]
w = [[0.3, 0.34, 0.4, 0.03],
[0.5, 0.2, 0.1, 0.22],
[0.23, 0.37, 0.4, 0.3 ],
[0.02, 0.14, 0.17, 0.05]]
b = [2, -0.3, 5, 7]
dot_product = np.dot(np.array(w), np.array(x_input)) + np.array(b)
print(dot_product)
# [5.37 2.8 9.67 8.39]
One thing that tripped me up at first: we did not actually write a transpose here, even though the math clearly said
. That is because NumPy treats a 1D array like x_input as a column vector of shape (4,). So np.dot(W, x) already computes, for each row of W (each neuron), the dot product with the input. The result is the same four numbers. When we move to batches of inputs later, the shapes change and the convention flips to X · W, which is why the reusable class below multiplies in that order. Shapes are the thing to keep your eye on in this whole exercise.
To prove there is no magic in np.dot, here is the same thing with plain Python loops:
layer_output = []
for weight_row, bias in zip(w, b):
output = 0
for n_input, weight in zip(x_input, weight_row):
output += n_input * weight # weighted sum
output += bias # add the bias
layer_output.append(output)
print(layer_output)
# [5.37, 2.8000000000000003, 9.67, 8.39]
Identical result. The loop makes the math obvious, but it is slow: its time complexity is roughly O(inputs times neurons), and it gets worse fast as layers grow. NumPy does the same arithmetic in optimized C under the hood, so we use it for everything from here on.
Finally, let us wrap the forward pass in a reusable class, the same idea as a Keras Dense layer:
class Layer_Dense:
def __init__(self, n_input, n_neurons):
# small random weights break symmetry so neurons learn different things
self.weights = 0.01 * np.random.randn(n_input, n_neurons)
self.bias = np.zeros((1, n_neurons))
def forward(self, inputs):
self.output = np.dot(inputs, self.weights) + self.bias
Two small but important choices. Weights start as small random numbers so that no two neurons are identical at the start, if they all started equal they would stay equal and the network could never learn. Biases start at zero, which is fine.
With this, the forward pass of one layer is done, and every hidden layer works the same way.
Activation functions: giving the network its power
If we only ever stacked linear steps ( ), the whole network would collapse into a single linear function, no matter how many layers. It could never learn a curve. Activation functions add the non-linearity that lets networks model complicated patterns.
ReLU, for hidden layers
The most common activation for hidden layers is ReLU, the Rectified Linear Unit. It is almost embarrassingly simple:
Anything negative becomes zero, anything positive passes through unchanged.
class Activation_ReLU:
def forward(self, inputs):
self.output = np.maximum(0, inputs)
Softmax, for the output layer
For the final layer of a multi-class classifier we want probabilities that sum to 1. That is what softmax gives us. It exponentiates each value and then normalizes:
There is one practical trick. Before exponentiating, we subtract the largest value in each row. This does not change the result mathematically (the ratio is identical), but it keeps the exponentials from blowing up to huge numbers and overflowing. This is called the numerical-stability trick.
class Activation_Softmax:
def forward(self, inputs):
# subtract the row max for numerical stability
exp_values = np.exp(inputs - np.max(inputs, axis=1, keepdims=True))
probabilities = exp_values / np.sum(exp_values, axis=1, keepdims=True)
self.output = probabilities
For binary classification you would use the sigmoid function instead, which squashes a single value into the range 0 to 1. Here we will stick with softmax for the multi-class case.
The loss: measuring how wrong we are
Once the network produces probabilities, we need one number that says how far off it was. For classification we use categorical cross-entropy. For a single sample it is:
where
is the probability the network assigned to the true class. If the network was confident and right (say 0.99), the loss is tiny. If it was confident and wrong (say 0.01 for the true class), the loss is large. That is exactly the behaviour we want.
class Loss:
def calculate(self, output, y):
sample_losses = self.forward(output, y)
return np.mean(sample_losses) # average loss over the batch
class Loss_CategoricalCrossentropy(Loss):
def forward(self, y_pred, y_true):
samples = len(y_pred)
# clip to avoid log(0), which is infinity
y_pred_clipped = np.clip(y_pred, 1e-7, 1 - 1e-7)
if len(y_true.shape) == 1: # sparse labels
correct_confidences = y_pred_clipped[range(samples), y_true]
elif len(y_true.shape) == 2: # one-hot labels
correct_confidences = np.sum(y_pred_clipped * y_true, axis=1)
return -np.log(correct_confidences)
What the shape == 1 vs shape == 2 check is doing. This is the part that confuses people, so let us be explicit. Labels can arrive in two formats:
-
Sparse (1D),
shape == 1: the labels are just the class indices, for example[0, 2, 1]meaning sample 0 is class 0, sample 1 is class 2, sample 2 is class 1. To grab the predicted probability of the correct class we index directly:y_pred_clipped[range(samples), y_true]picks one probability per row. -
One-hot (2D),
shape == 2: the labels are vectors like[[1,0,0], [0,0,1], [0,1,0]]. Here we multiply element-wise by the predictions and sum across the row, which keeps only the probability at the "1" position.
Both paths end up with the same thing, the probability the network gave to the right answer. We also clip the predictions to a tiny range so we never take log(0), which would be negative infinity.
Backpropagation: how the network learns
This is the step that felt most like magic to me, so let us go slowly.
Forward propagation gave us a loss. Backpropagation answers a single question: if I nudge each weight and bias a little, how does the loss change? That answer is the gradient. Once we know it, we move each parameter in the direction that reduces the loss. Repeat thousands of times and the network gets better.
The tool behind it is the chain rule from calculus, and honestly, once that clicked for me it stopped being scary. The loss depends on the output, the output depends on the activation, the activation depends on the weights, and so on down the chain. The chain rule just lets us multiply all those little local derivatives together and push the error backward, from the loss all the way to the first layer.
The guiding principle for the code is simple and worth stating clearly:
Every operation that has a forward pass also has a backward pass. If the forward step computed an output from an input, the backward step computes the gradient of the loss with respect to that input.
Let us work through each piece.
The dense layer, backward
For the linear step , calculus gives three gradients:
The gradient with respect to the weights tells us how to update them. The gradient with respect to the input is what we pass to the previous layer, so it can do its own update. The bias gradient is just the sum of the incoming gradients.
class Layer_Dense_v1:
def __init__(self, n_inputs, n_neurons):
self.weights = 0.01 * np.random.randn(n_inputs, n_neurons)
self.bias = np.zeros((1, n_neurons))
def forward(self, inputs):
self.inputs = inputs # remember inputs for the backward pass
self.output = np.dot(inputs, self.weights) + self.bias
def backward(self, dvalues):
self.dweights = np.dot(self.inputs.T, dvalues) # dL/dW
self.dbias = np.sum(dvalues, axis=0, keepdims=True) # dL/db
self.dinputs = np.dot(dvalues, self.weights.T) # dL/dX -> goes to previous layer
Notice we store self.inputs during the forward pass. The backward pass needs them, which is why frameworks keep a record of the forward computation.
ReLU, backward
The derivative of ReLU is beautifully simple: it is 1 wherever the input was positive, and 0 wherever it was negative or zero. So the backward pass just lets the gradient through where the input was positive and blocks it everywhere else.
class Activation_ReLU:
def forward(self, X):
self.input = X
self.output = np.maximum(0, X)
def backward(self, dvalue):
self.dinputs = dvalue.copy()
self.dinputs[self.input <= 0] = 0 # gradient is 0 where input was not positive
Softmax and cross-entropy, backward, together
You could compute the derivative of the loss with respect to the softmax output, and then the derivative of softmax with respect to its input, and multiply them. But that is expensive and messy, because the softmax derivative is a full matrix (the Jacobian). It turns out that when softmax is paired with cross-entropy loss, almost everything cancels, and the combined gradient collapses to something wonderfully clean:
That is it. The predicted probabilities minus the true labels. This is why the two are almost always implemented together.
class Activation_Softmax_Loss_CategoricalCrossentropy:
def __init__(self):
self.activation = Activation_Softmax()
self.loss = Loss_CategoricalCrossentropy()
def forward(self, inputs, y_true):
self.activation.forward(inputs)
self.output = self.activation.output
return self.loss.calculate(self.output, y_true)
def backward(self, dvalues, y_true):
samples = len(dvalues)
# if labels are one-hot, convert to class indices
if len(y_true.shape) == 2:
y_true = np.argmax(y_true, axis=1)
self.dinputs = dvalues.copy()
self.dinputs[range(samples), y_true] -= 1 # subtract 1 from the true class: (y_pred - y_true)
self.dinputs = self.dinputs / samples # average over the batch
Look at what self.dinputs[range(samples), y_true] -= 1 does. Starting from the predicted probabilities, it subtracts 1 from the entry of the correct class. That is exactly
in one line, because the one-hot true label is 1 only at the correct class. Then we divide by the number of samples to get the average gradient.
The optimizer: turning gradients into updates
We now have the gradients. The optimizer decides how to use them. The simplest one is plain gradient descent: move each parameter a small step in the opposite direction of its gradient.
where
is the learning rate, how big a step we take.
class Optimizer_GD:
def __init__(self, learning_rate=1.0):
self.learning_rate = learning_rate
def update_params(self, layer):
layer.weights -= self.learning_rate * layer.dweights
layer.bias -= self.learning_rate * layer.dbias
Adding learning-rate decay
A fixed learning rate is a compromise. Too large and you overshoot the minimum and bounce around. Too small and training crawls. Decay gives you both: start with big steps to cover ground quickly, then shrink the steps over time to settle precisely into the minimum.
where
is the iteration number.
class Optimizer_GD_with_decay:
def __init__(self, learning_rate=1.0, decay=0.0):
self.learning_rate = learning_rate
self.current_learning_rate = learning_rate
self.decay = decay
self.iterations = 0
def pre_update_params(self):
if self.decay:
self.iterations += 1
self.current_learning_rate = self.learning_rate * (1. / (1. + self.decay * self.iterations))
def update_params(self, layer):
layer.weights -= self.current_learning_rate * layer.dweights
layer.bias -= self.current_learning_rate * layer.dbias
The benefit over a plain learning rate is faster early progress with a more stable, precise finish. The shortcoming is that decay is blind to the actual terrain of the loss. It shrinks the step on a schedule whether or not that is the right thing to do at that moment, and it can slow down too early on hard problems.
Adding momentum
Momentum fixes some of that blindness by borrowing an idea from physics. Instead of stepping purely on the current gradient, it keeps a running velocity, a memory of the recent gradients, and moves along that. When gradients keep pointing the same way, momentum builds speed. When they oscillate, momentum cancels the wobble out.
where
is the momentum factor.
class Optimizer_GD_with_momentum:
def __init__(self, learning_rate=1.0, decay=0.0, momentum=0.0):
self.learning_rate = learning_rate
self.current_learning_rate = learning_rate
self.decay = decay
self.iterations = 0
self.momentum = momentum
def pre_update_params(self):
if self.decay:
self.iterations += 1
self.current_learning_rate = self.learning_rate * (1. / (1. + self.decay * self.iterations))
def update_params(self, layer):
if self.momentum:
if not hasattr(layer, 'weight_momentums'):
layer.weight_momentums = np.zeros_like(layer.weights)
layer.bias_momentums = np.zeros_like(layer.bias)
weight_updates = self.momentum * layer.weight_momentums - self.current_learning_rate * layer.dweights
layer.weight_momentums = weight_updates
bias_updates = self.momentum * layer.bias_momentums - self.current_learning_rate * layer.dbias
layer.bias_momentums = bias_updates
else:
weight_updates = -self.current_learning_rate * layer.dweights
bias_updates = -self.current_learning_rate * layer.dbias
layer.weights += weight_updates
layer.bias += bias_updates
The benefit over plain decay is that momentum accelerates through long, consistent slopes and dampens oscillation in narrow valleys, which usually means faster and steadier convergence. It can also carry the network across small bumps and saddle points where plain gradient descent would stall. The cost is one more hyperparameter to tune, and if you set momentum too high it can overshoot and become unstable. This same line of thinking, adaptive per-parameter steps, is what leads to optimizers like Adam.
Putting it all together
Here is the full training loop, using everything we built. Two dense layers, a ReLU in the middle, and the combined softmax-plus-cross-entropy at the end.
dense1 = Layer_Dense_v1(2, 64) # 2 inputs -> 64 neurons
activation1 = Activation_ReLU()
dense2 = Layer_Dense_v1(64, 3) # 64 -> 3 output classes
loss_act = Activation_Softmax_Loss_CategoricalCrossentropy()
optimizer = Optimizer_GD()
for epoch in range(10000):
# ---- forward pass ----
dense1.forward(X)
activation1.forward(dense1.output)
dense2.forward(activation1.output)
loss = loss_act.forward(dense2.output, y)
# ---- track accuracy ----
predictions = np.argmax(loss_act.output, axis=1)
if len(y.shape) == 2:
y = np.argmax(y, axis=1)
accuracy = np.mean(predictions == y)
if not epoch % 100:
print(f'epoch: {epoch}, acc: {accuracy:.3f}, loss: {loss:.3f}')
# ---- backward pass ----
loss_act.backward(loss_act.output, y)
dense2.backward(loss_act.dinputs)
activation1.backward(dense2.dinputs)
dense1.backward(activation1.dinputs)
# ---- update ----
optimizer.update_params(dense1)
optimizer.update_params(dense2)
Read the loop top to bottom and you can see the entire story of a neural network in one screen. Data flows forward through the layers to a loss. The gradient flows backward through the exact same layers in reverse order. The optimizer nudges every weight and bias. Then it all repeats, and with each pass the accuracy climbs and the loss falls.
What I took away from this
Building this from scratch changed how I see every deep learning model. A Dense layer is a weighted sum plus a bias. An activation is a small non-linear function with an equally small derivative. A loss is one honest number. Backpropagation is the chain rule applied carefully, backward. And an optimizer is just a rule for taking a step.
Once those five pieces click, the bigger architectures stop being intimidating. A CNN, an RNN, even a transformer are the same ideas rearranged. The box is no longer black.
If you want to run it yourself, the full implementation is here:
https://github.com/siva-rgb/Neural_Network_Scratch
And if you want the deeper theory and full implementation, Dr. Raj Dandekar's Video Lecture lecture series is well worth your time.
Top comments (0)