Deep Neural Networks

These notes are incomplete. They merely provide a summary. Please consult Chapters 6, 7, and 8 of Goodfellow-et-al-2016 for more details.

Deep neural networks as universal function approximators

Deep neural networks (DNN) are function approximators expressing hierarchically layered information. You can use DNNs to approximate a function of $d$ inputs to $q$ outputs using some parameters $\theta$. We, typicaly write $\mathbf{y} = f(\mathbf{x};\theta)$. Here $f(\mathbf{x};\theta)$ is the DNN, and $\theta$ are its parameters. We will clarify both these concepts below.

Mathematically, deep neural networks are compositions of simpler one-layer neural networks:

$$ f(\mathbf{x};\theta) = (f_L \circ f_{L-1} \circ \cdots \circ f_1)( \mathbf{x}). $$

In the simplest setting, the layers $f_l$s are a composition of an elementwise nonlinearity with a linear transformation:

$$ f_i ( \mathbf{z} ) = h^{(i)} ( \mathbf{W}^{(i)} \mathbf{z} + \mathbf{b}^{(i)} ), $$

where, $\mathbf{W}^{(i)}$ is a matrix of parameters, $\mathbf{b}^{(i)}$ is a vector of parameters, and $h^{(i)}$ is a nonlinear function applied in an elementwise fashion (i.e., applied separately to each one of the inputs that are provided to it). A DNN with this structure is called a fully-connected DNN.

In deep learning terminology, the matrix $\mathbf{W}^{(i)}$ is referred to as a weight matrix, the vector $\mathbf{b}^{(i)}$ is referred to as a bias. The function $h^{(i)}(\cdot)$ is called the activation function. It is typical for all but the last layer of a DNN to have the same activation function.

At the final layer, the output's dimensionality and the activation function's choice are dictated by constraints on the final output of the function $f$. For example:

1. If the output from $f$ is a real number with no constraints, the output dimensions are $d_L=1$ and $h^{(L)}(\mathbf{z}) = 1$.
2. If the output from $f$ is a positive real, $n^{(L)} = q = 1$ and $\sigma_L(x) = \exp(x)$.
3. If the output from $f$ is a probability mass function on $K$ categories, $n^{(L)} = q = K$ and $h^{(L)}(\mathbf{z}) = \frac{\exp(z_i)}{\sum_{j=1}^{K} \exp(z_j)}, i=1, 2, \dots, K$. Don't try to memorize this. We will revisit it in the following lecture.

Different ways of constructing the compositional structure of $f$ lead to different architectures such as fully connected networks (shown above), recurrent neural networks, convolutional neural networks, autoencoders, residual networks, etc.

Activation functions

The most common activation functions include the rectified Linear Units or ReLU (and variants), sigmoid functions, hyperbolic tangents, sinusoids, step functions, etc. We will visualize them in the hands-on activity.

Universal theorem for neural networks

The universal approximation theorem guarantees that DNNs are good function approximators. In plain English, the (original) theorem states that if you take any decent activation function and build a dense neural network with which you can approximate any continuous function (defined on a compact input domain) arbitrarily well, if you keep increasing the number of neurons you use. Recently, researchers have proven similar theorems for deep neural networks. In general, if you grow your network by adding neurons and layers, it can approximate almost anything you need. That's one of the reasons why deep neural networks have been (re)gaining momentum recently. In practice, however, it is a bit more difficult than that.

Training regression networks - Loss function view

Assume that you want to solve a regression problem. You have input data:

$$ \mathbf{x}_{1:n} = (\mathbf{x}_1,\dots,\mathbf{x}_n), $$

and output data:

$$ \mathbf{y}_{1:n} = (y_1,\dots,y_n). $$

You want to use them to find the map between $\mathbf{x}$ and $y$ using DNNs.

Well, you start by using a DNN $y=f(\mathbf{x};\theta)$ to represent the map from $\mathbf{x}$ to $y$. Here $\theta$ are the network parameters (the layers' weights and biases). Your problem is to fit $\theta$ to the available data.

The simplest way forward is to follow a least-squares approach. First, define a so-called loss function:

$$ L(\theta) = \frac{1}{n}\sum_{i=1}^n\left(y_i-f(\mathbf{x}_i;\theta)\right)^2. $$

This loss function is the sum of the squares of the prediction error of the DNN for a given $\theta$. Once you have the loss function, you can fit $\theta$ by minimizing it:

$$ \theta^* = \arg\min L(\theta). $$

However, this minimization problem does not have an analytical solution. Neither does it have a unique solution. It is a non-linear, non-convex optimization problem. It requires special treatment. We will talk about it in a while.

Training regression networks - Probabilistic view

Sometimes, it is not straightforward how to come up with loss functions. In such situations, we can employ a probabilistic view. We must develop a likelihood function that helps us connect the model to the observed data. So, in general, we need to come up with: $ p(y_{1:n}|\mathbf{x}_{1:n},\theta). $ Then, we can fit the parameters by maximizing the log-likelihood, which is the same as minimizing the "loss" function:

$$ L(\theta) = -\log p(y_{1:n}|\mathbf{x}_{1:n},\theta). $$

This approach is going to give you the same thing as the loss function approach under the following assumptions:

• The observations are independent (conditional on the model); and
• The measurement noise is Gaussian with a mean given by the DNN and a constant variance.

Let's show this. Take:

\begin{split} p(y_i|\mathbf{x}_i,\theta) &= N(y_i | f(\mathbf{x}_i;\theta), \sigma^2)\\ &= \frac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\frac{\left(y_i-f(\mathbf{x}_i;\theta)\right)^2}{2\sigma^2}\right\}, \end{split}

where $\sigma^2$ is the measurement noise variance. Then, from independence, we have:

$$ p(y_{1:n}|\mathbf{x}_{1:n},\theta) = \prod_{i=1}^np(y_i|\mathbf{x}_i,\theta). $$

So, we should be minimizing:

\begin{split} L'(\theta) &= -\log p(y_{1:n}|\mathbf{x}_{1:n},\theta)\\ &= -\sum_{i=1}^n\log p(y_i|\mathbf{x}_i,\theta)\\ &= \frac{1}{2\sigma^2}\sum_{i=1}^n\left(y_i-f(\mathbf{x}_i;\theta)\right)^2 + \text{const}. \end{split}

That's the same (up to an additive constant) as the $L(\theta)$ we had before. The benefit of the probabilistic approach is that it allows you to be more flexible with how you model the measurement process.

The minimization problem as a stochastic optimization problem

As mentioned, $L(\theta)$ is non-linear and non-convex. Classic, gradient-based optimization techniques do not work well on it. They tend to get trapped in bad local minima. Adding stochasticity to the optimization algorithm helps it avoid these bad local minima. Such stochastic optimization algorithms are still finding local minima, but they are better ones!

Another potential problem is that $L(\theta)$ may involve a summation over millions of observations (in the case of big data). In this regime, gradient-based optimization algorithms are also computationally inefficient. Stochastic optimization algorithms subsample the available data, allowing you to break them down into computationally digestible batches.

Let's first say what a stochastic optimization problem is. Then, we will show how to recast a typical $\min L(\theta)$ problem as a stochastic optimization problem. A stochastic optimization problem is a problem of the form:

$$ \min_\theta \mathbb{E}_Z[\ell(\theta;Z)], $$

where $\ell(\theta;Z)$ is some scalar function of $\theta$ and the random vector $Z$. The expectation is over $Z$. You want to minimize the expectation over $Z$ of $\ell(\theta;Z)$. That's it.

Okay. Back to our original problem. Take:

$$ L(\theta) = \frac{1}{n}\sum_{i=1}^n \left(y_i-f(\mathbf{x}_i;\theta)\right)^2. $$

We need to write this as an expectation of something. But an expectation of what? Well, it will be an expectation over randomly selected batches of the observed data. This is by no means the only choice. But it is a very useful choice. Let's see how we can do this.

First, let's visit the observations one by one. Take $I$ to be a Categorical random variable that picks with equal probability the index of one of the $n$ observations, i.e.,

$$ I \sim \operatorname{Categorical}\left(\frac{1}{n},\dots,\frac{1}{n}\right). $$

Take:

$$ \ell(\theta;I) = \left(y_I-f(\mathbf{x}_I;\theta)\right)^2. $$

So, here $Z = I$. Let's take the expectation over $I$ and see what it is going to give us:

\begin{split} \mathbb{E}_I[\ell(\theta;I)] &= \sum_{i=1}^np(I=i)\ell(\theta;i)\\ &= \sum_{i=1}^n\frac{1}{n}\left(y_i-f(\mathbf{x}_i;\theta)\right)^2\\ &= \frac{1}{n}\sum_{i=1}^n\left(y_i-f(\mathbf{x}_i;\theta)\right)^2\\ &= L(\theta). \end{split}

Great! Minimizing $L(\theta)$ is the same as minimizing $\mathbb{E}_I[\ell(\theta;I)]$.

Let's now do it again, using an $m$-sized randomly selected batch from the observed data. Take $I_1,I_2,\dots,I_m$ to be independent and identically distributed Categoricals that choose with equal probability an index from 1 to $n$. Then define:

$$ \ell_m(\theta;I_{1:m}) = \frac{1}{m}\sum_{j=1}^m\left(y_{I_j}-f(\mathbf{x}_{I_j};\theta)\right)^2. $$

So, here $Z = (I_1,\dots,I_m)$. Now take the expectation of this over the $I$'s:

\begin{split} \mathbb{E}[\ell_m(\theta;I_{1:m})] &= \mathbb{E}\left[\frac{1}{m}\sum_{j=1}^m\left(y_{I_j}-f(\mathbf{x}_{I_j};\theta)\right)^2\right]\\ &= \frac{1}{m}\sum_{j=1}^m\mathbb{E}\left[\left(y_{I_j}-f(\mathbf{x}_{I_j};\theta)\right)^2\right]\\ &= \frac{1}{m}\sum_{j=1}^m L(\theta)\\ &= \frac{m}{m}L(\theta)\\ &= L(\theta), \end{split}

where we have used that $\mathbb{E}\left[\left(y_{I_j}-f(\mathbf{x}_{I_j};\theta)\right)^2\right] = L(\theta)$ since it follows from our previous analysis. Therefore, minimizing $L(\theta)$ is the same as minimizing the expectation of $\ell_m(\theta;I_{1:m})$.

The Robbins-Monro algorithm

We reached the point where we can discuss the simplest variant of a stochastic optimization algorithm. It is the stochastic gradient descent or the Robbins-Monro algorithm, {cite:p}robbins-monro-1951. It goes as follows. Take the stochastic optimization problem:

$$ \min_\theta \mathbb{E}_Z[\ell(\theta;Z)]. $$

and consider the RM algorithm:

• initialize $\theta$ to $\theta_0$
• Iterate:

$$ \theta_{t+1} = \theta_t - \alpha_t \nabla_{\theta}\ell(\theta_t,z_t), $$

where $z_t$ are independent samples of $Z$.

In this algorithm, $\theta_t$ is gradually evolved following a noisy gradient signal. The sequence $\alpha_t$ is our choice, known as the learning rate. The Robbins-Monro theorem guarantees that the RM algorithm converges to a local minimum of the expectation $\mathbb{E}[\ell(\theta, Z)]$ if the learning rate satisfies the following properties:

$$ \sum_{t=1}^\infty \alpha_t = +\infty, $$

and

$$ \sum_{t=1}^\infty \alpha_t^2 < +\infty. $$

Intuitively, these properties say that the learning rate should converge to zero (an implication of the convergence of the second series) but not too fast (an implication of the divergence of the first series). Many sequences of learning rates satisfy these constraints. Here is a very commonly used one:

$$ \alpha_t = \frac{A}{(Bt + C)^\rho}, $$

with $\rho$ a number between $0.5$ and $1$ (exclusive).

Application of the Robbins-Monro algorithm to training regression networks

The algorithm for training regression networks becomes:

$$ \theta_{t+1} = \theta_t - \alpha_t\nabla_{\theta} \frac{1}{m}\sum_{j=1}^m\left(y_{i_{tj}}-f(\mathbf{x}_{i_{tj}};\theta_t)\right)^2, $$

where $i_{t1},\dots,i_{tm}$ are randomly selected indices of the observation data. Using properties of the gradient, you can also write this as:

$$ \theta_{t+1} = \theta_t - 2\alpha_t \frac{1}{m}\sum_{j=1}^m\left(y_{i_{tj}}-f(\mathbf{x}_{i_{tj}};\theta_t)\right)\nabla_{\theta}f(\mathbf{x}_{i_{tj}};\theta_t). $$

That's it.

Notice that to carry out the algorithm, we need to $\nabla_{\theta}f(\mathbf{x}_{i_{tj}};\theta_t)$, i.e., the gradient of the neural network output with respect to the parameters (weights and biases). This is done using the chain rule. The algorithm is known as the back-propagation algorithm. We are not going to cover it. Nowadays, you don't have to worry about derivatives. Software like PyTorch and JAX can find the derivatives for you. In the hands-on activity, I will introduce you to PyTorch.

Advanced variations of stochastic gradient descent

The RM algorithm is the simplest stochastic optimization algorithm I could explain in a lecture. It works, but it is only one of the most commonly used. More powerful algorithms like stochastic gradient descent with momentum, AdaGrad, or Adam (adaptive moment estimation) exist. I will show you in the hands-on activities how to use these algorithms as implemented in PyTorch, but I will not explain their details. If you want to know the details, please read Chapter 8 of {cite:p}Goodfellow-et-al-2016.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns
sns.set_context("paper")
sns.set_style("ticks")

import urllib.request
import os

def download(
    url : str,
    local_filename : str = None
):
    """Download a file from a url.
    
    Arguments
    url            -- The url we want to download.
    local_filename -- The filemame to write on. If not
                      specified 
    """
    if local_filename is None:
        local_filename = os.path.basename(url)
    urllib.request.urlretrieve(url, local_filename)

Regression with Deep Neural Networks

We want to understand the basics of PyTorch and set up and train regression DNNs with it. Some useful references on PyTorch are:

• Deep Learning with PyTorch: A 60-minute blitz and in particular:
  • What is PyTorch?
  • Autograd: Automatic differentation
  • Neural networks

What is PyTorch, and why are we using it?

• PyTorch is an alternative to Numpy that can harness the power of GPUs.
• PyTorch provides some core functionality for Neural Networks:
  • Some essential elements for building them up, like linear layers, activation functions, etc.
  • Automatic differentiation for getting the derivative of loss functions with respect to parameters.
  • Some stochastic optimization algorithms for minimizing loss functions
  • ...

I am not going to provide here a complete tutorial on PyTorch. You are advised to go over the first three topics of the Deep Learning with PyTorch: A 60-minute blitz before beginning this hands-on activity. Otherwise, it is unlikely that you understand the code that follows.

The Facebook AI Research Group (now Meta AI Research Group, I guess) developed PyTorch. Another powerful alternative Google Brain developed is TensorFlow. However, TensorFlow is gradually losing popularity. Instead of TensorFlow, Google Brain is now developing JAX. There are pros and cons to each of these frameworks. In research, I use either PyTorch or JAX depending on what we are doing. It is worth learning both of them, but we will stick to PyTorch in this course.

Making neural networks in PyTorch

PyTorch is pretty flexible in allowing you to make any neural network you like. You have absolute freedom on how your model looks like. However, it does provide a super easy way to make dense neural networks with a fixed activation function. That's what we are going to start with. First, import torch:

import torch

The submodule torch.nn is where the neural network building blocks reside:

import torch.nn as nn

First, let me show you how you can make a single linear layer:

$$ y = Wx + b. $$

The weights are selected randomly if not specified. Here you go:

layer = nn.Linear(1, 20)

This function now takes one-dimensional inputs and spits out 20-dimensional outputs. Here is how it works:

x = torch.rand(10, 1) # 10 randomly sampled one dimensinal inputs
print(x)

y = layer(x)
print(y)

print(y.shape)

So, this took us to 10, 20 dimensional outputs. Looks good.

But where are the weights and the bias term? Here they are:

layer.weight

layer.bias

You can directly change them if you wish. Notice the requires_grad=True flag. This is because PyTorch knows that these are parameters to be optimized.

There is a little bit of flexibility on nn.Linear. For example, you can altogether drop the bias if you wish. For the complete list of possibilities, you should always check the docs.

Now, let's get to the activation functions. There are a lot already in torch.nn. Here is the sigmoid:

h = nn.Sigmoid()

fig, ax = plt.subplots()
x = torch.linspace(-5, 5, 100)[:, None]
ax.plot(x, h(x))
ax.set_xlabel('$x$')
ax.set_ylabel('$z=h(x)$')
ax.set_title('Activation function: ' + str(h))
sns.despine(trim=True);

You could also implement the activation function by hand. The only restriction is that you should use' PyTorch' instead of numpy functions. Here is how we would do it for the sigmoid:

# Here is how you could do this by hand:
h_by_hand = lambda x: torch.exp(x) / (1.0 + torch.exp(x))

fig, ax = plt.subplots(dpi=100)
ax.plot(x, h_by_hand(x))
ax.set_xlabel('$x$')
ax.set_ylabel('$z=h(x)$')
ax.set_title('Activation function: Sigmoid (hand version)')
sns.despine(trim=True);

Here are now some of the most commonly used activation functions in torch.nn:

fig, ax = plt.subplots(dpi=100)
ax.plot(x, h_by_hand(x))

for Func in [nn.Sigmoid, nn.ReLU, nn.Tanh]:
    h = Func()
    ax.plot(x, h(x), label=str(h))
    
ax.set_xlabel('$x$')
ax.set_ylabel('$z=h(x)$')
plt.legend(loc='best', frameon=False)
sns.despine(trim=True);

Now that we have a linear layer and an activation function, here is how we can combine them to make a function that takes us from the input to the internal neurons:

h = nn.Sigmoid()
z_func = lambda x: h(layer(x))

This is pretty much it. And that's now a function:

z_func(x)

Now, for regression, we would like to bring this back to a scalar output. We need to add one more linear layer to take the 20 internal neurons back to one dimension to do this.

final_layer = nn.Linear(20, 1)
f = lambda x: final_layer(z_func(x))
print(f(x).shape)

Instead of doing this manually, we can use the class nn.Sequential of PyTorch:

f = nn.Sequential(layer, nn.Sigmoid(), final_layer)

This is recommended because nn.Sequential adds some additional functionality, which I will show you in a while. You can evaluate this as a function, and you can also plot it. But to plot it, you have to turn the output into a proper numpy array. This is because matplotlib doesn't like PyTorch tensors that depend on parameters. Here is what you need to do:

y = f(x).detach().numpy() # detach freezes the parameters to whatever they are
                          # numpy returns a proper numpy array
print(type(y))
print(y.shape)

And here is what it looks like (remember the weights are random):

fig, ax = plt.subplots(dpi=100)
ax.plot(x, f(x).detach().numpy())
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
sns.despine(trim=True);

The class nn.Sequential is convenient because it allows us to build deep networks quickly. Here is a 5-layer network that starts from one input, takes us through 3 layers with 20 neurons each, and ends on a single output:

f = nn.Sequential(nn.Linear(1, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 1))

Where are the parameters in an object created in this way? Here they are:

for theta in f.named_parameters():
    print(theta)

And that's why we love PyTorch. Because it does all the dirty work for us. Imagine having to keep track of all these parameters by hand.

For those who want to know what is happening inside nn.Linear, note that it is a special case of a PyTorch neural network module; see nn.Module. You would directly inherit the latter when writing your own class for a non-standard neural network. We will not cover it in this class, but you can find plenty of examples here.

Making a loss function

Let's now make the loss function that we want to minimize. It needs to be a PyTorch function as well. For regression problems, we can think of the loss as a function of the model predictions and the observed data. That depends on the parameters that come through the predictions. In this form, let's write down the mean square error (MSE) loss. It is:

$$ L_{\text{MSE}}(\theta) = L_{\text{MSE}}(y_{1:n}, f(x_{1:n};\theta)) = \frac{1}{n}\sum_{i=1}^n\left[y_i-f(x_i;\theta)\right]^2, $$

where $x_{1:n}$ are the observed inputs (features), $y_{1:n}$ are the observed outputs (targets), and $f(x_{1:n};\theta)$ contains the model predictions on the observed inputs.

You can implement the MSE loss like this:

mse_loss_ours = lambda y, f: torch.mean((y - f) ** 2)

Or we can use built-in PyTorch functionality:

mse_loss = nn.MSELoss()

Let's evaluate it for some random data:

# The number of fake observations
n = 20
# Some fake observed features
x_fake = torch.rand(n, 1)
# Some fake observed targets
y_fake = 4 * x_fake ** 2 - 5 * x_fake ** 3 + 0.1 * torch.rand(n, 1)
fig, ax = plt.subplots(dpi=100)
ax.plot(x_fake, y_fake, 'x')
sns.despine(trim=True);

Here is how to calculate the loss (for the random parameters that our net started with):

# Predict with the net:
y_pred = f(x_fake)
# Evaluate the loss
our_loss = mse_loss_ours(y_fake, y_pred)
built_in_loss = mse_loss(y_fake, y_pred)
print(our_loss)
print(built_in_loss)

Let's minimize the MSE loss for these fake data and see what fit we will get. Here is how you do this in PyTorch. Since I don't have a lot of data, I will use gradient descent - no random subsampling of the data. I will show you how to use stochastic gradient descent in the following example.

# Reinitialize the net:
f = nn.Sequential(nn.Linear(1, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 1))

# Initialize the optimizer - Notice that it needs to know about the 
# parameters it is optimizing
optimizer = torch.optim.SGD(f.parameters(), lr=0.01) # lr is the learning rate
# Some place to hold the training loss for visualizing it later
training_loss = []
# Iterate the optimizer. Let's just do 10 iterations.
for i in range(10000):
    # This is essential for the optimizer to keep
    # track of the gradients correctly
    # It is using some buffers internally that need to
    # be manually zeroed on each iteration.
    # This is because it doesn't know when you are done with the
    # calculation of the loss
    optimizer.zero_grad()
    # Make predictions
    y_pred = f(x_fake)
    # Evaluate the loss - That's what you are minimizing
    loss = mse_loss(y_fake, y_pred)
    # Evaluate the derivative of the loss with respect to
    # all parameters - It knows how to do it because of
    # PyTorch magick
    loss.backward()
    # And now you are ready to make a step
    optimizer.step()
    # Save the training loss of later visualization
    training_loss.append(loss.item())
    # Print the loss every one hundend iterations:
    if i % 1000 == 0:
        print('it = {0:d}: loss = {1:1.3f}'.format(i, loss.item()))

Let's plot the predictions of this model on the fake data:

fig, ax = plt.subplots(dpi=100)
ax.plot(x_fake, y_fake, 'x');
xx = torch.linspace(x_fake.min(), x_fake.max(), 100)[:, None]
yy = f(xx).detach().numpy()
ax.plot(xx, yy)
sns.despine(trim=True);

This may or may not work, depending on what random seed you start with. It may stack at some local minimum if you run it a few times. This algorithm is excellent if we do stochastic optimization, i.e., subsampling the data. Here is how the loss changes with each iteration:

fig, ax = plt.subplots(dpi=100)
ax.plot(training_loss)
ax.set_xlabel('Iteration')
ax.set_ylabel('Training loss')
sns.despine(trim=True);

The problem is the plateau we have at the beginning of the optimization.

Let's redo this thing with stochastic optimization. For stochastic optimization, we need to subsample the data during each iteration. We can either do this manually or use the PyTorch functionality. First, let's do it manually.

# Pick a subsampling batch size
m = 5

# Reinitialize the net:
f = nn.Sequential(nn.Linear(1, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 1))

# Reinitialize the optimizer
optimizer = torch.optim.SGD(f.parameters(), lr=0.01)
# Keep track of the training loss
training_loss_sgd = []
# Iterate the optimizer. Let's just do 10 iterations.
for i in range(10000):
    # Zero out the gradient buffers
    optimizer.zero_grad()
    # Sample m observation indices at random
    idx = np.random.randint(0, n, m)
    # Here is the subsample of the data
    x_batch = x_fake[idx]
    y_batch = y_fake[idx]
    # Make predictions
    y_pred = f(x_batch)
    # Evaluate the loss - That's what you are minimizing
    loss = mse_loss(y_batch, y_pred)
    # Evaluate the derivative of the loss with respect to
    # all parameters - It knows how to do it because of
    # PyTorch magick
    loss.backward()
    # And now you are ready to make a step
    optimizer.step()
    # Keep track of the training loss
    training_loss_sgd.append(loss.item())
    # Print the loss every one hundend iterations:
    if i % 1000 == 0:
        print('it = {0:d}: loss = {1:1.2e}'.format(i, loss.item()))

fig, ax = plt.subplots(dpi=100)
ax.plot(x_fake, y_fake, 'x');
xx = torch.linspace(x_fake.min(), x_fake.max(), 100)[:, None]
yy = f(xx).detach().numpy()
ax.plot(xx, yy)
sns.despine(trim=True);

This fit does look a little bit better. Let's now compare the training loss of stochastic gradient descent to the previous one:

fig, ax = plt.subplots(dpi=100)
ax.plot(training_loss_sgd, label='Gradient descent')
ax.plot(training_loss, label='Stochastic gradient descent')
ax.set_xlabel('Iteration')
ax.set_ylabel('Training loss')
plt.legend(loc='best', frameon=False)
sns.despine(trim=True);

This wiggly nature of stochastic gradient descent allows it to escape bad local minima.

Questions

• Rerun the stochastic gradient descent with one sample per iteration ($m=1$). Does it converge? Do you need fewer or more iterations? Is it more or less wiggly?
• Rerun the stochastic gradient descent with ten samples per iteration. How does it perform now?

Example - Motorcycle Data

Let's now use the motorcycle data to do regression with DNNs. This will help us demonstrate some best practices specifically:

• Standardizing the data
• Splitting in training and test subsets

First, start by loading the dataset:

# Load the data
data = np.loadtxt('motor.dat')
x = data[:, 0][:, None]
y = data[:, 1][:, None]

# Split into training and test datasets
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3)

# Visualize them
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(x_train, y_train, 'x', markeredgewidth=2, label='Training data')
ax.plot(x_test, y_test, 'o', markeredgewidth=2, label='Test data')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.legend(loc='best', frameon=False)
sns.despine(trim=True);

# Turn the data into torch tensors:
x_train = torch.tensor(x_train, dtype=torch.float)
y_train = torch.tensor(y_train, dtype=torch.float)
x_test = torch.tensor(x_test, dtype=torch.float)
y_test = torch.tensor(y_test, dtype=torch.float)

Please note that the dtype=torch.float specification is absolutely needed here. You need to include it so that the code is going to work. The problem is that the x_train, etc., are all numpy arrays and that numpy arrays have 64-bit floating point numbers by default. PyTorch is using 32-bit floating point numbers by default. We need, at some point, to make the two compatible.

Now, we are ready to train the network. Let's give it a shot. We will use the same architecture as before. The only difference is that I will print the validation loss instead of the training loss.

# The number of training samples
n = x_train.shape[0]

# Pick a subsampling batch size
m = 5

# Reinitialize the net:
f = nn.Sequential(nn.Linear(1, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 1))

# Reinitialize the optimizer
optimizer = torch.optim.SGD(f.parameters(), lr=0.01)
# Keep track of the training loss and the test loss
training_loss = []
test_loss = []
# Iterate the optimizer. Let's just do 10 iterations.
for i in range(10000):
    # Zero out the gradient buffers
    optimizer.zero_grad()
    # Sample m observation indices at random
    idx = np.random.randint(0, n, m)
    # Here is the subsample of the data
    x_batch = x_train[idx]
    y_batch = y_train[idx]
    # Make predictions
    y_pred = f(x_batch)
    # Evaluate the loss - That's what you are minimizing
    loss = mse_loss(y_batch, y_pred)
    training_loss.append(loss.item())
    # Evaluate the derivative of the loss with respect to
    # all parameters - It knows how to do it because of
    # PyTorch magick
    loss.backward()
    # And now you are ready to make a step
    optimizer.step()
    # Evaluate the test loss
    y_pred_test = f(x_test)
    ts_loss = mse_loss(y_test, y_pred_test)
    test_loss.append(ts_loss.item())
    # Print the loss every one hundend iterations:
    if i % 1000 == 0:
        print('it = {0:d}: loss = {1:1.2e}'.format(i, ts_loss.item()))

The above code may not work at all, giving you nan's. Or it may work and get you nowhere. The problem here is the scale of both the inputs and the outputs and the assumptions that have been made about them when we initialize the weights of the net and when we pick the learning rate of the optimization algorithm. The easiest way to overcome this problem is to standardize the data. This is achieved by subtracting the empirical mean and dividing the inputs and the outputs by the empirical standard deviation. By standardizing the data, we are making the default parameters (for weight initialization and stochastic gradient descent) valid again. Standardization is such a common process that it is already implemented in sklearn.preprocessing.StandardScaler. Here is how it works:

from sklearn.preprocessing import StandardScaler

feature_scaler = StandardScaler().fit(x)
target_scaler = StandardScaler().fit(y)

The feature_scaler.transform() is a function:

$$ x \rightarrow \frac{x-\mu}{\sigma}, $$

where $\mu$ and $\sigma$ are the features' empirical mean and standard deviation. Here they are:

# The mean:
feature_scaler.mean_

# The standard deviation:
feature_scaler.scale_

And here is how the scalers work:

x_scaled = feature_scaler.transform(x)
y_scaled = target_scaler.transform(y)

The data are now scaled, see this fig:

fig, ax = plt.subplots(dpi=100)
ax.plot(x_scaled, y_scaled, 'x')
ax.set_xlabel('$x$ (scaled)')
ax.set_ylabel('$y$ (scaled)')
sns.despine(trim=True);

We will train the net using x_scale and y_scaled. We can always go back to the original scales at the end. Let's see if it works.

# Split in training and test
x_s_train, x_s_test, y_s_train, y_s_test = train_test_split(x_scaled, y_scaled,
                                                            test_size=0.3)

# Turn the data into torch tensors:
x_s_train = torch.tensor(x_s_train, dtype=torch.float)
y_s_train = torch.tensor(y_s_train, dtype=torch.float)
x_s_test = torch.tensor(x_s_test, dtype=torch.float)
y_s_test = torch.tensor(y_s_test, dtype=torch.float)

# The number of training samples
n = x_train.shape[0]

# Pick a subsampling batch size
m = 5

# Reinitialize the net:
f = nn.Sequential(nn.Linear(1, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 20),
                  nn.ReLU(),
                  nn.Linear(20, 1))

# Reinitialize the optimizer
optimizer = torch.optim.SGD(f.parameters(), lr=0.01)
# Keep track of the training loss and the test loss
training_loss = []
test_loss = []
# Iterate the optimizer. Let's just do 10 iterations.
for i in range(10000):
    # Zero out the gradient buffers
    optimizer.zero_grad()
    # Sample m observation indices at random
    idx = np.random.randint(0, n, m)
    # Here is the subsample of the data
    x_batch = x_s_train[idx]
    y_batch = y_s_train[idx]
    # Make predictions
    y_pred = f(x_batch)
    # Evaluate the loss - That's what you are minimizing
    loss = mse_loss(y_batch, y_pred)
    training_loss.append(loss.item())
    # Evaluate the derivative of the loss with respect to
    # all parameters - It knows how to do it because of
    # PyTorch magick
    loss.backward()
    # And now you are ready to make a step
    optimizer.step()
    # Evaluate the test loss
    y_pred_test = f(x_s_test)
    ts_loss = mse_loss(y_s_test, y_pred_test)
    test_loss.append(ts_loss.item())
    # Print the loss every one hundend iterations:
    if i % 1000 == 0:
        print('it = {0:d}: loss = {1:1.2e}'.format(i, ts_loss.item()))

Let's visualize the fit:

xx_scaled = torch.linspace(x_scaled.min(), x_scaled.max(), 100)[:, None]
yy_scaled = f(xx_scaled).detach().numpy()
fig, ax = plt.subplots(dpi=100)
ax.plot(x_s_train, y_s_train, 'x', label='Training data')
ax.plot(x_s_test, y_s_test, 'o', label='Test data')
ax.plot(xx_scaled, yy_scaled, label='DNN fit')
ax.set_xlabel('$x$ (scaled)')
ax.set_ylabel('$y$ (scaled)')
plt.legend(loc='best', frameon=False)
sns.despine(trim=True);

Here is predictions-observations plot on the test data set:

y_pred_test = f(x_s_test).detach().numpy()
fig, ax = plt.subplots(dpi=100)
ax.plot(y_pred_test, y_s_test, '.')
yys = np.linspace(y_s_test.min(), y_s_test.max(), 10)
ax.plot(yys, yys, 'r')
ax.set_xlabel('Predictions')
ax.set_ylabel('Observations')
sns.despine(trim=True);

Also, if you wish, you can scale the predictions back to the original units:

xx = feature_scaler.inverse_transform(xx_scaled)
yy = target_scaler.inverse_transform(yy_scaled)
fig, ax = plt.subplots(dpi=100)
ax.plot(x_train, y_train, 'x', label='Training data')
ax.plot(x_test, y_test, 'o', label='Test data')
ax.plot(xx, yy, label='DNN fit')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.legend(loc='best', frameon=False)
sns.despine(trim=True);

It is instructive to observe how the training and test losses evolve as a function of the optimization iteration:

fig, ax = plt.subplots(dpi=100)
ax.plot(training_loss, label='Training loss')
ax.plot(test_loss, label='Test loss')
ax.set_xlabel('Iteration')
ax.set_ylabel('Loss')
plt.legend(loc='best', frameon=False)
sns.despine(trim=True);

The wiggliness is, of course, due to the stochastic nature of the optimization. The training error converges to a minimum as you keep iterating. This is a direct consequence of the Robbins-Monro theorem. You will reach a local minimum of the training error eventually. However, this is not the case for the test error. The test error will reach a minimum eventually and then increase! It will always do this when training networks by minimizing a loss function. What happens is that the algorithm will start overfitting the training data and will not be able to generalize correctly for the test data. There are ways around this. We will learn about the basic one in the next lecture (weight regularization and early stopping). There are some advanced ways to avoid overfitting (e.g., dropout, Bayesian neural networks), which we will not cover in the class.

Questions

• Change the activation function from nn.ReLU to nn.Tanh. Are you getting a better fit or a worse fit?
• Rerun the code above for 100,000 iterations. Does it start to overfit the training data? What happens to the test loss?
• Rerun the code for 5,000 iterations. What does the prediction look like now? Early stopping would stop at about this point.

Deep Neural Networks Continued

Loss functions for classification

Take some features $\mathbf{x}_{1:n}$ and some discrete targets $y_{1:n}$. Because the targets are discrete, we have a classification problem. What loss function should we use? Let's examine two cases: binary and multiclass classification.

Binary classification

In binary classification, we use a DNN $f(\mathbf{x};\theta)$ with parameters $\theta$ to model the probability that $y$ take the value $1$ by:

$$ p(y=1|\mathbf{x},\theta) = \operatorname{sigm}(f(\mathbf{x};\theta)) = \frac{\exp\{f(\mathbf{x};\theta)\}}{1 + \exp\{f(\mathbf{x};\theta)\}}. $$

Remember that the sigmoid function takes the scalar $f(\mathbf{x};\theta)$ and maps it on $[0,1]$ so that we get a probability. From the obvious rule of probability, we get that:

$$ p(y=0|\mathbf{x},\theta) = 1 - p(y=1|\mathbf{x},\theta) = 1 - \operatorname{sigm}(f(\mathbf{x};\theta)). $$

So, for an arbitrary $y$ (either 0 or 1), we can write:

$$ p(y|\mathbf{x},\theta) = \left[\operatorname{sigm}(f(\mathbf{x};\theta))\right]^y \left[1-\operatorname{sigm}(f(\mathbf{x};\theta))\right]^{1-y}. $$

This is a nice trick because it activates the right term based on what $y$ is.

Now that we have specified the likelihood of a single observation, the likelihood of the entire dataset is:

$$ p(y_{1:n}|\mathbf{x}_{1:n},\theta) = \prod_{i=1}^n p(y_i|\mathbf{x},\theta). $$

We are almost done. The idea is to train the network by maximizing the log-likelihood, which is the same as minimizing the following loss function:

\begin{split} L(\theta) &= -\log p(y_{1:n}|\mathbf{x}_{1:n},\theta)\\ &= -\sum_{i=1}^n \left\{y_i \log \operatorname{sigm}(f(\mathbf{x}_i;\theta)) + (1-y_i)\log \left[1-\operatorname{sigm}(f(\mathbf{x}_i;\theta))\right] \right\}. \end{split}

This loss function is known as the cross entropy loss.

Multiclass classification

Now assume that $y$ can take $K$ different values ranging from $0$ to $K-1$. We need to model the probability that $y=k$ given $\mathbf{x}$. To do this, we introduce a DNN $\mathbf{f}(\mathbf{x};\theta)$ with parameters $\theta$ and $K$ outputs:

$$ \mathbf{f}(\mathbf{x};\theta) = \left(f_0(\mathbf{x};\theta),\dots,f_{K-1}(\mathbf{x};\theta)\right). $$

However, $\mathbf{f}(\mathbf{x})$ is just a bunch of $K$ scalars. We need to turn it into a $K$-dimensional probability vector. To achieve this, we define:

$$ p(y=k|\mathbf{x},\theta) = \operatorname{softmax}_k(\mathbf{f}(\mathbf{x};\theta)) := \frac{\exp\left\{f_k(\mathbf{x};\theta)\right\}}{\sum_{k'=0}^{K-1}\exp\left\{f_{k'}(\mathbf{x};\theta)\right\}}. $$

So, the role of the softmax is dule. First, to turn the scalars to positive numbers and, second, to normalize them.

For an arbitrary $y$, we can just write:

$$ p(y|\mathbf{x},\theta)) = \operatorname{softmax}_y(\mathbf{f}(\mathbf{x};\theta)). $$

So, the likelihood of the dataset is:

$$ p(y_{1:n}|\mathbf{x}_{1:n},\theta) = \prod_{i=1}^n p(y_i|\mathbf{x},\theta). $$

Therefore, the loss function we should be minimizing is:

\begin{split} L(\theta) &= -\log p(y_{1:n}|\mathbf{x}_{1:n},\theta)\\ &= -\sum_{i=1}^n \log \left[\operatorname{softmax}_{y_i}(\mathbf{f}(\mathbf{x}_i;\theta))\right]. \end{split}

This is also called cross entropy loss. Sometimes, it is called the softmax cross entropy loss. Now, all these names are not important. What is important is to understand how these losses are derived from maximum likelihood.

Regularization

DNNs are incredibly flexible. However, this makes them prone to overfitting. As we have discussed, one way to avoid overfitting is to add regularization to the loss function. This works as follows. Instead of minimizing $L(\theta)$ you minimize:

$$ J(\theta) = L(\theta) + \lambda R(\theta), $$

where $\lambda$ is a positive parameter, and $R(\theta)$ is a term that penalizes very big weights. One possibility of $R(\theta)$ is the L2 regularization:

$$ R(\theta) = \sum_j \theta_j^2. $$

But there are many more. Also, you could add multiple regularizers - not just one. The bigger $\lambda$ is, the more you regularize, and the smaller it is, the less you regularize. So, $\lambda$ is a parameter you would like to tune - but more on this later.

Bayesian interpretation of regularization

As we showed during our regression lectures, the regularization can be interpreted as a maximum posterior approach to inference. To see this, assume that we have a data likelihood $p(y_{1:n}|\mathbf{x}_{1:n},\theta)$ for either a regression or a classification problem. Now, assume that you have a prior over all parameters, say $p(\theta)$. Then, write down the posterior of the parameters:

$$ p(\theta|\mathbf{x}_{1:n}, y_{1:n}) \propto p(y_{1:n}|\mathbf{x}_{1:n},\theta)p(\theta). $$

Instead of picking the parameters by maximizing the log-likelihood, let's pick the parameters by maximizing the log posterior. This is equivalent to minimizing the following loss function:

$$ J(\theta) = -\log p(y_{1:n}|\mathbf{x}_{1:n},\theta) - \log p(\theta). $$

Of course, the first term is just $L(\theta)$. The second term is $\lambda R(\theta)$. For a more concrete example, assume that the prior over $\theta$ is a zero-mean Gaussian with precision $\alpha$:

$$ p(\theta) = N(\theta|0, \alpha^{-1}I). $$

Then, we get:

$$ -\log p(\theta) = \frac{\alpha}{2\pi}\sum_j\theta_j^2 + \text{const}. $$

And this is how you get the L2 regularization term in a principled way.

Convolution neural networks

Convolutional neural networks try to mimic the operations carried out by the animal visual cortex. By construction, they satisfy certain properties (e.g., invariance to small translations and rotations) that make them particularly successful in image recognition tasks. Explaining the mathematical details of convolutional neural networks is beyond the scope. The details presented in the lecture videos should be enough for this course and for building standard image classifiers. If you want to know more, read Chapter 9 of Goodfellow-et-al-2016.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns
sns.set_context("paper")
sns.set_style("ticks")

import urllib.request
import os

def download(
    url : str,
    local_filename : str = None
):
    """Download a file from a url.
    
    Arguments
    url            -- The url we want to download.
    local_filename -- The filemame to write on. If not
                      specified 
    """
    if local_filename is None:
        local_filename = os.path.basename(url)
    urllib.request.urlretrieve(url, local_filename)

Classification with Deep Neural Networks

We implement image classification network in PyTorch. We add L2 regularization, convolutional layers, and hyper parameter tuning. You may find these tutorials useful:

• Deep Learning with PyTorch: A 60 minute blitz and in particular:
  • Training a Classifier - with which we use the same dataset in this hands-on activity.

The CIFAR10 dataset

We will demonstrate multiclass classification using the CIFAR10 dateset. The dataset consists of 60000 32x32 color images in 10 classes (plane, car, bird, cat, deer, dog, frog, horse, ship, and truck), with 6000 images per class. The dataset can be downloaded directly from PyTorch using the torchvision module.

You can think of the original images as 32x32x3 arrays. The first two dimensions correspond to the pixels. The third dimension corresponds to the color (red, green, blue). Of course, we must turn them into PyTorch tensors. Also, it is more convenient to scale them to be between $[-1,1]$. We will achieve this using a transformation. Don't worry about this now. We will explain it as we go.

import torch
import torchvision
import torchvision.transforms as transforms

# This is the transformation that we will apply to each image
transform = transforms.Compose(
    [transforms.ToTensor(),   # This turns the picture to a Tensor
     transforms.Normalize((0.5, 0.5, 0.5), (0.5, 0.5, 0.5))]) # This scales it to [-1, 1]

# Here is how you can download the training dataset
trainset = torchvision.datasets.CIFAR10(root='./data', train=True,
                                        download=True, transform=transform)

# And here is how to download the test dataset:
testset = torchvision.datasets.CIFAR10(root='./data', train=False,
                                       download=True, transform=transform)

# These are the class labels
classes = ('plane', 'car', 'bird', 'cat',
           'deer', 'dog', 'frog', 'horse', 'ship', 'truck')

Now, all these data went into the folder "./data." Here is what this folder contains:

!ls ./data/

The file cifar-10-python.tar.gz is a compressed file containing everything. The contents were automatically extracted and put in the folder cifar-10-batches-py. Let's look inside this folder:

!ls -lht data/cifar-10-batches-py

You see several files. The important ones are data_batch_1 to data_batch_5 and test_batch. Each of these contains 10,000 images in a binary format. The format is explained here. We can read them as follows:

def unpickle(file):
    import pickle
    with open(file, 'rb') as fo:
        dict = pickle.load(fo, encoding='bytes')
    return dict

data = unpickle('data/cifar-10-batches-py/data_batch_1')
# data is a dictionary
# Here are the keys
print(data.keys())

# One key has to do with the pictures
# It gives you a numpy array:
print(data[b'data'].shape)

# The first dimension correspond to differnt picture
# The second dimension is
32 * 32 * 3

# So this is the first picture:
img = data[b'data'][0, :].reshape((32, 32, 3), order='F')
# Here is the Red channel:
print(img[:, :, 0])

The numbers go from 0 (no red) to 255 (full red). Here is how to visualize it:

fig, ax = plt.subplots(figsize=(1, 1))
ax.imshow(np.transpose(img, (1, 0, 2)));

This is clearly a frog. Let's verify this:

classes[data[b'labels'][0]]

This is nice. And we could proceed manually like this. However, PyTorch offers some helpful functionality. Let's investigate the trainset that was returned by CIFAR10:

trainset

# Here are the classes:
trainset.classes

# Here is the correspondence between classes and discrete labels
trainset.class_to_idx

# Here are the images from all training batches
print(trainset.data.shape)

# Here are the labels
print(trainset.targets[:10])

Alright. Now, let's use PyTorch functionality for looping over the training and the test datasets. We need a DataLoader:

# One for the training data:
trainloader = torch.utils.data.DataLoader(trainset, batch_size=4,
                                          shuffle=True, num_workers=0)

# One for the test data:
testloader = torch.utils.data.DataLoader(testset, batch_size=4,
                                         shuffle=False, num_workers=0)

These objects work as follows:

# They help you loop over all the data in a random way (because we had shuffle=True)
for i, data in enumerate(trainloader, 0):
    inputs, labels = data
    # Here inputs are of size batch_size x (3 x 32 x 32)
    # Since we had specified, the batch_size to be 4
    # this essentially loads four images per iteration
    if i % 1000 == 0:
        print('Data point:', i, 'input size:', str(inputs.shape))

When you reach the end of the loop, you have visited all the images once. Notice that PyTorch has reshaped the images to 3 x 32 x 32 3D arrays. This is more convenient for the convolutional layers we will use later. Also, PyTorch uses the transformations we gave it to scale the data to array elements to $[-1, 1]$. Let me show you an example:

for i, data in enumerate(trainloader, 0):
    inputs, labels = data
    print(inputs[0])
    break

Training a classifier using dense DNNs

Let's train a classifier using a dense neural network. It could be better, but it is very easy to put together. We will start the network with 3 x 32 x 32 = 3072, followed by a few dense layers ending with ten outputs passed through softmax. However, for numerical stability reasons, we will not end with the softmax layer during training.

import torch.nn as nn

# The classifer - The dimensions of the layers have
# been picked to match those of the convolutional neural network
# that we are going to build later
# For now, just notice that we gradually take the 3072-dimensional input
# down to 10 dimensions (the number of classes we have)
# Also, notice that I do not add the softmax layer at this point
model_dense = nn.Sequential(nn.Linear(3072, 1176), nn.ReLU(),
                            nn.Linear(1176, 400), nn.ReLU(),
                            nn.Linear(400, 120), nn.ReLU(),
                            nn.Linear(120, 84), nn.ReLU(),
                            nn.Linear(84, 10))

# This is our loss function. 
# Read this: https://pytorch.org/docs/stable/generated/torch.nn.CrossEntropyLoss.html
criterion = nn.CrossEntropyLoss()
# The reason we did not add the Softmax layer at the end is because
# the loss function above is doing it internally.
# It expects that you provide "contain raw, unnormalized scores for each class"

# Here is the optimizer
import torch.optim as optim
optimizer = optim.SGD(model_dense.parameters(), lr=0.001, momentum=0.9)

Let's train the network. This is going to take a while (about 3 minutes on my desktop):

# How many times do you want to go over the entire dataset?
# Don't pick a very big number because you will overfit
num_epochs = 2

# Here is the main training algorithm
for epoch in range(num_epochs):  # loop over the dataset multiple times

    running_loss = 0.0
    for i, data in enumerate(trainloader, 0):
        # get the inputs; data is a list of [inputs, labels]
        inputs, labels = data

        # zero the parameter gradients
        optimizer.zero_grad()

        # forward + backward + optimize
        outputs = model_dense(inputs.reshape(4, 3 * 32 * 32))
        loss = criterion(outputs, labels)
        loss.backward()
        optimizer.step()

        # print statistics
        running_loss += loss.item()
        if i % 1000 == 999:    # print every 2000 mini-batches
            print('[%d, %5d] loss: %.3f' %
                  (epoch + 1, i + 1, running_loss / 1000))
            running_loss = 0.0

print('Finished Training')

Since training networks takes a while, it's a good idea to save it:

torch.save(model_dense.state_dict(), 'hands-on-25-model-dense.pth')

Here it is as a file:

!ls -lht hands-on-25-model-dense.pth

Now let's make some predictions:

# Get the first four images and their labels
dataiter = iter(testloader)
images, labels = next(dataiter)

print(labels)

# Make predictions with the net and pass them through 
# softmax to turn them into probabilities
st = nn.Softmax(dim=1)
predictions = st(model_dense(images.reshape(4, 3072)))

def imshow(img, ax):
    img = img / 2 + 0.5     # unnormalize
    npimg = img.numpy()
    ax.imshow(np.transpose(npimg, (1, 2, 0)))

# Plot the pictures and the predictions
for i in range(4):
    fig, ax = plt.subplots(figsize=(1,1))
    imshow(images[i], ax)
    fig2, ax2 = plt.subplots()
    ax2.bar(np.arange(10), predictions[i].detach().numpy())
    ax2.set_xticks(np.arange(10))
    ax2.set_xticklabels(classes)
    sns.despine(trim=True)

Now, let's do the same thing with a convolutional neural network. We are not going to use nn.Sequential this time. Instead, we will use nn.Module to manually create the network. The documentation is here. You need to inherit nn.Module, and implement __init__() and forward().

import torch.nn.functional as F

class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        # A convolutional layer:
        # 3 = input channels (colors),
        # 6 = output channels (features),
        # 5 = kernel size
        self.conv1 = nn.Conv2d(3, 6, 5)
        # A 2 x 2 max pooling layer - we are going to use it two times
        self.pool = nn.MaxPool2d(2, 2)
        # Another convolutional layer
        self.conv2 = nn.Conv2d(6, 16, 5)
        # Some linear layers
        self.fc1 = nn.Linear(16 * 5 * 5, 120)
        self.fc2 = nn.Linear(120, 84)
        self.fc3 = nn.Linear(84, 10)

    def forward(self, x):
        # This function implements your network output
        # Convolutional layer, followed by relu, followed by max pooling
        x = self.pool(F.relu(self.conv1(x)))
        # Same thing
        x = self.pool(F.relu(self.conv2(x)))
        # Flatting the output of the convolutional layers
        x = x.view(-1, 16 * 5 * 5)
        # Go throught the first dense linear layer followed by relu
        x = F.relu(self.fc1(x))
        # Through the second dense layer
        x = F.relu(self.fc2(x))
        # Finish up with a linear transformation
        x = self.fc3(x)
        return x


model_cnn = Net()

Here is the optimizer:

optimizer = optim.SGD(model_cnn.parameters(), lr=0.001, momentum=0.9)

And here is the training loop. Notice that this network trains faster.

# How many times do you want to go over the entire dataset?
# Don't pick a very big number because you will overfit
num_epochs = 5

# Here is the main training algorithm
for epoch in range(num_epochs):  # loop over the dataset multiple times

    running_loss = 0.0
    for i, data in enumerate(trainloader, 0):
        # get the inputs; data is a list of [inputs, labels]
        inputs, labels = data

        # zero the parameter gradients
        optimizer.zero_grad()

        # forward + backward + optimize
        outputs = model_cnn(inputs)
        loss = criterion(outputs, labels)
        loss.backward()
        optimizer.step()

        # print statistics
        running_loss += loss.item()
        if i % 1000 == 999:    # print every 2000 mini-batches
            print('[%d, %5d] loss: %.3f' %
                  (epoch + 1, i + 1, running_loss / 1000))
            running_loss = 0.0

print('Finished Training')

Make some predictions:

# Make predictions with the net and pass them through 
# softmax to turn them into probabilities
st = nn.Softmax(dim=1)
predictions = st(model_cnn(images))
for i in range(4):
    fig, ax = plt.subplots(figsize=(1,1))
    imshow(images[i], ax)
    fig2, ax2 = plt.subplots()
    ax2.bar(np.arange(10), predictions[i].detach().numpy())
    ax2.set_xticks(np.arange(10))
    ax2.set_xticklabels(classes)
    sns.despine(trim=True)

It doesn't work equally well for all classes. Here is some code from the PyTorch tutorial to get the accuracy for each class:

class_correct = list(0. for i in range(10))
class_total = list(0. for i in range(10))
with torch.no_grad():
    for data in testloader:
        images, labels = data
        outputs = model_cnn(images)
        _, predicted = torch.max(outputs, 1)
        c = (predicted == labels).squeeze()
        for i in range(4):
            label = labels[i]
            class_correct[label] += c[i].item()
            class_total[label] += 1


for i in range(10):
    print('Accuracy of %5s : %2d %%' % (
        classes[i], 100 * class_correct[i] / class_total[i]))

This could be better. There are several things that we can do. First, we would run this for more epochs. At least 50 epochs are needed to train it properly. Second, we could add data augmentation. This can be done through a transformation; see this. Third, we have to make the network bigger. Here is a list of large networks trained on CIFAR10. It is possible to reach an accuracy of 95%.

A general advice: When you train neural networks think big. Make them big. Add a lot of data. Train more (but not too much).

Questions

• Set the number of epochs for the CNN-based model to 40. How much better accuracy do you get? Make sure you do this before bed and look at it in the morning. It takes a while to train.