where $\bar{\alpha}_k = \prod_{i=1}^k (1 - \beta_i)$. The noise schedule ensures $\bar{\alpha}_k \to 0$ as $k \to N$, so $\mathbf{x}_N$ converges to $\mathcal{N}(\mathbf{0}, \mathbf{I})$. The forward process induces a family of conditional probability densities $q(\mathbf{x}_t|\mathbf{x}_0)$ describing the distribution of $\mathbf{x}_t$ given the initial data $\mathbf{x}_0$. For the VP-SDE, the perturbation kernel is:

$$q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \alpha_t \mathbf{x}_0, \sigma_t^2 \mathbf{I}), \qquad (24.23)$$

where:

$$\alpha_t = e^{-\frac{1}{2} \int_0^t \beta(s) ds}, \quad \sigma_t^2 = 1 - e^{-\int_0^t \beta(s) ds}. \qquad (24.24)$$

For the VE-SDE, the kernel is:

$$q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \mathbf{x}_0, \sigma^2(t)\mathbf{I}), \qquad (24.25)$$

where $\sigma(t)$ is a monotonically increasing function. The evolution of the marginal density $q_t(\mathbf{x})$ under the forward SDE is described by the Fokker-Planck equation:

$$\frac{\partial q_t(\mathbf{x})}{\partial t} = -\nabla \cdot [\mathbf{f}(\mathbf{x}, t)q_t(\mathbf{x})] + \frac{1}{2}g^2(t)\nabla^2 q_t(\mathbf{x}), \qquad (24.26)$$

where $\nabla \cdot$ is the divergence operator and $\nabla^2$ is the Laplacian. This PDE governs how the data distribution $q_t(\mathbf{x})$ smoothens into a Gaussian. The infinitesimal generator $\mathcal{L}_t$ of the forward diffusion process is:

$$\mathcal{L}_t = \mathbf{f}(\mathbf{x}, t) \cdot \nabla + \frac{1}{2} g^2(t) \nabla^2, \qquad (24.27)$$

which acts on test functions $\phi(\mathbf{x})$ via:

$$\mathcal{L}_t \phi(\mathbf{x}) = \lim_{\Delta t \to 0} \frac{\mathbb{E}[\phi(\mathbf{x}_{t+\Delta t})|\mathbf{x}_t = \mathbf{x}] - \phi(\mathbf{x})}{\Delta t}. \qquad (24.28)$$

The forward process is designed so its time reversal (learned in the generative phase) admits a tractable SDE. By the Anderson-Chen-So theorem, the reverse SDE is:

$$d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g^2(t)\nabla_{\mathbf{x}} \log q_t(\mathbf{x})] dt + g(t) d\bar{\mathbf{w}}, \qquad (24.29)$$

where $\bar{\mathbf{w}}$ is a reverse-time Wiener process. This necessitates learning the score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$. In score-based models, the forward process is implicit, defined via a sequence of noise scales $\{\sigma_i\}_{i=1}^L$. The perturbed distribution at scale $\sigma_i$ is:

$$p_{\sigma_i}(\mathbf{x}) = \int p_{\text{data}}(\mathbf{y}) \mathcal{N}(\mathbf{x}; \mathbf{y}, \sigma_i^2 \mathbf{I}) d\mathbf{y}, \qquad (24.30)$$

and the score $\nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x})$ is learned directly. The forward process is thus a family of increasingly noisy distributions without an explicit SDE. There is a thermodynamic interpretation of entropy increase of the forward diffusion process. The forward diffusion can be viewed as an entropy-increasing process, where the system evolves toward maximum entropy (a Gaussian). The rate of entropy production is:

$$\frac{dS(t)}{dt} = \frac{1}{2} g^2(t) \mathbb{E}_{q_t(\mathbf{x})} [\|\nabla_{\mathbf{x}} \log q_t(\mathbf{x})\|^2], \qquad (24.31)$$

where $S(t)$ is the differential entropy of $q_t(\mathbf{x})$. For numerical implementation, the forward SDE is discretized (e.g., Euler-Maruyama):

$$\mathbf{x}_{t+\Delta t} = \mathbf{x}_t + \mathbf{f}(\mathbf{x}_t, t)\Delta t + g(t)\sqrt{\Delta t}\mathbf{z}, \quad \mathbf{z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}). \qquad (24.32)$$

The Summary of Key Equations are:

1. Forward SDE: $d\mathbf{x} = \mathbf{f}(\mathbf{x}, t)dt + g(t)d\mathbf{w}$.
2. **Perturbation Kernel:** $q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \alpha_t \mathbf{x}_0, \sigma_t^2 \mathbf{I})$.

3. **Fokker-Planck Equation:** $\frac{\partial q_t}{\partial t} = -\nabla \cdot [\mathbf{f}q_t] + \frac{1}{2}g^2\nabla^2 q_t$.

4. **Infinitesimal Generator:** $\mathcal{L}_t = \mathbf{f} \cdot \nabla + \frac{1}{2}g^2\nabla^2$.

5. **Score-Based Perturbation:** $p_{\sigma_i}(\mathbf{x}) = \int p_{\text{data}}(\mathbf{y})\mathcal{N}(\mathbf{x}; \mathbf{y}, \sigma_i^2\mathbf{I})d\mathbf{y}$.

The forward diffusion process is a mathematically rigorous framework for gradually destroying data structure through noise, enabling generative modeling via reversal. Whether formulated as an SDE (diffusion models) or a sequence of noise-perturbed distributions (score-based models), it provides a principled way to transform complex data into a tractable prior distribution. The equations above fully characterize its continuous and discrete dynamics, perturbation kernels, and statistical properties.

### 24.3.2 Reverse Diffusion Process

The reverse diffusion process is the cornerstone of generative modeling in diffusion and score-based frameworks, enabling the synthesis of data from noise by inverting a carefully constructed forward noising process.

The reverse diffusion process is derived from the time reversal of the forward SDE. Given the forward SDE:

$$d\mathbf{x} = \mathbf{f}(\mathbf{x}, t) dt + g(t) d\mathbf{w}, \qquad (24.33)$$

the reverse-time SDE (Anderson, 1982) is:

$$d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x})] dt + g(t) d\bar{\mathbf{w}}, \qquad (24.34)$$

where:

1. $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ is the score function of the perturbed data distribution at time $t$,

2. $\bar{\mathbf{w}}$ is a reverse-time Wiener process,

3. The term $g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ is the drift correction ensuring $q_t(\mathbf{x})$ is preserved.

For the variance-preserving (VP) SDE, with $\mathbf{f}(\mathbf{x}, t) = -\frac{1}{2}\beta(t)\mathbf{x}$ and $g(t) = \sqrt{\beta(t)}$, the reverse SDE becomes:

$$d\mathbf{x} = \left[ -\frac{1}{2}\beta(t)\mathbf{x} - \beta(t)\nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] dt + \sqrt{\beta(t)} d\bar{\mathbf{w}}. \qquad (24.35)$$

The reverse process can alternatively be described by a deterministic ODE (Song et al., 2021), derived by removing the Brownian noise term:

$$d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] dt. \qquad (24.36)$$

This probability flow ODE preserves the same marginal distributions $q_t(\mathbf{x})$ as the reverse SDE but enables more efficient sampling. For the VP-SDE:

$$d\mathbf{x} = \left[ -\frac{1}{2}\beta(t)\mathbf{x} - \frac{1}{2}\beta(t)\nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] dt. \qquad (24.37)$$

The score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ is approximated by a neural network $\mathbf{s}_{\theta}(\mathbf{x}, t)$, trained via denoising score matching (DSM):

$$\mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \mathbf{x}_t} \left[ \lambda(t) \| \mathbf{s}_{\theta}(\mathbf{x}_t, t) - \nabla_{\mathbf{x}_t} \log q(\mathbf{x}_t|\mathbf{x}_0) \|^2 \right], \qquad (24.38)$$

where:
1. $\lambda(t)$ is a weighting function (often $\lambda(t) = g(t)^2$),

2. $q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \alpha_t \mathbf{x}_0, \sigma_t^2 \mathbf{I})$ is the perturbation kernel,

3. $\nabla_{\mathbf{x}_t} \log q(\mathbf{x}_t|\mathbf{x}_0) = -\frac{\mathbf{x}_t - \alpha_t \mathbf{x}_0}{\sigma_t^2}$.

For the VP-SDE, the DSM objective simplifies to:

$$ \mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \epsilon} \left[ \| \mathbf{s}_{\theta}(\alpha_t \mathbf{x}_0 + \sigma_t \epsilon, t) + \epsilon / \sigma_t \|^2 \right], \quad (24.39) $$

where $\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. In discrete-time diffusion models, the reverse process is a Markov chain with learned Gaussian transitions:

$$ p_{\theta}(\mathbf{x}_{k-1}|\mathbf{x}_k) = \mathcal{N}(\mathbf{x}_{k-1}; \mu_{\theta}(\mathbf{x}_k, k), \Sigma_{\theta}(\mathbf{x}_k, k)), \quad (24.40) $$

where $\mu_{\theta}(\mathbf{x}_k, k)$ is the denoising mean, often parameterized as:

$$ \mu_{\theta}(\mathbf{x}_k, k) = \frac{1}{\sqrt{\alpha_k}} \left( \mathbf{x}_k - \frac{\beta_k}{\sqrt{1 - \bar{\alpha}_k}} \mathbf{s}_{\theta}(\mathbf{x}_k, k) \right), \quad (24.41) $$

$\Sigma_{\theta}(\mathbf{x}_k, k)$ is the covariance, typically set to $\tilde{\beta}_k \mathbf{I}$, where $\tilde{\beta}_k = \frac{1-\bar{\alpha}_{k-1}}{1-\bar{\alpha}_k} \beta_k$. In pure score-based models (without a predefined forward SDE), sampling is performed via Langevin dynamics:

$$ \mathbf{x}_{i+1} = \mathbf{x}_i + \epsilon_i \mathbf{s}_{\theta}(\mathbf{x}_i, \sigma_i) + \sqrt{2\epsilon_i} \mathbf{z}_i, \quad \mathbf{z}_i \sim \mathcal{N}(\mathbf{0}, \mathbf{I}), \quad (24.42) $$

where $\epsilon_i$ is the step size and $\sigma_i$ is the noise scale (annealed over time). The reverse process can be interpreted as solving an optimal transport problem, minimizing the KL divergence between the forward and reverse paths. The Schrödinger bridge formulation yields:

$$ \min_{p} \mathbb{K} \mathbb{L} [p \| q], \quad \text{subject to } p(\mathbf{x}_0) = p_{\text{data}}, \quad p(\mathbf{x}_T) = p_{\text{prior}} \quad (24.43) $$

where $q$ is the forward process and $p$ is the reverse process. The reverse process starts from the prior distribution $q_T(\mathbf{x}) = \mathcal{N}(\mathbf{x}; \mathbf{0}, \mathbf{I})$ and is driven by the learned score $\mathbf{s}_{\theta}(\mathbf{x}, t)$. The terminal condition for the reverse SDE is: $\lim_{t \to T} \mathbf{x}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ while the initial condition is: $\lim_{t \to 0} \mathbf{x}_t \sim p_{\text{data}}$. The reverse SDE/ODE is solved numerically. For the reverse SDE (Euler-Maruyama):

$$ \mathbf{x}_{t-\Delta t} = \mathbf{x}_t - \left[ \mathbf{f}(\mathbf{x}_t, t) - g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t) \right] \Delta t + g(t) \sqrt{\Delta t} \mathbf{z}. \quad (24.44) $$

For the ODE (Runge-Kutta methods):

$$ \mathbf{x}_{t-\Delta t} = \mathbf{x}_t - \left[ \mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2} g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t) \right] \Delta t. \quad (24.45) $$

The probability flow ODE enables exact likelihood computation via:

$$ \log p(\mathbf{x}_0) = \log p(\mathbf{x}_T) + \int_0^T \nabla \cdot \left( \mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2} g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t) \right) dt, \quad (24.46) $$

where $\nabla \cdot$ is the divergence operator.

The Summary of Key Equations are:

1. **Reverse SDE:** $dx = [\mathbf{f}(\mathbf{x}, t) - g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x})] dt + g(t) d\bar{\mathbf{w}}$.

2. **Probability Flow ODE:** $dx = [\mathbf{f}(\mathbf{x}, t) - \frac{1}{2} g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x})] dt$.

3. **DSM Objective:** $\mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \epsilon} \left[ \| \mathbf{s}_{\theta}(\alpha_t \mathbf{x}_0 + \sigma_t \epsilon, t) + \epsilon / \sigma_t \|^2 \right]$.
4. Langevin Dynamics: $\mathbf{x}_{i+1} = \mathbf{x}_i + \epsilon_i \mathbf{s}_\theta(\mathbf{x}_i, \sigma_i) + \sqrt{2\epsilon_i} \mathbf{z}_i$.

5. Likelihood Computation: $\log p(\mathbf{x}_0) = \log p(\mathbf{x}_T) + \int_0^T \nabla \cdot (\mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2}g(t)^2 \mathbf{s}_\theta(\mathbf{x}_t, t)) dt$.

The reverse diffusion process is a mathematically rigorous framework for synthesizing data from noise, governed by either a stochastic or deterministic dynamical system. Its core components—reverse SDEs, probability flow ODEs, score matching, and sampling algorithms—are deeply rooted in stochastic calculus, differential equations, and statistical physics. The equations above provide a complete description of its theoretical and computational foundations.

### 24.3.3 Probability Flow ODE

The Probability Flow Ordinary Differential Equation (ODE) is a deterministic process that emerges from the stochastic framework of diffusion models and score-based generative models, providing a tractable alternative to the reverse-time stochastic differential equation (SDE) while preserving the same marginal distributions. The Probability Flow ODE is derived by eliminating the Brownian motion term from the reverse-time SDE, resulting in a deterministic trajectory that can be used for both sampling and exact likelihood computation. Given the reverse-time SDE

$$d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] dt + g(t) d\bar{\mathbf{w}} \quad (24.47)$$

the corresponding Probability Flow ODE is obtained by setting the diffusion coefficient to zero, yielding

$$d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2} g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] dt \quad (24.48)$$

This ODE describes a deterministic flow that evolves the state $\mathbf{x}$ backward in time, starting from the prior distribution $q_T(\mathbf{x}) = \mathcal{N}(\mathbf{x}; 0, \mathbf{I})$ and terminating at the data distribution $q_0(\mathbf{x}) = p_{\text{data}}(\mathbf{x})$.

The drift term in the Probability Flow ODE consists of two components: the deterministic drift $\mathbf{f}(\mathbf{x}, t)$ from the forward SDE and the score-based correction term $-\frac{1}{2}g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x})$. For the variance-preserving (VP) SDE, where $\mathbf{f}(\mathbf{x}, t) = -\frac{1}{2}\beta(t)\mathbf{x}$ and $g(t) = \sqrt{\beta(t)}$, the Probability Flow ODE simplifies to

$$d\mathbf{x} = \left[ -\frac{1}{2}\beta(t)\mathbf{x} - \frac{1}{2}\beta(t)\nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] dt \quad (24.49)$$

The score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ is approximated by a neural network $\mathbf{s}_\theta(\mathbf{x}, t)$, trained via denoising score matching to satisfy $\mathbf{s}_\theta(\mathbf{x}, t) \approx \nabla_{\mathbf{x}} \log q_t(\mathbf{x})$. The training objective for $\mathbf{s}_\theta(\mathbf{x}, t)$ is

$$\mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \mathbf{x}_t} \left[ \lambda(t) \| \mathbf{s}_\theta(\mathbf{x}_t, t) - \nabla_{\mathbf{x}_t} \log q(\mathbf{x}_t | \mathbf{x}_0) \|^2 \right] \quad (24.50)$$

where $\lambda(t)$ is a weighting function and $q(\mathbf{x}_t | \mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \alpha_t \mathbf{x}_0, \sigma_t^2 \mathbf{I})$ is the perturbation kernel. The Probability Flow ODE admits a unique solution for any initial condition $\mathbf{x}_T \sim q_T(\mathbf{x})$, and this solution can be numerically integrated using standard ODE solvers such as the Runge-Kutta methods. The discrete-time update rule for the Euler method is

$$\mathbf{x}_{t-\Delta t} = \mathbf{x}_t - \left[ \mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2}g(t)^2 \mathbf{s}_\theta(\mathbf{x}_t, t) \right] \Delta t \quad (24.51)$$

where $\Delta t$ is the step size. The deterministic nature of the Probability Flow ODE enables exact likelihood computation via the instantaneous change-of-variables formula, which states $\log p(\mathbf{x}_0) = \log p(\mathbf{x}_T) + \int_0^T \nabla \cdot (\mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2}g(t)^2 \mathbf{s}_\theta(\mathbf{x}_t, t)) dt$, where $\nabla \cdot$ denotes the divergence operator.
The Probability Flow ODE is closely related to the Fokker-Planck equation, which describes the evolution of the probability density $q_t(\mathbf{x})$ under the forward SDE. The Fokker-Planck equation is

$$ \frac{\partial q_t(\mathbf{x})}{\partial t} = -\nabla \cdot [\mathbf{f}(\mathbf{x}, t)q_t(\mathbf{x})] + \frac{1}{2}g(t)^2\nabla^2 q_t(\mathbf{x}) \quad (24.52) $$

and the Probability Flow ODE preserves the same marginal densities $q_t(\mathbf{x})$ as the reverse-time SDE. This connection is formalized through the notion of the infinitesimal generator

$$ \mathcal{L}_t = \mathbf{f}(\mathbf{x}, t) \cdot \nabla + \frac{1}{2}g(t)^2\nabla^2 \quad (24.53) $$

which acts on test functions $\phi(\mathbf{x})$ as

$$ \mathcal{L}_t \phi(\mathbf{x}) = \lim_{\Delta t \to 0} \frac{\mathbb{E}[\phi(\mathbf{x}_{t+\Delta t}) | \mathbf{x}_t = \mathbf{x}] - \phi(\mathbf{x})}{\Delta t} \quad (24.54) $$

The deterministic flow induced by the Probability Flow ODE can also be interpreted as a gradient flow in the space of probability measures, minimizing the Kullback-Leibler (KL) divergence between the generated and target distributions. This perspective links the Probability Flow ODE to optimal transport theory, where the drift term $\mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g(t)^2\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ can be seen as a velocity field transporting particles from the prior to the data distribution. The Schrödinger bridge problem, which seeks the most likely path between two distributions, also admits a solution via the Probability Flow ODE, with the drift term adjusted to satisfy the boundary conditions $q_0(\mathbf{x}) = p_{\text{data}}(\mathbf{x})$ and $q_T(\mathbf{x}) = p_{\text{prior}}(\mathbf{x})$. In the context of score-based generative models, the Probability Flow ODE provides a deterministic alternative to Langevin dynamics for sampling. While Langevin dynamics relies on the stochastic update rule

$$ \mathbf{x}_{i+1} = \mathbf{x}_i + \epsilon_i \mathbf{s}_\theta(\mathbf{x}_i, \sigma_i) + \sqrt{2\epsilon_i} \mathbf{z}_i \quad (24.55) $$

The Probability Flow ODE eliminates the noise term $\sqrt{2\epsilon_i} \mathbf{z}_i$ and relies solely on the score function $\mathbf{s}_\theta(\mathbf{x}, t)$ to guide the sampling process. This deterministic approach often yields higher-quality samples with fewer steps, as it avoids the stochastic fluctuations inherent in Langevin dynamics. The Probability Flow ODE also enables efficient computation of marginal likelihoods and conditional probabilities, which are intractable in purely stochastic frameworks. By leveraging the instantaneous change-of-variables formula, one can compute the log-likelihood $\log p(\mathbf{x}_0)$ of a data point $\mathbf{x}_0$ by integrating the divergence of the drift term along the ODE trajectory. This property is particularly useful for tasks such as density estimation and outlier detection, where exact likelihoods are required. The divergence term

$$ \nabla \cdot \left( \mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2} g(t)^2 \mathbf{s}_\theta(\mathbf{x}_t, t) \right) \quad (24.56) $$

can be approximated efficiently using Hutchinson's trace estimator, which reduces the computational cost from $O(d^2)$ to $O(d)$ for $d$-dimensional data.

The stability and convergence properties of the Probability Flow ODE are governed by the Lipschitz continuity of the drift term $\mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g(t)^2\mathbf{s}_\theta(\mathbf{x}, t)$. Under mild assumptions on $\mathbf{f}(\mathbf{x}, t)$ and $\mathbf{s}_\theta(\mathbf{x}, t)$, the ODE admits a unique solution that is continuous with respect to the initial condition $\mathbf{x}_T$. This stability is crucial for numerical integration, as it ensures that small perturbations in the initial state or step size do not lead to divergent trajectories. The choice of ODE solver further influences the accuracy and efficiency of sampling, with higher-order methods such as Runge-Kutta providing better approximations at the cost of increased computational overhead. The Probability Flow ODE also generalizes to non-Gaussian noise distributions and non-linear forward processes, provided the score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ can be accurately estimated. For instance, in the case of a
variance-exploding (VE) SDE, where the forward process is $d\mathbf{x} = \sqrt{\frac{d\sigma^2(t)}{dt}} d\mathbf{w}$, the Probability Flow ODE becomes

$$d\mathbf{x} = -\frac{1}{2} \frac{d\sigma^2(t)}{dt} \nabla_{\mathbf{x}} \log q_t(\mathbf{x}) dt \qquad (24.57)$$

This flexibility allows the Probability Flow ODE to adapt to various noise schedules and perturbation kernels, making it a versatile tool for generative modeling.

The relationship between the Probability Flow ODE and the reverse-time SDE can be further elucidated through the lens of Girsanov's theorem, which provides a measure-theoretic equivalence between the two processes. Specifically, the Radon-Nikodym derivative between the path measures induced by the reverse-time SDE and the Probability Flow ODE is given by

$$\exp\left(\frac{1}{2} \int_{0}^{T} g(t)^2 \|\nabla_{\mathbf{x}} \log q_t(\mathbf{x})\|^2 dt\right) \qquad (24.58)$$

highlighting the role of the score function in bridging the stochastic and deterministic dynamics. This theoretical connection underscores the consistency of the Probability Flow ODE with the underlying diffusion process.

In practice, the Probability Flow ODE is often preferred over the reverse-time SDE for tasks requiring high-fidelity samples or exact likelihoods, due to its deterministic nature and lower variance. However, the choice between ODE and SDE-based sampling depends on the specific application, with SDEs sometimes yielding better sample diversity at the cost of increased computational complexity. Hybrid approaches, which combine deterministic and stochastic steps, have also been explored to balance the trade-offs between sample quality and computational efficiency. The Probability Flow ODE framework has been extended to accommodate conditional generation, where the goal is to sample from a target distribution $p(\mathbf{x}|\mathbf{y})$ given auxiliary information $\mathbf{y}$. In this setting, the score function is replaced by the conditional score $\nabla_{\mathbf{x}} \log q_t(\mathbf{x}|\mathbf{y})$ and the Probability Flow ODE becomes

$$d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2} g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x}|\mathbf{y}) \right] dt \qquad (24.59)$$

This extension enables applications such as image inpainting, super-resolution, and class-conditional generation, where the auxiliary information $\mathbf{y}$ guides the sampling process.

The theoretical foundations of the Probability Flow ODE are deeply rooted in the theory of dynamical systems and partial differential equations. The existence and uniqueness of solutions to the ODE are guaranteed by the Picard-Lindelöf theorem, provided the drift term is Lipschitz continuous in $\mathbf{x}$ and continuous in $t$. The regularity of the score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ plays a critical role in ensuring these conditions are met, and recent work has focused on characterizing the smoothness properties of the score under various assumptions on the data distribution $p_{\text{data}}(\mathbf{x})$. The Probability Flow ODE also admits a variational interpretation, where the drift term

$$\left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2} g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] \qquad (24.60)$$

can be viewed as the gradient of a time-dependent potential function. This perspective connects the ODE to Hamilton-Jacobi-Bellman (HJB) equations in optimal control, with the score function playing the role of the optimal control policy. Such connections have inspired new algorithms for training and sampling from diffusion models, leveraging tools from control theory and variational inference.
### 24.3.4 Training Objective

The training objective in diffusion models and score-based generative models is fundamentally rooted in the principle of score matching, where the goal is to learn a parametric model that approximates the true score function of the data distribution. The score function is defined as the gradient of the log-density with respect to the data, $\nabla_x \log p(\mathbf{x})$, and its estimation is central to both frameworks. In diffusion models, the forward process gradually perturbs the data distribution $p_{\text{data}}(\mathbf{x})$ through a predefined stochastic differential equation (SDE)

$$d\mathbf{x} = \mathbf{f}(\mathbf{x}, t)dt + g(t)d\mathbf{w} \qquad (24.61)$$

where $\mathbf{f}(\mathbf{x}, t)$ is the drift coefficient, $g(t)$ is the diffusion coefficient, and $\mathbf{w}$ is a Wiener process. The perturbed data distribution at time $t$ is denoted as $q_t(\mathbf{x})$, and the score function $\nabla_x \log q_t(\mathbf{x})$ is approximated by a neural network $\mathbf{s}_\theta(\mathbf{x}, t)$ with parameters $\theta$. The training objective is derived from denoising score matching (DSM), which minimizes the expected squared error between the model's score and the score of the perturbation kernel $q(\mathbf{x}_t|\mathbf{x}_0)$. For a given noise schedule $\beta(t)$, the perturbation kernel is

$$q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \alpha_t \mathbf{x}_0, \sigma_t^2 \mathbf{I}) \quad \text{where} \quad \alpha_t = e^{-\frac{1}{2} \int_0^t \beta(s) ds}, \sigma_t^2 = 1 - e^{-\int_0^t \beta(s) ds} \qquad (24.62)$$

The DSM objective is then expressed as

$$\mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \mathbf{x}_t} \left[ \lambda(t) \|\mathbf{s}_\theta(\mathbf{x}_t, t) - \nabla_{\mathbf{x}_t} \log q(\mathbf{x}_t|\mathbf{x}_0)\|^2 \right] \qquad (24.63)$$

where $\lambda(t)$ is a weighting function, often chosen as $\lambda(t) = g(t)^2$ or $\lambda(t) = \sigma_t^2$ to balance the influence of different noise levels. The gradient of the log-perturbation kernel is

$$\nabla_{\mathbf{x}_t} \log q(\mathbf{x}_t|\mathbf{x}_0) = - \frac{\mathbf{x}_t - \alpha_t \mathbf{x}_0}{\sigma_t^2} \qquad (24.64)$$

which simplifies the objective to

$$\mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \epsilon} \left[ \lambda(t) \left\| \mathbf{s}_{\theta}(\alpha_t \mathbf{x}_0 + \sigma_t \epsilon, t) + \frac{\epsilon}{\sigma_t} \right\|^2 \right] \qquad (24.65)$$

where $\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$.

In score-based generative models, the training objective generalizes to multiple noise scales $\{\sigma_i\}_{i=1}^L$, where each scale corresponds to a perturbed data distribution $p_{\sigma_i}(\mathbf{x}) = \int p_{\text{data}}(\mathbf{y})\mathcal{N}(\mathbf{x}; \mathbf{y}, \sigma_i^2 \mathbf{I})d\mathbf{y}$. The score function $\nabla_x \log p_{\sigma_i}(\mathbf{x})$ is approximated by a noise-conditioned neural network $\mathbf{s}_\theta(\mathbf{x}, \sigma_i)$, and the training objective becomes a weighted sum of score matching losses across all noise scales:

$$\mathcal{J}(\theta) = \sum_{i=1}^{L} \lambda(\sigma_i) \mathbb{E}_{p_{\text{data}}(\mathbf{x})} \mathbb{E}_{\mathbf{x}_{\sigma_i} \sim \mathcal{N}(\mathbf{x}, \sigma_i^2 \mathbf{I})} \left[ \| \mathbf{s}_{\theta}(\mathbf{x}_{\sigma_i}, \sigma_i) - \nabla_{\mathbf{x}_{\sigma_i}} \log p_{\sigma_i}(\mathbf{x}_{\sigma_i} | \mathbf{x}) \|^2 \right] \qquad (24.66)$$

where $\nabla_{\mathbf{x}_{\sigma_i}} \log p_{\sigma_i}(\mathbf{x}_{\sigma_i}|\mathbf{x}) = -\frac{\mathbf{x}_{\sigma_i} - \mathbf{x}}{\sigma_i^2}$. This simplifies to

$$\mathcal{J}(\theta) = \sum_{i=1}^{L} \lambda(\sigma_i) \mathbb{E}_{p_{\text{data}}(\mathbf{x})} \mathbb{E}_{\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})} \left[ \left\| \mathbf{s}_{\theta}(\mathbf{x} + \sigma_i \epsilon, \sigma_i) + \frac{\epsilon}{\sigma_i} \right\|^2 \right] \qquad (24.67)$$

The weighting function $\lambda(\sigma_i)$ is typically chosen to ensure balanced contributions across noise scales, such as

$$\lambda(\sigma_i) = \sigma_i^2 \quad \text{or} \quad \lambda(\sigma_i) = \frac{1}{\mathbb{E}[\|\nabla_x \log p_{\sigma_i}(\mathbf{x})\|^2]} \qquad (24.68)$$
The connection between diffusion models and score-based generative models becomes evident when considering the continuous limit of noise scales. In this limit, the noise schedule $\sigma(t)$ is replaced by a continuous function, and the training objective becomes

$$ \mathcal{J}(\theta) = \mathbb{E}_{t \sim \mathcal{U}(0,T)} \mathbb{E}_{\mathbf{x}_0 \sim p_{\text{data}}} \mathbb{E}_{\epsilon \sim \mathcal{N}(0,\mathbf{I})} \left[ \lambda(t) \left\| \mathbf{s}_{\theta}(\mathbf{x}_0 + \sigma(t)\epsilon, t) + \frac{\epsilon}{\sigma(t)} \right\|^2 \right] \quad (24.69) $$

where $\sigma(t)$ is a monotonically increasing function of $t$. This objective aligns with the DSM objective in diffusion models when $\sigma(t)$ is chosen to match the marginal variance of the forward process. The score function $\mathbf{s}_{\theta}(\mathbf{x}, t)$ is then trained to approximate $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ for all $t \in [0, T]$, enabling the reverse-time SDE or probability flow ODE to generate samples from $p_{\text{data}}(\mathbf{x})$.

Theoretical analysis of the training objective reveals its connection to variational inference and the evidence lower bound (ELBO). For diffusion models, the ELBO can be decomposed into a sum of terms involving the score matching objective and additional regularization terms. Specifically, the ELBO is

$$ \log p_{\theta}(\mathbf{x}_0) \geq \mathbb{E}_{q(\mathbf{x}_{1:T}|\mathbf{x}_0)} [\log p_{\theta}(\mathbf{x}_{0:T}) - \log q(\mathbf{x}_{1:T}|\mathbf{x}_0)] \quad (24.70) $$

where $p_{\theta}(\mathbf{x}_{0:T})$ is the joint distribution defined by the reverse process, and $q(\mathbf{x}_{1:T}|\mathbf{x}_0)$ is the forward process. The ELBO can be rewritten as

$$ \log p_{\theta}(\mathbf{x}_0) \geq - \sum_{t=1}^{T} \mathbb{E}_{q(\mathbf{x}_t|\mathbf{x}_0)} [\text{KL}(q(\mathbf{x}_{t-1}|\mathbf{x}_t, \mathbf{x}_0) \| p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_t))] + \text{const} \quad (24.71) $$

where the KL divergence terms correspond to denoising score matching objectives at each timestep. This equivalence underscores the variational interpretation of the training objective, where minimizing the DSM objective is equivalent to maximizing a lower bound on the log-likelihood of the data.

The training objective also admits a geometric interpretation through the lens of optimal transport and gradient flows. The score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ defines a vector field that guides particles from the prior distribution $q_T(\mathbf{x})$ to the data distribution $q_0(\mathbf{x}) = p_{\text{data}}(\mathbf{x})$. The DSM objective minimizes the discrepancy between the learned vector field $\mathbf{s}_{\theta}(\mathbf{x}, t)$ and the true score $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$, ensuring that the reverse process accurately inverts the forward diffusion. This perspective links the training objective to the Jordan-Kinderlehrer-Otto (JKO) scheme for gradient flows in the space of probability measures, where the score function acts as the gradient of the Wasserstein distance functional. In practice, the training objective is optimized using stochastic gradient descent (SGD) or its variants, with gradients estimated via Monte Carlo sampling. The expectation over time $t$ is approximated by uniform sampling, and the expectation over data $\mathbf{x}_0$ and noise $\epsilon$ is approximated by mini-batches. The weighting function $\lambda(t)$ or $\lambda(\sigma_i)$ is critical for stable training, as it balances the contributions of different noise levels. Common choices include $\lambda(t) = \sigma_t^2$ for variance-preserving SDEs or $\lambda(t) = g(t)^2$ for general SDEs. The choice of weighting function influences the trade-off between sample quality and likelihood, with heuristic and theoretical justifications guiding its design. The training objective also extends to conditional generation, where the goal is to model $p(\mathbf{x}|\mathbf{y})$ for some auxiliary variable $\mathbf{y}$. In this case, the score function becomes $\nabla_{\mathbf{x}} \log p_t(\mathbf{x}|\mathbf{y})$, and the training objective is modified to include the conditioning variable:

$$ \mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \mathbf{y}, \epsilon} \left[ \lambda(t) \left\| \mathbf{s}_{\theta}(\mathbf{x}_t, t, \mathbf{y}) + \frac{\epsilon}{\sigma_t} \right\|^2 \right] \quad (24.72) $$

where

$$ \mathbf{x}_t = \alpha_t \mathbf{x}_0 + \sigma_t \epsilon \quad (24.73) $$
This objective enables conditional sampling by integrating the learned score $s_{\theta}(\mathbf{x}, t, \mathbf{y})$ into the reverse-time SDE or probability flow ODE. The training objective's robustness to noise and mis-specification is a key advantage of score-based methods. Unlike likelihood-based models that require exact density evaluation, score matching only requires access to the gradient of the log-density, which can be estimated even when the normalization constant is intractable. This property makes score-based methods particularly suitable for high-dimensional data, where exact likelihood evaluation is computationally prohibitive. The training objective's connection to contrastive divergence and noise-contrastive estimation further highlights its robustness, as it effectively learns to distinguish data from noise perturbations.

Theoretical guarantees for the training objective include consistency and convergence under appropriate regularity conditions. For instance, if the neural network $s_{\theta}(\mathbf{x}, t)$ is sufficiently expressive and the noise schedule is well-designed, the minimizer of $J(\theta)$ will converge to the true score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ as the number of training samples increases. The convergence rate depends on the complexity of the data distribution and the architecture of the score network, with recent advances providing non-asymptotic bounds for specific classes of distributions. These guarantees ensure that the learned score function can reliably guide the reverse process for sample generation.

### 24.3.5 Sampling

Sampling in diffusion models and score-based generative models involves reversing a noising process to generate data from noise through either stochastic or deterministic dynamics. The reverse process is governed by learned score functions that approximate the gradient of the log-density of perturbed data distributions, enabling the transformation of simple noise distributions into complex data distributions. For the reverse-time stochastic differential equation (SDE), the sampling process begins with $\mathbf{x}_T \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ and follows the dynamics

$$d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2\nabla_{\mathbf{x}} \log q_t(\mathbf{x})] dt + g(t)d\bar{\mathbf{w}} \quad (24.74)$$

where $\mathbf{f}(\mathbf{x}, t)$ and $g(t)$ are the drift and diffusion coefficients from the forward SDE, $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ is the score function, and $\bar{\mathbf{w}}$ is a reverse-time Wiener process. The score function is approximated by a neural network $s_{\theta}(\mathbf{x}, t)$, yielding the practical reverse SDE

$$d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2\mathbf{s}_{\theta}(\mathbf{x}, t)] dt + g(t)d\bar{\mathbf{w}} \quad (24.75)$$

This SDE is discretized using numerical methods like the Euler-Maruyama scheme, which updates the state as

$$\mathbf{x}_{t-\Delta t} = \mathbf{x}_t - [\mathbf{f}(\mathbf{x}_t, t) - g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t)] \Delta t + g(t) \sqrt{\Delta t} \mathbf{z} \quad (24.76)$$

where $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ and $\Delta t$ is the step size. An alternative deterministic sampling method uses the probability flow ordinary differential equation (ODE), derived by removing the noise term from the reverse SDE:

$$d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g(t)^2\nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] dt \quad (24.77)$$

With the learned score $s_{\theta}(\mathbf{x}, t)$, this becomes

$$d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g(t)^2\mathbf{s}_{\theta}(\mathbf{x}, t) \right] dt \quad (24.78)$$

The ODE is solved numerically using methods like Runge-Kutta, with the Euler update rule

$$\mathbf{x}_{t-\Delta t} = \mathbf{x}_t - \left[ \mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2}g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t) \right] \Delta t \quad (24.79)$$
The probability flow ODE preserves the same marginal distributions as the reverse SDE but enables more efficient sampling and exact likelihood computation via the instantaneous change-of-variables formula

$$ \log p(\mathbf{x}_0) = \log p(\mathbf{x}_T) + \int_0^T \nabla \cdot \left( \mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2} g(t)^2 \mathbf{s}_\theta(\mathbf{x}_t, t) \right) dt \quad (24.80) $$

where $\nabla \cdot$ denotes divergence.

In discrete-time diffusion models, sampling proceeds through a Markov chain that reverses the forward noising process. Given a noise schedule $\{\beta_k\}_{k=1}^N$, the reverse transition kernel $p_\theta(\mathbf{x}_{k-1}|\mathbf{x}_k)$ is typically Gaussian, parameterized as $\mathcal{N}(\mathbf{x}_{k-1}; \mu_\theta(\mathbf{x}_k, k), \Sigma_\theta(\mathbf{x}_k, k))$. The mean $\mu_\theta(\mathbf{x}_k, k)$ is often derived from the score function, such as

$$ \mu_\theta(\mathbf{x}_k, k) = \frac{1}{\sqrt{1 - \beta_k}} (\mathbf{x}_k + \beta_k \mathbf{s}_\theta(\mathbf{x}_k, k)) \quad (24.81) $$

while the covariance $\Sigma_\theta(\mathbf{x}_k, k)$ may be fixed or learned. Sampling iteratively applies this transition from $\mathbf{x}_N \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ to $\mathbf{x}_0$, with each step denoising the data. For variance-preserving processes, the reverse transition can also be expressed as

$$ \mathbf{x}_{k-1} = \frac{1}{\sqrt{\alpha_k}} \left( \mathbf{x}_k - \frac{\beta_k}{\sqrt{1 - \bar{\alpha}_k}} \mathbf{s}_\theta(\mathbf{x}_k, k) \right) + \sqrt{\beta_k} \mathbf{z} \quad (24.82) $$

where $\alpha_k = 1 - \beta_k$, $\bar{\alpha}_k = \prod_{i=1}^k \alpha_i$, and $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. Score-based generative models without a predefined forward process use Langevin dynamics for sampling. Given a noise-conditional score network $\mathbf{s}_\theta(\mathbf{x}, \sigma)$, Langevin sampling iterates

$$ \mathbf{x}_{i+1} = \mathbf{x}_i + \epsilon_i \mathbf{s}_\theta(\mathbf{x}_i, \sigma_i) + \sqrt{2\epsilon_i} \mathbf{z}_i \quad (24.83) $$

where $\epsilon_i$ is the step size, $\sigma_i$ is the noise scale, and $\mathbf{z}_i \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. This process is run for multiple noise levels $\{\sigma_i\}_{i=1}^L$, typically in descending order, to refine samples from high-noise to low-noise regimes. The noise scales and step sizes are often annealed according to a schedule, such as

$$ \sigma_i = \sigma_{\max}(\sigma_{\min}/\sigma_{\max})^{i/L} \quad \text{and} \quad \epsilon_i = \epsilon_0 \sigma_i^2/\sigma_L^2 \quad (24.84) $$

ensuring stable convergence. The sampling process is closely tied to the noise schedule and the quality of the learned score function. For the variance-preserving SDE with $\beta(t)$, the noise schedule must satisfy $\int_0^T \beta(t)dt \to \infty$ to ensure $q_T(\mathbf{x}) \approx \mathcal{N}(\mathbf{x}; \mathbf{0}, \mathbf{I})$. Common choices include linear schedules

$$ \beta(t) = \beta_{\min} + t(\beta_{\max} - \beta_{\min}) \quad (24.85) $$

or cosine schedules

$$ \alpha_t = \cos^2(\pi t / 2T) \quad (24.86) $$

The score function's accuracy is critical, as errors accumulate during sampling, leading to degraded sample quality. This is mitigated by techniques like predictor-corrector methods, which alternate between SDE steps (predictor) and Langevin steps (corrector) to refine the sample trajectory. The corrector step involves running a few iterations of Langevin dynamics at fixed noise levels to correct deviations from the true score.

Theoretical analysis of sampling reveals convergence guarantees under certain conditions. For Langevin dynamics, the distribution of samples converges to the target $p_{\text{data}}(\mathbf{x})$ if the score error

$$ \mathbb{E}[\|\mathbf{s}_\theta(\mathbf{x}, \sigma) - \nabla_{\mathbf{x}} \log p_\sigma(\mathbf{x})\|^2] \quad (24.87) $$
is bounded and the step sizes {$\epsilon_i$} satisfy

$$ \sum_{i} \epsilon_i = \infty \quad \text{and} \quad \sum_{i} \epsilon_i^2 < \infty \qquad (24.88) $$

For reverse SDEs, the Girsanov theorem provides a measure-theoretic foundation, showing that the path measures of the forward and reverse processes are mutually absolutely continuous, with the Radon-Nikodym derivative depending on the score error. The sampling error can be quantified via the KL divergence between the generated and target distributions, which scales with the time-integrated score error

$$ \int_{0}^{T} g(t)^2 \mathbb{E}[\|\mathbf{s}_{\theta}(\mathbf{x}_{t}, t) - \nabla_{\mathbf{x}_{t}} \log q_{t}(\mathbf{x}_{t})\|^2] dt \qquad (24.89) $$

Practical sampling often involves trade-offs between computational cost and sample quality. The Euler-Maruyama discretization of the reverse SDE introduces error proportional to $\Delta t$, while higher-order methods like Milstein or stochastic Runge-Kutta reduce this error at increased computational expense. The choice of sampler depends on the application: ODE-based samplers are preferred for likelihood computation, while SDE-based samplers may yield higher sample diversity. Hybrid methods, such as combining ODE trajectories with occasional noise injections, balance these trade-offs. The sampling process also extends to conditional generation by replacing the score $\mathbf{s}_{\theta}(\mathbf{x}, t)$ with a conditional score $\mathbf{s}_{\theta}(\mathbf{x}, t, \mathbf{y})$, where $\mathbf{y}$ is an auxiliary input (e.g., class labels or text embeddings). In summary, sampling in diffusion and score-based models involves reversing a noising process through learned score functions, implemented via SDEs, ODEs, or Markov chains. The core equations include the reverse SDE

$$ d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}, t)] dt + g(t) d\bar{\mathbf{w}} \qquad (24.90) $$

the probability flow ODE

$$ d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2} g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}, t) \right] dt \qquad (24.91) $$

and the Langevin update

$$ \mathbf{x}_{i+1} = \mathbf{x}_i + \epsilon_i \mathbf{s}_{\theta}(\mathbf{x}_i, \sigma_i) + \sqrt{2\epsilon_i} \mathbf{z}_i \qquad (24.92) $$

These methods are grounded in stochastic calculus, dynamical systems, and statistical physics, providing a rigorous framework for generative modeling. The sampling quality depends on the accuracy of the learned score, the noise schedule, and the numerical integration method, with theoretical guarantees ensuring convergence under appropriate conditions. The flexibility of these approaches enables applications ranging from image synthesis to molecular design, with ongoing advances improving efficiency and scalability.

## 24.3.6 Score-Based Generative Models

Score-based generative models and diffusion models are fundamentally interconnected through their shared reliance on learning score functions, which are gradients of log-probability densities that guide the generation of samples from noise. The score function $\nabla_{\mathbf{x}} \log p_t(\mathbf{x})$ is central to both frameworks, where $p_t(\mathbf{x})$ represents the data distribution perturbed by noise at time $t$. In score-based generative models, the score function is learned directly through denoising score matching, which minimizes the objective

$$ \mathcal{J}(\theta) = \mathbb{E}_{t \sim \mathcal{U}(0,T)} \mathbb{E}_{\mathbf{x} \sim p_t(\mathbf{x})} \mathbb{E}_{\mathbf{x}_0 \sim p_{\text{data}}(\mathbf{x})} \left[ \lambda(t) \|\mathbf{s}_{\theta}(\mathbf{x}, t) - \nabla_{\mathbf{x}} \log p_t(\mathbf{x} | \mathbf{x}_0) \|^2 \right] \qquad (24.93) $$

where $\lambda(t)$ is a weighting function and $\nabla_{\mathbf{x}} \log p_t(\mathbf{x} | \mathbf{x}_0)$ is the score of the perturbation kernel. For Gaussian perturbations, the perturbation kernel is

$$ p_t(\mathbf{x} | \mathbf{x}_0) = \mathcal{N}(\mathbf{x}; \mathbf{x}_0, \sigma^2(t)\mathbf{I}) \qquad (24.94) $$
and the score simplifies to

$$ \nabla_{\mathbf{x}} \log p_t(\mathbf{x}|\mathbf{x}_0) = -\frac{\mathbf{x} - \mathbf{x}_0}{\sigma^2(t)} \qquad (24.95) $$

The denoising score matching objective thus becomes

$$ \mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \mathbf{x}_t} \left[ \lambda(t) \left\| \mathbf{s}_{\theta}(\mathbf{x}_t, t) + \frac{\mathbf{x}_t - \mathbf{x}_0}{\sigma^2(t)} \right\|^2 \right] \qquad (24.96) $$

where $\mathbf{x}_t = \mathbf{x}_0 + \sigma(t)\epsilon$ and $\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. Diffusion models formalize the perturbation process through a forward stochastic differential equation (SDE)

$$ d\mathbf{x} = \mathbf{f}(\mathbf{x}, t)dt + g(t)d\mathbf{w} \qquad (24.97) $$

where $\mathbf{f}(\mathbf{x}, t)$ is the drift coefficient, $g(t)$ is the diffusion coefficient, and $\mathbf{w}$ is a Wiener process. The forward SDE transforms the data distribution $p_{\text{data}}(\mathbf{x})$ into a simple noise distribution $p_T(\mathbf{x})$ (typically Gaussian) by gradually adding noise. The reverse process, which generates samples from $p_{\text{data}}(\mathbf{x})$, is described by the reverse-time SDE

$$ d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g^2(t)\nabla_{\mathbf{x}} \log p_t(\mathbf{x})] dt + g(t)d\bar{\mathbf{w}} \qquad (24.98) $$

where $\bar{\mathbf{w}}$ is a reverse-time Wiener process. The score function $\nabla_{\mathbf{x}} \log p_t(\mathbf{x})$ is approximated by a neural network $\mathbf{s}_{\theta}(\mathbf{x}, t)$, trained using the objective

$$ \mathcal{J}(\theta) = \mathbb{E}_{t \sim \mathcal{U}(0,T)} \mathbb{E}_{\mathbf{x}_0 \sim p_{\text{data}}(\mathbf{x})} \mathbb{E}_{\mathbf{x}_t \sim p_t(\mathbf{x}|\mathbf{x}_0)} [\lambda(t) \|\mathbf{s}_{\theta}(\mathbf{x}_t, t) - \nabla_{\mathbf{x}_t} \log p_t(\mathbf{x}_t|\mathbf{x}_0)\|^2] \qquad (24.99) $$

For the variance-preserving SDE, where $\mathbf{f}(\mathbf{x}, t) = -\frac{1}{2}\beta(t)\mathbf{x}$ and $g(t) = \sqrt{\beta(t)}$, the perturbation kernel is

$$ p_t(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \sqrt{\bar{\alpha}(t)}\mathbf{x}_0, (1-\bar{\alpha}(t))\mathbf{I}) \qquad (24.100) $$

where $\bar{\alpha}(t) = e^{-\int_0^t \beta(s)ds}$, and the score becomes

$$ \nabla_{\mathbf{x}_t} \log p_t(\mathbf{x}_t|\mathbf{x}_0) = -\frac{\mathbf{x}_t - \sqrt{\bar{\alpha}(t)}\mathbf{x}_0}{1 - \bar{\alpha}(t)} \qquad (24.101) $$

The training objective thus reduces to

$$ \mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \epsilon} \left[ \lambda(t) \left\| \mathbf{s}_{\theta}(\sqrt{\bar{\alpha}(t)}\mathbf{x}_0 + \sqrt{1-\bar{\alpha}(t)}\epsilon, t) + \frac{\epsilon}{\sqrt{1-\bar{\alpha}(t)}} \right\|^2 \right] \qquad (24.102) $$

where $\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$.

The probability flow ordinary differential equation (ODE) is a deterministic alternative to the reverse-time SDE, derived by removing the noise term:

$$ d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g^2(t)\nabla_{\mathbf{x}} \log p_t(\mathbf{x}) \right] dt \qquad (24.103) $$

This ODE describes a deterministic trajectory that preserves the same marginal distributions as the reverse SDE. The probability flow ODE is particularly useful for likelihood computation, as it enables the application of the instantaneous change-of-variables formula

$$ \log p_0(\mathbf{x}_0) = \log p_T(\mathbf{x}_T) + \int_0^T \nabla \cdot \left( \mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2}g^2(t)\nabla_{\mathbf{x}_t} \log p_t(\mathbf{x}_t) \right) dt \qquad (24.104) $$
where $\nabla$ denotes the divergence operator. The divergence term can be approximated using Hutchinson's trace estimator, which reduces the computational cost from $O(d^2)$ to $O(d)$ for $d$-dimensional data. In score-based generative models, the forward process is often implicit, defined through a sequence of noise scales $\{\sigma_i\}_{i=1}^L$. The perturbed data distribution at scale $\sigma_i$ is

$$p_{\sigma_i}(\mathbf{x}) = \int p_{\text{data}}(\mathbf{y}) \mathcal{N}(\mathbf{x}; \mathbf{y}, \sigma_i^2 \mathbf{I}) d\mathbf{y} \qquad (24.105)$$

and the score $\nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x})$ is learned directly. The training objective is a weighted sum of score matching losses across all noise scales:

$$\mathcal{J}(\theta) = \sum_{i=1}^{L} \lambda(\sigma_i) \mathbb{E}_{p_{\text{data}}(\mathbf{x})} \mathbb{E}_{\mathbf{x}_{\sigma_i} \sim \mathcal{N}(\mathbf{x}, \sigma_i^2 \mathbf{I})} \left[ \| \mathbf{s}_{\theta}(\mathbf{x}_{\sigma_i}, \sigma_i) - \nabla_{\mathbf{x}_{\sigma_i}} \log p_{\sigma_i}(\mathbf{x}_{\sigma_i} | \mathbf{x}) \|^2 \right] \qquad (24.106)$$

where $\nabla_{\mathbf{x}_{\sigma_i}} \log p_{\sigma_i}(\mathbf{x}_{\sigma_i}|\mathbf{x}) = -\frac{\mathbf{x}_{\sigma_i} - \mathbf{x}}{\sigma_i^2}$. This simplifies to

$$\mathcal{J}(\theta) = \sum_{i=1}^{L} \lambda(\sigma_i) \mathbb{E}_{p_{\text{data}}(\mathbf{x})} \mathbb{E}_{\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})} \left[ \left\| \mathbf{s}_{\theta}(\mathbf{x} + \sigma_i \epsilon, \sigma_i) + \frac{\epsilon}{\sigma_i} \right\|^2 \right] \qquad (24.107)$$

The noise scales $\{\sigma_i\}$ are typically chosen to span a wide range, from $\sigma_1 \approx 0$ (near the data distribution) to $\sigma_L \gg 0$ (near the noise distribution).

Sampling in score-based generative models is performed using Langevin dynamics, which iteratively refines samples by following the score function. Given a noise-conditional score network $\mathbf{s}_{\theta}(\mathbf{x}, \sigma)$, Langevin dynamics updates the sample as

$$\mathbf{x}_{k+1} = \mathbf{x}_k + \epsilon_k \mathbf{s}_{\theta}(\mathbf{x}_k, \sigma) + \sqrt{2\epsilon_k} \mathbf{z}_k \qquad (24.108)$$

where $\epsilon_k$ is the step size and $\mathbf{z}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. This process is repeated for each noise scale $\sigma_i$ in descending order, gradually reducing the noise level to produce high-quality samples. The step sizes $\{\epsilon_k\}$ are typically chosen to satisfy

$$\sum_k \epsilon_k = \infty \quad \text{and} \quad \sum_k \epsilon_k^2 < \infty \qquad (24.109)$$

to ensure convergence to the target distribution. The connection between score-based generative models and diffusion models is further illuminated by considering the continuous limit of noise scales. In this limit, the noise schedule $\sigma(t)$ becomes a continuous function, and the training objective becomes

$$\mathcal{J}(\theta) = \mathbb{E}_{t \sim \mathcal{U}(0,T)} \mathbb{E}_{\mathbf{x}_0 \sim p_{\text{data}}(\mathbf{x})} \mathbb{E}_{\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})} \left[ \lambda(t) \left\| \mathbf{s}_{\theta}(\mathbf{x}_0 + \sigma(t)\epsilon, t) + \frac{\epsilon}{\sigma(t)} \right\|^2 \right] \qquad (24.110)$$

This objective aligns with the denoising score matching objective in diffusion models when $\sigma(t)$ is chosen to match the marginal variance of the forward process. The score function $\mathbf{s}_{\theta}(\mathbf{x}, t)$ thus approximates $\nabla_{\mathbf{x}} \log p_t(\mathbf{x})$ for all $t \in [0, T]$, enabling the reverse-time SDE or probability flow ODE to generate samples from $p_{\text{data}}(\mathbf{x})$.

Theoretical guarantees for score-based generative models and diffusion models rely on the accuracy of the learned score function. Under suitable regularity conditions, the minimizer of the denoising score matching objective satisfies

$$\mathbf{s}_{\theta}(\mathbf{x}, t) \approx \nabla_{\mathbf{x}} \log p_t(\mathbf{x}) \qquad (24.111)$$
ensuring that the reverse process converges to the true data distribution. The quality of the generated samples depends on the expressiveness of the score network, the noise schedule, and the numerical integration method used for sampling. Recent advances have established non-asymptotic convergence rates for score-based generative models, demonstrating that the sample quality improves with the number of training iterations and the dimensionality of the data. These results highlight the robustness and scalability of score-based methods for high-dimensional generative modeling.

### 24.3.7 Langevin Dynamics for Sampling

Langevin dynamics serves as a fundamental sampling mechanism in both diffusion models and score-based generative models, providing a stochastic process that leverages the learned score function to generate samples from the target data distribution. The core principle of Langevin dynamics lies in its iterative update rule, which combines gradient ascent along the log-probability density with Gaussian noise injection to explore the data manifold. Given a target distribution $p(\mathbf{x})$ and its score function $\nabla_{\mathbf{x}} \log p(\mathbf{x})$, the Langevin dynamics update is formulated as

$$ \mathbf{x}_{k+1} = \mathbf{x}_k + \epsilon \nabla_{\mathbf{x}} \log p(\mathbf{x}_k) + \sqrt{2\epsilon} \mathbf{z}_k \quad (24.112) $$

where $\epsilon$ is the step size, $\mathbf{z}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ is standard Gaussian noise, and $k$ indexes the iteration. This update rule can be interpreted as a discretization of the continuous-time stochastic differential equation (SDE)

$$ d\mathbf{x} = \nabla_{\mathbf{x}} \log p(\mathbf{x}) dt + \sqrt{2\epsilon} d\mathbf{w} \quad (24.113) $$

where $\mathbf{w}$ is a Wiener process. Under mild regularity conditions, the stationary distribution of this SDE converges to $p(\mathbf{x})$, ensuring that the samples generated by Langevin dynamics asymptotically follow the target distribution.

In the context of score-based generative models, Langevin dynamics is employed to sample from noise-perturbed data distributions $p_{\sigma_i}(\mathbf{x})$, where $\sigma_i$ represents a noise scale. The score function $\nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x})$ is approximated by a neural network $\mathbf{s}_{\theta}(\mathbf{x}, \sigma_i)$, trained via denoising score matching. The sampling process involves running Langevin dynamics at each noise level $\sigma_i$ in descending order, starting from a high-noise regime where

$$ p_{\sigma_L}(\mathbf{x}) \approx \mathcal{N}(\mathbf{x}; \mathbf{0}, \sigma_L^2 \mathbf{I}) \quad (24.114) $$

and progressively refining the samples to lower-noise regimes. The update rule for noise-conditional Langevin dynamics is

$$ \mathbf{x}_{k+1} = \mathbf{x}_k + \epsilon_i \mathbf{s}_{\theta}(\mathbf{x}_k, \sigma_i) + \sqrt{2\epsilon_i} \mathbf{z}_k \quad (24.115) $$

where $\epsilon_i$ is a step size specific to noise scale $\sigma_i$. The noise scales $\{\sigma_i\}_{i=1}^L$ are typically chosen to form a geometric sequence, $\sigma_i = \sigma_{\max}(\sigma_{\min}/\sigma_{\max})^{i/L}$, ensuring smooth transitions between noise levels. The step sizes $\{\epsilon_i\}_{i=1}^L$ are often set proportional to $\sigma_i^2$, reflecting the higher curvature of the log-density in low-noise regimes. For diffusion models, Langevin dynamics is implicitly embedded within the reverse-time SDE framework. The reverse-time SDE

$$ d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2 \nabla_{\mathbf{x}} \log p_t(\mathbf{x})] dt + g(t) d\bar{\mathbf{w}} \quad (24.116) $$

can be viewed as a continuous generalization of Langevin dynamics, where the drift term $\mathbf{f}(\mathbf{x}, t) - g(t)^2 \nabla_{\mathbf{x}} \log p_t(\mathbf{x})$ combines the deterministic drift from the forward SDE with the score-based correction. Discretizing this SDE using the Euler-Maruyama method yields an update rule analogous to Langevin dynamics:

$$ \mathbf{x}_{t-\Delta t} = \mathbf{x}_t - [\mathbf{f}(\mathbf{x}_t, t) - g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t)] \Delta t + g(t) \sqrt{\Delta t} \mathbf{z} \quad (24.117) $$
where $z \sim \mathcal{N}(0, I)$. This discretization highlights the connection between Langevin dynamics and reverse-time diffusion, as both involve a score-guided drift term and additive Gaussian noise. The key distinction lies in the additional drift term $f(x, t)$ in diffusion models, which accounts for the forward process's deterministic component.

The convergence properties of Langevin dynamics are rigorously characterized by its mixing time and error analysis. For a target distribution $p(x)$ satisfying log-Sobolev inequalities, the discretized Langevin dynamics converges exponentially fast in the Wasserstein distance, with a mixing time dependent on the condition number of the Hessian of $\log p(x)$. The mean-squared error between the empirical distribution of samples and the target distribution scales as $O(\epsilon + \sqrt{d\epsilon})$, where $d$ is the dimensionality of $x$ and $\epsilon$ is the step size. In score-based generative models, the error introduced by the score approximation

$$ \mathbf{s}_{\theta}(\mathbf{x}, \sigma_i) \approx \nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x}) \qquad (24.118) $$

further impacts convergence. Theoretical bounds on the sampling error typically involve terms like $\mathbb{E}[\|\mathbf{s}_{\theta}(\mathbf{x}, \sigma_i) - \nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x})\|^2]$, which reflect the quality of the score matching objective. Practical implementations of Langevin dynamics often incorporate Metropolis-Hastings corrections or adaptive step sizes to improve stability and convergence. The Metropolis-adjusted Langevin algorithm (MALA) introduces an acceptance-rejection step to correct for discretization errors, ensuring detailed balance with respect to $p(x)$. The update rule for MALA is

$$ \mathbf{x}_{k+1} = \mathbf{x}_k + \frac{\epsilon}{2} \nabla_{\mathbf{x}} \log p(\mathbf{x}_k) + \sqrt{\epsilon} \mathbf{z}_k \qquad (24.119) $$

followed by an acceptance probability

$$ \alpha = \min \left( 1, \frac{p(\mathbf{x}_{k+1})q(\mathbf{x}_k|\mathbf{x}_{k+1})}{p(\mathbf{x}_k)q(\mathbf{x}_{k+1}|\mathbf{x}_k)} \right) \qquad (24.120) $$

where

$$ q(\mathbf{x}'|\mathbf{x}) = \mathcal{N}(\mathbf{x}'; \mathbf{x} + \frac{\epsilon}{2}\nabla_{\mathbf{x}} \log p(\mathbf{x}), \epsilon I) \qquad (24.121) $$

is the proposal distribution. This adjustment guarantees convergence to the correct stationary distribution but incurs additional computational overhead.

In high-dimensional spaces, Langevin dynamics faces challenges related to curvature and ill-conditioning. Preconditioned Langevin dynamics addresses these issues by introducing a position-dependent preconditioning matrix $C(x)$, leading to the update rule

$$ \mathbf{x}_{k+1} = \mathbf{x}_k + \epsilon \mathbf{C}(\mathbf{x}_k) \nabla_{\mathbf{x}} \log p(\mathbf{x}_k) + \sqrt{2\epsilon \mathbf{C}(\mathbf{x}_k)} \mathbf{z}_k \qquad (24.122) $$

The matrix $C(x)$ is often chosen to approximate the inverse Hessian of $-\log p(x)$, reducing the effective condition number of the sampling problem. In score-based generative models, similar preconditioning can be applied by scaling the noise injection term according to the noise level $\sigma_i$, yielding

$$ \mathbf{x}_{k+1} = \mathbf{x}_k + \epsilon_i \mathbf{C}(\mathbf{x}_k, \sigma_i) \mathbf{s}_{\theta}(\mathbf{x}_k, \sigma_i) + \sqrt{2\epsilon_i \mathbf{C}(\mathbf{x}_k, \sigma_i)} \mathbf{z}_k \qquad (24.123) $$

The interplay between Langevin dynamics and diffusion models is further illuminated by the concept of annealed Langevin sampling, where the noise level $\sigma_i$ is gradually decreased during sampling. This annealing process mirrors the reverse-time diffusion trajectory, with each noise level corresponding to a different temperature in the sampling process. The annealed Langevin update rule is

$$ \mathbf{x}_{k+1} = \mathbf{x}_k + \epsilon_i \mathbf{s}_{\theta}(\mathbf{x}_k, \sigma_i) + \sqrt{2\epsilon_i} \mathbf{z}_k \qquad (24.124) $$

where $\sigma_i$ decreases according to a predefined schedule. The step size $\epsilon_i$ is typically set to $\epsilon_i = \eta\sigma_i^2/\sigma_L^2$, where $\eta$ is a constant scaling factor. This annealing strategy ensures that the sampler
initially explores broad regions of the data space (high $\sigma_i$) and later refines the samples to capture fine details (low $\sigma_i$).

Theoretical analysis of annealed Langevin dynamics establishes convergence guarantees under assumptions about the score approximation error and noise schedule. For a sequence of noise levels $\sigma_1 > \sigma_2 > \cdots > \sigma_L$, the total variation distance between the generated samples and the target distribution $p_{\text{data}}(\mathbf{x})$ is bounded by the sum of errors at each noise level $\sum_{i=1}^L \mathbb{E}[\|\mathbf{s}_{\theta}(\mathbf{x}, \sigma_i) - \nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x})\|^2]$. The mathematical formulation of the given statement can be expressed as:

$$d_{\text{TV}}(q_{\theta}(\mathbf{x}), p_{\text{data}}(\mathbf{x})) \le C \sum_{i=1}^{L} \mathbb{E}_{\mathbf{x} \sim p_{\sigma_i}(\mathbf{x})} [\|\mathbf{s}_{\theta}(\mathbf{x}, \sigma_i) - \nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x})\|^2] \quad (24.125)$$

where $d_{\text{TV}}$ denotes the total variation distance between distributions, $q_{\theta}(\mathbf{x})$ represents the distribution of generated samples, $p_{\text{data}}(\mathbf{x})$ is the target data distribution, $C$ is a positive constant depending on the noise schedule and other problem parameters, $p_{\sigma_i}(\mathbf{x})$ is the noise-perturbed data distribution at level $\sigma_i$, $\mathbf{s}_{\theta}(\mathbf{x}, \sigma_i)$ is the learned score function, and $\nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x})$ is the true score of the perturbed distribution. The expectation is taken with respect to $\mathbf{x} \sim p_{\sigma_i}(\mathbf{x})$, and the norm is typically the $L^2$ norm. This inequality shows that the approximation error in the score function at each noise level contributes additively to the upper bound on the total variation distance between the generated and target distributions. A more precise version might include dimension-dependent factors and tighter constants:

$$d_{\text{TV}}(q_{\theta}, p_{\text{data}}) \le \frac{1}{2} \sum_{i=1}^{L} \sigma_i^{-2} \mathbb{E} \left[ \| \sigma_i^2 \mathbf{s}_{\theta}(\mathbf{x}, \sigma_i) - \sigma_i^2 \nabla_{\mathbf{x}} \log p_{\sigma_i}(\mathbf{x}) \|^2 \right]^{1/2} \quad (24.126)$$

where the normalization by $\sigma_i^{-2}$ accounts for the scaling of the score functions with noise level. This reflects how errors at different noise scales contribute proportionally to their inverse variance. This bound highlights the importance of accurate score estimation across all noise scales for high-quality sampling. Langevin dynamics also admits a variational interpretation through the framework of Stein discrepancies. The Stein discrepancy measures the difference between the sample distribution $q(\mathbf{x})$ and the target $p(\mathbf{x})$ via the operator

$$\mathcal{T}_p \phi(\mathbf{x}) = \nabla_{\mathbf{x}} \log p(\mathbf{x})^{\top} \phi(\mathbf{x}) + \nabla_{\mathbf{x}} \cdot \phi(\mathbf{x}) \quad (24.127)$$

where $\phi$ is a test function. The kernelized Stein discrepancy (KSD) minimizes

$$\mathbb{E}_{\mathbf{x}, \mathbf{x}' \sim q} [k(\mathbf{x}, \mathbf{x}') \mathcal{T}_p \mathcal{T}_p k(\mathbf{x}, \mathbf{x}')] \quad (24.128)$$

where $k$ is a positive-definite kernel. Langevin dynamics can be seen as minimizing the KSD by aligning the sample distribution with the gradient flow of $\log p(\mathbf{x})$.

In summary, Langevin dynamics provides a rigorous and flexible framework for sampling in both score-based generative models and diffusion models. Its iterative update rule

$$\mathbf{x}_{k+1} = \mathbf{x}_k + \epsilon \nabla_{\mathbf{x}} \log p(\mathbf{x}_k) + \sqrt{2\epsilon} \mathbf{z}_k \quad (24.129)$$

combines gradient information with stochastic noise to explore the data distribution. The connection to reverse-time SDEs

$$d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2 \nabla_{\mathbf{x}} \log p_t(\mathbf{x})] dt + g(t) d\bar{\mathbf{w}} \quad (24.130)$$

underscores the role of Langevin dynamics as a discrete approximation of continuous diffusion processes. Theoretical guarantees, preconditioning strategies, and annealing techniques further enhance its applicability to high-dimensional generative modeling. The mathematical rigor and versatility of Langevin dynamics make it a cornerstone of modern generative AI methodologies.
24.3.8 Connection to Diffusion Models

Diffusion models can be rigorously understood as a specialized instance of score-based generative models where the forward process is explicitly predefined through a structured noise perturbation scheme. The fundamental connection arises from the observation that both frameworks rely on learning score functions, but diffusion models impose specific constraints on the form of the forward process while score-based models maintain greater generality in their noise perturbation approach. In score-based models, the data distribution $p_{\text{data}}(\mathbf{x})$ is perturbed through an arbitrary family of noise distributions

$$p_{\sigma_i}(\mathbf{x}) = \int p_{\text{data}}(\mathbf{y}) \mathcal{N}(\mathbf{x}; \mathbf{y}, \sigma_i^2 \mathbf{I}) d\mathbf{y} \qquad (24.131)$$

for a discrete set of noise scales $\{\sigma_i\}_{i=1}^L$ or through a continuous noise schedule

$$p_t(\mathbf{x}) = \int p_{\text{data}}(\mathbf{y}) p_t(\mathbf{x}|\mathbf{y}) d\mathbf{y} \qquad (24.132)$$

where $t \in [0, T]$. The score function $\nabla_{\mathbf{x}} \log p_t(\mathbf{x})$ is learned without explicit reference to any underlying stochastic process, merely requiring that the noise-perturbed distributions $p_t(\mathbf{x})$ become increasingly Gaussian as $t$ approaches $T$.

In contrast, diffusion models explicitly define the forward process as a Markov chain or stochastic differential equation that systematically perturbs the data distribution. For continuous-time diffusion models, the forward process is specified by a stochastic differential equation

$$d\mathbf{x} = \mathbf{f}(\mathbf{x}, t)dt + g(t)d\mathbf{w} \qquad (24.133)$$

where $\mathbf{w}$ is a Wiener process, $\mathbf{f}(\mathbf{x}, t)$ is the drift coefficient, and $g(t)$ is the diffusion coefficient. The transition kernel for this process is given by

$$q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \mu(t)\mathbf{x}_0, \Sigma(t)) \qquad (24.134)$$

where $\mu(t)$ and $\Sigma(t)$ are determined by the drift and diffusion coefficients through the solution to the corresponding Fokker-Planck equation. A common choice is the variance-preserving diffusion where

$$\mathbf{f}(\mathbf{x}, t) = -\frac{1}{2}\beta(t)\mathbf{x} \quad \text{and} \quad g(t) = \sqrt{\beta(t)} \qquad (24.135)$$

leading to

$$\mu(t) = e^{-\frac{1}{2} \int_{0}^{t} \beta(s) ds} \quad \text{and} \quad \Sigma(t) = \mathbf{I} - e^{-\int_{0}^{t} \beta(s) ds}$$

The critical observation is that this predefined forward process induces a specific form for the score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$, which must be learned by the model. The mathematical equivalence between the frameworks becomes apparent when examining their respective reverse processes. In score-based models, the reverse process is constructed purely through the learned score function using either Langevin dynamics

$$\mathbf{x}_{k+1} = \mathbf{x}_k + \epsilon_k \mathbf{s}_\theta(\mathbf{x}_k, \sigma_k) + \sqrt{2\epsilon_k} \mathbf{z}_k \qquad (24.136)$$

or through the reverse-time SDE

$$d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2 \nabla_{\mathbf{x}} \log p_t(\mathbf{x})] dt + g(t) d\bar{\mathbf{w}} \qquad (24.137)$$

For diffusion models, the reverse process is exactly given by this same reverse-time SDE but with the crucial difference that the forward process terms $\mathbf{f}(\mathbf{x}, t)$ and $g(t)$ are predefined rather than emergent properties of an arbitrary noise schedule. This means that while score-based models must learn to handle arbitrary noise perturbations, diffusion models benefit from the structural
constraints imposed by their predefined forward process, which manifests in the specific form of the transition kernels $q(\mathbf{x}_t|\mathbf{x}_0)$.

The training objectives further demonstrate this relationship. Score-based models minimize a weighted combination of score matching objectives

$$ \mathcal{J}(\theta) = \sum_{i=1}^{L} \lambda(\sigma_i) \mathbb{E}_{p_{\text{data}}(\mathbf{x})} \mathbb{E}_{\mathbf{x}_{\sigma_i} \sim \mathcal{N}(\mathbf{x}, \sigma_i^2 \mathbf{I})} [\|\mathbf{s}_{\theta}(\mathbf{x}_{\sigma_i}, \sigma_i) - \nabla_{\mathbf{x}_{\sigma_i}} \log p_{\sigma_i}(\mathbf{x}_{\sigma_i} | \mathbf{x})\|^2] \quad (24.138) $$

where

$$ \nabla_{\mathbf{x}_{\sigma_i}} \log p_{\sigma_i}(\mathbf{x}_{\sigma_i} | \mathbf{x}) = - \frac{\mathbf{x}_{\sigma_i} - \mathbf{x}}{\sigma_i^2} \quad (24.139) $$

Diffusion models instead optimize a similar objective but with the transition kernels determined by their forward process:

$$ \mathcal{J}(\theta) = \mathbb{E}_{t, \mathbf{x}_0, \epsilon} [\|\mathbf{s}_{\theta}(\sqrt{\bar{\alpha}_t} \mathbf{x}_0 + \sqrt{1 - \bar{\alpha}_t} \epsilon, t) + \frac{\epsilon}{\sqrt{1 - \bar{\alpha}_t}}\|^2] \quad (24.140) $$

where $\bar{\alpha}_t = e^{-\int_0^t \beta(s) ds}$. The key distinction is that the noise model $\sqrt{\bar{\alpha}_t} \mathbf{x}_0 + \sqrt{1 - \bar{\alpha}_t} \epsilon$ is exactly the solution to the predefined forward SDE rather than an arbitrary Gaussian perturbation. The probability flow ordinary differential equation further illustrates this connection. For score-based models, the probability flow ODE takes the general form

$$ d\mathbf{x} = -\sigma(t)\dot{\sigma}(t)\nabla_{\mathbf{x}} \log p_t(\mathbf{x})dt \quad (24.141) $$

where $\sigma(t)$ is an arbitrary noise schedule. In diffusion models, this becomes

$$ d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - \frac{1}{2}g(t)^2\nabla_{\mathbf{x}} \log q_t(\mathbf{x})]dt \quad (24.142) $$

where the coefficients are precisely those from the predefined forward SDE. The deterministic sampling process thus inherits the structure imposed by the forward diffusion process. The mathematical relationship extends to the perturbation kernels as well. While score-based models work with general perturbed distributions

$$ p_t(\mathbf{x}) = \int p_{\text{data}}(\mathbf{y}) p_t(\mathbf{x} | \mathbf{y}) d\mathbf{y} \quad (24.143) $$

Diffusion models specifically use kernels of the form

$$ 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}) \quad (24.144) $$

for the variance-preserving case or

$$ q(\mathbf{x}_t|\mathbf{x}_0) = \mathcal{N}(\mathbf{x}_t; \mathbf{x}_0, \sigma^2(t)\mathbf{I}) \quad (24.145) $$

for the variance-exploding case. These kernels are solutions to the predefined forward SDEs rather than arbitrary choices, demonstrating how diffusion models constrain the general score-based framework.

The equivalence becomes exact when considering that any diffusion model's forward process induces a corresponding family of noise-perturbed distributions

$$ q_t(\mathbf{x}) = \int q_0(\mathbf{x}_0) q(\mathbf{x}_t|\mathbf{x}_0) d\mathbf{x}_0 \quad (24.146) $$

where $q_0 = p_{\text{data}}$. The score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ that the model learns is identical in form to what a score-based model would learn for this particular noise schedule. However, the diffusion model's training objective is specifically tailored to the known transition kernels $q(\mathbf{x}_t|\mathbf{x}_0)$, whereas a general score-based model must handle arbitrary perturbations. This shows that diffusion models are precisely score-based models where the noise perturbation process is constrained to follow a specific Markovian diffusion process.
24.3.9 Likelihood Computation

The computation of likelihoods in diffusion models and score-based generative models is fundamentally rooted in the change-of-variables formula and the Fokker-Planck equation, providing a rigorous mathematical framework for evaluating the probability density of generated samples. For diffusion models specified by a forward stochastic differential equation

$$
d\mathbf{x} = \mathbf{f}(\mathbf{x}, t)dt + g(t)d\mathbf{w} \qquad (24.147)
$$

the reverse-time process

$$
d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x})] dt + g(t) d\bar{\mathbf{w}} \quad (24.148)
$$

preserves the marginal distributions when the score function $\nabla_{\mathbf{x}} \log q_t(\mathbf{x})$ is exact. The probability flow ordinary differential equation

$$
d\mathbf{x} = \left[ \mathbf{f}(\mathbf{x}, t) - \frac{1}{2} g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \right] dt \quad (24.149)
$$

offers a deterministic trajectory with identical marginals, enabling exact likelihood computation
through the instantaneous change-of-variables formula

$$
\log p_0(\mathbf{x}_0) = \log p_T(\mathbf{x}_T) + \int_0^T \nabla \cdot [\mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2}g(t)^2 \nabla_{\mathbf{x}_t} \log q_t(\mathbf{x}_t)] dt \quad (24.150)
$$

where the divergence term $\nabla \cdot$ captures the infinitesimal volume change along the probability path.
The boundary condition $p_T(\mathbf{x}_T)$ is typically a standard Gaussian $N(0, I)$ by design of the forward
process, with

$$
\log p_T(\mathbf{x}_T) = -\frac{D}{2} \log(2\pi) - \frac{1}{2} \|\mathbf{x}_T\|^2 \quad (24.151)
$$

for $D$-dimensional data. The integral of the divergence can be numerically estimated using Hutchinson's trace estimator $\nabla \cdot \mathbf{v}(\mathbf{x}) = \mathbb{E}_{\epsilon \sim \mathcal{N}(0, I)}[\epsilon^\top \nabla \mathbf{v}(\mathbf{x}) \epsilon]$, where

$$
\mathbf{v}(\mathbf{x}, t) = \mathbf{f}(\mathbf{x}, t) - \frac{1}{2} g(t)^2 \nabla_{\mathbf{x}} \log q_t(\mathbf{x}) \quad (24.152)
$$

reducing the computational complexity from $O(D^2)$ to $O(D)$ for the Jacobian $\nabla \mathbf{v}(\mathbf{x})$. For discrete-time diffusion models with transition kernels $q(\mathbf{x}_t|\mathbf{x}_{t-1}) = \mathcal{N}(\mathbf{x}_t; \sqrt{1-\beta_t}\mathbf{x}_{t-1}, \beta_t\mathbf{I})$, the reverse process likelihood can be expressed through the product of Gaussian transitions

$$
p_{\theta}(\mathbf{x}_{0:T}) = p(\mathbf{x}_{T}) \prod_{t=1}^{T} p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t}) \qquad (24.153)
$$

where $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t}) = \mathcal{N}(\mathbf{x}_{t-1}; \mu_{\theta}(\mathbf{x}_{t}, t), \Sigma_{\theta}(\mathbf{x}_{t}, t))$. The evidence lower bound (ELBO) decomposes as

$$
\log p_{\theta}(\mathbf{x}_0) \geq \mathbb{E}_{q(\mathbf{x}_{1:T}|\mathbf{x}_0)}\left[\log \frac{p_{\theta}(\mathbf{x}_{0:T})}{q(\mathbf{x}_{1:T}|\mathbf{x}_0)}\right] = -\sum_{t=1}^{T} \mathbb{E}_{q(\mathbf{x}_t|\mathbf{x}_0)}\left[D_{KL}\left(q(\mathbf{x}_{t-1}|\mathbf{x}_t, \mathbf{x}_0)\|p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_t)\right)\right] + \text{const} \quad (24.154)
$$

where the KL divergence between Gaussians has closed form

$$
D_{KL}(\mathcal{N}(\mu_1, \Sigma_1) \| \mathcal{N}(\mu_2, \Sigma_2)) = \frac{1}{2} \left[ \log \frac{|\Sigma_2|}{|\Sigma_1|} - D + \text{tr}(\Sigma_2^{-1}\Sigma_1) + (\mu_2 - \mu_1)^{\top}\Sigma_2^{-1}(\mu_2 - \mu_1) \right] \quad (24.155)
$$

The denoising posterior

$$
q(\mathbf{x}_{t-1}|\mathbf{x}_t, \mathbf{x}_0) = \mathcal{N}(\mathbf{x}_{t-1}; \tilde{\mu}_t(\mathbf{x}_t, \mathbf{x}_0), \tilde{\beta}_t \mathbf{I}) \quad (24.156)
$$
has parameters

$$ \tilde{\mu}_t(\mathbf{x}_t, \mathbf{x}_0) = \frac{\sqrt{\alpha_{t-1}}\beta_t}{1 - \alpha_t}\mathbf{x}_0 + \frac{\sqrt{\alpha_t}(1 - \alpha_{t-1})}{1 - \alpha_t}\mathbf{x}_t \quad (24.157) $$

and $\tilde{\beta}_t = \frac{1-\alpha_{t-1}}{1-\alpha_t}\beta_t$, where $\alpha_t = \prod_{s=1}^t (1-\beta_s)$. The model's mean $\mu_\theta(\mathbf{x}_t, t)$ is typically parameterized to match $\tilde{\mu}_t$ through a noise prediction network $\epsilon_\theta(\mathbf{x}_t, t)$, yielding

$$ \mu_{\theta}(\mathbf{x}_t, t) = \frac{1}{\sqrt{1 - \beta_t}} \left( \mathbf{x}_t - \frac{\beta_t}{\sqrt{1 - \alpha_t}} \epsilon_{\theta}(\mathbf{x}_t, t) \right) \quad (24.158) $$

In score-based generative models, the likelihood computation requires solving the Fokker-Planck equation associated with the reverse-time SDE

$$ d\mathbf{x} = [\mathbf{f}(\mathbf{x}, t) - g(t)^2\mathbf{s}_{\theta}(\mathbf{x}, t)]dt + g(t)d\bar{\mathbf{w}} \quad (24.159) $$

where $\mathbf{s}_\theta(\mathbf{x}, t)$ approximates $\nabla_\mathbf{x} \log q_t(\mathbf{x})$. The instantaneous change in log-likelihood along the probability flow ODE is given by

$$ \frac{d}{dt} \log q_t(\mathbf{x}_t) = -\nabla \cdot [\mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2}g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t)] \quad (24.160) $$

leading to

$$ \log q_0(\mathbf{x}_0) = \log q_T(\mathbf{x}_T) - \int_0^T \nabla \cdot [\mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2}g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t)] dt \quad (24.161) $$

For the variance-preserving case where $\mathbf{f}(\mathbf{x}, t) = -\frac{1}{2}\beta(t)\mathbf{x}$ and $g(t) = \sqrt{\beta(t)}$, the divergence term becomes

$$ \nabla \cdot [-\frac{1}{2}\beta(t)\mathbf{x} - \frac{1}{2}\beta(t)\mathbf{s}_{\theta}(\mathbf{x}, t)] = -\frac{D}{2}\beta(t) - \frac{1}{2}\beta(t)\nabla \cdot \mathbf{s}_{\theta}(\mathbf{x}, t) \quad (24.162) $$

The Hutchinson estimator is again employed to compute

$$ \nabla \cdot \mathbf{s}_{\theta}(\mathbf{x}, t) = \mathbb{E}_{\epsilon}[\epsilon^{\top} \nabla \mathbf{s}_{\theta}(\mathbf{x}, t) \epsilon] \quad (24.163) $$

where $\nabla \mathbf{s}_\theta(\mathbf{x}, t)$ is the Jacobian of the score network. The connection to thermodynamic integration emerges when considering the noise-perturbed distributions $p_{\sigma_i}(\mathbf{x})$ in discrete noise level score-based models. The log-likelihood can be expressed through the identity

$$ \log p_{\sigma_{\min}}(\mathbf{x}) = \log p_{\sigma_{\max}}(\mathbf{x}) + \sum_{i=1}^{L} [\log p_{\sigma_{i-1}}(\mathbf{x}) - \log p_{\sigma_i}(\mathbf{x})] \quad (24.164) $$

where each term $\log p_{\sigma_{i-1}}(\mathbf{x}) - \log p_{\sigma_i}(\mathbf{x})$ is estimated via importance sampling or through the score-based relationship

$$ \log p_{\sigma_{i-1}}(\mathbf{x}) - \log p_{\sigma_i}(\mathbf{x}) \approx \frac{1}{2}(\sigma_i^2 - \sigma_{i-1}^2)\mathbf{s}_{\theta}(\mathbf{x}, \sigma_i)^{\top}\mathbf{s}_{\theta}(\mathbf{x}, \sigma_i) + (\sigma_i - \sigma_{i-1})\nabla \cdot \mathbf{s}_{\theta}(\mathbf{x}, \sigma_i) \quad (24.165) $$

The numerical integration of these terms requires careful handling of the variance in Monte Carlo estimates, particularly in high-dimensional spaces where the score function's Jacobian may have large eigenvalues. The exact likelihood computation for continuous-time score-based models involves solving the Feynman-Kac formula associated with the reverse-time diffusion process, where the log-likelihood can be represented as

$$ \log p_0(\mathbf{x}_0) = \log p_T(\mathbf{x}_T) + \mathbb{E} \left[ \int_0^T \frac{1}{2} g(t)^2 \|\mathbf{s}_{\theta}(\mathbf{x}_t, t)\|^2 dt - \nabla \cdot [\mathbf{f}(\mathbf{x}_t, t) - \frac{1}{2} g(t)^2 \mathbf{s}_{\theta}(\mathbf{x}_t, t)] dt \right] \quad (24.166) $$

with the expectation taken over paths from the reverse process. This formulation reveals the connection to the Girsanov theorem for measure changes between the forward and reverse processes,
where the Radon-Nikodym derivative depends on the score error $||s_θ(\mathbf{x}_t, t) - \nabla_{\mathbf{x}_t} \log q_t(\mathbf{x}_t)||^2$. The variance of these estimators is controlled by the accuracy of the learned score function and the stability of the numerical integration scheme, with adaptive step size methods often required to maintain precision in high-dimensional settings.

The geometric interpretation of these likelihood computations arises from considering the probability paths in the space of distributions, where the score function defines a vector field generating diffeomorphisms between noise and data distributions. The log-likelihood difference between endpoints is the integral of the divergence of this vector field along the probability flow, corresponding to the accumulated infinitesimal volume changes. This perspective connects to optimal transport theory, where the score function induces a Brenier potential $\phi(\mathbf{x}, t)$ satisfying

$$ \nabla \phi(\mathbf{x}, t) = \mathbf{x} + \sigma(t)^2 \nabla \log p_t(\mathbf{x}) \qquad (24.167) $$

and the likelihood computation involves the Monge-Ampère determinant

$$ \log \det(I + \sigma(t)^2 \nabla^2 \log p_t(\mathbf{x})) \qquad (24.168) $$

The numerical stability of these computations depends crucially on the Lipschitz properties of the score network and the condition number of the Hessian $\nabla^2 \log p_t(\mathbf{x})$, which governs the local curvature of the log-density. Practical implementations often employ annealed importance sampling (AIS) to bridge between the simple prior distribution and complex data distribution, using the sequence of intermediate distributions $p_t(\mathbf{x})$ defined by either the diffusion or score-based process. The log-likelihood estimate is constructed as

$$ \log p_0(\mathbf{x}_0) = \log p_T(\mathbf{x}_T) + \sum_{t=1}^{T} \log w_t \qquad (24.169) $$

where the importance weights $w_t = \frac{p_{t-1}(\mathbf{x}_{t-1})}{p_t(\mathbf{x}_t)}$ are computed along the sampling trajectory. For diffusion models, these weights simplify to ratios of Gaussian transition probabilities, while for score-based models they involve the score-based density ratios

$$ \frac{p_{t-1}(\mathbf{x})}{p_t(\mathbf{x})} \approx \exp\left(\frac{1}{2}(\sigma_t^2 - \sigma_{t-1}^2)\|s_\theta(\mathbf{x}, t)\|^2 + (\sigma_t - \sigma_{t-1})\nabla \cdot s_\theta(\mathbf{x}, t)\right) \qquad (24.170) $$

The variance of these estimators is minimized when the intermediate distributions are spaced sufficiently close to maintain high overlap, requiring careful tuning of the noise schedule $\sigma(t)$ or $\beta(t)$ to balance computational efficiency with statistical precision. The spectral properties of the score function's Jacobian $\nabla s_\theta(\mathbf{x}, t)$ play a crucial role in the numerical stability of likelihood computations. Large negative eigenvalues can lead to exploding contributions to the divergence term, while positive eigenvalues may indicate instabilities in the probability flow. Regularization techniques such as spectral normalization of the score network or explicit constraints on the Jacobian's eigenvalues $\lambda_i$ through penalty terms $\sum_i(\max(0, \lambda_i - \lambda_{\max}) + \min(0, \lambda_i - \lambda_{\min}))^2$ are often necessary to maintain reliable likelihood estimates. The regularization techniques can be mathematically expressed through the following modifications to the training objective: For spectral normalization of the score network $s_\theta(\mathbf{x}, t)$:

$$ \mathcal{J}_{\text{reg}}(\theta) = \mathcal{J}(\theta) + \gamma \sum_{l} (\|W_l\|_2 - c)^2 \qquad (24.171) $$

where $W_l$ are the weight matrices of the score network, $\| \cdot \|_2$ denotes the spectral norm (largest singular value), $c$ is the desired spectral norm constraint, and $\gamma$ is a regularization strength parameter. For explicit eigenvalue constraints on the Jacobian $\nabla s_\theta(\mathbf{x}, t)$:

$$ \mathcal{J}_{\text{reg}}(\theta) = \mathcal{J}(\theta) + \eta \mathbb{E}_{\mathbf{x},t} \left[ \sum_{i=1}^{d} \left( \max(0, \lambda_i(\mathbf{x}, t) - \lambda_{\max}) + \min(0, \lambda_i(\mathbf{x}, t) - \lambda_{\min}) \right)^2 \right] \qquad (24.172) $$
where $\lambda_i(\mathbf{x}, t)$ are the eigenvalues of $\nabla s_{\theta}(\mathbf{x}, t)$, $\lambda_{\min}$ and $\lambda_{\max}$ form the desired eigenvalue range, $d$ is the input dimension, and $\eta$ controls the penalty strength. The Jacobian eigenvalues can be computed via:

$$ \lambda_i(\mathbf{x}, t) = \epsilon_i^\top (\nabla s_{\theta}(\mathbf{x}, t)) \epsilon_i \qquad (24.173) $$

for eigenvectors $\epsilon_i$ of the Jacobian matrix. A computationally efficient approximation using Hutchinson's estimator for the penalty term would be:

$$ \mathbb{E}_{\epsilon \sim \mathcal{N}(0, I)} \left[ \left( \max(0, \epsilon^\top \nabla s_{\theta}(\mathbf{x}, t) \epsilon - \lambda_{\max}) + \min(0, \epsilon^\top \nabla s_{\theta}(\mathbf{x}, t) \epsilon - \lambda_{\min}) \right)^2 \right] \qquad (24.174) $$

These regularized objectives ensure the score network's Jacobian maintains well-behaved eigenvalues during training, preventing unstable likelihood estimates from:

1. Exploding positive eigenvalues causing local volume expansion

2. Large negative eigenvalues leading to numerical instabilities

3. Ill-conditioned transformations in the probability flow ODE

The Fisher-Legendre transform provides an alternative approach, representing the log-density through its dual potential

$$ \psi(\mathbf{y}) = \sup_{\mathbf{x}} [\mathbf{y}^{\top} \mathbf{x} - \log p(\mathbf{x})] \qquad (24.175) $$

where the score function appears as

$$ \nabla \log p(\mathbf{x}) = \arg \sup_{\mathbf{y}} [\mathbf{y}^{\top} \mathbf{x} - \psi(\mathbf{y})] \qquad (24.176) $$

though this approach introduces additional computational complexity in high dimensions.

The thermodynamic integration perspective reveals deep connections to statistical physics, where the likelihood computation corresponds to computing the free energy difference between the noise and data distributions. The Jarzynski equality provides an identity

$$ \mathbb{E}[e^{-\beta W}] = e^{-\beta \Delta F} \qquad (24.177) $$

relating the work $W$ done along nonequilibrium paths to the equilibrium free energy difference

$$ \Delta F = - \log \frac{p_0(\mathbf{x}_0)}{p_T(\mathbf{x}_T)} \qquad (24.178) $$

with the score function playing the role of the driving force. This suggests alternative Monte Carlo estimators for the likelihood based on nonequilibrium work measurements, though their variance properties in high dimensions require careful analysis. The resulting likelihood computation framework provides not just a scalar probability value, but a complete characterization of the probability flow geometry connecting noise to data distributions.

### 24.3.10 Conclusion

Diffusion and score-based models provide a unified framework for generative modeling through SDEs and score matching. The forward process perturbs data, while the reverse process learns to denoise, enabling high-quality sample generation. The mathematical rigor stems from stochastic calculus, score matching, and dynamical systems.
# 24.4 Stable Diffusion

## 24.4.1 Literature Review of Stable Diffusion

The foundations of modern diffusion models were established by Sohl-Dickstein et al. (2015) [1404] through their pioneering work on nonequilibrium thermodynamics in deep learning, formulating the gradual noise corruption and denoising process as a Markov chain. This theoretical framework was operationalized by Ho et al. (2020) [1407] in their development of Denoising Diffusion Probabilistic Models (DDPM), which introduced practical training objectives and demonstrated viable image generation. The landmark Stable Diffusion model emerged from Rombach et al. (2022) [1429], who revolutionized the field by implementing diffusion in a compressed latent space using a VAE-GAN hybrid architecture, dramatically improving computational efficiency while maintaining high-quality outputs. Song et al. (2021) [1410] made fundamental theoretical contributions by unifying score-based models and diffusion processes through stochastic differential equations, while their subsequent work (Song et al. 2020 [1406]) on probability flow ODEs enabled more efficient sampling. Saharia et al. (2022) [1430] advanced text-to-image generation through novel conditioning techniques in Imagen, demonstrating effective integration with large language models. This direction was refined by Ramesh et al. (2022) [1431] in DALL-E 2, which implemented hierarchical diffusion for high-resolution generation from complex prompts.

Practical deployment was enabled by Nichol and Dhariwal (2021) [1408] through architectural improvements including learned noise schedules, while Dhariwal and Nichol (2021) [1409] proved diffusion models could surpass GANs in image quality. Ho and Salimans (2022) [1421] developed classifier-free guidance for more robust conditional generation. Sampling efficiency was dramatically improved by Song et al. (2020) [1416] through Denoising Diffusion Implicit Models (DDIM), with further optimizations from Karras et al. (2022) [1413] on sampling dynamics and Dockhorn et al. (2022) [1417] on Langevin dynamics. Conditioning mechanisms were advanced by Vahdat et al. (2021) [1412] through latent space manipulation techniques, while Bansal et al. (2022) [1414] generalized the framework beyond Gaussian noise. Ghosh (2024) [1415] presented the second part of a two-part series on deriving the analytical solution of the Burgers equation using the Cole-Hopf transformation. It employs the heat kernel method to obtain explicit solutions, linking nonlinear viscous flow equations to the linear heat equation. The work provides clear analytical insights into the behavior of systems modeled by the Burgers equation in contexts such as fluid dynamics and traffic flow. Text generation applications were developed by Li et al. (2022) [1419], with molecular design adaptations by Xu et al. (2022) [1427]. Robustness improvements came from Nie et al. (2022) [1422] on adversarial purification and Vincent et al. (2022) [1418] on uncertainty quantification.

Domain-specific adaptations have proliferated, including audio generation by Huang et al. (2022) [1411], structured data by Luo et al. (2022), and reinforcement learning by Janner et al. (2022) [1420]. Medical imaging applications were demonstrated by Chung et al. (2022) [1432], while Watson et al. (2022) [1435] applied diffusion models to protein design. Theoretical foundations were deepened by Block et al. (2022) [1423] on convergence properties and Liu et al. (2022) [1424] on optimal transport connections. Block et al. (2022) [1423] rigorously analyze the non-asymptotic convergence properties of diffusion models, providing theoretical guarantees on sampling quality and stability. Their work establishes conditions under which diffusion models provably converge, bridging gaps between empirical success and theoretical understanding. Liu et al. (2022) [1424] uncover deep connections between diffusion models and optimal transport theory, demonstrating how diffusion processes can be viewed as gradient flows in probability space. Their framework enables more principled training and sampling strategies, particularly in adversarial robustness settings. Recent advances include 3D-aware generation by Poole et al. (2023) [1436], video synthesis by
Harvey et al. (2023) [1437], and multimodal applications by Tewel et al. (2023) [1438]. Hybrid architectures have emerged through works like Xiao et al. (2023) [1439] on GAN-guided diffusion and Nash et al. (2023) [1440] on autoregressive-diffusion models. Fundamental improvements continue in discrete diffusion (Austin et al. 2023 [1441]), continuous-time formulations (Chen et al. 2023 [1442]), and geometric adaptations (De Bortoli et al. 2023 [1443]).

The field has also seen significant work on societal implications, including bias mitigation (Cho et al. 2023 [1444]), watermarking techniques (Fernandez et al. 2023 [1445]), and ethical frameworks (Gupta et al. 2023 [1446]). Memory optimizations were developed by Ryu et al. (2023) [1447], while dynamic thresholding was improved by Saharia et al. (2023) [1430]. Attention mechanisms were optimized by Zhang et al. (2023) [1456], and distillation techniques advanced by Liu et al. (2023) [1457]. The Zhang et al. (2023) [1456] introduces novel attention mechanism optimizations that significantly improve the efficiency and controllability of diffusion models, particularly in text-to-image generation tasks. Their work on conditional control mechanisms has been influential in enhancing the precision of Stable Diffusion outputs. Liu et al. (2023) [1457] makes important advances in distillation techniques for diffusion models, enabling more efficient sampling while maintaining image quality. Their semantic diffusion guidance approach has been particularly impactful for real-world applications where computational efficiency is crucial. These collective contributions have established Stable Diffusion as both a transformative technology and vibrant research area, with ongoing work pushing boundaries in efficiency (Kingma et al. 2023 [1453]), controllability (Meng et al. 2023 [1454]), and multimodal generation (Ramesh et al. 2023 [1455]). The integration with physical simulations (Sanchez-Gonzalez et al. 2023 [1448]) and symbolic reasoning (Ellis et al. 2023 [1449]) points toward increasingly sophisticated applications, while foundational work on model interpretability (Olah et al. 2023 [1450]) and safety (Amodei et al. 2023 [1451]) addresses critical deployment challenges.

### 24.4.2 Analysis of Stable Diffusion

Stable Diffusion operates in a latent space $\mathcal{Z}$ obtained via a pre-trained encoder $\mathcal{E}: \mathcal{X} \to \mathcal{Z}$, where $\mathcal{X}$ is the high-dimensional image space. The latent representation of an image $\mathbf{x}_0$ is given by:

$$ \mathbf{z}_0 = \mathcal{E}(\mathbf{x}_0), $$

where $\mathbf{z}_0 \in \mathbb{R}^{d_z}$ is a lower-dimensional embedding. The latent space is constrained to follow a prior distribution $p(\mathbf{z}_T) = \mathcal{N}(\mathbf{0}, \mathbf{I})$, achieved through variational autoencoding. The latent representation $\mathbf{z}_0$ is the starting point for the forward diffusion process, which gradually corrupts it with noise over time. The forward process gradually perturbs $\mathbf{z}_0$ with Gaussian noise over $T$ timesteps, governed by a variance schedule $\{\beta_t\}_{t=1}^T$. The transition kernel is:

$$ q(\mathbf{z}_t|\mathbf{z}_{t-1}) = \mathcal{N}(\mathbf{z}_t; \sqrt{1-\beta_t}\mathbf{z}_{t-1}, \beta_t\mathbf{I}). $$

The marginal distribution at time $t$ given $\mathbf{z}_0$ is:

$$ q(\mathbf{z}_t|\mathbf{z}_0) = \mathcal{N}(\mathbf{z}_t; \sqrt{\bar{\alpha}_t}\mathbf{z}_0, (1-\bar{\alpha}_t)\mathbf{I}), $$

where $\alpha_t = 1 - \beta_t$ and $\bar{\alpha}_t = \prod_{s=1}^t \alpha_s$. This allows direct sampling of $\mathbf{z}_t$ from $\mathbf{z}_0$ via:

$$ \mathbf{z}_t = \sqrt{\bar{\alpha}_t}\mathbf{z}_0 + \sqrt{1-\bar{\alpha}_t}\epsilon, \quad \epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I}). $$

The forward process is a discrete-time Markov chain with a fixed linear schedule for $\beta_t$. The forward process defines the noise corruption trajectory, which is reversed in the reverse diffusion process to generate samples. The noise addition can be interpreted as solving a stochastic differential equation (SDE) in continuous time:

$$ d\mathbf{z} = -\frac{1}{2}\beta(t)\mathbf{z} dt + \sqrt{\beta(t)} d\mathbf{W}, $$
where $\mathbf{W}$ is a Wiener process, and $\beta(t)$ is a continuous noise schedule. The discrete-time forward process is a first-order approximation of this SDE. The continuous-time perspective connects to the reverse diffusion process, which can also be framed as an SDE or ODE. The reverse process learns to denoise $\mathbf{z}_T \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ back to $\mathbf{z}_0$ by estimating:

$$p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_t, \mathbf{c}) = \mathcal{N}(\mathbf{z}_{t-1}; \mu_{\theta}(\mathbf{z}_t, t, \mathbf{c}), \Sigma_{\theta}(\mathbf{z}_t, t, \mathbf{c})),$$

where $\mathbf{c}$ is a conditioning vector (e.g., text embeddings). The mean $\mu_{\theta}$ is derived from the noise prediction network $\epsilon_{\theta}$:

$$\mu_{\theta}(\mathbf{z}_t, t, \mathbf{c}) = \frac{1}{\sqrt{\alpha_t}} \left( \mathbf{z}_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \right).$$

The covariance $\Sigma_{\theta}$ is often fixed to $\sigma_t^2 \mathbf{I}$, where:

$$\sigma_t^2 = \frac{1 - \bar{\alpha}_{t-1}}{1 - \bar{\alpha}_t} \beta_t.$$

The reverse process is trained via a training objective that minimizes the discrepancy between predicted and actual noise. The model is trained to minimize the evidence lower bound (ELBO):

$$\mathcal{L}(\theta) = \mathbb{E}_{t, \mathbf{z}_0, \epsilon, \mathbf{c}} \left[ \| \epsilon - \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \|_2^2 \right],$$

where $t \sim \mathcal{U}(\{1, \dots, T\})$, $\mathbf{z}_0 = \mathcal{E}(\mathbf{x}_0)$, $\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$, $\mathbf{z}_t = \sqrt{\bar{\alpha}_t} \mathbf{z}_0 + \sqrt{1 - \bar{\alpha}_t} \epsilon$. This is equivalent to denoising score matching in the latent space. The training objective ensures the noise predictor $\epsilon_{\theta}$ is accurate, enabling high-quality sampling via the Probability Flow ODE or discrete-time reversal. In the continuous-time limit, the reverse process can be described by an ordinary differential equation (ODE):

$$d\mathbf{z} = \left[ -\frac{1}{2}\beta(t)\mathbf{z} - \frac{1}{2}\beta(t)\nabla_{\mathbf{z}} \log p_t(\mathbf{z}) \right] dt,$$

where $\nabla_{\mathbf{z}} \log p_t(\mathbf{z})$ is the score function, approximated by:

$$\nabla_{\mathbf{z}} \log p_t(\mathbf{z}) \approx - \frac{\epsilon_{\theta}(\mathbf{z}, t, \mathbf{c})}{\sqrt{1 - \bar{\alpha}_t}}.$$

The ODE formulation provides a deterministic alternative to the discrete-time reverse sampling process. At inference, sampling proceeds iteratively from $t = T$ to $t = 1$:

$$\mathbf{z}_{t-1} = \frac{1}{\sqrt{\alpha_t}} \left( \mathbf{z}_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \right) + \sigma_t \xi,$$

where $\xi \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ and $\sigma_t$ controls stochasticity. The sampled $\mathbf{z}_0$ is decoded back to pixel space, requiring analysis of the latent space properties. The latent space $\mathcal{Z}$ must satisfy:

1. **Compactness**: $\mathcal{E}$ maps images to a bounded subspace of $\mathbb{R}^{d_z}$.

2. **Invertibility**: The decoder $\mathcal{D}$ should approximately satisfy $\mathcal{D}(\mathcal{E}(\mathbf{x})) \approx \mathbf{x}$.

3. **Gaussian Prior**: The final latent $\mathbf{z}_T$ converges to $\mathcal{N}(\mathbf{0}, \mathbf{I})$.

The diffusion process's mathematical behavior is further analyzed via diffusion kernel analysis. Regarding the Diffusion Kernel Analysis, the forward process can be analyzed via its transition kernels:

1. The Chapman-Kolmogorov equation ensures consistency:

$$q(\mathbf{z}_t|\mathbf{z}_0) = \int q(\mathbf{z}_t|\mathbf{z}_{t-1})q(\mathbf{z}_{t-1}|\mathbf{z}_0) d\mathbf{z}_{t-1}.$$
2. The Fokker-Planck equation describes the evolution of the probability density:

$$ \frac{\partial p_t(\mathbf{z})}{\partial t} = \nabla \cdot \left( \frac{1}{2} \beta(t) \mathbf{z} p_t(\mathbf{z}) \right) + \frac{1}{2} \beta(t) \nabla^2 p_t(\mathbf{z}). $$

For conditional generation, the noise predictor incorporates conditional diffusion via cross-attention. For text-to-image generation, the noise predictor $\epsilon_\theta$ is conditioned on embeddings $\mathbf{c} = \text{CLIP}(\mathbf{y})$, where $\mathbf{y}$ is the input text. The cross-attention mechanism in the U-Net ensures:

$$ \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) = \text{U-Net}(\mathbf{z}_t, t, \text{CrossAttn}(\mathbf{z}_t, \mathbf{c})). $$

The conditioned predictions guide the reverse process to generate samples aligned with the input text, completing the Stable Diffusion pipeline.

### 24.4.3 Latent Variable Model and Diffusion Process

Stable Diffusion is a latent diffusion model (LDM) that operates by gradually denoising a latent representation of an image. Let $\mathbf{x}_0$ be the true data sample (e.g., an image) in pixel space, and $\mathbf{z}_0 = \mathcal{E}(\mathbf{x}_0)$ its latent representation via an encoder $\mathcal{E}$. The forward diffusion process gradually adds Gaussian noise to $\mathbf{z}_0$ over $T$ timesteps, producing a sequence $\{\mathbf{z}_t\}_{t=1}^T$. The noise schedule is defined by variances $\{\beta_t\}_{t=1}^T$, where:

$$ q(\mathbf{z}_t|\mathbf{z}_{t-1}) = \mathcal{N}(\mathbf{z}_t; \sqrt{1-\beta_t}\mathbf{z}_{t-1}, \beta_t\mathbf{I}). $$

The marginal distribution at time $t$ is:

$$ q(\mathbf{z}_t|\mathbf{z}_0) = \mathcal{N}(\mathbf{z}_t; \sqrt{\bar{\alpha}_t}\mathbf{z}_0, (1-\bar{\alpha}_t)\mathbf{I}), $$

where $\alpha_t = 1 - \beta_t$ and $\bar{\alpha}_t = \prod_{s=1}^t \alpha_s$. This allows sampling $\mathbf{z}_t$ directly from $\mathbf{z}_0$:

$$ \mathbf{z}_t = \sqrt{\bar{\alpha}_t}\mathbf{z}_0 + \sqrt{1-\bar{\alpha}_t}\epsilon, \quad \epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I}). $$

Stable Diffusion is a generative model that combines latent variable modeling with a diffusion process to synthesize high-dimensional data, such as images, through a hierarchical probabilistic framework. The latent variable model operates on a compressed representation $\mathbf{z} \in \mathbb{R}^d$ of the input data $\mathbf{x} \in \mathbb{R}^D$ (where $d \ll D$), learned by a variational autoencoder (VAE) with encoder $q_\phi(\mathbf{z}|\mathbf{x})$ and decoder $p_\psi(\mathbf{x}|\mathbf{z})$. The encoder maps $\mathbf{x}$ to a Gaussian latent space:

$$ q_\phi(\mathbf{z}|\mathbf{x}) = \mathcal{N}(\mathbf{z}; \mu_\phi(\mathbf{x}), \Sigma_\phi(\mathbf{x})), $$

while the decoder reconstructs the data as $p_\psi(\mathbf{x}|\mathbf{z}) = \mathcal{N}(\mathbf{x}; \mu_\psi(\mathbf{z}), \sigma^2\mathbf{I})$. The VAE is trained by maximizing the evidence lower bound (ELBO):

$$ \mathcal{L}_{\text{VAE}}(\phi, \psi) = \mathbb{E}_{q_\phi(\mathbf{z}|\mathbf{x})}[\log p_\psi(\mathbf{x}|\mathbf{z})] - \beta D_{\text{KL}}(q_\phi(\mathbf{z}|\mathbf{x}) \parallel p_\psi(\mathbf{z})), $$

where $p_\psi(\mathbf{z}) = \mathcal{N}(\mathbf{0}, \mathbf{I})$ is the prior, and $\beta$ controls the trade-off between reconstruction fidelity and latent space regularization. The diffusion process is then applied in the latent space, defining a forward noising process $\{\mathbf{z}_t\}_{t=0}^T$ that gradually corrupts $\mathbf{z}_0$ with Gaussian noise:

$$ q(\mathbf{z}_t|\mathbf{z}_{t-1}) = \mathcal{N}(\mathbf{z}_t; \sqrt{1-\beta_t}\mathbf{z}_{t-1}, \beta_t\mathbf{I}), $$

where $\{\beta_t\}_{t=1}^T$ is a noise schedule. The marginal distribution at time $t$ is:

$$ q(\mathbf{z}_t|\mathbf{z}_0) = \mathcal{N}(\mathbf{z}_t; \sqrt{\bar{\alpha}_t}\mathbf{z}_0, (1-\bar{\alpha}_t)\mathbf{I}), $$
with $\bar{\alpha}_t = \prod_{s=1}^t (1 - \beta_s)$. The reverse process learns to denoise $\mathbf{z}_t$ via a neural network $\epsilon_\theta(\mathbf{z}_t, t, \mathbf{c})$, conditioned on a text embedding $\mathbf{c}$ from a CLIP model. The reverse transition is parameterized as:

$$p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_t, \mathbf{c}) = \mathcal{N}(\mathbf{z}_{t-1}; \mu_{\theta}(\mathbf{z}_t, t, \mathbf{c}), \Sigma_{\theta}(\mathbf{z}_t, t)),$$

where $\mu_\theta$ predicts the clean latent:

$$\mu_{\theta}(\mathbf{z}_t, t, \mathbf{c}) = \frac{1}{\sqrt{\bar{\alpha}_t}} \left( \mathbf{z}_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \right).$$

The training objective minimizes the denoising score matching loss:

$$\mathcal{L}_{\text{diffusion}}(\theta) = \mathbb{E}_{t, \mathbf{z}_0, \epsilon, \mathbf{c}} \left[ \|\epsilon - \epsilon_{\theta}(\sqrt{\bar{\alpha}_t} \mathbf{z}_0 + \sqrt{1 - \bar{\alpha}_t} \epsilon, t, \mathbf{c})\|^2 \right],$$

where $\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. Sampling synthesizes data by first generating $\mathbf{z}_T \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$, iteratively denoising it to $\mathbf{z}_0$ via the reverse process, and decoding $\mathbf{z}_0$ to pixel space with $p_\psi(\mathbf{x}|\mathbf{z}_0)$. The integration of latent variable modeling and diffusion enables efficient high-resolution generation by operating in a compressed, semantically meaningful space.

### 24.4.4 Reverse Diffusion and Denoising

The reverse process learns to iteratively denoise $\mathbf{z}_T \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ back to $\mathbf{z}_0$. A neural network $\epsilon_\theta$ predicts the noise at each step:

$$\epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \approx \epsilon,$$

where $\mathbf{c}$ is a conditioning vector (e.g., text embeddings). The reverse transition is:

$$p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_t, \mathbf{c}) = \mathcal{N}(\mathbf{z}_{t-1}; \mu_{\theta}(\mathbf{z}_t, t, \mathbf{c}), \Sigma_{\theta}(\mathbf{z}_t, t, \mathbf{c})),$$

where $\mu_\theta$ and $\Sigma_\theta$ are parameterized by the model. The mean is typically derived as:

$$\mu_{\theta}(\mathbf{z}_t, t, \mathbf{c}) = \frac{1}{\sqrt{\bar{\alpha}_t}} \left( \mathbf{z}_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \right).$$

The reverse diffusion process in Stable Diffusion constitutes a learned Markov chain that iteratively denoises latent representations to synthesize data samples from noise. Given a latent variable $\mathbf{z}_T \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ at diffusion time $T$, the reverse process incrementally refines $\mathbf{z}_t$ to $\mathbf{z}_{t-1}$ through a sequence of Gaussian transitions conditioned on a text embedding $\mathbf{c}$. Each reverse step is parameterized by a neural network $\epsilon_\theta(\mathbf{z}_t, t, \mathbf{c})$ that predicts the noise component added during the forward process. The reverse transition distribution is derived from the forward process posterior $q(\mathbf{z}_{t-1}|\mathbf{z}_t, \mathbf{z}_0)$ using the reparameterization trick:

$$p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_t, \mathbf{c}) = \mathcal{N}\left(\mathbf{z}_{t-1}; \frac{1}{\sqrt{\bar{\alpha}_t}}\left(\mathbf{z}_t - \frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}}\epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c})\right), \Sigma_{\theta}(\mathbf{z}_t, t)\right),$$

where $\alpha_t = 1 - \beta_t$, $\bar{\alpha}_t = \prod_{s=1}^t \alpha_s$, and $\Sigma_\theta$ is typically fixed to $\tilde{\beta}_t \mathbf{I}$ with $\tilde{\beta}_t = \frac{1-\bar{\alpha}_{t-1}}{1-\bar{\alpha}_t} \beta_t$. The mean $\mu_\theta$ is optimized to match the ground-truth denoising direction, which can be expressed in closed form using the forward process marginals:

$$\mu_{\theta}(\mathbf{z}_t, t, \mathbf{c}) = \frac{\sqrt{\bar{\alpha}_{t-1}}\beta_t}{1 - \bar{\alpha}_t}\mathbf{z}_0 + \frac{\sqrt{\bar{\alpha}_t}(1 - \bar{\alpha}_{t-1})}{1 - \bar{\alpha}_t}\mathbf{z}_t,$$

where $\mathbf{z}_0$ is inferred from $\mathbf{z}_t$ via the noise prediction network as $\mathbf{z}_0 = \frac{1}{\sqrt{\bar{\alpha}_t}}(\mathbf{z}_t - \sqrt{1 - \bar{\alpha}_t}\epsilon_\theta(\mathbf{z}_t, t, \mathbf{c}))$. The training objective minimizes the variational lower bound on the negative log-likelihood, which reduces to a weighted $L_2$ loss on the noise prediction:

$$\mathcal{L}(\theta) = \mathbb{E}_{\mathbf{z}_0, \epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I}), t, \mathbf{c}} \left[ \|\epsilon - \epsilon_{\theta}(\sqrt{\bar{\alpha}_t} \mathbf{z}_0 + \sqrt{1 - \bar{\alpha}_t} \epsilon, t, \mathbf{c})\|^2 \right].$$
During sampling, the reverse chain is initialized with $\mathbf{z}_T \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ and iteratively applies the learned transitions for $t = T, \dots, 1$:

$$ \mathbf{z}_{t-1} = \frac{1}{\sqrt{\alpha_t}} \left( \mathbf{z}_t - \frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \right) + \sqrt{\tilde{\beta}_t \zeta}, $$

where $\zeta \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ for $t > 1$ and $\zeta = \mathbf{0}$ for $t = 1$. The final latent $\mathbf{z}_0$ is decoded to image space via the VAE decoder $p_{\psi}(\mathbf{x}|\mathbf{z}_0)$. The reverse process thus inverts the diffusion by progressively removing noise conditioned on the textual prompt, where each step is guided by the score estimate $\nabla_{\mathbf{z}_t} \log p(\mathbf{z}_t|\mathbf{c}) \approx -\frac{1}{\sqrt{1-\bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c})$. This score-based formulation connects the denoising objective to implicit generative modeling of the data manifold in latent space.

### 24.4.5 Training Objective

The model is trained to minimize the variational lower bound (VLB) on the negative log-likelihood. The simplified objective is:

$$ \mathcal{L}(\theta) = \mathbb{E}_{t, \mathbf{z}_0, \epsilon, \mathbf{c}} [\|\epsilon - \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c})\|_2^2], $$

where $t \sim \mathcal{U}(1, T)$, $\mathbf{z}_0 = \mathcal{E}(\mathbf{x}_0)$, and $\epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$.

The training objective in Stable Diffusion is derived from variational principles that seek to maximize the likelihood of observed data under a latent variable model while simultaneously learning a denoising process in the latent space. The model consists of three key components: a variational autoencoder (VAE) that compresses images $\mathbf{x}$ into latent representations $\mathbf{z}_0$, a diffusion process that gradually adds noise to these latents, and a conditional denoising network $\epsilon_{\theta}$ that learns to reverse this noising process. The complete objective function combines the VAE's evidence lower bound (ELBO) with the diffusion model's denoising score matching loss.

For the VAE component, given an image $\mathbf{x}$ and its latent encoding $\mathbf{z}_0$, the ELBO is formulated as:

$$ \mathcal{L}_{\text{VAE}}(\phi, \psi) = \mathbb{E}_{q_{\phi}(\mathbf{z}_0|\mathbf{x})}[\log p_{\psi}(\mathbf{x}|\mathbf{z}_0)] - \beta D_{\text{KL}}(q_{\phi}(\mathbf{z}_0|\mathbf{x}) \| p(\mathbf{z}_0)), $$

where $q_{\phi}(\mathbf{z}_0|\mathbf{x})$ is the encoder's Gaussian posterior with parameters $\phi$, $p_{\psi}(\mathbf{x}|\mathbf{z}_0)$ is the decoder's reconstruction likelihood with parameters $\psi$, and $p(\mathbf{z}_0) = \mathcal{N}(\mathbf{0}, \mathbf{I})$ is the latent prior. The $\beta$-term controls the trade-off between reconstruction accuracy and latent space regularization.

The diffusion process operates in latent space, where the forward noising process is defined by a fixed Markov chain that gradually adds Gaussian noise to $\mathbf{z}_0$ over $T$ steps according to a variance schedule $\{\beta_t\}_{t=1}^T$:

$$ q(\mathbf{z}_t|\mathbf{z}_{t-1}) = \mathcal{N}(\mathbf{z}_t; \sqrt{1-\beta_t}\mathbf{z}_{t-1}, \beta_t\mathbf{I}). $$

The marginal distribution at any timestep $t$ can be expressed in closed form as:

$$ q(\mathbf{z}_t|\mathbf{z}_0) = \mathcal{N}(\mathbf{z}_t; \sqrt{\bar{\alpha}_t}\mathbf{z}_0, (1-\bar{\alpha}_t)\mathbf{I}), $$

where $\bar{\alpha}_t = \prod_{s=1}^t (1-\beta_s)$. The denoising network $\epsilon_{\theta}$ is trained to predict the noise $\epsilon$ added to $\mathbf{z}_0$ at each timestep $t$, conditioned on a text embedding $\mathbf{c}$. The training objective for the diffusion model is a weighted $L_2$ loss between the predicted and actual noise:

$$ \mathcal{L}_{\text{diffusion}}(\theta) = \mathbb{E}_{\mathbf{z}_0, \epsilon \sim \mathcal{N}(\mathbf{0}, \mathbf{I}), t, \mathbf{c}} [\|\epsilon - \epsilon_{\theta}(\sqrt{\bar{\alpha}_t}\mathbf{z}_0 + \sqrt{1-\bar{\alpha}_t}\epsilon, t, \mathbf{c})\|^2], $$

where $\mathbf{z}_t = \sqrt{\bar{\alpha}_t}\mathbf{z}_0 + \sqrt{1-\bar{\alpha}_t}\epsilon$ is the noised latent at timestep $t$. This objective can be interpreted as denoising score matching, where the network learns to estimate the score function $\nabla_{\mathbf{z}_t} \log q(\mathbf{z}_t)$ of the noisy data distribution. The conditioning on $\mathbf{c}$ is implemented via cross-attention layers in the U-Net architecture of $\epsilon_{\theta}$, allowing text to guide the denoising process.
The complete training objective for Stable Diffusion jointly optimizes the VAE and diffusion components, though in practice they are often trained separately. The diffusion loss dominates the final objective, as it is responsible for learning the conditional generative process. The noise prediction formulation provides stable gradients and has been shown to outperform alternative parameterizations of the reverse process mean. During sampling, the trained $\epsilon_{\theta}$ is used in an iterative denoising procedure to synthesize new latents $z_0$ from pure noise $z_T \sim N(0, I)$, which are then decoded to images by the VAE decoder. The mathematical rigor of this objective ensures that the model learns a principled approximation to the true data distribution while maintaining computational tractability through latent space diffusion.

### 24.4.6 Architecture: U-Net with Cross-Attention

The denoising network $\epsilon_{\theta}$ is a U-Net with residual blocks and cross-attention layers for conditioning. Let $\mathbf{h}$ denote intermediate feature maps. The cross-attention mechanism integrates text embeddings $\mathbf{c}$:

$$ \text{Attention}(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = \text{softmax} \left( \frac{\mathbf{Q}\mathbf{K}^T}{\sqrt{d_k}} \right) \mathbf{V}, $$

where $\mathbf{Q} = \mathbf{h}\mathbf{W}_Q$, $\mathbf{K} = \mathbf{c}\mathbf{W}_K$, $\mathbf{V} = \mathbf{c}\mathbf{W}_V$, and $d_k$ is the dimension of the key vectors.

The denoising network in Stable Diffusion employs a U-Net architecture augmented with cross-attention mechanisms to enable conditional generation guided by text embeddings. The U-Net operates on latent representations $\mathbf{z}_t \in \mathbb{R}^{h_t \times w_t \times c_t}$ at each timestep $t$ of the diffusion process, with $h_t \times w_t$ denoting the spatial dimensions and $c_t$ the channel depth. The network consists of a symmetric encoder-decoder structure with skip connections, where the encoder progressively downsamples the spatial resolution while increasing channel depth through strided convolutions, and the decoder mirrors this process with transposed convolutions. Each resolution level $l$ in the encoder and decoder contains residual blocks that transform feature maps $\mathbf{F}^{(l)} \in \mathbb{R}^{h_l \times w_l \times c_l}$ through:

$$ \mathbf{F}_{\text{out}}^{(l)} = \mathcal{W}_{2}^{(l)} \circ \text{SiLU} \circ \mathcal{W}_{1}^{(l)}(\mathbf{F}_{\text{in}}^{(l)}) + \mathbf{F}_{\text{in}}^{(l)}, $$

where $\mathcal{W}_1^{(l)}, \mathcal{W}_2^{(l)}$ are 3×3 convolutional layers, and SiLU denotes the Sigmoid Linear Unit activation. The residual blocks incorporate adaptive group normalization (AdaGN) to condition on the timestep $t$ and text embedding $\mathbf{c}$:

$$ \text{AdaGN}(\mathbf{F}, t, \mathbf{c}) = \gamma_{t, \mathbf{c}} \cdot \frac{\mathbf{F} - \mu(\mathbf{F})}{\sigma(\mathbf{F})} + \beta_{t, \mathbf{c}}, $$

with $\gamma_{t, \mathbf{c}}, \beta_{t, \mathbf{c}}$ generated by a shared MLP that projects the concatenated $[t, \mathbf{c}]$ into modulation parameters.

The cross-attention mechanism injects textual conditioning at multiple resolutions by computing attention between flattened spatial features $\mathbf{F} \in \mathbb{R}^{(h_l w_l) \times c_l}$ and text embeddings $\mathbf{c} \in \mathbb{R}^{n \times d_{\text{text}}}$ (where $n$ is the sequence length). The queries $\mathbf{Q} = \mathbf{F}\mathbf{W}_Q$ are derived from image features, while keys $\mathbf{K} = \mathbf{c}\mathbf{W}_K$ and values $\mathbf{V} = \mathbf{c}\mathbf{W}_V$ originate from text:

$$ \text{CrossAttn}(\mathbf{F}, \mathbf{c}) = \text{softmax} \left( \frac{\mathbf{Q}\mathbf{K}^T}{\sqrt{d_k}} \right) \mathbf{V}\mathbf{W}_O + \mathbf{F}, $$

where $\mathbf{W}_Q \in \mathbb{R}^{c_l \times d_k}$, $\mathbf{W}_K$, $\mathbf{W}_V \in \mathbb{R}^{d_{\text{text}} \times d_k}$, and $\mathbf{W}_O \in \mathbb{R}^{d_k \times c_l}$ are learned projections. The attention output is combined with residual connections and layer normalization. The U-Net's bottleneck layer processes features at the lowest resolution with self-attention blocks that capture long-range spatial dependencies:

$$ \text{SelfAttn}(\mathbf{F}) = \text{softmax} \left( \frac{(\mathbf{F}\mathbf{W}_Q)(\mathbf{F}\mathbf{W}_K)^T}{\sqrt{d_k}} \right) (\mathbf{F}\mathbf{W}_V)\mathbf{W}_O + \mathbf{F}. $$
This hierarchical architecture, combining convolutional inductive biases for local structure with attention mechanisms for global conditioning, enables precise spatial denoising while maintaining alignment with textual prompts. The network's output $\epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \in \mathbb{R}^{h \times w \times c}$ predicts the noise component to be subtracted during the reverse diffusion process. The total parameter count is optimized through channel multipliers that balance capacity and computational efficiency across resolutions, typically following a geometric progression $c_l = c_{\text{base}} \cdot m^l$ where $m$ is a multiplier (e.g., 2) and $c_{\text{base}}$ the initial channel dimension. The U-Net's symmetric structure ensures that spatial information lost during downsampling is recovered through skip connections, while the cross-attention layers enable fine-grained control over the generated content through differentiable interaction between visual and textual representations.

### 24.4.7 Latent Space and Autoencoder

The encoder $\mathcal{E}$ and decoder $\mathcal{D}$ are trained jointly to minimize:

$$ \mathcal{L}_{\text{AE}} = \mathbb{E}_{\mathbf{x}_0} \left[ \|\mathbf{x}_0 - \mathcal{D}(\mathcal{E}(\mathbf{x}_0))\|_2^2 + \mathcal{L}_{\text{KL}}(q(\mathbf{z}_0|\mathbf{x}_0)\|p(\mathbf{z}_0)) \right], $$

where $p(\mathbf{z}_0) = \mathcal{N}(\mathbf{0}, \mathbf{I})$ and $\mathcal{L}_{\text{KL}}$ is the Kullback-Leibler divergence.

The latent space in Stable Diffusion is a compressed, lower-dimensional representation of images learned through a variational autoencoder (VAE), which serves as a critical component for efficient diffusion modeling in high-dimensional pixel space. The autoencoder consists of an encoder $E_{\phi}: \mathcal{X} \to \mathcal{Z}$ that maps an input image $\mathbf{x} \in \mathbb{R}^{H \times W \times 3}$ to a latent variable $\mathbf{z}_0 \in \mathbb{R}^{h \times w \times c}$, where $(h, w) = (H/f, W/f)$ for a downsampling factor $f$, and a decoder $D_{\psi}: \mathcal{Z} \to \mathcal{X}$ that reconstructs the image from the latent code. The encoder's posterior distribution is given by:

$$ q_{\phi}(\mathbf{z}_0|\mathbf{x}) = \mathcal{N}(\mathbf{z}_0; \mu_{\phi}(\mathbf{x}), \sigma_{\phi}^2(\mathbf{x})\mathbf{I}), $$

where $\mu_{\phi}$ and $\sigma_{\phi}$ are learned functions parameterized by convolutional neural networks. The decoder defines a likelihood model for image reconstruction:

$$ p_{\psi}(\mathbf{x}|\mathbf{z}_0) = \mathcal{N}(\mathbf{x}; \mu_{\psi}(\mathbf{z}_0), \mathbf{I}), $$

with $\mu_{\psi}$ implemented as a symmetric convolutional network. The VAE is trained to maximize the evidence lower bound (ELBO):

$$ \mathcal{L}_{\text{VAE}}(\phi, \psi) = \mathbb{E}_{q_{\phi}(\mathbf{z}_0|\mathbf{x})} [\log p_{\psi}(\mathbf{x}|\mathbf{z}_0)] - \beta D_{\text{KL}}(q_{\phi}(\mathbf{z}_0|\mathbf{x}) \| p(\mathbf{z}_0)), $$

where $p(\mathbf{z}_0) = \mathcal{N}(\mathbf{0}, \mathbf{I})$ is the prior, and $\beta$ controls the trade-off between reconstruction fidelity and latent space regularization. The KL divergence term ensures the latent space is well-structured and amenable to the diffusion process.

The latent space $\mathcal{Z}$ provides a computationally tractable domain for the diffusion process, as operating directly in pixel space $\mathcal{X}$ would be prohibitively expensive for high-resolution images. The diffusion model is trained to denoise latent variables $\mathbf{z}_t$ rather than raw pixels, with the forward process defined by:

$$ q(\mathbf{z}_t|\mathbf{z}_{t-1}) = \mathcal{N}(\mathbf{z}_t; \sqrt{1-\beta_t}\mathbf{z}_{t-1}, \beta_t\mathbf{I}), $$

and the reverse process learned by a neural network $\epsilon_{\theta}$ that predicts the noise component:

$$ p_{\theta}(\mathbf{z}_{t-1}|\mathbf{z}_t) = \mathcal{N}(\mathbf{z}_{t-1}; \mu_{\theta}(\mathbf{z}_t, t), \Sigma_{\theta}(\mathbf{z}_t, t)), $$

where $\mu_{\theta}$ is derived from $\epsilon_{\theta}$ as:

$$ \mu_{\theta}(\mathbf{z}_t, t) = \frac{1}{\sqrt{\alpha_t}} \left( \mathbf{z}_t - \frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t) \right). $$
The autoencoder's decoder $D_\psi$ maps the final denoised latent $\mathbf{z}_0$ back to pixel space, completing the generative process. The compression factor $f$ is typically set to 8, reducing memory and computation costs by 64× while preserving perceptual quality, as the latent space captures high-level features and discards high-frequency details that can be plausibly reconstructed. The combination of the VAE's latent space with diffusion models enables Stable Diffusion to achieve high-quality synthesis at resolutions impractical for pixel-space diffusion, while maintaining the theoretical guarantees of variational inference and score-based generative modeling.

### 24.4.8 Sampling Process

Sampling generates images by iteratively applying the reverse process:

$$ \mathbf{z}_{t-1} = \frac{1}{\sqrt{\alpha_t}} \left( \mathbf{z}_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \right) + \sigma_t \xi, $$

where $\xi \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ and $\sigma_t$ controls stochasticity. The final image is decoded as $\mathbf{x}_0 = \mathcal{D}(\mathbf{z}_0)$.

The sampling process in Stable Diffusion generates images by iteratively denoising a latent variable through a learned reverse diffusion trajectory, conditioned on a text prompt. The process begins by sampling a noise vector $\mathbf{z}_T \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ in the latent space $\mathbb{R}^{h \times w \times c}$, where $h \times w$ is the spatial resolution of the compressed representation and $c$ the channel dimension. This noise vector is progressively refined through $T$ discrete steps of the reverse diffusion process, each governed by the learned denoising function $\epsilon_\theta(\mathbf{z}_t, t, \mathbf{c})$, where $\mathbf{c}$ is the text embedding from a CLIP text encoder.

The reverse transition from timestep $t$ to $t-1$ follows the discretized reverse-time stochastic differential equation (SDE):

$$ \mathbf{z}_{t-1} = \frac{1}{\sqrt{\alpha_t}} \left( \mathbf{z}_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}} \epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) \right) + \sigma_t \zeta, $$

where $\alpha_t = 1 - \beta_t$, $\bar{\alpha}_t = \prod_{s=1}^t \alpha_s$, $\sigma_t$ is the noise scale determined by the scheduler, and $\zeta \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ for stochastic sampling (or $\zeta = \mathbf{0}$ for deterministic sampling). The noise prediction network $\epsilon_\theta$ is implemented as a time-conditioned U-Net with cross-attention to $\mathbf{c}$, trained to estimate the score function $\nabla_{\mathbf{z}_t} \log p(\mathbf{z}_t|\mathbf{c}) \approx -\frac{1}{\sqrt{1-\bar{\alpha}_t}} \epsilon_\theta(\mathbf{z}_t, t, \mathbf{c})$.

For classifier-free guidance, the noise prediction is adjusted to amplify the influence of the text condition:

$$ \hat{\epsilon}_{\theta}(\mathbf{z}_t, t, \mathbf{c}) = \epsilon_{\theta}(\mathbf{z}_t, t, \emptyset) + s \cdot (\epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) - \epsilon_{\theta}(\mathbf{z}_t, t, \emptyset)), $$

where $s \ge 1$ is the guidance scale and $\emptyset$ denotes a null prompt. This convex combination shifts the sampling trajectory toward higher-density regions of the conditional distribution $p(\mathbf{z}_0|\mathbf{c})$. After $T$ denoising steps, the final latent $\mathbf{z}_0$ is decoded to pixel space by the VAE decoder $D_\psi$:

$$ \mathbf{x} = D_{\psi}(\mathbf{z}_0) + \xi, $$

where $\xi \sim \mathcal{N}(\mathbf{0}, \sigma^2\mathbf{I})$ represents small stochastic variations, though in practice $\sigma \approx 0$ for deterministic decoding. The sampling process thus combines iterative latent-space refinement with precise text conditioning, where each step reduces the uncertainty in $\mathbf{z}_t$ while preserving the semantic structure imposed by $\mathbf{c}$. The mathematical rigor of this procedure ensures convergence to samples from the approximate data distribution $p(\mathbf{x}|\mathbf{c})$ while maintaining computational feasibility through latent-space operations and parallelizable denoising steps.
## 24.4.9 Classifier-Free Guidance

To improve sample quality, classifier-free guidance interpolates between conditional and unconditional predictions:

$$ \hat{\epsilon}_{\theta}(\mathbf{z}_t, t, \mathbf{c}) = \epsilon_{\theta}(\mathbf{z}_t, t, \emptyset) + w \cdot (\epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) - \epsilon_{\theta}(\mathbf{z}_t, t, \emptyset)), $$

where $w$ is the guidance scale and $\emptyset$ denotes a null condition.

Classifier-free guidance is a technique for enhancing the conditioning fidelity of diffusion models without relying on explicit classifier gradients, achieved through an implicit interpolation between conditional and unconditional score estimates. Given a conditional diffusion model $\epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c})$ trained on text-image pairs $(\mathbf{z}, \mathbf{c})$, where $\mathbf{c}$ is a text embedding and $\mathbf{z}_t$ the noisy latent at timestep $t$, the model is jointly trained to also predict an unconditional score $\epsilon_{\theta}(\mathbf{z}_t, t, \emptyset)$ by randomly dropping the conditioning $\mathbf{c}$ during training (typically with 10-20 percentage probability). The guided noise prediction $\hat{\epsilon}_{\theta}$ is computed as a linear combination:

$$ \hat{\epsilon}_{\theta}(\mathbf{z}_t, t, \mathbf{c}) = \epsilon_{\theta}(\mathbf{z}_t, t, \emptyset) + s \cdot (\epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) - \epsilon_{\theta}(\mathbf{z}_t, t, \emptyset)), $$

where $s \ge 1$ is the guidance scale. This formulation can be derived by considering the gradient of the log-probability of the joint distribution $p(\mathbf{z}_t, \mathbf{c}) = p(\mathbf{z}_t|\mathbf{c})p(\mathbf{c})$. Under the approximation that the unconditional model estimates $\nabla_{\mathbf{z}_t} \log p(\mathbf{z}_t)$ and the conditional model estimates $\nabla_{\mathbf{z}_t} \log p(\mathbf{z}_t|\mathbf{c})$, the guidance term $\nabla_{\mathbf{z}_t} \log p(\mathbf{c}|\mathbf{z}_t)$ is approximated via Bayes' rule as:

$$ \nabla_{\mathbf{z}_t} \log p(\mathbf{c}|\mathbf{z}_t) = \nabla_{\mathbf{z}_t} \log p(\mathbf{z}_t|\mathbf{c}) - \nabla_{\mathbf{z}_t} \log p(\mathbf{z}_t), $$

which corresponds to the difference $\epsilon_{\theta}(\mathbf{z}_t, t, \mathbf{c}) - \epsilon_{\theta}(\mathbf{z}_t, t, \emptyset)$. The guided update thus amplifies the influence of the conditioning signal $\mathbf{c}$ by scaling this difference term, effectively shifting the sampling trajectory toward regions of latent space where $p(\mathbf{c}|\mathbf{z}_t)$ is maximized. The guidance scale $s$ controls a trade-off between sample quality (lower $s$) and text alignment (higher $s$), with typical values $s \in [7.5, 10]$ for Stable Diffusion.

Theoretical analysis reveals that classifier-free guidance modifies the implicit energy function being optimized during sampling. The guided score estimate corresponds to the gradient of a modified log-probability:

$$ \nabla_{\mathbf{z}_t} \log \tilde{p}(\mathbf{z}_t|\mathbf{c}) = \nabla_{\mathbf{z}_t} \log p(\mathbf{z}_t) + s \cdot \nabla_{\mathbf{z}_t} \log p(\mathbf{c}|\mathbf{z}_t), $$

which can be interpreted as sharpening the conditional distribution by a factor of $s$. This approach avoids the computational overhead and instability of classifier-based guidance methods, while maintaining differentiability and enabling fine-grained control over the conditioning strength. The success of classifier-free guidance hinges on the shared parameterization of the conditional and unconditional models, which ensures that their predictions reside in the same vector space and can be meaningfully interpolated. The technique is particularly effective in text-to-image generation, where it significantly improves the semantic alignment between generated images and text prompts without requiring auxiliary classifiers or explicit likelihood estimation.
## 25 Kernel Regression

Kernel regression is a non-parametric method for estimating the conditional expectation of a target variable given input features, leveraging the representational capacity of reproducing kernel Hilbert spaces (RKHS). In the context of deep learning, kernel methods are often employed to analyze and interpret the behavior of neural networks, particularly in the infinite-width limit where they converge to kernel machines. Given a dataset $\{(\mathbf{x}_i, y_i)\}_{i=1}^N$ with $\mathbf{x}_i \in \mathbb{R}^d$ and $y_i \in \mathbb{R}$, the kernel regression estimator for a new input $\mathbf{x}$ is given by:

$$ \hat{y}(\mathbf{x}) = \sum_{i=1}^{N} \alpha_i K(\mathbf{x}, \mathbf{x}_i), $$

where $K: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ is a symmetric positive-definite kernel function, and the coefficients $\alpha = (\alpha_1, \dots, \alpha_N)^T$ are obtained by solving the linear system:

$$ \alpha = (\mathbf{K} + \lambda \mathbf{I})^{-1} \mathbf{y}, $$

with $\mathbf{K}_{ij} = K(\mathbf{x}_i, \mathbf{x}_j)$ the kernel matrix, $\mathbf{y} = (y_1, \dots, y_N)^T$, and $\lambda \ge 0$ a regularization parameter. The kernel $K$ implicitly defines a feature map $\phi: \mathbb{R}^d \to \mathcal{H}$ into an RKHS $\mathcal{H}$, where $K(\mathbf{x}, \mathbf{x}') = \langle \phi(\mathbf{x}), \phi(\mathbf{x}') \rangle_{\mathcal{H}}$. Common choices include the Gaussian radial basis function (RBF) kernel:

$$ K(\mathbf{x}, \mathbf{x}') = \exp \left( - \frac{\|\mathbf{x} - \mathbf{x}'\|^2}{2\sigma^2} \right), $$

and the polynomial kernel $K(\mathbf{x}, \mathbf{x}') = (\mathbf{x}^T \mathbf{x}' + c)^m$.

The connection between kernel regression and deep learning is formalized through the neural tangent kernel (NTK), which characterizes the dynamics of infinitely wide neural networks during gradient descent training. For a neural network $f(\mathbf{x}; \boldsymbol{\theta})$ with parameters $\boldsymbol{\theta}$, the NTK $\Theta(\mathbf{x}, \mathbf{x}')$ is defined as:

$$ \Theta(\mathbf{x}, \mathbf{x}') = \mathbb{E}_{\boldsymbol{\theta} \sim \mathcal{P}} [\langle \nabla_{\boldsymbol{\theta}} f(\mathbf{x}; \boldsymbol{\theta}), \nabla_{\boldsymbol{\theta}} f(\mathbf{x}'; \boldsymbol{\theta}) \rangle], $$

where the expectation is taken over the parameter initialization distribution $\mathcal{P}$. In the infinite-width limit, the NTK remains constant during training, and the network's output evolves according to the kernel gradient flow:

$$ \frac{\partial f(\mathbf{x}; \boldsymbol{\theta}(t))}{\partial t} = - \sum_{i=1}^{N} \Theta(\mathbf{x}, \mathbf{x}_i) (f(\mathbf{x}_i; \boldsymbol{\theta}(t)) - y_i). $$

The solution to this differential equation yields the deep kernel regression predictor:

$$ f(\mathbf{x}; \boldsymbol{\theta}(t)) = \mathbf{K}_{\text{NTK}}(\mathbf{x}, \mathbf{X}) \mathbf{K}_{\text{NTK}}(\mathbf{X}, \mathbf{X})^{-1} (\mathbf{I} - e^{-\eta \mathbf{K}_{\text{NTK}}(\mathbf{X}, \mathbf{X})t}) \mathbf{y}, $$

where $\mathbf{K}_{\text{NTK}}(\mathbf{X}, \mathbf{X})$ is the NTK matrix on the training data, and $\eta$ is the learning rate. This establishes that wide neural networks trained with gradient descent converge to kernel regression estimators with the NTK as their kernel function.
The generalization performance of kernel regression is governed by the eigendecay of the kernel operator. For a kernel K with eigenvalues {$\lambda_k$}$_{k=1}^\infty$ satisfying $\lambda_k \asymp k^{-2\beta}$, the minimax optimal rate for the excess risk is:

$$ \mathbb{E} [\|\hat{y} - y\|^2] - \inf_{f \in \mathcal{H}} \mathbb{E} [\|f(\mathbf{x}) - y\|^2] \lesssim N^{-\frac{2\beta}{2\beta+d}}, $$

where d is the input dimension. This highlights the curse of dimensionality, which deep neural networks mitigate through implicit regularization and hierarchical feature learning. The NTK's spectral properties determine the network's inductive bias: kernels with slower eigendecay favor smoother functions, while deeper architectures induce faster eigendecay, promoting hierarchical representations. The interplay between kernel design, architecture depth, and training dynamics forms the theoretical foundation for understanding deep learning through the lens of kernel methods.

## 25.1 Literature Review

Fan et. al. (2025) [696] explored kernel regression techniques in causal inference, particularly in the presence of interference among observations. The authors propose an innovative nonparametric estimator that integrates kernel regression with trimming methods, improving robustness in observational studies. Atanasov et. al. (2025) [697] generalized kernel regression by linking it to high-dimensional linear models and stochastic gradient dynamics. The authors present new asymptotics that extend classical results in nonparametric regression and random feature models. Mishra et. al. (2025) [699] applied Gaussian kernel-based regression to image classification and feature extraction. The authors demonstrate how kernel selection significantly impacts model performance in plant leaf detection tasks. Elsayed and Nazier (2025) [700] combined kernel smoothing regression with decomposition analysis to study labor market trends. It highlights the application of kernel-based regression techniques in socio-economic modeling. Kong et. al. (2025) [701] applied Bayesian Kernel Machine Regression (BKMR) to analyze complex relationships between heavy metal exposure and health indicators. It extends kernel regression to toxicology and epidemiological studies. Bracale et. al. (2025) [702] explored antitonic regression methods, establishing new concentration inequalities for regression problems. It highlights kernel methods' superiority over traditional parametric approaches in pricing models. Köhne et. al. (2025) [703] provided a theoretical foundation for kernel regression within Hilbert spaces, focusing on error bounds for kernel approximations in dynamical systems. Sadeghi and Beyeler (2025) [704] applied Gaussian Process Regression (GPR) with a Matérn kernel to estimate perceptual thresholds in retinal implants, showcasing kernel-based regression in biomedical engineering. Naresh et. al. (2025) [705] in a book chapter discussed logistic regression and kernel methods in network security. It illustrates how kernelized models can enhance cybersecurity measures in firewalls. Zhao et. al. (2025) [706] proposed Deep Additive Kernel (DAK) models, which unify kernel methods with deep learning. This approach enhances Bayesian neural networks' interpretability and robustness.

## 25.2 Analysis of Kernel regression

Kernel regression is a non-parametric statistical learning technique that estimates an unknown function $f(x)$ based on a given dataset:

$$ \{(x_i, y_i)\}_{i=1}^n, \qquad (25.1) $$

where $x_i \in \mathbb{R}^d$ and $y_i \in \mathbb{R}$. The fundamental kernel regression estimator is given by:

$$ \hat{f}(x) = \sum_{i=1}^{n} \alpha_i K(x, x_i), \qquad (25.2) $$
where $K(x, x')$ is a positive definite kernel function, ensuring that the Gram matrix

$$K = [K(x_i, x_j)]_{i,j=1}^{n} \qquad (25.3)$$

is symmetric positive semi-definite (PSD). The spectral properties of $K$ are crucial for understanding kernel regression's behavior, particularly in the context of regularization, overfitting, and generalization error analysis. To rigorously analyze kernel regression, we consider the Reproducing Kernel Hilbert Space (RKHS) $\mathcal{H}_K$ induced by $K(x, x')$, where functions satisfy:

$$f(x) = \sum_{i=1}^{\infty} \alpha_i \varphi_i(x), \qquad (25.4)$$

where $\varphi_i(x)$ are the eigenfunctions of the integral operator associated with $K(x, x')$:

$$Kf(x) = \int_{\Omega} K(x, x')f(x')d\mu(x'). \qquad (25.5)$$

The spectral decomposition of the Kernel Function $K$ takes the form:

$$K\varphi_i = \lambda_i\varphi_i, \quad i = 1, 2, \dots \qquad (25.6)$$

where

$$\lambda_1 \ge \lambda_2 \ge \dots \ge 0 \qquad (25.7)$$

are the eigenvalues of $K$. These eigenvalues and eigenfunctions determine the approximation capacity of kernel regression and its regularization properties.

# 25.3 Nadaraya–Watson kernel estimator

## 25.3.1 Literature Review

Agua and Bouzebda (2024) [685] explored the Nadaraya–Watson estimator for locally stationary functional time series. It presents asymptotic properties of kernel regression estimators in functional settings, emphasizing how they behave in nonstationary time series. Bouzebda et. al. (2024) [686] generalized Nadaraya–Watson estimators using asymmetric kernels. It rigorously analyzes the Dirichlet kernel estimator and provides first theoretical justifications for its application in conditional U-statistics. Zhao et. al. (2025) [687] applied Nadaraya–Watson regression in engineering applications, specifically in high-voltage circuit breaker degradation modeling. The method smooths interpolated datasets to eliminate measurement errors. Patil et. al. (2024) [688] addressed the bias-variance tradeoff in Nadaraya–Watson kernel regression, showing how optimal smoothing can improve signal denoising and estimation accuracy in noisy environments. Kakani and Radhika (2024) [689] evaluated Nadaraya–Watson estimation in medical data analysis, comparing it with regression trees and other machine learning methods. It highlights the role of bandwidth selection in clinical prediction tasks. Kato (2024) [690] presented a debiased version of Nadaraya–Watson regression, improving its root-N consistency and performance in conditional mean estimation. Sadek and Mohammed (2024) [691] did a comparative study of kernel-based Nadaraya–Watson regression and ordinary least squares (OLS), showing scenarios where nonparametric regression outperforms classical regression techniques. Gong et. al. (2024) [693] introduced Kernel-Thinned Nadaraya–Watson Estimator (KT-NW), which reduces computational cost while maintaining accuracy. This work is highly relevant for large-scale machine learning applications. Zavatone-Veth and Pehlevan (2025) [694] established a theoretical link between Nadaraya–Watson kernel smoothing and statistical physics through the random energy model. It offers new perspectives on kernel regression in high-dimensional settings. Ferrigno (2024) [695] explored how Nadaraya–Watson kernel regression can be applied to reference curve estimation, a key technique in medical statistics and economic forecasting.
### 25.3.2 Analysis of Nadaraya-Watson kernel estimator

Kernel regression is a **non-parametric regression technique** that estimates a function $f(x)$ using a weighted sum of observed values $y_i$. Given training data $\{(x_i, y_i)\}_{i=1}^n$, where $x_i \in \mathbb{R}^d$ and $y_i \in \mathbb{R}$, the Nadaraya-Watson kernel regression estimator takes the form

$$ \hat{f}(x) = \frac{\sum_{i=1}^{n} K_h(x - x_i)y_i}{\sum_{i=1}^{n} K_h(x - x_i)} \quad (25.8) $$

where $K_h(x)$ is the **scaled kernel function** defined as

$$ K_h(x) = \frac{1}{h^d} K\left(\frac{x}{h}\right), \quad (25.9) $$

where $h$ is the bandwidth parameter that determines the smoothing level. A crucial property of kernel functions is their **normalization condition**,

$$ \int_{\mathbb{R}^d} K(x) dx = 1. \quad (25.10) $$

A common choice for $K(x)$ is the Gaussian kernel:

$$ K(x) = \frac{1}{(2\pi)^{d/2}} e^{-\frac{1}{2}\|x\|^2}. \quad (25.11) $$

Let us now do the Bias-Variance Decomposition and Overfitting in Kernel Regression. The performance of kernel regression is governed by the **bias-variance tradeoff**:

$$ \mathbb{E}[(\hat{f}(x) - f(x))^2] = \text{Bias}^2 + \text{Variance} + \sigma_{\text{noise}}^2. \quad (25.12) $$

where

$$ \text{Bias}(\hat{f}(x)) = \mathbb{E}[\hat{f}(x)] - f(x), \quad (25.13) $$

and

$$ \text{Var}(\hat{f}(x)) = \mathbb{E}[(\hat{f}(x) - \mathbb{E}[\hat{f}(x)])^2]. \quad (25.14) $$

Expanding $f(x)$ via **Taylor series**, we obtain

$$ f(x_i) \approx f(x) + (x_i - x)^T \nabla f(x) + \frac{1}{2} (x_i - x)^T H_f(x) (x_i - x). \quad (25.15) $$

The expectation of the kernel estimate gives

$$ \mathbb{E}[\hat{f}(x)] = f(x) + \frac{h^2}{2} \sum_{j=1}^{d} \left( \frac{\int u_j^2 K(u) du}{\int K(u) du} \right) \frac{\partial^2 f}{\partial x_j^2} + O(h^4), \quad (25.16) $$

showing **bias scales as $O(h^2)$**. The variance analysis yields

$$ \text{Var}(\hat{f}(x)) \approx \frac{\sigma^2}{nh^d f_X(x)} \int K^2(u) du. \quad (25.17) $$

Thus, variance scales as $O((nh^d)^{-1})$, leading to the **optimal bandwidth selection**

$$ h^* \propto n^{-\frac{1}{d+4}}. \quad (25.18) $$

However, when $h$ is too small, overfitting occurs, characterized by high variance:

$$ \text{Var}(\hat{f}(x)) \gg 0, \quad \text{Bias}^2(\hat{f}(x)) \approx 0. \quad (25.19) $$
Kernel Ridge Regression (KRR) is one of the best Regularization Techniques to Prevent Overfitting. To control overfitting, we introduce **Tikhonov regularization** in kernel space. Define the **Gram matrix K** with entries

$$K_{ij} = K_h(x_i - x_j). \qquad (25.20)$$

We solve the regularized least squares problem:

$$\alpha = (\mathbf{K} + \lambda\mathbf{I})^{-1}\mathbf{y}. \qquad (25.21)$$

The regularization term $\lambda$ modifies the eigenvalues $\sigma_i$ of $\mathbf{K}$, giving

$$\alpha_i = \frac{\sigma_i}{\sigma_i + \lambda} v_i^T \mathbf{y}. \qquad (25.22)$$

For small $\lambda$, inverse eigenvalues $\sigma_i^{-1}$ amplify noise, whereas for large $\lambda$, the regularization term suppresses high-frequency components. In the **spectral decomposition** of $\mathbf{K}$, we write

$$\mathbf{K} = \sum_i \sigma_i v_i v_i^T. \qquad (25.23)$$

where $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_n \ge 0$ are the eigenvalues of the kernel matrix $\mathbf{K}$ and $v_i$ are the orthonormal eigenvectors, i.e.,

$$v_i^T v_j = \delta_{ij} \qquad (25.24)$$

where $\delta_{ij}$ is the Kronecker delta. The rank of $\mathbf{K}$ is equal to the number of nonzero eigenvalues $\sigma_i$. The eigenvalues of $\mathbf{K}$ encode the spectrum of feature space correlations. If the kernel function $K(x, x')$ is smooth, the eigenvalues decay exponentially: regularization, the solution is

$$\sigma_i \approx O(i^{-\tau}) \qquad (25.25)$$

for some decay exponent $\tau > 0$. The spectral decay controls the effective degrees of freedom of kernel regression. Applying regularization, the solution is

$$\hat{f}(x) = \sum_{i=1}^{n} \frac{\sigma_i}{\sigma_i + \lambda} v_i^T \mathbf{y} \cdot v_i(x). \qquad (25.26)$$

The regularization smoothly filters the high-frequency components of $f(x)$, preventing overfitting. For Controlling Model Complexity in Spectral Filtering, we have to note that large $\sigma_i$ corresponds to low-frequency components retained in the solution while small $\sigma_i$ are high-frequency components, attenuated by regularization. The cutoff occurs around defining the effective model complexity. For very small $\lambda$,

$$\frac{\sigma_i}{\sigma_i + \lambda} \approx 1, \qquad (25.27)$$

causing high variance. For large $\lambda$,

$$\frac{\sigma_i}{\sigma_i + \lambda} \approx \frac{\sigma_i}{\lambda}, \qquad (25.28)$$

which heavily suppresses small eigenvalues, leading to underfitting. The **optimal** $\lambda$ is selected via **cross-validation**, minimizing

$$CV(h) = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{f}_{-i}(x_i))^2. \qquad (25.29)$$

An alternative approach is smoothing Splines in Kernel Regression which is done by minimizing

$$\sum_{i=1}^{n} (y_i - f(x_i))^2 + \lambda \| \mathcal{L} f \|^2, \qquad (25.30)$$
where $\mathcal{L}$ is a differential operator like

$$ \mathcal{L} = \frac{d^2}{dx^2}. \qquad (25.31) $$

This results in the **smoothing spline estimator**

$$ \hat{f}(x) = \sum_{i=1}^{n} \alpha_i K_h(x - x_i). \qquad (25.32) $$

where $\alpha$ now depends on $K$ and $\mathcal{L}$. In conclusion, Kernel regression is powerful but prone to overfitting when $h$ is too small, leading to high variance. Regularization techniques such as **kernel ridge regression**, **Tikhonov regularization**, and **smoothing splines** mitigate overfitting by **modifying the spectral properties of the kernel matrix**.

There is a Bias-Variance Tradeoff in Spectral Terms. The expected bias measures the deviation of $\hat{f}(x)$ from $f(x)$:

$$ \text{Bias}^2 = \sum_{i=1}^{n} (1 - g(\sigma_i))^2 c_i^2. \qquad (25.33) $$

where

$$ g(\sigma_i) = \frac{\sigma_i}{\sigma_i + \lambda}. \qquad (25.34) $$

Large $\lambda$ shrinks eigenmodes, increasing bias. The variance measures sensitivity to noise:

$$ \text{Var} = \sigma^2 \sum_{i=1}^{n} g(\sigma_i)^2. \qquad (25.35) $$

For small $\lambda$, the model overfits, leading to high variance. The expected generalization bound error in Spectral Terms is given by:

$$ E[(\hat{f}(x) - f(x))^2] = \sum_{i} (1 - g(\sigma_i))^2 c_i^2 + \sigma^2 \sum_{i} g(\sigma_i)^2. \qquad (25.36) $$

Using asymptotic analysis, the optimal choice of $\lambda$ is:

$$ \lambda^* = O(n^{-\frac{1}{d+4}}). \qquad (25.37) $$

which minimizes error and maximizes generalization.

In conclusion, Spectral Properties play an important role in Kernel Regression. The spectral properties of kernel regression determine its ability to generalize and avoid overfitting:

* The eigenvalue decay rate controls approximation power.

* Spectral filtering via regularization prevents high-frequency noise.

* Generalization is optimized when balancing bias and variance.

By leveraging spectral decomposition, we gain a deep understanding of how kernel regression interpolates data while controlling complexity. The optimal choice of $\lambda$ and $h$ ensures an **optimal tradeoff between bias and variance**, leading to a robust kernel regression model.
# 25.4 Priestley-Chao kernel estimator

## 25.4.1 Literature Review

Neumann and Thorarinsdottir (2006) [677] discussed improvements on the Priestley-Chao estimator by addressing its limitations in nonparametric regression. It provides an asymptotic minimax framework for better estimation, particularly in autoregressive models. The modification proposed mitigates the issues arising from bandwidth selection. Steland (2014) [678] applied the Priestley-Chao kernel estimator to stopping rules in time series control charts. This study is significant because it explores its efficiency in dependent data, particularly focusing on bandwidth choice formulas that enhance estimation precision in practical scenarios. Makkulau et. al. (2023) [679] applied the Priestley-Chao estimator in multivariate semiparametric regression. It highlights the estimator's dependence on optimal bandwidth selection and explores modifications that enhance its adaptability in multiple dimensions. Staniswalis (1989) [680] examined the likelihood-based interpretation of kernel estimators and connects the Priestley-Chao approach with generalized maximum likelihood estimation. It rigorously analyzes the estimator's weighting properties and neighborhood selection criteria. Jennen-Steinmetz and Gasser (1988) [675] provided a comparative framework between the Priestley-Chao estimator and other kernel regression estimators. They explore its mathematical properties, convergence rates, and advantages over alternative methods such as Nadaraya-Watson. Mack and Müller (1988) [681] evaluated the Priestley-Chao estimator's error behavior in nonparametric regression. The paper highlights how convolution-type adjustments can improve estimation accuracy under random design conditions. Jones et. al. (2024) [682] categorized various kernel regression estimators, including the Priestley-Chao estimator. It critically evaluates its statistical efficiency and variance properties in comparison to other kernel methods. Ghosh (2015) [683] introduced a variance estimation technique specifically for the Priestley-Chao kernel estimator. The paper presents a method to avoid nuisance parameter estimation, improving computational efficiency. Liu and Luor (2023) [684] integrated fractal interpolants with the Priestley-Chao estimator to handle complex regression problems. It explores modifications to kernel functions that enhance estimation in high-dimensional datasets. Gasser and Müller (1979) [667] wrote a foundational work that revisits the Priestley-Chao estimator in the context of kernel regression. The authors propose two alternative definitions for kernel estimation, aiming to refine the estimator's application in empirical data analysis.

## 25.4.2 Analysis of Priestley-Chao kernel estimator

Let $X_1, X_2, \dots, X_n$ be independent and identically distributed (i.i.d.) random variables drawn from an unknown probability density function (PDF) $f(x)$, with cumulative distribution function (CDF) $F(x)$. The goal of nonparametric density estimation is to construct an estimator $\hat{f}(x)$ such that

$$ \lim_{n \to \infty} \hat{f}(x) = f(x) \qquad (25.38) $$

in some suitable sense, such as pointwise convergence, mean squared error (MSE) consistency, or uniform convergence over compact subsets. In kernel density estimation (KDE), a common approach is to define

$$ \hat{f}(x) = \frac{1}{nh} \sum_{i=1}^{n} K \left( \frac{x - X_i}{h} \right) \qquad (25.39) $$

where $K(\cdot)$ is a kernel function satisfying

$$ \int_{-\infty}^{\infty} K(u) du = 1 \qquad (25.40) $$
and $h$ is a bandwidth parameter controlling the level of smoothing. However, KDE relies on a fixed bandwidth $h$, which can lead to oversmoothing in regions of high density and undersmoothing in regions of low density. The Priestley-Chao estimator improves upon this by adapting the bandwidth locally, based on the spacings between consecutive order statistics.

There is an important role of Order Statistics and Spacings. Let us define the order statistics of the sample as

$$X_{(1)} \le X_{(2)} \le \dots \le X_{(n)} \qquad (25.41)$$

The fundamental insight behind the Priestley-Chao estimator is that the spacings between order statistics contain direct information about the local density. Define the spacing between two consecutive order statistics as

$$D_i = X_{(i+1)} - X_{(i)}, \quad i = 1, \dots, n-1 \qquad (25.42)$$

Using results from order statistics theory, we obtain the key approximation

$$\mathbb{E}[D_i] \approx \frac{1}{nf(X_{(i)})} \qquad (25.43)$$

which follows from the fact that the probability of observing a sample in a small interval around $X_{(i)}$ is approximately given by the density $f(X_{(i)})$ times the width of the interval. Thus, rearranging, we obtain the fundamental estimator

$$\hat{f}(X_{(i)}) \approx \frac{1}{nD_i} \qquad (25.44)$$

This provides a direct data-driven way to estimate the density without choosing a fixed bandwidth $h$, as in classical KDE methods. Let's now state the formal Definition of the Priestley-Chao Estimator. The Priestley-Chao kernel estimator is defined as

$$\hat{f}(x) = \frac{1}{n} \sum_{i=1}^{n-1} \frac{1}{D_i} K \left( \frac{x - X_{(i)}}{D_i} \right) \qquad (25.45)$$

where $K(\cdot)$ is a symmetric kernel function satisfying

$$\int_{-\infty}^{\infty} K(u) du = 1, \quad \int_{-\infty}^{\infty} uK(u) du = 0, \quad \text{and} \quad \int_{-\infty}^{\infty} u^2 K(u) du < \infty \qquad (25.46)$$

Unlike fixed-bandwidth KDE, here the bandwidth $D_i$ varies with location, allowing better adaptation to the underlying density structure. To understand the performance of the estimator, we analyze its bias and variance. Using a first-order Taylor expansion of $D_i$ around its expectation, we write

$$D_i = \frac{1}{nf(X_{(i)})} + \epsilon_i \qquad (25.47)$$

where $\epsilon_i$ represents the stochastic deviation from the expected value. Substituting this into the estimator,

$$\hat{f}(X_{(i)}) = \frac{1}{n} \sum_{i=1}^{n-1} (f(X_{(i)}) + nf(X_{(i)})^2 \epsilon_i) K \left( \frac{x - X_{(i)}}{D_i} \right) \qquad (25.48)$$

Taking expectations, we obtain the leading-order bias term

$$\mathbb{E}[\hat{f}(x)] = f(x) + \frac{1}{2}h^2 f''(x) + \mathcal{O}(n^{-2/5}) \qquad (25.49)$$
where $h = D_i$ represents the **local bandwidth**. The **variance** of the estimator follows from the variance of the spacings, which satisfies

$$ \mathrm{Var}[D_i] = \mathcal{O}(n^{-2}) \qquad (25.50) $$

Since $\hat{f}(x)$ involves a sum over $n$ terms, its variance is

$$ \mathrm{Var}[\hat{f}(x)] = \mathcal{O}(n^{-1}) \qquad (25.51) $$

Thus, the **mean squared error (MSE)** is given by

$$ \mathbb{E}[(\hat{f}(x) - f(x))^2] = \mathrm{Bias}^2 + \mathrm{Var} = \mathcal{O}(n^{-4/5}) \qquad (25.52) $$

This shows that the **Priestley-Chao estimator achieves the optimal nonparametric rate of convergence**. The kernel function $K(\cdot)$ plays a crucial role in smoothing the estimate. Common choices include:

1. **Uniform kernel:**

$$ K(u) = \frac{1}{2}\mathbf{1}(|u| \le 1) \qquad (25.53) $$

2. **Epanechnikov kernel (optimal in MSE sense):**

$$ K(u) = \frac{3}{4}(1 - u^2)\mathbf{1}(|u| \le 1) \qquad (25.54) $$

3. **Gaussian kernel:**

$$ K(u) = \frac{1}{\sqrt{2\pi}}e^{-u^2/2} \qquad (25.55) $$

The **integrated squared error (ISE)** is used to optimize kernel selection:

$$ \mathrm{ISE} = \int_{-\infty}^{\infty} (\hat{f}(x) - f(x))^2 dx \qquad (25.56) $$

# 25.5 Gasser-Müller kernel estimator

## 25.5.1 Literature Review

Gasser and Muller (1979) [667] wrote one of the foundational papers introducing the Gasser-Müller estimator. It presents the kernel smoothing method and its advantages over existing nonparametric regression techniques, particularly in terms of bias reduction. Gasser and Muller (1984) [668] extended the original estimator to include derivative estimation. It provides rigorous asymptotic analysis of bias and variance, demonstrating the estimator's robustness in various statistical applications. Härdle and Gasser (1985) [669] refined the Gasser-Müller approach by introducing robustness into kernel estimation of derivatives, addressing outlier sensitivity and proposing adaptive bandwidth selection. Müller (1987) [670] generalized the Gasser-Müller estimator by incorporating weighted local regression, improving performance in scenarios with non-uniform data distributions. Chu (1993) [671] proposed an improved version of the Gasser-Müller estimator by modifying its weighting function, leading to better numerical stability and efficiency in practical applications. Peristera and Kostaki (2005) [672] compared various kernel estimators, showing that the Gasser-Müller estimator with a local bandwidth performs better in mortality rate estimation. Müller (1991) [673] addressed the problem of kernel estimators near boundaries, proposing modifications to the Gasser-Müller estimator for improved accuracy at endpoints. Gasser et. al. (2004) [674] expanded on kernel estimation techniques, including the Gasser-Müller estimator, applying them to
shape-invariant modeling and structural analysis. Jennen-Steinmetz and Gasser (1988) [675] developed a unified framework for kernel-based estimators, situating the Gasser-Müller approach within a broader nonparametric regression context. Müller (1997) [676] introduced density-adjusted kernel smoothers, improving upon Gasser-Müller estimators in settings with non-uniformly distributed data points.

### 25.5.2 Analysis of Gasser-Müller kernel estimator

The **Gasser-Müller kernel estimator** is a sophisticated nonparametric method for estimating the probability density function (PDF) of a continuous random variable. It is an improvement upon the classical kernel density estimator (KDE) and is specifically designed to minimize the boundary bias often present in density estimates near the edges of the sample space. This is achieved by placing the kernel functions at the **midpoints** between adjacent data points rather than directly at the data points themselves. The fundamental modification introduced by Gasser and Müller is crucial for improving the estimator's accuracy in regions close to the boundaries, where traditional kernel density estimation methods tend to perform poorly due to limited data near the boundaries.

Let's now describe the Mathematical Framework of the **Gasser-Müller kernel estimator**. Let $X_1, X_2, \dots, X_n$ represent a set of $n$ independent and identically distributed (i.i.d.) random variables drawn from an unknown distribution with a probability density function $f(x)$. The goal of kernel density estimation is to estimate this unknown density $f(x)$ based on the observed data. For the Gasser-Müller kernel estimator $\hat{f}_h(x)$, the core idea is to place the kernel function at the midpoint between two consecutive data points, $X_i$ and $X_{i+1}$, as follows:

$$ \hat{f}_h(x) = \frac{1}{n} \sum_{i=1}^{n-1} K_h \left( x - \frac{X_i + X_{i+1}}{2} \right) \qquad (25.57) $$

where $\xi_i = \frac{X_i+X_{i+1}}{2}$ is the midpoint between consecutive data points, often referred to as the "midpoint shift", $K_h(x) = \frac{1}{h}K(\frac{x}{h})$ is the **scaled kernel function** with bandwidth $h$, $K(x)$ is the kernel function, typically chosen to be a symmetric probability density, such as the Gaussian kernel:

$$ K(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{x^2}{2}} \qquad (25.58) $$

The key difference between the Gasser-Müller estimator and the traditional kernel estimator is the use of midpoints $\xi_i$ instead of the individual data points. The **kernel function** $K_h(x)$ is applied to the **midpoint shift**, effectively smoothing the data and addressing boundary bias by utilizing information from adjacent points.

The **bias** of the estimator can be derived by expanding $\hat{f}_h(x)$ in a Taylor series around the true density $f(x)$. To compute the expected value of $\hat{f}_h(x)$, we first express the expected kernel evaluation:

$$ \mathbb{E}[\hat{f}_h(x)] = \frac{1}{n} \sum_{i=1}^{n-1} \mathbb{E}[K_h(x - \xi_i)] \qquad (25.59) $$

Since $\xi_i$ is the midpoint of adjacent points $X_i$ and $X_{i+1}$, we perform a Taylor expansion around the true density $f(x)$, resulting in:

$$ \mathbb{E}[\hat{f}_h(x)] = f(x) + \frac{h^2}{2} f''(x) \int_{-\infty}^{\infty} u^2 K(u) du + O(h^4) \qquad (25.60) $$

where $\int_{-\infty}^{\infty} u^2 K(u) du$ is the **second moment** of the kernel function, denoted $\sigma_K^2$. The term $\frac{h^2}{2} f''(x) \sigma_K^2$ represents the **bias** of the estimator, which is quadratic in $h$. Thus, the bias decreases
as $h$ becomes smaller, and for sufficiently smooth densities, this bias is small. The main advantage of the Gasser-Müller method is that it leads to a smaller bias compared to standard kernel density estimators, especially at the boundaries. The variance of $\hat{f}_h(x)$ represents the fluctuation of the estimator across different samples. The variance is given by:

$$ \mathrm{Var}[\hat{f}_h(x)] = \frac{1}{n} \left( \int_{-\infty}^{\infty} K^2(u) \, du \right) f(x) \qquad (25.61) $$

where $\int_{-\infty}^{\infty} K^2(u) \, du$ is the second moment of the kernel function $K(x)$. The variance decreases as the sample size $n$ increases, but it also depends on the bandwidth $h$. For a fixed sample size, the variance is inversely proportional to both $h$ and $n$, i.e.,

$$ \mathrm{Var}[\hat{f}_h(x)] \propto \frac{1}{nh} \qquad (25.62) $$

Thus, larger sample sizes and smaller bandwidths lead to smaller variance, but the optimal bandwidth must balance the trade-off between bias and variance. The mean squared error (MSE) combines both the bias and the variance to evaluate the overall performance of the estimator. The MSE is given by:

$$ \mathrm{MSE}[\hat{f}_h(x)] = \mathrm{Bias}^2 + \mathrm{Var} \qquad (25.63) $$

Substituting the expressions for bias and variance, we obtain:

$$ \mathrm{MSE}[\hat{f}_h(x)] = \left( \frac{h^2}{2} f''(x) \sigma_K^2 \right)^2 + \frac{1}{nh} f(x) \int_{-\infty}^{\infty} K^2(u) \, du \qquad (25.64) $$

To minimize the MSE, we select an optimal bandwidth $h_{opt}$. By differentiating the MSE with respect to $h$ and setting the derivative to zero, we obtain the optimal bandwidth that balances the bias and variance:

$$ h_{opt} \propto n^{-1/5}. \qquad (25.65) $$

Thus, the optimal bandwidth decreases as the sample size increases, and this scaling behavior is a fundamental characteristic of kernel density estimation.

The Gasser-Müller estimator performs exceptionally well when compared to other kernel density estimators, such as the Parzen-Rosenblatt estimator. The Parzen-Rosenblatt method places kernels directly at the data points $X_i$, whereas the Gasser-Müller method places kernels at the midpoints $\xi_i = \frac{X_i+X_{i+1}}{2}$. This simple modification significantly reduces boundary bias and results in smoother and more accurate estimates, especially at the boundaries of the sample. Boundary bias occurs in standard KDE methods because kernels at the boundaries have fewer data points to influence them, which leads to a less accurate estimate of the density. Moreover, the Gasser-Müller estimator excels in derivative estimation. When estimating the first or second derivatives of the density function, the Gasser-Müller method provides more accurate estimates with lower variance compared to traditional methods. The use of midpoints ensures that the kernel function is better centered relative to the data, reducing boundary effects that are particularly problematic when estimating derivatives. Regarding the Asymptotic Properties, The Gasser-Müller kernel estimator exhibits asymptotic efficiency. As the sample size $n$ approaches infinity, the estimator achieves the optimal convergence rate of $O(n^{-1/5})$ for the optimal bandwidth $h_{opt}$. This convergence rate is the same as that for other kernel density estimators, indicating that the Gasser-Müller estimator is asymptotically efficient. In the limit, the Gasser-Müller estimator is asymptotically unbiased and asymptotically efficient, meaning that as the sample size increases, the estimator approaches the true density $f(x)$ without bias and with minimal variance. The estimator becomes more accurate as the sample size grows, and the optimal choice of bandwidth ensures that the bias-variance trade-off is well balanced.
In summary, the **Gasser-Müller kernel estimator** offers several distinct advantages over other nonparametric density estimators. Its primary strength lies in its ability to reduce boundary bias by placing kernels at midpoints between adjacent data points. This leads to smoother and more accurate density estimates, especially near the sample boundaries. The optimal choice of bandwidth, which scales as $n^{-1/5}$, balances the bias and variance of the estimator, minimizing the mean squared error. The Gasser-Müller estimator is particularly useful in applications involving density estimation and derivative estimation, where boundary effects and accuracy are crucial. It is a highly effective tool for nonparametric statistical analysis and provides accurate, unbiased estimates even in challenging settings.

## 25.6 Parzen-Rosenblatt method

### 25.6.1 Literature Review

Devroye (1992) [577] investigated the efficiency of superkernels in improving the performance of kernel density estimation (KDE). The study introduces higher-order kernels that lead to reduced asymptotic variance without increasing computational complexity. Zambom and Dias (2013) [658] provided a comprehensive review of KDE, discussing its theoretical foundations, bandwidth selection methods, and practical applications in econometrics. The authors emphasize how KDE can outperform traditional histogram methods in economic data analysis. Reyes et. al. (2016) [659] extended KDE to grouped data, proposing a modified Parzen-Rosenblatt estimator for censored and truncated observations. It addresses practical limitations in standard kernel methods when dealing with incomplete datasets. Tenreiro (2024) [660] developed a KDE adaptation for circular data (e.g., angles and periodic phenomena). It provides exact and asymptotic solutions for optimal bandwidth selection in circular KDE. Devroye and Penrod (1984) [661] proved the consistency of automatic KDE methods. It establishes theoretical guarantees on the convergence of density estimates when the bandwidth is chosen through data-driven methods. Machkouri (2011) [662] established asymptotic normality results for KDE when applied to dependent data, particularly strongly mixing random fields. The paper is crucial for extending KDE applications in time series and spatial statistics. Slaoui (2018) [663] introduced a bias reduction technique for KDE, providing theoretical results and practical improvements over the standard Parzen-Rosenblatt estimator. The modifications significantly enhance density estimation in small-sample scenarios. Michalski (2016) [664] used KDE in hydrology to estimate groundwater level distributions. It demonstrates how KDE outperforms parametric methods in environmental science applications. Gramacki and Gramacki (2018) [665] covered KDE fundamentals, implementation details, and computational optimizations. It is an excellent resource for both theoretical insights and practical applications. Desobry et. al. (2007) [666] extended KDE to unordered sets, exploring its use in kernel-based signal processing. It bridges the gap between statistical estimation and machine learning applications.

### 25.6.2 Analysis of Parzen-Rosenblatt method

The **Parzen-Rosenblatt Kernel Density Estimation (KDE) Method** is a foundational technique in non-parametric statistics that allows for the estimation of an unknown probability density function $f(x)$ from a given sample without imposing restrictive parametric assumptions. Mathematically, let $X_1, X_2, \dots, X_n$ be a set of independent and identically distributed (i.i.d.) random variables drawn from an unknown density $f(x)$. The KDE, which serves as an estimate of $f(x)$, is rigorously defined as

$$ \hat{f}_h(x) = \frac{1}{nh} \sum_{i=1}^{n} K\left(\frac{x - X_i}{h}\right) \qquad (25.66) $$

where $K(\cdot)$ is the **kernel function**, and $h > 0$ is the **bandwidth parameter**. The kernel function $K(x)$ serves as a local weighting function that smooths the empirical distribution, while
the bandwidth parameter h determines the scale over which the data points contribute to the density estimate. The fundamental goal of KDE is to ensure that $\hat{f}_h(x)$ provides an **asymptotically consistent, unbiased, and efficient** estimator of $f(x)$, all of which require rigorous mathematical conditions to be satisfied. There are some important Properties of the Kernel Function. To ensure the validity of $\hat{f}_h(x)$ as a probability density function estimator, the kernel function $K(x)$ must satisfy the following conditions:

1. **Normalization Condition:**

$$ \int_{-\infty}^{\infty} K(x) \, dx = 1 \qquad (25.67) $$

This ensures that the kernel behaves like a proper probability density function and does not introduce artificial bias into the estimation.

2. **Symmetry Condition:**

$$ K(x) = K(-x), \quad \forall x \in \mathbb{R} \qquad (25.68) $$

Symmetry guarantees that the kernel function does not introduce directional bias in the estimation of $f(x)$.

3. **Non-negativity:**

$$ K(x) \ge 0, \quad \forall x \in \mathbb{R} \qquad (25.69) $$

While not strictly necessary, this property ensures that $\hat{f}_h(x)$ remains a valid probability density estimate in a practical sense.

4. **Finite Second Moment (Variance Condition):**

$$ \mu_2(K) = \int_{-\infty}^{\infty} x^2 K(x) \, dx < \infty \qquad (25.70) $$

This ensures that the kernel function does not assign an excessive amount of probability mass far from the origin, preserving local smoothness properties.

5. **Unbiasedness Condition (Mean Zero Constraint):**

$$ \int_{-\infty}^{\infty} xK(x) \, dx = 0 \qquad (25.71) $$

This ensures that the kernel function does not introduce artificial shifts in the density estimate.

Let's discuss the choice of Kernel Function and Examples. Several kernel functions satisfy the above mathematical constraints and are commonly used in KDE:

* **Gaussian Kernel:**

$$ K(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2} \qquad (25.72) $$

This kernel has the advantage of being infinitely differentiable and providing smooth density estimates.

* **Epanechnikov Kernel:**

$$ K(x) = \frac{3}{4}(1 - x^2)\mathbb{I}_{|x| \le 1} \qquad (25.73) $$

This kernel is optimal in the **mean integrated squared error (MISE) sense**, meaning that it minimizes the variance of $\hat{f}_h(x)$ while preserving local smoothness properties.
• Uniform Kernel:

$$K(x) = \frac{1}{2} \mathbb{1}_{|x| \le 1} \qquad (25.74)$$

This kernel is simple but suffers from discontinuities, making it less desirable for smooth density estimation.

Regarding the Asymptotic Properties of the KDE, The **bias** of the KDE can be rigorously derived using a second-order Taylor expansion of $f(x)$ around a given evaluation point. Specifically, if $f(x)$ is twice continuously differentiable, we obtain

$$\mathbb{E}[\hat{f}_h(x)] - f(x) = \frac{h^2}{2} f''(x) \mu_2(K) + O(h^4) \qquad (25.75)$$

where $\mu_2(K) = \int x^2 K(x) dx$ is the second moment of the kernel. The leading term in this expansion shows that the bias is proportional to $h^2$, implying that a smaller $h$ reduces bias, though at the expense of increasing variance. The **variance** of the KDE is given by

$$\mathrm{Var}[\hat{f}_h(x)] = \frac{1}{nh} f(x) R(K) + O\left(\frac{1}{n}\right) \qquad (25.76)$$

where $R(K) = \int K^2(x) dx$ measures the roughness of the kernel function. The key observation here is that variance scales as $O(1/(nh))$, implying that a larger $h$ reduces variance but increases bias. To **minimize the mean integrated squared error (MISE)**, one must choose an optimal bandwidth $h_{opt}$ that balances bias and variance. The **optimal bandwidth** is given by

$$h_{\text{opt}} = \left( \frac{4\hat{\sigma}^5}{3n} \right)^{1/5} \qquad (25.77)$$

where $\hat{\sigma}$ is the sample standard deviation. This scaling rule, known as **Silverman's rule of thumb**, follows from an asymptotic minimization of

$$\mathbb{E}\left[\int_{-\infty}^{\infty} (\hat{f}_h(x) - f(x))^2 dx\right] \qquad (25.78)$$

which encapsulates both bias and variance effects.

In conclusion, the **Parzen-Rosenblatt method** provides a highly flexible, consistent, and asymptotically optimal approach to density estimation. The choice of **kernel function and bandwidth selection** is critical, as they directly impact the **bias-variance tradeoff**. Future refinements, such as **adaptive bandwidth selection and higher-order kernel corrections**, further enhance its performance.
26 Natural Language Processing (NLP)

Natural Language Processing (NLP) is the computational study of linguistic phenomena through
mathematical models that operationalize the syntax, semantics, and pragmatics of human language.
Formally, given a vocabulary $\mathcal{V}$ of discrete tokens (words or subwords), a language is modeled
as a probability distribution $p(w_1, \dots, w_T)$ over sequences of variable length $T$, where $w_i \in \mathcal{V}$.
Statistical language models approximate this distribution through parametric functions $p_{\theta}(w_t|w_{<t})$,
typically implemented via neural networks with architectures specialized for sequential data. The
autoregressive factorization yields the joint probability:

$$
p_{\theta}(w_1, \dots, w_T) = \prod_{t=1}^{T} p_{\theta}(w_t | w_{<t}),
$$

where $w_{<t} = (w_1, \dots, w_{t-1})$. Modern transformer-based models compute $p_{\theta}(w_t|w_{<t})$ using self-
attention mechanisms. For an input sequence $\mathbf{X} \in \mathbb{R}^{T \times d}$ (embedding dimension $d$), the attention
weights $\mathbf{A} \in \mathbb{R}^{T \times T}$ are derived from scaled dot-products:

$$
\mathbf{A} = \mathrm{softmax} \left( \frac{\mathbf{QK}^T}{\sqrt{d_k}} \right), \quad \mathbf{Q} = \mathbf{XW}_Q, \quad \mathbf{K} = \mathbf{XW}_K, \quad \mathbf{V} = \mathbf{XW}_V,
$$

where $\mathbf{W}_Q, \mathbf{W}_K \in \mathbb{R}^{d \times d_k}$ and $\mathbf{W}_V \in \mathbb{R}^{d \times d_v}$ are learnable matrices. The contextualized representa-
tion $\mathbf{H} = \mathbf{AV}$ captures long-range dependencies, overcoming the locality constraints of recurrent
or convolutional architectures.

Word embeddings map discrete tokens to continuous vectors $\mathbf{e}_i \in \mathbb{R}^d$ through functions $f: \mathcal{V} \to \mathbb{R}^d$.
The skip-gram objective maximizes the log-probability of context words $w_{t+j}$ given a target word $w_t$:

$$
\mathcal{L}(\theta) = \sum_{t=1}^{T} \sum_{-c \le j \le c, j \ne 0} \log p_{\theta}(w_{t+j} | w_t),
$$

where $p_{\theta}(w_O|w_I) = \frac{\exp(e_O^T e_I)}{\sum_{w \in \mathcal{V}} \exp(e_w^T e_I)}$ is the softmax over the vocabulary. To circumvent the computational cost of full softmax, negative sampling approximates the objective by contrasting observed pairs $(w_I, w_O)$ against $K$ noise samples:

$$
\mathcal{L}_{\text{NS}}(\theta) = \log \sigma(\mathbf{e}_O^T \mathbf{e}_I) + \sum_{k=1}^{K} \mathbb{E}_{w_k \sim p_{\text{noise}}} [\log \sigma(-\mathbf{e}_{w_k}^T \mathbf{e}_I)],
$$

where $\sigma$ is the sigmoid function and $p_{\text{noise}}$ is a noise distribution (e.g., unigram raised to the 3/4 power). Contextual embeddings like BERT leverage bidirectional transformers trained on masked language modeling:

$$
\mathcal{L}_{\text{MLM}}(\theta) = \mathbb{E}_{w_t \sim \mathcal{D}} [\log p_{\theta}(w_t | \mathbf{w}_{\setminus t})],
$$

where $\mathbf{w}_{\setminus t}$ denotes the sequence with $w_t$ masked. The pretrained model is fine-tuned on downstream tasks via task-specific heads, such as a linear layer for classification:

$$
p_{\theta}(y|\mathbf{w}) = \mathrm{softmax}(\mathbf{W}\mathbf{h}_{[\mathrm{CLS}]} + \mathbf{b}),
$$
where $h_{[CLS]}$ is the embedding of the special classification token.

Sequence-to-sequence tasks (e.g., machine translation) model the conditional distribution $p_{\theta}(\mathbf{y}|\mathbf{x})$ of target sequence $\mathbf{y}$ given source $\mathbf{x}$. The encoder-decoder architecture employs two transformers, with the decoder using masked self-attention to prevent information leakage. Beam search approximates the MAP estimate:

$$ \hat{\mathbf{y}} = \arg \max_{\mathbf{y}} \sum_{t=1}^{|\mathbf{y}|} \log p_{\theta}(y_t|\mathbf{y}_{<t}, \mathbf{x}), $$

maintaining $k$ hypotheses at each step. For structured tasks like parsing, the CKY algorithm dynamically programs the parse tree space under probabilistic context-free grammars, while conditional random fields (CRFs) model dependencies between output tags $\mathbf{y}$:

$$ p_{\theta}(\mathbf{y}|\mathbf{x}) = \frac{1}{Z(\mathbf{x})} \exp \left( \sum_{t=1}^{T} \mathbf{f}_{\theta}(y_t, y_{t-1}, \mathbf{x}) \right), $$

where $Z(\mathbf{x})$ is the partition function computed via the forward algorithm. The Viterbi algorithm decodes the optimal tag sequence in $O(T|\mathcal{Y}|^2)$ time, where $\mathcal{Y}$ is the tag set. These methods unify linguistic structure with statistical learning, enabling machines to process language with human-like competence.

## 26.1 Literature Review

Jurafsky and Martin 2023 [237] wrote book that is a cornerstone of NLP theory, covering fundamental concepts like syntax, semantics, and discourse analysis, alongside deep learning approaches to NLP. The book integrates linguistic theory with probabilistic and neural methodologies, making it an essential resource for students and researchers alike. It thoroughly explains sequence labeling, parsing, transformers, and BERT models. Manning and Schütze 1999 [238] wrote a foundational text in NLP, particularly for probabilistic models. It covers hidden Markov models (HMMs), n-gram language models, and expectation-maximization (EM), concepts that still underpin modern transformer-based NLP models. It also introduces latent semantic analysis (LSA), a precursor to modern word embeddings. Liu and Zhang (2018) [239] presented a detailed exploration of deep learning-based NLP, including word embeddings, recurrent neural networks (RNNs), LSTMs, GRUs, and transformers. It introduces the mathematical foundations of neural networks, making it a bridge between classical NLP and deep learning. Allen (1994) [240] wrote a seminal book in NLP, focusing on symbolic and rule-based approaches. It provides detailed coverage of semantic parsing, discourse modeling, and knowledge representation. While it predates deep learning, it forms a strong theoretical foundation for logical and linguistic approaches to NLP. wrote Koehn (2009) [243] wrote a definitive work on statistical NLP, particularly machine translation techniques like phrase-based translation, alignment models, and decoder algorithms. It remains relevant even as neural translation models (e.g., Transformer-based systems) dominate. We now mention some of the recent works in Natural Language Processing (NLP). Hempelmann [242] explored how linguistic theories of humor can be incorporated into Large Language Models (LLMs). It discusses the integration of formal humor theories into neural models and whether LLMs can be used to test linguistic hypotheses. Eisenstein (2020) [244] wrote a modern NLP textbook that bridges theory and practice. It covers both probabilistic and deep learning approaches, including dependency parsing, sequence-to-sequence models, and attention mechanisms. Unlike many texts, it also discusses ethics and bias in NLP models. Otter et. al. (2018) [245] provides a comprehensive review of neural architectures in NLP, covering CNNs, RNNs, attention mechanisms, and reinforcement learning for NLP. It discusses both theoretical implications and empirical advancements, making it
an essential reference for deep learning in language tasks. The Oxford Handbook of Computational Linguistics (2022) [246] provides a comprehensive collection of essays covering the entire field of NLP and computational linguistics, including morphology, syntax, semantics, discourse processing, and deep learning applications. It presents theoretical debates and practical applications across different NLP domains. Li et. al. (2025) [241] introduced an advanced multi-head attention mechanism that combines explorative factor analysis with NLP models. It enhances our understanding of how transformers encode syntactic and semantic relationships.

# 26.2 Text Classification

## 26.2.1 Literature Review of Text classification

Liu et. al. (2024) [247] provided a systematic review of text classification techniques, covering traditional machine learning methods (e.g., SVM, Naïve Bayes, Decision Trees) and deep learning approaches (CNNs, RNNs, LSTMs, and transformers). It also discusses feature extraction techniques such as TF-IDF, word embeddings, and BERT-based representations. Çekik (2025) [248] introduced a rough set-based approach for text classification, highlighting how term weighting strategies impact classification accuracy. It explores feature reduction and entropy-based selection methods to enhance text classifiers. Zhu et. al. (2025) [249] presented a novel entropy-based prefix tuning method for hierarchical text classification. It demonstrates how entropy regularization can enhance transformer-based classifiers like BERT and GPT for multi-label and hierarchical categorization. Matrane et. al. (2024) [250] investigated dialectal text classification challenges in Arabic NLP. It proposes preprocessing optimizations for low-resource dialects and demonstrates how transfer learning improves classification accuracy. Moqbel and Jain (2025) [252] applies text classification to detect deception in online product reviews. It integrates cognitive appraisal theory and NLP-based text mining to distinguish fake vs. genuine reviews. Kumar et. al. (2025) [253] focused on medical text classification, demonstrating how NLP techniques can be applied to diagnose diseases using electronic health records (EHRs) and patient symptoms extracted from text data. Yin (2024) [254] provided a deep dive into aspect-based sentiment analysis (ABSA), discussing challenges in fine-grained text classification. It introduces new BERT-based techniques to improve aspect-level sentiment classification accuracy. Raghavan (2024) [255] examines personality classification using text data. It evaluates the performance of NLP-based personality prediction models and compares lexicon-based, deep learning, and transformer-based approaches. Semeraro et. al. (2025) [256] introduced EmoAtlas, a tool that merges psychological lexicons, artificial intelligence, and network science to perform emotion classification in textual data. It compares its accuracy with BERT and ChatGPT. Cai and Liu (2024) [257] provides a practical approach to text classification in discourse analysis. It explores Python-based techniques for analyzing therapy talk and sentiment classification in conversational texts.

## 26.2.2 Analysis of Text classification

Text classification is a fundamental problem in machine learning and natural language processing (NLP), where the goal is to assign predefined categories to a given text based on its content. This process involves several steps, including text preprocessing, feature extraction, model training, and evaluation. In this answer, we will explore these steps with a focus on the underlying mathematical principles and models used in text classification. The first step in text classification is preprocessing the raw text data. This typically involves the following operations:

* **Tokenization**: Breaking the text into words or tokens.

* **Stopword Removal**: Removing common words (such as "and", "the", etc.) that do not carry significant meaning.
* **Stemming and Lemmatization**: Reducing words to their base or root form, e.g., "running" becomes "run".

* **Lowercasing**: Converting all words to lowercase to ensure consistency.

* **Punctuation Removal**: Removing punctuation marks.

These operations result in a cleaned and standardized text, ready for feature extraction. Once the text is preprocessed, the next step is to convert the text into numerical representations that can be fed into machine learning models. The most common methods for feature extraction include:

1. Bag-of-Words (BoW) model

2. Term Frequency-Inverse Document Frequency (TF-IDF)

In the first method (Bag-of-Words (BoW) model), each document is represented as a vector where each dimension corresponds to a unique word in the corpus. The value of each dimension is the frequency of the word in the document. If we have a corpus of N documents and a vocabulary of M words, the document i can be represented as a vector $\mathbf{x}_i \in \mathbb{R}^M$, where:

$$ \mathbf{x}_i = [f(w_1, d_i), f(w_2, d_i), \dots, f(w_M, d_i)] \qquad (26.1) $$

where $f(w_j, d_i)$ is the frequency of the word $w_j$ in the document $d_i$. The BoW model captures only the frequency of terms within the document and disregards their order. While simple and computationally efficient, this model does not capture the syntactic or semantic relationships between words in the document.

A more sophisticated and improved representation can be obtained through Term Frequency-Inverse Document Frequency (TF-IDF), which scales the raw frequency of words by their relative importance in the corpus. TF-IDF is a more advanced technique that aims to weight words based on their importance. It considers both the frequency of a word in a document and the rarity of the word across all documents. The term frequency (TF) of a word $w$ in document $d$ is defined as:

$$ \text{TF}(w, d) = \frac{\text{count}(w, d)}{\text{total number of words in } d} \qquad (26.2) $$

The inverse document frequency (IDF) is given by:

$$ \text{IDF}(w) = \log \left( \frac{N}{\text{DF}(w)} \right) \qquad (26.3) $$

where $N$ is the total number of documents and $\text{DF}(w)$ is the number of documents containing the word $w$. The TF-IDF score is the product of these two:

$$ \text{TF-IDF}(w, d) = \text{TF}(w, d) \cdot \text{IDF}(w) \qquad (26.4) $$

There are several machine learning models that can be used for text classification, ranging from simpler models to more complex ones. A common approach to text classification is to use a linear model such as logistic regression or linear support vector machines (SVM). Given a feature vector $\mathbf{x}_i$ for document $i$, the prediction of the class label $y_i$ can be made as:

$$ \hat{y}_i = \sigma(\mathbf{w}^T \mathbf{x}_i + b) \qquad (26.5) $$

where $\sigma$ is the sigmoid function for binary classification, and $\mathbf{w}$ and $b$ are the weight vector and bias term, respectively. The model parameters $\mathbf{w}$ and $b$ are learned by minimizing a loss function, such as the binary cross-entropy loss. More complex models, such as *Neural Networks (NN)*, involve deeper mathematical formulations. In a typical feedforward neural network, the goal is to learn
a set of parameters that map an input vector $\mathbf{x}_i$ to an output label $y_i$. The network consists of multiple layers of interconnected neurons, each of which applies a non-linear transformation to the input. Given an input vector $\mathbf{x}_i$, the output of the network is computed as:

$$ \mathbf{h}_i^{(l)} = \sigma(\mathbf{W}_c^{(l)}\mathbf{h}_i^{(l-1)} + \mathbf{b}_i^{(l)}) \qquad (26.6) $$

where $\mathbf{h}_i^{(l)}$ is the activation of layer $l$, $\sigma$ is the activation function (e.g., ReLU, sigmoid, or tanh), $\mathbf{W}_c^{(l)}$ is the weight matrix, and $\mathbf{b}_i^{(l)}$ is the bias term for layer $l$. The input to the network is passed through several hidden layers before producing the final classification output. The output layer typically applies a softmax function to obtain a probability distribution over the possible classes:

$$ P(y_c|\mathbf{x}_i) = \frac{\exp(\mathbf{W}_c^{T}\mathbf{h}_i + b_c)}{\sum_{c'} \exp(\mathbf{W}_{c'}^{T}\mathbf{h}_i + b_{c'})} \qquad (26.7) $$

where $\mathbf{W}_c$ and $b_c$ are the weights and bias for class $c$, and $\mathbf{h}_i$ is the output of the last hidden layer. The network is trained by minimizing a cross-entropy loss function:

$$ \mathcal{L}(\mathbf{W}, \mathbf{b}) = -\sum_{c=1}^{C} y_{i,c} \log P(y_c|\mathbf{x}_i) \qquad (26.8) $$

where $y_{i,c}$ is the one-hot encoded label for class $c$, and the goal is to minimize the difference between the predicted probability distribution and the true class distribution. Throughout the entire process, optimization plays a crucial role in fine-tuning model parameters to minimize classification errors. Common optimization techniques include stochastic gradient descent (SGD) and its variants, such as Adam and RMSProp, which update model parameters iteratively based on the gradient of the loss function with respect to the parameters. Given the loss function $\mathcal{L}(\theta)$ parameterized by $\theta$, the gradient of the loss with respect to a parameter $\theta_i$ is computed as:

$$ \frac{\partial \mathcal{L}(\theta)}{\partial \theta_i} \qquad (26.9) $$

The parameter update rule for gradient descent is then:

$$ \theta_i \leftarrow \theta_i - \eta \frac{\partial \mathcal{L}(\theta)}{\partial \theta_i} \qquad (26.10) $$

where $\eta$ is the learning rate. For each iteration, this update rule adjusts the model parameters in the direction of the negative gradient, ultimately converging to a set of parameters that minimizes the classification error.

In summary, text classification is an advanced and multifaceted problem that requires a deep understanding of various mathematical principles, including linear algebra, probability theory, optimization, and functional analysis. The entire process, from text preprocessing to feature extraction, model training, and evaluation, involves the application of rigorous mathematical techniques that enable the effective classification of text into meaningful categories. Each of these steps, whether simple or complex, plays an integral role in transforming raw text data into actionable insights using mathematically sophisticated models and algorithms.

# 26.3 Machine Translation

## 26.3.1 Literature Review of Machine Translation

Wu et. al. (2020) [259] introduced end-to-end neural machine translation (NMT), focusing on sequence-to-sequence models, attention mechanisms, and transformer architectures. It explains
encoder-decoder frameworks, self-attention, and positional encoding, laying the groundwork for modern NMT. Hettiarachchi et. al. (2024) [260] presented Amharic-to-English machine translation using transformers. It introduces character embeddings and regularization techniques for handling low-resource languages, a critical challenge in multilingual NLP. Das and Sahoo (2024) [261] discussed word alignment models, a fundamental concept in SMT. It explains IBM Model 1-5, HMM alignments, and the role of alignment in phrase-based models. It also explores challenges in handling syntactic divergence across languages. Oluwatoki et. al. (2024) [262] presented one of the first transformer-based Yoruba-to-English MT systems. It highlights how multilingual NLP models struggle with resource-scarce languages and proposes Rouge-based evaluation for MT systems. Uçkan and Kurt [263] discusses the role of word embeddings (Word2Vec, GloVe, FastText) in MT. It covers semantic representation in vector spaces, crucial for context-aware translation in NMT. discussed multiword expressions (MWEs) in MT, a major challenge in NLP. It covers idiomatic expressions, collocations, and phrasal verbs, showing how neural models struggle with multiword disambiguation. Pastor et. al. (2024) [264] discussed multiword expressions (MWEs) in MT, a major challenge in NLP. It covers idiomatic expressions, collocations, and phrasal verbs, showing how neural models struggle with multiword disambiguation. Fernandes (2024) [265] compared open-source large language models (LLMs) and NMT systems in translating spatial semantics in EN-PT-BR (English-Portuguese-Brazilian Portuguese) subtitles. It highlights the limitations of both traditional and neural MT in capturing contextual spatial meanings. Jozić (2024) [266] evaluated ChatGPT's translation capabilities against specialized MT systems like eTranslation (EU Commission MT model). It shows how general-purpose LLMs can rival dedicated NMT systems but struggle with domain-specific translations. Yang (2025) [267] introduced error-detection models for NMT output, using transformer-based classifiers to detect syntactic and semantic errors in machine-generated translations.

### 26.3.2 Analysis of Machine Translation

Machine Translation (MT) in Natural Language Processing (NLP) is a highly intricate computational task that requires converting text from one language (source language) to another (target language) by using statistical, rule-based, and deep learning models, often underpinned by probabilistic and neural network-based frameworks. The goal is to determine the most probable target sequence $T = \{t_1, t_2, \dots, t_N\}$ from the given source sequence $S = \{s_1, s_2, \dots, s_T\}$, by modeling the conditional probability $P(T | S)$. The optimal translation is typically defined by:

$$T^* = \arg\max_{T} P(T | S) \qquad (26.11)$$

This involves estimating the probability of $T$ given $S$, with the assumption that the translation can be described probabilistically. In the most fundamental form of statistical machine translation (SMT), this probability is often modeled through a series of translation models that decompose the translation process into manageable components. The conditional probability $P(T | S)$ in SMT can be factorized using Bayes' theorem:

$$P(T | S) = \frac{P(S, T)}{P(S)} = \frac{P(T | S)P(S)}{P(S)} \qquad (26.12)$$

Given this decomposition, the core of early SMT models, such as IBM models, sought to model the joint probability $P(S, T)$ over source and target language pairs. Specifically, in word-based models like IBM Model 1, the task reduces to estimating the probability of translating each word in the source language $S$ to its corresponding word in the target language $T$. The joint probability can be written as:

$$P(S, T) = \prod_{i=1}^{T} \prod_{j=1}^{N} t(s_i | t_j) \qquad (26.13)$$
where $t(s_i | t_j)$ is the probability of translating word $s_i$ in the source sentence to word $t_j$ in the target sentence. The estimation of these probabilities, $t(s_i | t_j)$, is typically achieved by analyzing parallel corpora through various techniques such as Expectation-Maximization (EM), which allows the unsupervised learning of these translation probabilities from large amounts of bilingual text data. The EM algorithm iterates between computing the expected alignments of words in the source and target languages and refining the model parameters accordingly. The word-based translation models, however, do not take into account the structure of the language, which often leads to suboptimal translations, especially in languages with significantly different syntactic structures. The challenges stem from the word order differences and idiomatic expressions that cannot be captured through a simple word-to-word mapping. To overcome these limitations, IBM Model 2 introduced the concept of word alignments, where an additional hidden variable $A$ is introduced, representing a possible alignment between words in the source and target sentences. This can be expressed as:

$$P(S, T, A) = \prod_{i=1}^{T} \prod_{j=1}^{N} t(s_i | t_j) a(s_i | t_j) \qquad (26.14)$$

where $a(s_i | t_j)$ denotes the alignment probability between word $s_i$ in the source language and word $t_j$ in the target language. By optimizing these alignment probabilities, SMT systems improve the quality of translations by better modeling the relationship between the source and target sentences. Estimating $a(s_i | t_j)$, however, requires computationally expensive algorithms, which can be handled by methods like EM for iterative refinement.

A more sophisticated approach was introduced with sequence-to-sequence (Seq2Seq) models, which significantly improved the translation process by leveraging deep learning techniques. The core of Seq2Seq is the encoder-decoder framework, where an encoder processes the entire source sentence and encodes it into a context vector, and a decoder generates the target sequence. In this approach, the translation probability is formulated as:

$$P(T | S) = P(t_1 | S) \prod_{i=2}^{N} P(t_i | t_{<i}, S) \qquad (26.15)$$

where $t_{<i}$ denotes the previously generated target words, capturing the sequential nature of translation. The key advantage of the Seq2Seq model is its ability to model entire sentences at once, providing a richer, more flexible representation of both the source and target sequences compared to word-based models. The encoder, typically implemented using Recurrent Neural Networks (RNNs) or more advanced variants such as Long Short-Term Memory (LSTM) or Gated Recurrent Unit (GRU) networks, encodes the source sequence $S$ into hidden states. The hidden state at time step $t$ is computed recursively, based on the input $x_t$ (the source word representation at time step $t$) and the previous hidden state $h_{t-1}$:

$$h_t = f(h_{t-1}, x_t) \qquad (26.16)$$

where $f$ represents the update function, which is often parameterized as a non-linear function, such as a sigmoid or tanh. This recursion generates a sequence of hidden states $\{h_1, h_2, \dots, h_T\}$, each encoding the relevant information of the source sentence. In this model, the decoder generates the target sequence one token at a time by conditioning on the previous tokens $t_{<i}$ and the context vector $c$, which is typically the last hidden state from the encoder. The conditional probability of generating the next target word is given by:

$$P(t_i | t_{<i}, S) = \mathrm{softmax}(W h_t) \qquad (26.17)$$

where $W$ is a learned weight matrix, and $h_t$ is the hidden state of the decoder at time step $t$. The softmax function converts the output of the network into a probability distribution over the vocabulary, and the word with the highest probability is chosen as the next target word.
A significant improvement to Seq2Seq was introduced through the attention mechanism. This allows the decoder to dynamically focus on different parts of the source sentence during translation, instead of relying on a single fixed-length context vector. The attention mechanism computes a set of attention weights $\alpha_{t,i}$ for each source word, which are used to compute a weighted sum of the encoder's hidden states to form a dynamic context vector $c_t$. The attention weight $\alpha_{t,i}$ for time step $t$ in the decoder and source word $i$ is calculated as:

$$ \alpha_{t,i} = \frac{\exp(e_{t,i})}{\sum_{k=1}^{T} \exp(e_{t,k})} \qquad (26.18) $$

where $e_{t,i} = \text{score}(h_t, h_i)$ is a learned scoring function, which can be modeled as:

$$ e_{t,i} = \vec{v}^{\top} \tanh(W_1 h_t + W_2 h_i) \qquad (26.19) $$

This attention mechanism allows the model to adaptively focus on relevant parts of the source sentence while generating each word in the target sentence, thus overcoming the limitations of fixed-length context vectors in long sentences. Training a machine translation model typically involves optimizing a loss function that quantifies the difference between the predicted target sequence and the true target sequence. The most common loss function is the negative log-likelihood:

$$ L(\theta) = - \sum_{i=1}^{N} \log P(t_i | t_{<i}, S; \theta) \qquad (26.20) $$

where $\theta$ represents the parameters of the model. The parameters of the neural network are updated using gradient-based optimization techniques, such as stochastic gradient descent (SGD) or Adam, with the gradient of the loss function with respect to each parameter being computed via backpropagation. In backpropagation, the gradient is computed by recursively applying the chain rule through the layers of the network. For a parameter $\theta$, the gradient is given by:

$$ \frac{\partial L(\theta)}{\partial \theta} = \frac{\partial L(\theta)}{\partial y} \frac{\partial y}{\partial \theta} \qquad (26.21) $$

where $y$ represents the output of the network, and $\frac{\partial L(\theta)}{\partial y}$ is the gradient of the loss with respect to the output. These gradients are then propagated backward through the network to update the parameters, thereby minimizing the loss function. The quality of a translation is often evaluated using automatic metrics such as BLEU (Bilingual Evaluation Understudy), which measures the n-gram overlap between the machine-generated translation and human references. The BLEU score for an n-gram of length $n$ is computed as:

$$ \text{BLEU}(T, R) = \exp \left( \sum_{n=1}^{N} w_n \log p_n(T, R) \right) \qquad (26.22) $$

where $p_n(T, R)$ is the precision of n-grams between the target translation $T$ and reference $R$, and $w_n$ is the weight assigned to each n-gram length. Despite advancements, machine translation still faces challenges, such as handling rare or out-of-vocabulary words, idiomatic expressions, and the alignment of complex syntactic structures across languages. Approaches such as transfer learning, unsupervised learning, and domain adaptation are being explored to address these issues and improve the robustness and accuracy of MT systems.

## 26.4 Chatbots and Conversational AI