From Variational Bayes to Diffusion Models: The Journey of Generative Mathematics

1. Introduction

Denoising Diffusion Probabilistic Models feel, at first glance, like a strange piece of engineering: corrupt an image with a thousand tiny doses of Gaussian noise, then train a neural network to walk that corruption backward, one timestep at a time. There is nothing in that recipe that obviously screams "Bayesian statistics." And yet, as we will see over the course of this post, a DDPM is not a new idea bolted onto deep learning — it is the direct, almost inevitable descendant of a single 18th-century theorem and a 21st-century trick for turning intractable integrals into optimization problems.

The trick is called Variational Inference (VI), and its story is older and more general than generative image models. VI was born out of a very mundane frustration: Bayes' theorem tells you exactly how to update your beliefs given data, but the formula requires computing an integral that, for any model with more than a handful of latent variables, simply cannot be evaluated. Statisticians spent decades building sampling-based workarounds (Markov Chain Monte Carlo) before realizing that the posterior distribution could instead be approximated by turning the whole problem into an optimization — searching over a family of simple, tractable distributions for the one member that best resembles the true, intractable posterior.

That single idea — replace an impossible integral with a solvable optimization — is the seed from which an entire lineage of generative models grows. Wrap the approximate posterior in an encoder network and you get a Variational Autoencoder (VAE). Stack several layers of latent variables and you get a Hierarchical VAE (HVAE). Freeze the encoder into a fixed, hand-designed noise schedule instead of learning it, and — almost by accident — you get a Denoising Diffusion Probabilistic Model. Push the number of noising steps to infinity and the whole discrete construction melts into a Stochastic Differential Equation, quietly revealing that the noise-prediction network trained by a simple mean-squared-error loss was secretly estimating a score function the entire time.

This post walks that exact path, stop by stop, keeping the mathematics honest at every step rather than skipping straight to the familiar $L_{simple}$ loss. The goal is not merely to state that "DDPM is a special case of HVAE" — a sentence that appears, unproven, in a great many blog posts — but to actually show why it is true, and to show exactly which piece of freedom Ho et al. (2020) removed from the general HVAE picture to make the mathematics collapse so beautifully.

Below is a summary of the upcoming sections:

• Section 2 - First stop: We state the Bayesian inference problem precisely, show exactly where it becomes intractable, and derive the Evidence Lower Bound (ELBO) from scratch — the single inequality that everything else in this post is built on top of.
• Section 3 - Second stop: We parameterize the variational family with neural networks, arriving at the Variational Autoencoder, and diagnose precisely why a single Gaussian bottleneck tends to produce blurry reconstructions.
• Section 4 - The transformation: We stack the single latent variable into a chain of latents, obtaining the Hierarchical VAE, and show exactly where its ELBO refuses to simplify — a tension that sets up everything that follows.
• Section 5 - The Leap Forward: We show that DDPM resolves that tension by making one dramatic simplification: the "encoder" direction of the hierarchy is no longer trained at all, but fixed to a known linear-Gaussian noise schedule.
• Section 6 - The Elegance of Mathematics: We use the fixed encoder from Section 5 to derive a closed-form posterior via a Bayes'-rule conditioning trick, and watch the resulting KL divergence collapse into the simple mean-squared-error loss between true and predicted noise.
• Section 7 - From discrete to continuous: We take the number of diffusion steps to infinity, arriving at the Stochastic Differential Equation formulation, and show that the noise-prediction network was, all along, secretly a score estimator.
• Section 8 - Conclusion: We look back at the full journey and draw out the single recurring principle that ties Bayes' theorem to modern generative AI together.

2. First stop: The pain of Bayes' Theorem and the Variational Inference Solution

2.1 The Bayesian Inference Problem

Every Bayesian model begins with the same question. Given a set of observed data $\mathbf{x}$ and a set of latent (unobserved) variables $\mathbf{z}$, what is the posterior distribution $p(\mathbf{z}|\mathbf{x})$ — that is, what do we believe about $\mathbf{z}$ once we have seen $\mathbf{x}$? Bayes' theorem answers this in one line: $$p(\mathbf{z}|\mathbf{x}) = \frac{p(\mathbf{x},\mathbf{z})}{p(\mathbf{x})} = \frac{p(\mathbf{x}|\mathbf{z})p(\mathbf{z})}{p(\mathbf{x})}$$ Both the likelihood $p(\mathbf{x}|\mathbf{z})$ and the prior $p(\mathbf{z})$ are, in most models of interest, things we can write down and evaluate directly. The entire difficulty of Bayesian inference is hiding inside the denominator.

2.2 The Intractable Evidence

The denominator $p(\mathbf{x})$, called the evidence or the marginal likelihood, is obtained by integrating the joint distribution over every possible configuration of the latent variables: $$p(\mathbf{x}) = \int_{\mathbb{R}^m} p(\mathbf{x},\mathbf{z})\,d\mathbf{z}$$ For a model with $m$ latent dimensions, this is an $m$-dimensional integral with, in general, no closed form. Even a modest latent space renders it numerically infeasible — this is precisely the "curse of dimensionality" that made exact Bayesian inference impractical for anything beyond the simplest conjugate models for most of the 20th century, and it is exactly as intractable whether the latent variable is a handful of cluster assignments or a 3-channel image's worth of Gaussian noise.

2.3 Variational Inference: From Inference to Optimization

Variational Inference sidesteps the integral entirely, rather than trying to compute it more cleverly. Instead of insisting on the exact posterior, we pick a family of tractable distributions $\mathcal{Q}$ — say, all Gaussians, or all distributions that factorize across dimensions — and search within that family for the member $q^*(\mathbf{z})$ that lies closest to the true posterior, where "closest" is measured by the Kullback-Leibler (KL) divergence:

This is the philosophical pivot the entire rest of this post rests on: a problem of statistical inference — "what is the true posterior?" — has just been rewritten as a problem of numerical optimization — "which member of an easy family minimizes this divergence?" Optimization problems can be attacked with gradient descent; intractable integrals cannot.

2.4 Deriving the Evidence Lower Bound

There is an obvious objection: the objective in the box above still contains $p(\mathbf{z}|\mathbf{x})$, the very quantity we cannot compute. A short algebraic rearrangement dissolves this problem entirely. Expanding the KL divergence and substituting Bayes' theorem, $$\mathrm{KL}(q(\mathbf{z})\|p(\mathbf{z}|\mathbf{x})) = \mathbb{E}_q[\log q(\mathbf{z})] - \mathbb{E}_q[\log p(\mathbf{z}|\mathbf{x})] = \mathbb{E}_q[\log q(\mathbf{z})] - \mathbb{E}_q\left[\log\frac{p(\mathbf{x},\mathbf{z})}{p(\mathbf{x})}\right]$$ $$= \mathbb{E}_q[\log q(\mathbf{z})] - \mathbb{E}_q[\log p(\mathbf{x},\mathbf{z})] + \log p(\mathbf{x})$$ The last step uses the fact that $\log p(\mathbf{x})$ does not depend on $\mathbf{z}$, so $\mathbb{E}_q[\log p(\mathbf{x})] = \log p(\mathbf{x})\int q(\mathbf{z})\,d\mathbf{z} = \log p(\mathbf{x})$. This single line is the whole trick: the untouchable evidence $p(\mathbf{x})$ has been isolated as an additive constant that does not depend on $q$, so it plays no role in the minimization over $\mathcal{Q}$. Define the Evidence Lower Bound (ELBO) as the remaining, fully tractable piece:

Notice what has happened: minimizing the KL divergence to the intractable posterior is now exactly equivalent to maximizing a quantity, $\mathrm{ELBO}(q)$, that requires only the joint distribution $p(\mathbf{x},\mathbf{z})$ — something we always know how to write down — and the variational distribution $q(\mathbf{z})$, which we chose to be tractable in the first place. The evidence $p(\mathbf{x})$ never needs to be computed. This is the engine underneath every model discussed in this post; only the choice of $\mathcal{Q}$ and the parameterization of $q$ will change from here on.

The classical choice of $\mathcal{Q}$ is the mean-field family, where the latent variables are assumed mutually independent, $q(\mathbf{z}) = \prod_j q_j(z_j)$, which yields the Coordinate Ascent Variational Inference (CAVI) algorithm (Blei et al., 2017). We will not need CAVI's coordinate updates in this post — deep generative models replace mean-field factorization with neural-network parameterizations instead — but it is worth keeping in mind that the ELBO derived above is completely agnostic to how $\mathcal{Q}$ is chosen. That freedom is exactly what the rest of this post exploits.

3. Second stop: Variational Autoencoder (VAE) and Bottleneck

3.1 Parameterizing the ELBO with Neural Networks

The mean-field family is a reasonable choice when $q_j(z_j)$ can be updated in closed form, as in the Gaussian-mixture and LDA examples worked out in the companion technical note. It becomes a poor choice the moment the likelihood $p(\mathbf{x}|\mathbf{z})$ is itself a deep neural network with no conjugate structure to exploit. Kingma and Welling (2013) proposed a different route: instead of restricting $\mathcal{Q}$ by factorization, restrict it by parameterization. Let an encoder network with parameters $\theta$ output the parameters of a Gaussian, $$q_\theta(\mathbf{z}|\mathbf{x}) = \mathcal{N}(\mathbf{z};\boldsymbol\mu_\theta(\mathbf{x}), \boldsymbol\sigma_\theta^2(\mathbf{x})\mathbf{I})$$ and let a decoder network with parameters $\phi$ define the likelihood $p_\phi(\mathbf{x}|\mathbf{z})$. Plugging these into the ELBO of Section 2.4 and using the standard Gaussian prior $p(\mathbf{z}) = \mathcal{N}(0,\mathbf{I})$ gives the familiar Variational Autoencoder objective: $$\mathcal{L}_{VAE}(\theta,\phi;\mathbf{x}) = \underbrace{\mathbb{E}_{q_\theta(\mathbf{z}|\mathbf{x})}\big[\log p_\phi(\mathbf{x}|\mathbf{z})\big]}_{\text{reconstruction}} - \underbrace{\mathrm{KL}\big(q_\theta(\mathbf{z}|\mathbf{x})\,\|\,p(\mathbf{z})\big)}_{\text{regularization}}$$ The first term rewards the decoder for reconstructing $\mathbf{x}$ from a sampled $\mathbf{z}$; the second penalizes the encoder for straying too far from the simple, samplable prior $\mathcal{N}(0,\mathbf{I})$ — this is exactly what makes the model generative, since we can later sample $\mathbf{z}\sim\mathcal{N}(0,\mathbf{I})$ directly and skip the encoder altogether.

3.2 The Reparameterization Trick

Maximizing $\mathcal{L}_{VAE}$ by gradient descent requires differentiating through a sampling operation, $\mathbf{z}\sim q_\theta(\mathbf{z}|\mathbf{x})$, which is not directly differentiable. The reparameterization trick moves the randomness outside the computation graph: $$\mathbf{z} = \boldsymbol\mu_\theta(\mathbf{x}) + \boldsymbol\sigma_\theta(\mathbf{x})\odot\boldsymbol\epsilon, \qquad \boldsymbol\epsilon\sim\mathcal{N}(0,\mathbf{I})$$ Now $\mathbf{z}$ is a deterministic, differentiable function of $\theta$ and an independent noise source $\boldsymbol\epsilon$, so ordinary backpropagation applies to both the encoder and the decoder simultaneously — end-to-end training of an ELBO, at last, in a single gradient step.

3.3 Why VAEs Blur

VAE samples are notorious for looking soft and blurry rather than crisp. The reason is visible directly in the reconstruction term. If the decoder likelihood is modeled as an isotropic Gaussian, $p_\phi(\mathbf{x}|\mathbf{z}) = \mathcal{N}(\mathbf{x}; D_\phi(\mathbf{z}), \mathbf{I})$, then maximizing $\mathbb{E}_q[\log p_\phi(\mathbf{x}|\mathbf{z})]$ is equivalent to minimizing $\mathbb{E}_q\|\mathbf{x} - D_\phi(\mathbf{z})\|^2$ — an ordinary mean-squared-error regression loss. But many different plausible images $\mathbf{x}$ can be compatible with the same compressed code $\mathbf{z}$ (a bottleneck that has thrown away high-frequency detail cannot distinguish between them), and the squared-error-minimizing prediction for a multimodal target is not any one of those plausible images — it is their average, $\mathbb{E}[\mathbf{x}|\mathbf{z}]$. Averaging over many plausible sharp images produces a single blurry one. This single-step, single-bottleneck limitation is exactly what the next section tries to fix.

4. The Transformation: Hierarchical VAE (HVAE)

4.1 Stacking the Bottleneck

If one Gaussian bottleneck loses too much information, a natural remedy is to use several, arranged as a chain rather than a single layer. This is the Hierarchical VAE (HVAE) (Sønderby et al., 2016). From here on we index the latents to match the notation used throughout the rest of this post — and, not coincidentally, throughout Diffusion Models — writing $\mathbf{x}_0 := \mathbf{x}$ for the data itself and $\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_T$ for a chain of $T$ increasingly abstract latent layers.

The encoder (inference model) runs bottom-up, from data toward abstraction, and is assumed Markov — each layer depends only on the one below it: $$q_\phi(\mathbf{x}_{1:T}|\mathbf{x}_0) = \prod_{t=1}^{T} q_\phi(\mathbf{x}_t|\mathbf{x}_{t-1})$$ The decoder (generative model) runs top-down, from abstraction back to data, and is likewise Markov: $$p_\theta(\mathbf{x}_{0:T}) = p(\mathbf{x}_T)\prod_{t=1}^{T} p_\theta(\mathbf{x}_{t-1}|\mathbf{x}_t)$$ Both chains are, in a standard HVAE, fully learned neural networks — one for every layer of the encoder, one for every layer of the decoder.

4.2 The Hierarchical ELBO

Plugging these two factorizations into the definition of the ELBO from Section 2.4 and expanding the logarithms of the products into sums, $$\log p_\theta(\mathbf{x}_{0:T}) - \log q_\phi(\mathbf{x}_{1:T}|\mathbf{x}_0) = \log p(\mathbf{x}_T) + \sum_{t=1}^{T}\log p_\theta(\mathbf{x}_{t-1}|\mathbf{x}_t) - \sum_{t=1}^{T}\log q_\phi(\mathbf{x}_t|\mathbf{x}_{t-1})$$ Peel the $t=1$ term off each sum — it is the only term touching the observed data $\mathbf{x}_0$ directly — and the rest reorganizes as: $$\mathrm{ELBO} = \mathbb{E}_q\Big[\underbrace{\log p_\theta(\mathbf{x}_0|\mathbf{x}_1) - \log q_\phi(\mathbf{x}_1|\mathbf{x}_0)}_{\text{data layer}} + \log p(\mathbf{x}_T) + \sum_{t=2}^{T}\Big(\log p_\theta(\mathbf{x}_{t-1}|\mathbf{x}_t) - \log q_\phi(\mathbf{x}_t|\mathbf{x}_{t-1})\Big)\Big]$$ Here is the tension this section wants the reader to notice. In the plain single-layer VAE ELBO of Section 3.1, the reconstruction term and the KL term were clean and separate because both $q_\theta(\mathbf{z}|\mathbf{x})$ and the target prior $p(\mathbf{z})$ shared the same argument, $\mathbf{z}$. Here, no such luck: inside the sum, $p_\theta(\mathbf{x}_{t-1}|\mathbf{x}_t)$ is a distribution over $\mathbf{x}_{t-1}$, while $q_\phi(\mathbf{x}_t|\mathbf{x}_{t-1})$ is a distribution over $\mathbf{x}_t$ — different random variables. They cannot be paired into a clean $\mathrm{KL}(q\|p)$ term the way the single-layer case could. In practice, this forces real HVAE architectures (such as Ladder VAEs) to build top-down inference networks, where $q_\phi(\mathbf{x}_t|\mathbf{x}_{t-1})$ is additionally corrected using information flowing down from the generative path, purely so that the two sides of each layer line up into a matched, tractable KL divergence.

There is, however, a more drastic fix available — one that does not patch the inference network at all, but instead removes its need to be learned in the first place. That fix is exactly what turns an HVAE into a diffusion model.

5. The Leap Forward: DDPM - A Special Form of HVAE

5.1 A Fixed, Untrained Encoder

Denoising Diffusion Probabilistic Models (Ho et al., 2020) make a single, decisive choice: the encoder chain $q(\mathbf{x}_t|\mathbf{x}_{t-1})$ from Section 4.1 is no longer a trainable network at all. It is fixed in advance to a simple linear-Gaussian Markov chain governed by a hand-picked noise schedule $\{\beta_t\}_{t=1}^{T}$: $$q(\mathbf{x}_t|\mathbf{x}_{t-1}) := \mathcal{N}\big(\mathbf{x}_t;\sqrt{1-\beta_t}\,\mathbf{x}_{t-1},\,\beta_t\mathbf{I}\big)$$ There is no $\phi$ anywhere in this expression — nothing here is learned. This chain has no encoder network to train, no amortization gap to worry about, and, because it is linear-Gaussian, it composes: applying the reparameterization trick recursively across all $t$ steps gives a closed-form expression for jumping directly from $\mathbf{x}_0$ to any $\mathbf{x}_t$ in one shot, with $\alpha_t := 1-\beta_t$ and $\bar\alpha_t := \prod_{s=1}^{t}\alpha_s$: $$q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}\big(\mathbf{x}_t;\sqrt{\bar\alpha_t}\,\mathbf{x}_0,\,(1-\bar\alpha_t)\mathbf{I}\big), \qquad \mathbf{x}_t = \sqrt{\bar\alpha_t}\,\mathbf{x}_0 + \sqrt{1-\bar\alpha_t}\,\boldsymbol\epsilon$$ This closed form is precisely what will let us resolve the mismatched-argument problem left open at the end of Section 4.

Notice also a second, quieter difference from the classical HVAE picture: every $\mathbf{x}_t$ lives in the same dimensional space as $\mathbf{x}_0$. There is no genuine compression happening as $t$ grows — only the gradual replacement of signal with noise. It would be more accurate to call this a "noise-neck" than a bottleneck.

5.2 The Learned Decoder

Everything that is learned lives in the reverse, generative direction, exactly as in a standard HVAE: $$p_\theta(\mathbf{x}_{t-1}|\mathbf{x}_t) = \mathcal{N}\big(\mathbf{x}_{t-1};\boldsymbol\mu_\theta(\mathbf{x}_t,t),\Sigma_\theta(\mathbf{x}_t,t)\big)$$ A single neural network, shared across all $T$ layers and conditioned on a timestep embedding $t$, plays the role that would otherwise have required $T$ separate decoder networks — a practical simplification on top of the theoretical one.

5.3 DDPM as a Constrained HVAE

This single restriction — fixing rather than learning the encoder — is the entire conceptual leap from Section 4 to Section 5. Its mathematical payoff, worked out next, is enormous: because $q$ is now a known, explicit Gaussian rather than an opaque neural network, we can finally compute $q(\mathbf{x}_{t-1}|\mathbf{x}_t,\mathbf{x}_0)$ in closed form, which is exactly the missing ingredient needed to pair every term of the HVAE ELBO into a clean, matched KL divergence.

6. The Elegance of Mathematics: Transforming ELBO into a Simple MSE Function

6.1 The Full Diffusion ELBO

Take the negative of the hierarchical ELBO from Section 4.2 and specialize it to DDPM's fixed forward process. Writing the loss as $L := -\mathrm{ELBO}$, and temporarily leaving the $t=1$ and $t=T$ boundary terms aside, we are left with the raw, unresolved sum from Section 4.2: $$L = \Big[-\log p_\theta(\mathbf{x}_0|\mathbf{x}_1) + \log q(\mathbf{x}_1|\mathbf{x}_0)\Big] - \log p(\mathbf{x}_T) + \sum_{t=2}^{T}\Big[\log q(\mathbf{x}_t|\mathbf{x}_{t-1}) - \log p_\theta(\mathbf{x}_{t-1}|\mathbf{x}_t)\Big]$$ Every term inside the sum still has the mismatched-argument problem flagged in Section 4.2: $q(\mathbf{x}_t|\mathbf{x}_{t-1})$ is over $\mathbf{x}_t$, while $p_\theta(\mathbf{x}_{t-1}|\mathbf{x}_t)$ is over $\mathbf{x}_{t-1}$. This is exactly where the fixed encoder from Section 5.1 earns its keep.

6.2 The Conditioning Trick

Because the forward chain is Markov, conditioning it additionally on $\mathbf{x}_0$ changes nothing once $\mathbf{x}_{t-1}$ is already known: $q(\mathbf{x}_t|\mathbf{x}_{t-1}) = q(\mathbf{x}_t|\mathbf{x}_{t-1},\mathbf{x}_0)$. This lets us apply Bayes' rule to flip the direction of conditioning, using the closed forms $q(\mathbf{x}_t|\mathbf{x}_0)$ and $q(\mathbf{x}_{t-1}|\mathbf{x}_0)$ derived in Section 5.1: $$q(\mathbf{x}_t|\mathbf{x}_{t-1}) = q(\mathbf{x}_t|\mathbf{x}_{t-1},\mathbf{x}_0) = \frac{q(\mathbf{x}_{t-1}|\mathbf{x}_t,\mathbf{x}_0)\,q(\mathbf{x}_t|\mathbf{x}_0)}{q(\mathbf{x}_{t-1}|\mathbf{x}_0)}$$ Substituting this into the sum from Section 6.1, the ratio $q(\mathbf{x}_t|\mathbf{x}_0)/q(\mathbf{x}_{t-1}|\mathbf{x}_0)$ telescopes across $t=2,\dots,T$ — every intermediate value of $q(\mathbf{x}_t|\mathbf{x}_0)$ cancels against the next term's $q(\mathbf{x}_{t-1}|\mathbf{x}_0)$, leaving only the two endpoints, $q(\mathbf{x}_T|\mathbf{x}_0)$ and $q(\mathbf{x}_1|\mathbf{x}_0)$. The latter exactly cancels the boundary term we set aside in Section 6.1. What remains reorganizes, after this bookkeeping, into a clean sum of matched-argument KL divergences: $$L = \underbrace{\mathrm{KL}\big(q(\mathbf{x}_T|\mathbf{x}_0)\,\|\,p(\mathbf{x}_T)\big)}_{L_T,\ \approx 0\text{ by design}} + \sum_{t=2}^{T}\underbrace{\mathrm{KL}\big(q(\mathbf{x}_{t-1}|\mathbf{x}_t,\mathbf{x}_0)\,\|\,p_\theta(\mathbf{x}_{t-1}|\mathbf{x}_t)\big)}_{L_{t-1}} \underbrace{- \log p_\theta(\mathbf{x}_0|\mathbf{x}_1)}_{L_0}$$ Both sides of every $L_{t-1}$ term are now distributions over the same random variable, $\mathbf{x}_{t-1}$ — precisely the resolution promised at the end of Section 5.3. $L_T$ has no trainable parameters (the schedule is fixed so that $\bar\alpha_T\approx 0$, making $q(\mathbf{x}_T|\mathbf{x}_0)$ almost exactly the prior $\mathcal{N}(0,\mathbf{I})$), so training reduces to minimizing the $T-1$ terms $L_{t-1}$.

6.3 A Closed-Form Posterior

We still need $q(\mathbf{x}_{t-1}|\mathbf{x}_t,\mathbf{x}_0)$ explicitly. Because it is a ratio of Gaussian densities (Bayes' rule above), its log is quadratic in $\mathbf{x}_{t-1}$, so the distribution itself must be Gaussian — the same completing-the-square argument used throughout the mean-field derivations in Section 2.

This is exactly what "special form of HVAE" was buying us: a fully explicit, closed-form target distribution for $\mathbf{x}_{t-1}$ at every single layer, obtained without ever training an encoder — impossible in a general HVAE, where $q_\phi(\mathbf{x}_t|\mathbf{x}_{t-1})$ is an opaque neural network with no such conjugate structure.

6.4 From Mean-Matching to Noise-Matching

Set the decoder's covariance to match the fixed target, $\Sigma_\theta(\mathbf{x}_t,t) := \tilde\beta_t\mathbf{I}$. The KL divergence between two Gaussians that share an identical, non-trainable covariance collapses to a plain squared distance between their means (the covariance terms in the KL formula cancel exactly): $$L_{t-1} = \mathbb{E}_q\left[\frac{1}{2\tilde\beta_t}\big\|\tilde{\boldsymbol\mu}_t(\mathbf{x}_t,\mathbf{x}_0) - \boldsymbol\mu_\theta(\mathbf{x}_t,t)\big\|^2\right] + C$$ Now substitute the forward-process reparameterization $\mathbf{x}_t = \sqrt{\bar\alpha_t}\mathbf{x}_0 + \sqrt{1-\bar\alpha_t}\boldsymbol\epsilon$ (so $\mathbf{x}_0 = (\mathbf{x}_t - \sqrt{1-\bar\alpha_t}\boldsymbol\epsilon)/\sqrt{\bar\alpha_t}$) into $\tilde{\boldsymbol\mu}_t$ from Section 6.3. After the algebra simplifies, the target mean depends on $\mathbf{x}_t$ and the specific noise $\boldsymbol\epsilon$ that produced it — nothing else: $$\tilde{\boldsymbol\mu}_t(\mathbf{x}_t,\mathbf{x}_0) = \frac{1}{\sqrt{\alpha_t}}\left(\mathbf{x}_t - \frac{\beta_t}{\sqrt{1-\bar\alpha_t}}\boldsymbol\epsilon\right)$$ The natural move is to give $\boldsymbol\mu_\theta$ the exact same functional form, but with $\boldsymbol\epsilon$ replaced by a network's prediction $\boldsymbol\epsilon_\theta(\mathbf{x}_t,t)$: $$\boldsymbol\mu_\theta(\mathbf{x}_t,t) = \frac{1}{\sqrt{\alpha_t}}\left(\mathbf{x}_t - \frac{\beta_t}{\sqrt{1-\bar\alpha_t}}\boldsymbol\epsilon_\theta(\mathbf{x}_t,t)\right)$$ Substituting both expressions back into $L_{t-1}$, every $\mathbf{x}_t$-dependent term cancels identically, leaving a plain weighted noise-matching loss: $$L_{t-1} = \mathbb{E}_{\mathbf{x}_0,\boldsymbol\epsilon}\left[\frac{\beta_t^2}{2\tilde\beta_t\,\alpha_t(1-\bar\alpha_t)}\big\|\boldsymbol\epsilon - \boldsymbol\epsilon_\theta(\mathbf{x}_t,t)\big\|^2\right]$$ Ho et al. (2020) observed empirically that dropping this timestep-dependent weighting entirely — treating every timestep equally — trains more stably and produces higher sample quality. The result is the celebrated simplified training objective, which the second post in this series introduced without derivation: $$L_{simple}(\theta) = \mathbb{E}_{t\sim\mathcal{U}\{1,\dots,T\},\;\mathbf{x}_0,\;\boldsymbol\epsilon\sim\mathcal{N}(0,\mathbf{I})}\Big[\big\|\boldsymbol\epsilon - \boldsymbol\epsilon_\theta(\mathbf{x}_t,t)\big\|^2\Big]$$ Trace the chain backward and the entire, seemingly ad hoc "predict-the-noise" recipe is revealed as nothing more than the Evidence Lower Bound of Section 2.4, specialized layer by layer through Sections 4 and 5, and simplified by a single conditioning trick. Bayes' theorem, an intractable integral, and a mean-squared-error loss turn out to be three views of the exact same object.

7. From discrete to continuous: Stochastic Differential Equations (SDEs)

7.1 The Continuum Limit

Everything so far has treated $T$ as a large but finite number of discrete layers. Song et al. (2021) showed that letting $T\to\infty$, while shrinking each individual noise increment $\beta_t\to 0$ accordingly, turns the discrete Markov chain of Section 5.1 into a continuous-time process. Rescale the discrete index to a continuous time variable $t\in[0,1]$ and define a noise-rate function $\beta(t)$ such that $\beta_t \approx \beta(t)\,\Delta t$ for step size $\Delta t = 1/T$. The forward update $\mathbf{x}_t = \sqrt{1-\beta_t}\,\mathbf{x}_{t-1} + \sqrt{\beta_t}\,\boldsymbol\epsilon_t$ then converges, in this limit, to the Variance-Preserving Stochastic Differential Equation: $$d\mathbf{x} = -\frac{1}{2}\beta(t)\,\mathbf{x}\,dt + \sqrt{\beta(t)}\,d\mathbf{w}$$ where $\mathbf{w}$ is a standard Wiener process. Every one of the $T$ separate Gaussian transition kernels from Section 5.1 has been absorbed into a single drift-and-diffusion pair, $\mathbf{f}(\mathbf{x},t) = -\frac{1}{2}\beta(t)\mathbf{x}$ and $g(t) = \sqrt{\beta(t)}$ — the exact SDE form used throughout "Advanced Diffusion Models".

7.2 Forward and Reverse-Time SDEs

A classical result of Anderson (1982) states that any such forward SDE has a corresponding reverse-time SDE with the same marginal densities $p_t(\mathbf{x})$ at every $t$, running backward from $t=1$ to $t=0$: $$d\mathbf{x} = \Big[-\frac{1}{2}\beta(t)\,\mathbf{x} - \beta(t)\,\nabla_{\mathbf{x}}\log p_t(\mathbf{x})\Big]dt + \sqrt{\beta(t)}\,d\bar{\mathbf{w}}$$ Everything about this reverse process is known and computable — except for a single missing ingredient, the score function $\nabla_{\mathbf{x}}\log p_t(\mathbf{x})$. Simulating this equation backward from pure noise at $t=1$ down to $t=0$, if the score were known exactly, would produce a genuine sample from $p_{data}$.

7.3 The Score Was There All Along

This is where the entire journey folds back on itself. Because $q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t;\sqrt{\bar\alpha_t}\mathbf{x}_0,(1-\bar\alpha_t)\mathbf{I})$ is a known Gaussian (Section 5.1), its own score has a closed form, $\nabla_{\mathbf{x}_t}\log q(\mathbf{x}_t|\mathbf{x}_0) = -\boldsymbol\epsilon/\sqrt{1-\bar\alpha_t}$. Denoising score matching shows that training a network to recover this conditional score, averaged over the data distribution, recovers the true marginal score $\nabla_{\mathbf{x}_t}\log p_t(\mathbf{x}_t)$ as well — and so: $$\nabla_{\mathbf{x}_t}\log p_t(\mathbf{x}_t) \approx -\frac{\boldsymbol\epsilon_\theta(\mathbf{x}_t,t)}{\sqrt{1-\bar\alpha_t}}$$ The exact same $\boldsymbol\epsilon_\theta$ trained by the plain mean-squared-error loss $L_{simple}$ derived purely from Variational Inference in Section 6.4 is, up to this fixed rescaling, precisely the score estimator that the reverse-time SDE of Section 7.2 needs. Nobody set out to train a score function — it fell out of maximizing an Evidence Lower Bound. This is the third of the "triple perspective" on diffusion models introduced informally in "Overview of Diffusion Models" (nonequilibrium thermodynamics, hierarchical variational inference, and score-based SDEs) meeting concretely, in one equation, for the first time.

8. Conclusion

We started this post with the single most mundane complaint in Bayesian statistics: the evidence $p(\mathbf{x})$ is an intractable integral. Variational Inference turned that complaint into an optimization problem via the ELBO. Wrapping the ELBO's variational family in an encoder network gave a VAE; a Gaussian bottleneck that was too narrow gave way to a stack of bottlenecks, the HVAE; and a single, decisive simplification — freezing that stack's encoder into a fixed noise schedule instead of learning it — gave DDPM. A conditioning trick made possible only by that frozen encoder collapsed a tangle of mismatched-argument KL divergences into the simple, almost embarrassingly ordinary mean-squared-error loss between real and predicted noise. And taking the number of noising steps to infinity dissolved the whole discrete machinery into a Stochastic Differential Equation, revealing that this simple regression loss had, from its very first derivation, been secretly estimating a score function.

The recurring lesson is not any one equation, but the shape of the whole journey: every major advance covered here — VAE, HVAE, DDPM, and the SDE limit — is the exact same optimization principle from Section 2.4, applied again and again with a progressively cleverer choice of variational family $\mathcal{Q}$. Nothing about the neural network architecture changed the underlying mathematics; what changed, every time, was which constraints we were willing to place on $q$. From here, the natural next step is the one taken in "Advanced Diffusion Models": having arrived at the score function and the SDE, one can ask whether the stochastic path from noise to data was ever necessary at all — and discover Flow Matching, Optimal Transport, and Schrödinger Bridges waiting on the other side of that same question.