Let me preface this blog post by saying that the topic is somewhat complicated. Even a moderate explanation would take many pages, and so my summary will leave out a lot of detail.

A new technique for anomaly detection is to use a variational autoencoder and compute a metric called reconstruction probability. The idea is related to, but significantly different from anomaly detection using a standard autoencoder, and is also significantly differently from using a VAE with reconstruction error.

At a high level, a VAE computes an internal representation of the probability distribution of the source data in the form of a mean and a standard deviation (actually the log of the variance — one of many tricky details). A VAE outputs the distribution mean and log-variance, and a reconstructed version of the input. The idea of reconstruction probability anomaly detection is to compute a second probability distribution and then use it to calculate the likelihood that an input item came from the distribution. Data items with a low reconstruction probability are not likely to have come from the distribution, and so are anomalous in some way.

The source research paper is “Variational Autoencoder Based Anomaly Detection Using Reconstruction Probability” (2015), written by Jinwon An, and Sungzoon Cho. The paper is simultaneously excellent and . . . not so excellent. The idea of reconstruction probability is very clever, and the background and motivation are clearly explained. But the exact definition of “reconstruction probability” is never given, and the hints about how to implement it are somewhat contradictory.

I spent several days exploring the ideas of VAE reconstruction probability and finally got a possible demo running. I say “possible demo” because the ideas are so complex, there are dozens of places where I might have gone wrong.

The algorithm in the research paper is shown in the image below:

A close inspection of the algorithm reveals that it has ambiguities. For example, the g() function is the decoder component of a VAE, but g() in a standard VAE architecture does not return a mean and standard deviation as indicated by the algorithm. This could mean the researchers used a modified VAE architecture (where the output is a second mean and standard deviation — probably, or it could mean something else. The bottom line is that VAE reconstruction probability is not rigidly defined.

For this experiment, I set up a VAE with the architecture shown in the diagram below. The source X determines a u1 mean and s1 standard deviation which define a Normal distribution. That distribution determines a reconstructed x. To compute a reconstruction probability, I use the reconstructed X as a stand-in for the mean of a second Normal distribution. I use a dummy standard deviation of all 1s for the second distribution.

My preliminary PyTorch implementation of reconstruction probability is:

import torch as T
import numpy as np
import scipy.stats as sps

def recon_prob(model, xi, n_samples=10):
  # xi is one Tensor item, not a batch
  # assumes model.eval() has been set
  
  with T.no_grad():
    (u, logvar) = model.encode(xi)
  u = u.numpy()               # mean
  v = np.exp(logvar.numpy())  # variance

  samples = sps.multivariate_normal.rvs(u, \
    np.diag(v), size=n_samples) 
  samples = T.tensor(samples, \
    dtype=T.float32).to(device)  

  with T.no_grad():
    x_computeds = model.decode(samples)
  est_means = T.mean(x_computeds, dim=0)
  est_vars = T.ones(len(est_means), \
    dtype=T.float32).to(device)
  est_prob = sps.multivariate_normal.pdf( \
    xi.numpy(), est_means.numpy(), \
    np.diag(est_vars.numpy()))

  return est_prob

The recon_prob() function accepts a trained VAE model and a single data item. The encoder computes a mean and a log-var. These are used to construct a multivariate Normal distribtuion, which in turn is used to generate a set of samples. These samples are fed to the decoder component of the VAE. The outputs are averaged to produce a mean. I make a simplifying assumption and use a diagonal matrix of 1-values to stand in for the covariance matrix. The mean and covariance matrix produce a new distribution, and then the distribution’s pdf() function (probability density function) emits the result probability. My implementation gives a single probability, unlike the research paper which indicates it produces L probabilities and then averages them.

I coded up a demo where I use 240 dummy Employee data items. Each data item has a sex, age, city, income, and job-type. My VAE reconstruction probability implementation appears to be working, but there are many hours of exploration ahead before I’ll be ready to say I have an implementation that’s ready for posting. I’m looking forward to traveling the journey of investigation.

Three 1960s United Airlines travel posters by artist Stan Galli (1912-2009). I like his style a lot. These posters emit a sense of joy with probability = 1.