{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Normalizing Flows\n",
"\n",
"The VAE was our first example of a generative model that is capable of sampling from $P(x)$. The VAE has two disadvantages though. The first that you cannot numerically evaluate $P(x)$, the probability of a single point. The second disadvantage is that training the VAE is difficult, especially because you're assuming the latent space should be normal. **Generative adversarial networks** (GANs) are similar to VAEs and have the same two disadvantages.\n",
"\n",
"A **normalizing flow** is similar to a VAE in that we try to build up $P(x)$ by starting from a simple known distribution $P(z)$. We use functions, like the decoder from a VAE, to go from $x$ to $z$. However, we make sure that the functions we choose keep the probability mass normalized ($\\sum P(x) = 1$) and can be used forward (to sample from x) and backward (to compute $P(x)$). We call these functions **bijectors** because they are bijective (1 to 1, onto). An example of a bijector is an element-wise cosine $y_i = \\cos x_i$ (assuming $x_i$ is between 0 and $\\pi$) and non-bijective function would be a reduction $y = \\sum_i x_i$. Any function which changes the number of elements is automatically not bijective. A consequence of using only bijectors in constructing our normalizing flow is that the size of the latent space must be equal to the size of the feature space. Remember the VAE used a smaller latent space than the feature space. \n",
"\n",
"\n",
"```{admonition} Audience & Objectives\n",
"This chapter builds on {doc}`VAE` and assumes the same background of probability theory. This chapter is an introduction to the key ideas, but is not fully developed yet. Some knowledge of vector calculus (Jacobians) is assumed as well. After completing it, you should be able to \n",
"\n",
" * Understand the trade-offs between a VAE, GAN, and normalizing flow. \n",
" * Identify a bijector and construct a bijector chain\n",
" * Construct a normalizing flow using common bijectors types and train it\n",
" * Sample from a normalizing flow and compute sample probabilities \n",
"```\n",
"\n",
"You can find a recent review of normalizing [flows here](https://arxiv.org/pdf/1908.09257.pdf) {cite}`kobyzev2020normalizing` and [here](https://arxiv.org/abs/1912.02762){cite}`papamakarios2019normalizing`. Although generating images and sound is the most popular application of normalizing flows, some of their biggest scientific impact has been on more efficient sampling from posteriors or likelihoods and other complex probability distributions {cite}`papamakarios2019sequential`. You find details on how to do normalizing flows on categorical (discrete) data in Hoogeboom et al. {cite}`hoogeboom2021argmax`.\n",
"\n",
"## Flow Equation\n",
"\n",
"Recall for the VAE decoder, we had an explicit formula for $p(x | z)$. This allowed us to compute $p(x) = \\int\\,dz p(x | z)p(z)$ which is the quantity of interest. The VAE decoder is a conditional probability density function. In the normalizing flow, we do not use probability density functions. We use bijective functions. So we cannot just compute an integral to change variables. We can use the change of variable formula. Consider our normalizing flow to be defined by our bijector $x = f(z)$, its inverse $z = g(x)$, and the starting probability distribution $P_z(z)$. Then the formula for probability of $x$ is\n",
"\n",
"\\begin{equation}\n",
"P(x) = P_z\\left(g(x)\\right) \\,\\left| \\textrm{det}\\left[\\mathbf{J}_g\\right]\\right|\n",
"\\end{equation}\n",
"\n",
"where the term on the right is the absolute value of the determinant of the Jacobian of $g$. [Jacobians](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) are matrices that describe how infinitesimal changes in each domain dimension change each range dimension. This term corrects for the volume change of the distribution. For example, if $f(z) = 2z$, then $g(x) = x / 2$, and the Jacobian determinant is $1 / 2$. The intuition is that we are stretching out $z$ by 2, so we need to account for the increase in volume to keep the probability normalized. You can read more about the change of variable formula for [probability distributions here](https://cranmer.github.io/stats-ds-book/distributions/change-of-variables.html)\n",
"\n",
"\n",
"## Bijectors\n",
"\n",
"A bijector is a function that is [injective](https://en.wikipedia.org/wiki/Injective_function) (1 to 1) and [surjective](https://en.wikipedia.org/wiki/Surjective_function) (onto). An equivalent way to view a bijective function is if it has an inverse. For example, a sum reduction has no inverse and is thus not bijective. $\\sum [1,0] = 1$ and $\\sum [-1, 2] = 1$. Multiplying by a matrix which has an inverse is bijective. $y = x^2$ is not bijective, since $y = 4$ has two solutions. \n",
"\n",
"Remember that we must compute the determinant of the bijector Jacobian. If the Jacobian is dense (all output elements depend on all input elements), computing this quantity will be $O\\left(|x|_0^3\\right)$ where $|x|_0$ is the number of dimensions of $x$ because a determinant scales by $O(n^3)$. This would make computing normalizing flows impractical in high-dimensions. However, in practice we restrict ourselves to bijectors that have easy to calculate Jacobians. For example, if the bijector is $x_i = \\cos z_i$ then the Jacobian will be diagonal. Typically, the trick that is done is to make the Jacobian triangular. Then $x_0$ only depends on $z_0$, $z_1$ depends on $z_0, Z_1$, and $x_2$ depends on $z_0, z_1, z_2$, etc. The matrix determinant is then computed in linear time with respect to the number of dimensions.\n",
"\n",
"### Bijector Chains\n",
"\n",
"Just like in deep neural networks, multiple bijectors are chained together to increase how complex of the final fit distribution $\\hat{P}(x)$ can be. The change of variable equation can be repeatedly applied:\n",
"\n",
"\\begin{equation}\n",
"P(x) = P_z\\left[g_1\\left(g_0(x)\\right)\\right] \\,\\left| \\textrm{det}\\left[\\mathbf{J}_{g_1}\\right]\\right| \\left|\\textrm{det}\\left[\\mathbf{J}_{g_0}\\right]\\right|\n",
"\\end{equation}\n",
"\n",
"where we would compute $x$ with $f_0\\left(f_1(z)\\right)$. One critical point is that you should also include a **permute bijector** that swaps the order of dimensions. Since the bijectors typically have triangular Jacobians, certain output dimensions will depend on many input dimensions and others will only depend on a single one. By applying a permutation, you allow each dimension to influence each other.\n",
"\n",
"\n",
"## Training\n",
"\n",
"At this point, you may be wondering how you could possibly train a normalizing flow. The trainable parameters appear in the bijectors. They have adjustable parameters. The loss equation is quite simple: the negative log-likelihood (negative to make it minimization). Explicitly:\n",
"\n",
"\\begin{equation}\n",
"\\mathcal{l} = -\\log P_z\\left[g_1\\left(g_0(x)\\right)\\right] - \\sum_i \\log\\left| \\textrm{det}\\left[\\mathbf{J}_{g_i}\\right]\\right|\n",
"\\end{equation}\n",
"\n",
"where $x$ is the training point and when you take the gradient of the loss, it is with respect to the parameters of the bijectors. \n",
"\n",
"\n",
"## Common Bijectors\n",
"\n",
"The choice of bijector functions is a fast changing area. I will thus only mention a few. You can of course use any bijective function or matrix, but these become inefficient at high-dimension due to the Jacobian calculation. One class of efficient bijectors are autoregressive bijectors. These have triangular Jacobians because each output dimension can only depend on the dimensions with a lower index. There are two variants: masked autoregressive flows (MAF){cite}`papamakarios2017masked` and inverse autoregressive flows (IAF) {cite}`kingma2016improved`. MAFs are efficient at training and computing probabilities, but are slow for sampling from $P(x)$. IAFs are slow at training and computing probabilities but efficient for sampling. Wavenets combine the advantages of both {cite}`kim2018flowavenet`. I'll mention one other common bijector which is not autoregressive: real non-volume preserving (RealNVPs) {cite}`dinh2016density`. RealNVPs are less expressive than IAFs/MAFs, meaning they have trouble replicating complex distributions, but are efficient at all three tasks: training, sampling, and computing probabilities. Another interesting variant is the Glow bijector,which is able to expand the rank of the normalizing flow, for example going from a matrix to an RGB image {cite}`das2019dimensionality`. What are the equations for all these bijectors? Most are variants of standard neural network layers but with special rules about which outputs depend on which inputs. \n",
"\n",
"\n",
"```{warning}\n",
"Remember to add permute bijectors between autoregressive bijectors to ensure the dependence between dimensions is well-mixed.\n",
"```\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Running This Notebook\n",
"\n",
"\n",
"Click the above to launch this page as an interactive Google Colab. See details below on installing packages.\n",
"\n",
"````{tip} My title\n",
":class: dropdown\n",
"To install packages, execute this code in a new cell. \n",
"\n",
"```\n",
"!pip install dmol-book\n",
"```\n",
"\n",
"If you find install problems, you can get the latest working versions of packages used in [this book here](https://github.com/whitead/dmol-book/blob/master/package/requirements.txt)\n",
"\n",
"````"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The hidden code below imports the tensorflow probability package and other necessary packages. Note that the tensorflow probability package (`tfp`) is further broken into distributions (`tfd`) and bijectors (`tfb`). "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import tensorflow as tf\n",
"import tensorflow_probability as tfp\n",
"import matplotlib.pyplot as plt\n",
"import sklearn.datasets as datasets\n",
"import numpy as np\n",
"import dmol"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(0)\n",
"tf.random.set_seed(0)\n",
"\n",
"\n",
"tfd = tfp.distributions\n",
"tfb = tfp.bijectors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Moon Example\n",
"\n",
"We'll start with a basic 2D example to learn the two moons distribution with a normalizing flow. Two moons is a common example dataset that is hard to cluster and model as a probability distribution. \n",
"\n",
"When doing normalizing flows you have two options to implement them. You can do all the Jacobians, inverses, and likelihood calculations analytically and implement them in a normal ML framework like Jax, PyTorch, or TensorFlow. This is actually most common. The second option is to utilize a probability library that knows how to use bijectors and distributions. The packages for that are PYMC, TensorFlow Probability (which has a non-tensorflow JAX version confusingly), and Pyro (Pytorch). We'll use TensorFlow Probability for this work.\n",
"\n",
"### Generating Data\n",
"\n",
"In the code below, I set-up my imports and sample points which will be used for training. Remember, this code has nothing to do with normalizing flows -- it's just to generate data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"moon_n = 10000\n",
"ndim = 2\n",
"data, _ = datasets.make_moons(moon_n, noise=0.05)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(data[:, 0], data[:, 1], \".\", alpha=0.8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Z Distribution\n",
"\n",
"Our Z distribution should always be as simple as possible. I'll create a 2D Gaussian with unit variance, no covariance, and centered at 0. I'll be using the tensorflow probability package for this example. The key new concept is that we organize our tensors that were *sampled* from a probability distribution in a specific way. We, by convention, make the first axes be the **sample** shape, the second axes be the **batch** shape, and the final axes be the **event** shape. The sample shape is the number of times we sampled from our distribution. It is a *shape* and not a single dimension because occasionally you'll want to organize your samples into some shape. The batch shape is a result of possibly multiple distributions batched together. For example, you might have 2 Gaussians, instead of a single 2D Gaussian. Finally, the event shape is the shape of a single sample from the distribution. This is overly complicated for our example, but you should be able to read information about the distribution now by understanding this nomenclature. You can find a tutorial on these [shapes here](https://www.tensorflow.org/probability/examples/Understanding_TensorFlow_Distributions_Shapes) and more tutorials on [tensorflow probability here](https://www.tensorflow.org/probability/examples/A_Tour_of_TensorFlow_Probability). "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"zdist = tfd.MultivariateNormalDiag(loc=[0.0] * ndim)\n",
"zdist"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With our new understanding of shapes, you can see that this distribution has no `batch_shape` because there is only one set of parameters and the `event_shape` is `[2]` because it's a 2D Gaussian. Let's now sample from this distribution and view it"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"zsamples = zdist.sample(moon_n)\n",
"plt.plot(zsamples[:, 0], zsamples[:, 1], \".\", alpha=0.8)\n",
"plt.xlim(-4, 4)\n",
"plt.ylim(-4, 4)\n",
"plt.gca().set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As expected, our starting distribution looks nothing like are target distribution. Let's demonstrate a bijector now. We'll implement the following bijector:\n",
"\n",
"$$\n",
"x = \\vec{z} \\times (1, 0.5)^T + (0.5, 0.25)\n",
"$$\n",
"\n",
"This is bijective because the operations are element-wise and invertible. Rather than just write this out using operations like `@` or `*`, we'll use the built-in bijectors from TensorFlow probability. The reason we do this is that they have their inverses and Jacobian determinants already defined. We first create a bijector that scales by `[1, 0.5]` and one then one that shifts by `[0.5,0.25]`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"shift_bij = tfb.Shift([0.5, 0.25])\n",
"scale_bij = tfb.Scale([1, 0.5])\n",
"# make composite via function convolution\n",
"b = shift_bij(scale_bij)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To now apply the change of variable formula, we create a **transformed distribution**. What is important about this choice is that we can compute likelihoods, probabilities, and sample from it."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"td = tfd.TransformedDistribution(zdist, bijector=b)\n",
"td"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"zsamples = td.sample(moon_n)\n",
"plt.plot(zsamples[:, 0], zsamples[:, 1], \".\", alpha=0.8)\n",
"plt.xlim(-4, 4)\n",
"plt.ylim(-4, 4)\n",
"plt.gca().set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We show above the sampling from this new distribution. We can also plot it's probability, which is impossible for a VAE-like model!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# make points for grid\n",
"zpoints = np.linspace(-4, 4, 150)\n",
"(\n",
" z1,\n",
" z2,\n",
") = np.meshgrid(zpoints, zpoints)\n",
"zgrid = np.concatenate((z1.reshape(-1, 1), z2.reshape(-1, 1)), axis=1)\n",
"# compute P(x)\n",
"p = np.exp(td.log_prob(zgrid))\n",
"fig = plt.figure()\n",
"# plot and set axes limits\n",
"plt.imshow(p.reshape(z1.shape), aspect=\"equal\", extent=[-4, 4, -4, 4])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Normalizing Flow\n",
"\n",
"Now we will build bijectors that are expressive enough to capture the moon distribution. I will use 3 sets of a MAF and permutation for 6 total bijectors. MAF's have dense neural network layers in them, so I will also set the usual parameters for a neural network: dimension of hidden layer and activation. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"num_layers = 3\n",
"my_bijects = []\n",
"# loop over desired bijectors and put into list\n",
"for i in range(num_layers):\n",
" # Syntax to make a MAF\n",
" anet = tfb.AutoregressiveNetwork(\n",
" params=ndim, hidden_units=[128, 128], activation=\"relu\"\n",
" )\n",
" ab = tfb.MaskedAutoregressiveFlow(anet)\n",
" # Add bijector to list\n",
" my_bijects.append(ab)\n",
" # Now permuate (!important!)\n",
" permute = tfb.Permute([1, 0])\n",
" my_bijects.append(permute)\n",
"# put all bijectors into one \"chain bijector\"\n",
"# that looks like one\n",
"big_bijector = tfb.Chain(my_bijects)\n",
"# make transformed dist\n",
"td = tfd.TransformedDistribution(zdist, bijector=big_bijector)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this point, we have not actually trained but we can still view our distribution. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"zpoints = np.linspace(-4, 4, 150)\n",
"(\n",
" z1,\n",
" z2,\n",
") = np.meshgrid(zpoints, zpoints)\n",
"zgrid = np.concatenate((z1.reshape(-1, 1), z2.reshape(-1, 1)), axis=1)\n",
"p = np.exp(td.log_prob(zgrid))\n",
"fig = plt.figure()\n",
"plt.imshow(p.reshape(z1.shape), aspect=\"equal\", extent=[-4, 4, -4, 4])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can already see that the distribution looks more complex than a Gaussian. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Training\n",
"To train, we'll use TensorFlow Keras, which just handles computing derivatives and the optimizer. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# declare the feature dimension\n",
"x = tf.keras.Input(shape=(2,), dtype=tf.float32)\n",
"# create a \"placeholder\" function\n",
"# that will be model output\n",
"log_prob = td.log_prob(x)\n",
"# use input (feature) and output (log prob)\n",
"# to make model\n",
"model = tf.keras.Model(x, log_prob)\n",
"# define a loss\n",
"def neg_loglik(yhat, log_prob):\n",
" # losses always take in label, prediction\n",
" # in keras. We do not have labels,\n",
" # but we still need to accept the arg\n",
" # to comply with Keras format\n",
" return -log_prob\n",
"\n",
"\n",
"# now we prepare model for training\n",
"model.compile(optimizer=tf.optimizers.Adam(1e-3), loss=neg_loglik)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One detail is that we have to create fake labels (zeros) because Keras expects there to always be training labels. Thus our loss we defined above (negative log-likelihood) takes in the labels but does nothing with them. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"result = model.fit(x=data, y=np.zeros(moon_n), epochs=100, verbose=0)\n",
"plt.plot(result.history[\"loss\"])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Training looks reasonable. Let's now see our distribution. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"zpoints = np.linspace(-2.5, 2.5, 200)\n",
"(\n",
" z1,\n",
" z2,\n",
") = np.meshgrid(zpoints, zpoints)\n",
"zgrid = np.concatenate((z1.reshape(-1, 1), z2.reshape(-1, 1)), axis=1)\n",
"p = np.exp(td.log_prob(zgrid))\n",
"fig = plt.figure()\n",
"plt.imshow(\n",
" p.reshape(z1.shape), aspect=\"equal\", origin=\"lower\", extent=[-2.5, 2.5, -2.5, 2.5]\n",
")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wow! We now can compute the probability of any point in this distribution. You can see there are some oddities that could be fixed with further training. One issue that cannot be overcome is the connection between the two curves -- it is not possible to get fully disconnected densities. This is because of our requirement that the bijectors are invertible and volume preserving -- you can only squeeze volume so far but cannot completely disconnect. Some work has been done on addressing this issue by adding sampling to the flow and this gives more expressive normalizing flows {cite}`wu2020stochastic`.\n",
"\n",
"Finally, we'll sample from our model just to show that indeed it is generative."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"zsamples = td.sample(moon_n)\n",
"plt.plot(zsamples[:, 0], zsamples[:, 1], \".\", alpha=0.2, markeredgewidth=0.0)\n",
"plt.xlim(-2.5, 2.5)\n",
"plt.ylim(-2.5, 2.5)\n",
"plt.gca().set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Relevant Videos\n",
"\n",
"### Normalizing Flow for Molecular Conformation\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Chapter Summary \n",
"\n",
"* A normalizing flow builds up a probability distribution of $x$ by starting from a known distribution on $z$. Bijective functions are used to go from $z$ to $x$.\n",
"* Bijectors are functions that keep the probability mass normalized and are used to go forward and backward (because they have well-defined inverses).\n",
"* To find the probability distribution of $x$ we use the change of variable formula, which requires a function inverse and Jacobian.\n",
"* The bijector function has trainable parameters, which can be trained using a negative log-likelihood function. \n",
"* Multiple bijectors can be chained together, but typically must include a permute bijector to swap the order of dimensions. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cited References\n",
"\n",
"```{bibliography}\n",
":style: unsrtalpha\n",
":filter: docname in docnames\n",
"```"
]
}
],
"metadata": {
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}