Normalizing Flow

A normalizing flow is a generative model built from an invertible neural network. The network maps an image $x$ to a code $z$ of the same size, and a simple distribution, usually a standard Gaussian, is placed on $z$. Because the map can be inverted, the same network generates images, by drawing $z$ and running it backwards, and scores them: the probability density of any image follows exactly from the density of its code and the change of volume the network applies, the change-of-variables formula. Flows are therefore trained by plain maximum likelihood, with no bound, as in the variational autoencoder, and no adversary, as in the generative adversarial network. They are also called flow-based generative models, invertible generative models or reversible generative models.

The price is architectural. Every layer must be invertible, keep the dimension of its input, and have a Jacobian determinant that is cheap to compute. The layers that meet these constraints, coupling layers, autoregressive layers, invertible $1 \times 1$ convolutions and continuous-time flows, are the subject of this article, together with the way flows are fitted to discrete pixel values (dequantization) and scored (bits per dimension). The generative model article compares flows with the other families.

Core Idea

Write $f$ for the map from an image to its code, $z = f(x)$, and $g = f^{-1}$ for the generator, $x = g(z)$. Probability mass is conserved: a small region around $z$ and its image around $x$ hold the same probability. If $g$ stretches a region, the same mass is spread over a larger volume and the density falls there; if it compresses a region, the density rises. The determinant of the Jacobian measures this local change of volume (Papamakarios et al., 2021, §2.1). Rezende and Mohamed (2015, §3) describe the picture that gives the method its name: the initial density “flows” through a sequence of invertible maps, and the result is still a valid, normalized distribution. They credit the principle to Tabak and Vanden-Eijnden (2010) and Tabak and Turner (2013); Papamakarios et al. (2021, §2.4) date the term normalizing flow to Tabak and Turner and describe Rezende and Mohamed as carrying its idea and language into variational inference, after Rippel and Adams (2013) and NICE had parameterized such maps with neural networks.

A single invertible layer with a tractable determinant is not very expressive, but invertible maps compose. A deep flow stacks many simple layers, each easy to invert and with a cheap log-determinant, and the log-determinant of the stack is the sum of the layers’ terms. Most of the design effort goes into layers whose Jacobian is triangular, so that the determinant is the product of the diagonal (Kingma & Dhariwal, 2018, §2).

Mathematical Formulation

Change of variables

Let $f: \mathbb{R}^D \to \mathbb{R}^D$ be invertible and differentiable with a differentiable inverse, and let $p_Z$ be the base density. The density of $x$ is

$$\log p_X(x) = \log p_Z\big(f(x)\big) + \log \left| \det \frac{\partial f}{\partial x} \right|.$$

For a composition $f = f_K \circ \dots \circ f_1$ with intermediate values $h_0 = x$, $h_k = f_k(h_{k-1})$ and $h_K = z$,

$$\log p_X(x) = \log p_Z(z) + \sum_{k=1}^{K} \log \left| \det \frac{\partial h_k}{\partial h_{k-1}} \right|,$$

which is Glow’s Eq. 7 (Kingma & Dhariwal, 2018). The requirement that $z$ have the same dimension $D$ as $x$ follows from invertibility (Papamakarios et al., 2021, §2.1). For a layer whose Jacobian is triangular, the log-determinant is the sum of the logs of the diagonal entries, a cost linear in $D$; a general $D \times D$ determinant costs $O(D^3)$, which Papamakarios et al. (2021, §3) note can be intractable for large $D$.

Training by maximum likelihood

Given training images $x^{(1)}, \dots, x^{(N)}$, the parameters are chosen to minimize the average negative log-likelihood

$$\mathcal{L} = -\frac{1}{N} \sum_{n=1}^{N} \Big[ \log p_Z\big(f(x^{(n)})\big) + \log \left| \det \frac{\partial f}{\partial x}(x^{(n)}) \right| \Big],$$

which is, up to a constant, the forward KL divergence from the data distribution to the model (Papamakarios et al., 2021, §2.3.1, §6.1). Both terms pull in opposite directions: the first wants every image mapped to a code near the origin, the second penalizes squeezing the data into a small volume.

Expressive power

Papamakarios et al. (2021, §2.2) show that flows can in principle represent any distribution with a positive density and differentiable conditional distribution functions. The construction maps each coordinate through the cumulative distribution function of its conditional given the previous coordinates; that map is triangular and turns any such density into the uniform distribution on the unit cube. In practice the result depends on how much a finite stack of restricted layers can approximate.

Building Blocks

Coupling layers: NICE and RealNVP

A coupling layer splits its input into two parts, $x = (x_a, x_b)$, passes one part unchanged and transforms the other elementwise with parameters computed from the unchanged part. In the affine form used by Glow (Kingma & Dhariwal, 2018, Table 1),

$$y_a = s \odot x_a + t, \qquad y_b = x_b, \qquad (\log s, t) = \mathrm{NN}(x_b),$$

and the inverse is $x_b = y_b$, $x_a = (y_a - t) / s$, with the same network evaluated on $y_b$. The network $\mathrm{NN}$ can be any function, a convolutional network for images, because it is never inverted. The Jacobian is triangular with an identity block and the diagonal $s$, so the log-determinant is $\sum \log |s|$ (Papamakarios et al., 2021, §3.1.2, Eq. 45). Both directions cost one network pass, which is why, in the review’s account, coupling layers are among the most popular ways to build flows and are common in generative models of images, audio and video.

NICE (Dinh, Krueger & Bengio, 2015) used the additive special case, $s = 1$, which preserves volume, so the log-determinant is zero (Rezende & Mohamed, 2015, §5; Papamakarios, Pavlakou & Murray, 2017, §3.3). RealNVP (Dinh, Sohl-Dickstein & Bengio, 2017), whose NVP stands for non-volume preserving, added the scale (Papamakarios et al., 2017, §3.3). In the account of Papamakarios et al. (2021, §6.1.2), RealNVP was the first flow to produce sharp, visually compelling full-color images, with samples on CelebA, ImageNet at $64 \times 64$, CIFAR-10 and LSUN. Neither paper could be read here; their contents are given as Glow, MAF, IAF and the 2021 review describe them.

A coupling layer changes only half of its input, so successive layers must alternate or permute which half is passed through, otherwise part of the input is never transformed (Papamakarios et al., 2021, §3.2). For images RealNVP used two kinds of split: a spatial checkerboard pattern and a split along the channels; NICE and RealNVP reversed the channel order between layers (Kingma & Dhariwal, 2018, §3.3).

Multi-scale architecture

Because every layer keeps the full dimension, a deep flow over a large image is expensive. RealNVP’s multi-scale architecture factors out part of the variables at regular intervals: they are sent straight to the final code, and later layers act only on the rest (Papamakarios et al., 2021, §3.4). Glow uses the same design with $L$ levels, each a “squeeze”, $K$ steps of flow and a “split” (Kingma & Dhariwal, 2018, Figure 2), and refers to RealNVP for its details (§3). According to Papamakarios et al. (2021, §3.4), RealNVP showed that the dimensions that pass through every layer carry more global, semantic information, and those factored out early carry lower-level detail.

Autoregressive layers: MAF and IAF

An autoregressive model with Gaussian conditionals generates $x_i = u_i \exp(\alpha_i) + \mu_i$, with $\mu_i$ and $\alpha_i$ computed from $x_1, \dots, x_{i-1}$ and $u_i \sim \mathcal{N}(0, 1)$. Read this way, the model is an invertible map from the noise $u$ to the data $x$, with a triangular Jacobian whose log-determinant is $-\sum_i \alpha_i$ in the direction from data to noise (Papamakarios et al., 2017, §3.1, Eq. 5); Papamakarios et al. credit the observation to Kingma et al. (2016).

Masked Autoregressive Flow (MAF; Papamakarios, Pavlakou & Murray, 2017) stacks such layers, each implemented as a MADE network whose weight matrices are masked so that output $i$ sees only inputs before $i$. The density of a data point then takes one pass, but sampling takes $D$ sequential passes. Inverse Autoregressive Flow (IAF; Kingma et al., 2016) computes $\mu_i$ and $\alpha_i$ from the noise variables $u_{1:i-1}$ instead, so sampling and the density of its own samples take one pass, but the density of an arbitrary data point takes $D$ passes (Papamakarios et al., 2017, §3.2). That is why IAF was built for approximate posteriors in variational inference, where only the density of its own samples is needed, and MAF for density estimation. Papamakarios et al. also show that an affine coupling layer is a special case of both: set $\mu_i = \alpha_i = 0$ for the first $d$ variables and let the rest depend only on those $d$ (§3.3). Coupling layers give up flexibility per layer for one-pass sampling and density evaluation in both directions; Glow left MAF out of its comparisons because its sampling cannot be parallelized (Kingma & Dhariwal, 2018, §4).

Glow: actnorm and invertible 1 × 1 convolutions

Glow (Kingma & Dhariwal, 2018) builds each step of flow from three layers: actnorm, an invertible $1 \times 1$ convolution and an affine coupling layer (Figure 2, Table 1). For an $h \times w \times c$ tensor:

• Actnorm applies a per-channel scale and bias, $y_{i,j} = s \odot x_{i,j} + b$, with log-determinant $h \cdot w \cdot \sum \log |s|$. RealNVP had used batch normalization, whose noise grows as the per-device minibatch shrinks, and Glow trained large images with one image per device. Actnorm’s parameters are initialized so that the activations of an initial minibatch have zero mean and unit variance per channel, then trained as ordinary parameters independent of the data (§3.1). The normalization layers article covers the batch-dependent alternatives.
• Invertible 1 × 1 convolution replaces the fixed channel permutation with a learned $c \times c$ matrix $W$ applied at every pixel, $y_{i,j} = W x_{i,j}$, of which a permutation is a special case. Its log-determinant is $h \cdot w \cdot \log |\det W|$. $W$ starts as a random rotation, with log-determinant 0. Computing $\det W$ costs $O(c^3)$; parameterizing $W = P L (U + \mathrm{diag}(s))$, with $P$ a fixed permutation, $L$ unit lower triangular and $U$ strictly upper triangular, reduces the log-determinant to $\sum \log |s|$, though the authors did not measure a large difference in wall-clock time for their networks (§3.2).
• Affine coupling as above, with the last convolution of each coupling network initialized to zero so that every coupling layer starts as the identity, which the authors found helps train very deep networks; Glow splits only along channels, not with checkerboards (§3.3).

On CIFAR-10, with $K = 32$ steps and $L = 3$ levels and everything else fixed, the $1 \times 1$ convolution gave a lower negative log-likelihood and faster convergence than channel reversal or a fixed random shuffle, with additive and with affine coupling, averaged over three seeds; it added 0.2% parameters and about 7% wall-clock time (§5, Figure 3).

Residual flows

A residual layer $z' = z + g_\phi(z)$ is invertible when $g_\phi$ is contractive, which can be enforced by spectral normalization of each linear layer; the inverse is found by fixed-point iteration (Papamakarios et al., 2021, §3.3.1). Such layers have a dense Jacobian, so every variable can affect every other, but no efficient exact log-determinant is known. Training uses an unbiased estimate from a power series of traces, each trace estimated with Hutchinson’s estimator. Planar and radial flows, which Rezende and Mohamed (2015, §4.1) introduced for variational inference, are residual layers whose determinant follows from the matrix determinant lemma in $O(D)$ time.

Continuous-time flows

Letting the number of layers go to infinity gives a flow defined by an ordinary differential equation, $dz(t)/dt = g(z(t), t)$. Chen et al. (2018, Theorem 1) showed that the log-density then changes as

$$\frac{\partial \log p\big(z(t)\big)}{\partial t} = -\operatorname{tr}\!\left( \frac{\partial g}{\partial z(t)} \right),$$

the instantaneous change of variables, which needs a trace instead of a determinant, and they called these models continuous normalizing flows. The function $g$ needs only to be uniformly Lipschitz in $z$ and continuous in $t$; it does not itself have to be invertible, because a unique solution of the ODE makes the map from start to end automatically bijective (Chen et al., 2018, §4). Sampling and density evaluation each require solving the ODE numerically, and both directions cost the same (Papamakarios et al., 2021, §4.1). Computing the trace exactly still takes $O(D)$ backward passes; FFJORD (Grathwohl et al., 2019) estimates it with Hutchinson’s estimator, $\operatorname{tr}(J) \approx v^\top J v$ for a random vector $v$ with zero mean and unit covariance, which the review describes as about $D$ times more efficient than computing the trace exactly (Papamakarios et al., 2021, §4.1). Chen et al. (2018) also train through the solver with the adjoint method, a continuous-time analogue of backpropagation that does not store the solver’s intermediate steps. In our reading, the probability-flow ODE of diffusion models, which gives those models exact likelihoods, is a continuous-time flow of this kind whose vector field comes from a learned score.

How Data Flows

Scoring an image. The image enters the data-to-code direction $f$. Each step of flow returns its output and its log-determinant; the log-determinants are summed, the code is scored under the base Gaussian, and the two are added. In a multi-scale model, the factored-out variables are scored as well, so the whole code has $D$ dimensions. One pass gives the exact log-likelihood.

Generating an image. A code is drawn from the base distribution and passed through the inverse layers in reverse order: coupling layers subtract and divide instead of adding and multiplying, $1 \times 1$ convolutions apply $W^{-1}$, actnorm undoes its scale and bias, and squeezes are reversed. Glow reports generating a $256 \times 256$ image in about 130 ms on one NVIDIA GTX 1080 Ti (footnote 2). Glow and earlier likelihood models often gave better-looking samples at a reduced temperature $T$, sampling from a distribution proportional to $p_\theta(x)^{T^2}$; with additive coupling this amounts to multiplying the standard deviation of the base distribution by $T$, and the authors chose $T = 0.7$ (§6).

Encoding. The map from image to code is exact, unlike a VAE encoder, so every image has one point in the flow’s latent space. Glow interpolated linearly between the codes of two real images and edited images along the difference of mean codes of images with and without an attribute, computed after training from CelebA’s binary labels (§6).

Training Considerations

Dequantization

Pixel values are integers, but a flow is a density on continuous space. Fitting a continuous density to discrete points has a degenerate optimum: narrow spikes of density at the observed values give arbitrarily high likelihoods. Uria, Murray and Larochelle (2013, §4.2) avoided this “cheating” solution in RNADE by adding noise uniformly distributed between 0 and 1 to each pixel value before modelling image patches. This uniform dequantization became the usual practice (Ho et al., 2019, §3.1.1). Glow writes the objective for data with discretization level $a$ as

$$\mathcal{L} \simeq \frac{1}{N} \sum_{n} -\log p_\theta\big(x^{(n)} + u\big) + c, \qquad u \sim \mathcal{U}(0, a), \quad c = -M \log a,$$

with $M$ the number of dimensions (Kingma & Dhariwal, 2018, Eq. 2). The continuous model defines a discrete one, $P_\text{model}(x) = \int_{[0,1)^D} p_\text{model}(x + u)\, du$, and by Jensen’s inequality the expected log-density of dequantized data is a lower bound on the expected log-probability of that discrete model, an argument Ho et al. (2019, §3.1.1) attribute to Theis et al. (2016). The continuous objective therefore cannot be pushed past the discrete model’s likelihood by collapsing onto the data.

Uniform noise asks the model to spread density evenly over each unit cube around a data point, which Ho et al. (2019, §3.1.2) argue is hard for smooth function approximators. Their variational dequantization learns the noise distribution $q(u \mid x)$ with a conditional flow and maximizes $\mathbb{E}_{u \sim q}\big[\log p_\text{model}(x + u) - \log q(u \mid x)\big]$, a bound whose gap is the KL divergence between $q$ and the model’s own distribution of $u$ given $x$; uniform dequantization is the special case of a uniform $q$. In their CIFAR-10 ablations after 400 epochs, before the models had converged, switching back to uniform dequantization cost about 0.127 bits per dimension, more than replacing their coupling layers or removing self-attention (about 0.03 each), and the train–test gap was about 0.06 bits per dimension with uniform noise against about 0.02 with variational dequantization (§4.2, Table 2).

Bits per dimension

Flows report the negative log-likelihood in bits per dimension: $-\log_2 P_\text{model}(x) / D$, with $D = H \times W \times C$, averaged over the test set. For 8-bit images rescaled to $[0, 1)$, the bin width is $a = 1/256$ and the constant $c = -M \log a$ of Glow’s Eq. 2 adds $\log_2 256 = 8$ bits per dimension to the continuous negative log-density in bits. A model that assigns equal probability to all 256 values of every dimension scores exactly 8, and the uniform distribution is listed at 8.00 bits per dimension on CIFAR-10 by Kingma et al. (2016, Table 2). The dequantized estimate is an upper bound on the discrete model’s bits per dimension; Flow++’s CIFAR-10 result is a tighter value estimated by importance sampling (Ho et al., 2019, §4.1). Some flows are trained in a logit space: MAF, following RealNVP’s preprocessing, rescales the dequantized values to $[0, 1]$ and maps them through $\mathrm{logit}(\lambda + (1 - 2\lambda) x)$, and reports log-likelihoods in that space, converting to bits per pixel in its supplementary material (Papamakarios et al., 2017, §4.3). Numbers are comparable across papers only when they are converted back to the discrete pixel space.

Other practicalities

• Normalization. Batch normalization is an elementwise affine map with a diagonal Jacobian and can be inserted between flow layers; during training its statistics depend on the batch, so the formula for its log-determinant is an approximation that often suffices for large batches. Actnorm avoids the batch dependence (Papamakarios et al., 2021, §3.4).
• Initialization at the identity. Zero-initialized coupling outputs (Glow, §3.3) and a data-dependent actnorm initialization start the flow close to an identity map with normalized activations.
• Memory. Computing gradients in reversible networks needs an amount of memory that is constant rather than linear in depth (Kingma & Dhariwal, 2018, §1, citing RevNets); Glow itself used gradient checkpointing at high resolution and named the reversible approach as future work (§6).
• Depth. In the review’s account, Glow composed as many as 320 sub-transformations, trained across 40 GPUs (Papamakarios et al., 2021, §3.4).

Example

A coupling flow on $6 \times 6$ images with four grey levels. Each image is one of 82 clean pictures, a square of side 2 or 3 with intensity 2 or 3 on a background of 0, after which every pixel is, with probability 0.1, replaced by a level drawn uniformly from $\{0, 1, 2, 3\}$. Because the generator is known, the exact probability of any image is a sum over the 82 clean pictures, and the true entropy of the image distribution can be estimated on test images: no discrete model can average fewer bits than that. As a second reference, a model of independent pixels, fitted by counting each pixel’s levels in the training images, gives exact discrete bits per dimension. A small RealNVP-style flow, eight affine coupling layers with alternating checkerboard masks, each conditioned by a fully connected network with two hidden layers of 128 units, maps an image (with the constant 2 subtracted) to a 36-dimensional Gaussian code. The log-scale is bounded by a tanh times a learned factor, and the last layer of each conditioner is zero-initialized, as in Glow, so every coupling starts as the identity. Pixel values stay on the integer grid, so the bin width is $a = 1$ and Glow’s constant $c$ is 0: bits per dimension are $-\log_2$ of the density at the evaluated points divided by 36.

Three copies of the flow are trained identically on 20,000 images except for their inputs: uniformly dequantized values $x + u$; the integer values moved to the centres of their cells, $x + 0.5$; and the integer values themselves. For each, the script prints the reported number, the negative log-density at the kind of input the model was trained on, and the probability its density gives to the discrete image, $P(x) = \int_{[0,1)^{36}} p(x + u)\, du$, estimated from one or from 256 uniform draws of $u$. It then checks the change-of-variables log-determinant of the dequantized model against a brute-force determinant of its full $36 \times 36$ Jacobian, computed by automatic differentiation, and inverts the flow in both directions. It uses PyTorch 2.14.1 and NumPy 2.5.4 in double precision on one thread, and runs in about four minutes. The budget of 150 epochs per model was not part of the original plan: in a pilot at the planned 40 epochs the dequantized flow did not beat the independent-pixel model, so one further pilot at 150 epochs was run with a criterion fixed in advance, that the dequantized model’s 256-draw score fall below the independent-pixel model’s, and the criterion was met. After the first run, the script was extended to print the reported bits per dimension of the two undequantized models on the test images every 10 epochs; the evaluation draws no random numbers, and the rest of the output is unchanged.

import math
import numpy as np
import torch

torch.set_num_threads(1)
torch.set_default_dtype(torch.float64)
S, K, EPS = 6, 4, 0.1          # 6x6 images, 4 grey levels, pixel noise rate
D = S * S
rng = np.random.default_rng(0)

# 82 clean images: a square of side 2 or 3, intensity 2 or 3, on background 0.
clean = []
for side in (2, 3):
    for r in range(S - side + 1):
        for c in range(S - side + 1):
            for v in (2, 3):
                im = np.zeros((S, S), int); im[r:r + side, c:c + side] = v
                clean.append(im.ravel())
clean = np.array(clean)

def sample(n):                  # a clean image, then each pixel resampled with prob. EPS
    x = clean[rng.integers(0, len(clean), n)].copy()
    flip = rng.random((n, D)) < EPS
    x[flip] = rng.integers(0, K, flip.sum())
    return x

def log2_true(x):               # exact log2 P(x) under the generator
    same = x[:, None, :] == clean[None, :, :]
    lp = np.where(same, np.log(1 - EPS + EPS / K), np.log(EPS / K)).sum(2)
    return (np.logaddexp.reduce(lp, axis=1) - np.log(len(clean))) / np.log(2)

train, test = sample(20000), sample(5000)
H = -log2_true(test).mean() / D
m = np.stack([(train == k).mean(0) for k in range(K)])      # independent pixels
indep = -np.log2(m[test, np.arange(D)]).sum(1).mean() / D
print(f"true entropy {H:.3f} bits/dim; independent-pixel model {indep:.3f} bits/dim")

class Coupling(torch.nn.Module):    # affine coupling: y = x*exp(s) + t on unmasked pixels
    def __init__(self, mask):
        super().__init__()
        self.register_buffer("mask", mask)
        self.net = torch.nn.Sequential(torch.nn.Linear(D, 128), torch.nn.ReLU(),
                                       torch.nn.Linear(128, 128), torch.nn.ReLU(),
                                       torch.nn.Linear(128, 2 * D))
        torch.nn.init.zeros_(self.net[-1].weight); torch.nn.init.zeros_(self.net[-1].bias)
        self.scale = torch.nn.Parameter(torch.ones(1))   # zero init: starts as identity
    def st(self, xa):
        h = self.net(xa)
        s = self.scale * torch.tanh(h[:, :D])
        return s * (1 - self.mask), h[:, D:] * (1 - self.mask)
    def forward(self, x):           # image -> code, with log|det J|
        xa = x * self.mask; s, t = self.st(xa)
        return xa + (1 - self.mask) * (x * torch.exp(s) + t), s.sum(1)
    def inverse(self, z):
        za = z * self.mask; s, t = self.st(za)
        return za + (1 - self.mask) * ((z - t) * torch.exp(-s))

board = torch.tensor(((np.arange(S)[:, None] + np.arange(S)) % 2).ravel(), dtype=torch.float64)

class Flow(torch.nn.Module):        # 8 couplings with alternating checkerboard masks
    def __init__(self, L=8):
        super().__init__()
        self.layers = torch.nn.ModuleList(
            [Coupling(board if i % 2 == 0 else 1 - board) for i in range(L)])
    def forward(self, x):
        x, logdet = x - 2.0, torch.zeros(len(x))
        for f in self.layers:
            x, ld = f(x); logdet = logdet + ld
        return x, logdet
    def inverse(self, z):
        for f in reversed(self.layers): z = f.inverse(z)
        return z + 2.0
    def log_prob(self, x):          # natural-log density by change of variables
        z, ld = self(x)
        return -0.5 * (z ** 2).sum(1) - 0.5 * D * math.log(2 * math.pi) + ld

# Three training inputs: dequantized x+u; cell centres x+0.5; grid points x.
SHIFT = {"dequantized": None, "cell centres": 0.5, "grid points": 0.0}
def fit(kind, epochs=150):
    torch.manual_seed(1)
    flow = Flow(); opt = torch.optim.Adam(flow.parameters(), lr=1e-3)
    X = torch.tensor(train, dtype=torch.float64)
    trace = []                      # undequantized models: reported bits/dim every 10 epochs
    for ep in range(epochs):
        perm = torch.randperm(len(X))
        for i in range(0, len(X), 256):
            xb = X[perm[i:i + 256]]
            xb = xb + (torch.rand_like(xb) if SHIFT[kind] is None else SHIFT[kind])
            loss = -flow.log_prob(xb).mean()
            opt.zero_grad(); loss.backward(); opt.step()
        if SHIFT[kind] is not None and ep % 10 == 9:
            with torch.no_grad(): trace.append(bits(flow.log_prob(Xt + SHIFT[kind])))
    return flow, trace

Xt = torch.tensor(test, dtype=torch.float64)
bits = lambda lp: -lp.mean().item() / D / math.log(2)      # nats -> bits per dimension
print(f"{'trained on':14s}{'reported':>10s}{'discrete, S=1':>15s}{'discrete, S=256':>17s}")
flows, traces = {}, {}
for kind in SHIFT:
    flow, traces[kind] = fit(kind); flows[kind] = flow
    with torch.no_grad():
        y = Xt + (torch.rand_like(Xt) if SHIFT[kind] is None else SHIFT[kind])
        g = torch.Generator().manual_seed(2)
        lps = torch.stack([flow.log_prob(Xt + torch.rand(Xt.shape, generator=g))
                           for _ in range(256)])           # p(x+u), u ~ U[0,1)^D
    print(f"{kind:14s}{bits(flow.log_prob(y).detach()):10.3f}{bits(lps[0]):15.3f}"
          f"{bits(torch.logsumexp(lps, 0) - math.log(256)):17.3f}")
print("reported bits/dim of the undequantized models during training")
print(f"{'epoch':>5s}{'cell centres':>14s}{'grid points':>13s}")
for i, (a, b) in enumerate(zip(traces["cell centres"], traces["grid points"])):
    print(f"{10 * (i + 1):5d}{a:14.3f}{b:13.3f}")

f = flows["dequantized"]
x8 = Xt[:8] + 0.3               # brute force: full 36x36 Jacobian, then its determinant
J = torch.stack([torch.autograd.functional.jacobian(lambda v: f(v[None])[0][0], xi) for xi in x8])
_, ld_brute = torch.linalg.slogdet(J)
with torch.no_grad():
    _, ld = f(x8)
    print(f"log|det J|: coupling sum vs brute force, max difference {(ld - ld_brute).abs().max().item():.1e}")
    z, _ = f(Xt + 0.3)
    print(f"image -> code -> image, max error {(f.inverse(z) - (Xt + 0.3)).abs().max().item():.1e}")
    zs = torch.randn(1000, D, generator=torch.Generator().manual_seed(3))
    print(f"code -> image -> code, max error {(f(f.inverse(zs))[0] - zs).abs().max().item():.1e}")

Output:

true entropy 0.684 bits/dim; independent-pixel model 1.058 bits/dim
trained on      reported  discrete, S=1  discrete, S=256
dequantized        1.116          1.118            0.979
cell centres       2.042         28.572            2.258
grid points       -0.003          3.051            2.113
reported bits/dim of the undequantized models during training
epoch  cell centres  grid points
   10        -0.242        0.015
   20        -0.251        0.062
   30        -0.486        0.085
   40        -0.312        1.429
   50         1.317        0.419
   60        -0.603       -0.134
   70         0.708       -0.238
   80         0.345        0.011
   90         0.744       -0.039
  100         0.510       -0.227
  110        -0.141        0.148
  120         0.003        1.184
  130        -0.150        0.158
  140         0.076       -0.383
  150         2.042       -0.003
log|det J|: coupling sum vs brute force, max difference 1.1e-14
image -> code -> image, max error 1.3e-14
code -> image -> code, max error 1.3e-14

The exactness checks agree to about $10^{-14}$, the round-off of double precision: the sum of the coupling layers’ log-scales is the log-determinant of the whole network’s Jacobian, and the inverse layers undo the forward ones in both directions, on test images and on random codes. That is what the architecture guarantees, whatever the training did.

The dequantized model behaves as the bound says. Its reported 1.116 bits per dimension and the one-draw estimate 1.118 are the same quantity with different noise; averaging the density over 256 draws before taking the log gives 0.979 bits per dimension, tighter by 0.14. Since the logarithm of an average of draws underestimates the log of the expectation on average, 0.979 overstates the discrete model’s bits per dimension on average, and the model is at least that good. That is better than the independent-pixel model (1.058), which cannot represent a square, and above the true entropy (0.684), as any discrete model’s score must be on average. By the identity Ho et al. (2019, §3.1.2) give, the gap between the one-draw bound and the discrete model’s score is the KL divergence between uniform noise and the model’s own distribution of $u$ within each cell, so uniform dequantization costs this model, on average, at least about the 0.14 bits per dimension seen here: it does not spread its density evenly over each cell.

The undequantized models show why the reported number means nothing without dequantization. Their reported values never settle: from one 10-epoch checkpoint to the next they jump by up to about 2 bits per dimension, and most checkpoints of both models lie below the true entropy, many below zero, values no discrete model can reach, because the change-of-variables formula lets a continuous density grow without limit near the training points. Where training stops decides what is reported: this run ends with the grid-point model just below zero and the cell-centre model at its highest value, above the entropy, after a checkpoint near zero. As models of the discrete images, both score worse than independent pixels at the end of training (2.113 and 2.258 bits per dimension). The cell-centre model’s one-draw estimate of 28.572 against 2.258 with 256 draws shows that it puts its mass on a small part of each cell, in our reading near the centres it was trained on. Neither undequantized model’s reported number says how well it models the images, which is the point of RNADE’s added noise and of Theis et al.’s bound as Ho et al. (2019) present it.

Variants

• Planar and radial flows (Rezende & Mohamed, 2015): residual layers with $O(D)$ determinants, used to make approximate posteriors in variational inference more flexible than a diagonal Gaussian.
• IAF (Kingma et al., 2016): autoregressive posterior flows inside a VAE; their ResNet VAE with IAF reached 3.11 bits per dimension on CIFAR-10, a likelihood estimated by importance sampling, and sampled in about 0.05 s per image against 52.0 s for a naively sampled PixelCNN on a Titan X (§6.2, Table 2). The variational autoencoder article places IAF among VAE variants.
• Conditional flows: a condition $y$ is fed to every layer’s conditioner network. Conditional MAF on CIFAR-10 with class labels reached 2.98 bits per pixel in the main text’s account (Papamakarios et al., 2017, §4.3).
• Discrete flows act on integer-valued variables directly. They need no Jacobian term, but a bijection between discrete sets can only permute probabilities, which limits what a factorized base distribution can represent (Papamakarios et al., 2021, §5.3).
• Residual and continuous flows (above) trade one-pass exact evaluation for dense Jacobians and free-form networks.

Use in Computer Vision

Image generation and density estimation. Coupling flows give exact likelihoods and one-pass sampling on natural images, from RealNVP to Glow’s $256 \times 256$ faces trained on 5-bit CelebA-HQ images (Kingma & Dhariwal, 2018, §6). Their samples trailed those of GANs and later diffusion models in FID; the generative model article gives the comparison, and image generation covers the task.

Anomaly detection. A flow fitted to normal data gives a likelihood that can flag unusual inputs. In visual anomaly detection and industrial visual inspection the flow is fitted to features of a pretrained network rather than to pixels: DifferNet (Rudolph, Wandt & Rosenhahn, 2021) maps pooled multi-scale features to a normal distribution and scores images by their likelihood, and CFLOW-AD (Gudovskiy, Ishizaka & Kozuka, 2022) models the features at every position with conditional flows for pixel-level localization. The anomaly detection article gives their architectures and MVTec AD results. Papamakarios et al. (2021, §6.1.1) trace the idea back to Gaussianization used for one-class detection of urban areas in satellite images.

Pose and keypoint likelihoods. Residual log-likelihood estimation (Li et al., 2021) uses a small RealNVP flow over two-dimensional keypoint positions to learn the shape of the regression error density for human pose estimation; the predicted mean is the output and the flow is not run at test time, so inference costs the same as plain regression (§1, §3.3). The heatmap regression article gives its COCO results. Li et al. also note earlier work that used flows as priors over 3D human poses and body-model joint angles (§2).

Variational inference. Flows make the approximate posteriors of latent-variable models more flexible, the use for which Rezende and Mohamed (2015) applied them and for which IAF was designed. Flow++ (Ho et al., 2019) uses a conditional flow as the dequantization distribution itself.

Limitations

• Same dimension everywhere. The code is as large as the image, so a flow cannot compress to a low-dimensional latent the way an autoencoder does, and every layer of a deep flow acts on the full resolution unless variables are factored out (Papamakarios et al., 2021, §3.4).
• Restricted layers. Coupling layers buy one-pass evaluation and sampling with layers that are individually weak: a single coupling layer is not a universal approximator, and whether fewer than $O(D)$ coupling layers can be is an open problem (Papamakarios et al., 2021, §3.1.2). Autoregressive layers are more flexible but slow in one direction; residual and continuous flows lose exact one-pass evaluation.
• Likelihood is not sample quality. Flow likelihoods on CIFAR-10 came close to autoregressive models (3.08 bits per dimension for Flow++ against 2.92 for PixelCNN++ in Ho et al., 2019, Table 1), but good likelihoods neither require nor guarantee good samples; the generative model article discusses why.
• Dequantization is part of the model. The reported number depends on the noise used and how the bound is estimated, and a continuous density evaluated on undequantized integers is not a likelihood at all, as the example shows.

Important Models

Model	Year	Layer type	Key change
Planar and radial flows (Rezende & Mohamed)	2015	Residual	Flows for variational posteriors, $O(D)$ determinants
NICE (Dinh, Krueger & Bengio)	2015	Additive coupling	Volume-preserving coupling layers
IAF (Kingma et al.)	2016	Inverse autoregressive	Scalable posterior flows in a VAE
RealNVP (Dinh, Sohl-Dickstein & Bengio)	2017	Affine coupling	Scale term, checkerboard and channel masks, multi-scale
MAF (Papamakarios, Pavlakou & Murray)	2017	Masked autoregressive	Stacked MADEs for density estimation
Glow (Kingma & Dhariwal)	2018	Affine coupling	Actnorm, invertible $1 \times 1$ convolution, high-resolution faces
Neural ODE CNF (Chen et al.)	2018	Continuous time	Instantaneous change of variables
FFJORD (Grathwohl et al.)	2019	Continuous time	Hutchinson trace estimator for free-form dynamics
Flow++ (Ho et al.)	2019	Logistic-mixture coupling	Variational dequantization, self-attention

Published unconditional results in bits per dimension on test data (lower is better):

Model	CIFAR-10	ImageNet 32 × 32	ImageNet 64 × 64	Source
NICE	4.48	—	—	Kingma et al. (2016), Table 2
RealNVP	3.49	4.28	3.98	Kingma & Dhariwal (2018), Table 2
Glow	3.35	4.09	3.81	Kingma & Dhariwal (2018), Table 2
Flow++	3.08	3.86	3.69	Ho et al. (2019), Table 1

Glow’s Table 2 also gives LSUN at $64 \times 64$ (random crops of images downsampled to $96 \times 96$): 2.38, 2.46 and 2.67 bits per dimension for bedroom, tower and church outdoor, against 2.72, 2.81 and 3.08 for RealNVP. Its text says the table also includes the bits per dimension of the $256 \times 256$ CelebA-HQ model, but the table in the NeurIPS version has no such column. Flow++’s CIFAR-10 figure is an importance-sampled estimate with 16,384 samples; with one sample, as in the dequantized bound, it is 3.12 (§4.1). Flow++ lists RealNVP’s ImageNet $64 \times 64$ entry as not available, while Glow’s table gives 3.98.