2.3.1.1 The standard basis and the identity transformation.

2.3.1.2 Example: a \( 90^\circ\) rotation.

2.3.1.3 Beyond rotations: general linear transformations.

2.3.2.1 Matrix–vector multiplication: the column-combination view.

2.3.2.2 Matrix–vector multiplication: the row-by-column view.

2.3.2.3 Extending to matrix–matrix multiplication.

Given a square matrix \( A \in \mathbb{R}^{d \times d}\) , a nonzero vector \( v \in \mathbb{R}^d\) is called an eigenvector of \( A\) if

for some scalar \( \lambda \in \mathbb{R}\) , called the corresponding eigenvalue. This equation tells us that when \( A\) acts on \( v\) , the output remains along the same line spanned by \( v\) ; its magnitude is scaled by \( |\lambda|\) , and if \( \lambda<0\) the orientation is reversed.

For a general vector, multiplication by \( A\) may rotate, stretch, shear, or otherwise distort the direction. Eigenvectors therefore identify directions that are not mixed with others under the transformation: the output remains along the same line, and differs only by a scalar factor \( \lambda\) . Eigenvectors can therefore be viewed as directions of pure scaling under the transformation defined by \( A\) .

If a matrix \( A \in \mathbb{R}^{d \times d}\) has \( d\) linearly independent eigenvectors, we may assemble them as the columns of a matrix \( V\) , and place the corresponding eigenvalues along the diagonal of a matrix \( \Lambda\) .

Concretely, if \( v_1, v_2, \dots, v_d\) are eigenvectors of \( A\) with corresponding eigenvalues \( \lambda_1, \lambda_2, \dots, \lambda_d\) , then

In compact matrix form, this relationship can be written as

where

Rearranging this expression gives the eigen-decomposition

For a general square matrix, an eigen-decomposition may fail to exist over the reals, or may involve complex eigenvalues or an insufficient set of linearly independent eigenvectors. However, when \( A\) is a real symmetric matrix (\( A = A^\top\) ), the situation becomes much simpler.

In this case, several important properties hold:

• All eigenvalues of \( A\) are real.

• There exists an orthogonal matrix \( U\) , whose columns are eigenvectors of \( A\) , such that

  \[ A = U \Lambda U^\top, \]

  where the columns of \( U\) have unit length and are mutually orthogonal, with \( U^\top U = I\) .

This special structure means that symmetric matrices can always be diagonalised using an orthonormal change of basis, so the change of coordinates itself preserves lengths and angles. As we will see shortly, this property is crucial for PCA, since the empirical covariance matrix is always symmetric.

With these tools in place, we are now ready to formulate and solve the PCA problem.

2.5.1.1 Centring the data and defining projected variance.

2.5.1.2 Covariance matrix form.

The covariance matrix \( \Sigma\) is symmetric by construction. As discussed in the previous section, symmetric matrices admit an orthogonal eigen-decomposition:

where the columns of \( U\) are orthonormal eigenvectors and the eigenvalues are ordered as

To solve the optimisation problem, we express the objective in this eigenvector coordinate system.

2.5.2.1 Step 1: Change of basis to the eigenvector coordinates.

2.5.2.2 Step 2: Diagonalising the objective.

2.5.2.3 Step 3: Bounding the variance.

2.5.2.4 Step 4: Optimal direction.

2.5.2.5 Subsequent principal components.

Once the principal components have been identified, PCA provides a concrete procedure for dimensionality reduction by projecting data onto a low-dimensional subspace and, if desired, reconstructing it back in the original space.

Let

denote the matrix formed by the top \( k\) principal components, that is, the eigenvectors of the covariance matrix \( \Sigma\) corresponding to the largest eigenvalues \( \lambda_1 \ge \lambda_2 \ge … \ge \lambda_k\) . PCA reduces the dimensionality of a data point \( x \in \mathbb{R}^d\) by representing it using only these \( k\) directions that capture the largest variance.

2.5.3.1 Projection (encoding).

2.5.3.2 Reconstruction (decoding).

2.5.3.3 Explained variance ratio.

2.6.1.1 Trace of the covariance matrix.

2.6.1.2 Similarity transformations.

2.6.1.3 Trace invariance under change of basis.

We now revisit PCA from the perspective of a change of basis. For an individual data point \( x \in \mathbb{R}^d\) , the low-dimensional representation is obtained by projection onto the principal component directions, \( z = U_k^\top x.\) At the level of the entire dataset, with centred data points collected row-wise in a matrix \( X\) , applying PCA corresponds to the linear transformation \( Z = XU,\) where \( U\) is the orthogonal matrix whose columns are the eigenvectors of the covariance matrix \( \Sigma\) .

We now examine how this change of basis affects the covariance structure of the data.

2.6.2.1 Covariance of the transformed data.

2.6.2.2 Consequences of the similarity transformation.

The discussion above shows that PCA may be understood as a change of basis that diagonalises a covariance matrix and reorganises variance along orthogonal directions. This naturally leads to a broader question:

What if the matrix we want to understand is not a covariance matrix, not symmetric, or not even square?

The Singular Value Decomposition (SVD) provides the corresponding general framework [5].

Eigen-decomposition is based on the relation \( A v = \lambda v,\) which identifies directions that are mapped onto themselves up to scaling. Along such directions, the transformation is especially simple: the output remains on the same line, scaled by \( \lambda\) .

This viewpoint is powerful, but limited. It applies only to square matrices, and additional conditions are generally needed to guarantee a full and well-behaved real decomposition.

SVD generalises the idea. Instead of asking for directions that are preserved by the matrix, it asks a more flexible question: which orthogonal input directions are mapped to which orthogonal output directions, and by how much are they stretched?

2.6.3.1 The SVD viewpoint.

Let \( U_k = [\,u_1,\dots,u_k\,] \in \mathbb{R}^{d\times k}\) denote the matrix whose columns are the top \( k\) eigenvectors of the covariance matrix. For a centred data point \( x_i \in \mathbb{R}^d\) , PCA produces the low-dimensional representation

which records the coordinates of \( x_i\) along the principal directions. Reconstructing from this representation gives

which is precisely the orthogonal projection of \( x_i\) onto the \( k\) -dimensional subspace spanned by the columns of \( U_k\) .

Once both the original point \( x_i\) and its reconstruction \( \hat{x}_i\) are available, the natural measure of information loss is their difference. This leads to the squared reconstruction error

Summing across all data points gives the total reconstruction loss

From this viewpoint, PCA may be formulated as

which asks:

Which orthonormal \( k\) -dimensional subspace minimises the distortion caused by compressing the data and then reconstructing it?

This reconstruction formulation is fully equivalent to the earlier variance-maximisation view of PCA, and is closely related to the classical best low-rank approximation perspective [6]. Although the objective is written differently, both formulations recover the same optimal \( k\) -dimensional subspace spanned by the leading eigenvectors of the covariance matrix.

This perspective is especially valuable because it recasts PCA as a pair of maps:

That is, PCA may be viewed as an encoder followed by a decoder, both linear and analytically determined.

2.7.1.1 Autoencoders as learned encoder–decoder models.

The connection between linear autoencoders and PCA is classical [7]. By itself, it does not guarantee that PCA will be recovered. To recover PCA exactly, two structural ingredients matter.

2.7.2.1 1. Linearity.

2.7.2.2 2. Subspace versus basis.

2.7.2.3 Subspaces and column spaces.

2.7.2.4 Subspace recovery under reconstruction loss.

The true power of autoencoders appears once the encoder and decoder are allowed to be nonlinear. In the linear case, if the decoder has the form

then all reconstructions lie in the single linear subspace

No matter how the latent code \( z\) varies, the output remains confined to one flat \( k\) -dimensional surface inside \( \mathbb{R}^d\) .

If the decoder becomes nonlinear, this geometry changes fundamentally. The set of possible reconstructions becomes

which is no longer restricted to a linear subspace. Instead, it may form a curved low-dimensional surface embedded in the ambient space. Such a surface is called a manifold.

A nonlinear autoencoder therefore generalises PCA in a profound way. Rather than learning a single flat subspace, it can learn a curved manifold that adapts more flexibly to the intrinsic structure of the data.

2.7.3.1 Nonlinear encoder and decoder.

2.7.3.2 Why deterministic autoencoders are not yet generative models.

Before introducing the model formally, it is worth pausing on the term latent variable. The word latent simply means hidden: it refers to a variable that is not directly observed, but is assumed to lie behind the data and explain part of its structure. In PPCA, the observed vector \( \mathbf{x}\) is the data point we see, while \( \mathbf{z}\) is a lower-dimensional hidden representation that we do not observe directly. The model assumes that \( \mathbf{x}\) is generated from \( \mathbf{z}\) , rather than treating \( \mathbf{x}\) as an isolated object. This shift—from observed data alone to observed data generated from hidden causes—is one of the defining ideas of latent-variable generative modelling.

The PPCA generative process consists of two steps:

1. Sample a latent vector

   \[ \mathbf{z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}), \]

   from a standard Gaussian prior.

2. Generate an observation by applying a linear transformation, adding a mean offset, and then injecting Gaussian noise \( \boldsymbol{\epsilon}\) :

   \[ \mathbf{x} = \mathbf{W}\mathbf{z} + \boldsymbol{\mu} + \boldsymbol{\epsilon}, ~~ \boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \alpha^2 \mathbf{I}). \]

Here,

• \( \mathbf{z} \in \mathbb{R}^k\) is the latent variable,
• \( \mathbf{x} \in \mathbb{R}^d\) is the observed data point,
• \( \mathbf{W} \in \mathbb{R}^{d\times k}\) maps latent coordinates into the observation space,
• \( \boldsymbol{\mu} \in \mathbb{R}^d\) is the mean of the data,
• and \( \alpha^2\) controls the isotropic observation noise.

This model already contains the essential ingredients of many later generative models: a prior over latent variables, a rule that maps latent variables to observations, and a probability distribution over possible outputs. The key novelty relative to PCA is that the latent representation no longer determines a single point exactly. Instead, it determines the centre of a distribution of plausible observations.

Although the generative equation

looks similar to a deterministic mapping, its meaning is fundamentally probabilistic because of the random noise term. Conditioned on a latent vector \( \mathbf{z}\) , the observation is Gaussian:

This expression is worth pausing over. Once \( \mathbf{z}\) is fixed, the mean of \( \mathbf{x}\) is fixed at \( \mathbf{W}\mathbf{z} + \boldsymbol{\mu}\) , but the realised observation is not. Instead, it is sampled from a Gaussian cloud centred at that mean. Thus, a single latent point corresponds to many possible observations, all distributed around the same linear reconstruction.

This is the central probabilistic shift. In PCA or a deterministic linear autoencoder, a latent coordinate gives one reconstruction. In PPCA, the same latent coordinate defines a distribution of possible reconstructions.

To make the geometry concrete, consider the simplest nontrivial case: a one-dimensional latent variable \( z \sim \mathcal{N}(0,1)\) and a two-dimensional observation \( \mathbf{x} \in \mathbb{R}^2\) . Then \( \mathbf{W}\) is a \( 2\times 1\) vector and the model becomes

Without noise, the term \( \mathbf{W}z + \boldsymbol{\mu}\) traces out a line in the two-dimensional observation space. Each value of \( z\) would map to a single point on that line, exactly as in a deterministic linear latent representation. Once Gaussian noise is added, however, each point on the line becomes the centre of a small Gaussian cloud of possible observations.

Geometrically, PPCA therefore describes the data as being generated near a low-dimensional latent subspace, but not exactly on it. The latent variable controls where we are along the subspace, while the noise term allows observations to deviate from it. This is already much closer to a realistic generative story: real data rarely lie perfectly on an ideal low-dimensional surface.

At this point it is natural to ask how PPCA relates back to classical PCA. The answer is most clearly seen through the zero-noise limit. Recall the PPCA model:

As the observation noise tends to zero \( \alpha^2 \to 0,\) , the conditional distribution

becomes increasingly concentrated around its mean. In the limit, the randomness disappears and each latent coordinate \( \mathbf{z}\) maps to a single deterministic point:

In this sense, PCA may be understood as the zero-noise limit of PPCA. The latent subspace remains, but the probabilistic spread around it vanishes. PPCA therefore contains PCA as a special limiting case while adding something essentially new: uncertainty around the reconstruction.

In the next section, we will see that PPCA remains unusually tractable: because the model is linear and Gaussian, the marginal likelihood \(p(\mathbf{x})\) can still be derived analytically.

Suppose we observe a dataset

The log-likelihood of the dataset under PPCA is

It is often convenient to divide by \(N\) and write this as an average log-likelihood,

which may also be expressed as an expectation under the empirical data distribution:

Here, \(q(\mathbf{x})\) simply means the empirical distribution that places equal mass on each observed data point. Thus, taking an expectation under \(q(\mathbf{x})\) is nothing more than averaging over the training samples, a notation commonly appearing in more general probabilistic setting used later in the book.

The PPCA objective is therefore

Here the marginal likelihood depends on the model parameters \((\mathbf{W}, \boldsymbol{\mu}, \alpha)\). To avoid overly heavy notation, however, we will often suppress this dependence and write simply \(p(\mathbf{x})\) when the meaning is clear.

The key quantity here is the marginal likelihood \(p(\mathbf{x})\), which, by the law of total probability, is obtained by integrating out the latent variable:

At first glance, this marginal likelihood may appear difficult to compute. Substituting the PPCA model gives

Since \(\mathbf{z}\) ranges over a continuous latent space, it may seem that evaluating this integral exactly should be intractable. In most latent-variable models, that intuition is correct. PPCA is special precisely because its linear-Gaussian structure keeps this marginalisation analytically tractable.

These facts allow us to identify the distribution of \(\mathbf{x}\) without explicitly computing the integral by brute force.

Recall the PPCA generative equation:

with the latent variable \(\mathbf{z}\) and the noise term \(\boldsymbol{\epsilon}\) assumed to be independent. This shows that \(\mathbf{x}\) is constructed as an affine transformation of a Gaussian latent variable together with an independent Gaussian noise term. Because affine transformations of Gaussians remain Gaussian, and sums of independent Gaussians are again Gaussian, it follows that \(\mathbf{x}\) itself must also be Gaussian.

The same marginal can also be derived by expanding the joint density

and then integrating over \(\mathbf{z}\) using the so-called completing the square trick. This algebraic route is reviewed in Appendix A, where the same idea is used to show how Gaussian joint densities lead to Gaussian marginal and posterior distributions. For now, we use the simpler distributional argument based on affine transformations of Gaussian variables.

We can now work out its mean and covariance directly.

3.4.2.1 Mean.

3.4.2.2 Covariance.

3.4.2.3 Interpreting the covariance structure.

Because the marginal distribution of \(\mathbf{x}\) is Gaussian, PPCA admits an unusually elegant route to learning. In principle, one may write down the log-likelihood explicitly, differentiate it with respect to the parameters, set the resulting derivatives to zero, and solve for the stationary point:

This is the classical maximum-likelihood strategy. Its attraction is clear: if the equations can be solved exactly, then the model parameters can be expressed directly in terms of the empirical covariance structure of the data.

For PPCA, this route succeeds because the model remains linear and Gaussian all the way through. The final solution can be written in closed form and is closely tied to the eigenstructure of the sample covariance matrix. This is one reason PPCA is such an important pedagogical model: it is a genuine latent-variable generative model that still admits an exact maximum-likelihood solution.

At the same time, this route also reveals a limitation. Although elegant, direct analytic optimisation depends strongly on the model being simple enough for the derivatives to be solved exactly. As soon as the latent-variable model becomes nonlinear or non-Gaussian, this closed-form route typically breaks down. Moreover, even when an analytic solution exists, obtaining it may still require computationally expensive operations such as matrix inversion or eigen-decomposition, especially in high dimensions.

A second possibility is to treat the log-likelihood simply as an objective function and optimise it numerically by gradient ascent. In this view, the parameters are updated iteratively according to

where \( \theta = \{\mathbf{W}, \boldsymbol{\mu}, \alpha\},\) and \(\eta\) is a step size.

This is the most general optimisation viewpoint. Unlike direct analytic maximum likelihood, it does not require that the optimal parameters themselves be solvable in closed form. Instead, it treats PPCA like any other differentiable probabilistic model and uses gradients to climb the likelihood surface iteratively.

In PPCA, the marginal likelihood \(p(\mathbf{x})\) is available explicitly, so its gradient can also be computed explicitly. The difference from the analytic route is therefore not that gradients avoid the likelihood itself, but that they avoid having to solve the optimisation problem exactly in closed form. More generally, even when no exact closed-form solution for the parameters exists, gradient-based methods can still be used whenever the objective or its gradient can be computed—or at least approximated.

This perspective is worth emphasising because it foreshadows how most modern generative models are trained. In later chapters, when closed-form solutions are no longer available, gradient-based optimisation will become the default computational tool.

However, in the particular case of PPCA, gradient ascent is not the most illuminating route. It can certainly be applied, but it does not exploit the special latent-variable structure of the model. Nor does it reveal why latent-variable models are often approached through posterior inference rather than direct optimisation alone. For that, we need a third viewpoint.

The third route is the Expectation–Maximisation (EM) algorithm [10]. Unlike direct analytic optimisation, EM does not try to solve the likelihood in a single step. Unlike generic gradient ascent, it does not treat the latent variables as merely hidden nuisance variables inside a black-box objective. Instead, it works with the latent-variable structure of the model.

The basic idea is simple. Although the marginal likelihood

may be difficult to optimise directly, the joint distribution

often has a much simpler form. EM exploits this by alternating between two tasks:

• estimating the latent variables under the current parameters;
• updating the parameters using those estimates.

In PPCA, EM is not strictly necessary, because the model can already be solved analytically. But it is still extremely valuable. It reveals a general optimisation principle for latent-variable models, and that principle survives even when exact solutions do not.

Seen from a broader perspective, PPCA is especially valuable because it is one of the few latent-variable models in which all of the central probabilistic quantities can be written down explicitly: the prior \(p(\mathbf{z})\), the conditional distribution \(p(\mathbf{x}\mid\mathbf{z})\), the marginal likelihood \(p(\mathbf{x})\), and the posterior \(p(\mathbf{z}\mid\mathbf{x})\). In more complex models, these same quantities remain central, but they no longer stay analytically tractable. In particular, learning is driven by the marginal likelihood \(p(\mathbf{x})\), while inference depends on the posterior \(p(\mathbf{z}\mid\mathbf{x})\). This is precisely why EM is so important: it provides a structured way to alternate between these two aspects of the latent-variable problem.

More importantly for what follows, EM naturally leads to two deeper ideas: the introduction of an auxiliary distribution over latent variables, and the replacement of a difficult log-likelihood by a lower bound that can be optimised iteratively. These ideas are not unique to PPCA. They belong to the general theory of latent-variable modelling, and they will reappear in a broader form when we later study variational inference.

For this reason, we now turn to the mathematical foundation underlying EM: Jensen’s inequality and the Evidence Lower Bound (ELBO).

The latent-variable structure of PPCA is not unique to PPCA itself. More generally, many generative models assume that data are produced in two stages:

1. sample a latent variable \(\mathbf{z}\) from a prior \(p(\mathbf{z})\);
2. generate an observation \(\mathbf{x}\) from a conditional distribution \(p(\mathbf{x}\mid\mathbf{z};\theta)\).

The corresponding marginal likelihood is

For a dataset \(\{\mathbf{x}^{(i)}\}\), maximum likelihood learning becomes

The difficulty lies in the form of this expression: the latent variable must be integrated out, and the logarithm sits outside the integral. In PPCA, the linear-Gaussian structure allowed us to bypass this obstacle. In general, we need a more systematic tool.

That tool begins with Jensen’s inequality. Recall that a convex function curves upward, so that the straight line joining any two points on its graph lies above the graph itself. A concave function curves downward, so that the straight line joining any two points lies below the graph. For a random variable \(X\) and a convex function \(f\), Jensen’s inequality states

If \(f\) is concave, the direction of the inequality reverses:

In the present context, the function \(f\) we care about is specifically the logarithm. Since \(\log\) is concave,

This is precisely the form we need. In latent-variable models, the problematic quantity is typically

or equivalently \(\log \mathbb{E}[\cdot]\), where the logarithm sits outside an integral or expectation. Jensen’s inequality allows us to move this logarithm inside the expectation, trading exactness for a lower bound that is much easier to work with.

To apply Jensen’s inequality, we first rewrite the marginal likelihood as an expectation. Let \(Q(\mathbf{z})\) be any probability distribution over the latent variable, so that

Then

Taking logarithms gives

At first glance this may look like a purely algebraic trick, but it is an essential step. The distribution \(Q(\mathbf{z})\) is a freely chosen auxiliary distribution, introduced so that the marginal likelihood can be expressed as an expectation. Once this has been done, Jensen’s inequality becomes available.

We are now ready to apply Jensen’s inequality to the rewritten marginal likelihood:

Since the logarithm is concave, Jensen’s inequality gives

Hence,

This lower bound is called the Evidence Lower Bound, or ELBO:

Expanding the logarithm of the fraction gives the equivalent form

This expression is much easier to work with because the logarithm now acts directly on the joint density inside an expectation, rather than sitting outside an intractable integral. As a result, it becomes much more amenable to iterative optimisation.

This is the central result of the section. Instead of maximising the generally intractable quantity \(\log p(\mathbf{x};\theta)\) directly, we may instead maximise a tractable lower bound that depends on two objects: the model parameters \(\theta\) and the auxiliary distribution \(Q(\mathbf{z})\).

A lower bound is useful only if we understand when it is tight. For Jensen’s inequality, equality holds when the argument of the concave function—in our case, the logarithm—is constant with respect to the averaging variable. Here, that averaging variable is \(\mathbf{z}\), so equality holds precisely when

for some constant \(c\) independent of \(\mathbf{z}\).

Rearranging this expression gives \( Q(\mathbf{z}) \propto p(\mathbf{x},\mathbf{z};\theta).\) To determine the normalising constant, integrate both sides with respect to \(\mathbf{z}\). Since \(Q(\mathbf{z})\) is a valid probability distribution, \( \int Q(\mathbf{z})\,d\mathbf{z} = 1.\)

Therefore,

So the constant must be

Substituting this back into the equality condition yields

which is exactly the true posterior distribution of the latent variable.

Hence the ELBO becomes equal to the true log-likelihood exactly when \( Q(\mathbf{z})=p(\mathbf{z}\mid\mathbf{x};\theta),\) that is, when the auxiliary distribution matches the true posterior over the latent variable. This is a crucial point. The ELBO is not merely an arbitrary surrogate objective. Its tightness is controlled by how well \(Q(\mathbf{z})\) approximates the true posterior.

This perspective also clarifies the role of PPCA in the chapter. PPCA is not special because it avoids the latent-variable problem, but because it solves it exactly. Its posterior \(p(\mathbf{z}\mid\mathbf{x})\) is tractable, and its marginal likelihood \(p(\mathbf{x})\) can also be computed analytically. For that reason, PPCA does not strictly require the ELBO.

However, the ELBO framework reveals the more general principle that PPCA exemplifies. The same latent-variable structure persists in more complex models, but the exact posterior and exact marginal likelihood are no longer available. The ELBO is the device that allows us to keep going when exact Gaussian calculations break down.

We are now in a position to introduce the Expectation–Maximisation algorithm. The ELBO depends on two coupled ingredients: the auxiliary distribution \(Q(\mathbf{z})\) and the model parameters \(\theta\). This naturally suggests an alternating strategy: first improve \(Q(\mathbf{z})\) so that the bound becomes as tight as possible, and then improve \(\theta\) so that the bound itself becomes larger.

That alternating logic is exactly what EM formalises. In the next section, we will show how EM arises as an iterative procedure for maximising the ELBO, and how, in the special case of PPCA, it recovers the same solution that we previously obtained through exact Gaussian calculations.

Each iteration of EM consists of two stages.

3.7.1.1 E-step (Expectation step).

3.7.1.2 M-step (Maximisation step).

3.7.1.3 Why EM does not decrease the likelihood.

In PPCA, the latent posterior is tractable, so EM can be carried out explicitly. This makes PPCA an especially useful example: it shows how a general ELBO-based algorithm reduces, in a simple linear-Gaussian model, to closed-form updates.

For PPCA, the posterior distribution of the latent variable is Gaussian. This follows from the same linear-Gaussian structure used to derive the marginal distribution, and the completing-the-square derivation is given in Appendix A. The posterior is

From this posterior, the E-step provides the expectations

and

The M-step then substitutes these expectations into the expected complete-data log-likelihood and maximises with respect to the parameters \((\mathbf{W},\boldsymbol{\mu},\alpha)\). This leads to closed-form update equations for the PPCA parameters.

The important point here is not merely that PPCA can be fitted by EM, but that the exact PPCA solution may be reinterpreted as the fixed point of an EM procedure. In other words, EM does not give a different model; it gives a different view of learning the same model.

More broadly, EM teaches us that learning in latent-variable models can be understood as an alternation between

• inference, where we estimate hidden variables given the data;
• learning, where we update the model parameters using those inferred latent variables.

Although PPCA introduces a prior over latent variables, \( p(\mathbf{z})=\mathcal{N}(\mathbf{0},\mathbf{I}),\) and defines a valid generative model, it does not yet provide the kind of flexible and regularised inference framework that will become central in modern latent-variable models such as the VAE.

For any individual observation \(\mathbf{x}\), the posterior \( p(\mathbf{z}\mid \mathbf{x})\) describes which latent coordinates are plausible explanations for that data point under the model. If we average these posteriors over the empirical data distribution, we obtain the aggregate posterior: \( \mathbb{E}_{q(\mathbf{x})}\left[p(\mathbf{z}\mid \mathbf{x})\right].\) This distribution describes how the dataset as a whole is represented in latent space.

In PPCA, this posterior can be computed exactly because the model is linear and Gaussian. This is a strength of PPCA, but it also reveals the source of its limitation: the clean mathematics depends on strong assumptions. The latent-to-observed mapping is linear, the observation noise is Gaussian, and posterior inference can be worked out analytically.

This distinction becomes important when we move beyond the linear-Gaussian world. Modern generative models aim to keep the same latent-variable idea while making the mapping from latent variables to observations much more flexible. Once this flexibility is introduced, the simple closed-form calculations that made PPCA so transparent are usually no longer available. We therefore need a new framework that can learn expressive generative mappings, infer latent representations from data, and keep the latent space organised in a way that supports generation.

This is the role played by the Variational Autoencoder. The VAE can be viewed as a nonlinear and variational extension of the PPCA story: it preserves the idea of latent variables generating observations, but replaces exact Gaussian calculations with a learnable inference-and-generation framework. The details of this framework, including the precise role of the variational posterior and the latent-space regularisation term, are developed in the next chapter.

Thus, the transition from PPCA to the VAE is not a move from a non-generative model to a generative one. PPCA is already generative. Rather, the transition is from an analytically tractable linear-Gaussian latent-variable model to a more flexible nonlinear latent-variable model in which inference and generation must be learned together.

In the next chapter, the Variational Autoencoder develops this idea fully by combining reconstruction, probabilistic inference, and latent-space regularisation within a single variational objective.

In a deterministic autoencoder, the encoder is a function

For each observation \(\mathbf{x}\), it returns exactly one latent vector \(\mathbf{z}\).

In a VAE, the encoder instead defines a tractable conditional distribution over latent variables,

which serves as a learned approximation to the true posterior. This distribution is often called the variational posterior or approximate posterior. Its role is to describe which latent variables are plausible given the observed data point \(\mathbf{x}\).

The word posterior is used deliberately. In the previous PPCA chapter, we saw that in a latent-variable model the true posterior is

the distribution over latent variables after observing \(\mathbf{x}\). In the VAE, this true posterior is generally intractable. The main reason is that the encoder and decoder are usually modelled by neural networks with multiple nonlinear layers, so the posterior no longer retains the simple linear-Gaussian form that made PPCA analytically solvable. We therefore introduce the tractable approximation \(q_{\phi}(\mathbf{z}\mid \mathbf{x})\) instead. For now, the key structural change is simply this: the encoder no longer outputs a single vector, but the parameters of a distribution.

In practice, the encoder is typically implemented by a neural network that maps \(\mathbf{x}\) to the parameters of a Gaussian distribution. Most commonly, we assume

or, in a simplified isotropic presentation,

Thus the encoder network does not produce \(\mathbf{z}\) directly. Instead, it produces the mean and variance parameters that specify a Gaussian distribution in latent space.

Conceptually, the encoder now has the form

A latent code is then obtained by sampling

This makes the encoder stochastic. Even for the same input \(\mathbf{x}\), different samples of \(\mathbf{z}\) may be drawn from the same latent distribution. The representation is therefore no longer a single fixed point, but a small region of latent space together with a probability law over that region.

A similar shift occurs on the decoding side. In a deterministic autoencoder, the decoder is a function

which produces a single reconstruction \(\hat{\mathbf{x}}\) for each latent input.

In a VAE, the decoder instead defines a conditional distribution

which describes how an observation is generated from a latent variable. This is the generative part of the model.

Again, the decoder is typically parameterised by a neural network. Given \(\mathbf{z}\), the decoder outputs the parameters of a distribution in data space. For continuous data, a common choice is a Gaussian distribution such as

while for binary data one often uses a Bernoulli likelihood. The exact choice depends on the type of data being modelled, but the underlying principle is always the same:

the decoder does not output a single deterministic reconstruction, but a distribution over possible observations given the latent variable.

In practice, one often uses the mean or mode of this distribution as the model’s representative reconstruction, but conceptually the decoder defines a full conditional probability law rather than a single point prediction. So, whereas the encoder maps data to a distribution in latent space, the decoder maps latent variables to a distribution in data space.

4.3.3.1 The latent prior.

4.3.3.2 Why this creates a new difficulty.

In Chapter 2, the ELBO was derived by introducing an auxiliary distribution \(Q(\mathbf{z})\) over the latent variable. In the VAE, that auxiliary distribution now becomes a learned, parameterised approximation to the true posterior:

where \(q_{\phi}(\mathbf{z}\mid \mathbf{x})\) is the encoder network introduced in the previous section.

This substitution is the key step that turns the general ELBO framework into the VAE objective. Instead of trying to compute the intractable true posterior \(p_{\theta}(\mathbf{z}\mid \mathbf{x})\), we use the tractable approximation \(q_{\phi}(\mathbf{z}\mid \mathbf{x})\) and optimise a lower bound on \(\log p_{\theta}(\mathbf{x})\).

Applying the ELBO construction gives

Thus the VAE objective for a single data point is

This is the ELBO specialised to the VAE setting. It depends on both sets of parameters:

• \(\theta\), the parameters of the decoder network \(p_{\theta}(\mathbf{x}\mid \mathbf{z})\);
• \(\phi\), the parameters of the encoder network \(q_{\phi}(\mathbf{z}\mid \mathbf{x})\).

Training the VAE means maximising this lower bound jointly with respect to both parts of the model.

We will now rewrite the same objective in three equivalent forms. Each form highlights a different aspect of the learning problem: the latent-variable structure, the gap to the exact likelihood, and the practical balance between reconstruction and regularisation.

4.4.2.1 First form of the ELBO: joint-density form.

4.4.2.2 Second form of the ELBO: exact likelihood minus posterior mismatch.

4.4.2.3 Third form of the ELBO: reconstruction plus regularisation.

4.4.2.4 What the ELBO achieves.

The key idea is to separate the source of randomness from the dependence on the encoder parameters. For the common case of a Gaussian approximate posterior with diagonal covariance,

we can generate a latent sample by first drawing noise from a standard normal distribution

and then defining

This is the reparameterisation trick [13, 14].

The important point is that the randomness now lies entirely in the auxiliary noise variable \(\boldsymbol{\epsilon}\), which is independent of the parameters. Once \(\boldsymbol{\epsilon}\) has been sampled, the latent variable \(\mathbf{z}\) becomes a deterministic differentiable function of \(\mathbf{x}\), \(\phi\), and \(\boldsymbol{\epsilon}\):

This allows gradients to pass through the encoder in the usual backpropagation framework.

Conceptually, the reparameterisation trick does not remove stochasticity. Rather, it relocates it. Instead of sampling \(\mathbf{z}\) directly from a parameter-dependent distribution, we sample parameter-free noise and then transform it deterministically. This is what makes stochastic latent-variable learning compatible with standard gradient-based deep learning.

With the reparameterisation trick in place, the VAE can be trained end-to-end by maximising the ELBO jointly with respect to both parameter sets:

In practice, for each input the encoder produces \(\boldsymbol{\mu}_{\phi}(\mathbf{x})\) and \(\boldsymbol{\sigma}_{\phi}(\mathbf{x})\), a latent sample is formed via reparameterisation, the ELBO is evaluated, and gradients are then backpropagated through both encoder and decoder. This optimisation is typically carried out using minibatch stochastic gradient methods or one of their modern variants.

This is one of the most important conceptual achievements of the VAE. It turns approximate posterior inference and generative learning into a single differentiable computation graph. From this point onward, the encoder and decoder are no longer trained as separate objects: they are coupled through the ELBO and learned together.

In this sense, training a VAE is an instance of variational inference. Rather than computing the true posterior \(p_{\theta}(\mathbf{z}\mid\mathbf{x})\) exactly, we learn a tractable approximation \(q_{\phi}(\mathbf{z}\mid\mathbf{x})\) and optimise the ELBO jointly with respect to both the decoder parameters \(\theta\) and the encoder parameters \(\phi\).

The balance between reconstruction and regularisation can be adjusted in different ways, and one important extension is the \(\beta\)-VAE [15], which introduces a weighting factor on the KL term:

When \(\beta>1\), the latent space is more strongly regularised, often leading to more structured or disentangled representations, but usually at the cost of reconstruction quality.

This extension makes especially clear a central tension in variational latent-variable models:

the latent space must be regular enough to support generation, but informative enough to support useful representation and reconstruction.

Seen from this perspective, the VAE marks an important step in the development of modern generative modelling. PPCA introduced a tractable probabilistic latent-variable model. The VAE carries that same latent-variable logic into the nonlinear neural-network setting by replacing exact inference with a learned approximation and by turning the ELBO into a practical training objective.

These ideas will continue to matter in later chapters. In particular, although diffusion models are often introduced from other viewpoints, ELBO-based reasoning and latent-variable perspectives also provide a useful framework for approaching them. In the next chapter, we will see how these same probabilistic ideas can be extended further in one of the most important modern generative AI frameworks.

5.3.1.1 A step-by-step Gaussian noising process.

5.3.1.2 A closed-form link between \( \mathbf{x}_t\) and \( \mathbf{x}_0\) .

Because the diffusion process introduces a whole sequence of latent states, the full forward trajectory is described by the joint distribution

By the chain rule of probability,

The crucial simplifying assumption is that the forward process is Markovian. This means that each new state depends only on the immediately preceding one:

Hence the full forward distribution factorises as

This makes the forward process mathematically tractable and easy to analyse. But it is not the main object of learning. Once the data sample \(\mathbf{x}_0\) and the noise schedule \(\{\beta_t\}_{t=1}^{T}\) are fixed, the forward process \(q\) is fully known.

What ultimately matters for generation, however, is the reverse or denoising process: starting from a random Gaussian point and gradually reconstructing a data-like sample. That reverse process is not known analytically and must be learned, typically by a neural network. To understand it, we will now bring back the tools of ELBO and variational inference developed earlier and adapt them to this time-dependent latent sequence, so that the reverse process can be expressed, analysed, and eventually trained.

To describe generation probabilistically, we introduce a parameterised reverse process. By the general chain rule of probability, the joint distribution over the whole reverse trajectory could in principle be written as

This is the fully general expansion: each earlier state may depend on all later states.

In diffusion models, however, we make the same simplifying choice as in the forward process and assume that the reverse chain is also Markovian. That is,

Under this first-order Markov assumption, the reverse joint simplifies to

This expression should be read from right to left. We begin from a simple prior at the final time step, typically

and then apply a sequence of learned reverse transitions \(p_\theta(\mathbf{x}_{t-1}\mid \mathbf{x}_t)\) until we eventually obtain a sample at time \(0\).

This is the key unknown object in the model. Unlike the forward process \(q\), which was fully specified once the noise schedule \(\{\beta_t\}\) was chosen, the reverse transitions are not known in closed form and must therefore be learned from data. In practice, the parameters \(\theta\) belong to a neural network that will be trained to model these reverse denoising steps.

At first sight, one might imagine simply reversing the forward rule. But this is not straightforward. Although the forward transition

is Gaussian and analytically defined, the corresponding reverse transition

is not directly available in a usable form.

To see why, recall Bayes’ rule:

The forward conditional \(q(\mathbf{x}_t\mid \mathbf{x}_{t-1})\) is known by construction. The difficulty lies in the marginal distributions \(q(\mathbf{x}_{t-1})\) and \(q(\mathbf{x}_t)\). These are not simple fixed Gaussians chosen by hand; rather, they are induced by taking the unknown data distribution \(q(\mathbf{x}_0)\) and pushing it through the forward noising chain. So, although the forward step itself is simple, the true reverse conditional depends implicitly on the underlying data distribution and is therefore not available as a simple rule that we can directly evaluate or sample from.

A second difficulty appears when we consider likelihood-based learning. The model defines a joint distribution over the whole reverse trajectory \(\mathbf{x}_{0:T}\), but what we ultimately care about is the marginal likelihood of the observed clean data:

This is the diffusion-model analogue of the marginal likelihoods we encountered in earlier latent-variable models. The observed variable is the clean sample \(\mathbf{x}_0\), while the whole trajectory \(\mathbf{x}_{1:T}\) plays the role of latent variables and must be integrated out.

This quantity \(p_\theta(\mathbf{x}_0)\) is the key object we would like to learn and optimise, because it measures how much probability the model assigns to a real data sample. If a generative model assigns high probability to the observed data, it means it has learned a good distribution over the data domain. But in diffusion, this marginal likelihood is defined only indirectly through the reverse chain

So learning the reverse transitions \(p_\theta(\mathbf{x}_{t-1}\mid \mathbf{x}_t)\) well is precisely how we hope to build a model that assigns high likelihood to \(\mathbf{x}_0\).

The difficulty, however, is that computing \(p_\theta(\mathbf{x}_0)\) requires integrating over the entire latent trajectory \(\mathbf{x}_{1:T}\). Since this trajectory is high-dimensional and time-indexed, the integral is generally intractable.

So, just as in the VAE, the central learning problem is:

how can we optimise the model when the exact marginal likelihood is not available in closed form?

The answer is again to use an Evidence Lower Bound (ELBO).

5.4.3.1 Recasting the ELBO in the diffusion setting

5.4.3.2 The diffusion ELBO.

5.5.1.1 Substituting the forward and reverse chains.

5.5.1.2 What should the learned denoising process approximate?

To bring this true reverse posterior into the ELBO, we apply Bayes’ rule to the forward transition:

Substituting this into the denominator of the ELBO gives

Now split the logarithm of the product into a sum:

The second sum telescopes. To see this explicitly,

because all intermediate terms cancel pairwise. The term \(q(\mathbf{x}_0\mid \mathbf{x}_0)\) is formally a degenerate distribution concentrated at the observed value \(\mathbf{x}_0\), corresponding to a Dirac delta rather than an ordinary density. Since it contains no learnable parameters, it contributes only an additive constant to the objective and may be ignored for optimisation purposes. So the ELBO becomes

We now reorganise the ELBO into the standard step-wise decomposition. There are two main pieces to handle:

• the boundary term involving the endpoint \(\mathbf{x}_T\);
• the time-step terms involving the learned denoising transitions \(p_\theta(\mathbf{x}_{t-1}\mid \mathbf{x}_t)\).

Starting from the expression obtained above, we first isolate the part involving only the endpoint \(\mathbf{x}_T\):

Since this integrand depends only on \(\mathbf{x}_T\), all the other latent variables \(\mathbf{x}_1,\dots,\mathbf{x}_{T-1}\) can be marginalised out under the forward chain, leaving only the marginal \(q(\mathbf{x}_T\mid \mathbf{x}_0)\):

This is the prior-matching term: it compares the distribution reached by the forward process at its final time step with the simple Gaussian prior from which the learned reverse process begins.

We now turn to the sum over \(t\). Here one must handle the inner expectation \(\mathbb{E}_{q(\mathbf{x}_{1:T}\mid \mathbf{x}_0)}\) carefully. The whole sum sits inside the expectation, but each summand

depends only on the local variables \(\mathbf{x}_{t-1}\) and \(\mathbf{x}_t\), together with the observed sample \(\mathbf{x}_0\). So for each fixed \(t\), all the other latent variables can again be marginalised out.

For the special case \(t=1\), this gives

Since \(q(\mathbf{x}_0\mid \mathbf{x}_1,\mathbf{x}_0)\) is degenerate once \(\mathbf{x}_0\) is conditioned on, it corresponds to a point mass at the observed value and contains no learnable parameters. It therefore contributes only an additive constant, leaving the reconstruction term

For the remaining terms \(t=2,\dots,T\), we do the same localisation of the inner expectation:

Now recognise the inner expectation as the negative of a KL divergence:

Therefore each intermediate time-step contribution becomes

Putting the boundary term and the time-step terms together, the negative ELBO can be written as

This decomposition shows exactly where the learning problem comes from:

• the prior-matching term ensures that the endpoint of the forward process is consistent with the Gaussian prior from which generation begins;
• the step-wise denoising terms compare the learned reverse transition to the true reverse posterior at each intermediate time step;
• the reconstruction term handles the final recovery of the clean data from the first noisy state.

The key consequence of this decomposition is that learning reduces to matching, at each time step, the learned denoising transition \( p_\theta(\mathbf{x}_{t-1}\mid \mathbf{x}_t)\) to the true reverse posterior \( q(\mathbf{x}_{t-1}\mid \mathbf{x}_t,\mathbf{x}_0).\) The next step is therefore to work out the analytical form of this posterior.

From the previous section, we have the identity

This is useful because every term on the right-hand side is determined by the forward process, which is already known. In particular,

where \(\alpha_t=1-\beta_t\). From the earlier closed-form forward relation, we also know

and

So the true reverse posterior is expressed entirely in terms of Gaussian quantities from the forward diffusion process. Because the terms above are Gaussian, the posterior in \(\mathbf{x}_{t-1}\) must also be Gaussian. So we may write

for some mean \(\tilde{\boldsymbol{\mu}}_t(\mathbf{x}_t,\mathbf{x}_0)\) and variance coefficient \(\tilde{\beta}_t\).

Using the same completing-the-square idea introduced earlier, and derived in Appendix B, the true reverse posterior takes the following Gaussian form: with

and

This result is worth interpreting. The posterior mean \(\tilde{\boldsymbol{\mu}}_t(\mathbf{x}_t,\mathbf{x}_0)\) is a weighted combination of two pieces of information: the clean sample \(\mathbf{x}_0\), which tells us where the trajectory began, and the current noisy state \(\mathbf{x}_t\), which tells us where we are at time \(t\). So the true reverse posterior identifies the most plausible slightly cleaner state by combining knowledge of both the clean origin and the current noisy observation.

The posterior variance \(\tilde{\beta}_t\mathbf{I}\) is also fixed by the forward diffusion schedule. That is, once the forward noising process is chosen, the uncertainty of the corresponding reverse step is determined as well. So, from the ELBO viewpoint, the reverse model is not being asked to match an arbitrary denoising distribution, but a specific Gaussian target induced by the known forward process.

Although this posterior is exact, it still depends explicitly on the clean data sample \(\mathbf{x}_0\). That is acceptable during training, when \(\mathbf{x}_0\) is known, but it is not usable directly at generation time, when \(\mathbf{x}_0\) is precisely the unknown object we wish to create.

The key observation is that the forward process already provides a direct relationship between the clean sample \(\mathbf{x}_0\), the noisy sample \(\mathbf{x}_t\), and a standard Gaussian noise variable:

Solving this for \(\mathbf{x}_0\) gives

We now substitute this into the posterior mean \(\tilde{\boldsymbol{\mu}}_t(\mathbf{x}_t,\mathbf{x}_0)\) in order to remove its explicit dependence on \(\mathbf{x}_0\).

Starting from

substituting the expression for \(\mathbf{x}_0\) yields

After simplifying the coefficients, this becomes

So the true reverse posterior may be rewritten as

This form is extremely important. It shows that the reverse step can be characterised through the current noisy sample \(\mathbf{x}_t\) and the aggregate noise \(\boldsymbol{\epsilon}\) that produced it. So instead of trying to recover \(\mathbf{x}_0\) directly, the model can aim to predict the underlying noise.

This leads to the standard DDPM parameterisation. Instead of using the true noise \(\boldsymbol{\epsilon}\), which is unknown to the model, we introduce a neural network

to predict it from the noisy sample \(\mathbf{x}_t\) and the time step \(t\).

The learned reverse transition is then defined as

with mean

This has a clear interpretation: if the network predicts the true noise accurately, then the learned reverse mean matches the exact reverse posterior mean, and so the learned denoising process matches the correct local reverse step.

Now recall that the central trainable terms in the ELBO are

Now both distributions are Gaussian with the same variance \(\tilde{\beta}_t\mathbf{I}\). So their KL divergence reduces to a squared difference between means:

Substituting the two mean expressions gives

So, up to a time-dependent scaling factor, the KL term becomes a mean-squared error in noise prediction:

This is the central simplification in DDPMs. A complicated latent-variable ELBO ultimately reduces to a regression problem:

given a noisy sample \(\mathbf{x}_t\) and time step \(t\), predict the Gaussian noise that was used to produce it.

The commonly used practical DDPM objective [17] is therefore

where \( \mathbf{x}_t = \sqrt{\bar{\alpha}_t}\,\mathbf{x}_0 + \sqrt{1-\bar{\alpha}_t}\,\boldsymbol{\epsilon}.\)

The training loop now becomes conceptually simple. For each clean data sample \(\mathbf{x}_0\), we:

1. sample a time step \(t\);
2. sample Gaussian noise \(\boldsymbol{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I})\);

3. construct the noisy state

   \[ \mathbf{x}_t = \sqrt{\bar{\alpha}_t}\,\mathbf{x}_0 + \sqrt{1-\bar{\alpha}_t}\,\boldsymbol{\epsilon}; \]
4. train the network to predict \(\boldsymbol{\epsilon}\) from \((\mathbf{x}_t,t)\).

Because \(\mathbf{x}_t\) can be sampled directly from \(\mathbf{x}_0\), training does not require simulating the full forward chain step by step.

Once trained, the model generates data by running the learned reverse process iteratively. We begin from

and then for \(t=T,T-1,\dots,1\) sample

At each step, the model uses its predicted noise \(\boldsymbol{\epsilon}_\theta(\mathbf{x}_t,t)\) to compute the reverse mean and remove a small amount of corruption. After repeating this over all time steps, the final state \(\mathbf{x}_0\) is obtained as the generated sample.

This completes the main probabilistic story of DDPMs in discrete time. The forward process is fixed, the reverse process is learned, the ELBO reduces to local denoising targets, and the final practical objective becomes noise prediction.

We begin by introducing the basic language used to describe motion in continuous time. The aim is to make precise what it means for a state to evolve smoothly, how its local change is measured, and how discrete update rules can be reinterpreted in infinitesimal form.

Let \(x(t)\) denote the position of a quantity at time \(t\). In one dimension, \(x(t)\) is simply a scalar-valued function of time. In multiple dimensions, it becomes a vector whose coordinates all evolve with time. The first derivative, \(\frac{dx}{dt}\), measures the instantaneous rate of change of \(x\) with respect to time, and is commonly interpreted as velocity. It is also often written in dot notation as \(\dot{x}(t)\). Similarly, the second derivative, \(\frac{d^2x}{dt^2}=\ddot{x}(t)\), measures the rate of change of velocity, and is commonly interpreted as acceleration.

The symbols \(dx\) and \(dt\) are often read informally as an infinitesimal change in position and an infinitesimal change in time. The basic intuition is that if \(\frac{dx}{dt}\) is the rate of change, then over an extremely small time interval, the corresponding change in position is represented by

So, in infinitesimal form, the local change is again understood as

Once we write equations describing how the position changes over time, we arrive at differential equations. In the simplest case, such an equation specifies the velocity of the state at each possible location as a function of its current position, and sometimes of time as well. In higher dimensions, this local motion rule can be viewed geometrically as assigning to each point in space a direction of motion. One may therefore think of a differential equation as a rule of motion, and of a vector field as the picture formed by all the local arrows indicating how a particle would move if placed at different locations.

These ideas are especially useful here because they allow us to reinterpret iterative update rules as continuous motion laws. Over a small time interval \(\Delta t\), a first-order approximation gives

This says that the value at a slightly later time is approximately equal to the current value plus velocity times a small time increment. It is therefore the continuous-time analogue of a discrete update rule.

We first consider deterministic motion, where the future evolution is fully specified once the current state is known. The simplest example is contraction: a motion in which the state is continuously pulled back toward the origin. This begins in one dimension and then extends naturally to multiple dimensions, where different coordinates may contract at different rates.

6.3.2.1 One-dimensional contraction.

6.3.2.2 Extension to multiple dimensions.

The discussion so far has been entirely deterministic: once the initial state is fixed, the future trajectory is completely determined by the governing differential equation. Many continuous-time generative models, however, also involve randomness. At each small time interval, the state may receive a random perturbation, and the resulting motion is no longer fully predictable from its initial condition alone.

The basic model of such random continuous-time fluctuation is Brownian motion, also called Wiener motion [20]. It is usually denoted by \(\{W_t\}_{t\ge 0}\), and may be characterised by three basic properties:

1. \(W_0=0\);
2. increments over disjoint time intervals are independent;

3. for any \(\Delta t>0\),

   \[ W_{t+\Delta t}-W_t \sim \mathcal{N}(0,\Delta t). \]

The third property is especially important. It says that over a short time interval of length \(\Delta t\), the random increment has mean zero and variance \(\Delta t\). So Brownian motion introduces fluctuations that are unbiased on average, but whose uncertainty grows with elapsed time.

The Gaussian assumption is also natural. If we imagine the motion over a short interval as the accumulated effect of many tiny, approximately independent random perturbations, then their total is well modelled by a Gaussian random variable. This is closely related to the intuition behind the central limit theorem: when many small independent effects are added together, their sum tends toward a Gaussian form.

A useful way to understand this is to divide the interval \([0,t]\) into many small pieces and write \(W_t\) as a sum of independent Gaussian increments. If we divide \([0,t]\) into \(n\) equal pieces of length \( \Delta t=\frac{t}{n},\) then

Since \(W_0=0\), this becomes

So at each fixed time \(t\), Brownian motion is centred at zero, but its spread grows linearly with time. This means that Brownian motion does not pull the state toward any preferred location; instead, it describes random wandering whose uncertainty increases as time passes.

A useful heuristic way to think about Brownian motion is through the approximation

This makes the scaling especially clear. Because the increment has variance \(\Delta t\), its typical size is of order \(\sqrt{\Delta t}\), not \(\Delta t\). This square-root scaling is one of the most characteristic features of stochastic continuous-time motion and will later reappear in diffusion models.

Brownian motion is often written informally in differential form as \( dW_t,\) which represents an infinitesimal random increment. Although \(dW_t\) is not an ordinary differential in the deterministic sense, this notation is extremely useful because it allows stochastic motion to be written in a form parallel to ordinary differential equations.

A particularly important example of stochastic continuous-time motion is the Ornstein–Uhlenbeck (OU) process. It combines deterministic contraction with Brownian noise:

This equation contains two parts. The drift term \( -\alpha X_t\,dt\) pulls the state back toward the origin, just as in the one-dimensional contraction model introduced earlier. The noise term \( \sigma\,dW_t\) injects random perturbations driven by Brownian motion. So the OU process balances two competing tendencies: deterministic pull toward the origin and stochastic fluctuations away from it.

Over a small time interval \(\Delta t\), one may write the corresponding update heuristically as

This discrete-time view is very useful. It shows that the next state is obtained from the current one by adding two pieces:

• a deterministic correction of size \(\Delta t\), and
• a random perturbation of size \(\sqrt{\Delta t}\).

So the OU process already has the basic drift-plus-noise structure that will later appear in Langevin dynamics and diffusion models.

The OU process also has an explicit solution:

This formula reveals its structure clearly. The first term, \( X_0e^{-\alpha t},\) is exactly the deterministic exponential contraction of the initial state. The second term accumulates all past random perturbations, but each perturbation is weighted by the decay factor \(e^{-\alpha(t-s)}\). So older noise contributions gradually lose influence, while more recent ones matter more.

This is why the OU process is often described as a mean-reverting stochastic process. The OU process shows how a system can remain random while still having a stable overall tendency, with the drift-plus-noise structure that will be one of the main templates for the continuous-time generative models developed later.

So far, however, we have still described motion only at the level of individual particles. The next step is to understand how such deterministic and stochastic motions reshape an entire probability density over time.

We begin with the simplest density principle: probability is conserved under deterministic transport. Here this means that particles are moved by the dynamics from one location to another in a fully determined way, while the associated probability mass is simply carried with them.

Suppose that at time \(0\), a particle has position \(x_0\), distributed according to density \(p_0(x_0)\). Suppose further that the deterministic dynamics transport this initial point to a later position

where \(\phi_t\) denotes the trajectory or flow map that moves the particle from time \(0\) to time \(t\).

To understand how densities change, consider a very small interval around \(x_0\), of length \(dx_0\). The probability mass contained in that tiny interval is approximately

After evolving to time \(t\), that same small piece of probability mass is transported to a new interval around \(x_t=\phi_t(x_0)\). Let the length of that transported interval be \(dx_t\). Since deterministic motion merely transports particles and does not create or destroy probability mass, the total probability must remain the same:

The key question is therefore how the small interval \(dx_0\) is transformed under the flow map. In one dimension, this local stretching or compression is measured by the derivative of the map. The derivation is as follows.

This is the one-dimensional change-of-variable formula induced by a deterministic flow. Its interpretation is very intuitive. If the transformation stretches an interval, then the same amount of probability is spread over a larger region, so the density decreases. If it compresses an interval, then the same probability mass is squeezed into a smaller region, so the density increases. Thus the derivative \(\frac{dx_t}{dx_0}\) measures how much local length has changed, and its reciprocal tells us how the density rescales.

The same idea extends naturally to multiple dimensions. In one dimension, we tracked how a small length changes under the map. In two dimensions, we must track how a small area changes; in three dimensions, a small volume; and in general, a small \(d\)-dimensional volume element. The quantity that measures this local scaling is no longer an ordinary derivative, but the determinant of the Jacobian matrix of the transformation. Thus the multidimensional change-of-variable formula becomes

So the determinant plays, in higher dimensions, the role that \(\frac{dx_t}{dx_0}\) plays in one dimension, measuring how much a tiny local region has been stretched or compressed by the flow.

It is important to note that this formula describes a finite-time deterministic transformation. It tells us how density changes between time \(0\) and time \(t\), but it does not yet describe the infinitesimal rate of density change. For that, we next derive a local-in-time continuity equation.

The finite-time change-of-variable formula tells us how a density transforms between time \(0\) and time \(t\). In continuous-time modelling, however, it is often even more important to understand how density changes over an infinitesimal time interval, because continuous dynamics is built from such local updates. This leads to the Liouville equation, also called the continuity equation in the deterministic setting.

Suppose particles move according to the one-dimensional ordinary differential equation

where \(f(x,t)\) is the local velocity field. To understand how the density changes, we now adopt a fixed-space viewpoint: rather than following one particle, we stand at a small interval in space and watch probability mass flow through it.

Consider a short interval \([x,x+\Delta x]\). The total probability mass inside this interval at time \(t\) is approximately

Over an infinitesimal time interval, this local mass can change only because probability flows in through the left boundary or flows out through the right boundary.

The relevant flow quantity is the probability flux, namely density times velocity:

This tells us how much probability mass crosses a location per unit time. Thus the inflow at the left boundary is

while the outflow at the right boundary is

Conservation of probability therefore says that the rate of change of probability mass inside the interval must equal inflow minus outflow:

Dividing both sides by \(\Delta x\) gives

Rearranging,

Now let \(\Delta x\to 0\). The right-hand side becomes a spatial derivative, so we obtain

This is the one-dimensional Liouville equation.

The equation says that density changes according to the net outflow of probability flux. If more probability flows out of a region than into it, the density decreases there. If more flows in than out, the density increases. It is often useful to expand the right-hand side using the product rule:

This reveals two distinct contributions.

The term \( -f\,\partial_x p\) depends on the spatial slope of the density itself. It therefore describes transport: an existing density profile is being carried through space by the velocity field. If the flow simply moves particles from one place to another without changing their relative spacing, then the main effect is that the whole density profile shifts.

By contrast, the term \( -p\,\partial_x f\) depends on the spatial slope of the velocity field. It therefore describes deformation: nearby particles no longer move at exactly the same speed, so they either spread apart or compress together. If \(\partial_x f>0\), the velocity increases with position, so nearby particles separate and the density decreases. If \(\partial_x f<0\), the velocity decreases with position, so nearby particles compress and the density increases.

The Liouville equation describes how density changes at a fixed location in space. But there is another viewpoint that is also very useful: how does the density change along the trajectory of a moving particle? This viewpoint shifts attention from what happens at a fixed position to what is experienced by a particle as it moves with the flow.

Suppose a particle follows the deterministic dynamics

Now consider the density evaluated along that particle path:

Because both the position and time are changing, we differentiate this quantity using the chain rule. Equivalently, this is the total derivative: it combines the explicit change with respect to time and the implicit change coming from the fact that the particle itself is moving through space. Thus

Since \(\dot{x}(t)=f(x(t),t)\), this becomes

We now substitute the Liouville equation

This gives

The transport terms cancel, leaving

Dividing both sides by \(p(x(t),t)\) gives

Hence

This is an important result. In the previous subsection, the density evolution was decomposed into a transport part and a deformation part. Here we see that, once we move with the particle itself, the transport part disappears: from the particle’s own point of view, it is no longer observing the profile being carried through space, because it is travelling with that flow. What remains is only the local expansion or compression effect.

Thus, if \(\partial_x f>0\), nearby particles tend to separate and the density along the trajectory decreases. If \(\partial_x f<0\), nearby particles tend to compress and the density increases. So the fixed-space and moving-particle viewpoints emphasise different aspects of the same dynamics. The Liouville equation describes how the whole density field changes across space, whereas the formula above tells us how density changes along an individual trajectory.

We now extend density evolution from one dimension to a full multidimensional state space. The main new feature is that both the motion and the probability flux become vector-valued. In one dimension, the scalar derivative \(\partial_x\) was enough to describe how probability flowed in or out of a small interval. In multiple dimensions, the corresponding quantity is the divergence, which measures the net outflow of a vector field from a tiny region of space. This section shows how divergence arises naturally from local volume change, and how it leads to the multidimensional Liouville equation.

Let \(x(t)\in\mathbb{R}^d\) and suppose particles evolve according to

where the velocity field is

As in the one-dimensional case, the basic idea is still local conservation of probability. But before we can write down the density equation, we first need to understand how a tiny region of space changes under the flow. In one dimension, this was measured by the derivative of the map, which told us how a small interval was stretched or compressed. In multiple dimensions, the corresponding objects are the Jacobian matrix and its determinant.

Over a very small time step \(t\to t+dt\), the motion is approximated by

So the infinitesimal transformation is

To understand how a tiny volume changes, we ask how the updated point depends on its starting location. That dependence is captured by differentiating the map with respect to \(x\). This gives the Jacobian matrix of the transformation:

If we denote the Jacobian matrix of the velocity field by

then

The Jacobian tells us how a small displacement changes under the flow. To pass from displacements to volumes, we take its determinant.

This gives the key geometric interpretation of divergence: it is the instantaneous rate at which a tiny volume expands or contracts under the flow. If \(\nabla\cdot v>0\), the flow locally expands volume. If \(\nabla\cdot v<0\), it locally compresses volume.

With that geometric picture in place, we can now derive the density-evolution equation. The probability mass inside a small moving volume must remain the same:

Substituting the volume scaling relation,

we obtain

Next, apply a multivariable first-order Taylor expansion to the density.

This is the multidimensional Liouville equation. It is the natural generalisation of the one-dimensional continuity equation. Instead of the scalar derivative \(\partial_x(pf)\), we now have the divergence of the vector-valued probability flux \(pv\). So the density changes according to the net outflow of probability from a small region of space.

The trajectory-wise viewpoint also extends naturally. Along a path satisfying

the log-density evolves according to

So in multiple dimensions, just as in one dimension, once we move with the particle itself, the transport part disappears and the remaining density change is governed by local expansion or compression of the flow.

To build intuition, consider first pure Brownian diffusion in one dimension:

Over a small interval \(\Delta t\), the increment is Gaussian:

Thus a particle previously at location \(y\) has some probability of arriving near \(x\) after this short time step, and that probability is described by a Gaussian kernel centred at \(y\). Therefore the updated density is obtained by collecting contributions from all nearby previous locations:

In other words, the density at \(x\) after a short time is a weighted average of the previous density around \(x\), where the weights are given by the short-time Gaussian transition kernel. This is a convolution with a narrow Gaussian kernel. Intuitively, probability mass diffuses into \(x\) from nearby locations, while some of the probability previously near \(x\) diffuses away. The net result is smoothing and spreading of the density.

To make this more convenient to analyse, let the Gaussian increment be denoted by \(\eta\), where

Introduce the change of variable

Then the density update may be written as

where

is the Gaussian density of the increment. Since

this can also be viewed as an expectation:

We now expand \(p(x-\eta,t)\) around \(x\) using Taylor series. In general, for a sufficiently smooth function of one variable, \( f(x+h)=\sum_{k=0}^{\infty}\frac{f^{(k)}(x)}{k!}h^k.\) So, taking \(f(\cdot)=p(\cdot,t)\) and \(h=-\eta\), we obtain

Taking expectation with respect to \(\eta\) gives

Because the Gaussian increment has mean zero,

the first-order term vanishes. Its variance is

so the second-order term survives:

Subtracting \(p(x,t)\) from both sides, dividing by \(\Delta t\), and taking \(\Delta t\to 0\), we obtain

Thus pure Brownian motion gives the pure diffusion equation.

This is the density-level signature of random spreading. Unlike the Liouville equation, which moves density through deterministic transport, the diffusion equation smooths and spreads the density through space.

Now suppose we also include deterministic drift:

Over a short interval, the particle undergoes two effects:

• a deterministic displacement of order \(\Delta t\) due to the drift \(f\),
• a random Gaussian displacement of order \(\sqrt{\Delta t}\) due to the Brownian term.

The resulting density evolution combines the deterministic Liouville transport term with the diffusion term derived above:

This is the one-dimensional Fokker–Planck equation for constant diffusion coefficient \(\sigma\).

In multiple dimensions, the deterministic part becomes

while the diffusion term becomes

where the Laplacian is

Thus, for constant isotropic diffusion coefficient \(\sigma\), the multidimensional Fokker–Planck equation takes the form

This equation is one of the most important density-evolution laws in continuous-time probabilistic modelling. The Liouville equation tells us how densities move under deterministic flow. The Fokker–Planck equation extends that picture by adding stochastic diffusion. Together, they provide the main bridge from elementary continuous-time motion to the density dynamics that underlie fully continuous diffusion models.

In the next chapter, we will build directly on this foundation. The key next step is to understand how density evolution can be linked to the gradient of log-density, and how that viewpoint leads naturally to score-based learning and fully continuous diffusion.

Let \(p(x)\) be a probability density on \(\mathbb{R}^d\). Its score function is defined as

If the density depends on time as well, then the time-dependent score is written as

At first sight, it may seem slightly surprising that the gradient of the log-density should be so important. The reason is that the score captures local directional information in a particularly useful form. The density \(p(x)\) tells us how much probability mass is concentrated near \(x\), whereas the score tells us in which direction the log-density increases most rapidly. It therefore provides a local vector field pointing toward regions of higher probability.

In one dimension, the score becomes

This expression is already revealing. The numerator \(p'(x)\) measures the local change in density, but dividing by \(p(x)\) normalises that change relative to the local scale of the density itself. The score therefore measures a relative slope rather than just an absolute slope. This makes it more meaningful as a field of probabilistic motion.

A simple Gaussian example makes this interpretation concrete. Consider the one-dimensional Gaussian density

Taking logarithms gives

Differentiating with respect to \(x\), we obtain

Hence the score is

This Gaussian example immediately shows the geometric meaning of the score. If \(x>\mu\), then the score is negative, so it points back toward the mean. If \(x<\mu\), then the score is positive, again pointing back toward the mean. The further \(x\) is from the mean, the larger the magnitude of the score. Thus, for a Gaussian density, the score behaves like a restoring field that directs particles toward the high-density region.

In the multivariate Gaussian case,

the score becomes

So again, the score is a vector field pointing toward the high-density centre, but now shaped by the covariance structure. Directions with smaller variance produce stronger restoring forces, while directions with larger variance produce weaker ones. This connects naturally to the contraction dynamics discussed earlier: in the Gaussian case, the drift toward the centre is directly given by the score of the density.

This suggests a natural idea. If we define a deterministic flow by

then particles move in the direction of steepest increase of log-density. Since the logarithm is monotone, this also moves particles toward regions of higher density. In this sense, the score may be interpreted as a probabilistic drift field: it tells us how to move if we wish to guide particles toward more likely regions under the target distribution.

This viewpoint is extremely useful, but it also reveals an important limitation. A purely deterministic score flow tends to move particles uphill toward modes. It may therefore collapse trajectories into concentrated high-density regions rather than explore the full target distribution. This is useful for optimisation, but not yet sufficient for sampling. To obtain a genuine sampling dynamics, we must add noise. This leads directly to Langevin dynamics.

Langevin dynamics combines two ingredients:

1. a drift term given by the score, which pulls particles toward higher-density regions;
2. a Brownian noise term, which injects randomness and maintains exploration.

In continuous time, the Langevin stochastic differential equation takes the form

where \(p^\ast(x)\) denotes the target density of interest. The notation \(p^\ast\) is used here to emphasise that this is the density we would like to sample from.

In discrete time, this corresponds heuristically to the update

where \(\epsilon>0\) is a small step size.

This update has a direct interpretation. The current state \(x_k\) is first moved by the score term

which nudges the particle toward regions of greater probability under the target density. Then a Gaussian perturbation

is added. This random term prevents simple collapse into a mode and allows the particle to keep exploring the state space. In the continuous-time limit, the deterministic score drift and the Brownian perturbation combine into the Langevin SDE above.

The structure should look familiar from the previous chapter. Langevin dynamics is a drift-plus-noise system, and therefore its density evolves according to a Fokker–Planck equation. For the Langevin dynamics

the drift field is

and the diffusion coefficient is \(\sigma=\sqrt{2}\). Since the Fokker–Planck equation for constant isotropic diffusion is

substituting \(v(x)=\nabla_x\log p^\ast(x)\) and \(\sigma^2=2\) gives

This equation combines two effects. The first term, \( -\nabla\cdot\bigl(p\nabla_x \log p^\ast\bigr),\) describes deterministic transport induced by the score field. The second term, \( \Delta p,\) describes diffusion caused by Brownian noise. The drift term tends to concentrate probability toward high-density regions of the target, while the diffusion term tends to spread probability outward. The important fact is that, with exactly this combination, the target density becomes stationary. This is the key reason Langevin dynamics is useful as a sampling mechanism.

It is also useful to compare this with the Ornstein–Uhlenbeck process introduced earlier. The OU process has a linear drift toward the origin and a Gaussian stationary distribution. Langevin dynamics generalises the same principle: instead of using a fixed linear restoring force, it uses the score of a chosen target density as the restoring field. In this sense, the OU process may be viewed as a Gaussian special case of a broader score-based sampling mechanism.

To see this explicitly, suppose the target density is the standard Gaussian \(p^\ast(x)=\mathcal{N}(0,I)\). Then

and Langevin dynamics becomes

which is precisely an Ornstein–Uhlenbeck process with a particular noise scaling. Thus the Langevin equation recovers the OU process when the target density is Gaussian, and extends it to non-Gaussian targets through the score \(\nabla_x\log p^\ast(x)\).

We now explain why Langevin dynamics is not merely a formal equation, but a meaningful sampling mechanism.

Suppose a density \(p(x,t)\) evolves according to some density evolution equation

A density \(p^\ast(x)\) is called a stationary distribution if, once the system has density \(p^\ast\), it remains there. In PDE language, this means

Equivalently,

This concept is important in generative modelling. If the target density is stationary for a dynamics, then the dynamics is aligned with that target distribution. If, in addition, trajectories or densities started elsewhere approach this stationary law over time, then the process becomes a mechanism for sampling from the target distribution. A full convergence theory requires additional assumptions and is beyond the present tutorial, but the stationarity calculation already reveals the central idea.

Let us now verify stationarity for Langevin dynamics. Recall that its Fokker–Planck equation is

We now substitute

Then the right-hand side becomes

Now use the identity

because

Therefore the previous expression simplifies to

But by definition,

Hence

So indeed,

Therefore \(p^\ast\) is a stationary distribution of the Langevin dynamics

This is one of the central conceptual results behind score-based generative modelling. It shows that if we use the score of a target density as the drift field, and combine it with the appropriate Brownian noise, then the resulting density evolution preserves that target density. The score is therefore not merely an interesting derivative of the density. It is precisely the drift field that makes the desired density stationary under a natural stochastic dynamics.

At the same time, it is important to be precise about what has been shown. The argument above establishes stationarity, not a full proof that every initial density converges to \(p^\ast\). The latter is a deeper question involving the long-time behaviour of the stochastic process. Nevertheless, the stationarity result closes an important conceptual gap. It explains why score-based drift, together with Brownian noise, provides a principled route toward sampling from distributions of interest.

If the Langevin process is allowed to evolve long enough under suitable conditions, its distribution may approach the stationary law \(p^\ast\). At the density level, this means that probability mass progressively takes on the shape of the desired target distribution. At the particle level, it means that repeatedly updating samples through Langevin dynamics can produce particles whose statistical behaviour resembles draws from the target law.

Thus the density-level and particle-level viewpoints are consistent with one another. The Fokker–Planck equation explains how the distribution evolves, while the Langevin update explains how individual samples move. The score links the two: it acts locally as a probabilistic drift field, and globally as the ingredient that makes the target density stationary.

We begin by expanding the squared norm:

The final term,

does not depend on \(\theta\), and therefore does not affect the optimisation. Thus minimising \(\mathcal{J}(\theta)\) is equivalent to minimising

The first term is directly computable from samples, because it only requires evaluating the model \(s_\theta(x)\) at data points. The challenge lies in the cross term

because it still involves the unknown true score.

The key idea of score matching is to rewrite this cross term using integration by parts. We first do this carefully in one dimension, where the calculation is most transparent, and then extend the result to multiple dimensions.

To make the main idea clear, suppose a one-dimensional scenario, with the true data score

Substituting this into the cross term gives

At this stage, the derivative acts on the unknown density \(p_{\text{data}}(x)\). Our goal is to transfer this derivative onto the model \(s_\theta(x)\) instead, using the integration by parts.

We now apply this identity by choosing

Then

Under the standard assumption that \(p_{\text{data}}(x)\) decays sufficiently fast as \(|x|\to\infty\), and that \(s_\theta(x)\) does not grow too rapidly, the product \(s_\theta(x)p_{\text{data}}(x)\) vanishes at the boundary. Hence the boundary contribution is zero. We therefore obtain

Recognising the right-hand side as an expectation under \(p_{\text{data}}\), we conclude that

The derivative has therefore been transferred from the unknown density \(p_{\text{data}}(x)\) to the model \(s_\theta(x)\), requiring no longer explicit knowledge of the true data score. Substituting this result into the reduced Fisher divergence gives

Since multiplying an objective by a positive constant does not change its minimiser, it is common to write the one-dimensional score matching objective as

We now extend the same argument to the general case \(x\in\mathbb{R}^d\). The true data score is

and the cross term in the Fisher divergence becomes

The role played by the ordinary derivative in one dimension is now played by the divergence. For a vector field

the divergence is

As introduced in the previous chapter, divergence measures the local net outflow of a vector field. Here it appears because the model score \(s_\theta(x)\) is itself a vector field.

The multidimensional analogue of integration by parts follows from the divergence theorem. The relevant product rule is

Rearranging gives

Integrating over \(\mathbb{R}^d\), and assuming suitable boundary conditions so that the divergence term contributes no boundary flux, we obtain

Therefore,

This is the direct multidimensional counterpart of the one-dimensional result

The only structural change is that the derivative of the scalar score is replaced by the divergence of the vector-valued score field.

Substituting the cross-term identity into the reduced Fisher divergence gives

Equivalently, after multiplying by \(1/2\), the general score matching objective is

This is the classical score matching objective. It was derived from the Fisher divergence between the learned score and the true data score, but the final expression no longer contains \(s_{\text{data}}(x)\) explicitly. Once \(s_\theta(x)\) approximates \(s_{\text{data}}(x)\), it can be inserted into the Langevin dynamics introduced earlier. Under suitable conditions, this learned dynamics can then be used to generate samples whose distribution approximates the data distribution.

We begin by perturbing the data with additive Gaussian noise. Given a clean sample

we construct a noisy observation

Equivalently, the conditional distribution of the corrupted sample given the clean sample is

The marginal distribution of \(\tilde{x}\) is obtained by integrating out the clean variable \(x\):

This is a convolution of the original data distribution with a Gaussian kernel. The effect is smoothing: probability mass is spread into nearby regions, sharp structures are softened, and areas that previously had negligible density may receive nonzero probability under the corrupted distribution.

Rather than learning the score of the original data distribution, we now aim to learn the score of this smoothed density:

We next examine the structure of this noisy score. By the logarithmic gradient identity,

The gradient of the smoothed density is

where we differentiate under the integral sign under standard regularity conditions. Therefore,

This expression has a useful probabilistic interpretation. By Bayes’ rule, the posterior distribution of the clean variable given the corrupted observation is

Using this posterior, the noisy score can be written as

or equivalently,

For Gaussian corruption,

Its logarithm is, up to constants independent of \(\tilde{x}\),

Hence

Substituting this into the conditional expectation gives

This identity gives an intuitive interpretation of the smoothed score. At a corrupted point \(\tilde{x}\), the score points toward the conditional average clean point that could have produced it. In other words, the noisy score gives a local denoising direction: it tells us how \(\tilde{x}\) should be shifted, on average, in order to move back toward the data distribution.

The identity above is conceptually very useful, but it does not yet give us a directly computable training target. The forward corruption process is simple: we can sample

Thus we can generate pairs \((x,\tilde{x})\) without difficulty.

The problem is that the noisy score involves an expectation under the posterior distribution \(p(x\mid\tilde{x})\):

By Bayes’ rule,

where

Although the likelihood \(p(\tilde{x}\mid x)\) is known, the data distribution \(p_{\text{data}}(x)\) is not available in closed form. We only have samples from it. Therefore the marginal density \(p_\sigma(\tilde{x})\) and the posterior \(p(x\mid\tilde{x})\) are not directly available.

This prevents us from directly minimising the noisy Fisher divergence

because evaluating this objective would require explicit access to the true noisy score \(s_\sigma(\tilde{x})\). Thus, although the forward corruption process is tractable, the corresponding score remains inaccessible in closed form.

We therefore face a familiar situation: we have a well-defined score estimation problem, but its target cannot be computed directly. The crucial question is whether the known corruption process can be used to avoid evaluating the posterior.

The key observation is that the corruption process defines a tractable joint distribution over \((x,\tilde{x})\):

We may not know the marginal density \(p_\sigma(\tilde{x})\) or the posterior \(p(x\mid\tilde{x})\) explicitly, but we can sample from this joint distribution exactly: draw a clean data point \(x\), add Gaussian noise, and obtain \(\tilde{x}\).

Denoising score matching exploits this joint sampling procedure through a basic property of least-squares regression.

We now apply this identity to the denoising setting. Let the input variable be the corrupted observation

and let the regression target be the scaled residual

Now consider the objective

By the least-squares regression identity, the optimal function satisfies

But this conditional expectation is exactly the noisy score derived earlier:

Therefore,

This is the central insight of denoising score matching. Instead of attempting to compute the posterior \(p(x\mid\tilde{x})\) or the noisy score \(s_\sigma(\tilde{x})\) explicitly, we generate pairs \((x,\tilde{x})\) from the known forward corruption process and train the model to predict the scaled residual

The optimal predictor under squared error loss is precisely the score of the smoothed density.

Equivalently, since

we have

Thus the denoising score matching objective may also be written as

This form makes the learning target especially concrete: the model is trained to predict the negative noise direction, scaled by \(1/\sigma^2\).

Denoising score matching therefore replaces an intractable score estimation problem with a supervised regression problem defined entirely through a known corruption process. The objective depends only on clean data samples and Gaussian noise. It does not require evaluating \(p_{\text{data}}(x)\), \(p_\sigma(\tilde{x})\), or the posterior \(p(x\mid\tilde{x})\).

It is important to note that, for a fixed \(\sigma\), the model learns the score of the corrupted density \(p_\sigma\), not directly the score of \(p_{\text{data}}\). The noise level therefore controls the distribution whose geometry is being learned. However, learning the score at a single noise level raises a further question: how should the noise scale \(\sigma\) be chosen?

The previous section showed that, for a fixed noise level \(\sigma\), denoising score matching learns the score of a smoothed data distribution \(p_\sigma\). We corrupt data by

and train

The corresponding denoising target is

However, a single noise level gives only one smoothed view of the data distribution. If \(\sigma\) is large, the density is smooth and the score field provides stable global guidance toward the broad data region, but fine details are blurred away. If \(\sigma\) is small, the density is closer to the original data distribution and the score field captures sharper local detail, but it may be difficult to estimate reliably away from the data support. Thus one noise level cannot simultaneously provide robust global guidance and precise local refinement.

For example, suppose the data lie near two thin rings in \(\mathbb{R}^2\). With large Gaussian noise, the two rings are heavily blurred and may appear as one broad cloud. With an intermediate noise level, the two regions begin to separate. With small noise, the thin circular structures become visible again. The corresponding score fields therefore move from global guidance, to mode separation, to local refinement near the data manifold.

This suggests that the noise scale should not be chosen once and for all. Instead, we should learn a family of score fields across a sequence of noise levels.

Let

be a decreasing sequence of noise scales. For each \(\sigma_k\), define the corrupted density

Instead of learning one score field, we now train a noise-conditional model

The denoising score matching objective extends directly from the single-scale case:

Equivalently, writing \(\tilde{x}=x+\varepsilon\) with \(\varepsilon\sim\mathcal{N}(0,\sigma_k^2I)\), this becomes

For each fixed noise level, the same regression argument as before applies. The optimal predictor recovers the score of the corresponding smoothed density:

The difference is that the model now learns a family of score fields rather than a single one. Instead of modelling one density, we therefore model a hierarchy of smoothed densities, ranging from highly noisy and simple to lightly noisy and detailed.

Once scores have been learned at multiple noise levels, they can be used for sampling. This leads to annealed Langevin dynamics, which applies Langevin updates from large noise scales to small noise scales.

Suppose

with

We initialise from a broad noise distribution, for example

and then, for each scale \(\sigma_k\), perform several updates of the form

After several steps at scale \(\sigma_k\), the procedure moves to the next smaller scale \(\sigma_{k+1}\).

The logic is simple. At large noise levels, the score gives coarse guidance toward the broad data region. At smaller noise levels, the score gives sharper local corrections. Sampling therefore proceeds as progressive refinement:

Annealed Langevin dynamics turns multi-scale score learning into an explicit generative procedure. It also suggests the next conceptual step. If the number of noise levels becomes very large and the gap between neighbouring levels becomes very small, this discrete refinement process begins to resemble a continuous-time diffusion process.

Annealed Langevin dynamics is already a generative procedure, but it is not yet a diffusion model in the DDPM sense. It uses a collection of noisy densities \(p_{\sigma_k}\) and applies Langevin refinement at each scale. A diffusion model organises the intermediate noisy distributions as the marginals of a single forward noising process. In DDPMs, this process is a discrete Markov chain; in the continuous-time score-based view, it is described by a stochastic differential equation.

The multi-scale construction uses finitely many noise levels to approximate a smooth progression from data to noise. When the number of levels becomes large and neighbouring levels are close together, it is natural to replace the discrete index \(\sigma_k\) by a continuous time variable \(t\). Instead of thinking of separate smoothed densities

we consider a continuous family of densities

This continuous-time description should be read as the analogue of the discrete noising path introduced in the DDPM chapter: the finite sequence of intermediate noisy variables is replaced by a continuously indexed family of noisy distributions.

At one end, \(p_0\) is the data distribution, or an idealised version of it. At the other end, \(p_T\) is designed to be a simple noise-like distribution, often approximately Gaussian.

The score model then changes from a noise-conditional field

to a time-dependent field

Thus the multi-scale score-learning picture becomes

This is the point at which multi-scale score learning becomes tied to a diffusion model. Learning scores at many noise levels gives us a family of vector fields, but by itself it does not yet specify a diffusion model. A diffusion model is obtained when the noisy densities are organised as the time marginals of a single forward noising process. In other words, there must be a stochastic dynamics whose evolving distribution is \(p_t\). The learned score \(s_\theta(x,t)\) is then used not merely as an isolated denoising field, but as the correction term needed to run that noising process backward.

A continuous noising process can be described by a stochastic differential equation of the form

The term \(f(X_t,t)\,dt\) is the deterministic drift, while \(g(t)\,dW_t\) injects Brownian noise. The function \(g(t)\) controls the strength of the noise at time \(t\).

This is exactly the drift-plus-noise structure developed in the previous chapter. If \(p_t(x)\) denotes the density of \(X_t\), then its evolution is governed by the Fokker–Planck equation

for the scalar, isotropic diffusion coefficient \(g(t)\). The drift term transports probability mass through space, while the diffusion term spreads and smooths the density.

In generative modelling, the forward process is not used to generate new data directly. Instead, its role is to define a controlled path from data to noise:

The generative problem is then to run this process backward: start from noise and recover a structured sample.

The forward diffusion process is easy to simulate. Starting from a data sample, we add noise according to the chosen stochastic differential equation and obtain increasingly noisy states. But generation requires the opposite direction. We want to start from a simple noise sample and move backward toward the data distribution.

At first sight, reversing a noisy process may seem impossible, because diffusion appears to destroy information. However, under suitable regularity conditions, the reverse-time evolution of a diffusion process can also be described by a stochastic differential equation. The crucial fact is that the reverse drift depends on the score of the intermediate density \(p_t\).

For the forward process

in the scalar-diffusion setting considered here, the reverse-time dynamics has the form

Here the equation is written in the conventional reverse-time notation: the process is simulated from \(t=T\) down to \(t=0\), so the time increment is understood in the reverse direction. The Brownian motion \(\bar{W}_t\) denotes the noise associated with this reversed process.

This equation is the point where the learned score becomes operational for generation. The forward SDE defines the family of noisy densities \(p_t\), but the reverse SDE tells us how to move from \(p_T\) back toward \(p_0\). Since the true score \(\nabla_x\log p_t(x)\) is unknown, it is replaced in practice by the learned score model \(s_\theta(x,t)\). Thus the learned score field becomes the driving correction used during sample generation.

The most important part of the reverse-time expression is the score-dependent correction

This term contains the local information needed to reverse the effect of diffusion. The forward process spreads probability mass by injecting noise. The reverse process must guide samples back toward regions that are likely under the intermediate density \(p_t\). The score \( \nabla_x\log p_t(x)\) provides exactly this local direction.

This is why continuous-time diffusion models require a time-dependent score model:

The model does not learn only the score of the data distribution. It learns the scores of all intermediate noisy distributions along the path from data to noise.

Once these scores are known or well approximated, generation can be carried out by simulating the reverse-time dynamics. One begins with

where \(p_T\) is a simple noise distribution, and then evolves backward from \(T\) to \(0\) using the learned score \(s_\theta(x,t)\). The final state is then treated as an approximate sample from the data distribution.

This completes the conceptual transition from score learning to score-based diffusion. Denoising score matching explains how to learn scores of noisy distributions. Multi-scale score learning extends this across many noise levels. The forward diffusion process organises these noisy distributions into a single time-indexed path, and the reverse-time SDE shows how learned scores can be used to travel back along that path.

Thus, a diffusion model is not simply a score estimator, and score matching alone is not yet a complete generative model. Rather, the forward diffusion process defines a structured path from data to noise, while score learning provides the local information needed to reverse that path. Generation is obtained by combining the two: start from a simple noise distribution and follow a learned, score-guided reverse dynamics back toward the data distribution.

Recall that the DDPM forward process gradually corrupts a clean data sample \(x_0\) by adding Gaussian noise. A key property of this process is that the noisy variable \(x_t\) can be sampled directly from \(x_0\) as

where \(\bar{\alpha}_t\) controls the remaining signal strength at time \(t\). Equivalently,

From the score-based viewpoint, the relevant object is the gradient of a log-density with respect to the noisy variable \(x_t\). For the conditional Gaussian \(q(x_t\mid x_0)\), we have

Differentiating with respect to \(x_t\) gives

Using the reparameterisation

this becomes

This identity explains why noise prediction and score estimation are closely connected. If a model predicts the added noise \(\varepsilon\), then, up to the known scaling factor

it also predicts the conditional score of the noisy Gaussian distribution \(q(x_t\mid x_0)\).

It is important to be precise about what this means. The expression above is the conditional score, where the clean sample \(x_0\) is treated as known. The score needed by the reverse process is the marginal score

where \(q_t(x_t)\) is the full distribution of noisy samples at time \(t\). As in denoising score matching, the marginal score is obtained by averaging over the possible clean samples \(x_0\) that could have produced the same noisy observation \(x_t\).

Thus DDPM noise prediction provides a practical regression route to the same kind of information required by score-based models: the local direction that tells a noisy sample how it should be corrected in order to become more data-like.

We can now place the different viewpoints side by side.

From the DDPM latent-variable perspective, a diffusion model is built from a discrete forward noising process

together with a learned reverse process

Training can be expressed through a variational objective, and in the standard parameterisation this leads to a neural network that predicts the noise added at each timestep.

From the score-based perspective, the same sequence of noisy variables defines a family of intermediate noisy distributions

At each timestep, the reverse process needs local information about how to move a noisy sample toward regions of higher probability under the corresponding intermediate distribution. This local information is encoded by the score

The two views are therefore complementary. The DDPM view emphasises the probabilistic latent-variable structure and the discrete reverse Markov chain. The score-based view emphasises the geometry of the intermediate noisy densities and the local denoising directions needed to reverse the corruption.

In continuous time, this geometric picture becomes especially clear. The discrete family

is replaced by a continuous path of densities

and the reverse dynamics explicitly contains the score

Thus the score-based continuous-time formulation can be viewed as a smooth analogue of the same denoising principle that appears in DDPMs.

The main conceptual message is not that all formulations are identical in their details. DDPMs, annealed Langevin dynamics, and continuous-time score-based diffusion organise the generative process in different ways. But they share the same central idea:

and use those learned corrections to transform noise into data.

In this sense, diffusion models are structured forms of time-dependent denoising. The forward process creates a hierarchy of noisy distributions. The learning problem estimates the information needed to move backward through that hierarchy. The sampling procedure starts from simple noise and repeatedly applies learned corrections until a data-like sample is obtained.

This completes the connection between the discrete DDPM viewpoint and the score-based viewpoint developed in this chapter. The DDPM chapter showed how diffusion models can be derived as latent-variable generative models. The present chapter shows why their denoising behaviour can also be understood geometrically: they learn score-like information across noise levels and use it to guide the reverse transformation from noise to data.

Although the change-of-variables formula was introduced in Chapter 5 from the more general viewpoint of deterministic density transport, it is worth revisiting it here because it is the modelling principle of normalising flows.

Suppose \(X\) has a known density \(p_X(x)\), and we apply a differentiable, invertible function \( Y=f(X).\) Our goal is to determine the density of \(Y\), denoted \(p_Y(y)\). Intuitively, the density of \(Y\) needs to reflect two things.

First, it must account for where probability mass came from. Since \(y=f(x)\), the likelihood of observing \(y\) depends on how likely its corresponding pre-image \( x=f^{-1}(y)\) was under the original distribution.

Second, it must account for how the transformation stretches or compresses space. A transformation that spreads points further apart should reduce the density, while one that compresses points should increase it. Density is always probability per unit length, area, or volume, so local spacing matters.

In one dimension, these ideas combine in the change-of-variables formula

The first term tells us which value of \(x\) generated \(y\). The second term is a scaling factor that accounts for how stretched or compressed the neighbourhood around \(x\) becomes when mapped to \(y\).

8.3.1.1 A simple geometric example.

8.3.1.2 The higher-dimensional case: from stretching lengths to stretching volumes.

8.3.1.3 The Jacobian as a local linear lens.

The examples above were intentionally simple. They showed how the change-of-variables formula tracks density under linear stretching and volume scaling. However, linear transformations have limited expressive power. They can stretch, rotate, or shear space, but cannot bend or warp it in rich nonlinear ways. To approximate complex distributions, we need transformations that introduce carefully controlled nonlinearity while preserving the two core requirements of normalising flows: invertibility and a tractable Jacobian determinant.

Planar flows provide a classic example of this trade-off, originally introduced as part of the normalising-flow framework of Rezende and Mohamed [27]. A foundational nonlinear flow layer is defined as

where \(x\in\mathbb{R}^d\), \(u,w\in\mathbb{R}^d\), \(b\in\mathbb{R}\), and

is a smooth activation function such as \(\tanh\) or sigmoid. Here \(x\) denotes the input to this individual flow layer; in a stacked flow it would usually be one of the intermediate variables.

The term

acts as an affine projection of the input, reducing the high-dimensional vector \(x\) to a single scalar. On its own, this projection would contribute only a linear quantity. To introduce curvature and richer behaviour, we pass this scalar through a nonlinear activation \(h(\cdot)\). Functions such as \(\tanh\) and sigmoid squash their inputs into bounded ranges, so the output of \(h\) can be interpreted as a soft gating value that controls how strongly the transformation should modify the input.

Multiplying this gated scalar by the vector \(u\) produces a directionally scaled perturbation. The deformation always occurs along the direction of \(u\), but its magnitude depends nonlinearly on the projection \(w^\top x+b\). Adding this perturbation back to \(x\) ensures that the output remains in the full \(d\)-dimensional space while introducing a flexible nonlinear adjustment.

8.3.2.1 Jacobian derivation.

8.3.2.2 Invertibility condition.

Affine coupling layers, introduced in coupling-based flow models such as NICE and RealNVP [28, 29], provide one of the most practical ways to build expressive flows with tractable inverses and log-determinants.

An affine coupling layer begins by splitting the input vector into two parts. Let \(x\in\mathbb{R}^d\) be the input to a flow layer, and write

where \(x_1\in\mathbb{R}^{d_1}\), \(x_2\in\mathbb{R}^{d_2}\), and \(d_1+d_2=d\). The affine coupling transformation is defined by

where \(s(\cdot)\) and \(t(\cdot)\) are functions, usually neural networks, and \(\odot\) denotes elementwise multiplication.

The first part of the input is left unchanged. The second part is transformed by an elementwise affine map whose scale and shift are predicted from the first part. Thus \(x_1\) acts as a context variable, while \(x_2\) is scaled and translated according to that context.

Although the transformation is affine in \(x_2\), the functions \(s(x_1)\) and \(t(x_1)\) may be highly nonlinear. This allows the overall transformation to be expressive while preserving a simple algebraic structure. The layer is therefore deliberately asymmetric: it sacrifices full simultaneous mixing inside one layer in exchange for exact inversion and an efficient determinant.

8.4.1.1 Jacobian determinant.

8.4.1.2 Invertibility.

A single flow layer can only produce a limited deformation of space. Planar flows perturb the identity in one learned direction, while affine coupling layers transform only part of the input at a time. The expressive power of normalising flows therefore comes from composition: many simple, invertible transformations are stacked so that their small deformations accumulate into a rich global mapping.

Let

be a sequence of differentiable and invertible transformations. Starting from a base variable

we define intermediate variables

and the final output

Equivalently,

Because each layer is invertible, the full composition is also invertible. This means we can move in both directions:

or

This two-way structure gives the two basic operations of a flow model: likelihood evaluation and sampling.

8.4.2.1 Likelihood evaluation: data to base space.

8.4.2.2 Sampling: base space to data space.

8.4.2.3 Composition, mixing, and expressiveness.

Let

be a \(d\)-dimensional random vector. The chain rule of probability follows from repeated application of the identity

Starting from the full joint density, we can first separate out the final variable:

Applying the same identity recursively to the remaining joint term gives

Reordering the factors yields the standard autoregressive factorisation:

where

and for \(i=1\), the term \(p(x_1\mid x_{<1})\) is simply interpreted as \(p(x_1)\).

This identity is exact. It makes no approximation and imposes no modelling assumption by itself. Any joint distribution can be decomposed in this way once an ordering of the variables has been chosen.

The modelling step comes from parameterising the conditional densities:

If each conditional density is normalised, then the product defines a properly normalised joint density. This is why autoregressive models give exact likelihoods by construction.

The chain-rule factorisation becomes a learnable generative model once we specify how each conditional density is represented. In general, we write

where \(\theta\) denotes model parameters.

An important point is that these conditionals are usually not represented by \(d\) unrelated models. Instead, a single parameterised function, often a neural network, is used to output the parameters of all conditional distributions while respecting the autoregressive ordering. The same parameters are reused across the factorisation; what changes from one conditional to the next is the conditioning context.

For continuous variables, a common simple choice is a Gaussian conditional:

where the mean and variance are functions of the previous variables. More flexible conditionals may use mixtures of Gaussians, discretised logistic mixtures, or other expressive one-dimensional distributions.

For discrete variables, the model may output probabilities for a categorical distribution:

In both cases, the model predicts a full conditional distribution, not merely a point estimate.

The same model parameters \(\theta\) are reused throughout. What changes from one conditional distribution to the next is not the model itself, but the amount of conditioning information available at that step. This resolves a common misconception: the autoregressive factorisation lists a sequence of conditional distributions, but this does not mean that we train a separate unrelated model for each factor. In modern autoregressive models, the conditionals are usually represented together by one shared architecture with restricted information flow.

The same factorisation also gives a direct sampling procedure. Given

we may generate a sample sequentially:

then

then

and so on until

After all components have been generated, the resulting vector

is a sample from the autoregressive model.

This procedure should be read as a full generative mechanism. The model does not first sample a latent variable and then decode it. Nor does it run a stochastic denoising process. It directly constructs the data point one component at a time, using previously generated components as context for later ones.

Autoregressive models have an important computational asymmetry. Likelihood evaluation during training can often be parallelised, while sampling is inherently sequential.

Given an observed data point

the log-likelihood is

During training, all components of \(x\) are known. Therefore, every conditioning context \(x_{<i}\) can be read directly from the observed data. No sampling is required. The model evaluates each conditional term using the ground-truth previous components.

This is why autoregressive likelihood training is exact and stable. The model is not asked to generate \(x_1\), then use the generated value to evaluate \(x_2\), and so on. Instead, it evaluates all conditional probabilities against the observed data:

Depending on the architecture, many or all of these conditional terms can be computed in parallel.

Sampling is different. When generating a new sample, future components are unknown because they have not yet been produced. We must therefore proceed in order:

Each newly sampled component becomes part of the conditioning context for the next one. This makes sampling sequential by construction.

To implement the autoregressive factorisation in a neural network, we must also ensure that the parameters of the \(i\)-th conditional distribution depend only on \(x_{<i}\), and not on \(x_i\) itself or on future variables. This restriction is commonly enforced through masking. A mask blocks forbidden connections inside the model, ensuring that information flows only from earlier variables to later variables.

Thus, even if the network processes the whole vector \(x\) at once during training, the output parameters for \(x_i\) cannot depend on \(x_i,x_{i+1},\dots,x_d\). The computation may be parallel, but the dependency structure remains ordered.

Abstractly, masking imposes a lower-triangular dependency pattern. If we define a dependency matrix \(D\) whose entry \(D_{ij}\) indicates whether the conditional for \(x_i\) is allowed to depend on \(x_j\), then autoregressive ordering requires

For example,

where \(\ast\) denotes a permitted dependency.

This triangular dependency is a key reason the model remains tractable in practice. It avoids cyclic relationships: the distribution of \(x_i\) may depend on earlier components, but earlier components do not simultaneously depend on \(x_i\). As a result, the high-dimensional joint likelihood decomposes into a sum of conditional log-likelihoods.

In this sense, triangular dependency in autoregressive models plays a role analogous to triangular Jacobians in coupling-based normalising flows. In both cases, dependencies are organised in one direction so that a high-dimensional probability calculation factorises into simple terms.

With this principle in place, many well-known architectures can be understood at a structural level. MADE uses masking to enforce autoregressive structure in feedforward networks [30]; PixelRNN and PixelCNN apply autoregressive factorisation to images [31]; and WaveNet uses causal convolutions for raw audio generation [32].

Autoregressive models therefore provide a second major route to exact density modelling. Normalising flows achieve exact likelihoods through invertible transformations and tractable Jacobian determinants. Autoregressive models achieve exact likelihoods through ordered conditional factorisation. Both show how high-dimensional probability models can remain computable when their dependencies are carefully structured.

A GAN begins with a generator, a parameterised function

where \(z\) is a latent variable drawn from a simple fixed distribution, such as a standard Gaussian or a uniform distribution. The generator \(G_\theta\), often implemented as a neural network, maps this simple source of randomness into the data space.

The mapping

transforms simple noise into structured, data-like samples. By sampling \(z\sim p(z)\) and passing it through \(G_\theta\), the generator induces a distribution over \(x\), denoted \(p_\theta(x)\). Crucially, this distribution is defined only implicitly: we can draw samples from it, but we cannot usually evaluate its density at a given point.

This distinguishes GANs from the likelihood-based models studied earlier. Unlike VAEs, the latent variable \(z\) is not inferred from data and no posterior distribution is learned. Unlike diffusion models, \(z\) does not arise from a forward noising process. Unlike normalising flows or autoregressive models, the generator need not be invertible, and no Jacobian determinant or explicit normalisation of \(p_\theta(x)\) is required. The generator therefore behaves like a deterministic simulator: it specifies how to generate data, but not how likely a given data point is.

Once training is complete, sampling from a GAN is straightforward. We draw \( z\sim p(z)\) from the latent noise distribution and compute \( x=G_\theta(z).\) The discriminator is not needed at generation time; it serves only as a training signal that shapes the generator. Thus GANs provide direct sampling, but not direct likelihood evaluation.

If a model produces only samples, a fundamental question immediately arises:

In likelihood-based models, this question is answered by probability: samples are good if they receive high likelihood under the model. In GANs, likelihoods are unavailable. There is no tractable function \(p_\theta(x)\) to evaluate or maximise.

To resolve this, GANs introduce a second learned component: the discriminator. Rather than measuring probability, the discriminator provides a data-driven comparison between two sources of samples:

The first source consists of real data samples; the second consists of generated samples.

Formally, the discriminator is a function

parameterised by \(\phi\). For an input \(x\), the value \(D_\phi(x)\) is interpreted as the discriminator’s confidence that \(x\) came from the real data distribution rather than from the generator. Thus the discriminator is trained as a binary classifier:

The discriminator is not an end in itself. It does not define a probability density, nor is it the final object we want to use after training. Its role is to provide a moving, adaptive notion of discrepancy between real and generated samples. The generator learns only through this feedback, adjusting its parameters so that generated samples become harder to distinguish from real ones.

This interaction is formalised by the minimax objective

The discriminator tries to maximise this objective by assigning high scores to real samples and low scores to generated samples. The generator tries to minimise it by producing samples for which \(D_\phi(G_\theta(z))\) becomes large, so that the discriminator struggles to classify them as fake.

Although the objective resembles a classification loss, its role is fundamentally generative. The discriminator defines a learned comparison between two distributions, while the generator reshapes its output distribution so that this comparison gradually disappears.

To understand adversarial training at the distributional level, consider an idealised setting in which the generator \(G_\theta\) is fixed. Sampling \(z\sim p(z)\) and mapping \(x=G_\theta(z)\) induces an implicit generator distribution \(p_\theta(x)\). With the generator fixed, the discriminator’s task reduces to binary classification between

The discriminator objective can then be written as

Writing these expectations as integrals gives

Because this objective decomposes over \(x\), the optimal discriminator can be derived pointwise. For each fixed \(x\), we maximise

with respect to the scalar \(D(x)\).

Differentiating with respect to \(D(x)\) gives

Setting this derivative to zero yields

Solving for \(D(x)\), we obtain the optimal discriminator

This formula gives an important interpretation. At optimality, the discriminator does not merely label samples as real or fake. It estimates the relative amount of probability mass assigned by the two distributions at \(x\). Where \(p_{\text{data}}(x)\) dominates, \(D^*(x)\) approaches one. Where \(p_\theta(x)\) dominates, \(D^*(x)\) approaches zero. Thus the discriminator acts as a local density-ratio estimator, even though neither density is explicitly available to the model.

Having derived the optimal discriminator for a fixed generator, we can ask what the adversarial game represents at the level of entire distributions. To do this, we write the discriminator objective as a value function, which scores a generator–discriminator pair by how well the discriminator separates real samples from generated ones.

We define

For fixed \(G_\theta\), maximising \(V(G_\theta,D)\) over \(D\) gives the best discriminator against that generator. If we then evaluate the value function at this best discriminator, we can identify the distributional discrepancy that the generator is implicitly minimising.

Using

after algebraic simplification, we obtain

where \(\mathrm{JS}(\cdot\|\cdot)\) denotes the Jensen–Shannon divergence.

The constant term \(-\log 4\) does not depend on the generator, and therefore does not affect the optimisation. Thus, in the idealised setting where the discriminator is optimal, GAN training corresponds to minimising the Jensen–Shannon divergence between the data distribution and the generator distribution.

The expression above introduces a new divergence between probability distributions: the Jensen–Shannon divergence. To understand it, it is helpful to recall the Kullback–Leibler divergence, which appeared earlier in the VAE chapter as part of the evidence lower bound. The KL divergence compares two distributions by

It measures how costly it is to use \(q\) as an approximation to \(p\). However, KL divergence is asymmetric:

and it is sensitive to support mismatch. If \(q(x)=0\) in a region where \(p(x)>0\), then \(\mathrm{KL}(p\,\|\,q)\) becomes infinite.

The Jensen–Shannon divergence modifies this comparison by comparing both distributions to a shared intermediate distribution. Given two distributions \(p\) and \(q\), define their mixture

This mixture can be understood as the distribution obtained by first choosing between \(p\) and \(q\) with equal probability, and then sampling from the chosen distribution.

The Jensen–Shannon divergence is then defined as

Thus, rather than asking one distribution to explain the other directly, JS divergence asks how far each distribution is from their shared mixture. This makes it symmetric:

and bounded:

This is precisely the divergence that appears when the optimal discriminator is substituted back into the original GAN value function. In the idealised setting where the discriminator is optimal, the generator is therefore not maximising likelihood. Instead, it is adjusting its implicit distribution \(p_\theta\) so as to reduce

This reveals the distributional meaning of adversarial training. GANs do not construct an explicit density; they perform distribution matching through comparison. The discriminator estimates how distinguishable the real and generated distributions are, and the generator learns by trying to reduce that distinguishability.

The global optimum is achieved when

In that case,

everywhere, because real and generated samples are indistinguishable.

Before moving to Wasserstein GANs, we first discuss a common practical workaround within the original GAN framework: changing the generator’s loss to improve the gradient signal during training. The original GAN objective is written as the minimax problem

For the discriminator, the interpretation is straightforward: it is rewarded for assigning high values to real samples and low values to generated samples. For the generator, the situation is more subtle.

If we fix the discriminator and minimise the same value function with respect to the generator parameters \(\theta\), the only term that depends on \(\theta\) is

This objective asks the generator to reduce the discriminator’s confidence that generated samples are fake.

The difficulty appears early in training. At this stage, the discriminator can often identify generated samples confidently:

Then

The generator objective is therefore close to flat in this region, and the gradient passed back to the generator can become very small even though its samples are still poor.

For this reason, GANs are commonly trained with a different generator objective, called the non-saturating loss:

This loss is not obtained by algebraically simplifying the previous expression. Instead, it changes the generator’s goal from reducing the discriminator’s “fake” confidence to increasing its “real” confidence.

The two objectives are different, but they have the same desired equilibrium: generated samples should become indistinguishable from real data. Their difference lies in the gradient behaviour. When

the non-saturating term

is large, so the generator receives a much stronger learning signal. This is why the non-saturating loss is commonly preferred in practice, even though the original minimax objective provides the cleaner theoretical link to Jensen–Shannon divergence. This optimisation behaviour helps explain both the strength and the weakness of GANs. On the positive side, the discriminator penalises features that make generated samples look unrealistic, so adversarial training can encourage sharp, visually convincing outputs rather than blurry averages. On the negative side, the mechanism can give rise to mode collapse. Since the generator is trained to fool the discriminator, it may model only a subset of the data distribution convincingly. If that subset is enough to confuse the discriminator, other modes may be ignored.

Likelihood-based models are penalised when they fail to cover observed modes, because likelihood depends on assigning probability mass to the data. GANs, by contrast, optimise distinguishability rather than explicit probability coverage. Thus the non-saturating loss improves the generator’s gradients, but it does not change the underlying geometry of the adversarial comparison. When two distributions have little overlap, a classification-based discrepancy may still say mainly that they are distinguishable, not how one should move toward the other. This motivates a more geometric notion of distributional discrepancy.

To develop this idea, it is helpful to think of probability distributions geometrically. Rather than viewing a distribution as a static density, imagine its probability mass as something that can be moved through space. If transforming one distribution into another requires only small displacements of mass, the distributions should be considered close. If large amounts of mass must be transported over long distances, they should be considered far apart.

The Wasserstein distance formalises this intuition. In the form used by Wasserstein GANs, the relevant object is the \(1\)-Wasserstein distance, also called the Earth Mover’s distance. It measures the minimal cost of transporting probability mass from one distribution to another, with cost proportional to the distance moved.

A key theoretical result, the Kantorovich–Rubinstein duality, gives a practical way to express this distance. For two probability distributions \(p\) and \(q\) on \(\mathbb{R}^d\), the \(1\)-Wasserstein distance can be written as

Here, \(f\) is a real-valued scoring function, and the constraint

means that \(f\) must be \(1\)-Lipschitz: its value cannot change faster than unit rate with respect to distance in the input space.

The expression

measures how strongly the scoring function separates the two distributions on average. The supremum means that we choose the scoring function that separates them as much as possible, subject to the Lipschitz constraint.

This constraint is essential. Without it, the scoring function could assign arbitrarily large values to samples from one distribution and arbitrarily small values to samples from the other, making the comparison meaningless. The Lipschitz constraint forces the score difference to reflect the geometry of the data space: large score differences are only allowed when the corresponding points are far apart.

Note this is different from the score in score-based diffusion models, where the score means \(\nabla_x\log p(x)\), the gradient of a log-density; in WGANs, the critic output \(f_\phi(x)\) is a real-valued scalar used for comparison, rather than a probability; the critic is constrained through its Lipschitz behaviour, not through probabilistic normalisation.

This example shows why the Wasserstein distance remains informative when distributions do not overlap. An overlap-based divergence can detect that the supports are disjoint, but may not reflect whether the distributions are one unit apart or ten units apart. The Wasserstein distance, by contrast, directly measures spatial separation.

The dual characterisation of the Wasserstein distance leads directly to a practical adversarial objective. If the generator induces the distribution \(p_\theta\), then the Wasserstein distance between the data distribution and the generator distribution can be written as

Since samples from \(p_\theta\) are produced by

we can write

Wasserstein GANs represent the scoring function \(f\) by a neural network \(f_\phi\), and optimise

where \(\mathcal{L}\) denotes the set of parameter values for which \(f_\phi\) satisfies the required Lipschitz constraint.

The maximisation over \(\phi\) seeks a scoring function that separates real and generated samples as strongly as possible while respecting the geometry of the data space. The minimisation over \(\theta\) then updates the generator so that this separation becomes smaller.

At this point, the terminology changes. In the original GAN, the sample-evaluating network is called a discriminator. Its output

is interpreted as the probability that \(x\) is real. In a Wasserstein GAN, the corresponding network is called a critic. It outputs a real-valued score

not a probability. There is no sigmoid output, no binary classification interpretation, and no requirement that the scores be normalised. Only differences in critic values matter.

Seen in this way, the critic is a neural-network approximation to the optimal scoring function in the Wasserstein distance. Its purpose is not to make a hard real–fake decision, but to assign scores whose average difference reflects the geometric separation between the two distributions. Real samples should receive higher scores on average than generated samples, subject to the constraint that scores cannot change arbitrarily fast across the data space. Since the Wasserstein distance changes continuously as the generator distribution moves, the generator can still receive useful gradients even when \(p_\theta\) and \(p_{\text{data}}\) do not overlap, so optimisation is guided by geometric proximity rather than overlap alone. As in the original GAN, the critic is used only during training. Once the generator has been learned, sampling is again performed directly by drawing \(z\sim p(z)\) and returning \(G_\theta(z)\).

The remaining challenge is practical: how can the Lipschitz constraint on the critic be enforced during training? Recall that this constraint ensures the critic’s scores change at a controlled rate, so that score differences reflect geometric distance rather than arbitrary scaling. Neural networks do not satisfy this property automatically.

The original WGAN enforced the constraint through weight clipping, which restricts each parameter of the critic to a fixed interval. Although simple, this is a crude way to control the Lipschitz constant. It can reduce the capacity of the critic and may lead to unstable optimisation.

A more effective alternative is the gradient penalty, introduced in WGAN-GP [37]. Instead of directly constraining the parameters, the gradient penalty constrains the critic’s input–output behaviour. The idea is to control how rapidly the critic’s score changes in response to small changes in its input.

A penalty term of the form

is added to the critic’s training objective, where \(\hat{x}\) is sampled along straight-line interpolations between real samples and generated samples.

The gradient

measures how sensitive the critic’s score is to infinitesimal changes in the input. Its Euclidean norm is therefore the local rate of change, or local slope, of the critic at \(\hat{x}\). Penalising deviations of \( \|\nabla_{\hat{x}} f_\phi(\hat{x})\|_2\) from one encourages the critic to change smoothly and at a controlled rate in the regions most relevant for comparing real and generated samples.

The interpolation points \(\hat{x}\) lie between real and generated samples, precisely where the critic must learn how to compare the two distributions. Enforcing the constraint along these paths stabilises training while preserving the geometric interpretation of the Wasserstein distance. The coefficient \(\lambda\) controls the strength of the penalty: larger values enforce the constraint more strongly, while smaller values allow greater flexibility.

By enforcing the Lipschitz constraint in this way, WGANs with gradient penalty preserve the transport-based meaning of the Wasserstein distance while achieving more reliable training dynamics. More broadly, they show that generative learning need not rely on explicit likelihoods or normalised densities. It can instead be guided by a learned real-valued function whose geometry pushes generated samples toward the data distribution.

This provides a natural bridge to energy-based models. Like WGANs, energy-based models use real-valued functions to compare, rank, or guide samples. But instead of framing this through an adversarial game, they define distributions through unnormalised energy landscapes. We turn to this perspective next.

The key modelling primitive in an energy-based model is the scalar energy function

Unlike normalising flows or autoregressive models, this function does not need to provide an exact density. Unlike VAEs, it does not require an explicit latent posterior. Unlike GANs, it does not require a separate discriminator trained in an adversarial game. The model is specified directly through a scalar landscape.

Only relative energy differences matter. If we add a constant \(c\) to the energy function,

then every configuration is shifted by the same amount. The ordering of points from low to high energy is unchanged. Since the model’s behaviour depends on which configurations are preferred relative to others, this constant shift has no meaningful effect.

This relative nature is important. An energy function is not trying to assign an absolute probability value to every point. Instead, it encodes preferences: which configurations should be encouraged, and which should be discouraged. In this sense, energy-based modelling asks a different question from explicit density modelling:

The answer is learned from data.

Although an energy function is not itself a probability distribution, it can be connected to one through the Boltzmann form

where

is the normalising constant, or partition function.

This construction has a simple interpretation. The exponential transformation turns energies into positive weights. Because of the negative sign, low-energy configurations receive large weights, while high-energy configurations receive small weights. The partition function then normalises these weights so that they integrate to one.

Thus, if the partition function were available, the energy function would define a probability density:

The difficulty is that \(Z_\theta\) is usually intractable in high-dimensional spaces. Computing it requires integrating over all possible configurations:

For images, audio, text representations, or other high-dimensional data, this integral is generally impossible to evaluate exactly.

This is the central trade-off of energy-based modelling. By avoiding the structural constraints required by exact density models, energy functions can be very flexible. But this flexibility comes at the cost of normalisation: the model usually defines an unnormalised density

rather than a tractable normalised likelihood.

This distinguishes energy-based models from the exact density models studied in the previous chapter. Normalising flows and autoregressive models are carefully designed so that \(p_\theta(x)\) can be evaluated exactly. Energy-based models instead place the emphasis on the energy landscape itself. Relative energy differences and energy gradients can be useful even when the normalising constant is unknown.

Learning in an energy-based model means shaping the energy function so that it reflects the structure of the data. Suppose we observe samples

A desirable energy function should assign low energy to these observed samples:

But this condition alone is not enough. If we simply lowered the energy everywhere, the model would not distinguish data-like configurations from implausible ones. Learning must therefore be contrastive: it must lower the energy on data samples while raising the energy elsewhere.

Informally, the desired behaviour is

This contrastive structure is what makes energy-based learning meaningful. The model must not merely reward observed data; it must also learn which alternative configurations should be discouraged.

A useful way to see this is through the energy gradient

This gradient describes how the energy changes locally in data space. Moving in the direction

decreases energy, guiding a configuration toward regions that the model considers more plausible. Thus the energy landscape defines not only a scalar preference value, but also a local direction of improvement.

This connects energy-based models to the geometric viewpoint developed earlier for GANs and score-based models. In a Wasserstein GAN, the critic assigns real-valued scores whose differences guide the generator. In an energy-based model, the energy function assigns scalar values whose gradients define how configurations should move through the data space. The two ideas are related, but not identical: a WGAN critic is trained to separate two distributions, whereas an energy function is intended to define a landscape of plausibility.

Although energy-based models are often discussed as unnormalised models, the Boltzmann form still provides a principled learning objective. If

then maximum likelihood seeks to maximise

Substituting the energy-based density gives

and therefore

The first term has an intuitive effect: maximising it lowers the energy of data samples. The second term prevents the trivial solution of lowering energy everywhere, because the partition function increases when too much of the space is assigned low energy. Thus maximum likelihood balances two forces:

This becomes clearer by differentiating the objective with respect to \(\theta\). The gradient of the data term is

For the partition function term, we use

Differentiating gives

Since the log-likelihood contains \(-\log Z_\theta\), its contribution to the gradient is

Combining the two terms, the maximum-likelihood gradient becomes

This expression is central to energy-based learning. It has two parts:

The data term lowers the energy of observed samples. The model term raises the energy of samples drawn from the current model distribution. This is often described as a contrast between a positive phase and a negative phase: the positive phase makes data more preferred, while the negative phase pushes back against configurations that the model itself currently assigns too much probability mass.

This contrastive gradient exposes the central practical difficulty of energy-based models. The expectation over the data distribution is easy to approximate using the training set. The expectation over the model distribution

is much harder, because drawing samples from \(p_\theta(x)\propto \exp(-E_\theta(x))\) is itself a nontrivial problem.

This creates a coupling between learning and sampling. To improve the energy function, we need samples from the current model. But to draw good samples, we need an energy landscape that has already been shaped well. This circular dependence is one of the main reasons energy-based models are powerful but difficult to train.

The next section develops this connection more explicitly. Since an energy-based model does not provide a direct sampler, generation is usually formulated as movement on the energy landscape. This leads naturally to energy-gradient dynamics, Langevin sampling, and the connection between energy-based modelling and the score-based diffusion methods studied earlier.

An energy-based model assigns low energy to plausible configurations, but it does not directly tell us how to draw independent samples from the corresponding distribution. Unlike a normalising flow, there is no invertible map from a base distribution to data space. Unlike an autoregressive model, there is no left-to-right product of conditionals that can be sampled sequentially. Instead, sampling is usually formulated as an iterative dynamical process.

The basic idea is to initialise a configuration \(x_0\) and repeatedly update it so that it moves toward lower-energy regions. A purely deterministic update would take the form

where \(\eta>0\) is a step size. This resembles gradient descent, but now the descent is performed in the data space rather than in the parameter space. The configuration \(x_t\) itself is being moved toward lower energy.

However, deterministic descent alone is not sufficient for sampling. It may collapse into a local minimum of the energy landscape and fail to explore multiple plausible regions. To obtain a sampling procedure, randomness is added:

This is the discrete Langevin-style update associated with the energy function.

The deterministic term

moves the current configuration toward lower energy. The noise term

injects randomness, preventing the dynamics from simply collapsing to a single minimum and allowing it to explore the broader energy landscape. Under suitable conditions and with an appropriate continuous-time interpretation, this gradient-plus-noise dynamics can have the energy-based distribution as its stationary distribution.

This is the same Langevin idea studied earlier in the score-based chapter. The only change is how the drift field is written. There, the drift was expressed using the score \(\nabla_x\log p(x)\). Here, it is expressed using the negative energy gradient \(-\nabla_xE_\theta(x)\). The next subsection makes this connection explicit.

Recall that the score of a density \(p(x)\) is

For an energy-based model,

Taking logarithms gives

Now differentiate with respect to \(x\). Since the partition function \(Z_\theta\) depends on \(\theta\) but not on \(x\), we have

Therefore,

This identity is the bridge between energy-based modelling and score-based modelling. The score field of an energy-based distribution is exactly the negative gradient of its energy function. Thus the local direction used to move toward higher probability is the same as the direction used to move toward lower energy.

The difference lies in how the object is represented. In a classical energy-based model, we parameterise a scalar function \(E_\theta(x)\), and the driving vector field is obtained by differentiating it:

In a score-based model, we often parameterise the vector field directly:

Thus, an energy model starts from a scalar landscape and derives a vector field, whereas a score model often learns the vector field itself.

This distinction is useful. A scalar energy automatically gives a conservative gradient field: the vector field comes from the gradient of one underlying function. A directly learned score network is more flexible as a vector field, but it may not always be explicitly represented as the gradient of a scalar energy unless additional structure is imposed. For the purposes of generative modelling, however, the conceptual link remains strong: both approaches learn local directions that guide samples toward regions of higher data density.

The previous section showed that maximum likelihood for energy-based models leads to a contrastive gradient involving samples from the model distribution. This creates a practical difficulty because sampling from the model is itself hard. Score matching [23] provides a different route: instead of matching normalised probabilities, it matches local score fields.

In the score-based chapter, we derived score matching in detail using the Fisher divergence and integration by parts. Here we only recall the consequence needed for the energy-based view. The goal is to learn a model score

For an energy-based model,

Thus score matching can be interpreted as learning an energy landscape whose negative gradient agrees with the data score.

The classical score-matching objective can be written, up to constants independent of \(\theta\):

Substituting

gives the energy-form objective

where

is the Laplacian of the energy function.

The important point is not the detailed derivation, which was already developed earlier, but the change in what is being matched. Maximum likelihood compares normalised probability distributions and produces a model expectation term. Score matching compares local gradients of log-density and avoids the partition function, because

In energy terms, it learns the shape of the energy landscape near data samples through derivatives, rather than requiring exact evaluation of the normalised density.

This is why score matching is so natural from an energy-based perspective. It does not ask for the probability of a data point. It asks for the local direction in which probability should increase, or equivalently the direction in which energy should decrease.

The relationship between energy-based models and score-based diffusion can now be stated clearly. An energy-based model defines a scalar landscape \(E_\theta(x)\). If this landscape is associated with the unnormalised density

then its score is

Thus the local direction used in Langevin sampling is the same object that score-based models try to learn directly: a direction pointing toward higher probability, or equivalently toward lower energy.

From this perspective, score-based diffusion models [25, 26] can be viewed as learning time-dependent energy-gradient-like directions across a hierarchy of noisy distributions. The difference is therefore not the local geometry itself, but how that geometry is represented and learned. Energy-based models parameterise a scalar energy and derive a vector field by differentiation. Score-based models usually parameterise the vector field directly. Diffusion models go one step further by learning a time-dependent score field

where \(p_t(x)\) is the distribution obtained after corrupting the data with noise up to time \(t\).

From an energy-based viewpoint, each noisy distribution \(p_t\) can be associated with an implicit energy landscape

where \(C_t\) is an arbitrary constant. Since only gradients matter, this constant has no effect, and

Therefore, a diffusion model can be interpreted as learning

the negative gradients of a whole family of energy landscapes indexed by noise level.

This viewpoint is useful because the forward noising process creates a multi-scale hierarchy of distributions. At large noise levels, \(p_t\) is smooth and close to a simple distribution, so the corresponding energy landscape is broad and regular. At small noise levels, \(p_t\) is closer to the data distribution, so the energy landscape contains finer structure. Diffusion training therefore learns the local geometry of this sequence:

In this sense, diffusion models may be viewed as a practical and stable form of multi-scale energy-gradient learning. Instead of learning one sharply structured energy landscape for the data distribution and sampling from it directly, they spread the problem across many smoothed distributions. The high-noise levels provide global guidance, while the low-noise levels provide detailed local refinement.

This also explains why diffusion models avoid some of the traditional difficulties of energy-based learning. Classical EBM maximum likelihood requires samples from the current model distribution in order to estimate the negative phase. Diffusion models instead create noisy training pairs through a known forward corruption process and learn denoising directions across noise levels. The result remains deeply connected to energy gradients, but the training procedure is far more stable.

Energy-based models therefore provide the broader landscape view of this part of generative modelling. They describe distributions through scalar compatibility functions, while score-based diffusion can be seen as one particularly effective modern way to learn and use the corresponding energy-gradient directions across a multi-scale path from noise to data.