Evidence Lower Bound (ELBO)

The evidence lower bound (ELBO) is a quantity that can be computed and maximized in place of the log-likelihood $\log p_\theta(x)$ of an image when that likelihood is an intractable integral over hidden variables. It is the training objective of variational autoencoders and the objective that diffusion models are derived from. The “evidence” is the probability of the observed data, $p_\theta(x)$, and the ELBO falls short of its logarithm by exactly the KL divergence between an approximate posterior and the true one.

The generative model article places the ELBO among the model families; this article derives the bound, explains when it is tight, and checks each claim on small digit images whose likelihood is exact.

Definition

Let a latent-variable model generate an image $x$ by first drawing a latent code $z$ from a prior $p(z)$ and then drawing $x$ from a decoder $p_\theta(x \mid z)$. For any distribution $q_\phi(z \mid x)$ over codes, called the approximate posterior or encoder, the evidence lower bound is

$$\mathcal{L}(\theta, \phi; x) = \mathbb{E}_{q_\phi(z \mid x)}\big[\log p_\theta(x \mid z)\big] - \mathrm{KL}\big(q_\phi(z \mid x) \,\Vert\, p(z)\big) \le \log p_\theta(x).$$

In a variational autoencoder, the prior is $p(z) = \mathcal{N}(0, I)$, and both $q_\phi$ and $p_\theta$ are neural networks.

Intuition

To score an image under a latent-variable model, one would have to average $p_\theta(x \mid z)$ over every code the prior could produce. For an image of a handwritten “3”, almost all codes explain it terribly; only a small set, the codes of 3s in roughly the right place, matter. The posterior $p_\theta(z \mid x)$ is exactly the distribution over those relevant codes, but computing it requires the same intractable integral.

The ELBO replaces it with a guess $q_\phi(z \mid x)$, an encoder that proposes codes for the image. The bound rewards proposals that reconstruct the image well and penalizes proposals that stray from the prior. A perfect guess makes the bound exact; an imperfect one costs its KL divergence from the true posterior.

Formal Definition

Derivation by Jensen’s inequality

Write the marginal likelihood as an expectation under $q_\phi$ by multiplying and dividing by it, and move the logarithm inside, which Jensen’s inequality allows because $\log$ is concave:

$$\log p_\theta(x) = \log \mathbb{E}_{q_\phi(z \mid x)}\!\left[\frac{p_\theta(x, z)}{q_\phi(z \mid x)}\right] \ge \mathbb{E}_{q_\phi(z \mid x)}\big[\log p_\theta(x, z) - \log q_\phi(z \mid x)\big] = \mathcal{L}(\theta, \phi; x).$$

Jordan et al. (1999, §6, Eq. 43) give this derivation for graphical models with discrete hidden nodes, and Burda, Grosse and Salakhutdinov (2016, Eq. 3, arXiv v4) for the variational autoencoder.

Derivation from the posterior

The second derivation gives the gap exactly. Since $p_\theta(x, z) = p_\theta(z \mid x)\, p_\theta(x)$, and $\log p_\theta(x)$ does not depend on $z$,

$$\log p_\theta(x) - \mathcal{L}(\theta, \phi; x) = \mathbb{E}_{q_\phi(z \mid x)}\big[\log q_\phi(z \mid x) - \log p_\theta(z \mid x)\big] = \mathrm{KL}\big(q_\phi(z \mid x) \,\Vert\, p_\theta(z \mid x)\big),$$

so

$$\mathcal{L}(\theta, \phi; x) = \log p_\theta(x) - \mathrm{KL}\big(q_\phi(z \mid x) \,\Vert\, p_\theta(z \mid x)\big).$$

Kingma and Welling (2014, §2.2, Eq. 1, arXiv v11) start from this identity: the KL divergence is non-negative (cross-entropy loss covers its properties), so $\mathcal{L}$ is a lower bound. Splitting $\log p_\theta(x, z) = \log p_\theta(x \mid z) + \log p(z)$ in the Jensen form gives the reconstruction-minus-KL form of the definition (their Eq. 3).

Tightness and the gap

The bound is tight, $\mathcal{L} = \log p_\theta(x)$, exactly when $q_\phi(z \mid x) = p_\theta(z \mid x)$. In the Jensen derivation this is the equality condition: the ratio $p_\theta(x, z) / q_\phi(z \mid x)$ must be constant in $z$, which holds only for the posterior.

• With $\theta$ fixed, maximizing $\mathcal{L}$ over $\phi$ is the same as minimizing $\mathrm{KL}(q_\phi \Vert p_\theta(z \mid x))$, because $\log p_\theta(x)$ does not depend on $\phi$. If the family of $q_\phi$ cannot represent the posterior, for example because it factorizes over latent variables that the posterior couples, a gap remains at the optimum: the bound is the best available within the family, and exact only if $q$ may range over all distributions (Jordan et al., 1999, §6.1).
• Maximizing $\mathcal{L}$ over $\theta$ either raises $\log p_\theta(x)$ or shrinks the gap. Burda et al. (2016, §1–2) point out that this trade-off pushes the model toward representations whose posteriors are easy to approximate, which helps inference but limits how expressive the model becomes.

When $q$ may be any distribution and $\theta$ is updated with $q$ fixed, alternating the two maximizations is the expectation–maximization (EM) algorithm: Jordan et al. (1999, §6.2), following Neal and Hinton, describe EM as coordinate ascent on this bound.

Properties

Reconstruction and regularization

The two terms of the definition have the reading of an autoencoder’s loss (Kingma & Welling, 2014, §2.3). The first, $\mathbb{E}_{q}[\log p_\theta(x \mid z)]$, is a negative reconstruction error: with a Bernoulli decoder over binary pixels it is the negative binary cross-entropy summed over pixels, and with a Gaussian decoder of fixed variance it is a negative squared error up to constants. The second, $\mathrm{KL}(q_\phi(z \mid x) \Vert p(z))$, keeps each image’s codes close to the prior, so that codes drawn from the prior at generation time decode to plausible images. Rezende, Mohamed and Wierstra (2014, §4.1) minimize the negative bound, which they call the free energy, and read its terms the same way.

The Gaussian KL in closed form

With a diagonal Gaussian encoder $q_\phi(z \mid x) = \mathcal{N}(\mu, \operatorname{diag}(\sigma^2))$ in $J$ dimensions, where $\mu$ and $\sigma$ are outputs of the encoder network, and the prior $p(z) = \mathcal{N}(0, I)$, the KL term needs no sampling (Kingma & Welling, 2014, Appendix B, arXiv v11):

$$\mathrm{KL}\big(q_\phi(z \mid x) \,\Vert\, p(z)\big) = \frac{1}{2} \sum_{j=1}^{J} \left( \mu_j^2 + \sigma_j^2 - 1 - \log \sigma_j^2 \right).$$

Only the reconstruction term must then be estimated by sampling, an estimator that the authors note typically has lower variance than sampling the whole bound (§2.3).

The reparameterization trick

The reconstruction term is an expectation over a distribution that depends on $\phi$, so its gradient with respect to $\phi$ cannot be taken by differentiating samples directly; the generic score-function estimator works but has high variance (Kingma & Welling, 2014, §2.2). The reparameterization trick writes the sample as a differentiable function of $\phi$ and parameter-free noise:

$$z = \mu_\phi(x) + \sigma_\phi(x) \odot \epsilon, \qquad \epsilon \sim \mathcal{N}(0, I).$$

The expectation becomes one over $\epsilon$, a single sample gives an unbiased estimate of $\mathcal{L}$ and of its gradient, and backpropagation carries the gradient through $z$ into the encoder (§2.4). Kingma and Welling found one sample per image enough with minibatches of 100, so the bound is trained with ordinary stochastic gradient descent. Rezende et al. (2014, §3) derive gradient rules for Gaussian latents, which they call stochastic backpropagation. The trick needs a continuous latent; the example below enumerates its discrete latent instead.

Importance-weighted bounds

Averaging several importance weights inside the logarithm gives a family of tighter bounds. With $z_1, \dots, z_K$ drawn independently from $q_\phi(z \mid x)$, Burda et al. (2016, Eq. 8) define

$$\mathcal{L}_K(x) = \mathbb{E}_{z_1, \dots, z_K \sim q_\phi(z \mid x)}\!\left[\log \frac{1}{K} \sum_{k=1}^{K} \frac{p_\theta(x, z_k)}{q_\phi(z_k \mid x)}\right].$$

The average of the weights is an unbiased estimate of $p_\theta(x)$, so Jensen’s inequality again gives $\mathcal{L}_K \le \log p_\theta(x)$. $\mathcal{L}_1$ is the ELBO, and their Theorem 1 shows that $\log p_\theta(x) \ge \mathcal{L}_{K+1} \ge \mathcal{L}_K$, with $\mathcal{L}_K \to \log p_\theta(x)$ as $K \to \infty$ when the weights are bounded. Training on $\mathcal{L}_K$ gives the importance weighted autoencoder (IWAE), which has the same networks as a VAE. On binarized MNIST with two stochastic layers and $K = 50$, they report a test negative log-likelihood of 82.90 nats for the IWAE against 84.78 for a VAE trained with the same $K$, both estimated with $\mathcal{L}_{5000}$ (Table 1), so the log-likelihoods behind these numbers are themselves stochastic lower bounds, as the authors note. The IWAE also used more of its latent dimensions.

The ELBO in diffusion models

Diffusion models are latent-variable models whose latents are the noisy images $x_1, \dots, x_T$. Sohl-Dickstein et al. (2015, §2.4, arXiv v8) trained them on a Jensen lower bound whose derivation, they note, parallels variational Bayes. Ho, Jain and Abbeel (2020, arXiv v2) start from the same bound (Eq. 3), with the forward noising process $q(x_{1:T} \mid x_0)$ as a fixed approximate posterior without learnable parameters, and rewrite it as a decoder term plus KL divergences between Gaussians, computable in closed form (Eq. 5). Their best samples came from a simplified noise-prediction loss (Eq. 14), which they describe as a weighted variational bound that down-weights the terms for small noise levels. The diffusion model and DDPM articles cover this in detail.

History

The bound predates deep learning as the basis of variational methods for approximate inference; Jordan et al. (1999) gave a tutorial account for graphical models, with the Jensen derivation, the KL gap and the link to EM. In 2014, Kingma and Welling, and Rezende, Mohamed and Wierstra, made it the objective of deep generative models: a recognition network amortizes the approximate posterior across images, and reparameterized gradients train it jointly with the decoder. Importance weighting (Burda et al., 2016) and diffusion models followed.

Examples

Exact ELBO, gap and IWAE on digit glyphs. The images are $8 \times 8$ binary pictures of a $3 \times 5$ digit glyph. The latent $z$ is the digit and the glyph’s position (10 digits × 24 positions = 240 states, uniform prior, a discrete stand-in for the Gaussian prior of a VAE), and the decoder turns on the glyph’s pixels with probability 0.9 and the others with probability 0.1. With 240 states, $\log p_\theta(x)$, the posterior and every expectation under $q$ are exact sums. The script compares three choices of $q$: the prior, a mean-field $q(d)\,q(\text{position})$ fitted per image by coordinate ascent on the ELBO, and the exact posterior. It then estimates $\mathcal{L}_K$ by sampling, over 20 repetitions with the standard error shown.

import sys
import numpy as np
from scipy.special import logsumexp

rng = np.random.default_rng(int(sys.argv[1]) if len(sys.argv) > 1 else 0)

# 3x5 digit glyphs. Latent z = (digit d, top-left position r, c) of the glyph on
# an 8x8 canvas; every pixel then flips with probability 0.1.
FONT = ["111101101101111", "010110010010111", "111001111100111", "111001111001111",
        "101101111001001", "111100111001111", "111100111101111", "111001001001001",
        "111101111101111", "111101111001111"]
H = W = 8; FLIP = 0.1
glyphs = np.array([[int(b) for b in s] for s in FONT]).reshape(10, 5, 3)
templ = []                                  # one clean image per latent state
for d in range(10):
    for r in range(H - 5 + 1):
        for c in range(W - 3 + 1):
            img = np.zeros((H, W)); img[r:r + 5, c:c + 3] = glyphs[d]
            templ.append(img.ravel())
templ = np.array(templ)                     # (Z, 64): Z = 10 digits x 4 x 6 positions
Z = len(templ); n_pos = Z // 10
mu = np.where(templ == 1, 1 - FLIP, FLIP)   # decoder p(x|z): independent Bernoulli pixels
log_prior = np.full(Z, -np.log(Z))          # p(z) uniform over the 240 states

x = (rng.random((2000, 64)) < mu[rng.integers(0, Z, 2000)]).astype(float)
log_lik = x @ np.log(mu).T + (1 - x) @ np.log(1 - mu).T   # log p(x|z), (N, Z)
log_joint = log_lik + log_prior
log_px = logsumexp(log_joint, axis=1)                       # exact log p(x)
log_post = log_joint - log_px[:, None]                      # exact log p(z|x)

gaps = {}
def report(name, log_q):
    q = np.exp(log_q)
    recon = (q * log_lik).sum(1)                            # E_q[log p(x|z)]
    kl_prior = (q * (log_q - log_prior)).sum(1)             # KL(q || p(z))
    elbo = recon - kl_prior
    kl_post = (q * (log_q - log_post)).sum(1)               # KL(q || p(z|x))
    err = np.abs(log_px - elbo - kl_post).max()
    gaps[name] = log_px - elbo
    print(f"{name:22s}{recon.mean():9.3f}{kl_prior.mean():8.3f}{elbo.mean():9.3f}"
          f"{round((log_px - elbo).mean(), 3) + 0.0:8.3f}{kl_post.mean():8.3f}{err:10.1e}")

def mean_field(n_iter=50):
    """Fit q(d) q(pos) by coordinate ascent on the ELBO: each factor is set to
    the exponentiated expected log joint under the other factor."""
    lj = log_joint.reshape(-1, 10, n_pos)
    lq_pos = np.full((len(x), n_pos), -np.log(n_pos))
    for _ in range(n_iter):
        lq_d = (lj * np.exp(lq_pos)[:, None, :]).sum(2)
        lq_d -= logsumexp(lq_d, axis=1, keepdims=True)
        lq_pos = (lj * np.exp(lq_d)[:, :, None]).sum(1)
        lq_pos -= logsumexp(lq_pos, axis=1, keepdims=True)
    return (lq_d[:, :, None] + lq_pos[:, None, :]).reshape(-1, Z)

print(f"mean exact log p(x): {log_px.mean():.3f} nats over {len(x)} images, {Z} latent states")
print(f"{'q(z|x)':22s}{'recon':>9s}{'KL(q|p)':>8s}{'ELBO':>9s}{'gap':>8s}{'KL post':>8s}{'max err':>10s}")
log_q_prior = np.broadcast_to(log_prior, log_joint.shape)
log_q_mf = mean_field()
report("prior p(z)", log_q_prior)
report("mean-field q(d)q(pos)", log_q_mf)
report("true posterior", log_post)
g = gaps["mean-field q(d)q(pos)"]
print(f"mean-field gap above 0.1 nats on {np.mean(g > 0.1):.1%} of images, "
      f"which hold {g[g > 0.1].sum() / g.sum():.1%} of the total gap")

# Importance-weighted bound L_K = E[log (1/K) sum_k p(x,z_k)/q(z_k|x)], z_k ~ q.
def iwae(log_q, K, reps=20):
    rows = np.arange(len(x))[:, None]
    cdf = (np.exp(log_q).cumsum(1) / np.exp(log_q).sum(1, keepdims=True) + rows).ravel()
    out = np.empty(reps)
    for i in range(reps):                   # inverse-CDF sampling of K latents per image
        u = rng.random((len(x), K)) + rows
        z = np.minimum(np.searchsorted(cdf, u) - rows * Z, Z - 1)
        log_w = np.take_along_axis(log_joint - log_q, z, axis=1)
        out[i] = (logsumexp(log_w, axis=1) - np.log(K)).mean()
    return out.mean(), out.std() / np.sqrt(reps)

print(f"{'K':>5s}{'L_K, q = prior':>22s}{'L_K, q = mean-field':>24s}")
for K in [1, 10, 100, 1000]:
    a, sa = iwae(log_q_prior, K); b, sb = iwae(log_q_mf, K)
    print(f"{K:5d}{a:14.3f} ± {sa:5.3f}{b:16.3f} ± {sb:5.3f}")

Output:

mean exact log p(x): -25.906 nats over 2000 images, 240 latent states
q(z|x)                    recon KL(q|p)     ELBO     gap KL post   max err
prior p(z)              -50.143   0.000  -50.143  24.236  24.236   2.1e-14
mean-field q(d)q(pos)   -21.701   4.853  -26.553   0.647   0.647   2.2e-13
true posterior          -21.057   4.850  -25.906   0.000   0.000   1.4e-13
mean-field gap above 0.1 nats on 11.7% of images, which hold 99.7% of the total gap
    K        L_K, q = prior     L_K, q = mean-field
    1       -50.126 ± 0.034         -26.553 ± 0.001
   10       -36.747 ± 0.038         -26.549 ± 0.000
  100       -28.091 ± 0.013         -26.546 ± 0.000
 1000       -26.023 ± 0.003         -26.539 ± 0.001

In every row the gap $\log p_\theta(x) - \mathcal{L}$ equals the KL divergence to the true posterior, image by image, up to rounding error ($2.2 \times 10^{-13}$ nats at most). Using the prior as $q$ pays nothing in the KL term but reconstructs badly, because most of its codes put the wrong digit in the wrong place; its gap is 24.236 nats. The mean-field $q$ gets within 0.647 nats, and 99.7% of that gap sits on 11.7% of the images; in our reading these are images where noise makes digit and position jointly ambiguous, a coupling that a product of two factors cannot express. The exact posterior closes the gap completely even though its KL from the prior is 4.850 nats: the KL term of the objective is not the gap.

Importance weighting tightens the prior’s bound from $-50.126$ to $-26.023$ nats at $K = 1000$, within 0.117 nats of the exact value, and the $K = 1$ estimate matches the exact ELBO within its standard error. For the mean-field $q$ it gains only 0.014 nats. In our reading, the reason is that the mean-field fit is concentrated on one interpretation of each ambiguous image and rarely samples the other, so extra samples bring little new information; with $K = 1000$ the broad prior becomes the better proposal. Run with seeds 0 to 8 as its command-line argument (seed 0 is the default), the script gives a mean-field gap of 0.47 to 0.67 nats, a remaining prior gap at $K = 1000$ of 0.105 to 0.117 nats, and a mean-field IWAE gain of 0.014 to 0.028 nats. The script runs in about 5 seconds.

Common Misconceptions

• “The ELBO is the log-likelihood.” It is a lower bound, and the gap depends on the encoder as well as the model. Comparing two models by their ELBOs compares their bounds; a tighter estimate such as $\mathcal{L}_K$ with large $K$ is used when likelihoods are reported.
• “The KL term in the ELBO measures how loose the bound is.” The objective contains the KL divergence from $q$ to the prior; the gap is the KL divergence from $q$ to the posterior, as the example’s exact posterior shows.

Where It Is Used

The ELBO is the training objective of the variational autoencoder and its hierarchical and importance-weighted variants, the source of the likelihood bounds reported for VAEs and diffusion models, and the objective from which the DDPM loss is derived. Its encoder–decoder reading connects it to the autoencoder, and its reconstruction term to the loss functions of image models.