NON-DENOISING FORWARD-TIME DIFFUSIONS
Year: 2021

Abstract:
The scope of this paper is generative modeling through diffusion processes. An approach falling within this paradigm is the work of Song et al. (2021), which relies on a time-reversal argument to construct a diffusion process targeting the desired data distribution. We show that the time-reversal argument, common to all denoising diffusion probabilistic modeling proposals, is not necessary. We obtain diffusion processes targeting the desired data distribution by taking appropriate mixtures of diffusion bridges. The resulting transport is exact by construction, allows for greater flexibility in choosing the dynamics of the underlying diffusion, and can be approximated by means of a neural network via novel training objectives. We develop a unifying view of the drift adjustments corresponding to our and to time-reversal approaches and make use of this representation to inspect the inner workings of diffusion-based generative models. Finally, we leverage on scalable simulation and inference techniques common in spatial statistics to move beyond fully factorial distributions in the underlying diffusion dynamics. The methodological advances contained in this work contribute toward establishing a general framework for generative modeling based on diffusion processes.

heading: 1 INTRODUCTION
text: Denoising diffusion probabilistic modeling (DDPM) (Sohl-Dickstein et al., 2015; Ho et al., 2020; Song et al., 2021) is a recent generative modeling paradigm exhibiting strong empirical performance. Consider a dataset of N samples D = {x(n)}Nn=1 with empirical distribution PD. The unifying key steps underlying DDPM approaches are: (i) the definition of a stochastic process with initial distribution PD, whose forward-time (noising) dynamics progressively transform PD toward a simple data-independent distribution PZ ; (ii) the derivation of the backward-time (denoising / sampling) dynamics transforming PZ toward PD; (iii) the approximation of the backward-time transitions by means of a neural network. Following the training step (iii), a sample whose distribution approximates PD is drawn by (iv) simulating from the approximated backward-time transitions starting with a sample from PZ . Both discrete-time (Ho et al., 2020) and continuous-time (Song et al., 2021) formulations of DDPM have been pursued. This work focuses on the latter case, to which we refer as diffusion time-reversal transport (DTRT). As in DTRT, dynamics are specified through diffusion processes, i.e. solutions to stochastic differential equations (SDE) with associated drift f( ·) and diffusion g( ·) coefficients. A number of approximations are involved in the aforementioned steps. Firstly, as the dynamics are defined on a finite time interval, a dependency from PD is retained through the noising process. Hence, starting with a sample from the data-independent distribution PZ in (iv) introduces an approximation. Secondly, while the backward-time dynamics of (ii) are directly available for diffusions, they are approximated by means of a neural network in (iii). Thirdly, sampling in (iv) is achieved through a discretization on a time-grid, which introduces a discretization error. De Bortoli et al. (2021, Theorem 1) links these approximations to the total variation distance between the distribution of the generated samples from (iv) and PD. In our first methodological contribution we develop a procedure for constructing diffusion processes targeting PD without relying on time-reversal arguments. The proposed transport (coupling) between PZ and PD is achieved by: (1) specifying a diffusion process X on [0, τ ] starting from a generic x0; (2) conditioning X on hitting a generic xτ at time τ , thus obtaining a diffusion bridge; (3) taking a bivariate mixture Π0,τ of diffusion bridges over (x0, xτ ) with marginals Π0 = PZ and Πτ = PD, obtaining a mixture process M ; (4) matching the marginal distribution of M over [0, τ ] with a diffusion process, resulting in a diffusion with initial distribution PZ and terminal distribution PD.
The realized diffusion bridge mixture transport (DBMT) between PZ and PD is exact by construction. We thus sidestep the approximation common to all DDPM approaches due to the dependency from PD retained through the noising process. Moreover, the DBMT can be realized for almost arbitrary PZ , f( ·) and g( ·). This increased flexibility is a departure from the DTRT where f( ·) and g( ·) need to be chosen to obtain convergence toward a simple distribution PZ . Similarly to the DTRT, achieving the DBMT requires the computation of a drift adjustment term which depends onD. For a SDE class of interest, we develop a unified and interpretable representation of DTRT and DBMT drift adjustments as simple transformations of conditional expectations over D. This novel result provides insights on the target mapping that we aim to approximate and on the quality of approximation achieved by the trained score models of Song et al. (2021). Having defined for the DBMT a Fisher divergence objective similarly to Song et al. (2021), we leverage on this unified representation to define two additional training objectives featuring appealing properties.
In our last methodological contribution we extend the class of SDEs that can be realistically employed in computer vision applications. Specifically, computational considerations have so far restricted the transitions of the stochastic processes employed in DDPM to be fully factorials. We view images at a given resolution as taking values over a 2D lattice which discretizes the continuous coordinate system [0, 1]2 representing heights and widths. Diffusion processes are viewed as spatio-temporal processes with spatial support [0, 1]2. Doing so, it is possible to leverage on scalable simulation and inference techniques from spatial statistics and consider more realistic diffusion transitions.
This paper is structured as follows. In Section 2 we review the DTRT of Song et al. (2021) and in Section 3 we introduce the DBMT. In order to implement the DTRT and the DBMT it is necessary to specify the underlying SDE, i.e. the coefficients f( ·) and g( ·). We study a class of interest in Section 4. The unified view of drift adjustments is introduced in Section 5. Section 6 develops the training objectives and Section 7 reviews the obtained results and finalizes the DBMT construction. In Section 8 we establish the connection with spatio-temporal processes. We conclude in Section 9. Appendices A to D contain the theoretical framework, assumptions, proofs, and additional material.
Notation and conventions: we use uppercase notation for probability distributions (measures, laws) and lowercase notation for densities; each probability distribution, and corresponding density, is uniquely identified by its associated letter not by its arguments (which are muted); for example P (dx) is a distribution, p(x) is its corresponding density; random elements are always uppercase (an exception is made for times, always lowercase for typographical reasons); if P is the distribution of a stochastic process, we use subscript notation to refer to its finite dimensional distributions (densities with p), conditional or not, for some collection of times; for example pt′|t denotes a transition density, which is understood to be a function of four arguments pt′|t(y|x) = f(t, t′, x, y); δx is the delta distribution at x and ⊗ is used for product distributions; we refer directly to a given SDE instead of referring to the diffusion process satisfying such SDE when no ambiguity arises; we use [a]i and [A]i,j for vector and matrix indexing, A> for matrix transposition.

heading: 2 DIFFUSION TIME-REVERSAL TRANSPORT
text: The starting point of Song et al. (2021) is a diffusion process Y satisfying a generic D-dimensional time-inhomogenous SDE with initial distribution Y0 ∼ PD
dYr = f(Yr, r)dr + g(Yr, r)dWr, (1) over noising time r ∈ [0, τ ]. Thorough this paper we denote with Q the law of the diffusion solving (1) and with q the corresponding densities. Thus, let qr′|r(y|x), 0 ≤ r < r′ ≤ τ , be the transition density of (1), and let qr(y), 0 < r ≤ τ , be the marginal density of (1). As Y0 ∼ PD, we have
qr(y) = 1
N N∑ n=1 qr|0(y|x(n)). (2)
The dynamics of (1) over the reversed, i.e. sampling, time t = τ − r, t ∈ [0, τ ], are given by (Anderson, 1982; Haussmann & Pardoux, 1986; Millet et al., 1989)
dXt = [−f(Xt, r) +∇ ·G(Xt, r) +G(Xt, r)∇Xt ln qr(Xt)] dt+ g(Xt, r)dWt, (3) where r = τ − t is the remaining sampling time, G(x, r) = g(x, r)g(x, r)> and the D-dimensional vector ∇ ·G(x, r) is defined by [∇ ·G(x, r)]i = ∑D j=1∇xj [G(x, r)]i,j . That is the processes Xt
and Yr = Yτ−t have the same distribution. Approximating the terminal distribution Qτ of (1), i.e. the initial distribution of (3), with PZ , X0 is sampled from PZ and (3) is discretized and integrated over t to produce a sample Xτ approximately distributed as PD. The computation of the multiplicative drift adjustment ∇y ln qr(y) entering (3), i.e. the score of the marginal density (2), requires in principle O(N) operations. Let sφ(y, r) be a neural network for which we would like sφ(y, r) ≈ ∇y ln qr(y). It remains to find a suitable training objective for which unbiased gradients with respect to φ can be obtained at O(1) cost with respect to the dataset size N . As qr(y) has a mixture representation, the identity of Vincent (2011) for Fisher divergences provides us with the desired objective for a fixed r ∈ (0, τ ]
LFD,DTRT(φ, r) = E Yr∼Qr [∥∥∇Yr ln qr(Yr)− sφ(Yr, r)∥∥2] = E
(Y0,Yr)∼Q0,r [∥∥∇Yr ln qr|0(Yr|Y0)− sφ(Yr, r)∥∥2]. (4) The key point is that an unbiased, O(1) with respect to N , mini-batch Monte Carlo (MC) estimator for the expectation (4) can be trivially obtained by sampling a batch Y0 ∼ PD, Yr ∼ Qr|0(dyr|Y0), and evaluating the average loss over the batch. In order to achieve a global approximation over the whole time interval (0, τ ], Song et al. (2021) proposes uniform sampling of time r
LFD,DTRT(φ) = E r∼U(0,τ ],(Y0,Yr)∼Q0,r
[ Rr ∥∥∇Yr ln qr|0(Yr|Y0)− sφ(Yr, r)∥∥2], (5)
where Rr = E[‖∇Yr ln qr|0(Yr|Y0)‖2] −1 is a regularization term. A MC estimator for (5) is constructed by augmenting the MC estimator for (4) with the additional sampling step r ∼ U(0, τ ].

heading: 3 DIFFUSION BRIDGE MIXTURE TRANSPORT
text: Our starting point is a generic D-dimensional time-inhomogenous SDE which, in contrast to Song et al. (2021), is directly defined on the sampling time t ∈ [0, τ ]
dXt = f(Xt, t)dt+ g(Xt, t)dWt. (6) We reserve P · |0( · |x0) to denote the law of the diffusion solving (6) for a given starting value x0 and p · |0( · |x0) to denote the corresponding densities.

heading: 3.1 DIFFUSION BRIDGES
text: Diffusion bridges are central to the proposed methodology, in this Section we cover their basic theory. A diffusion bridge is a diffusion process starting from a given value which is conditioned on hitting a terminal value. It is a deep result, and consequence of Doob h-transforms (Särkkä & Solin (2019, Chapter 7.9), Rogers & Williams (2000, Chapter IV.6.39)), that a diffusion processes pinned down on both ends is still a diffusion process. In particular the Markov property is preserved. More precisely, (6) with initial value x0 conditioned on hitting a terminal value xτ at time τ is characterized the following SDE on [0, τ ] with initial value x0 (Särkkä & Solin, 2019, Theorem 7.11)
dXt = [ f(Xt, t) +G(Xt, t)∇Xt ln pτ |t(xτ |Xt) ] dt+ g(Xt, t)dWt, (7)
where G(x, t) = g(x, t)g(x, t)>. The multiplicative adjustment factor ∇xt ln pτ |t(xτ |xt) forces the process to hit xτ at time τ and the diffusion process solving (7) is known as the diffusion bridge from (x0, 0) to (xτ , τ). As previously noted, pτ |t(xτ |xt) in (7) refers to the transition density of (6).

heading: 3.2 DIFFUSION MIXTURES
text: The proposed transport construction relies on a representation result for diffusion mixtures. We present here an informal version and report the precise statement, the required assumptions, and the proof in Appendix A. Theorem 1 (Diffusion mixture representation — informal). Let {Xλ}, λ ∈ Λ be a collection of diffusions with associated SDEs {dXλt } and marginal densities {πλt }. Let L be a mixing distribution on Λ, πt be the L-mixture of {πλt }. Then there exists a diffusion process X with marginal πt. X follows a SDE whose drift and diffusion coefficients are weighted averages of the corresponding coefficients in {dXλt }, where the weights are proportional to {πλt } and to the mixing density.
Theorem 1 is first established in Brigo (2002, Corollary 1.3) limitedly to finite mixtures and 1- dimensional diffusions. The proof of Theorem 1 in Appendix A is more direct and extends the result to the required multivariate setting. In Section 3.1 we introduced diffusion bridges mapping arbitrary initial values x0 to arbitrary final values xτ . Let Π0,τ denote a generic bivariate distribution on RD × RD with marginals Π0,Πτ . We define the diffusion mixture M as the mixture of diffusion bridges corresponding to (X0, Xτ ) ∼ Π0,τ . That is, we apply Theorem 1 to the collection of diffusion bridges (7) indexed by their initial and terminal values, λ = (x0, xτ ), Λ = RD × RD, with mixing distribution L(dλ) = Π0,τ (dx0, dxτ ). By Theorem 1 the following SDE on [0, τ ] with initial distribution Π0 has the same marginal distribution as M , in particular its terminal distribution is Πτ
dXt = µ(Xt, t)dt+ g(Xt, t)dWt,
µ(xt, t) = f(xt, t) +G(xt, t) ∫ ∇xt ln pτ |t(xτ |xt)
pt|0,τ (xt|x0, xτ ) πt(xt)
Π0,τ (dx0, dxτ )︸ ︷︷ ︸ A(xt,t) ,
πt(xt) = ∫ pt|0,τ (xt|x0, xτ )Π0,τ (dx0, dxτ ).
(8)
In (8), A(xt, t) gives the multiplicative drift adjustment factor for (6). A case of particular interest occurs when Π0 puts all the mass on a single value x0. In the following we refer to A(xt, t, x0) in stance of A(xt, t), and to πt|0(xt|x0) in stance of πt(xt), when it is necessary to distinguish this specific case. We also extend the scope of Π to indicate the law of M . Indeed, we already denoted with Π0,τ its initial-terminal distribution, and with πt its marginal density. Accordingly, A(xt, t) = EXτ∼Π(dxτ |xt)[∇xt ln pτ |t(Xτ |xt)]. The transport from PZ to PD is then achieved by Πτ = PD and PZ = Π0 (PZ can be arbitrarily defined). As PD is an empirical distribution, the integral in (8) with respect to xτ reduces to averages over D. In summary, the diffusion X solution of (8) realizes the proposed transport from PZ to PD by matching the marginal distribution of M .

heading: 4 SDE CLASS
text: The starting point of the proposed transport is the unconstrained SDE (6). In this Section we define SDEs which are realized through a time-change of simpler SDEs and which are general enough to subsume the SDEs introduced in Song et al. (2021). Consider the D-dimensional SDEs
dZt = Γ 1/2dWt, (9)
dZt = αtZtdt+ Γ 1/2dWt, (10)
where αt 6= 0 is a scalar function and G(Xt, t) = Γ introduces an arbitrary covariance structure. (9) is the SDE of a correlated and scaled Brownian motion and (10) is the SDE of an Ornstein-Uhlenbeck process driven by a correlated and scaled Brownian motion. The transition densities of (9) and (10) are Gaussian (Appendix B). We denote both with p̃t′|t, informally (9) is a special case of (10) with αt = 0. SDE (10) in the time-homogenous case αt = −1/2 has stationary distribution ND(0,Γ). We now introduce the time-change. Let βt > 0 be a continuous function on [0, τ ]. Then bt = ∫ t 0 βudu defines a monotonically (strictly) increasing function bt : [0, τ ]→ [0, bτ ]. The following SDEs on [0, τ ] represent the class of dynamics for (1) and (6) on which we focus on the rest of this paper
dXt = √ βtΓ 1/2dWt, (11)
dXt = αtβtXtdt+ √ βtΓ
1/2dWt, (12) and G(x, t) = βtΓ. The standard time-change result for diffusions (Øksendal, 2003, Theorem 8.5.1) establishes that the processes Xt respectively from (11) and (12) are equivalent in law to their timescaled counterparts Zbt from (9) and (10). That is, SDEs (11) and (12) correspond to the evolution of the simpler SDEs (9) and (10) under a non-linear time wrapping where time flows with instantaneous intensity βt. For both (11) and (12) the time-change argument yields pτ |t(y|x) = p̃bτ |bt(y|x) for the transition density of (6), and equivalently for the transition density qτ |t of (1). We thus obtain
pτ |t(xτ |xt) = ND(xτ ; xta(t, τ), Γv(t, τ)) (13) for appropriate scalar functions a(t, τ), v(t, τ) with v(t, τ) > 0 (Appendix B). By direct computation
∇xt ln pτ |t(xτ |xt) = Γ−1 (
xτ a(t, τ)
− xt ) a2(t, τ)
v(t, τ) , (14)
∇xτ ln pτ |t(xτ |xt) = Γ−1 ( xta(t, τ)− xτ ) 1
v(t, τ) . (15)
From Bayes theorem and the Markov property we have
pt|0,τ (xt|x0, xτ ) = ND (xt; x0abr(0, t, τ) + xτabr(0, t, τ), Γvbr(0, t, τ)) , (16)
where once again abr(0, t, τ), abr(0, t, τ) and vbr(0, t, τ) > 0 are scalar functions given in Appendix B. Finally, by direct computation
∇xt ln pt|0,τ (xt|x0, xτ ) = Γ−1 x0abr(0, t, τ) + xτabr(0, t, τ)− xt
vbr(0, t, τ) . (17)
These results provide all the analytical formulas required for the computation of the adjustment factors A(xt, t), A(xt, t, x0) and of the training objectives used to approximate them (Section 6).

heading: 4.1 INTERPRETATION OF DENOISING TIME-REVERSED SDES
text: Song et al. (2021) introduces two specifications of (1), named VESDE and VPSDE, which are respectively given by
dYr = √ βve,rdWr, (18)
dYr = − 1
2 βvp,rYrdr +
√ βvp,rdWr. (19)
See Appendix B for the functional form of βve,r and βvp,r. We thus recover (18) and (19) from (11) and (12) with Γ = I and αt = −1/2. That is, VESDE and VPSDE correspond to a time change of the much simpler SDEs for the standard Brownian motion and for the standard Langevin SDE
dZr = dWr,
dZr = − 1
2 Zrdr + dWr.

heading: 5 UNIFIED VIEW OF DRIFT ADJUSTMENTS
text: The linearity of SDEs (11) and (12), underlying our and Song et al. (2021) works, has the important consequence that (14) and (15) are linear in xt. This in turn allow us to derive an alternative representation for the drift adjustment in (8). Indeed, substituting (14) in (8) gives (Appendix A)
G(x, t)A(x, t) = βt
( 1
a(t, τ) E Xτ∼Πτ|t(dxτ |x) [Xτ ]− x
) a2(t, τ)
v(t, τ) . (20)
Similarly, for the time-reversal drift adjustment term in (3) we have (Appendix A) G(x, r)∇x ln qr(x) = βr ( a(0, r) E
Xτ∼Q0|r(dxτ |x) [Xτ ]− x
) 1
v(0, r) . (21)
The relations (20) and (21) provide a unified view of the inner workings of the DTRT and of the DBMT targeting PD. In the following we always refer to sampling time t. Remember that r = τ − t is the remaining sampling time. For ease of exposition we assume βt = 1, as shown in Section 4 the term βt corresponds to a time-warping. The terms a(t, τ), a(0, r) are “integrated scalings”. They are equal to 1 for (11) and the same holds for (12) as r → 0. The terms v(t, τ), v(0, r) are “integrated variances”. They are equal to r for (11) and the same holds for (12) as r → 0. We commonly refer to E[Xτ |x, t] for expectation terms in (20) and (21). Both drift adjustments (20) and (21) are thus essentially of the form (E[Xτ |x, t]− x)v−1r where the term v−1r diverges as r → 0.
The expectations E[Xτ |x, t] are convex linear combinations of the samples x(n) from D. Explicitly, E[Xτ |x, t] = ∑N n=1 ω(x, t)
(n)x(n), where the weights ω(x, t)(n) are the probabilities, under the distributions Q (time-reversal sampling process (3)) and Π (mixture of diffusions process M from Section 3.2), of reaching each state x(n) at terminal time τ from x at time t. By construction, the initial weights entering expectation (20) are all equal to 1/N when X starts from a fixed value x0, and are so on average when X0 is stochastic. The initial weights entering expectation (21) are
on average approximately equal to 1/N, depending on the quality of the approximation PZ ≈ Qτ . Thus, E[Xτ |X0, 0] is an averaging of many samples x(n). As time progresses, changes in Xt correspond to changes in E[Xτ |Xt, t] through changes in the weights ω(x, t)(n). Eventually all mass concentrates on a single weight ω(x, t)(∗) corresponding to a dataset sample x(∗). Ultimately, the attractor dynamics implied by (E[Xτ |x, t]− x)v−1r drive Xt to x(∗). We provide an inspection in Figure 1, where D(CIFAR) stands for the training portion of the CIFAR10 dataset, and Euler(T) corresponds to the Euler scheme (Kloeden & Platen, 1992) applied with T discretization steps.
In the VESDE and VPSDE of Song et al. (2021) we have G(x, r) = βrI and reversing (21) gives
E Xτ∼Q0|r(dxτ |x)
[Xτ ] = v(0, r)∇x ln qr(x) + x
a(0, r) , (22)
where ∇y ln qr(y) is the true score. We can thus take a trained score model sφ(x, r) ≈ ∇x ln qr(y), plug it in (22), and verity the extent to which E[Xτ |x, t] has been approximated, see Figure 1.
We pause for a moment to review the findings of Figure 1 (see Appendix D for additional related plots). Firstly, we can classify the dynamics of E[Xτ |Xt, t] and of the associated weights in three stages. In the 1st stage the weights do not move much. During the 2nd stage, roughly t ∈ [0.4, 0.6], the weights’ mass gets distributed over a limited number of samples. Interestingly, the weights’ dynamics are not monotonic. As time progresses the weights’ mass shifts between different objects from different classes. From the beginning of the 3rd stage all mass gets allocated to a single weight, the terminal image is decided well in advance of the terminal time. These dynamics are suboptimal. We would like to shorten the 1st stage, but it is associated with large values of βt (i.e. quick time passing) which are required to decouple Qτ from PD. This is an intrinsic limitation of time-reversal approaches. It is also dubious that (partially) sampling multiple objects over t is beneficial for efficient generative modeling when we make use only of the terminal sample. This issue applies to trained models as well, as the 3rd row of Figure 1 shows. An interesting open question is how to obtain more suitable dynamics, where perhaps class transitions happen rarely. Secondly, E[Xτ |Xt, t] provides a denoised representation of Xt across the whole 3rd stage. An alternative to the noise removal step applied to Xτ in Song et al. (2021) is to consider E[Xτ |Xt, t] as the sampling process instead. Thirdly, Figure 1 makes it clear that lowering the number of discretization steps affects generative
sampling in multiple ways. On the one hand the terminal sample Xτ is more noisy. This is not very surprising: close to τ the drift adjustment is approximately (x(∗) − Xt)v−1r , which is the drift of a Brownian bridge. Bridge sampling is notoriously problematic (Bladt et al., 2016). On the other hand larger discretization errors also significantly affect the dynamics of E[Xτ |Xt, t] resulting in less coherent samples. We remark that none of these insights could have been gained by observing Xt alone, i.e. the even cells’ rows of Figure 1. To conclude, (20) and (21) give an additional meaning to “denoising”. Neural network approximators need to map from a noisy input Xt to an adjustment toward a smoother superimposition of samples. The desire to minimize the discrepancy between the smoothness properties of Xt and that of E[Xτ |Xt, t] motivates the developments of Section 8.

heading: 6 TRANSPORTS APPROXIMATION
text: As in Song et al. (2021), computing the multiplicative drift adjustment A(xt, t) requires O(N) operations. In this Section we introduce three training objectives for which unbiased and scalable, i.e. O(1) with respect to N , MC estimators can be immediately derived. The first training objective applies only to A(xt, t, x0). It relies on the identity (Appendix A)
A(xt, t, x0) = ∇xt lnπt|0(xt|x0)−∇xt ln pt|0(xt|x0). (23)
It is advantageous to consider the right-hand side of (23) because from (8) we know that πt|0(xt|x0) has mixture representation. As in Song et al. (2021), we can rely on Vincent (2011) to obtain a scalable objective to train a neural network approximator sφ(xt, t) ≈ ∇xt lnπt|0(xt|x0), i.e.
LFD,DBMT(φ) = E t∼U(0,τ),Xt∼Πt|0
[ Jt ∥∥∇Xt lnπt|0(Xt|x0)− sφ(Xt, t)∥∥2]
= E t∼U(0,τ),(Xt,Xτ )∼Πt,τ|0
[ Jt ∥∥∇Xt ln pt|0,τ (Xt|x0, Xτ )− sφ(Xt, t)∥∥2], (24)
where Jt = E[‖∇Xt ln pt|0,τ (Xt|x0, Xτ )‖2] −1 is a regularization term.
The remaining training objectives rely on the identities (20) and (21). The goal is directly approximate the expectations of (20) and (21) which, as in Section 5, we denote with a generic E[Xτ |x, t]. That is, we aim to train a neural network approximator sφ(x, t) ≈ E[Xτ |x, t]. As conditional expectations are mean squared error minimizers, suitable objectives for the expectation terms of (20) and (21) are
LCE,DBMT(φ) = E t∼U [0,τ),(Xt,Xτ )∼Πt,τ [∥∥Xτ − sφ(Xt, t)∥∥2], (25) LCE,DTRT(φ) = E
r∼U [0,τ),(Y0,Yr)∼Q0,r [∥∥Y0 − sφ(Yr, r)∥∥2]. (26) In Table 1 we summarize the operations needed to implement the plain MC estimators for the four objectives considered in this work. We reference where to find the required quantities for SDEs (11) and (12). The MC estimators for the Fisher divergence losses LFD,∗ involve multiplications by Γ−1 (by (15) and (17)). Moreover, computing the drift adjustment at generation time requires multiplications by Γ. In Section 8 we discuss how to manage the computational burden. An appealing property of LCE,∗ is that computing the drift adjustment only requires the application of simple scalar functions (see (20) and (21)), and that their MC estimators only requires sampling operations. A further advantage of LCE,∗ is that no regularization is required. In contrast, in the absence of regularization terms, LFD,∗ are divergent for t ≈ τ due to the term v−1r (Section 5).

heading: 7 DBMT OVERVIEW AND NUMERICAL EXPERIMENT
text: In this section we finalize the DBMT construction, putting together the results of Sections 3, 4 and 6. The unconstrained SDE follows (11) or (12). It remains to choose the mixing distribution Π0,τ . The marginal Πτ needs to match PD, but there is flexibility in the choice of Π0|τ . Song et al. (2021) derived a random ordinary differential equation (RODE) matching the marginal distribution of a generative SDE, leading to faster sampling and to likelihood evaluation. RODE-matching requires Π0 to have density. A natural implementation is given by the factorial distribution Π0,τ = PZ ⊗ PD with PZ = ND(0,Γ) and the unconstrained SDE following (12) with αt = 1/2 which preserves PZ . If instead the DBMT starts from a fixed value x0, i.e. Π0,τ = δx0 ⊗ PD, we can choose x0 = 1/N ∑N n=1 x
(n)a(0, τ)−1 to remove the drift adjustment at t = 0 (see (20)) and reduce the work required to transport x0 to PD. Finally, the use of non-factorial distributions can lead to a more efficient implementation, by linking the initial distribution to PD. The training steps for the simplest objective (25) of Section 6 are reported in Algorithm 1. Batch size is assumed to be 1 to ease the description. It is also assumed that Π0,τ is factorial, otherwise the obvious modification applies to line 2 (Algorithm 2 is unaffected) where the endpoints are sampled. At line 3 a random central time and the corresponding state are sampled. The function optimizationstep implements a step of stochastic gradient descent update based on the loss L. The corresponding sampling algorithm is reported in Algorithm 2 where the Euler(T) discretization is assumed in line 6. Pt|0,τ (dxt|X0, Xτ ), a(t, τ), v(t, τ), βt are defined in Section 4. Section 8 shows how to sample efficiently from Pt|0,τ (dxt|X0, Xτ ) and ND(0,Γ) in computer vision applications.
Algorithm 1 DBMT training (LCE,DBMT) Input: PD, PZ , SDE (11) or (12), NN sφ(x, t) Output: trained sφ(x, t)
1: repeat 2: Xτ ∼ PD, X0 ∼ PZ 3: t ∼ U [0, τ), Xt ∼ Pt|0,τ (dxt|X0, Xτ ) 4: L ←
∥∥Xτ − sφ(Xt, t)∥∥2 5: φ← optimizationstep(φ,L) 6: until convergence
Algorithm 2 DBMT sampling (LCE,DBMT) Input: PZ , SDE (11) or (12), trained sφ(x, t) Output: Discretized path X0:T
1: X0 ∼ PZ 2: for s = 1, . . . , T do 3: t← (s− 1) τT , x← Xs−1 4: us ← βt ( 1 a(t,τ)sφ(x, t)− x ) a2(t,τ) v(t,τ) 5: Es ∼ ND(0,Γ) 6: Xs ← x+(f(x, t)+us) τT +g(x, t) √ τ T Es 7: end for
We consider a toy numerical example with PZ = PD = 1/3(δ−2 + δ0 + δ2), D = τ = 1. The unconstrained SDE follows the standard Brownian motion. We consider two mixing distributions: independent mixing Π⊥⊥0,1 where X0 and X1 are independent and fully dependent mixing Π = 0,1 where X0 = X1. The results are reported in Figure 2. For both couplings the correct terminal distribution PD is recovered, as can be seen by taking the row-wise sum of the transition matrices. The initialterminal distribution of X solving (8), which realizes the DBMT, is different from the corresponding mixing distribution Π0,1 which is realized by the mixture process M . Π⊥⊥0,1 results in a transition matrix of equal entries 1/9, Π⊥⊥0,1 results in a diagonal transition matrix of equal diagonal entries 1/3.

heading: 8 NON-DENOISING DIFFUSIONS
text: In computer vision applications, images of resolution H×W corresponds to D = 3HW . The use of an arbitrary covariance matrix Γ in (11) and (12) requires its Cholesky (or equivalent) decomposition with cost O(D3). As the resolution increases the computational burden gets intractable very quickly. Indeed, to the best of the authors’ knowledge, all prior DDPM literature only considers independent transitions, that is Γ = I . We suggest to view SDEs (11) and (12) as corresponding to the space discretization on an H×W grid of a spatio-temporal process defined over the spatial domain [0, 1]2. Consider the Euler discretization of (12): Xt+∆t = Xt+αtβtXt∆t+ √ ∆tEt, where Et ∼ ND(0,Γ),
the idea is to adopt a functional perspective: X(t+ ∆t, s) = X(t, s) +αtβtX(t, s)∆t+ √
∆tE(t, s) where s ∈ [0, 1]2 defines space coordinates. That is, both X(t) and E(t) at each time t are random processes over [0, 1]2. We assume the innovations E(t) to be a Gaussian process (GP) for each t. As the GPs E(t) are defined on a 2D domain we can leverage on scalable inference techniques from spatial statistics. As an example, we consider the circulant embedding method (CEM) (Wood & Chan, 1994; Dietrich & Newsam, 1997) which exploits a connection with the fast Fourier transform (FFT). See Appendix C for a cursory review of the CEM. Consider an H×W uniform grid S of size S = HW discretizing [0, 1]2, i.e. the support of images. For a stationary covariance function the CEM samples E(t) on S with cost O(D ln(D)). This is close to the O(D) cost of sampling from a pure white-noise process, and compares very favorably to the O(D3) cost of a Cholesky decomposition. One limitation of CEM is that generated samples, while always Gaussian, might not have the correct covariances. Whether this happens, and in that case the quality of the approximation, depends on the covariance function. In Appendix C we select and fit an isotropic GP to the microscale properties of D(CIFAR). For this estimated GP sampling is exact. Figure 3 shows samples from a pure white-noise GP, i.e. Γ = I , (1st row) and from the fitted GP using CEM (2nd row). As noted in Section 6, sampling is enough to implement the MC estimators for LCE,∗, but the MC estimators and drift adjustments for LCE,∗ involve additional matrix multiplications by Γ and Γ−1. The CEM allows to compute these at the same O(D ln(D)) cost if we define the GP E(t) on a 2D torus (Rue & Held, 2005, Chapter 2.1). This corresponds to introducing dependencies between “opposing” boundaries of [0, 1]2. Figure 3 (3rd row) shows some samples, in the highlighted patch the opposing-boundaries dependency is evident. Either way, all samples from the 2nd and 3rd rows of Figure 3 match the smoothness properties of D(CIFAR).

heading: 9 CONCLUSIONS
text: The DBMT construction of Section 3 is exact. The SDE class of Section 4 is tractable as it results in linear diffusion bridges. The time-space factorization of the diffusion coefficient g(x, t) = √ βtΓ 1/2 separates modeling concerns: βt corresponds to a time-wrapping, Γ can be efficiently modeled by fitting the microscale properties of PD. Availability of GPU-accelerated FFT implementations motivates our focus on the CEM. Alternative scalable approaches abound, from Gaussian Markov Random Fields (Rue & Tjelmeland, 2002; Rue, 2001) to Karhunen–Loève expansions (Betz et al., 2014). It remains to apply the results of this work to perform an empirical benchmarking. Section 6 develops three novel training objectives, two of which with desirable properties compared to the objective of Song et al. (2021), especially for non-factorial transitions. We remark the simplicity of the proposed DBMT approach (Algorithms 1 and 2) compared to alternatives grounded in the Schrödinger bridge problem (De Bortoli et al., 2021; Wang et al., 2021; Vargas et al., 2021). The understanding of the target mappings ((20) and (21)) can guide the development of neural networks more closely matching the target structure compared to the U-Net default choice. https://github.com/? links to the code accompanying this paper which is made available under the MIT license.

heading: A THEORETICAL FRAMEWORK

heading: A.1 ASSUMPTIONS
text: Assumption 1 (SDE solution). A given D-dimensional SDE(f, g) with associated initial distribution V0 and integration interval [0, τ ] admits a unique strong solution on [0, τ ].
Assumption 1 can be checked through the application of the standard existence and uniqueness theorems for SDE solutions. Of particular relevance to our setting is the formulation of Krylov (1995, Chapter 5, Theorem 1) that limits the monotonic requirement to, informally speaking, drifts that pull the process toward infinities.
Assumption 2 (SDE density). A given D-dimensional SDE(f, g) with associated initial distribution V0 and integration interval [0, τ ] admits a marginal / transition density on (0, τ) with respect to the D-dimensional Lebesgue measure that uniquely satisfies the Fokker-Plank / Kolmogorov-forward partial differential equation (PDE).
We refer to Särkkä & Solin (2019, Chapter 5) and to Karatzas & Shreve (1996, Chapter 5.7) for connections between SDEs and PDEs.
All theoretical results of this work rely on simple algebraic manipulations and re-arrangements of quantities of interest. The main complication stems from the need to justify differentiation and integration exchanges, i.e. exchange of limits.
Assumption 3 (exchange of limits). We assume that limits exchanges are justified in the steps marked with (?) and (??).
Similarly, various steps in the derivations involve considering fractional quantities with densities appearing in the denominators.
Assumption 4 (positivity). For a given stochastic process, all finite-dimensional densities, conditional or not, are strictly positive.
Assumption 4 is easy to verify. We resorted to the practical but somewhat unsatisfactory formulation of Assumption 3 because in full generality it is complicated to give easy to check conditions. We just note that when ΠT = PD, the limit exchange marked with (?) is always justified. So are the limits exchanges marked with (??) when in addition Π0 puts all the mass to a fixed initial value, or (by direct verification) when Π0 is Gaussian for the SDE class of Section 4. Thorough this paper, both in the main text and in the proofs that follow, it is supposed that Assumptions 1, 2 and 4 are satisfied by SDEs (1), (6) and (7). This is the case for the SDE class of Section 4, i.e. (11) and (12), for any Π0 with finite variance.
Remark: For ease of exposition it is assumed thorough this paper that all diffusions take values in the state space RD. There is no impediment in extending the presented results to the case of diffusions taking values in a subset X ⊂ RD. The obvious changes to Assumptions 1 to 4 apply, the proofs carry over without substantial modifications. This extension could be of practical interest as images are often represented as floating point values in [0, 1].

heading: A.2 STATEMENT AND PROOF OF DIFFUSION MIXTURE REPRESENTATION THEOREM
text: Theorem 2 (Diffusion mixture representation). Consider the family of D-dimensional SDEs on t ∈ [0, τ ] indexed by λ ∈ Λ
dXλt = µ λ(Xλt , t)dt+ σ λ(Xλt , t)dW λ t , Xλ0 ∼ Vλ0 , (27)
where the initial distributions Vλ0 and the BMs Wλt are all independent. Let νλt , t ∈ (0, τ) denote the marginal density of Xλt . For a generic mixing distribution L on Λ, define the mixture marginal density νt for t ∈ (0, τ) and the mixture initial distribution V0 by
νt(x) = ∫ Λ νλt (x)L(dλ), V0(dx) = ∫ Λ Vλ0 (dx)L(dλ). (28)
Consider the D-dimensional SDE on t ∈ [0, τ ] defined by
µ(x, t) =
∫ Λ µλ(x, t)νλt (x)L(dλ)
νt(x) ,
σ(x, t) =
∫ Λ σλ(x, t)νλt (x)L(dλ)
νt(x) ,
dXt = µ(Xt, t)dt+ σ(Xt, t)dWt,
Y0 ∼ V0.
(29)
It is assumed that all diffusion processes Xλ and the diffusion process X solving (29) satisfy the regularity assumptions Assumptions 1, 2 and 4 and that Assumption 3 holds. Then the marginal distribution of the diffusion X is νt.
Proof of Theorem 2. We start by establishing that the law of X is indeed given by the solution of (29). In this proof we make use of the following notation: for f scalar-valued (f)t = d dtf , for a
vector-valued (a)x = ∑D i=1 d dxi a, for A matrix-valued (A)xx = ∑D i,j=1 d2 dxidxj A. This notation allows for a compact representation of PDEs reminiscent of the 1-dimensional setting. Then, for 0 < t < τ we have that
(ν(x, t))t = (∫ Λ νλ(x, t)L(dλ) ) t
= ∫ Λ ( νλ(x, t) ) t L(dλ) (??)
= ∫ Λ ( µλ(x, t)νλ(x, t) ) x + 1 2 ( σλ(x, t)νλ(x, t) ) xx L(dλ)
= ∫ Λ ( µλ(x, t)νλ(x, t) ν(x, t) ν(x, t) ) x + 1 2 ( σλ(x, t)νλ(x, t) ν(x, t) ν(x, t) ) xx L(dλ)
= (∫ Λ µλ(x, t)νλ(x, t) ν(x, t) L(dλ)ν(x, t) ) x + 1 2 (∫ Λ σλ(x, t)νλ(x, t) ν(x, t) L(dλ)ν(x, t) ) xx . (??)
The second line is an exchange of limits, the third line is the application of the Fokker-Plank PDEs for the collection of processes Xλ, the fourth line is a rewriting in terms of ν(y, t), the last line is another exchange of limits. The result follows by noticing that the last line gives the Fokker-Plank representation of (29).

heading: A.3 DRIFT ADJUSTMENTS
text: Limitedly to this section, we lighten the notation by removing subscripts from probability measures and densities. The missing time points can be inferred without ambiguity from the variables.

heading: A.3.1 DRIFT ADJUSTMENT IDENTITIES FOR CONSTANT INITIAL VALUE
text: First identity: ∇xt ln ∫
p(xτ |xt) p(xτ |x0) Π(dxτ )
= ∫ ∇xtp(xτ |xt) p(xτ |x0) p(xt|x0)Π(dxτ ) /∫ p(xτ |xt) p(xτ |x0) p(xt|x0)Π(dxτ ) (?)
= ∫ ∇xt ln p(xτ |xt)
p(xτ , xt|x0) p(xτ |x0) Π(dxτ ) /∫ p(xτ , xt|x0) p(xτ |x0) Π(dxτ )
= ∫ ∇xt ln p(xτ |xt)p(xt|x0, xτ )Π(dxτ ) / π(xt|x0)
= A(xt, t, x0).
Second identity: ∇xt ln ∫
p(xτ |xt) p(xτ |x0) Π(dxτ )
= ∇xt ln ∫
p(xτ |xt) p(xτ |x0) p(xt|x0)Π(dxτ )−∇xt ln p(xt|x0)
= ∇xt lnπ(xt|x0)−∇xt ln p(xt|x0).

heading: A.3.2 DRIFT ADJUSTMENTS AS EXPECTATIONS
text: To establish (20) notice that from (14) we have
G(xt, t)A(xt, t) = βtΓΓ −1 ∫ (
xτ a(t, τ)
− xt ) a2(t, τ)
v(t, τ) p(xt|x0, xτ ) π(xt) Π0,τ (dx0, dxτ )
= βt
( 1
a(t, τ)
∫ xτ π(xt|x0, xτ )
π(xt) Π0,τ (dx0, dxτ )− xt
) a2(t, τ)
v(t, τ)
= βt
( 1
a(t, τ) E Xτ∼Π(dxτ |xt) [Xτ ]− xt
) a2(t, τ)
v(t, τ) .
To establish (21) note that from (15) we have G(yr, r)∇yr ln q(yr) = Γ ∫ ∇yr ln q(yr|y0)
q(yr|y0) q(yr) PD(dy0)
= βrΓΓ −1 ∫ (
a(0, r)y0 − yr v(0, r) ) q(yr|y0) q(yr) PD(dy0)
= βr
( a(0, r) ∫ y0 q(yr|y0) q(yr) PD(dy0)− yr ) 1 v(0, r)
= βr
( a(0, r) E
Y0∼Q(dx0|yr) [Y0]− yr
) 1
v(0, r) .

heading: B SDES CLASS FORMULAS
text: The transition densities of (9) and (10) are given respectively by
p̃bm,τ |t(zτ |zt) = ND (zτ ; zt, Γ(τ − t)) , p̃ou,τ |t(zτ |zt) = ND ( zτ ; zte αt:τ (τ−t), Γ ( 1
2αt e2αt:τ (τ−t) − 1 2ατ
)) .
Here we used the notation f t:τ = 1 τ−t ∫ τ t fudu, i.e. f t:τ is the average value of a function fu on the interval [t, τ ]. The time-homogenous case of (10), where αt:τ = α, is thus immediately recovered.
The scalar functions abm(t, τ), aou(t, τ), vbm(t, τ) and vou(t, τ) are given by
abm(t, τ) = 1, vbm(t, τ) = bτ − bt,
aou(t, τ) = e αbt:bτ (bτ−bt), vou(t, τ) =
1 2αbt e2αbt:bτ (bτ−bt) − 1 2αbτ .
The scalar functions vbr(0, t, τ), abr(0, t, τ) and abr(0, t, τ) are given by
vbr(0, t, τ) = v(0, t)v(t, τ)
v(0, t)a2(t, τ) + v(t, τ) ,
abr(0, t, τ) = v(t, τ)a(0, t)
v(0, t)a2(t, τ) + v(t, τ) ,
abr(0, t, τ) = v(0, t)a(t, τ)
v(0, t)a2(t, τ) + v(t, τ) .
The scalar functions βve,r and βvp,r are given by
βve,r = σ 2 min ( σmax σmin )2r 2 log σmax σmin , βvp,r = ( β̄min + r ( β̄max − β̄min )) .
The constants σmin, σmax, β̄min, β̄max depend in part on the dataset considered, but are consistently chosen to have βve,r, βvp,r small for r ≈ 0 and large for r ≈ τ .
The approximating distributions PZ in Song et al. (2021) are PveZ = ND(0, Iσ2max) for VESDE, PvpZ = ND(0, I) for VPSDE.

heading: C ADDITIONAL MATERIAL

heading: C.1 CLOSELY RELATED WORK
text: A work closely related to the present paper is that of Wang et al. (2021) as it similarly avoids the time-reversal construction. Wang et al. (2021) construct a 2-stages diffusion process from a constant initial value x0 to PD by relying on the theory of Schrödinger bridges. The most notable differences with respect to the DBMT transport are: (i) the dynamics considered in Wang et al. (2021) are less general, in our notation they correspond to f( ·) = 0, g( ·) = σI for a fixed scalar σ; (ii) the transport proposed in Wang et al. (2021) necessarily starts from x0, the general result of (8) allows for (almost) arbitrary initial distributions and initial-terminal dependencies. For an initial x0, i.e. for case of A(xt, tx0) in Section 3.2 and for the more limited dynamics considered in Wang et al. (2021), the achieved transport is the same. In this sense the DBMT generalizes the first stage diffusion of Wang et al. (2021). It is interesting to note that in the case of a constant x0 the DBMT can also be obtained by an application of Doob h-transforms as we show in the following section.
We now review two additional works grounded in the Schrödinger bridge problem: De Bortoli et al. (2021); Vargas et al. (2021). Both works rely on the Iterative Proportional Fitting (IPF) procedure to solve the (dynamic) Schrödinger bridge problem. Both works leverage on time-reversal results to carry out the alternated Schrödinger half-bridge IPF iterations. The main difference between the two works is that De Bortoli et al. (2021) estimates the optimal SDE drifts via neural network approximations and score-matching, while Vargas et al. (2021) relies on Gaussian Processes and maximum likelihood fitting. The work of De Bortoli et al. (2021) can be seen as an extension of Song et al. (2021), and similarly to our work allows the use of shorter time intervals. Compared to our proposal, it solves a harder problem but also presents additional difficulties. Training is more involved as all the neural network approximations, one for each IPF iterate, need to converge. Moreover, there is limited guidance on how to optimally choose the number of integration steps over the number of IPF iterates.

heading: C.2 CONNECTION WITH DOOB H-TRANSFORMS
text: The previously established identity A(xt, t, x0) = ∇xt ln ∫
pτ |t(xτ |xt) pτ |0(xτ |x0) Πτ (dxτ ),
shows that the drift adjustment can be equivalently expressed as µ(xt, t) = f(xt, t) +G(xt, t)∇xth(xt, t), h(xt, t) = ln ∫
pτ |t(xτ |xt) pτ |0(xτ |x0) Πτ (dxτ ),
as x0 is a constant. It can be verified that the h function satisfies the required space-time regularity property (Särkkä & Solin, 2019, Eq. (7.73)). As such, it is a genuine Doob h-transform. That pht′|t(xt′ |xt) = pt′|t(xt′ |xt)h(xt′ , t
′)/h(xt, t) is the transition density of the DBMT transport from δx0 to Πτ follows by direct computation.

heading: C.3 GP MODELLING ON CIFAR10
text: For simplicity, we assume a factorial distribution over the channels and an isotopic stationary covariance function. We rely on the semivariogram approach (Cressie, 1993) to compare how
different covariance functions fit D(CIFAR). A semivariogram is a measure of dependency across space. In the case of an isotropic stationary covariance it simplifies to a scalar function of the Euclidean distance between points: γ(‖∆s‖) = E[(∆xs)2]/2 with ∆xs = xs+∆s − xs. The rate of decrease of γ(‖∆s‖) toward 0 as ‖∆s‖ → 0 gives a measure of the infinitesimal spatial dependency, i.e. the smoothness of the spatial process. Semivariograms corresponding to different covariance functions are here fitted to their empirical counterparts via a weighted minimum-leastsquares procedure (Cressie, 1993). Figure 4 illustrates the exponential and RBF semivariogram fits for two images of D(CIFAR). The exponential covariance, which corresponds to rougher paths, provides a much better fit to the shown samples. This result is consistent across D(CIFAR). We remark that Γ = I corresponds to a pure white-noise process with a perfectly flat semivariogram which would clearly result in a very poor fit to the empirical semivariograms shown in Figure 4. Based on these findings, we model the innovations of each image channel as a GP with exponential covariance function with length-scale θ = 0.205, the median estimated value (Figure 4 (right)). We match the marginal variance to that of D(CIFAR), σ2 = 0.063.

heading: C.4 CIRCULANT EMBEDDING METHOD
text: We start by providing a cursory explanation leading to efficient sampling in the 1D case, before giving the intuition behind the extension to the 2D case. We refer to Wood & Chan (1994); Dietrich & Newsam (1997) for a complete explanation and to Rue & Held (2005) for the results underlying efficient density (likelihood) computation. Let [0, 1] be the spatial domain of interest. Let S = {si}Mi=1 be a uniform grid (regular lattice) discretizing [0, 1], where the points si are assumed to be ordered. The first key observation is that for a stationary covariance function ρ( · , ·) the covariance matrix C with entries Ci,j = ρ(sj , sj) is symmetric and Toeplitz, i.e. with constant-diagonals. See Figure 5 (leftmost) for an example where M = 9. A property of symmetric Toeplitz matrices is that they can always be embedded in larger symmetric circulant matrices. A circulant matrix of size M ′×M ′ is defined by the property that all its rows (and columns) are obtained by cycling through the same M ′-dimensional vector. Circulant matrices correspond to covariance matrices of GPs defined on a (here 1D) torus (Rue & Held, 2005, Chapter 2.1). The circulant embedding matrix just introduced corresponds to an artificial enlargement of the spatial domain [0, 1] to a larger interval leading to a torus. See Figure 5 (2nd from left) for a circulant embedding of C. The second key observation is that a circulant matrix is diagonalized by the 1D FFT matrix. Having obtained the eigenvalues of C, efficient sampling on the enlarged domain is achieved by the 1D FFT applied to
complex standard random numbers multiplied by the (square root of the) eigenvalues. The real and imaginary part of the generated samples are independent. The main issue with the CEM is that the circulant matrix embedding might fail to be positive definite. The issue can be avoided by considering progressively larger embeddings, see the theoretical and empirical findings of Dietrich & Newsam (1997). Otherwise, a level of approximation can be accepted by modifying the covariance function or by truncating the eigenvalues to be positive.
The development of the 2D CEM follows very similar steps. The domain of interest is now [0, 1]2, the uniform grid is S = {si,j}Mi,j=1 and the points si,j are assumed to be lexicographically ordered. The stationarity of the covariance functions results in a symmetric block-Toeplitz covariance matrix, as shown in Figure 5 (3rd from left). Again, symmetric block-Toeplitz matrices can be embedded in symmetric block-circulant matrices, as exemplified by Figure 5 (rightmost, zooming might be required to see the block structure). Block-circulant matrices can be shown to be diagonalized by the 2D FFT matrix, and efficient sampling follows from similar steps to the ones seen in the 1D case.

heading: D ADDITIONAL FIGURES

Review summaries:
The paper studies the Brownian bridge formulation for diffusion-based generative models, derives theory for it, and shows connections to earlier works. Some spatial developments are also discussed. This is a theoretical paper that only has rudimentary experiments. The overall idea of Brownian Bridge diffusion is interesting, but the idea’s merits are not shown. The paper suffers from high technical complexity making it difficult to digest.  

The paper introduces a methodological framework for generative modeling through diffusion processes without time-reversal arguments. By utilizing diffusion bridges the authors consider processes that start and end in pre-determined points $x_0$ and $a_{\tau}$. Then by extending this to mixtures of diffusion bridges they show how to transport a distribution $\Pi_0 = P_D$ to $\Pi_{\tau} = P_Z$ without time-reversal of the diffusion process. For specific classes of SDEs by drawing connections to the paradigm introduced by Song et al 2021 and provide a unified view on the drift adjustment in forward and backward SDEs. This allows for defining a time and state-dependent class probability function (given by conditional expectation) providing further insight into the inner-workings of DTRT and DBMT. The authors then introduce two other training objectives $L_{\text{FD}},L_{\text{CE}}$ and discuss their favorable properties. Finally, methods for moving beyond fully factorial noise models are presented by extending the process and the noise term from functions of time to be functions of both time and space. This allows for incorporating "priors" with better spatial dependency models more suitable for specific datasets. 

In this paper the authors introduce a new model for diffusion-based generative modelling.  Instead of relying on time-reversed processes as in existing works, they propose an approach based on the mixing of diffusion bridges. They show that a mixing of diffusion bridges remains a diffusion with closed-form drift and volatility. Similarly to existing score-matching based generative modelling works the drift in this new formulation can be expressed as a conditional expectation and therefore can be seen as the minimizer of some loss function, allowing the use of neural-network based approximations. In addition, the authors propose some unification of the SDE classes used in [1] and to use non-identity volatility matrices.  [1] Song, Sohl-Dickstein, Kingma, Kumar, Ermon, Poole -- Score-based Generative Modeling through Stochastic Differential Equations