Neural Network From Scratch

This post reads best on my site, where inline equations are properly typeset and every reference to a numbered figure or equation is a link: Neural Network from Scratch.

In the previous post in my Machine Learning from Scratch series, I covered logistic regression. Now I will finally jump into the star of modern machine learning: neural networks. As the name suggests, neural networks are loosely inspired by the structure of the human brain. They use artificial models of neurons laced together into a web to learn complex patterns.

In this post, I will build and train a neural network to solve the MNIST handwritten digit classification problem. This is something of a “Hello World” task in machine learning. Specifically, I will create a submission for Kaggle’s Digit Recognizer competition. This has turned into the longest and most involved post yet, so buckle in. Let’s create a neural network from scratch.

Math Primer

Here are some formulas and definitions that will be used throughout the post, shown both in vectorized form and index form when appropriate. There is a lot of ‘index juggling’ in the derivations later on, as well as switching between index and vector notation. It can get confusing, so refer back to these as needed.

Matrix-vector product

Matrix-matrix product

Vector outer product

Hadamard product

Matrix transpose

Kronecker delta

The sifting property of the Kronecker delta allows it to pick out or ‘sift’ a specific element from a summation:

Multivariable chain rule

Let

Then

Artificial Neuron

The building block of neural networks is, unsurprisingly, an artificial neuron model.

A human neuron receives input signals from other neurons through dendrites, processes them in the cell body, and then sends an output signal along the axon.

Figure 1: Neuron diagram. Courtesy of Clker-Free-Vector-Images on Pixabay.

The diagram below shows an (extremely) simplified model of a neuron which can have any number of real-valued inputs. Each input is multiplied by a number called a weight, and the results are then summed together, along with another number called a bias. Without the bias term, if all the inputs to the neuron were zero, their weighted sum would be zero no matter what. This would limit the neuron’s ability to fit data. Lastly, the summed result is passed through an activation function to produce the final output.

Figure 2: Artificial neuron. The bias is shown as an extra input, marked in red, with a weight of 1.

Let us suppose the inputs are x₁, x₂, …, xₙ, the corresponding weights are w₁, w₂, …, wₙ, and the bias is b. The sum of the weighted inputs and the bias is then:

or

The final output is the result of applying the activation function to this sum.

If we collect the inputs and weights into vectors x and w, we can restate the neuron’s computations more compactly as:

Maybe some of this looks familiar if you’ve read my previous posts. I’ll let you in on a secret: everything we’ve covered in the series so far is actually a special case of an artificial neuron!

Linear regression can be modeled by a neuron whose activation function is the identity f(z) = z. Logistic regression can be modeled by a neuron whose activation function is the sigmoid function

If just a single neuron can do all that, imagine what a network of them is capable of…

Fully-Connected Layer

Imagine we have inputs x₁, …, xₙ and feed every input simultaneously into m neurons. This is called a fully-connected layer.

We will denote the weights of the ith neuron by wᵢ₁, wᵢ₂, and so on, meaning the jth weight of the ith neuron is wᵢⱼ. The output of the ith neuron before the activation function is applied (aka the unactivated output) is

Suppose every neuron has the same activation function f. The final output of the ith neuron is then

Figure 3: Fully-connected layer with 3 inputs and 4 outputs.

The inputs to a fully-connected layer could represent many things such as the pixels in an image or data about the passengers of the Titanic. As before, we can collect the inputs into a vector x. Let’s call an input vector a sample. If it is used to train the neural network, it’s a training sample.

We could just feed samples into a layer one-by-one. However, there is a more efficient method called batching. We can arrange sample vectors so that each one is a column of a matrix.

Figure 4: Batching visualized.

We create an input matrix X, a batch of samples where each column is an input vector. In NumPy terms its shape would be (num_inputs, num_samples) where num_inputs is the length of each vector and num_samples is the number of samples in the batch.

In index form, xᵢⱼ is the ith element of the jth sample. The unactivated output of the ith neuron for sample j is

or more compactly,

and the corresponding final output is

In vectorized form, these equations can be written as:

and

In Equation 1, 1 is a vector of 1s of length num_samples. The term b1ᵀ may look a bit confusing but all it does is create a matrix of biases where each column is identical, which is what we want since the same bias vector is used for every sample.

Figure 5: Bias matrix visualized.

Equation 1 and Equation 2 provide a very compact way of describing a fully-connected layer. Such layers will form the main building block of our neural network.

Feedforward Neural Network

We are now ready to create a feedforward neural network. We do this by connecting multiple fully-connected layers in sequence, using the output of one layer as the input to the next. The number of neurons and activation function can vary between layers. The initial point where data is fed into the network is called the input layer. The final layer is called the output layer, and in between are hidden layers.

Figure 6: A feedforward neural network. The arrow shows the direction of data flow from input to output.

Let X be a batch of input data. Using superscripts to index which layer we are referring to, we can write the computation of a feedforward neural network as follows:

where L is the total number of layers.

More compactly,

for i = 1 up to i = L.

In index form,

Elementwise Assumption

The second line of Equation 4 only holds if the activation function is applied elementwise to each neuron’s output. This matches the earlier description of a fully-connected layer, but later we will encounter an activation function called softmax which works differently. Fortunately, it will only be used in the output layer so the formula still applies to the rest of the network.

There is a mathematical result called the Universal Approximation Theorem which shows that a feedforward neural network having at least one hidden layer with a suitable activation function can approximate any continuous function. What this means is that in principle a neural network can model arbitrarily complex relationships between inputs and outputs. The question is how to obtain a network which models a phenomenon of interest. This is where the training process comes in.

Backpropagation

In the post on logistic regression, we discussed gradient descent as a way to minimize a loss function expressing the error between predicted and target values. If you aren’t clear on this, I would recommend taking a look at that post before reading further. We can use gradient descent to train a neural network, but we need to find the gradient of the loss function with respect to all the weights and biases in the network. This is done using a technique called backpropagation, which computes these gradients by applying the chain rule of calculus.

Backpropagation is probably the most important algorithm in modern machine learning. It is what allows neural networks to do everything from recognizing cat pictures to playing Go, predicting how proteins fold, and generating humanlike dialog. As with gradient descent, the basic idea behind it is surprisingly simple, though unfortunately the nitty-gritty math details are not.

A feedforward neural network applies a sequence of operations to input data to produce its output. We can plug this final output into a loss function like Loss(Yᴸ, T) where T is the target output. Assuming we know how to take the derivative of Loss, we can then compute its gradient w.r.t. Yᴸ:

That’s great and all, but we really need the gradients of Loss with respect to the weights and biases in every layer of the neural network. After all, these are the tunable parameters of the network which must be updated for it to learn. How do we obtain these? Using the chain rule.

Consider that the output of a given layer, Yⁱ, influences the loss only through the rest of the sequence of operations in the network. Let ℒⁱ represent the rest of the neural network after layer i, including the loss function. So

By Equation 3,

so

Now suppose we already know the gradient of loss w.r.t. the output of layer i. This is ∇_(Yⁱ)Loss, or in index form,

Via the multivariable chain rule,

By Equation 4, yₘₙⁱ only depends on zⱼₖⁱ when m = j and n = k. In all other cases,

Thus all but one term of the double sum vanishes and we are left with

Now we can find the derivatives of the loss w.r.t. the weights and biases of the layer. Applying the multivariable chain rule to the first line of Equation 4,

Note that if m ≠ j then zₘₙⁱ does not depend on wⱼₖⁱ and thus

Consequently, the expression reduces to

Similarly for the bias,

Finally we can get the derivatives of loss w.r.t. the input to layer i, which is the output of layer i − 1, by once again making use of the chain rule and Equation 4:

Note that yⱼₖ^(i − 1) only occurs in the computation of zₘₙⁱ if k = n (that is, if they belong to the same training sample). Thus the equation simplifies to

We have come full circle. If we know the gradient of the loss w.r.t. the output of layer i we can not only obtain the gradients for that layer’s weights and biases, but also the gradient w.r.t. the output of layer i − 1. This enables us to repeat the whole process for the previous layer, and the one before that, and so on, propagating the gradients back throughout the entire network.

In order to perform the calculations, we need the derivatives

and

The first can easily be obtained from Equation 4:

For the second, reindexing Equation 4 gives

Thus

And for the bias,

Finally, reindexing Equation 4 again,

so

At this point we have everything required to perform backpropagation through a feedforward neural network. However, the index forms of these equations are neither the easiest to read nor the most efficient to implement in code. We can vectorize them to obtain more compact forms.

Vectorized Backprop

By substituting Equation 9 into Equation 5 we obtain

In matrix form this is

where ⊙ is the Hadamard product.

Next, substituting Equation 10 into Equation 6 gives

Vectorized:

Substituting Equation 11 into Equation 7 gives

Vectorized:

since matrix multiplication with a vector of 1s is equivalent to summing each row of the matrix.

Substituting Equation 12 into Equation 8 leads to

Vectorized:

There we go! These four matrix equations allow backprop to be implemented more or less directly in NumPy using efficient linear algebra operations.

Figure 7: Backpropagation through a layer.

The MNIST Dataset

Now that we know how to train a neural network, let’s take a look at the data we will train it on.

Get Perrin Alaine-Sedano’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

The dataset consists of grayscale images of handwritten digits, 28x28 pixels in size. Each image is flattened into a vector of 784 numbers in the range 0–255. The training data is provided as a CSV file, which we load into a pandas DataFrame with 785 columns. The first column, the label, indicates the correct digit for the image. The following columns contain the pixel values. The test data is stored similarly but omits the labels. Let’s take a look.

train_data = pd.read_csv('/kaggle/input/competitions/digit-recognizer/train.csv')
train_data.head()

Figure 8: MNIST training data.

train_data.shape

(42000, 785)

test_data = pd.read_csv('/kaggle/input/competitions/digit-recognizer/test.csv')
test_data.head()

Figure 9: MNIST test data.

test_data.shape

(28000, 784)

Architecture of the Network

With knowledge of the MNIST dataset, we can design a neural network for the digit classification task. I chose to create a network with one hidden layer. It takes 784 inputs since this is the number of pixels in each image. The hidden layer contains 4096 neurons, and its activation function is the Rectified Linear Unit or ReLU function, defined by f(z) = max (0, z). Why this function in particular? A simple answer is that it’s very cheap to compute, and has been shown to work well in practice. Another reason is that in order for the Universal Approximation Theorem to apply, the hidden layer must use a non-polynomial activation function and ReLU fits the bill.

Figure 10: Graph of the ReLU function.

The network has 10 neurons in its output layer, one for each of the digits 0–9. It will be trained to output a probability distribution over digits: the output of the first neuron corresponds to digit 0, the output of the second to digit 1, and so on. The network’s classification for an image will be the digit with the highest estimated probability.

In order for the outputs to form a valid probability distribution, each one must be nonnegative and the sum of all the outputs must equal 1. The output layer will use the softmax activation function to ensure this.

We now need to define a loss function to see how well our network is performing. For this we will use something called cross-entropy loss. It provides an indication of how closely two probability distributions match. To calculate the loss for an individual training sample, we first pass it through the network and obtain the output yᴸ representing the predicted probability distribution. We need to compare it to a target (correct) probability distribution. What should this be? Easy: since we know with certainty the correct digit for each training sample, the target distribution should assign probability 1 to the correct digit and probability 0 to all other digits. This is called a one-hot representation. Each digit label in train_data will be converted to one-hot, as shown in the diagram below.

Figure 11: One-hot representation of digits.

The following diagram summarizes the architecture of the neural network.

Figure 12: Architecture of the MNIST classifier.

The two parts which are still unexplained are the softmax function and cross-entropy loss. Let’s take a look at them now.

Softmax

In the post on logistic regression, we used the sigmoid function to map the real-valued output of a linear classifier to a probability. A single number was sufficient for binary classification, but now we need to map a vector of 10 real numbers in the range (−∞, ∞) into probabilities for each class. We want the relative sizes of the network’s outputs to be preserved. That is, if the unactivated output for class i is larger than that for class j, then the probability of class i should be larger than that of class j.

How can we achieve these objectives? Let’s work backwards starting from a vector of class probabilities p = [p₁, p₂, …, pₙ] such that

For any pair of probabilities pᵢ and pⱼ, we can compute the odds pᵢ/pⱼ. We take the logarithm of this and obtain log ( pᵢ/pⱼ) = log pᵢ − log pⱼ.

Now we will define a vector z = [z₁, z₂, …, zₙ] where for all i, j ∈ {1, 2, …, n},

The elements of this vector are called logits in reference to the term log-odds.

Obviously such a vector can be obtained by setting zᵢ = log pᵢ, but that is not the only possibility. For an arbitrary constant C, zᵢ = log pᵢ + C also satisfies Equation 17.

This allows us to go from a vector of class probabilities p to a vector of logits z whose elements can range from −∞ to ∞. Equation 17 implies that if zᵢ > zⱼ, then pᵢ > pⱼ as we wanted. However, the choice of z is not unique since it depends on the arbitrary constant C.

Now imagine we are given only z without being told the value of C, and want to recover the probabilities. At first glance this might seem impossible. However, we can actually solve for C using the fact that the probabilities must sum to 1.

implies that

Thus

Therefore,

Equation 18 can be used to convert any vector z of real numbers into a vector p of probabilities. The probabilities will satisfy Equation 17, which ensures that the larger the value of zᵢ relative to the other logits, the larger the probability pᵢ will be.

In fact, Equation 18 is the definition of the softmax function:

Batched Softmax

So far we have discussed softmax applied to single vectors. However, it can be extended in a straightforward way to handle a batch of logit vectors.

For a batch Z of shape (num_classes, num_samples), we just need to compute softmax over every column (sample vector) separately. Different columns have no bearing on each other. Thus the formula becomes

Cross-Entropy Loss

Now we can move on to explaining the loss function. Let the network’s output for a single training sample be yᴸ (a probability distribution over digits) and the target be t (the target distribution, in this case a one-hot vector). The cross-entropy loss for that sample is given by

For a batch, we simply take the average over the training samples:

where N is the number of samples in the batch.

Surprise!

If you are not already familiar with the concepts, this formula probably doesn’t mean much. I said earlier that it relates to the similarity between two probability distributions. Now I want to provide some motivation for the formula and share an interesting perspective in the process.

In probability and information theory, there is a quantity called surprise, which is defined as the negative log of probability, −log p. This formula matches up with our intuition in an important way. Generally speaking, an event is more surprising if it has a lower probability of occurring. For example, it would not be surprising to flip a fair coin and have it land on tails, but to do so ten times and have it come up tails every time would be very surprising. The quantity −log p behaves analogously since it grows as pshrinks. For p = 1, −log p = 0, and as p → 0, −log p → ∞.

Surprise also has a fascinating and surprising relationship to information. Imagine generating a length-k string of random bits. The probability of getting any specific bitstring is 2⁻ᵏ. The surprise associated with that bitstring is then −log 2⁻ᵏ, which equals k if we use a base-2 logarithm. Thus, surprise can be interpreted as a measurement of the number of bits of information received from an observation.

The entropy of a probability distribution p over some set of observations represented by a random variable X can be defined as the expected surprise for a sample drawn according to that distribution:

This means entropy can also be interpreted as the amount of information you would expect to gain from observing said sample.

Cross-entropy is defined as the expected surprise when an observation is drawn from one probability distribution p, but surprise is calculated using a different probability distribution q:

In the context of machine learning, x is the data, p is the true (target) probability distribution over the data (e.g. tᵢⱼ), and q is the model’s estimated probability distribution (e.g. yᵢⱼᴸ). Cross-entropy thus quantifies the expected surprise of the model according to its own distribution q(x) for a sample drawn from the true distribution p(x).

The more accurately the model can predict the data, i.e. the more similar q(x) is to p(x), the less surprising it will find the observed data and the lower cross-entropy will be. Another interpretation is that the more knowledge the model already possesses, the less new information it will gain from observing the data.

Hopefully this digression provided some insight into why minimizing cross-entropy is a useful objective in machine learning. Now we can move on to obtaining the gradient of the loss function.

Computing Gradients for Backprop

In order to get started with backpropagation, we need the gradient of the loss w.r.t. the output of the neural network. Rather than using the softmax output Yᴸ, it will be more convenient to consider the gradient of the loss w.r.t. the unactivated output Zᴸ. This has a particularly simple form, as we will see.

Derivatives of Loss w.r.t. Yᴸ

By Equation 20,

Note that

if m = i and n = j. Otherwise it is 0. Thus all but one of the summation terms vanish and we are left with

Softmax Derivatives

We also need the partial derivatives of the softmax activation function.

Recall Equation 19. Since Yᴸ = softmax(Zᴸ),

Note that

since if i = m the derivative is equal to e^(zᵢⱼᴸ), and it’s 0 otherwise.

Thus

Final Loss Gradient

Now we can compute the gradient of the loss w.r.t. the logits via the multivariable chain rule.

By Equation 22,

vanishes unless n = j, so the above reduces to

By Equation 21,

so

Via ∑ⱼδᵢⱼxⱼ = xᵢ this becomes

Since each column of T is a one-hot vector, ∑ₘtₘⱼ = 1. Thus the above reduces to

In vectorized form, the gradient of the loss w.r.t. the logits is

An incredibly simple result! It’s just the difference between the predicted probabilities and the target probabilities divided by the number of samples.

ReLU Derivative

The last piece we need for backprop is the derivative of the ReLU activation function which is used in the hidden layer. For this we will use

Observant readers may notice that the derivative of ReLU is actually undefined at z = 0 since there is a sharp “corner” in the graph where the slope abruptly transitions from 0 to 1. However, in practice the derivative is set to 0 at this location as above.

Implementation

At last we have all the ingredients, and it’s time to code the neural network.

It will use a modular design informed by Figure 12. We will need several different types of layers.

The first type will compute the unactivated output of a fully-connected layer, i.e. Equation 1. This corresponds to the blue rectangles in Figure 12. It will be called Linear because it implements a linear function of the input, and also to follow PyTorch’s naming convention.

rng = np.random.default_rng()


class Linear:
    """
    A neural network layer which computes a linear function of the input using weights 
    and biases that are trainable by backpropagation with gradient descent.

    Parameters
    ----------
    n_inputs : int
        Number of inputs to the layer.
    n_outputs : int
        Number of outputs from the layer.
    """

    def __init__(self, n_inputs, n_outputs):
        # Weights are Kaiming initialized, drawn from a normal distribution with mean 0 and
        # variance 2 / n_inputs.
        self.weights = rng.normal(0, np.sqrt(2 / n_inputs), (n_outputs, n_inputs))
        self.bias = np.zeros((n_outputs, 1))
        self._last_input = None

As you can see, a Linear layer is defined by the number of inputs and number of outputs. It then creates a weight matrix and bias vector with the corresponding shapes. The weights are initialized using a method called Kaiming initialization which has been shown to improve the performance of neural networks using the ReLU activation function.

The class has a forward method which performs the computation of Equation 1.

    def forward(self, x):
        """
        Computes the output and saves the input for a future backward() call.
        
        Parameters
        ----------
        x : ndarray
            Input array.

        Returns
        -------
        ndarray
            Output array.
        """
        self._last_input = x
        
        return self.weights @ x + self.bias

Note that NumPy’s array broadcasting rules render the multiplication of the bias by 1ᵀ unnecessary as long as the bias has shape (n_outputs, 1).

The class also has a backward method which takes a gradient as input. Using this and self._last_input stored from the most recent forwardcall, it computes grad_wrt_weights, grad_wrt_bias, and grad_wrt_input via Equation 14, Equation 15, and Equation 16 respectively. It then performs a step of gradient descent by subtracting grad_wrt_weights and grad_wrt_bias from the weights and biases, nudging them in a direction that reduces loss for the training batch.

    def backward(self, gradient):
        """
        Updates weights and biases via gradient descent and returns a gradient for
        backpropagation to a previous layer.

        Parameters
        ----------
        gradient : ndarray of float
            Gradient of loss w.r.t. the output of this layer.

        Returns
        -------
        ndarray of float
            Gradient of loss w.r.t. the input to this layer.
        """
        if self._last_input is None:
            raise RuntimeError("forward(input) must be called before backward(gradient)")

        grad_wrt_weights = gradient @ self._last_input.T
        grad_wrt_bias = np.sum(gradient, axis=1, keepdims=True)
        grad_wrt_input = self.weights.T @ gradient
        
        self.weights -= grad_wrt_weights
        self.bias -= grad_wrt_bias
        
        return grad_wrt_input

Note that np.sum(gradient, axis=1, keepdims=True) produces an equivalent result to Equation 15. Due to keepdims=True, the shape of grad_wrt_bias is (n_outputs, 1) rather than just (n_outputs,). This means grad_wrt_bias has the same shape as self.bias, which is necessary for self.bias -= grad_wrt_bias to work as intended.

Next we will create a ReLU layer, which computes Equation 2.

class ReLU:
    """
    A neural network layer which applies the ReLU activation function to its input.
    """
    
    def __init__(self):
        self._last_input = None

    def forward(self, x):
        """
        Computes the output and saves its input for a future backward() call.

        Parameters
        ----------
        x : ndarray
            Input array.

        Returns
        -------
        ndarray
            Result of applying ReLU to each element of the input array.
        """
        self._last_input = x
        
        return np.maximum(0, x)

    def backward(self, gradient):
        """
        Computes the gradient of loss w.r.t. this layer's input.

        Parameters
        ----------
        gradient : ndarray of float
            Gradient of loss w.r.t. the output of this layer.

        Returns
        -------
        ndarray of float
            Gradient of loss w.r.t. the input to this layer.
        """
        if self._last_input is None:
            raise RuntimeError("forward(input) must be called before backward(gradient)")

        relu_deriv = (self._last_input > 0).astype(float)
        
        return gradient * relu_deriv

Like the Linear layer, it has forward and backward methods, and stores the input to forward for later use in backward. backward uses Equation 13 and Equation 24 to compute the backpropagated gradient.

Together, Linear and ReLU can be combined to form a fully-connected layer. You may be wondering why I didn’t just combine them into a single class. This was my original design, but splitting them up has several advantages. It simplifies the gradient computations in each class, and enables a Linear layer to easily be coupled with different activation functions. It is also the method used by PyTorch.

Next we implement the softmax function:

def softmax(logits, axis=0):
    """
    Parameters
    ----------
    logits : ndarray of float
        Array of logits.
    axis : int
        The axis over which to apply softmax.

    Returns
    -------
    ndarray of float
        An array of probabilities, the result of applying softmax to the logits over the
        specified axis.
    """
    shifted = logits - np.max(logits, axis=axis, keepdims=True)
    exp = np.exp(shifted)
    
    return exp / np.sum(exp, axis=axis, keepdims=True)

Subtracting the maximum logit from all the logits with shifted = logits - np.max(logits, axis=axis, keepdims=True) does not alter the resulting probabilities, and ensures that the largest value in the array is 0. This prevents numerical overflow caused by exponentiating a large number.

Next we can implement cross-entropy loss.

class CELoss:
    """
    Cross-entropy loss on softmax(logits), averaged over samples.
    """

    def __init__(self):
        self._last_target = None
        self._last_probs = None

Its forward method takes in raw logits and the target output, computes probabilities from the logits using softmax, and then computes the loss via Equation 20.

    def forward(self, logits, target):
        """
        Computes the cross-entropy loss for the given logits and targets. Saves targets and
        probabilities for a future backward() call.

        Parameters
        ----------
        logits : ndarray of float of shape (num_classes, num_samples)
            Array of logits representing a probability distribution over classes for each
            sample.
        target : ndarray of shape (num_classes, num_samples)
            The correct probability distributions over the classes, often one-hot
            (probability 1 for the correct class and 0 for the incorrect classes).

        Returns
        -------
        float
            The cross-entropy loss averaged over the samples.
        """
        probs = softmax(logits)

        self._last_target = target
        self._last_probs = probs

        num_samples = logits.shape[1]
        loss = -np.sum(target * np.log(probs)) / num_samples
        
        return loss

The backward method uses Equation 23 to get the gradient of the loss w.r.t. the logits.

    def backward(self):
        """
        Computes the gradient of the loss w.r.t. the logits from the most recent forward()
        call.
        
        Returns
        -------
        ndarray of float of shape (num_classes, num_samples)
            Gradient of loss w.r.t. the logits.
        """
        if self._last_probs is None:
            raise RuntimeError(
                "forward(logits, target) must be called before backward()"
            )
            
        num_samples = self._last_probs.shape[1]
        
        return (self._last_probs - self._last_target) / num_samples

Now we can create a class to represent a feedforward neural network. I have tried to make it maximally flexible, allowing it to contain any number of layers of any type as long as they have forward and backward methods.

class FeedforwardNN:
    """
    A feedforward neural network with customizable layers, loss function and learning rate.

    Parameters
    ----------
    layers : list of object
        A list of layers in the neural network. Each layer should have a forward method
        which computes its output and a backward method which computes gradients for
        backpropagation and updates the layer's parameters by gradient descent if
        applicable.
    loss : object
        The loss to use for the neural network. Requires a backward method to compute the
        gradient of loss w.r.t. the output of the neural network.
    learning_rate : float
        Parameter controlling how large an update is made to the network's parameters
        during each step of gradient descent.
    """
    
    def __init__(self, layers, loss, learning_rate):
        self.layers = layers
        self._loss = loss
        self.learning_rate = learning_rate
        self._output = None

    def forward(self, x):
        """
        Passes input through all the layers of the network, computing the final output.
        Saves the output for a future loss() call.

        Parameters
        ----------
        x : ndarray
            Array of input to the network.

        Returns
        -------
        ndarray
            The output of the network for the given input.
        """
        self._output = x
        
        for layer in self.layers:
            self._output = layer.forward(self._output)
        
        return self._output

    def loss(self, target):
        """
        Computes the loss from the last forward() call.

        Parameters
        ----------
        target : ndarray
            The correct outputs which the network should have returned.

        Returns
        -------
        float
            Loss for the most recent forward() call.
        """
        if self._output is None:
            raise RuntimeError("forward(input) must be called before loss(target)")
            
        return self._loss.forward(self._output, target)

    def backward(self):
        """
        Performs backpropagation with gradient descent to train the neural network.
        Will only work if forward() and loss() have already been called.
        """
        gradient = self.learning_rate * self._loss.backward()

        for layer in reversed(self.layers):
            gradient = layer.backward(gradient)

Note that the backward method scales the gradient by self.learning_rate before passing it back through the layers. This allows the learning rate to control the size of the gradient descent steps each Linear layer takes. A more advanced implementation could enable per-layer tuning of the learning rate.

Training the Neural Network

Let’s write a function that evaluates a network’s performance on a set of samples and targets.

def evaluate(model, X, Y):
    """
    Gets mean loss and accuracy over (X, Y), without updating the model.

    Parameters
    ----------
    model : FeedforwardNN
        The neural network.
    X : ndarray of float of shape (num_samples, num_inputs)
        The array of inputs to be used for the evaluation.
    Y : ndarray of shape (num_samples, num_classes)
        The array of targets (correct classes) for the corresponding inputs.

    Returns
    -------
    tuple
        The model's average loss and proportion of correct answers.
    """
    num_samples = X.shape[0]

    # transpose to comply with input format expected by the model.
    X = X.T # (num_samples, num_inputs) -> (num_inputs, num_samples)
    Y = Y.T # (num_samples, num_classes) -> (num_classes, num_samples)
    
    logits = model.forward(X)
    loss = model.loss(Y)
    num_correct = np.sum(np.argmax(logits, axis=0) == np.argmax(Y, axis=0))

    return loss, num_correct / num_samples

Next we can write the core training loop.

def train(model, X, Y, epochs, batch_size, X_eval, Y_eval):
    """
    Trains a neural network on the given inputs and targets for a certain number of epochs.
    Reports training and evaluation data for the training run.   

    Parameters
    ----------
    model : FeedforwardNN
        The neural network to be trained.
    X : ndarray of float of shape (num_samples, num_inputs)
        The array of inputs to be used for training.
    Y : ndarray of shape (num_samples, num_classes)
        The array of targets (correct classes) for the corresponding inputs.
    epochs : int
        Number of training epochs.
    batch_size : int
        How many inputs to feed the network at once during training.
    X_eval : ndarray of float of shape (num_samples, num_inputs)
        An array of inputs to be used for evaluation, separate from the training set.
    Y_eval : ndarray of shape (num_samples, num_classes)
        The array of targets for the corresponding evaluation inputs.

    Returns
    -------
    dict
        A history dict of per-epoch metrics. Epoch 0 performance is measured before
        training. After that, training metrics are averaged over each epoch's batches as
        the weights update, while eval metrics are measured once the epoch has finished.
    """
    num_samples = X.shape[0]
    history = {"train_loss": [], "train_acc": [], "eval_loss": [], "eval_acc": []}

    def record(epoch, train_loss, train_acc):
        history["train_loss"].append(train_loss)
        history["train_acc"].append(train_acc)
        message = f"epoch {epoch}: train_loss={train_loss:.4f} train_acc={train_acc:.4f}"

        eval_loss, eval_acc = evaluate(model, X_eval, Y_eval)
        history["eval_loss"].append(eval_loss)
        history["eval_acc"].append(eval_acc)
        message += f" eval_loss={eval_loss:.4f} eval_acc={eval_acc:.4f}"

        print(message)

    record(0, *evaluate(model, X, Y))
    
    for epoch in range(1, epochs + 1):
        perm = rng.permutation(num_samples)
        X, Y = X[perm], Y[perm]

        total_loss = 0.0
        num_correct = 0
        
        for start in range(0, num_samples, batch_size):
            # transpose to comply with input format expected by the model.
            batch_X = X[start:start + batch_size].T
            batch_Y = Y[start:start + batch_size].T

            logits = model.forward(batch_X)
            total_loss += model.loss(batch_Y) * batch_X.shape[1]
            num_correct += np.sum(
                np.argmax(logits, axis=0) == np.argmax(batch_Y, axis=0)
            )
            model.backward()

        record(epoch, total_loss / num_samples, num_correct / num_samples)

    return history

The model is trained on training samples X and targets Y, and evaluated on X_eval and Y_eval. Evaluating on separate data is important because it checks whether the model is overfitting to training data.

The helper function record appends training loss and accuracy to the history dict and then runs evaluate on the model. Before any training has occurred the model’s performance is assessed with record(0, *evaluate(model, X, Y)). Then the main training loop runs: for each training epoch, X, Y = X[perm], Y[perm] shuffles the training data and targets. We step through the shuffled data in chunks of size batch_size(the last chunk will be shorter if the data length is not divisible by batch_size), run the model on each batch, compute the loss and accuracy, and train the model via model.backward(). At the end of each epoch, we call record, passing it the training loss and accuracy for that epoch.

We can now preprocess the data, define the model, and train it!

EPOCHS = 30

if __name__ == '__main__':
    # divide by 255 to normalize pixel values
    # makes training more numerically well-behaved
    X = train_data.drop('label', axis=1).values / 255.0
    labels = train_data['label'].values
    Y = np.eye(10)[labels] # convert labels to one-hot representation
    
    # Hold out 10% of the labelled data for evaluation
    perm = rng.permutation(X.shape[0])
    num_eval = X.shape[0] // 10
    X_eval, Y_eval = X[perm[:num_eval]], Y[perm[:num_eval]]
    X_train, Y_train = X[perm[num_eval:]], Y[perm[num_eval:]]
    
    model = FeedforwardNN(
        layers=[Linear(784, 4096), ReLU(), Linear(4096, 10)],
        loss=CELoss(),
        learning_rate=0.01,
    )
    
    history = train(
        model,
        X_train,
        Y_train,
        epochs=EPOCHS,
        batch_size=64,
        X_eval=X_eval,
        Y_eval=Y_eval
    )

The line Y = np.eye(10)[labels] deserves some extra explanation. np.eye(10) creates a 10x10 identity matrix. Then np.eye(10)[labels]takes the labels (a vector of digits), and uses each digit as an index for this identity matrix. The result is a new matrix whose rows are one-hot vectors for each digit, since the ith row of the identity matrix has a 1 in the ith position and zeros in all other positions.

A FeedforwardNN is created based on Figure 12. One difference you may notice is that the network’s output is the logits and the final softmax step occurs only within CELoss. This is fine for classification since the largest logit will correspond to the largest probability, so the predicted class is just the argmax of the logits. Hence the use of np.argmax in the evaluate and train functions.

Finally, the model is trained for 30 epochs with a batch size of 64.

Kaggle Submission

Below is the code used to submit the model’s predictions for the Digit Recognizer competition.

X_test = test_data.values / 255.0
# Transpose X_test to comply with model input shape.
predictions = model.forward(X_test.T).argmax(axis=0)
submission = pd.DataFrame({
    'ImageId': np.arange(1, len(predictions) + 1),
    'Label': predictions,
})
submission.to_csv('submission.csv', index=False)

Now let’s look at how well the network completes the task of classifying handwritten digits.

Results

Below is a summary of the model’s training and evaluation performance at intervals of 5 epochs.

epoch 0: train_loss=2.3755 train_acc=0.1045 eval_loss=2.3752 eval_acc=0.1029
epoch 5: train_loss=0.2543 train_acc=0.9306 eval_loss=0.2575 eval_acc=0.9274
epoch 10: train_loss=0.1968 train_acc=0.9470 eval_loss=0.2131 eval_acc=0.9374
epoch 15: train_loss=0.1636 train_acc=0.9560 eval_loss=0.1847 eval_acc=0.9457
epoch 20: train_loss=0.1397 train_acc=0.9639 eval_loss=0.1659 eval_acc=0.9502
epoch 25: train_loss=0.1213 train_acc=0.9699 eval_loss=0.1519 eval_acc=0.9538
epoch 30: train_loss=0.1066 train_acc=0.9737 eval_loss=0.1425 eval_acc=0.9579

Here are graphs of accuracy and loss for the same training run.

Figure 13: Average accuracy per epoch.

Figure 14: Average loss per epoch.

As you can see, the training accuracy is slightly higher than the evaluation accuracy, which is to be expected since the model is never trained on evaluation data. However, the discrepancy is small.

The model’s best competition score so far is 0.95857.

Conclusion

At last, we have come to the end of the journey. This has definitely been the longest, most technically dense, and most difficult post to write so far. I didn’t expect it to be such a beast, but it was worth it. After performing all the derivations for backprop by hand and writing an MNIST classifier from scratch, I have gained a much better understanding of neural network fundamentals. I hope that readers have learned a lot as well. Be sure to tune in next time!

The code is available on Kaggle.