Variational Autoencoder in PyTorch, commented and annotated.

I have recently become fascinated with (Variational) Autoencoders and with PyTorch.

Kevin Frans has a beautiful blog post online explaining variational autoencoders, with examples in TensorFlow and, importantly, with cat pictures. Jaan Altosaar’s blog post takes an even deeper look at VAEs from both the deep learning perspective and the perspective of graphical models. Both of these posts, as well as Diederik Kingma’s original 2014 paper Auto-Encoding Variational Bayes, are more than worth your time.

In my case, I wanted to understand VAEs from the perspective of a PyTorch implementation. I started with the VAE example on the PyTorch github, adding explanatory comments and Python type annotations as I was working my way through it.

This post summarises my understanding, and contains my commented and annotated version of the PyTorch VAE example. I hope it helps!

What is PyTorch?

PyTorch is FAIR’s (that’s Facebook AI Research) Python dynamic deep learning / neural network library. The way that FAIR has managed to make neural network experimentation so dynamic and so natural is nothing short of miraculous. Read this post by fast.ai to find out more about their reasons for excitement, many of which I share.

What is an autoencoder?

The general idea of the autoencoder (AE) is to squeeze information through a narrow bottleneck between the mirrored encoder (input) and decoder (output) parts of a neural network. (see the diagram below)

Because the network achitecture and loss function are setup so that the output tries to emulate the input, the network has to learn how to encode input data on the very limited space represented by the bottleneck.

What is a variational autoencoder?

Variational Autoencoders, or VAEs, are an extension of AEs that additionally force the network to ensure that samples are normally distributed over the space represented by the bottleneck.

They do this by having the encoder output two n-dimensional (where n is the number of dimensions in the latent space) vectors representing the mean and the standard devation. These Gaussians are sampled, and the samples are sent through the decoder. This is the reparameterization step, also see my comments in the reparameterize() function.

What a fabulous trick!

The loss function has a term for input-output similarity, and, importantly, it has a second term that uses the KullbackÔÇôLeibler divergence to test how close the learned Gaussians are to unit Gaussians.

In other words, this extension to AEs enables us to derive Gaussian distributed latent spaces from arbitrary data. Given for example a large set of shapes, the latest space would be a high-dimensional space where each shape is represented by a single point, and the points would be normally distributed over all dimensions. With this one can represent existing shapes, but one can also synthesise completely new and plausible shapes by sampling points in latent space.

Results using MNIST

Below you see 64 random samples of a two-dimensional latent space of MNIST digits that I made with the example below, with ZDIMS=2.

pytorch-vae-sample-z2-epoch10.png

Next is the reconstruction of 8 random unseen test digits via a more reasonable 20-dimensional latent space. Keep in mind that the VAE has learned a 20-dimensional normal distribution for any input digit, from which samples are drawn that reconstruct via the decoder to output that appear similar to the input.

pytorch-vae-reconstruction-z10-epoch10.png

A diagram of a simple VAE

An example VAE, incidentally also the one implemented in the PyTorch code below, looks like this:

pytorch-vae-arch-2.png

A simple VAE implemented using PyTorch

I used PyCharm in remote interpreter mode, with the interpreter running on a machine with a CUDA-capable GPU to explore the code below. PyCharm parses the type annotations, which helps with code completion. I also made extensive use of the debugger to better understand logic flow and variable contents. (Debuggability is one of PyTorch’s strong points.)

Let me know in the comments to this post if you have any suggestions on how the code comments could be further improved.

# example from https://github.com/pytorch/examples/blob/master/vae/main.py
# commented and type annotated by Charl Botha <cpbotha@vxlabs.com>

import os
import torch
import torch.utils.data
from torch import nn, optim
from torch.autograd import Variable
from torch.nn import functional as F
from torchvision import datasets, transforms
from torchvision.utils import save_image

# changed configuration to this instead of argparse for easier interaction
CUDA = True
SEED = 1
BATCH_SIZE = 128
LOG_INTERVAL = 10
EPOCHS = 10

# connections through the autoencoder bottleneck
# in the pytorch VAE example, this is 20
ZDIMS = 20

# I do this so that the MNIST dataset is downloaded where I want it
os.chdir("/home/cpbotha/Downloads/pytorch-vae")

torch.manual_seed(SEED)
if CUDA:
    torch.cuda.manual_seed(SEED)

# DataLoader instances will load tensors directly into GPU memory
kwargs = {'num_workers': 1, 'pin_memory': True} if CUDA else {}

# Download or load downloaded MNIST dataset
# shuffle data at every epoch
train_loader = torch.utils.data.DataLoader(
    datasets.MNIST('data', train=True, download=True,
                   transform=transforms.ToTensor()),
    batch_size=BATCH_SIZE, shuffle=True, **kwargs)

# Same for test data
test_loader = torch.utils.data.DataLoader(
    datasets.MNIST('data', train=False, transform=transforms.ToTensor()),
    batch_size=BATCH_SIZE, shuffle=True, **kwargs)


class VAE(nn.Module):
    def __init__(self):
        super(VAE, self).__init__()

        # ENCODER
        # 28 x 28 pixels = 784 input pixels, 400 outputs
        self.fc1 = nn.Linear(784, 400)
        # rectified linear unit layer from 400 to 400
        # max(0, x)
        self.relu = nn.ReLU()
        self.fc21 = nn.Linear(400, ZDIMS)  # mu layer
        self.fc22 = nn.Linear(400, ZDIMS)  # logvariance layer
        # this last layer bottlenecks through ZDIMS connections

        # DECODER
        # from bottleneck to hidden 400
        self.fc3 = nn.Linear(ZDIMS, 400)
        # from hidden 400 to 784 outputs
        self.fc4 = nn.Linear(400, 784)
        self.sigmoid = nn.Sigmoid()

    def encode(self, x: Variable) -> (Variable, Variable):
        """Input vector x -> fully connected 1 -> ReLU -> (fully connected
        21, fully connected 22)

        Parameters
        ----------
        x : [128, 784] matrix; 128 digits of 28x28 pixels each

        Returns
        -------

        (mu, logvar) : ZDIMS mean units one for each latent dimension, ZDIMS
            variance units one for each latent dimension

        """

        # h1 is [128, 400]
        h1 = self.relu(self.fc1(x))  # type: Variable
        return self.fc21(h1), self.fc22(h1)

    def reparameterize(self, mu: Variable, logvar: Variable) -> Variable:
        """THE REPARAMETERIZATION IDEA:

        For each training sample (we get 128 batched at a time)

        - take the current learned mu, stddev for each of the ZDIMS
          dimensions and draw a random sample from that distribution
        - the whole network is trained so that these randomly drawn
          samples decode to output that looks like the input
        - which will mean that the std, mu will be learned
          *distributions* that correctly encode the inputs
        - due to the additional KLD term (see loss_function() below)
          the distribution will tend to unit Gaussians

        Parameters
        ----------
        mu : [128, ZDIMS] mean matrix
        logvar : [128, ZDIMS] variance matrix

        Returns
        -------

        During training random sample from the learned ZDIMS-dimensional
        normal distribution; during inference its mean.

        """

        if self.training:
            # multiply log variance with 0.5, then in-place exponent
            # yielding the standard deviation
            std = logvar.mul(0.5).exp_()  # type: Variable
            # - std.data is the [128,ZDIMS] tensor that is wrapped by std
            # - so eps is [128,ZDIMS] with all elements drawn from a mean 0
            #   and stddev 1 normal distribution that is 128 samples
            #   of random ZDIMS-float vectors
            eps = Variable(std.data.new(std.size()).normal_())
            # - sample from a normal distribution with standard
            #   deviation = std and mean = mu by multiplying mean 0
            #   stddev 1 sample with desired std and mu, see
            #   https://stats.stackexchange.com/a/16338
            # - so we have 128 sets (the batch) of random ZDIMS-float
            #   vectors sampled from normal distribution with learned
            #   std and mu for the current input
            return eps.mul(std).add_(mu)

        else:
            # During inference, we simply spit out the mean of the
            # learned distribution for the current input.  We could
            # use a random sample from the distribution, but mu of
            # course has the highest probability.
            return mu

    def decode(self, z: Variable) -> Variable:
        h3 = self.relu(self.fc3(z))
        return self.sigmoid(self.fc4(h3))

    def forward(self, x: Variable) -> (Variable, Variable, Variable):
        mu, logvar = self.encode(x.view(-1, 784))
        z = self.reparameterize(mu, logvar)
        return self.decode(z), mu, logvar


model = VAE()
if CUDA:
    model.cuda()


def loss_function(recon_x, x, mu, logvar) -> Variable:
    # how well do input x and output recon_x agree?
    BCE = F.binary_cross_entropy(recon_x, x.view(-1, 784))

    # KLD is Kullback–Leibler divergence -- how much does one learned
    # distribution deviate from another, in this specific case the
    # learned distribution from the unit Gaussian

    # see Appendix B from VAE paper:
    # Kingma and Welling. Auto-Encoding Variational Bayes. ICLR, 2014
    # https://arxiv.org/abs/1312.6114
    # - D_{KL} = 0.5 * sum(1 + log(sigma^2) - mu^2 - sigma^2)
    # note the negative D_{KL} in appendix B of the paper
    KLD = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
    # Normalise by same number of elements as in reconstruction
    KLD /= BATCH_SIZE * 784

    # BCE tries to make our reconstruction as accurate as possible
    # KLD tries to push the distributions as close as possible to unit Gaussian
    return BCE + KLD

# Dr Diederik Kingma: as if VAEs weren't enough, he also gave us Adam!
optimizer = optim.Adam(model.parameters(), lr=1e-3)


def train(epoch):
    # toggle model to train mode
    model.train()
    train_loss = 0
    # in the case of MNIST, len(train_loader.dataset) is 60000
    # each `data` is of BATCH_SIZE samples and has shape [128, 1, 28, 28]
    for batch_idx, (data, _) in enumerate(train_loader):
        data = Variable(data)
        if CUDA:
            data = data.cuda()
        optimizer.zero_grad()

        # push whole batch of data through VAE.forward() to get recon_loss
        recon_batch, mu, logvar = model(data)
        # calculate scalar loss
        loss = loss_function(recon_batch, data, mu, logvar)
        # calculate the gradient of the loss w.r.t. the graph leaves
        # i.e. input variables -- by the power of pytorch!
        loss.backward()
        train_loss += loss.data[0]
        optimizer.step()
        if batch_idx % LOG_INTERVAL == 0:
            print('Train Epoch: {} [{}/{} ({:.0f}%)]\tLoss: {:.6f}'.format(
                epoch, batch_idx * len(data), len(train_loader.dataset),
                100. * batch_idx / len(train_loader),
                loss.data[0] / len(data)))

    print('====> Epoch: {} Average loss: {:.4f}'.format(
          epoch, train_loss / len(train_loader.dataset)))


def test(epoch):
    # toggle model to test / inference mode
    model.eval()
    test_loss = 0

    # each data is of BATCH_SIZE (default 128) samples
    for i, (data, _) in enumerate(test_loader):
        if CUDA:
            # make sure this lives on the GPU
            data = data.cuda()

        # we're only going to infer, so no autograd at all required: volatile=True
        data = Variable(data, volatile=True)
        recon_batch, mu, logvar = model(data)
        test_loss += loss_function(recon_batch, data, mu, logvar).data[0]
        if i == 0:
          n = min(data.size(0), 8)
          # for the first 128 batch of the epoch, show the first 8 input digits
          # with right below them the reconstructed output digits
          comparison = torch.cat([data[:n],
                                  recon_batch.view(BATCH_SIZE, 1, 28, 28)[:n]])
          save_image(comparison.data.cpu(),
                     'results/reconstruction_' + str(epoch) + '.png', nrow=n)

    test_loss /= len(test_loader.dataset)
    print('====> Test set loss: {:.4f}'.format(test_loss))


for epoch in range(1, EPOCHS + 1):
    train(epoch)
    test(epoch)

    # 64 sets of random ZDIMS-float vectors, i.e. 64 locations / MNIST
    # digits in latent space
    sample = Variable(torch.randn(64, ZDIMS))
    if CUDA:
        sample = sample.cuda()
    sample = model.decode(sample).cpu()

    # save out as an 8x8 matrix of MNIST digits
    # this will give you a visual idea of how well latent space can generate things
    # that look like digits
    save_image(sample.data.view(64, 1, 28, 28),
               'results/sample_' + str(epoch) + '.png')

Miniconda3, TensorFlow, Keras on Google Compute Engine GPU instance: The step-by-step guide.

Google recently announced the availability of GPUs on Google Compute Engine instances. For my deep learning experiments, I often need more beefy GPUs than the puny GTX 750Ti in my desktop workstation, so this was good news. To make the GCE offering even more attractive, their GPU instances are also available in their EU datacenters, which is in terms of latency a big plus for me here on the Southern tip of the African continent.

Last night I had some time to try this out, and in this post I would like to share with you all the steps I took to:

  1. Get a GCE instance with GPU up and running with miniconda, TensorFlow and Keras
  2. Create a reusable disk image with all software pre-installed so that I could bring up new instances ready-to-roll at the drop of a hat.
  3. Apply the pre-trained Resnet50 deep neural network on images from the web, as a demonstration that the above works. Thanks to Keras, this step is fun and fantastically straight-forward.

Pre-requisites

I started by creating a project for this work. On the Compute Engine console, check that this project is active at the top.

Before I was able to allocate GPUs to my instance, I had to fill in the “request quote increase” form available from the Compute Engine quotas page. My request for two GPUs in the EU region was approved within minutes.

I installed my client workstation’s id_rsa.pub public SSH key as a project-wide SSH key via the metadata screen.

Start an instance for the first time

I configured my GPU instance as shown in the following screenshot:

gce-create-instance.png

  • Under Machine type switch to Customize to be able to select a GPU.
  • I selected an Ubuntu 16.04 image, and changed the persistent disk to SSD.
  • I selected the europe-west1-b zone. Choose whatever is closest for you. The interface will warn you if the selection does NOT support GPUs.

After this, click on the Create button and wait for your instance to become ready.

Once it’s up and running, you’ll be able to ssh to the displayed public IP. I used the ssh on my client workstation, but of course you could opt for the Google-supplied web-based versions.

Install NVIDIA drivers and CUDA

I used the following handy script from the relevant GCE documentation:

#!/bin/bash
echo "Checking for CUDA and installing."
# Check for CUDA and try to install.
if ! dpkg-query -W cuda; then
  # The 16.04 installer works with 16.10.
  curl -O http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1604/x86_64/cuda-repo-ubuntu1604_8.0.61-1_amd64.deb
  dpkg -i ./cuda-repo-ubuntu1604_8.0.61-1_amd64.deb
  apt-get update
  apt-get install cuda -y
fi

After this, download the CUDNN debs from the NVIDIA download site using your developer account. Install the two debs using dpkg -i.

To confirm that the drivers have been installed, run the nvidia-smi command:

gce-nvidia-smi.png

Install miniconda, tensorflow and keras

I usually download the 64bit Linux miniconda installer from conda.io and then install it into ~/miniconda3 by running the downloaded .sh script.

After this, I installed TensorFlow 1.0.1 and Keras 2.0.1 into a new conda environment by doing:

conda create -n ml python=3.6
conda install jupyter pandas numpy scipy scikit-image
pip install --ignore-installed --upgrade https://storage.googleapis.com/tensorflow/linux/gpu/tensorflow_gpu-1.0.1-cp36-cp36m-linux_x86_64.whl
pip install keras h5py

The keras package also installed theano, which I then uninstalled using pip uninstall theano in the active ml environment.

To test, run ipython and then type import keras. It should look like this:

gce-import-keras.png

Note that it’s picking up the TensorFlow backend, and successfully loading all of the CUDA librarias, including CUDNN.

Save your disk as an image for later

You will get billed for each minute that the instance is running. You also get billed for persistent disks that are still around, even if they are not used by any instance.

Creating a reusable disk image will enable you to delete instances and disks, and later to restart an instance with all of your software already installed.

To do this, follow the steps in the documentation, which I paraphrase and extend here:

  1. Stop the instance.
  2. In the instance list, click on the instance name itself; this will take you to the edit screen.
  3. Click the edit button, and then uncheck Delete boot disk when instance is deleted.
  4. Click the save button.
  5. Delete the instance, but double-check that delete boot disk is unchecked in the confirmation dialog.
  6. Now go to the Images screen and select Create Image with the boot disk as source.

Next time, go to the Images screen, select your image and then select Create Instance. That instance will come with all of your goodies ready to go!

Apply the ResNet50 neural network on images from the interwebs

After connecting to the instance with an SSH port redirect:

ssh -L 8889:localhost:8888 cpbotha@EXTERNAL_IP

… and then starting a jupyter notebook on the GCE instance:

cpbotha@instance-1:~$ source ~/miniconda3/bin/activate ml
(ml) cpbotha@instance-1:~$ jupyter notebook

So that I can connect to the notebook on my localhost:8889, I enter and execute the following code (adapted from the Keras documentation) in a cell:

from keras.applications.resnet50 import ResNet50
from keras.preprocessing import image
from keras.applications.resnet50 import preprocess_input, decode_predictions
import numpy as np
from PIL import Image

model = ResNet50(weights='imagenet')

# adapted from https://github.com/fchollet/deep-learning-models
# to accept also a PIL image
def load_and_predict_image(img_or_path):
    target_size = (224,224)
    if type(img_or_path) is str:
        img = image.load_img(img_or_path, target_size=target_size)

    else:
        img = img_or_path.resize(target_size)

    x = image.img_to_array(img)
    x = np.expand_dims(x, axis=0)
    x = preprocess_input(x)

    preds = model.predict(x)
    # decode the results into a list of tuples (class, description, probability)
    # (one such list for each sample in the batch)
    print('Predicted:', decode_predictions(preds, top=3)[0])

In the next cell, I do:

from PIL import Image
import urllib.request

url1 = "https://upload.wikimedia.org/wikipedia/commons/thumb/c/c8/KuduKr%C3%BCger.jpg/1920px-KuduKr%C3%BCger.jpg"
url2 = "https://upload.wikimedia.org/wikipedia/commons/thumb/9/9d/Struthio_camelus_-_Etosha_2014_%283%29.jpg/800px-Struthio_camelus_-_Etosha_2014_%283%29.jpg"
im = Image.open(urllib.request.urlopen(url2))

load_and_predict_image(im)

To be greeted with the following results:

gce-keras-resnet.png

The pre-trained ResNet50 network identified the Kudu I gave it initially as an ostrich, so I decided to make it a bit easier for the poor network by actually giving it an ostrich, which it did identify with a 99.98% probability.

Looking at the photos, the former taken in the Kruger National Park and the second in Etosha, I can image that the network could identify the former as the latter due to similar background and foreground colouring, and clearly having not been trained on Kudu.

Let’s file that under future work!