We observe a degraded, noisy measurement \( {\bf {y}} \in \mathbb{R}^m\) of the clean image \( {\bf {x}} \in \mathbb{R}^n\) ,

where \( \mathbf{A}\) is called a forward operator and \( \mathbf{n}\) denotes observation noise.

According to Bayes’ rule, the posterior distribution is:

In the Bayesian formulation of inverse problems, there are two different major solution paradigms: Maximum A Posteriori (MAP) estimation and generative posterior sampling. Both approaches operate on the exact same posterior distribution \( p(\mathbf{x}|\mathbf{y})\) , but MAP estimation is an optimization problem, whereas posterior sampling is a simulation problem.

Point Estimation

Posterior sampling

Given a semantic text prompt or conditioning context \( c\) , the goal of conditional generation is to sample from the conditional distribution \( p_t(\mathbf{x}_t|c)\) along the generative trajectory. Applying Bayes’ rule yields:

where \( t \in [0, T]\) denotes the noise level index (diffusion time step) along a generative trajectory.

In practice, computing the classifier score \( \nabla_{\mathbf{x}_t} \log p_t(c|\mathbf{x}_t)\) directly via an external, noise-aware image classifier is notoriously unstable and prone to adversarial artifacts. To bypass this, modern foundation models employ Classifier-Free Guidance (CFG). By jointly training a single network to estimate both the conditional score and the unconditional score (via random condition dropping), the generation process can extrapolate between them during sampling using a guidance scale \( w\) :

This algebraic formulation allows the model to amplify the conditioning signal and trade off sample diversity for visual fidelity without ever explicitly evaluating a classifier.

Even though both reduce to evaluating a conditional score function, they exhibit different behaviors: In inverse problems, the likelihood term \( p(\mathbf{y}|\mathbf{x})\) is defined by physics and noise statistics. If we know the forward operator \( \mathbf{A}\) and assume Gaussian noise, \( p(\mathbf{y}|\mathbf{x}) = \mathcal{N}(\mathbf{y}; \mathbf{Ax}, \sigma^2 \mathbf{I})\) , the gradient

can be computed exactly by linear algebra. In T2I, the term \( p(c|\mathbf{x})\) describes the probability that a human would write the text description \( c\) given the image \( \mathbf{x}\) . This is a semantic mapping that cannot be expressed as a closed-form equation. It must be learned from data pairs using alignment networks (e.g., CLIP [2]).

Understanding that both problems are conditional sampling tasks allows us to combine them into a unified Text-Guided Image Restoration framework: conditional sampling from \( p_t(\mathbf{x}_t|\mathbf{y}, c)\) along a generative trajectory.

When leveraging a Foundation Image Model to solve an inverse problem, we assume the physical measurement \( \mathbf{y}\) and the semantic context \( c\) are conditionally independent given the true image state. This allows us to decompose the total time-dependent score function into a linear combination of vector fields:

Because the physical likelihood \( p(\mathbf{y}|\mathbf{x}_0)\) is defined in the clean image domain, the gradient of the data-fidelity term (the likelihood) cannot be evaluated directly on the noisy latent \( \mathbf{x}_t\) . Instead, implementations, such as Diffusion Posterior Sampling, utilize the foundation model to estimate the clean image \( \hat{\mathbf{x}}_0(\mathbf{x}_t)\) at each step, passing the physical residual back through the network via the chain rule to enforce consistency.

In summary, imaging inverse problems and text-to-image generation leverage complementary conditioning terms grounded on the exact same probabilistic framework to navigate the natural image manifold.

Early image modeling, in the pre-1980 era, assumed image intensity varies smoothly over the spatial coordinates. These global smoothness constraints came in the form of deterministic or stochastic models.

Tikhonov regularization [3], also referred to as ridge regression in statistics, is the classic deterministic approach to stabilizing ill-posed inverse problems by a variational optimization formulation, where the regularizer \( R(\mathbf{x})\) is defined as the squared (\( L_2\) ) norm of the image gradient,

Tikhonov regularization is mathematically equivalent to modeling images as a realization of a homogeneous (stationary) constant-mean Gaussian random process, with the prior

which makes the score linear: \( \nabla_\mathbf{x} \log p(\mathbf{x}) = -2\lambda \mathbf{\nabla^2 x}\) , where \( \nabla^2=\nabla \cdot \nabla\) denotes the Laplacian.

While mathematically elegant and computationally efficient due to linearity of its score—leading to closed-form solutions or simple iterative gradients—this global smoothness prior results in "oversmoothing" of sharp edges because the \( L_2\) penalty treats all variations with equal weight regardless of their structural significance. The stationary Gaussian prior interpretation, which allows complete description of an image by its power spectral density, treating the image as a uniform texture, provided the theoretical foundation for early statistical image processing, such as Wiener filtering in the frequency domain. However, it fails to account for the non-stationary nature of real-world images. Consequently, stationary models cannot adapt to abrupt changes in local image statistics and produce ringing artifacts near edges.

The non-stationary-mean non-stationary-variance (NMNV) Gaussian model [4, 5] is a stochastic framework that characterizes the local statistics of natural images. It is considered a cornerstone in the transition from global (stationary) filters based on global smoothness priors to the spatially adaptive "edge-directed" filters in image processing. Unlike earlier "stationary" models, the NMNV model allows the mean and variance to vary spatially to better capture local structures. The image \( x(i, j)\) is modeled as:

where \( \mu_x(i, j)\) is the local average intensity (non-stationary mean), capturing low-frequency structure, \( \sigma_x(i, j)\) (non-stationary variance) capturing local activity, e.g., edges, textures, and high-frequency details, and \( w(i, j)\) is a zero-mean, unit-variance white Gaussian noise. It states that if we subtract the local mean \( \mu_x(i, j)\) from the image \( x(i, j)\) , the difference can be modeled as spatially-varying scaled white noise.

This model has been the basis for the Kuan filter [5], which is a minimum mean square error (MMSE) estimator for image restoration based on local image statistics. By allowing \( \mu\) and \( \sigma\) to change at every pixel, the model applies stronger smoothing to suppress noise in low-variance regions, while preserving the signal in high-variance regions (edges and texture) to avoid blurring critical details.

Markov Random Field (MRF) is a probabilistic framework to model the spatial context of an image by assuming the state of a pixel depends only on its spatial neighborhood. The application of MRFs to image modeling was revolutionized by the Hammersley-Clifford theorem that establishes equivalence between MRFs and Gibbs distributions. It allowed expressing the probability of an image \( \mathbf{x}\) as a Gibbs measure:

where \( Z\) is the partition function, \( \mathcal{C}\) denotes the set of local cliques (e.g., adjacent pixel pairs), and \( V_c(\cdot)\)  are the handcrafted potential functions that penalize local variations. The seminal work of Geman and Geman [6] demonstrated the power of MRFs in image restoration and segmentation tasks. However, MRFs rely on handcrafted local clique potentials to model piece-wise smoothness while preserving boundaries and struggle to capture the long-range dependencies in high-resolution natural images.

The Total Variation (TV) prior assumes images consist of piecewise-constant areas of intensity separated by sharp edges. Unlike Tikhonov regularization, which penalizes the \( L_2\)  norm of image derivatives, the isotropic TV prior [7] penalizes the spatial \( L_1\) norm of Euclidean gradient magnitudes at each pixel (also referred to as un-squared mixed \( \ell_{1,2}\) norm):

or in discrete notation for a vector \( \mathbf{x}\) as

This is equivalent to assuming an isotropic bivariate Laplace distribution on image gradients:

Because the Euclidean gradient magnitude \( \sqrt{g_1^2 + g_2^2}\) has a non-differentiable conical point at the origin, the Stein score function does not formally exist at zero. To mitigate this, the original ROF paper [7] replaced the gradient magnitude with a smoothed approximation, now known as the Charbonnier penalty:

where \( 0 < \varepsilon \ll 1\) . Alternatively, it is possible to use the piecewise-smooth Huber loss. Note that proximal algorithms (like Chambolle’s Projection Algorithm or the Split Bregman/ADMM methods) to handle non-differentiable norms had not been invented yet when the ROF paper was published in 1992.

The TV prior penalizes the "total amount" of change in pixel intensities. Because of its spatial \( L_1\) structure, it finds it "cheaper" to have one large jump (an edge) than many small random variations (noise). Hence, it flattens noise into smooth regions while leaving important edges sharp and intact. However, modeling images as piecewise-constant blocks with sharp transitions generates artifacts known as the "staircase effect," where regions with soft linear intensity variations appear as step-like stairs.

To overcome the limitations of the ROF model, several generalized TV frameworks were introduced. These include \( L_1\) -TV for impulsive noise, vectorial TV for multi-channel color imaging, non-local TV for texture preservation, and Total Generalized Variation (TGV) for mitigating staircasing artifacts. For a comprehensive review of these architectural variations and their numerical solvers, we refer the reader to Chambolle et al. [8] and the space-variant TV under a Bayesian framework by Pragliola et al. [9].

While natural images are dense in pixel space, they are "sparse" in some linear transform space. The sparse representation paradigm models an image (or patches) as a linear combination of a small number of atoms, chosen from an overcomplete dictionary matrix. Given a fixed dictionary \( \mathbf{D} \in \mathbb{R}^{m \times n}\) (where \( m \ll n\) ), sparse coding maps an observed signal \( \mathbf{y} \in \mathbb{R}^m\) into a latent coefficient vector \( \boldsymbol{\alpha} \in \mathbb{R}^n\) via \( \mathbf{y} = \mathbf{D}\boldsymbol{\alpha}\) . The "sparse" part of the definition requires that most entries in \( \mathbf{\alpha}\) are zero. There exist different mathematical formulations:

1. The Penalized Form [10]:

where \( p\) defines the nature of sparsity constraint, and \( \lambda\) balances reconstruction fidelity vs. sparsity.

2. The Sparsity-Constrained Form:

where \( T_0\) is a strict, predefined budget of allowed non-zero coefficients. This form is natively solved via greedy pursuit metrics like Orthogonal Matching Pursuit (OMP) [11].

3. The Error-Constrained Form (Basis Pursuit Denoising) [10]:

where \( \epsilon\) represents a bounded tolerance for noise or representation error.

If \( p=1\) , the objective functions and constraints are convex and these three forms are weakly equivalent, meaning they map out the same Pareto frontier. If \( p=0\) , the \( L_0\) "norm" counts non-zero coefficients, making the optimization problem non-convex and NP-hard, and the equivalence breaks down.

These optimization formulations can be generalized through a hierarchical Bayesian framework. Let’s treat the latent coefficients \( \mathbf{\alpha}\) as a random variable and the modeling error as additive noise \( n\) : \( \mathbf{y} = \mathbf{D\alpha} + \mathbf{n}\) . Under the MAP estimation objective, finding the optimal representation corresponds to minimizing:

Assuming an independent and identically distributed (i.i.d.) zero-mean Gaussian noise \( \mathbf{n} \sim \mathcal{N}(0, \sigma^2 \mathbf{I})\) , the conditional likelihood function takes the exponential form \( p(\mathbf{y}|\boldsymbol{\alpha}) \propto \exp\left(-\frac{\|\mathbf{y} - \mathbf{D}\boldsymbol{\alpha}\|_2^2}{2\sigma^2}\right)\) . Its negative logarithm generates the standard \( L_2\) data fidelity term:

Assuming that \( \alpha_i\) are i.i.d. and follow a heavy-tailed Laplace distribution \( p(\alpha_i) = \frac{1}{2b}\exp\left(-\frac{|\alpha_i|}{b}\right)\) , the log-prior yields the convex \( L_1\) -norm regularization seen in the LASSO framework:

Dictionary Learning Strategies

Low-Rank Models

Traditional image priors, such as NMNV, TV, and sparse modeling, are local: they assume a pixel is most like its immediate neighbors. Non-local (NL) priors are built on the observation that natural images often contain patches with repeating patterns. If there is a specific texture on a building’s brick or a pattern in a fabric, chances are that exact same pattern exists elsewhere in the image.

Non-Local Means (NL-Means) proposed by Buades et al. [13] was the first practical implementation that demonstrated "non-local averaging" across the image was far superior to "local averaging." A non-local estimate \( \hat{x}(i)\) at pixel \( i\) is a weighted average of all pixels \( j\) in the image:

The weights \( w_{\text{NLM}}(i, j)\) are determined by the structural similarity between the patches \( \mathcal{P}_i\) and \( \mathcal{P}_j\) centered at pixels \( i\) and \( j\) , usually computed using a Gaussian kernel of the Euclidean distance:

where \( \|\cdot\|_{2, a}^2\) denotes a Gaussian-weighted Euclidean distance across the patch footprint, \( h\) acts as a filtering parameter that controls the "strictness" of the similarity, and \( Z(i) = \sum_{j \in \Omega} w_{\text{NLM}}(i, j)\) is a normalizing partition function.

The conceptual precursor to this paradigm is the Bilateral Filter, which smooths images by combining spatial domain closeness with radiometric (intensity) similarity:

where \( i\) and \( j\) represent pixel coordinates, \( y\) represents pixel intensities, \( \sigma_s\) controls spatial decay, and \( \sigma_r\) controls radiometric decay.

NL priors are renowned for restoring highly complex, periodic textures. They are equally potent when regularizing edges. A local estimator can only smooth along an immediate neighborhood of an edge. In contrast, an NL prior identifies structurally similar edge-segments across non-adjacent regions, effectively averaging noise along the entire directional trajectory of the edge without inducing spatial blur.

While a global search over the entire image domain \( \Omega\) is theoretically optimal for infinite ensembles, practical implementations restrict the search space to a localized bounding window \( \mathcal{S}_i\) of size \( S \times S\) centered around the target pixel \( i\) . This bounded search strategy delivers two critical advantages:
1) Mitigation of Structural Hallucinations: Searching the entire image introduces a high probability of encountering "spurious matches"—unrelated distant patches that exhibit low Euclidean distance due to noise realizations. Restricting the search space to \( \mathcal{S}_i\) ensures that matches remain anchored to contextually relevant structures, yielding superior peak signal-to-noise ratio (PSNR) and fewer visual artifacts.
2) Computational Complexity Reduction: An exhaustive global search over an image with \( N\) pixels requires an intractable computational complexity of \( \mathcal{O}(N^2 \cdot |\mathcal{P}|)\) , where \( |\mathcal{P}|\) is the patch cardinality. By bounding the search to a localized window of size \( |\mathcal{S}|\) , the complexity scales linearly with the image size to \( \mathcal{O}(N \cdot |\mathcal{S}| \cdot |\mathcal{P}|)\) , making the framework viable for high-resolution processing pipelines.

Bayesian Formulation of Non-Local Regularization: We can formalize the non-local prior within a MAP-based Bayesian framework. Let \( x\) represent the latent clean image. We treat the non-local weights \( w(i,j)\) as structural hyperparameters that define a non-local graph manifold over the image. The non-local energy functional can be cast as a Graph Dirichlet Form, which functions as the log-prior distribution over the image space:

By applying the negative logarithm, this prior seamlessly pairs with a standard Gaussian noise likelihood model (\( \mathbf{y} = \mathbf{x} + \mathbf{n}\) ), yielding the variational optimization problem:

The Non-Local Laplacian Operator: Taking the derivative of the non-local log-prior with respect to \( x_i\) yields the functional’s score, which corresponds to the Non-Local Laplacian operator (\( \Delta_{\text{NL}}\) ):

The corresponding Euler-Lagrange equation for the MAP inference steps satisfies:

This Bayesian formulation reveals that the non-local prior acts as an adaptive, data-driven diffusion process. In homogeneous regions or along consistent geometric edges, the weights \( w(i,j)\) remain large, driving strong non-local smoothing. Conversely, at isolated structures where no matching patches exist within \( \mathcal{S}_i\) , the weights collapse to zero, completely deactivating the regularizer to preserve the rare detail.

A Taxonomy of Non-Local Frameworks: Algorithms based on NL image priors were the top performers before the deep learning era. The field of non-local image processing has expanded beyond basic patch-averaging filters into sophisticated sparse and low-rank structural group models. Examples include:

1) Block-Matching and 3D Filtering (BM3D) [14], which combines patch self-similarity with transform-domain sparsity, was the "gold standard" of image denoising for over a decade. It consists of three stages: (a) Grouping: Find similar patches (Block-Matching) and stack them into a 3D array (a "3D block").
(b) Collaborative Filtering: Apply a 3D transform (like DCT or Wavelets) to the stack. Because the patches are similar, the signal is highly sparse in 3D, while the noise is spread out. (c) Shrinkage and Inversion: Apply 3D hard-thresholding and empirical Wiener filtering and transform back to the 2D image.

2) Non-Local Total Variation (NL-TV) replaces the local spatial gradient vector with a non-local gradient operator defined as \( \nabla_{\text{NL}} x(i, j) = (x_j - x_i)\sqrt{w(i,j)}\) . Minimizing the \( \ell_1\) -norm of this non-local gradient matrix allows the optimizer to enforce piecewise flatness across non-contiguous but structurally similar regions, eliminating the localized staircasing artifacts typical of standard ROF models.

3) Weighted Nuclear Norm Minimization (WNNM) [15] treats the stack of similar patches as a low-rank matrix and minimizes its nuclear norm (sum of singular values) to enforce structural consistency. WNNM assigns inversely proportional weights to different singular values. Larger singular values (which represent dominant structural components) are penalized less, while smaller singular values (which represent noise and residual artifacts) are penalized more, yielding state-of-the-art preservation of subtle textures.

The word "proximal" comes from "proximity." These algorithms replace "gradient step" (in the direction of the negative gradient) with a generalized projection operator ("proximal operator") jumping to a point nearby an input vector \( \mathbf{v}\) that minimizes:

Parikh and Boyd [16] is a highly-cited reference for Proximal Algorithms. Methods foundational to proximal optimization include Proximal Gradient Descent (PGD), and proximal splitting techniques such as the Alternating Direction Method of Multipliers (ADMM) and Half-Quadratic Splitting (HQS).

The Proximal Gradient Descent algorithm alternates between two simple steps:

Step 1: The Gradient Step (Handling the smooth part) First, we ignore the non-smooth part \( R(\mathbf{x})\) and take a standard gradient descent step using only \( D(\mathbf{x}, \mathbf{y})\) , which moves the solution closer to the data:

Step 2: The Proximal Step (Handling the non-smooth part) Now, we apply the "Proximal Operator" of \( R(\mathbf{x})\) to project the result \( \mathbf{v}_k\) onto the desired structure (e.g., clipping small values to zero for sparsity)

A popular example of a proximal operator is Soft Thresholding. If \( R(\mathbf{x}) = \Vert{}\mathbf{x}\Vert{}_1\) (\( L_1\) norm), which is used for sparse coding or Lasso, the proximal operator \( \text{prox}_{\lambda R}(\mathbf{v})\) applies soft thresholding component-wise:

If the magnitude of an entry \( \vert{}v_i\vert{}\) is smaller than \( \lambda\) (noise), it is clipped to 0; if it is larger (signal), it is shrunk toward 0 by \( \lambda\) .

Iterative Shrinkage-Thresholding (ISTA) is a specific instance of PGD. Because the proximal operator of the \( L_1\) norm is the soft-thresholding (or shrinkage) operator, PGD applied to the Lasso problem, where \( R(x)\) is an \( L_1\) -regularizer, is called ISTA. The connection between "soft-thresholding" and Lasso problems was formalized by Daubechies et al. [17]. Fast ISTA (FISTA) [18] is an accelerated version of ISTA that converges at a rate of \( O(1/k^2)\) using a Nesterov-like momentum term, while ISTA converges at a rate of \( O(1/k)\) , which means the difference between the objective value at iteration \( k\) and the true minimum decreases at a rate proportional to \( 1/k^2\) and \( 1/k\) , respectively. FISTA has been used for sparse coding and image restoration for over a decade. number of iterations.

HQS and ADMM are advanced "Variable Splitting" algorithms. They are used to solve optimization problems when the argument of the "spiky" function involves a linear transformation (like a blur kernel or a matrix), where simple proximal gradient descent cannot be used.

Recall that ISTA solves \( D(\mathbf{x},\mathbf{y}) + R(\mathbf{x})\) . But real-world imaging problems often look like:

where \( \mathbf{L}\) is a linear operator, which maps the image \( \mathbf{x}\) into a domain where it becomes sparse (e.g., derivative operator for TV). Because \( \mathbf{L}\) is inside the function \( R\) , we can’t just threshold \( \mathbf{x}\) to compute the proximal operator, because the sparsity constraint applies to \( \mathbf{Lx}\) , not \( \mathbf{x}\) itself.

Both HQS and ADMM invent a new variable \( \mathbf{z}\) to replace \( \mathbf{Lx}\) . We rewrite the problem as:

where two variables \( \mathbf{x}\) and \( \mathbf{z}\) are linked by a constraint. The solution consists of two steps:
(i) \( \mathbf{x}\) -step: Solves an inversion problem with quadratic loss (usually has closed-form solution).
(ii) \( \mathbf{z}\) -step: Solves the non-smooth problem using a Proximal Operator (easy because \( \mathbf{L}\) is gone).

HQS was introduced to handle non-convex regularizers in image restoration. It relaxes the strict equality constraint \( \mathbf{Lx} = \mathbf{z}\) by adding a penalty if they are different. The additive form of HQS minimizes [19]

The HQS Algorithm:
- Update \( \mathbf{x}\) : Given \( \mathbf{z}_k\) , solve a quadratic inverse problem (hence, it is called "half-quadratic"):

- Update \( \mathbf{z}\) : Apply the Proximal Operator of \( R\) . Minimizing \( R(\mathbf{z}) + \frac{\mu}{2}\Vert{}\mathbf{Lx}_{k+1} - \mathbf{z}\Vert{}^2\) yields:

- Increase \( \mu\) : Make the penalty stricter in the next iteration.

HQS is simple and fast. However, to get an exact solution, \( \mu\) must go to infinity, which can cause numerical instability. ADMM uses an augmented Lagrangian. It keeps the penalty term from HQS, but adds a "Lagrange multiplier" (\( \mathbf{u}\) ) that tracks the error between \( \mathbf{Lx}\) and \( \mathbf{z}\) .
The ADMM Algorithm (3 Steps):
- Update \( \mathbf{x}\) (Inversion): Minimize \( D(\mathbf{x},\mathbf{y})\) while trying to match \( \mathbf{z} - \mathbf{u}\) .

- Update \( \mathbf{z}\) (Denoising/Proximal Step): Apply Proximal Operator to \( \mathbf{Lx} + \mathbf{u}\) .

- Update \( \mathbf{u}\) (Dual Correction): Update the multiplier based on the error.

Unlike HQS, ADMM guarantees convergence to the exact solution without needing the penalty parameter to go to infinity. Boyd et al. [20] is a highly-cited guide to ADMM.

The Plug-and-Play (PnP) framework allows solving complex inverse problems by using an off-the-shelf denoiser as a "prior," even if that denoiser was never explicitly designed for the task at hand.

In the proximal step, we normally solve :

This specific optimization is what a denoiser does. Venkatakrishnan et al. [21] coined the term "Plug-and-Play Priors," demonstrating how in the "z-step" of ADMM, the sophisticated proximal operator of a handcrafted prior (like TV) can be replaced with a simple, state-of-the-art denoiser, \( \mathbf{x}_{k+1} = \text{Denoiser}(\mathbf{v}_k)\) . They used the BM3D (Block-Matching and 3D Filtering) denoiser. Following their work, ADMM and HQS have become practical methods to solve inverse problems using PnP priors.

While [21] demonstrated that the PnP framework worked empirically, they did not provide a proof of convergence for the algorithm. Chan et al. [22] provided the first theoretical grounding for PnP. They established conditions (such as bounded denoiser error) under which the PnP-ADMM algorithm is guaranteed to converge to a fixed point and implicitly contains the prior of the image distribution.

As deep learning began to dominate denoising, researchers realized they could "plug in" a learned denoiser into the PnP framework to get the best of both worlds: the physics-based accuracy of traditional optimization and the speed and power of deep learning. Zhang et al. has [23] demonstrated that a set of CNN denoisers (DnCNN), trained at different noise levels, could be "plugged" into HQS or ADMM outperforming traditional PnP methods while being significantly faster.

While PnP originated as an engineering heuristic, RED [24], introduced by Romano, Elad, and Milanfar, aimed at presenting a mathematical framework to justify using a denoiser inside the optimization loop. The original RED formulation claims to construct an explicit regularization functional \( R(\mathbf{x})\) directly from a user-defined, black-box denoising function \( f(\mathbf{x})\) . The proposed regularizer takes the form:

This functional penalizes the difference between the latent image \( \mathbf{x}\) and its denoised counterpart \( f(\mathbf{x})\) , acting as a structural mismatch penalty. When paired with a data-fidelity term under the MAP inference objective, the optimization problem is stated as:

where \( D(\mathbf{x}, \mathbf{y})\) is the log-likelihood data fidelity term (e.g., \( D(\mathbf{x}, \mathbf{y})=\frac{1}{2}\|\mathbf{y} - \mathbf{Ax}\|_2^2\) for Gaussian noise). To minimize this energy functional via gradient descent or proximal splitting methods, Romano et al. asserted that the gradient of the regularizer simplifies to the standard denoising residual:

This expression allowed RED to bypass the complex derivatives of the underlying denoiser, leading to highly efficient algorithms such as RED-Gradient Descent, RED-Fixed Point, and RED-ADMM.

Despite the empirical success of RED, Reehorst and Schniter [25] subsequently demonstrated that the interpretation of RED as minimizing an explicit regularization function is, in general, mathematically untenable. They showed that the exact gradient of the proposed regularizer \( R(\mathbf{x}) = \frac{1}{2}\mathbf{x}^T \mathbf{x} - \frac{1}{2}\mathbf{x}^T f(\mathbf{x})\) by applying standard vector calculus is actually:

where \( J_f(\mathbf{x}) \in \mathbb{R}^{n \times n}\) represents the Jacobian matrix of the denoiser \( f(\mathbf{x})\) , defined element-wise as \( [J_f(\mathbf{x})]_{i,j} = \frac{\partial f_i(\mathbf{x})}{\partial x_j}\) . For the original RED gradient assumption \( \nabla R(\mathbf{x}) = \mathbf{x} - f(\mathbf{x})\) to hold, the denoiser must satisfy two restrictive mathematical conditions simultaneously: 1) The Jacobian matrix must be perfectly symmetric everywhere, meaning \( J_f(\mathbf{x}) = [J_f(\mathbf{x})]^T\) . 2) The denoiser must satisfy the condition \( J_f(\mathbf{x})\mathbf{x} = f(\mathbf{x})\) , meaning \( f(\mathbf{x})\) behaves as a locally homogeneous function of degree one.

In optimization theory, these conditions restrict \( f(\mathbf{x})\) to behaving identically to a proximal operator. Reehorst and Schniter [25] showed that practical denoisers do not possess a symmetric Jacobian. The patch-matching weights in algorithms such as NL-Means and BM3D depend non-linearly on the input image \( \mathbf{x}\) . Taking the derivative of the patch weights with respect to the input intensities breaks Jacobian symmetry completely. Likewise, convolutional neural networks, such as DnCNN, UNet, utilize non-linear activation functions (e.g., ReLU, GeLU) and data-dependent normalizations. The composition of these layers yields a highly asymmetric Jacobian matrix \( J_f(\mathbf{x}) \neq [J_f(\mathbf{x})]^T\) .

Because practical denoisers violate Jacobian symmetry, the pseudo-gradient step \( \mathbf{x} - f(\mathbf{x})\) in the RED literature is not the gradient of an actual scalar objective function. Consequently, running a RED algorithm does not trace a path down an energy landscape; instead, it executes an iterative path through a non-conservative vector field, meaning proofs of variational optimization convergence do not apply. Yet the RED algorithms remain highly powerful. To resolve this paradox, the literature has abandoned the explicit MAP framework in favor of alternative sound mathematical interpretations of RED.

Reehorst and Schniter [25] proposed replacing explicit regularization with Score-Matching by Denoising (SMD). Drawing from Tweedie’s formula and kernel density estimation, SMD proves that the expression \( \mathbf{x} - f(\mathbf{x})\) represents the score of the data distribution under a MMSE denoiser. Thus, RED algorithms do not minimize an explicit cost function; rather, they perform a root-finding task that matches the score of the underlying data manifold.

The overarching theme of the evolution of image modeling is the transition from early handcrafted image priors defined analytically via regularizers \( \mathcal{R}(\mathbf{x})\) used in point estimation to implicit (data-driven) priors defined by \( \mathcal{R}(\mathbf{x}) = -\log p_\sigma(\mathbf{x})\) , whose gradients are computed operationally by an algorithmic denoiser \( f(\mathbf{x})\) , where \( p_\sigma(\mathbf{x})\) denotes the noise-perturbed distribution, leveraging Tweedie’s formula.

We must emphasize that finding a point MAP estimate via optimization vs. Posterior Sampling via generative/foundation image models are fundamentally different solutions. The "unification" does not mean they will lead to the same solution. Rather, they are unified because they share the same mathematical engine to represent the image prior. Once we have a minimum mean square error (MMSE) denoiser to estimate the Stein score \( \nabla_\mathbf{x} \log p_\sigma(\mathbf{x})\) , the unification lies in the fact that we can use this same neural network to either optimize (MAP) or simulate (sample).

All generative image models (posterior sampling) address learning high-dimensional probability densities, \( p(\mathbf{x})\) , to generate new data. However, they have different architectures and different training objectives. VAEs and GANs perform sampling in one step, while DDPM, score- and flow-based models sample along a trajectory. The VAE consists of a learned encoder and a learned generative decoder, synthesized via a low-variance reparameterized bottleneck to maximize a normalized density lower bound (ELBO). The Diffusion Model is a multi-step, hierarchical latent variable model that uses score matching for optimization. Its architecture is defined by a pre-defined forward process and a single, learned recursive component for the reverse (sampling) process. Generative image models include:

1. DDPM and Score-Based SDEs: While early diffusion models (DDPM by Ho et al.) were described by variational inference and Markov chains, Y. Song et al. unified them via Score-Based Generative Modeling through SDEs. A single network is trained across a continuous spectrum of noise levels \( t \in [0, T]\) to estimate the time-dependent score field \( \nabla_{\mathbf{x}_t} \log p_t(\mathbf{x}_t)\) . By adding controlled noise and stepping along the score vector field using Langevin Dynamics, we can start from pure white noise at \( t=T\) and walk systematically down into a highly realistic, high-probability image at \( t=0\) .

Pan/Zoom

2. Flow Matching and Rectified Flows: While diffusion models navigate the score field using curved, stochastic paths, Flow Matching and Rectified Flows parameterize the problem as a deterministic velocity vector field. Instead of modeling the score directly, they learn a vector field that generates straight, linear trajectories between noise and data, resulting in significantly faster and more stable sampling.

3. Variational Autoencoders (VAE): VAEs avoid modeling the complex image-space score field directly by mapping the data to a low-dimensional, latent space \( \mathbf{z}\) governed by an explicit prior like a standard Gaussian \( p(\mathbf{z})\) . However, both the generative model (the decoder, \( p(\mathbf{x}|\mathbf{z})\) ) and the optimal inference process (the posterior, \( p(\mathbf{z}|\mathbf{x})\) ) are intractable unknowns. VAEs minimize ELBO with the reparameterization trick to find low-variance gradients because it is the only computationally stable way to learn a decoder \( p(\mathbf{x}|\mathbf{z})\) , where the optimal \( \mathbf{z}\) itself is a moving target that must also be learned. However, the gradient of the VAE’s reconstruction loss with respect to the input image acts as an empirical approximation of the score field, pushing corrupted inputs back toward the manifold determined by the latent bottleneck.

4. Latent Diffusion Models (LDM): LDMs are the ultimate synthesis of Variational Inference and Denoising Score Matching. A VAE is trained first to learn a normalized \( p(\mathbf{x}|\mathbf{z})\) with tractable ELBO. Then, a diffusion model is trained within the latent space \( \mathbf{z}\) using score matching to learn the prior \( p(\mathbf{z})\) .

5. Generative Adversarial Networks (GANs): GANs bypass probability density estimation, framing generation as a two-player game. The gradient of the log-odds ratio: \( \nabla_\mathbf{x} \log \left( \frac{D(\mathbf{x})}{1-D(\mathbf{x})} \right) = \nabla_\mathbf{x} \log p_{\text{real}}(\mathbf{x}) - \nabla_\mathbf{x} \log p_{\text{fake}}(\mathbf{x})\) , where \( D(\mathbf{x})\) is the discriminator, pushes generated images toward the high-density "real" image manifold.

The various stages of the evolution of image modeling from handcrafted image priors in pre 1980’s to priors learned from big data is summarized in Figure 2. Clearly, the emergence of today’s diffusion-based foundation models is not a disjointed leap, but rather a continuous evolution centered around learning and traversing the score of the image manifold.

VAEs were almost simultaneously introduced by Kingma et al. [28] and Rezende et al. [29]. Variational inference picks a simple distribution \( q(\mathbf{z})\) (such as a Gaussian), and poses minimization of the difference between our guess \( q(\mathbf{z})\) and the true posterior \( p(\mathbf{z}|\mathbf{x})\) as a variational problem. The measure of difference between distributions is the Kullback-Leibler (KL) Divergence \( D_{KL}(q(\mathbf{z}) || p(\mathbf{z}|\mathbf{x}))\) .

Let’s start with the definition of KL Divergence. It is the expectation (average) of the log difference:

Since \( p(\mathbf{z}|\mathbf{x}) = \frac{p(\mathbf{x}, \mathbf{z})}{p(\mathbf{x})}\) , we have \( \log p(\mathbf{z}|\mathbf{x}) = \log p(\mathbf{x}, \mathbf{z}) - \log p(\mathbf{x})\) . Substituting this back into above:

Notice that \( \log p(\mathbf{x})\) does not depend on \( \mathbf{z}\) . Therefore, it can be pulled out of the expectation: Since we want to solve for \( \log p(\mathbf{x})\) , let’s move everything else to the other side:

Since the KL divergence is always non-negative, dropping the KL divergence term yields an inequality:

i.e., the ELBO is strictly lower than or equal to the Evidence. Since the Evidence \( \log p(\mathbf{x})\) is fixed (the data is given), maximizing the ELBO means minimizing the KL Divergence (the error). From basic probability, \( p(\mathbf{x}, \mathbf{z}) = p(\mathbf{x}|\mathbf{z}) \cdot p(\mathbf{z})\) . Hence, \( \log p(\mathbf{x}, \mathbf{z}) = \log p(\mathbf{x}|\mathbf{z}) + \log p(\mathbf{z})\) . Substituting this into the definition of ELBO, we have

where we regroup the terms to separate the "Likelihood" part from the "Latent" part. In order to express the latent part as KL Divergence, recall: \( D_{KL}(q || p) = E_q [ \log q(\mathbf{z}) - \log p(\mathbf{z}) ]\) and observe that the latent part is the negative of the KL Divergence:

The second term is regularization loss, which forces  the latent distribution to be close to a prior (usually \( \mathcal{N}(0, \textbf{I})\) ), preventing the model from “cheating” by assigning a unique point to every input. Hence, the latent space of VAEs is regularized and continuous allowing learning the underlying probability distribution of a dataset for new sample generation.

In a VAE, the encoder doesn’t output a latent vector; it outputs parameters of a distribution: a mean \( \mu_i\) and a variance \( \sigma_i^2\) for each latent dimension. We then sample a latent vector \( \mathbf{z}\) from this distribution: \( \mathbf{z} \sim \mathcal{N}\left(\mathbf{\mu}, \text{diag}(\mathbf{\sigma}^2)\right)\) where \( \mathbf{\mu} = [\mu_1, \mu_2, \dots, \mu_d]^T\) is the mean vector and \( \mathbf{\sigma}^2 = [\sigma_1^2, \sigma_2^2, \dots, \sigma_d^2]^T\) is the variance vector. Sampling is a stochastic process and has no gradient. Hence, we cannot use the gradients to improve \( \mathbf{\mu}\) and \( \mathbf{\sigma}^2\) estimates. The reparameterization trick expresses \( \mathbf{z}\) as a deterministic function of another random variable \( \mathbf{\epsilon}\) with a standard normal distribution \( \mathbf{\epsilon} \sim \mathcal{N}(0, \textbf{I})\) as [28]:

This way, \( \mathbf{z}\) becomes a differentiable function of \( \mathbf{\mu}\) and \( \mathbf{\sigma}\) , and the gradients can flow from the loss function, through the decoder, through the \( \mathbf{z}\) calculation, and back into the encoder’s \( \mathbf{\mu}\) and \( \mathbf{\sigma}\) layers.

Given a training set of images \( \mathbf{X} = \{\mathbf{x}^{(1)}, \mathbf{x}^{(2)}, \dots, \mathbf{x}^{(N)}\}\) , the workflow to train an Encoder (\( q_\phi\) ) and a Decoder (\( p_\theta\) ) simultaneously for a mini-batch of \( B\) training images \( \{\mathbf{x}^{(b)}\}_{b=1}^B\) is as follows:

5.3.1 Encode

5.3.2 Reparameterization Trick (Sampling)

5.3.3 Decode (Reconstruction)

5.3.4 Compute the Batch ELBO Objective

5.3.5 Backpropagation and Optimization

Once trained, the VAE is used as a generative model to synthesize new samples by following the steps:

1. Discard the Encoder. The inference network \( q_\phi(\mathbf{z}|\mathbf{x})\) is no longer needed.
2. Sample a latent vector \( \mathbf{z}\) from the standard normal distribution: \( \mathbf{z} \sim p(\mathbf{z}) = \mathcal{N}(\mathbf{0}, \mathbf{I})\)
3. Pass \( \mathbf{z}\) through the generative model \( p_\theta(\mathbf{x}|\mathbf{z})\) to yield a generated sample: \( \hat{\mathbf{x}}_{\text{generated}} = \text{Decoder}(\mathbf{z})\)

VAEs provide a principled way to perform non-linear dimensionality reduction while ensuring the latent space is well-behaved for generation. By balancing reconstruction and regularization losses, they allow for the synthesis of new, plausible data samples.

The architecture of a GAN consists of two distinct sub-models. The generator \( G(\mathbf{z}; \mathbf{\theta}_g)\) maps latent noise \( \mathbf{z}\) to a data sample, while the discriminator \( D(\mathbf{x}; \mathbf{\theta}_d)\) maps a sample \( \mathbf{x}\) to a scalar value representing the probability that \( \mathbf{x}\) is real. Both models are implemented as deep neural networks.

The operation of a GAN is as follows: The generator learns to generate an image that is indistinguishable from a real image and the discriminator learns to detect whether a presented image is genuine (from the training set) or a fake (produced by the generator). Initially, the generator is inept, producing outputs that resemble random noise. However, as the training progresses, the discriminator provides feedback to guide the generator toward the real data distribution. This process continues until the generator’s outputs are so realistic that the discriminator can no longer distinguish between them, resulting in an equilibrium where the discriminator’s accuracy falls to 50%, akin to random guessing.

The objective function of a GAN is derived from the principles of game theory and information theory. We first analyze the discriminator’s task as a standalone binary classifier. In a binary classification task, for a given sample \( \mathbf{x}\) and a true label \( y \in \{0, 1\}\) , the binary cross-entropy loss (BCE) is defined as:

where \( \hat{y}\) is the probability predicted by the discriminator \( D(\mathbf{x})\) . For real samples \( \mathbf{x}\) from the true distribution \( p_{data}\) , the label is \( y = 1\) . Substituting this into the BCE formula yields: \( L_{real} = -\log(D(\mathbf{x}))\) . For generated samples \( G(\mathbf{z})\) from the noise distribution \( p_z\) , the label is \( y = 0\) . Substituting this yields: \( L_{fake} = -\log(1 - D(G(\mathbf{z})))\) . To maximize its classification accuracy, the discriminator seeks to minimize the total loss, or equivalently, to maximize the log-probabilities. By taking the expectation over the data and noise distributions, we arrive at the full discriminator objective:

The generator’s goal is the opposite of that of the discriminator. It seeks to minimize the same value function \( V(D, G)\) , effectively maximizing the probability that the discriminator misclassifies its outputs.

This leads to the minimax game:

where the discriminator \( D\) seeks the maximum value by correctly identifying real and fake data, while the generator \( G\) seeks the minimum value by confusing the discriminator. The term “minimax” refers to the search for a saddle point in the parameter space, which corresponds to a Nash Equilibrium in game theory—a state where neither player can improve their outcome by changing their strategy independently.

For any fixed generator \( G\) , there exists a unique optimal discriminator \( D^*_G(\mathbf{x})\) given by

which represents the posterior probability of the sample being “real” given \( \mathbf{x}\) , where \( p_{data}(\mathbf{x})\) is the real data density and \( p_G(\mathbf{x})\) is the probability density induced by the generator. That is, the discriminator is a density ratio estimator that indicates how much more likely a sample is from \( p_{data}(\mathbf{x})\) than from \( p_G(\mathbf{x})\) .

Substituting the expression for the optimal discriminator \( D^*_G\) into the value function \( V(G, D^*_G)\) yields the training criterion \( C(G)\) that the generator attempts to minimize during the minimax game:

which can be shown to be equivalent to

where \( D_{JSD}(P || Q)\) the Jensen-Shannon Divergence (JSD) between the distributions \( P\) and \( Q\) . Hence, the adversarial process is mathematically equivalent to minimizing the JSD between the model’s induced distribution and the true data distribution. The JSD is a symmetric and non-negative measure, which is equal to zero if and only if the two distributions are identical (\( p_{data} = p_G\) ). Therefore, the global minimum of the generator’s training criterion is achieved when the generator has perfectly replicated the data distribution. At this point, \( C(G)=-\log(4)\) , which corresponds to \( D(x) = 0.5\) everywhere.

Several architectural refinements have been developed to improve the performance of GANs. The Deep Convolutional GAN (DCGAN) [31] is perhaps the most significant one for high-resolution image synthesis. It replaced traditional pooling layers with strided convolutions and introduced Batch Normalization into the GAN architecture. Most importantly, it showed that the latent space of a GAN could be used for "vector arithmetic" (e.g., Man with glasses - Man + Woman = Woman with glasses).

The simultaneous training of two networks leads to the following practical difficulties:

The Vanishing Gradient Problem

Mode Collapse and Lack of Diversity

WGAN replaced the traditional Jensen-Shannon Divergence (JSD) cost function with the Wasserstein distance (Earth Mover’s distance), which measures the minimum cost of transporting mass from one distribution to another to solve these issues. Unlike JSD, the Wasserstein distance is continuous and differentiable even when distributions are disjoint, providing a meaningful gradient that guides the generator [32]. That is, instead of a "Discriminator" that tells if an image is real or fake (0 or 1), WGAN uses a "Critic" that provides a continuous score, ensuring the generator always has a gradient to follow.

While WGAN was theoretically sound, its method of enforcing the "1-Lipschitz" constraint (weight clipping) led to optimization difficulties. Gulrajani et al. [33] introduced the Gradient Penalty (GP), which penalizes the model if the gradient norm of the critic deviates from 1. This became the standard implementation for almost all GAN research that followed, including the famous StyleGAN series.

Phase 1: Causal Convolutions (RNN/CNN)
PixelRNN and PixelCNN [34] use Recurrent Neural Networks (RNNs) and causal masking, respectively, to maintain "causality"—ensuring a pixel only "sees" pixels that came before it. PixelCNN was presented as a faster alternative to PixelRNN. While PixelRNN was more accurate, PixelCNN replaced the slow, sequential LSTMs with masked convolutions. This allowed the model to be trained in parallel, significantly reducing training time. Gated PixelCNN [35] replaced simple ReLU activations with gated activation units (similar to LSTMs) and introduced a "stack" of two convolutions (horizontal and vertical) to eliminate the "blind spot." While PixelCNN [34] proved parallel training was possible for AR models, the Gated PixelCNN [35] proved that these models could be "conditioned" (e.g., "Generate a picture of a 5"), which laid the groundwork for modern text-to-image systems. The follow-up work PixelCNN++ [36] is a refinement that improved efficiency and quality by using a discretized logistic mixture likelihood.

Phase 2: Self-Attention (The Transformer)
As Transformers began to dominate NLP, researchers also applied self-attention to images. This allowed models to capture long-range dependencies (e.g., ensuring the left eye matches the right eye) that were difficult for local convolutional windows to manage. Some examples include: Image Transformer [37] is the first major application of the Transformer architecture to image synthesis, replacing convolutions with self-attention. iGPT (Image GPT) [38] by OpenAI demonstrated that a massive Transformer trained solely to predict the next pixel could learn powerful visual representations without any labels.

Phase 3: The Discrete Latent Space
Generating high-resolution images pixel-by-pixel is painfully slow. Vector Quantized Variational Autoencoder (VQ-VAE) [39] turns images into a discrete sequence of latent codes, making AR modeling practical for high-resolution image and video generation. The VQ-VAE architecture is a standard Encoder-Decoder structure. Encoding takes the image and maps it to a grid of discrete indices (tokens). Once the VQ-VAE is trained, it compresses an image into a grid of integers (e.g., \( 32 \times 32\) codebook indices). Decoding takes this grid of indices and reconstructs the image in parallel. This is where the AR model comes in: An AR model (typically a PixelCNN or a Transformer) is trained to predict the latent codes one by one. The original DALL-E [40] utilized a VQ-VAE and a massive AR Transformer to generate images from text, proving that the AR paradigm could scale to multi-modal foundation models.

While diffusion models have largely superseded AR models for pure image synthesis speed, the AR paradigm remains the backbone of multimodal intelligence in Foundation Image Models. It allows for a unified "vocabulary" where text, images, and audio are all just tokens in a single, long sequence.

Sohl-Dickstein et al. [41] introduced the concept of using a Markov chain to slowly destroy and then restore structure in a data distribution, inspired by statistical physics. However, the implementation was dense and lacked the visual fidelity required for applications. Ho et al. [42] took the physics-inspired foundation and re-engineered it into a DDPM denoising pipeline and showed that the model performed significantly better if it predicted the added noise (\( \epsilon\) ) at each step.

To unify DDPM, SGM and FM paradigms, we define a continuous time variable \( t \in [0, T]\) . Without loss of generality, let \( t=0\) index the clean data state, \( \mathbf{x}_0 \sim p_{\text{data}}(\mathbf{x})\) , and let \( t=T\) index the fully perturbed noise state, \( \mathbf{x}_T \sim p_{\text{prior}}(\mathbf{x})\) . The forward process, which systematically degrades the underlying structure of the data manifold, can be formalized as an Itô Stochastic Differential Equation (SDE):

where \( f(t)\) represents the deterministic drift coefficient, \( g(t)\) is the scalar diffusion coefficient controlling the injection rate of stochastic noise, and \( \mathbf{w}\) is a standard Brownian motion (Wiener process).

A foundational result in stochastic calculus established by Anderson [43] dictates that the reverse-time process—moving from the unstructured prior \( p_{\text{prior}}\) back to the data manifold \( p_{\text{data}}\) —is also a diffusion process running backward in time. It satisfies the reverse-time SDE:

where \( d\bar{\mathbf{w}}\) is a backward Wiener process running from \( t=T\) to \( t=0\) , and \( p_t(\mathbf{x})\) represents the marginal probability density of \( \mathbf{x}\) at time \( t\) . The term \( \nabla_\mathbf{x} \log p_t(\mathbf{x})\) is the Stein score function—a time-dependent vector field pointing toward regions of higher probability density at noise scale \( t\) .

The Probability Flow ODE

DDPMs and SGMs were developed independently using discrete-time formulations. The continuous-time SDE framework unifies them as specific parameterizations of the scalar drift \( f(t)\) and diffusion \( g(t)\) coefficients. However, understanding the operational mechanics of this unification requires grounding the framework in the explicit variational lower bound (ELBO) established by Ho et al. [42].

The Discrete ELBO Starting Point in Ho’s Formulation

Transition to Continuous-Time Score Matching

While score-based diffusion models provide a rigorous framework, their trajectories are inherently anchored to stochastic processes. This constraints the paths to be highly non-linear, which increases numerical complexity during ODE integration and necessitates a large number of sampling steps.

Flow Matching (FM) (Albergo and Vanden-Eijnden [48], Lipman et al. [49]) circumvents this limitation by abandoning the underlying SDE framework entirely. Instead, FM models generative trajectories strictly through a deterministic fluid-dynamics perspective using a time-dependent velocity vector field \( \mathbf{v}_t(x)\) . The generative process is governed by a simple, non-stochastic ODE:

Conditional Flow Matching (CFM)

Optimal Transport and Rectified Flows

The profound realization linking these three frameworks is that the Probability Flow ODE of score-based diffusion is identical to a Flow Matching ODE. By equating the deterministic velocity field \( \mathbf{v}_t(\mathbf{x})\) directly to the drift and score-diffusion terms of the PF-ODE, we obtain the universal trajectory identity:

Through this synthesis, the evolution of trajectory-based models can be viewed as a continuous optimization of path geometry: transitioning from stochastic paths (SDEs), to deterministic curves (Scores and PF-ODEs), and ultimately to straight vector lines (Flow Matching), systematically refining sampling efficiency while preserving the expressive capacity of the underlying data manifold.

Image-to-Image Schrödinger Bridge (SB) [51] generalizes the Flow Matching (FM) framework by formulating image generation as an Entropic Optimal Transport (EOT) problem that brings back the benefits of stochasticity. Technically, FM is the "zero-noise limit" of the SB, which recovers the optimal path that minimizes the KL divergence to a reference prior SDE, while satisfying exactly the marginal distributions at both the initial and final times. In image modeling, SBs can seamlessly construct optimal paths between domains, such as corrupted images and clean images, or sketches and photorealistic outputs.

In order to generate a \( 1024 \times 1024\) image, denoising diffusion models have to calculate over a million pixels at every timestep. This is computationally expensive. To handle high-resolution images efficiently, Rombach et al. [52] apply the diffusion process in the latent space of a pre-trained autoencoder.

A Latent Diffusion Model (LDM) consists of three blocks: (i) The VAE: An encoder that turns a high-resolution image into a compressed latent representation (\( z\) ). VAEs are often labeled as f8 or f16, which refers to the downsampling factor: e.g., for \( 512\times512\) input image, f8 refers to \( 512 \div 8\) , i.e., \( 64\times64\) latents, which was used for Stable Diffusion v1 and v2. f4 provides more detail but is much slower. f16/f32 are very fast, but the VAE "forgets" details. (ii) Diffusion Process: The U-Net has to denoise the smaller latent grid (e.g., \( 64 \times 64\) ). (iii) The Decoder: Once the latents are "clean," the VAE decoder expands them back up into a full-resolution, high-fidelity image.

To ensure the latent space is well-behaved for a diffusion model, Rombach et al. [52] explored two types of regularization to train the VAE: (i) Kullback-Leibler: Similar to a standard VAE, it imposes a penalty on the learned latents to keep them close to a standard normal distribution. (ii) Vector Quantization: It uses a codebook layer (like VQ-VAE) but the quantization is often integrated into the decoder. This makes the latent space discrete. Stable Diffusion primarily uses the KL-regularized VAE because it results in smoother, continuous latent spaces that are easier for the Diffusion U-Net to model. Esser et al. [53] developed the VQ-GAN (Vector Quantized GAN), which provided perceptual compression necessary to represent images as compact codes for transformer based denoising while maintaining sharp visual details.

During the training of the Diffusion Model (the U-Net), the VAE is usually frozen. We only use it to encode images into latents once, and then the U-Net is trained to learn how to denoise those latents.

In the evolution of image modeling, LDM [52] serves as a bridge between diffusion models and Foundation Image Models. The significance of LDM lies in: (i) Efficiency: Reducing the search space from \( H \times W \times 3\) to a latent representation (typically 1/8th the size) allowed researchers to train models on billions of images. (ii) Conditioning: It introduced Cross-Attention layers into the diffusion backbone, providing a flexible way to inject text prompts, depth maps, or semantic masks into the modeling process.

While the U-Net has been the workhorse of original LDMs, the transition to FIMs was triggered by the realization that Transformers, which revolutionized natural language processing (NLP), could be tamed for images. By replacing the U-Net with a Diffusion Transformer (DiT) architecture, proposed by Peebles and Xie [54], see Figure 3, researchers could finally apply Scaling Laws to image processing and vision. DGMs leverage convolutional frameworks (like U-Nets) which carry a strong structural inductive bias assuming localized pixel correlations. FIMs employ patch-based DiTs, utilizing global self-attention to capture long-range dependencies. By treating the latent representations of an image as a sequence of tokens ("Patchification")—analogous to words in a sentence—DiTs apply the same highly scalable architectures that power Large Language Models (LLMs) to model the continuous visual world. Furthermore, AR foundation models, which predict the next image token in a quantized latent space offer unified frameworks for both text and image generation. However, a model does not become a "Foundation Model" just because it uses a Transformer backbone. Other factors distinguishing DGMs from FIMs are:

Parameter and Dataset Scale: DGMs are typically smaller networks trained on specialized, narrow datasets (e.g., human face sets). FIMs scale to billions of parameters and undergo pre-training on diverse, web-scale datasets spanning billions of image-text pairs.

Pan/Zoom

Multiple Modalities: The integration of multiple modalities (text, images, video, and sequence semantics) within shared latent and embedding spaces is a fundamental requirement of FIM. Traditional DGMs are largely unimodal; they focus strictly on mapping noise vectors to pixel spaces based on visual statistics. FIMs, however, use joint cross-modal alignment mechanisms (such as CLIP) to structurally intertwine textual semantics with high-dimensional visual features. This multimodal alignment guides the model’s reverse sampling mechanics directly to a highly specific, semantically correct sub-manifold, eliminating the structural "hallucinations" common in older DGM configurations. Task Generalization (Narrow vs. Zero-Shot): DGMs are bound to narrow task scopes and struggle with out-of-distribution inputs. FIMs exhibit robust zero-shot generalization, enabling them to successfully execute tasks for which they were never explicitly or narrowly trained. In terms of versatility as an image prior, DGMs are fragile priors; if deployed to solve an inverse problem outside their training distribution, they can force incorrect features into the reconstruction. FIMs act as zero-shot universal image priors that can be plugged into arbitrary iterative solvers to fix any degradation without retraining the base network.

Foundation image models (FIM) can be broadly classified into the following categories:

Representation Models (Vision-language alignment)

Promptable Generative FIMs (Synthesizing pixels)

Promptable Discriminative Foundation Vision Models (Parsing and segmenting)

We can classify methods to solving inverse problems using FIM priors as training-time, supervised finetuning approaches vs. inference-time zero-shot physics-guided optimization on frozen FIM manifolds.

A. Training-time Task-Specific, Supervised Fine-Tuning of FIMs via Conditional Adapters

B. Inference-Time Zero-Shot, Physics-Guided Optimization on Frozen FIM Manifolds

When an INR is initialized with random weights, it doesn’t produce white noise. Surprisingly, it produces smooth, coherent, low-frequency patterns. That is, the structure of a neural network (convolutions or fully connected layers) is naturally biased toward generating smooth images rather than high-frequency noise. This is closely related to the Deep Image Prior (DIP) concept, where the structure of the network itself regularizes the signal. If we train an INR to replicate a noisy image, the network will learn the "clean" image structure first (because it’s easy to represent) and will struggle to learn the random noise (because it’s hard to represent). This makes INRs powerful unsupervised tools for tasks like Denoising (stop training before it learns the noise).

A major challenge with standard MLPs is that they are notoriously bad at learning high-frequency details (sharp edges and fine textures). This is known as Spectral Bias—networks converge on low-frequency components much faster. To make INRs viable for high-fidelity imaging, Sitzmann et al. [60] introduced Positional Encoding (similar to Transformers) or sinusoidal activation functions (SIREN):

By mapping coordinates into a higher-dimensional Fourier feature space, the network can tune specific frequencies, allowing it to capture both the smooth gradient of sky and the sharp edge of a building.