Schrodinger Bridge Flow for Unpaired Data Translation

Authors: valentin  de bortoli
google  deepmind
iryna  korshunova
andriy  mnih
arnaud  doucet
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 17:51:00

Abstract: Mass transport problems arise in many areas of machine learning whereby one wants to compute a map transporting one distribution to another. Generative modeling techniques like Generative Adversarial Networks (GANs) and Denoising Diffusion Models (DMMs) have been successfully adapted to solve such transport problems, resulting in CycleGAN and Bridge Matching respectively. However, these methods do not approximate Optimal Transport (OT) maps, which are known to have desirable properties. Existing techniques approximating OT maps for high-dimensional data-rich problems, including DDMs-based Rectified Flow and Schrodinger bridge procedures, require fully training a DDM-type model at each iteration, or use mini-batch techniques which can introduce significant errors. We propose a novel algorithm to compute the Schrodinger bridge, a dynamic entropy-regularized version of OT, that eliminates the need to train multiple DDMs-like models. This algorithm corresponds to a discretization of a flow of path measures, referred to as the Schrodinger Bridge Flow, whose only stationary point is the Schrodinger bridge. We demonstrate the performance of our algorithm on a variety of unpaired data translation tasks.

Main paper: 
 introduction: 
 The problem of finding a map to transport one probability distribution to another one has numerous applications in machine learning. In particular, it is at the core of generative modeling where the idea is to transform a noise distribution into the data distribution, and is also central to transfer learning tasks such as image-to-image translation. For discrete probability distributions, it is possible to compute the Optimal Transport (OT) map but this is computationally expensive (Peyré et al., 2019) . By showing that an entropy-regularised version of OT, the Entropic OT (EOT), could be computed much more efficiently using the Sinkhorn algorithm, Cuturi (2013) has enabled transport ideas to be used in numerous applications (Ge et al., 2021; Zhou et al., 2022) . However, the computational complexity of Sinkhorn algorithm is quadratic in the sample size, which makes its application to very large datasets impractical. Mini-batch versions have been proposed, see e.g. (Genevay et al., 2018) , but tend to introduce significant errors in high dimensions (Sommerfeld et al., 2019) .In the context of generative modeling, Denoising Diffusion Models (DDMs) (Song et al., 2021a; Ho et al., 2020) have shown impressive performance in a variety of domains. DDMs define a forward process progressively noising the data, and sample generation is achieved by approximating the timereversal of this diffusion. In order to leverage the iterative refinement properties of DDMs in the OT setting, methods exploiting the equivalence between the static versions of (E)OT and their dynamic counterparts (Benamou and Brenier, 2000; Léonard, 2014) have been developed. A procedure to approximate the dynamic OT is considered by Liu et al. (2023b) , while techniques to approximate the dynamic equivalent to EOT, the Schrödinger Bridge (SB), have been proposed in (De Bortoli et al., 2021; Vargas et al., 2021; Chen et al., 2022; Peluchetti, 2023; Shi et al., 2023) . These techniques are expensive however, as they require training multiple DDM-type models. Mini-batch versions of OT and Sinkhorn (Pooladian et al., 2023; Tong et al., 2024b) combined with bridge or flow matching have also been proposed to approximate the OT path and SB, but they optimise a minibatch OT objective that can introduce significant errors in high dimensions: the error in Wasserstein-1 distance is of order O(B -1/(2d) ), where d is the dimension of the problem and B the minibatch size, see (Sommerfeld et al., 2019, Corollary 1) .In this paper, we propose a novel approach to computing the SB. Similarly to Iterative Markovian Fitting (IMF) and its practical implementation, Diffusion Schrödinger Bridge Matching (DSBM) (Shi et al., 2023; Peluchetti, 2023) , it leverages the fact that the SB is the only Markov process with prescribed marginals at the endpoints which is in the reciprocal class of the Brownian motion, i.e. it has the same bridge as the Brownian motion (Léonard, 2014) ; see Section 2 for more details on Markov processes and the reciprocal class. Compared to DSBM, our approach is easier to implement as it does not require caching samples, alternating between optimising two different losses, and, optionally, uses one neural network instead of two. In Section 3, we start by introducing a flow of path measures whose time-discretisation yields a family of algorithms called α-IMF and presented in Section 4. Notably, we show that α-IMF converges to the Schrödinger Bridge for any α ∈ (0, 1]. Additionally, for a special value of the discretisation stepsize α = 1, we recover the IMF procedure (Peluchetti, 2023; Shi et al., 2023) , while α < 1 corresponds to online versions of IMF. We implement a parametric version of the α-IMF as an online DSBM procedure, called α-DSBM. We illustrate the efficiency of our approach in unpaired image-to-image translation settings in Section 6.Notation. We denote the space of path measures by P(C), i.e. P(C) = P(C([0, 1], R d )), where C([0, 1], R d ) is the space of continuous functions from [0, 1] to R d . The subset of Markov path measures associated with a diffusion of the form dX t = v t (X t )dt + σ t dB t , with σ, v locally Lipschitz, is denoted M. For Q induced by ( √ εB t ) t∈[0,1] , with ε > 0 and (B t ) t≥0 a d-dimensional Brownian motion, the reciprocal class of Q is denoted R(Q), see Definition 2.1. For any P ∈ P(C), we denote by P t its marginal distribution at time t, P s,t the joint distribution at times s, t, P s|t the conditional distribution at time s given the state at time t, and P |0,1 ∈ P(C) the distribution of the path on time interval (0, 1) given its endpoints; e.g. Q |0,1 is a scaled Brownian bridge. Unless specified otherwise, all gradient operators ∇ are w.r.t. the variable x t with time index t. Given probability spaces (X, X ) and (Y, Y), a Markov kernel K : X × Y → [0, 1], and a probability measure µ defined on X , we write µK for the probability measure on Y such that for any A ∈ Y we have µK(A) = X K(x, A)dµ(x). In particular, for any joint distribution Π 0,1 over R d × R d , we denote the mixture of bridges measure as Π = Π 0,1 P |0,1 ∈ P(C), which is short for Π(•) = R d ×R d P |0,1 (•|x 0 , x 1 )dΠ 0,1 (x 0 , x 1 ). Finally, we define the Kullback-Leibler (KL) divergence between two probability measures π 0 , π 1 ∈ P(X) as KL(π 0 |π 1 ) = X log((dπ 0 /dπ 1 )(x))dπ 0 (x) if π 0 is absolutely continuous w.r.t. π 1 and KL(π 0 |π 1 ) = +∞ otherwise. 
optimal transport and schrödinger bridge: 
 Unpaired Transfer and Optimal Transport. Given unpaired data samples from π 0 and π 1 , where π 0 , π 1 are two distributions on R d , we are interested in designing a transport map from π 0 to π 1 . This corresponds to an unpaired data transfer task. We can formulate this problem as finding a distribution Π on R d × R d with marginals Π 0 = π 0 and Π 1 = π 1 so that if X 0 ∼ π 0 then X 1 |X 0 ∼ Π 1|0 (•|X 0 ) satisfies X 1 ∼ π 1 . Among an infinite number of such so-called coupling distributions Π, we are here interested in finding the Entropic Optimal Transport (EOT) coupling Π ⋆ defined asEQUATIONwhere H(Π) is the differential entropy of Π and ε > 0 is a regularisation hyperparameter (Peyré et al., 2019) . For ε = 0, we recover the standard OT.In order to leverage the recent advances in generative modeling, and in particular the concept of iterative refinement central to DDMs, we turn to a dynamic formulation of EOT known as the Schrödinger Bridge problem (Léonard, 2014) . It is defined as follows: find P ⋆ ∈ P(C) such that P ⋆ = argmin P∈P(C) {KL(P|Q) ; P 0 = π 0 ,EQUATIONwith Q ∈ P(C) induced by a scaled d-dimensional Brownian motion ( √ εB t ) t∈ [0, 1] . The term dynamic here refers to the fact that (2) is defined on path measures, i.e. on (stochastic) processes, in contrast to the static problem (1) which is defined on measures on the space R d × R d . In Section 3, we show that solving (2) is equivalent to optimising the vector field of a stochastic process using objectives similar to the ones of bridge matching (Peluchetti, 2021; Albergo and Vanden-Eijnden, 2023; Lipman et al., 2023; Liu et al., 2023a) . Under mild assumptions, it can be shown that P ⋆ 0,1 = Π ⋆ , see e.g. (Léonard, 2014; Pavon et al., 2021) . Hence solving (1) reduces to solving (2). Once we have found P ⋆ associated with (X ⋆ t ) t∈[0,1] , we can sample from P ⋆ by first sampling X ⋆ 0 ∼ π 0 and then sampling the trajectory (X ⋆ t ) t∈(0,1] which yields (X ⋆ 0 , X ⋆ 1 ) ∼ Π ⋆ .Reciprocal and Markov projections. To introduce our methodology, it is necessary to recall the notions of reciprocal and Markov projections. We refer to Shi et al. (2023) for more details. For practitioners, a more intuitive explanation of these projections is given in Appendix E.Definition 2.1 (Reciprocal projection):P ∈ P(C) is in the reciprocal class R(Q) of Q if P = P 0,1 Q |0,1 .We define the reciprocal projection of P ∈ P(C) as P ⋆ = proj R(Q) (P) = P 0,1 Q |0,1 . We will write proj R instead of proj R(Q) to simplify notation.In other words, P is in the reciprocal class of Q if the conditional distribution of a path given its endpoints is identical under P and Q, see (Roelly, 2013) . Sampling from the reciprocal projection of P can be achieved by sampling a path (X t ) t∈[0,1] from P, keeping only the values of the endpoints, say X 0 , X 1 , and then sampling a new value for the bridge (X t ) t∈(0,1) from Q |0,1 .Definition 2.2 (Markov projection): Assume that Q is induced by ( √ εB t ) t∈[0,1] for ε > 0. Then, when it is well-defined, for any P ∈ R(Q), the Markovian projection M = proj M (P) ∈ M is the path measure induced by the diffusion (X ⋆ t ) t∈[0,1] with for any t ∈ [0, 1]dX ⋆ t = v ⋆ t (X ⋆ t )dt + √ εdB t , v ⋆ t (x t ) = E P 1|t [X 1 | X t = x t ] -x t /(1 -t), X ⋆ 0 ∼ P 0 .In practice, implementing a Markovian projection requires solving a regression problem to approximate E P 1|t [X 1 | X t = x t ], similar to the one appearing in bridge matching and flow matching.One key property of the Markovian projection is that M t = P t for all t ∈ [0, 1], i.e. the Markovian projection preserves the marginals; see (Peluchetti, 2021) for instance.Iterative Markovian Fitting. Leveraging the reciprocal and Markovian projections, Peluchetti (2023) and Shi et al. (2023) concurrently introduced IMF. Starting from P0 = (π 0 ⊗ π 1 )Q |0,1 , a measure where endpoints are sampled independently from π 0 and π 1 and then interpolated using a (scaled) Brownian bridge, they define a sequence of path measures (P n , Pn ) n∈N where P n = proj M ( Pn ) and Pn+1 = proj R (P n ). This ensures that P n 0 = π 0 , P n 1 = π 1 for all n, and it can be shown that the sequence (P n ) n∈N converges to the SB, see (Peluchetti, 2023, Theorem 2) . The practical implementation of this algorithm proposed by Shi et al. (2023) is called DSBM. Implementing DSBM poses challenges, as each Markovian projection requires training a neural network to approximate the relevant conditional expectations by minimising a bridge matching loss. Furthermore, in practice, generated model samples are stored in a cache in order to train the next iterations of DSBM. This introduces additional hyperparameters that require tuning. In Section 3 we propose α-IMF, an algorithm which can be interpreted as the discretisation of a flow of path measures. This leads to α-DSBM, an algorithm that is computationally much more efficient than DSBM as it does not rely on a Markovian projection at each step. 
schrödinger bridge flow: 
 We will now introduce a flow of path measures (P s ) s≥0 , and show that the time-discretisation of this flow with an appropriate stepsize α ∈ (0, 1] yields a family of procedures called α-IMF, which all converge to the Schrödinger Bridge. While α = 1 yields the classical IMF, α ∈ (0, 1) yields an incremental version of IMF. In Section 4 we show that α-IMF can be implemented as an online version of DSBM.Let (P s , Ps ) s≥0 be a flow of path measures defined for any s ≥ 0 by which we assume is well-defined. Note that for any s ≥ 0, P s is Markov while Ps is in the reciprocal class of Q. Crucially, the only fixed point of (3) is the SB. Indeed, let P be a fixed point of (P s ) s≥0 in (3). Then, we have that P = proj R (proj M ( P)). Hence, we get P = proj R (proj M (. . . (proj R (proj M ( P))) . . . )).EQUATIONHence, under mild assumptions, P is a limit point of IMF and therefore P is the SB P ⋆ given by (2), see (Peluchetti, 2023, Theorem 2) .Next, for any α ∈ (0, 1], we define the following discretisation of (3) called α-IMF:EQUATION)and P n = proj M ( Pn ). Note that for any n ∈ N, Pn ∈ R(Q). This recovers the IMF procedure (Shi et al., 2023; Peluchetti, 2023) when α = 1. Using the definition of the sequence ( Pn ) n∈N , it is possible to analyse the sequence (P n ) n∈N using the properties of the KL divergence as well as the Pythagorean identities derived in (Shi et al., 2023; Peluchetti, 2023) . We first introduce some assumptions on the Schrödinger Bridge problem. We recall that the differential entropy of a probability measure π is given byH(π) = - R d log((dπ/dLeb)(x))dπ(x),if π admits a density with respect to the Lebesgue measure and +∞ otherwise. Recall that Q is associated with ( √ εB t ) t∈[0,1] and assume that Q0 = Leb. Let π 0 , π 1 ∈ P(R d ) such that R d ∥x∥ 2 dπ i (x) < +∞, H(π i ) < +∞,for i ∈ {0, 1}. Under these assumptions, we can use the characterisation of the SB as the only path measure that preserves π 0 , π 1 , and is both Markov and in the reciprocal class of Q (see e.g. (Léonard, 2014, Theorem 2.12) ). We get the following result.Theorem 3.1 (Convergence of α-IMF): Let α ∈ (0, 1] and (P n , Pn ) n∈N defined by (4). Under mild assumptions, we have that lim n→+∞ P n = P ⋆ , where P ⋆ is the solution of the Schrödinger Bridge problem (2).We show here that α-IMF is associated with an incremental version of DSBM for α ∈ (0, 1).v : [0, 1] × R d → R d , we introduce the loss function L(v, P) = 1 0 L t (v t , P)dt = 1 0 (R d ) 3 v t (x t ) - x 1 -x t 1 -t 2 dP 0,1 (x 0 , x 1 )dQ t|0,1 (x t |x 0 , x 1 )dt,(5) where we recall that Q is induced by ( √ εB t ) t∈[0,1] for some ε > 0. This loss was already considered in (Peluchetti, 2021; Lipman et al., 2023; Liu et al., 2023a; Liu, 2022; Shi et al., 2023) . We also define the path measure P v ∈ P(C) associated withEQUATIONConsider first the sequence (v n ) n∈N defined byEQUATIONUsing Definition 2.2, we have that P v n+1 = proj M (proj R (P v n )), which corresponds to P n+1 in the IMF sequence. Therefore we have that lim n→+∞ P v n = P ⋆ under mild assumptions (Peluchetti, 2023, Theorem 2) .Functional gradient descent. We now introduce a relaxation of ( 7), where, instead of considering the argmin, we update the vector field with one gradient step. To define this relaxation, we recall that for a functional F : F → R, where F is an appropriate function space, its functional derivative (Courant and Hilbert, 2008) with reference measure µ is denoted ∇ µ F and is given for any ϕ ∈ F, when it exists, bylim γ→0 (F (f + γϕ) -F (f ))/γ = ⟨∇ µ F (f )(x), ϕ(x)⟩dµ(x).Initialised with v 0 t (x) = (E P0 1|t [X 1 | X t = x] -x)/(1 -t), where P0 = (π 0 ⊗ π 1 )Q |0,EQUATIONwith δ n > 0 and µ n ∈ P(C). The parameters (δ n , µ n ) n∈N will be made explicit in Proposition 3.2. We emphasize that, in contrast to the IMF procedure, in the online update (8) we do not need to solve a Markovian projection problem at every step; instead we simply take a gradient step on the loss (5).Connection with α-IMF. The following proposition shows that (P v n ) n∈N defined by ( 8) is associated with α-IMF defined in (4).Proposition 3.2 (Non-parametric updates are α-IMF): Let α ∈ (0, 1], (P n , Pn ) n∈N as in (4),δ n = α and µ n = (1 -α) Pn + αproj R (P n ).Then, under mild assumptions, we have P v n = P n for all n ∈ N.Combining Theorem 3.1 to Proposition 3.2, we get that lim n→+∞ P v n = P ⋆ , i.e. the non-parametric procedure converges to the SB. 
α-diffusion schrödinger bridge matching: 
 From DSBM to α-DSBM. In Section 3, we introduced α-IMF, a scheme which defines a sequence of path measures converging to the SB for all α ∈ (0, 1]. For α = 1, this corresponds to the IMF, whose practical DSBM implementation (Shi et al., 2023) requires repeatedly solving an expensive minimisation problem (7). In contrast, for α < 1 we are only required to take one (non-parametric) gradient step to update the vector field, see (8). This suggests the following practical implementation of α-IMF, called α-DSBM: First, pretrain a bridge matching model so that for t ∈ [0, 1] andx ∈ R d , v θ t (x) = (E P0 1|t [X 1 | X t = x] -x)/(1 -t), where P0 = (π 0 ⊗ π 1 )Q |0,1 .Then, perform the parametric version of the update (8):EQUATIONwhere Pθ is a stop-gradient version of P v θ . In Appendix D.2, we give a theoretical justification for this parametric equivalent of ( 5) and ( 8) by showing that, as α → 0, the update on the velocity fields v θ given by (9) corresponds to a direction of descent for the non-parametric loss (8) on average. Once again, we emphasize that if we replace the gradient step in (9) with the minimisation θ ← argmin θ L(θ, Pθ), we recover DSBM.Bidirectional online procedure. As with DSBM, directly implementing (9) leads to error quickly accumulating, see Appendix I for details. One way to circumvent this error accumulation issue is to consider a bidirectional procedure, in which we train both a forward and a backward model. This is possible because the Markovian projection coincides for forward and backward path measures, see (Shi et al., 2023, Proposition 9) . This suggests considering the loss L(v , v , P , P ) = 1 0 L t (v t , v t , P , P )dt, which is an extension of (5), whereEQUATION+ (R d ) 3 v 1-t (x t ) - x 0 -x t t 2 dP 0,1 (x 0 , x 1 )dQ t|0,1 (x t |x 0 , x 1 ).Similarly to (6), we define P v , P v , associated with (X t ) t∈[0,1] and (Y 1-t ) t∈[0,1] respectively, which are defined by forward and backward SDEs (fwd):dX t = v t (X t )dt+ √ εdB t , X 0 ∼ π 0 , (bwd): dY t = v t (Y t )dt+ √ εdB t , Y 0 ∼ π 1 . (11)Similarly to (8), we define non-parametric updates for any n ∈ N, t ∈ [0, 1] andEQUATION, P v n, , P v n, (x). We have the following proposition which ensures our bidirectional procedure is still valid and that the results of Proposition 3.2 still hold.Proposition 4.1 (Bidirectional updates): Let α ∈ (0, 1]. For any n ∈ N, define (P n , Pn ) n∈N by (4). Then, under mild assumption and assuming that δ n = α and µ n = (1 -α) Pn + αproj R (P n ), we have that for any n ∈ N, P v n, = P v n, = P n .In Appendix I, we show that in the Gaussian setting the bidirectional procedure (4.1) does not accumulate error when the vector field is approximated, while the unidirectional one (8) does.Vector field parameterisation. Contrary to existing procedures (Shi et al., 2023; Peluchetti, 2023; Liu, 2022) , we do not parameterise v and v using two separate networks. Instead, we consider an additional input s ∈ {0, 1} such that v θ (1, •) ≈ v and v θ (0, •) ≈ v . This allows us to substantially reduce the number of parameters in the model. The conditioning on s in the network is detailed in Appendix K. Before stating our full algorithm in Algorithm 1, we introduce a batched parametric version of (10). For ease of notation, we write Interp t for the operation corresponding to sampling from Q t|0,1 , i.e.EQUATIONWe are now ready to introduce the batched parametric version of (10). For a given batch of inputs X 1:B 0 and X 1:B 1 , timesteps t ∼ Unif([0, 1]) ⊗B , and X t = Interp t (X 0 , X 1 , Z) with Z ∼ N (0, Id) ⊗B , we compute the empirical forward and backward losses asEQUATION) ℓ (θ; t, X 0 , X t ) = 1 B B i=1 v θ 0, 1 -t i , X i t -X i 0 -X i t /t i 2 .We present the resulting α-DSBM in Algorithm 1. Note that in this algorithm, we maintain an Exponential Moving Average (EMA) of model parameters, as is common in diffusion models (Nichol and Dhariwal, 2021) . During the finetuning stage, when we generate samples to use as model's inputs, we then have a choice of sampling using the EMA or non-EMA parameters. At test time, we always sample using the EMA parameters, as it is known to improve the visual quality (Song and Ermon, 2020) . In Algorithm 1, we specify α ∈ (0, 1] as a stepsize parameter. In practice, we use Adam (Kingma and Ba, 2015) for optimization, thus the choice of α is implicit and adaptive throughout the training. To emphasize the importance of the parameter α, we sweep over its value with an explicit solver SGD in a toy setting, see Appendix K.2. We refer to Appendix K for more details on our experimental setup. Sample (X 0 , X 1 ) ∼ (π 0 ⊗ π 1 ) ⊗B 4: Sample t ∼ Unif([0, 1]) ⊗B and Z ∼ N (0, Id) ⊗B and compute X t = Interp t (X 0 , X 1 , Z)Update θ with a gradient step on 1 2 ℓ t 1:b , X1:b 1 , X 1:b t + ℓ t b+1:B , X b+1:B 0 , X b+1:B t 6:Update EMA parameters: θ EMA = γθ EMA + (1 -γ)θ 7: end for 8: for n ∈ {1, . . . , N finetuning } do9: Sample (X 0 , X 1 ) ∼ (π 0 ⊗ π 1 ) ⊗b 10: Sample X1 solving forward SDE (11)-(fwd) with v θ EMA (1, •) or v θ (1, •) starting from X 0 11: Sample X0 solving backward SDE (11)-(bwd) with v θ EMA (0, •) or v θ (0, •) starting from X 1 12: Sample t ∼ Unif([0, 1]) ⊗b and Z ∼ N (0, Id) ⊗b and compute X t = Interp t ( X0 , X 1 , Z ) 13: Sample t ∼ Unif([0, 1]) ⊗b and Z ∼ N (0, Id) ⊗b and compute X t = Interp t (X 0 , X1 , Z ) 14:Update θ with a gradient step on (Schrödinger, 1932) have been thoroughly studied through the lens of probability theory (Léonard, 2014) and stochastic control (Dai Pra, 1991; Chen et al., 2021) . They recently found applications in generative modeling and related fields leveraging recent advances in diffusion models (De Bortoli et al., 2021; Vargas et al., 2021; Chen et al., 2022) . Extensions of these methods to other machine learning problems and modalities were studied in (Shi et al., 2022; Thornton et al., 2022; Liu et al., 2022; Chen et al., 2023; Tamir et al., 2023) . Shi et al. (2023) ; Peluchetti (2023) concurrently introduced the DSBM algorithm which relies on a new procedure called IMF, while the DSB algorithm introduced in (De Bortoli et al., 2021) is based on the standard Iterative Proportional Fitting (IPF) scheme. Neklyudov et al. (2023a,b) ; Liu et al. (2022) generalise DSBM to arbitrary cost functions, albeit at the expense of having to learn the reciprocal projection which is no longer given by a Brownian bridge. These new methodologies translate to improved numerics when compared to their IPF counterparts, but they remain reliant on alternating between the optimisation of two losses. Finally, we note that the Schrödinger Bridge flow and the α-IMF procedure can be linked to the Sinkhorn flow recently introduced by Karimi et al. (2024) , see Appendix H.1 for a detailed discussion.Sampling-free methodologies. Sampling-free methodologies have been proposed to solve OT related objectives. In (Liu et al., 2023a; Somnath et al., 2023; Diefenbacher et al., 2024; Cao et al., 2024) , the authors perform one step of DSBM, i.e. only consider the pretraining stage of our algorithm. While the obtained bridge might enjoy transport properties, it does not solve an OT problem. In another line of work, Pooladian et al. (2023); Tong et al. (2024a,b) ; Eyring et al. (2024) have proposed simulation-free methods to minimise OT objectives. However, they target not the OT problem, but a minibatch version of it which coincides with OT only in the limit of infinite batch size, see (Pooladian et al., 2023, Theorem 4.2) . Other sampling-free methods to solve the Schrödinger Bridge problem include Kim et al. (2024) ; Gushchin et al. (2024b) both of which rely on adversarial losses to solve the OT problem. In (De Bortoli et al., 2021; Vargas et al., 2021; Liu et al., 2022; Shi et al., 2023; Peluchetti, 2023) the adversarial objective is dropped and instead the procedure requires alternating objectives during training and is not sampling-free. We also highlight the line of work of Korotin et al. (2024) ; Gushchin et al. (2024a) in which the Schrödinger Bridge potentials are parameterised with mixtures of Gaussians, allowing for fast training in small dimensions. Finally, recently Deng et al. (2024) introduced a variation on Schrödinger Bridge for generative modeling, which while still not sampling-free, does not require learning a forward process.  
experiments: 
 In this section, we illustrate the efficiency of α-DSBM on different tasks. In Section 6.1, we compare α-DSBM to DSBM in a Gaussian setting where the EOT coupling is tractable and show that α-DSBM recovers the solution faster than DSBM. In Section 6.2, we illustrate the scalability of our method through a range of unpaired image translation experiments.We compare α-DSBM to DSBM in the Gaussian setting where π 0 = N (0, σ 2 0 Id), π 1 = N (0, σ 2 1 Id) and Q is associated with ( √ εB t ) t∈[0,1] with √ ε = 0.5. In this case, the EOT coupling is N (0, Σ ⋆ ), with Σ ⋆ given byΣ ⋆ = σ 2 0 Id σ 2 ⋆ Id σ 2 ⋆ Id σ 2 1 Id, whereσ 2 ⋆ = (1/2)((σ 2 0 σ 2 1 + ε 2 ) 1/2 -ε),with Id being a d × d identity matrix. We consider d = 50, σ 0 = σ 1 = 1, resulting in σ 2 ⋆ ≈ 0.88. To showcase the robustness of α-DSBM, we consider the initial coupling P 0,1 , where (X 0 , X 1 ) ∼ P 0,1 , X 0 ∼ N (0, Id), X 1 = -X 0 , and let P0 = P 0,1 Q |0,1 . In this setting, the base model, i.e. bridge matching, significantly underestimates the true covariance σ 2 ⋆ , as shown in Section 6.1. Additionally, the figure illustrates that online finetuning approaches the true solution faster than the original iterative DSBM finetuning. For the latter, we can set how often we alternate between updating the forward and backward networks, and as this frequency increases, the behaviour approaches that of the online finetuning.Full covariance Gaussian case.Let π 0 = N (µ 0 , Σ 0 ), π 1 = N (µ 1 , Σ 1 ) with Σ i = Id + 1 2 Z i Z ⊤ ifor i ∈ {0, 1} and Z 0 , Z 1 independent d × d matrices with unit Gaussian entries. We also set µ 0 = µ 1 = 0. We consider the Entropic Optimal Transport (EOT) with regularization σ = 0.5 and d = 3, given byΠ = N (µ ⋆ , Σ ⋆ ), Σ ⋆ = Σ 0 C ⋆ C ⊤ ⋆ Σ 1 , with C ⋆ = 1 2 [Σ 1/2 0 D ⋆ Σ -1/2 0 -σ 2 Id], with D ⋆ = (4Σ 1/2 0 Σ 1 Σ 1/2 0 + σ 4 Id) 1/2. Let normFrob = ∥A -B∥ Fro /∥A∥ Fro be the normalized Frobenius distance between matrices A and B. The results are presented in Section 6.1 and confirm those presented in the original manuscript considering a diagonal covariance. Similarly to Shi et al. (2023) , we apply our method to image translation problems, such as MNIST digits to EMNIST letters (LeCun and Cortes, 2010; Cohen et al., 2017) , Wild to Cat domains from the Animal Faces-HQ (AFHQ) dataset (Choi et al., 2020) , downsampled to 64 × 64 and 256 × 256 resolutions and CelebA 64 × 64.The whole training procedure can be framed as a two-stage process: first, we train a base model on the true data samples, performing bridge matching (Peluchetti, 2021; Albergo and Vanden-Eijnden, 2023; Lipman et al., 2023; Liu et al., 2023a) , and then we finetune this model. We compare models that combine different vector field parameterisations (two networks vs. one bidirectional net), finetuning methods (iterative vs. online), and sample generation strategies during the finetuning stage.Following the established practice (Choi et al., 2020) , we evaluate our models using FID (Heusel et al., 2017) for visual quality, and mean squared distance (MSD) or LPIPS (Zhang et al., 2018) for alignment. It is important to note that for image translation tasks at hand, FID scores are not ideal, as FID was designed for natural RGB images, which is not the case for MNIST. It is also not well suited for small sample sizes as it is the case with AFHQ, where the test set in each domain has fewer than 500 examples. Thus quantitative results in Table 1 should be interpreted cautiously, and we recommend a visual inspection of samples to complement these quantitative measures, especially for the AFHQ models. Samples from the models along with the training and evaluation protocols are given in Appendix K.Compared to the iterative DSBM, our online finetuning α-DSBM reduces the number of tunable hyperparameters, i.e. inner and outer iterations, refresh rate and the size of the cache for storing generated samples. This simplifies implementation and makes the algorithm more practical. The primary remaining hyperparameter, the variance of a Brownian motion ε, requires careful tuning as it influences the trade-off between the visual quality and alignment, as was also observed in Shi et al. (2023) . An appropriate ε needs to balance the two: setting ε too low results in poor visual quality, while high values of ε cause poorly aligned and oversmoothed samples. Figure 3 illustrates how FID and MSD metrics vary with ε for the case of MNIST. Additionally, it demonstrates the impact of ε on the generated samples for the AFHQ-64 model.We run α-DSBM on CelebA with image size 64 × 64 with σ = 2.0. We do not change the training hyper-parameters compared to AFHQ. Visual results are reported in Figure 5 and Figure 6 . In Figure 5 , we show the influence of σ during the pretraining. The visual quality of the transfer is much lower for σ = 0 than for σ = 2.0. The case σ = 0 corresponds to the first step of Rectified Flow (i.e. Flow Matching). Given the poor quality of the samples, we do not perform finetuning with σ = 0. In Figure 6 , we compare the visual quality and alignment of DSBM and α-DSBM after 4000 training steps, corresponding to two outer DSBM iterations. In this case DSBM is trained with a bidirectional network and both procedures consist of finetuning the pretrained model obtained with σ = 2.0. We note that the alignment is better in the case of α-DSBM.  
discussion: 
 In this paper we have introduced α-Diffusion Schrödinger Bridge Matching (α-DSBM), a new methodology to solve Entropic Optimal Transport problems. α-DSBM is an improved version of DSBM, which does not require training multiple DDM-type models. We have shown that a non-parametric version of this method recovers the Schrödinger Bridge (SB). In addition, α-DSBM is easier to implement than existing SB methodologies while exhibiting similar performance. We illustrated the efficiency of our algorithm on a variety of unpaired transfer tasks.While α-DSBM solves one of the most critical limitations of DBSM, namely the alternative optimisation, several issues remain to be addressed in order for the method to scale comparably to generative DDMs. In particular, the method is not sampling-free, as during training it requires sampling from the model from the previous iteration to obtain the training data for the current iteration. While it seems difficult to derive a completely sampling-free method to solve SB problems without resorting to the Minibatch OT approximation, there is still room for improvement. The supplementary material is organised as follows. First in Appendix B, we analyze an Euclidean counterpart to the α-IMF sequence identified in Section 3 and the associated flow. In Appendix C, we show that the Markovian projection can be recovered as the parameterisation of the vector field that minimises the accumulation of errors, extending the results of Chen et al. (2024) . Theoretical results are gathered in Appendix D. In particular in Appendix D.1 we show that the proposed non-parametric method coincides with the α-IMF and prove the convergence of the α-IMF. In Appendix D.2, we show the connection between the non-parametric and the parametric updates. In Appendix E, we provide more background on DSBM and propose an extension of the DSBM methodology. Consistency losses similar to (Daras et al., 2023b; De Bortoli et al., 2024) are proposed in Appendix F. Model stitching procedures are described in Appendix G. We comment on extended related work in Appendix H. In particular we draw connections with Sinkhorn flows (Karimi et al., 2024) , Reinforcement Learning policies, Expectation-Maximisation schemes following (Brekelmans and Neklyudov, 2023) and comment on finetuning of diffusion models. In Appendix I, we investigate the accumulation of bias in a Gaussian setting and compare forward-forward and forward-backward methods. In Appendix J, we derive the preconditioning of loss following the principles of (Karras et al., 2022) in the case of bridge matching. Additional results and experimental details are presented in Appendix K.In this section, we study a simplified counterpart of the Schrödinger flow and of DSBM in a Euclidean setting. The goal of this section is to draw some conclusions in the Euclidean case which also remain true empirically when analyzing the Schrödinger Bridge problem.We consider the setA 1 = {(x, y) ∈ R 2 : y ≥ x} and the set A 2 = {(x, y) ∈ R 2 : y ≤ 0}.Loosely speaking, one can identify A 1 with the reciprocal class R(Q) and A 2 with the set of Markov path measures proj M . In that case, we have that for any(x, y) ∈ R 2 , proj A1 ((x, y)) = ((x + y)/2, (x + y)/2) if (x, y) / ∈ A 1 and otherwise, proj A1 ((x, y)) = (x, y).In addition, we have that for any (x, y) ∈ R 2 , proj A2 ((x, y)) = (x, 0) if (x, y) / ∈ A 2 and proj A2 ((x, y)) = (x, y) otherwise. We consider the following flow (x t , y t ) t≥0 given by ∂ t (x t , y t ) = proj A1 (proj A2 ((x t , y t ))) -(x t , y t ).Let (x 0 , y 0 ) / ∈ A 1 and (x 0 , y 0 ) / ∈ A 2 . Denote T the explosion time of (x t , y t ), i.e. for any t ≥ T we have that (x t , y t ) = ∞, where R 2 ∪ {∞} is the one-point compactification of R 2 . Finally, denote τ ≤ T such that for any t ∈ [0, τ ], (x t , y t ) ̸ ∈ A 1 and (x t , y t ) ̸ ∈ A 2 . Then, we have∂ t (x t , y t ) = (-x t /2, x t /2 -y t ).Hence, we have thatx t = x 0 exp[-t/2] for any t ∈ [0, τ ] and y t = x 0 exp[-t/2] + (x 0 exp[-t/2]) 2 (y 0 -x 0 )/x 2 0 .Therefore, we get that τ = T = +∞ and we have that for any t ≥ 0x t = x 0 exp[-t/2], y t = x t + x 2 t (y 0 -x 0 )/x 2 0 .Hence, ((x t , y t )) t≥0 converges exponentially fast to (0, 0) with rate 1/2.We now investigate the rate of convergence of the alternate projection scheme, i.e. the Euclidean equivalent of DSBM. We define ((x n , y n )) n∈N such that for any n ∈ N,(x n+1 , y n+1 ) = proj A1 (proj A2 ((x n , y n ))) = (x n /2, 0).Hence, we get that x n = x 0 2 -n and therefore ((x n , y n )) n∈N converges exponentially fast to (0, 0). Note that this procedure corresponds to a discretisation of the flow ((x t , y t )) n∈N with stepsize α = 1.More generally, we define for any α ∈ (0, 1], ((x α n , y α n )) n∈N such that for any n ∈ N,(x α n+1 , y α n+1 ) = αproj A1 (proj A2 ((x α n , y α n ))) + (1 -α)(x α n , y α n ).Hence, we get that x n = x 0 2 -n and therefore ((x n , y n )) n∈N converges exponentially fast to (0, 0). It can be shown that for any n ∈ N, x α n = x α 0 (1 -α/2) n and in addition, we have thaty α n = (1 -α) n y α 0 + αx α 0 n-1 k=0 (1 -α) k (1 -α/2) n-k .Therefore, we get thaty n = (1 -α) n y α 0 + 2(1 -(1 -α/(2 -α)) n )(1 -α/2) n x α 0 .We now analyse the complexity of the different discretisations assuming that the cost of discretising the flow with stepsize α ∈ (0, 1] is C α . In that case in order to reach the threshold value ε, i.e.|x α n | ≤ ε, we get a total cost C α n = O(log(1/ε)C α / log(1/(1 -α/2))), where we have neglected the terms that do not depend on log(1/ε). Hence, if C α is constant then the choice α = 1 is the best possible one in the range α ∈ (0, 1]. Otherwise, one has to consider the ratio C α / log(1/(1 -α/2)), where the lower is the better. The flow procedure and the iterative based one are presented in Figure 1 .Based on this simplified Euclidean experiment, we draw some conclusions which also remain true in our setting, see Appendix K for more experimental details. First, we have that different discretisations of the flow yield different convergence rates. Large stepsizes incur faster convergence. This suggests to choose α = 1. However, if the cost of choosing α = 1 is too high then one might turn to alternative schemes with α ∈ (0, 1) assuming that C α < C 1 in that case. To draw a parallel with our setting, in the case of DSBM (case α = 1), we need to solve the projection subproblem at each step which incurs a great cost. On the other hand, one step of the online algorithm only requires sampling once from the model and performing one gradient step.For a given non-Markovian (stochastic) interpolant process (see definition below), there exist an infinite number of Markov processes admitting the same marginals (Albergo and Vanden-Eijnden, 2023) . In this section, when it is well-defined, we show that the Markovian projection corresponds to the process which minimises an error measure (defined further) in case one has access to the oracle ofx t → E[X 1 | X t = x t ].Stochastic Interpolant. We first start by recalling the framework of Albergo and Vanden-Eijnden (2023) . Consider a coupling Π between π 0 and π 1 , one builds a (stochastic) flow between π 0 and π 1 using the following interpolation procedureX t = Interp t (X 0 , X 1 , Z) = α t X 0 + β t X 1 + γ t Z, (X 0 , X 1 ) ∼ Π, Z ∼ N (0, Id),where α 1 = β 0 = γ 0 = γ 1 = 0 and α 0 = β 1 = 1. This defines a non-Markovian process. We denote by π t the induced unconditional distribution of X t . Let us now consider the Markov process (X ε t ) given byEQUATIONwhere 14) can be rewritten as(B t ) t∈[0,1] is a d-dimensional) t∈[0,1] is a (stochastic) flow mapping π 0 onto π 1 . Note that (X ε t ) t∈[0,1] in (EQUATIONIn the specific case where α t = 1 -t, β t = t and γ t = σ 0 t(1 -t) then (15) becomesX t = Interp t (X 0 , X 1 , Z) = (1 -t)X 0 + tX 1 + σ 0 t(1 -t)Z.This corresponds to the marginal distribution of the bridge associated with (σ 0 B t ) t∈ [0, 1] . In this case, (15) becomesdX ε t = E [(X 1 -X t )/(1 -t)|X t = X ε t ] dt + √ 2σ 0 dB t for ε 2 t = (2γ t )( γt -γ t αt /α t ) = 2σ 2 0 .In Proposition C.1, we will show that this choice of (ε t ) t∈[0,1] is optimal in some sense.Consider ( Xε t ) t∈[0,1] given by d Xε t = ( αt /α t ) Xε t +( βt -β t αt /α t )E[X 1 | X t = Xε t ]+( γt -γ t αt /α t -ε 2 t /(2γ t )) Ẑ(t, Xε t )+ε t dB t ,with Xε 0 ∼ π 0 and where Ẑ(t, x) is an approximation of E[Z|X t = x]. We have the following result.Proposition C.1 (Optimality and stochastic interpolant): Denote P ε , respectively Pε the path measures associated with (X ε t ) t∈[0,1] and ( Xt ) t∈[0,1] respectively. Consider ℓ(ε) = KL(P ε | Pε ). Let ε ⋆ = argmin ε ℓ(ε). Then we have(ε ⋆ t ) 2 = 2γ t γt -2γ 2 t αt /α t .In particular, if α t = t, β t = 1 -t and γ t = σ 0 t(1 -t), then ε t = √ 2σ 0 . The value (ε ⋆ t ) t∈[0,1] corresponds to Markovian projection when it is well defined.Proof. We have that for anyKL(P ε | Pε ) = 1 0 1 ε 2 t ( γt -γ t αt /α t -ε 2 t /(2γ t )) 2 E[∆ t ]dt,where∆ t = ∥ Ẑ(t, X ε t ) -E[Z | X t = X ε t ]∥ 2and the expectation is w.r.t. P ε . We have thatKL(P ε | Pε ) = 0 if ε = ε ⋆ , which concludes the proof.Proposition C.1 is related to (Chen et al., 2024, Section 3.4 ). Therein it is noticed that, in the case of Augmented Bridge matching (De Bortoli et al., 2023) , the choice of ε t does not affect the joint distribution of (X ε 0 , X ε 1 ). The authors then select (ε t ) so as to minimise an approximation error. They show that, in that case, they recover the Föllmer process.We now show that Proposition C.1 can be further strengthened to establish that ε ⋆ is also the optimal value if we interpolate between π s and π 1 , or π 0 and π s , for any s ∈ [0, 1] and π s the distribution of X s . Consider in this context for any s, t ∈ [0, 1] with t ≥ s, γ t /γ s ≥ α t ≥ α s the following interpolation model.EQUATIONwhere X s ∼ π s , X 1 ∼ π 1 and Z ∼ N (0, Id). Assume that α t = 1-t, β t = t and γ t = σ 0 t(1 -t) for any t ∈ [0, 1] then (16) corresponds to the Brownian bridge associated with (σ 0 B t ) t∈[0,1] with endpoints X s at time s and X 1 at time 1. We have the following proposition.Proposition C.2 (Stochastic interpolant with intermediate time points): Define (X ε t,s ) t∈[s,1] given byEQUATION+ E ( γt,s -γ t,s αt /α t -ε 2 t,s /(2γ t,s ))Z | X t = X ε t,s + ε t,s dB t , X ε s,s ∼ π s , with γ t,s = γ 2 t -α 2 t γ 2 s /α 2 s .Then for any t ∈ [s, 1], X ε t,s and X t defined by (16) have the same distribution.Proof. We let s ∈ [0, 1] and X 1 , X s ∈ R d . From ( 16), we have directly that for any t ∈ [s, 1]dX t = [( αt /α s )X s + ( βt -αt β s /α s )X 1 + γt,s Z]dt,where Z ∼ N (0, Id). In addition, rearranging (16), we also have thatX s = (α s /α t )X t -(α s β t /α t -β s )X 1 -γ t,s (α s /α t )Z.Hence, by combining these two expressions, we get thatdX t = [( αt /α t )X t + ( βt -αt β t /α s )X 1 + ( γt,s -γ t,s ( αt /α t ))Z]dt.It follows that (X t,s ) t∈ [s,1] given bydX t,s = ( αt /α t )X t,s + ( βt -β t αt /α t )E[X 1 | X t = X t,s ] + ( γt,s -γ t,s αt /α t )E[Z | X t = X t,s ], X s,s ∼ π s ,is such that for any t ∈ [s, 1] the same distribution as X t defined by ( 16). Then, we conclude similarly to (Albergo et al., 2023, Theorem 2.8, Corollary 2.10) .We now consider the following approximate version of ( 17) d Xε t,s = ( αt /α t ) Xε t,s + ( βt -β t αt /α t )E X 1 | X t = Xε t,s + ( γt,s -γ t,s αt /α t -ε 2 t,s /(2γ t,s )) Ẑ(t, X ε t,s ) + ε t,s dB t , X ε s,s ∼ π s ,(ε ⋆ t ) 2 = 2γ t γt -2γ 2 t αt /α t .In particular, ε ⋆ does not depend on s ∈ [0, 1] and for everys 1 , s 2 ∈ [0, 1] with s 1 ≤ s 2 , we have that (X ε t,s1 ) t∈[s2,1] and (X ε t,s2 ) t∈[s2,1] coincide.Proof. Similarly to Proposition C.1, we get first that for any s, tEQUATIONSecond, we have that for any s, tEQUATIONThird, we have thatγ 2 t,s αt /α t = γ 2 t αt /α t -αt α t γ 2 s /α 2 s .(20) Combining ( 18), ( 19) and ( 20), we can conclude.In this section, we prove the main theoretical results of the paper. In Appendix D.1, we first prove the convergence of the α-IMF sequence, i.e. we prove Theorem 3.1. Second, we show that the non-parametric updates (8) correspond to the α-IMF sequence, i.e. we prove Proposition 3.2. In Appendix D.2, we link the non-parametric updates to the parametric updates.Let Q ∈ P(C) be associated with (√ εB t ) t∈[0,1], where (B t ) t∈[0,1] is a d-dimensional Brownian motion and ε > 0. In this section, we abuse notation and denote P(C) the set of path measures which are not necessarily probability path measures. In particular, we will consider Q ∈ P(C) associated with ( √ εB t ) t∈[0,1] with Q 0 = Leb. In that case, the Kullback-Leibler divergence is still well-defined and we refer to (Léonard, 2014) for more details. We recall that we have defined (P n , Pn ) n∈N for any n ∈ N and α ∈ (0, 1] byP n = proj M ( Pn ), Pn+1 = (1 -α) Pn + αproj R (proj M ( Pn )).In addition, for any n ∈ N, t ∈ [0, 1) and x ∈ R d we have definedv n+1 t (x) = v n t (x) -δ n ∇ µ n L t (v n t , P v n )(x), where L t (v t , P) = 1 2 (R d ) 3 v t (x t ) - x 1 -x t 1 -t 2 dP 0,1 (x 0 , x 1 )dQ t|0,1 (x t |x 0 , x 1 ) (21) = 1 2 (R d ) 3 v t (x t ) - x 1 -x t 1 -t 2 dproj R (P) 1,t .We define (P v n ) n∈N associated with ( 27), where for any suitable vector field v, P v is associated withdX t = v t (X t )dt + √ εdB t ,where(B t ) t∈[0,1] is a d-dimensional Brownian motion.In order to rigorously prove Proposition D.1 detailed further, we introduce P 2 (C), such that P ∈ P 2 (C) if P ∈ P(C) and for(R d ) 2 {∥x 0 ∥ 2 + ∥x 1 ∥ 2 }dP 0,1 (x 0 , x 1 ) < +∞.Note that if P ∈ P 2 (C) then we have that for any t ∈ [0, 1]R d ∥x t ∥ 2 dproj R (P) t < +∞.In addition, we recall that ϕ ∈ L 2 (µ)for µ ∈ P(R d ) if ϕ : R d → R d and R d ∥ϕ(x)∥ 2 dµ(x) < +∞.Finally, we define A 2 = {(ϕ, P) :P ∈ P 2 (C), ϕ ∈ L 2 (P)}.Then for any t ∈ [0, 1), we define L t : A 2 → R given for any (v, P) ∈ A 2 by ( 21).Proposition D.1 (Non-parametric updates are α-IMF): Let α ∈ (0, 1], (P n , Pn ) n∈N as in (4), δ n = α and µ n = (1 -α) Pn + αproj R (P n ). Assume that for any n ∈ N, P v n is well-defined.Then, for any n ∈ N, P v n = P n .Proof. First, we have that for any t ∈ [0, 1), v, P ∈ A 2 and ϕ ∈ L 2 (P t ) we haveL t (v t + εϕ, P) = L t (v t , P) + ε (R d ) 3 ⟨ϕ(x t ), v t (x t ) -x1-x1 1-t ⟩dP 0,1 (x 0 , x 1 )dQ t|0,1 (x t |x 0 , x 1 ) + (ε 2 /2) (R d ) 3 ∥ϕ(x t )∥ 2 dP 0,1 (x 0 , x 1 )dQ t|0,1 (x t |x 0 , x 1 ) = L t (v t , P) + ε (R d ) 2 ⟨ϕ(x t ), v t (x t ) - R d x 1 dproj R (P) 1|t (x 1 |x t ) -x t /(1 -t)⟩dproj R (P) t (x t ) + (ε 2 /2) R d ∥ϕ(x t )∥ 2 proj R (P) t (x t ).Hence, we have thatEQUATIONAssume that for some n ∈ N we have that for any t ∈ [0, 1) andx t ∈ R d , we have v k t (x t ) = E Pk [X 1 | X t = x t ] -x t /(1 -t).We are going to show that for any t ∈ [0, 1) andx t ∈ R d , we have v n+1 t (x t ) = E Pn+1 [X 1 | X t = x t ] -x t /(1 -t). For any t ∈ [0, 1) and x t ∈ R d , we denote δn t (x t ) = δ n (dproj R (P n ) t /dµ n t )(x t ).Since we have that δ n = α and µ n = (1 -α) Pn + αproj R (P n ), we obtain for any t ∈ [0, 1] andEQUATIONso thatEQUATIONTherefore, combining ( 8) with ( 23), ( 24), ( 22), we get that for any t ∈ [0, 1) andx t ∈ R d v n+1 t (x t ) = (1 -δn t (x t ))v n t (x t ) + δn t (x t ) E proj R (P n ) [X 1 | X t = x t ] -x t /(1 -t) = (1 -δn t (x t )) E Pn [X 1 | X t = x t ] -x t /(1 -t) + δn t (x t ) E proj R (P n ) [X 1 | X t = x t ] -x t /(1 -t) = (1 -δn t (x t )) R d x 1 d Pn 1|t (x 1 |x t ) -x t /(1 -t) + δn t (x t ) R d x 1 dproj R (P n ) 1|t (x 1 |x t ) -x t /(1 -t) = R d (x 1 -x t )/(1 -t)d[(1 -δn t (x t )) Pn 1|t + δn t (x t )proj R (P n ) 1|t ](x 1 |x t ) = R d x 1 d[(1 -α) Pn + αproj R (P n )] 1|t (x 1 |x t ).Hence, we have that for any t ∈ [0, 1) andx t ∈ R d , v n+1 t (x t ) = E Pn+1 [X 1 | X t = x t ] -x t /(1 - t). Since, for any t ∈ [0, 1) and x t ∈ R d , v 0 t (x t ) = E P0 [X 1 | X t = x t ] -x t (1 -t)by definition, we get that for any n ∈ N, t ∈ [0, 1) andx t ∈ R d , v n t (x t ) = E Pn [X 1 | X t = x t ] -x t /(1 -t).Using, Definition 2.2, we get that P v n = proj M ( Pn ), which concludes the proof.Before stating our convergence theorem, we show the following result which is a direct consequence of (Léonard, 2014, Theorem 2.12 ) and (Léonard et al., 2014, Theorem 2.14) . We recall that the differential entropy of a probability measure π is given byH(π) = - R d log((dπ/dLeb)(x))dπ(x),if π admits a density with respect to the Lebesgue measure and +∞ otherwise.Recall that Q is associated with ( √ εB t ) t∈[0,1] and assume that Q 0 = Leb. Let π 0 , π 1 ∈ P(R d ) such that R d ∥x∥ 2 dπ i (x) < +∞, H(π i ) < +∞,for i ∈ {0, 1}. Let P ⋆ such that P ⋆ is Markov, P ⋆ ∈ R(Q), P ⋆ 0 = π 0 and P ⋆ 1 = π 1 . Then P ⋆ is the Schrödinger Bridge, i.e. the unique solution to (2).Proof. First, we have that Q 0,1 is equivalent to Leb ⊗ Leb. Indeed, we have that for anyx 0 , x 1 ∈ R d (dQ 0,1 /d(Leb ⊗ Leb))(x 0 , x 1 ) = (2πε) -d/2 exp[-∥x 0 -x 1 ∥ 2 /(2ε)].Similarly, we have that for any t ∈ (0, 1) andx t ∈ R d , Q 0,1|t (•|x t ) is equivalent to Leb ⊗ Leb.Indeed, we have that for any t ∈ (0, 1) andx 0 , x t , x 1 ∈ R d (dQ 0,1|t (•|x t )/d(Leb ⊗ Leb))(x 0 , x 1 ) = (2πεt) -d/2 exp[-∥x 0 -x t ∥ 2 /(2εt)] × (2πε(1 -t)) -d/2 exp[-∥x t -x 1 ∥ 2 /(2ε(1 -t))].Hence, for any t ∈ (0, 1) andx t ∈ R d , Q 0,1|t (•|x t ) is equivalent to Q 0,1. Since P ⋆ is Markov and P ⋆ ∈ R(Q) we get that there exist φ • 0 and φ ⋆ 1 which are Lebesgue measurable such that for any ω ∈ C we have thatEQUATION)Second we verify that the conditions (i)-(vii) of (Léonard, 2014 , Theorem 2.12) are satisfied. First, Q is Markov and hence (i) is satisfied. Then, (ii) is satisfied since for any t ∈ (0, 1) andx t ∈ R d , Q 0,1|t (•|x t ) is equivalent to Q 0,1 . We have that Q 0 = Q 1 = Leb and (iii) is satisfied. We have that for any x 0 , x 1 ∈ R d (dQ 0,1 /d(Leb ⊗ Leb))(x 0 , x 1 ) = (2πε) -d/2 exp[-∥x 0 -x 1 ∥ 2 /(2ε)] ≥ (2πε) -d/2 exp[-∥x 0 ∥ 2 /ε -∥x 1 ∥ 2 /ε].Hence, (iv) is satisfied and we let A : R d → R + be given for any x ∈ R d by A(x) = ∥x∥ 2 /ε. In addition, we have that for anyx 0 , x 1 ∈ R d (R d ) 2 exp[-∥x 0 ∥ 2 /ε -∥x 1 ∥ 2 /ε]dQ 0,1 (x 0 , x 1 ) < +∞.Hence, (v) is satisfied and we let B : R d → R + given for any x ∈ R d by B(x) = ∥x∥ 2 /ε. By assumption (vi) and (vii) are satisfied. We conclude the proof upon using (Léonard, 2014 , Theorem 2.12-(b)) and ( 25).We are now ready to state our main convergence result. Proposition D.3 (Convergence of α-IMF): Let α ∈ (0,∈ N KL( Pn+1 |P ⋆ ) = KL((1 -α) Pn + αproj R (proj M ( Pn ))|P ⋆ ) ≤ (1 -α)KL( Pn |P ⋆ ) + αKL(proj R (proj M ( Pn ))|P ⋆ ) ≤ (1 -α)KL( Pn |P ⋆ ) + αKL(proj M ( Pn ) 0,1 |P ⋆ 0,1 ) ≤ (1 -α)KL( Pn |P ⋆ ) + αKL(proj M ( Pn )|P ⋆ ) ≤ (1 -α)KL( Pn |P ⋆ ) + αKL( Pn |P ⋆ ) -αKL( Pn |proj M ( Pn )).(26) Therefore, we get that αKL(Pn |proj M ( Pn )) ≤ KL( Pn |P ⋆ ) -KL( Pn+1 |P ⋆ ). Hence, it follows that n∈N KL( Pn |proj M ( Pn )) ≤ 2KL( P0 |P ⋆ ) < +∞.So we obtain lim n→+∞ KL( Pn |proj M ( Pn )) = 0. In addition, using (26) we have that KL( Pn |P ⋆ ) ≤ KL( P0 |P ⋆ ) < +∞ for all n ∈ N. Using (Shi et al., 2023, Lemma 6) , we also get that KL(proj M ( Pn )|P ⋆ ) ≤ KL( P0 |P ⋆ ) < +∞ for any n ∈ N. Hence both the sequences ( Pn ) n∈N and (P n ) n∈N = (proj M ( Pn )) n∈N are relatively compact in P(C). Let P ∈ P(C) be an adherent point to the sequence ( Pn ) n∈N and φ : N → N increasing such that lim n→+∞ P φ(n) = P. Similarly, let ϕ : N → N increasing such that (ϕ(n)) n∈N is a subsequence of (φ(n)) n∈N such that lim n→+∞ proj M ( Pn ) = P′ , with P′ and adherent point to the sequence (proj M ( Pn )) n∈N . Using the lower semi-continuity of the Kullback-Leibler divergence in both arguments (Dupuis and Ellis, 2011) , we get thatKL( P| P′ ) ≤ lim inf n→+∞ KL( Pϕ(n) |proj M ( Pϕ(n) )) = 0.Since the set of Markov measures and the set of reciprocal measures w.r.t. Q are both closed, we have that P is Markov and in the reciprocal class of Q. Since we also have that P0 = π 0 and P1 = π 1 , we get that P = P ⋆ using Appendix D.1. Since every adherent point of ( Pn ) n∈N is P ⋆ , we have that lim n→+∞ Pn = P ⋆ . Similarly, using that lim n→+∞ KL( Pn |proj M ( Pn )) = 0 and again the lower semi-continuity of the Kullback-Leibler divergence in both arguments, we get that every adherent point of (proj M ( Pn )) n∈N is P ⋆ . Hence, we have that lim n→+∞ P n = P ⋆ , which concludes the proof.In this section, we show that the parametric updates considered in (9) are a preconditioned version of the non-parametric updates considered in (8). We first recall the non-parametric lossL(v, P) = 1 0 L t (v t , P)dt = 1 2 1 0 (R d ) 2 v t (x t ) - x 1 -x t 1 -t 2 dproj R (P) t,1 (x t , x 1 )dtand the parametric lossL(θ, P) = 1 2 1 0 R d ×R d v θ t (x t ) - x 1 -x t 1 -t 2 dproj R(Q) (P) t,1 (x t , x 1 )dt.The non-parametric sequence (v n ) n∈N is given by ( 27), i.e. we have for any n ∈ N, t ∈ [0, 1] andEQUATION)Similarly the sequence of parametric updates is given for any n ∈ N, t ∈ [0, 1] and x ∈ R d byθ n+1 = θ n -α∇ θ L(θ n , P θn ).We recall that P θn is a stop gradient version of P v θn . We are going to show that on average the parametric algorithm yields a direction of descent for the non-parametric loss. We assume that the set of parameters Θ is an open subset of R p for some p ∈ N. For any t ∈ [0, 1] and x ∈ R d we assume that θ → v θ t (x) is twice continuously differentiable and denote D θ v θ t (x) ∈ R d×p its Jacobian and D 2 θ v θ t (x) its Hessian. For any θ ∈ Θ, we denoteh θ = ∇ θ L(θ, P θ ).We show the following result.Proposition D.4 (Velocity field parametric update): Assume that there exists C > 0 such that for any θ ∈ Θ andEQUATIONwhere h θ = ∇ θ L(θ, P θ ). We have that for any n ∈ N, t ∈ [0, 1) andx ∈ R d v θn+1 t (x) = v θn t (x) -αD θ v θn t (x) 1 0 R d D θ v θ s (x) ⊤ ∇ µ n L s (v θn , P v θn )(x)dµ n s (x)ds + o(α),where µ n = (1 -α)P n + αP v θn .Proof. First, we have that for any µ ∈ P(C)∇ θ L(θ, P θ ) = 1 0 R d D θ v θ s (x s ) ⊤ (v θ s (x s ) -(x 1 -x s )/(1 -s))dproj R (P θ ) s,1 (x s , x 1 )ds.Therefore, using ( 22), we get that∇ θ L(θ, P θ ) = 1 0 R d D θ v θ s (x s ) ⊤ ∇ µ L s (v θ s , P θ )(x s )dµ s (x s )ds.Let θ ∈ Θ and denote θ ′ = θ -α∇ θ L(θ, P θ ). Using a Taylor expansion, we get that for any θ ∈ Θ, we have thatv θ ′ t (x) = v θ t (x) -αD θ v θ t (x) 1 0 R d D θ v θ s (x s ) ⊤ ∇ µ L s (v θ s , P θ )(x s )dµ s (x s )ds + α 2 1 0 (1 -s)D 2 θ v θ-αsh θ s (x)(h θ , h θ )ds.Since the functional gradient is not applied on the second coordinate, we can drop the stop gradient operator and therefore we have for any θ ∈ Θv θ ′ t (x) = v θ t (x) -αD θ v θ t (x) 1 0 R d D θ v θ s (x s ) ⊤ ∇ µ L s (v θ s , P v θ )(x s )dµ s (x s )ds + α 2 1 0 (1 -s)D 2 θ v θ-αsh θ s (x)(h θ , h θ )ds.Combining this result with (28), we conclude the proof.The corresponding update on the velocity field is given for any n ∈ N, t ∈ [0, 1] and x ∈ R d byd n t (x) = -αD θ v θn t (x) 1 0 R d D θ v θ s (x) ⊤ ∇ µ n L s (v θn , P v θn )(x)dµ n s (x) + o(α).We immediately have the following corollary.Proposition D.5 (Parametric direction of descent): For any n ∈ N, if1 0 R d D θ v θ s (x) ⊤ ∇ µ n L s (v θn , P v θn )(x)dµ n s (x) ̸ = 0, then we have lim α→0 1 0 R d ⟨∇ µ n L t (v θn , P v θn )(x), d n t (x)⟩dµ n t (x) ≤ 0.In this section we recall some basics on Markovian and reciprocal projections in Appendix E.1. We explain the link between the concept of iterative refinement and Schrödinger Bridges in Appendix E.2. Then, we briefly present Diffusion Schrödinger Bridge Matching (DSBM) (Shi et al., 2023) in Appendix E.3 and propose some new extensions in Appendix E.4.In this section, we recall the definition of the reciprocal and Markov projection. We provide more details on how these different projections can be performed and illustrate them on simple examples.Markov projection. First, we recall the definition of the Markovian projection.Definition E.1 (Markov projection): Assume that Q is induced by ( √ εB t ) t∈[0,1] for ε > 0.Then, when it is well-defined, for any P ∈ R(Q), the Markovian projection M = proj M (P) ∈ M is the path measure induced by the diffusiondX ⋆ t = v ⋆ t (X ⋆ t )dt + √ εdB t , v ⋆ t (x t ) = E P 1|t [X 1 | X t = x t ] -x t /(1 -t), X ⋆ 0 ∼ P 0 .In Figure 7 and Figure 8 , we illustrate the effect of the Markovian projection, following the example of (Liu, 2022) . We consider two distributions π 0 and π 1 such thatπ 0 = 1 2 N ([-2, -2], Id) + 1 2 N ([-2, 2], Id), π 1 = 1 2 N ([2, -2], Id) + 1 2 N ([2, 2], Id).In Figure 7 , we display samples from the distributions π 0 and π 1 as well as trajectories from the path measure P = (π 0 ⊗ π 1 )Q |0,1 . Practically, this means that we sample X 0 ∼ π 0 and X 1 ∼ π 1 independently and then consider a Brownian bridge between X 0 and X 1 . The SDE associated with the Brownian bridge with scale ε > 0 is given for any t ∈ [0, 1] byEQUATIONNote that the measure P = (π 0 ⊗ π 1 )Q |0,1 is in the reciprocal class, i.e. P ∈ R(Q). Next in Figure 8 , we display samples from the distributions π 0 and π 1 as well as trajectories from the path measure P ⋆ = proj M (P). Note that in Figure 8 , contrary to Figure 7 , we observe less crossings between the trajectories. Indeed in the limit case where ε → 0 the Markov measures P ⋆ is an ODE with regular coefficients and therefore admits a unique solution for every starting point in the space so no crossing is possible. In particular, note that most of the trajectories starting from the upper-left Gaussian end at the upper-right Gaussian. Similarly, most of the trajectories starting from the lower-left Gaussian end at the lower-right Gaussian.o Figure 8 : Samples from the original distributions are shown in red, while sample paths from M = proj M (P) are shown in blue.In practice, computing the Markov projection involves finding the optimal drift v ⋆ t . This optimal drift is the minimizer of a regression problem, see (Shi et al., 2023) for more details. Hence, computing the Markovian projection requires training a neural network to define a vector field.Reciprocal projection. First, we recall the definition the reciprocal projection.∈ P(C) is in the reciprocal class R(Q) of Q if P = P 0,1 Q |0,1 .We define the reciprocal projection of P ∈ P(C) as P ⋆ = proj R(Q) (P) = P 0,1 Q |0,1 . We will write proj R instead of proj R(Q) to simplify notation.To sample from P ⋆ = proj R (P), we only need to sample (X 0 , X 1 ) ∼ P 0,1 and then to sample from the Brownian bridge conditioned on (X 0 , X 1 ). This means that in order to sample X ⋆ t ∼ P ⋆ t , we only need to sample (X 0 , X 1 ) ∼ P 0,1 and then computeEQUATIONwith Z ∼ N (0, Id). In particular, sampling from P ⋆ = proj R (P) does not require training any neural network. However, in practice, in order to obtain samples (X 0 , X 1 ) ∼ P, we have that P is associated with an SDE and therefore obtaining (X 0 , X 1 ) requires unrolling the SDE associated with P. In Algorithm 1, the measure P is associated with an SDE with parametric drift v θ .In Figure 9 , we continue our study of the example of (Liu, 2022) that we used to explain the concept of Markovian projection. We consider the path measure M obtained as the Markov projection of P = (π 0 ⊗ π 1 )Q |0,1 . In Figure 9 , we display samples from the distributions π 0 and π 1 as well as trajectories from the path measure P ⋆ = M 0,1 Q |0,1 . In order to sample from P ⋆ we first sample (X 0 , X 1 ) ∼ M 0,1 . This involves unrolling the SDE associated with M. Once we have access to samples (X 0 , X 1 ), we draw trajectories from the Brownian bridge following the SDE (29). We can also sample from any time t without having to unroll the SDE (29) by simply sampling from (30). This is what is done in Algorithm 1. (Peluchetti, 2023; Shi et al., 2023) . We highlight that in the case where ε → 0 then DSBM is equivalent to the Rectified Flow algorithm (Liu et al., 2023b) . One of the main limitation of those previously introduced procedures which provably converge to the solution of the Schrödinger Bridge problem is that they rely on these expensive iterative solvers and requires to consider two networks, one parameterising the forward process π 0 → π 1 and one parameterising the backward π 1 → π 0 .Diffusion Schrödinger Bridge Matching corresponds to the practical implementation of the Iterative Markovian Fitting procedure proposed in Shi et al. (2023) ; Peluchetti (2023) . The IMF procedure alternates between projecting on the Markov class M and the reciprocal class R Q . In what follows, we denote M n+1 = P 2n+1 ∈ M and Π n = P 2n ∈ R(Q). We also recall that Q is a (rescaled) Brownian motion associated with (σ 0 B t ) t∈[0,1] and that therefore sampling from Q |0,1 (•|x 0 , x 1 ) corresponds to sampling fromdX t = (x 1 -X t )/(1 -t)dt + σ 0 dB t , X 0 = x 0 .We recall that the main computational bottleneck of the DSBM lies in the approximation of the Markovian projection. Indeed, using (Shi et al., 2023 , Definition 1, Proposition 2), we have thatM ⋆ = proj M (Π) is associated with the process dX t = (E Π [X 1 | X t ] -X t )/(1 -t)dt + σ 0 dB t X 0 ∼ π 0 .We also have using (Shi et al., 2023, Proposition 2 ) that M ⋆ can be approximated using M θ ⋆ given byEQUATIONEQUATIONwhere {v θ : θ ∈ Θ} is a parametric family of functions, usually given by a neural network.Hence, since we can approximate proj R(Q) (M) and proj M (Π) we can approximate the IMF procedure. This is the DSBM algorithm introduced in Shi et al. ( 2023); Peluchetti (2023) . We describe the first few iterations. LetΠ 0 = Π 0 0,1 Q |0,1 where Π 0 0 = π 0 , Π 0 1 = π 1 . Learn M 1 ≈ proj M (Π 0 )given by ( 31) with v θ ⋆ given by ( 32). Next, sample fromΠ 1 = proj R(Q) (M 1 ) = M 1 0,1 Q |0,1by sampling from M 1 0,1 and reconstructing the bridge Q |0,1 . Upon iterating the previous procedure, we obtain a sequence (Π n , M n+1 ) n∈N . To mitigate the bias accumulation problem caused by approximating only the forward process, we alternate between a forward Markovian projection and a backward Markovian projection. We give more details on the advantage of using a forward-backward parameterisation instead of a forward-forward in Appendix I. This procedure is valid using (Shi et al., 2023, Proposition 9) . The optimal backward process is approximated withEQUATIONEQUATIONWe recall the full DSBM algorithm in Algorithm 2.Algorithm 2 Diffusion Schrödinger Bridge Matching1: Input: Joint distribution Π 0 0,1 , tractable bridge Q |0,1 , number of outer iterations N ∈ N. 2: Let Π 0 = Π 0 0,1 Q |0,1 . 3: for n ∈ {0, . . . , N -1} do 4:Learn v ϕ ⋆ using (34) with Π = Π 2n .Let M 2n+1 be given by (33).Let Π 2n+1 = M 2n+1 0,1 Q |0,1 .Learn v θ ⋆ using ( 5) with Π = Π 2n+1 .8:Let M 2n+2 be given by (31). 9:Let Π 2n+2 = M 2n+2 0,1 Q |0,1 . 10: end for 11: Output: v θ ⋆ , v ϕ ⋆First, we consider a reflection-projection method similar to the one investigated in Bauschke and Kruk (2004) . We recall that the DSBM algorithm is associated with a sequence (P n ) n∈N such that for any n ∈ N P n+1/2 = proj M (P n ),P n+1 = proj R(Q) (P n+1/2). In a reflection-projection scheme, one of the projection is replaced by a reflection. As noted in Bauschke and Kruk (2004) , this can yield faster convergence rates in practice. We consider the sequence (P n ) n∈N such that for any n ∈ N P n+1/2 = proj M (P n ),EQUATION)In what follows, we make the assumption that 2P n+1/2 -P n is a probability measure, even though it is not clear if this path measure is non-negative. However, even by making this strong assumption, we show that we can recover DSBM in Algorithm 4. By considering a relaxation of the reflectionprojection scheme.First, note that for any n ∈ N, P n |0 is Markov, see De Bortoli et al. (2023) for instance. Hence, we assume that P n is associated withEQUATION)Estimating P n+1 . First, we compute v n+1 t,0 assuming that we can sample from P n and P n+1/2 . Since P n+1 is in the reciprocal class, we have that P n+1 is associated withdX t = (E P n+1 1|0,t [X 1 | X t , X 0 ] -X t )/(1 -t)dt + σ 0 dB t .We refer to De Bortoli et al. ( 2023) for a proof of this fact. Hence, using (35), we have thatEQUATIONNext, we turn to the estimation of P n+3/2 .Estimating P n+3/2 . Next, we assume that for any n ∈ N, P n+1/2 is associated withdX n t = v n t (X t )dt + σ 0 dB t , X 0 ∼ π 0 .Using (36), we have that v n+1 t is given byv n+1 t = argmin v 1 0 R d ×R d ∥v t (t, x t ) -v n+1 t,0 (t, x t , x 0 )∥ 2 dP n+1 0,t (x 0 , x t ).We note also that using (Shi et al., 2023 , Proposition 2) and the fact P n is in the reciprocal class of Q, we also have thatEQUATION)Hence, assuming that we can sample from P n and P n+1/2 then we can estimate v n+1 t,0 and v n+1 t , i.e. sample from P n+1 and P n+3/2 . Note that the losses ( 37) and ( 38) only differ by the conditioning with respect to the initial condition x 0 and therefore the optimisation can be conducted in parallel.We are now able to propose the following projection-reflection algorithm, see Algorithm 3.1: Input: Vector field and conditional vector field v 0 t and v 0 t,0 , noise level σ 0 and associated bridge Q |0,1 , number of outer iterations N ∈ N, batch size B 2: for n ∈ {0, . . . , N -1} do 3:while not converged do4: Sample X 1:B 0 ∼ π ⊗B 0 5: Sample X 1:B 1 using dX 1:B t = v n t (X 1:B t )dt + σ 0 dB t 6: Sample X1:B 1 using d X1:B t = v n t,0 ( X1:B t , X 1:B 0 )dt + σ 0 dB t 7: L = 1 0 [ B k=1 ∥v t (X k t ) - X k 1 -X k t 1-t ∥ 2 -(1/2) B k=1 ∥v t (X k t ) - X k 1 -X k t 1-t ∥ 2 ]dt 8: L 0 = 1 0 [ B k=1 ∥v t,0 (X k t , X k 0 ) - X k 1 -X k t 1-t ∥ 2 -(1/2) B k=1 ∥v t,0 (X k t , X k 0 ) - X k 1 -X k t 1-t ∥ 2 ]dt 9: v n+1 t = Gradientstep(L) 10: v n+1 t,0 = Gradientstep(L 0 ) 11:end while 12: end for 13: Output:v N +1 t , v N +1 t,0Note that in Algorithm 3 we only consider the optimisation of a forward process but similarly to Algorithm 2, one can construct a forward backward extension to alleviate some of the bias accumulation problems. Finally, we can interpolate between DSBM and this new reflection algorithm and DSBM by introducing an hyperparameter α ≥ 0 and consider the following extension given in Algorithm 4Algorithm 4 Reflection Diffusion Schrödinger Bridge Matching 1: Input: Vector field and conditional vector field v 0 t and v 0 t,0 , noise level σ 0 and associated bridge Q |0,1 , number of outer iterations N ∈ N, batch size B 2: for n ∈ {0, . . . , N -1} do 3:while not converged do4: Sample X 1:B 0 ∼ π ⊗B 0 5: Sample X 1:B 1 using dX 1:B t = v n t (X 1:B t )dt + σ 0 dB t 6: Sample X1:B 1 using d X1:B t = v n t,0 ( X1:B t , X 1:B 0 )dt + σ 0 dB t 7: L = 1 0 [ B k=1 ∥v t (X k t ) - X k 1 -X k t 1-t ∥ 2 -α B k=1 ∥v t (X k t ) - X k 1 -X k t 1-t ∥ 2 ]dt 8: L 0 = 1 0 [ B k=1 ∥v t,0 (X k t , X k 0 ) - X k 1 -X k t 1-t ∥ 2 -α B k=1 ∥v t,0 (X k t , X k 0 ) - X k 1 -X k t 1-t ∥ 2 ]dt 9: v n+1 t = Gradientstep(L) 10: v n+1 t,0 = Gradientstep(L 0 ) 11:end while 12: end for 13: Output:v N +1 t , v N +1 t,0Using different values of α ≥ 0 in Algorithm 4, we recover different existing algorithms. If α = 1, we recover DSBM Shi et al. (2023) . Finally, if α = 1/2, we recover the reflection algorithm Algorithm 3.The idea of training both the forward and the backward jointly was mentioned in (Shi et al., 2023, Section G) . However, it was still assumed that, while being trained jointly, the forward and backward vector fields were obtained using an argmin operation, see (Shi et al., 2023, Equation (43) , ( 44)). In addition, in (Shi et al., 2023 , Section G) a consistency loss was proposed in order to enforce that the forward and backward processes match, see (Shi et al., 2023, Equation (49) ). In this section, we leverage new results from Daras et al. (2023a); De Bortoli et al. (2024) in order to enforce the internal consistency of the model. First, note that for any (X t ) t∈[0,1] associated with P ∈ R(Q) we have for any0 ≤ t 0 ≤ t ≤ t 1 ≤ 1 that X t = t -t 0 t 1 -t 0 X t1 + t 1 -t t 1 -t 0 X t0 + σ t0,t,t1 Z, Z ∼ N (0, Id),whereσ t0,t,t1 = (t -t 0 )(t 1 -t) t 1 -t 0 .Let p t be the density of X t with respect to the Lebesgue measure, we have that for any0 ≤ t 0 ≤ t ≤ t 1 ≤ 1 and x t ∈ R d p t (x t ) = R d ×R d (2πσ 2 t0,t,t1 ) -d/2 exp - ∥x t -t-t0 t1-t0 x t1 -t1-t t1-t0 x t0 ∥ 2 2σ 2 t0,t,t1 p t0 (x t0 )p t1 (x t1 )dx t0 dx t1 .Using the change of variable x t0 → x t0 + x t and x t1 → x t1 + x t we get that for anyEQUATIONThis identity for the score has already been presented in a bridge matching context in (De Bortoli et al., 2024, Section 3.3) . Let P ∈ R(Q) then we have that proj M (P) is such that for any t ∈ [0, 1], proj M (P) t = P t , see (Shi et al., 2023, Proposition 2) . We have that for any t ∈ [0, 1],v t (x) + v 1-t (x) = σ 2 0 ∇ log p t (x).Combining this result and (39), this suggests considering the following consistency lossEQUATION-v θ (t 0 , 1, X t ) -v θ (1 -t 0 , 0, X t0 ) -v θ (t 1 , 1, X t1 ) -v θ (1 -t 1 , 0, X t1 )∥ 2 ].Similarly to (13), we can consider an empirical version of (40).In Algorithm 1, the finetuning stage requires a pretrained bridge matching model interpolating between π 0 and π 1 (lines 2-7). However, for large datasets with complex distributions π 0 and π 1 , e.g. ImageNet, training this bridge model from scratch can be computationally expensive. To improve efficiency, we can leverage existing diffusion models targeting π 0 and π 1 . Specifically, we assume access to generative models transferring between N (0, Id) and π 0 , and between N (0, Id) and π 1 . In the rest of this section, we show how one can adapt Algorithm 1 to this setting. We then comment on the link between the proposed algorithm and Dual Diffusion Implicit Bridges (Su et al., 2023) .Setting. For simplicity, assume that we have two pretrained diffusion models for π 0 and π 1 . We describe our procedure for π 0 . Consider a forward process of the form X t = X 0 + σ t Z, with Z ∼ N (0, Id), where σ t is a hyperparameter. Note that we could have considered an interpolant of the form X t = α t X 0 + σ t Z instead, see Song et al. (2021a) for instance.We assume that the model X t = X 0 + σ t Z is associated with the forward diffusion modelEQUATIONwhere we assume that g t ≥ 0 for all t ∈ [0, 1]. Note that we have that for any t ∈ [0, 1], σ 2 t = t 0 g 2 s ds. In particular, we have that for any s, t ∈ [0, 1] with s ≤ tX t = X s + σ 2 t -σ 2 s Z, Z ∼ N (0, Id).Our goal is to solve the following Entropic Optimal Transport problemEQUATIONwhere ε > 0 is some entropic regularisation. We assume that ε > 0 is fixed and assume that there exists t ′ ∈ [0, 1] such that σ 2 t ′ = ε/2. We now consider a dynamic version of (42) withEQUATION)where Q is associated with (X t ) t∈[0,t ′ ] (41). Note that contrary to the setting presented in the main paper, here we do not consider the integration between time 0 and 1 but between time 0 and t ′ . It can be shown that for any t ∈ [0, t ′ ], (X t ) t∈[0,t ′ ] associated with Q t|0,t ′ is given byX t = Interp t (X 0 , X t ′ , Z) = 1 - σ 2 t σ 2 t ′ X 0 + σ 2 t σ 2 t ′ X t ′ + σ t 1 - σ 2 t σ 2 t ′ Z, Z ∼ N (0, Id).Solving ( 43) is equivalent to solving (42). We now propose an algorithm to solve (43). It corresponds to the finetuning stage of Algorithm 1 with a specific initialisation, similar to DSBM-IPF in Shi et al. (2023) .By v ϕ , we denote a DDM model associated with π 1 :EQUATIONSimilarly, v θ denotes a diffusion model associated with π 0 :EQUATIONIn analogy to Equation ( 11), the two equations above correspond to the forward and backward SDEs.For a given batch of inputs X 1:B 0 and X 1:B 1 , timesteps t ∼ Unif([0, t ′ ]) ⊗B , and interpolations X t and X t , we compute the empirical forward and backward losses as the following modification of Equation ( 13):ℓ (ϕ; t, X 1 , X t ) = 1 B B i=1 v ϕ t i , X i t -X i 1 -X i t /σ i t 2 , ℓ (θ; t, X 0 , X t ) = 1 B B i=1 ∥v θ t i , X i t -X i 0 -X i t / σ 2 t ′ -σ 2 t i )∥ 2 .Algorithm 5 corresponds to an online version of DSBM-IPF (Shi et al., 2023) with the initialisation given by two generative models. In Algorithm 5, we finetune the trained vector fields to solve the interpolation task. At inference time, the SDE associated with vector field v θ interpolates between π 1 → π 0 , while the SDE associated with the vector field v ϕ interpolates between π 0 → π 1 .Algorithm 5 α-Diffusion Schrödinger Bridge Matching for DDM finetuningWe highlight links between our proposed flow and Sinkhorn flows in Appendix H.1. We draw connection between our practical approach and Reinforcement Learning in Appendix H.2. We discuss how α-IMF is related to (incremental) Expectation-Maximisation in Appendix H.3. Finally, we discuss how our algorithm can be seen as an instance of continual learning in Appendix H.5.In this section, we discuss the links between our approach and the Sinkhorn flow introduced by Karimi et al. (2024) . We start by recalling how Sinkhorn flows are defined and then discuss how they are related to our approach.γ-Sinkhorn and Sinkhorn flows. We first consider the static EOT problem and recall the Sinkhorn procedure, also called Iterative Proportional Fitting. We define a sequence of coupling ( Πn , Π n ) n∈N , i.e. for any n ∈ N,Π n ∈ P(R d × R d ).We let Π 0 = Q 0,1 and we consider for any n ∈ N,Π n = argmin{KL(Π | Πn ) : Π ∈ P(R d × R d ), Π 0 = π 0 }, (46) Πn+1 = argmin{KL(Π | Π n ) : Π ∈ P(R d × R d ), Π 1 = π 1 },In Karimi et al. (2024) , the authors generalise ( 46) by introducing an extra hyperparameter γ ∈ (0, 1] and definingΠ n = argmin{KL(Π | Πn ) : Π ∈ P(R d × R d ), Π 1 = π 1 }, (47) Πn+1 = argmin{γKL(Π | Π n ) + (1 -γ)KL(Π | Πn ) : Π ∈ P(R d × R d ), Π 0 = π 0 },Using (Karimi et al., 2024 , Lemma 2), we have that for any γ ∈ (0, 1], any n ∈ N and anyEQUATIONwith f 0 γ = g 0 γ = 0 and for any n ∈ N, γ ∈ (0, 1] andEQUATIONIn addition, using (Karimi et al., 2024 , Equation ( 9)) we have that for any n ∈ N, γ ∈ (0, 1] andx 0 ∈ R d f n γ (x 0 ) = -log R d exp[g n γ (x 1 ) -(1/(2ε))∥x 0 -x 1 ∥ 2 ]dπ 1 (x 1 ) .When letting γ → 0, (48) and ( 49) suggest to consider for any s ≥ 0,x 0 , x 1 ∈ R d (d Πs /dQ 0,1 )(x 0 , x 1 ) = exp[f s (x 0 ) + g s (x 1 )], Π s = argmin{KL(Π | Πs ) : Π, Π 1 = π 1 },where for any s ≥ 0,x 1 ∈ R d ∂ s g s (x 1 ) = -log(d Πs 1 /dπ 1 )(x 1 ), ∂ s f s (x 0 ) = R d log(d Πs 1 /dπ 1 )(x 1 )d Πs (x 1 |x 0 ).Comparison with Schrödinger Bridge flows. In order to compare our approach with the one of Karimi et al. (2024) , we start by rewriting the γ-Sinkhorn algorithm defined by (47). To do so, we introduce the projection on the measures with fixed marginal.Definition H.1 (Projection on marginals): Let Π ∈ P(R d × R d ) and π 0 ∈ P(R d ), we define proj 0,π0 (Π) as followsproj 0,π0 (Π) = argmin{KL( Π | Π) : Π ∈ P(R d × R d ), Π0 = π 0 }.Similarly, for any Π ∈ P(R d × R d ) and π 1 ∈ P(R d ), we define proj 1,π1 (Π) as follows 2 : Comparison between γ-Sinkhorn and γ-IMF. N ∈ N is the size of the buffer equipped with two functions Add and Sample. We have that Addproj 1,π1 (Π) = argmin{KL( Π | Π) : Π ∈ P(R d × R d ), Π1 = π 1 }. γ-Sinkhorn γ-IMF Loss function KL(Π |proj 1,π1 ( Π)) 1 0 E proj R(Q) (P) [∥v t (X t ) -X1-Xt 1-t ∥ 2 ]dt regularisation KL(Π | Π) 1 0 R d ∥f t (x t ) -ft (x t )∥ 2 dµ t (x t )dt Update Implicit Explicit Table: Ω × k∈N (R 2d ) k × (R 2d ) N → (R 2d ) N ,where Ω is a probability space. In practice Add takes a random number (the function can be stochastic), any number of proposed samples as well as the current buffer. As an output Add returns the updated buffer. We also define Sample : Ω × N × (R 2d ) N → k∈N (R 2d ) k . This function takes a random number (the function can be stochastic), a natural number k representing the number of samples to return as well as the current buffer. As an output Sample returns a batch of k samples from the buffer.Algorithm 6 Replay Buffer Diffusion Schrödinger Bridge Matching 1: Input: π 0 , π 1 , ε (entropic regularisation), N pretraining (number of pretraining steps), N finetuning (number of finetuning steps), B (batch size), γ (EMA parameter), θ (initial parameters), B fwd (forward buffer), B bwd (backward buffer) 2: θ = θ 3: for n ∈ {0, . . . , N pretraining } do4: Sample (X 1:B 0 , X 1:B 1 ) ∼ (π 0 ⊗ π 1 ) ⊗B , t ∼ Unif([0, 1]), Z 1:B ∼ N (0, Id) ⊗B 5:Definition H.2 (Projections and maximisations): Let A be a subset of P(C). Then, for any P ∈ P(C), when it is well-defined, we define its E-projection on A as P ⋆ = argmin Q∈A KL(Q | P).Similarly, for any P ∈ P(C), when it is well-defined, we define its M-projection on A as P ⋆ = argmin Q∈A KL(P | Q).In Brekelmans and Neklyudov (2023) , the authors choose M-projection because this corresponds to the Maximisation step in an EM algorithm while the E-projection corresponds to the expectation step in the EM algorithm. In Brekelmans and Neklyudov (2023) , the authors highlight that the Iterative Proportional Fitting procedure is a Expectation-Expectation procedure, i.e. the alternating projections are both E-projections. In contrast, the Iterative Markovian Fitting procedure is a Maximisation-Maximisation procedure, i.e. the alternating projections are both M-projections. In particular, we can define the following sequence of path measures (P n ) n∈N , where for any n ∈ N we haveP n+1/2 = argmin P∈proj M KL(P n |P), P n+1 = argmin P∈R(Q) KL(P n+1/2 |P).In addition, we have thatEQUATIONH.4 Links with finetuning of diffusion models Algorithm 1 can be seen as a method to finetune bridge matching. Finetuning of diffusion models and flow matching procedures is an active research area. Most of the existing methodologies optimise for an external cost after a pretraining phase. These procedures rely on Reinforcement Learning strategies (Lee et al., 2023; Black et al., 2023; Fan et al., 2024) . Recently Direct Preference optimisation (DPO) (Rafailov et al., 2024) has been applied to the finetuning of diffusion models in (Yang et al., 2023; Rafailov et al., 2024) . Our approach departs from these works as the objective we minimise is given by the EOT cost. However all of these approaches involve some level of self-play, i.e. are not simulation free. (2023) . In (Masip et al., 2023) , the authors consider a weighted loss between a diffusion model loss and a distillation loss which ensures some consistency between the model being trained and the previous task model. Similarly to our approach this distillation loss is not simulation-free but, contrary to our loss, the clean samples are not obtained by unrolling the diffusion model but by applying a one-step prediction operator. In (Zaj ąc et al., 2023) , consider different replay buffer techniques to train continual diffusion models and observe that experience replay with a small coefficient can bring improvements. Finally, in (Smith et al., 2023) , the authors consider the continual training of a text-to-image diffusion model with LoRA (Hu et al., 2021) .I Forward-Forward, Forward-Backward and accumulation of errorIn this section, we investigate how error accumulates in the context of DSBM. In practice, we observe similar conclusions in the case of the online version of DSBM. We compare two methods: one which only trains a forward model and one which trains a forward and a backward model.In what follows, we assume that π 0 = π 1 = N (0, Id), we also assume that Q is associated with 1] . We recall that for any t ∈ [0, 1], we have that( √ 2B t ) t∈[0,X t = (1 -t)X 0 + tX 1 + 2t(1 -t)Z, Z ∼ N (0, Id).We are going to consider to approximate schemes to implement IMF.Forward-forward. First, we consider the following sequence of path measures (P n ) n∈N . We set P 0 = (π 0 ⊗ π 1 )Q |0,1 . For any n ∈ N, we define P 2n+2 = P 2n+1 0,1 Q |0,1 , i.e. P 2n+2 = proj R (P 2n+1 ). In addition, we define P 2n+1 = proj ε, M (P 2n ) such that P 2n+1 is associated with (X t ) t∈[0,1] where for any t ∈ [0, 1]EQUATIONwith ε ∈ R. Recall that if we define P2n+1 = proj M (P 2n ) we have that for any t ∈ [0, 1], P2n+1 is associated with (X t ) t∈ [0, 1] where for any t ∈ [0, 1]dX t = (E P 2n 1|t [X 1 | X t ] -X t )/(1 -t)dt + √ 2dB t .Hence, proj ε, M corresponds to making an error of order x → εx on the estimated velocity field. Doing so, we now longer have that for any n ∈ N, P n 1 = π 1 . In what follows, we are going to show how the error accumulates for the sequence P n 0,1 . Before stating Proposition I.1, we introduce f : R 4 → R such that for any c 0,0 , c 1,1 , c 0,1 > 0 and t ∈ [0, 1]f (c 0,0 , c 1,1 , c 0,1 , t) = [-(1 -t)c 0,0 + tc 1,1 + (1 -2t)c 0,1 -2t] /[(1 -t) 2 c 0,0 + t 2 c 1,1 + 2t(1 -t)c 0,1 + 2t(1 -t)].We define F (c 0,0 , c 1,1 , c 0,1 , ε, t) = 2 t 0 f (c 0,0 , c 1,1 , c 0,1 , s)ds + 2εt. Finally, we define f cov (c 0,0 , c 1,1 , c 0,1 , ε) = exp[ 1 2 F (c 0,0 , c 1,1 , c 0,1 , ε, 1)], as well asf var (c 0,0 , c 1,1 , c 0,1 , ε) = exp[ 1 2 F (c 0,0 , c 1,1 , c 0,1 , ε, 1)](1 + 2 1 0exp[-F (c 0,0 , c 1,1 , c 0,1 , ε, s)]ds).Proposition I.1 (Forward-Forward updates): For any n ∈ N, we have that P 2n+1 0,1 = N (0, Σ n+1 Id) whereΣ n+1 = Id c n+1 0,1 Id c n+1 0,1 Id c n+1 1,1 Id, and for any n ∈ Nc n+1 1,1 = f var (1, c n 1,1 , c n 0,1 , ε), c n+1 0,1 = f cov (1, c n 1,1 , c n 0,1 , ε).Proof. Let P = (P 0,1 )Q |0,1 where P 0,1 is a Gaussian random variable with zero mean and covariance matrix Σ ∈ R 2d×2d such thatΣ = Id c 0,1 Id c 0,1 Id c 1,1 Id ,where Id is the d-dimensional identity matrix and we assume that c 0,1 , c 1,1 > 0. We denote P ⋆ = proj ε, M (P). We have that P 1|t is a Gaussian random variable with zero mean. We now compute its covariance matrix. First, we have thatE[X t X ⊤ 1 ] = (1 -t)E[X 0 X ⊤ 1 ] + tE[X 1 X ⊤ 1 ] = [(1 -t)c 0,1 + tc 1,1 ]Id.We also have thatE[X t X ⊤ t ] = (1 -t) 2 E[X 0 X ⊤ 0 ] + t(1 -t)(E[X 1 X ⊤ 0 ] + E[X 1 X ⊤ 0 ]) + t 2 E[X 1 X ⊤ 1 ] + 2t(1 -t)Id = [(1 -t) 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 + 2t(1 -t)]Id = [1 -t 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 ]Id.Therefore, we get that for any t ∈ [0, 1] andx t ∈ R d E P 1|t [X 1 | X t = x t ] = ([(1 -t)c 0,1 + tc 1,1 ]/[1 -t 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 ])x t .Hence, we have that for any t ∈ [0, 1] andx t ∈ R d E P 1|t [X 1 | X t = x t ] -x t = ([(1 -t)c 0,1 + tc 1,1 ]/[1 -t 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 ] -1)x t = ([(1 -t)c 0,1 + tc 1,1 -1 + t 2 -t 2 c 1,1 -2t(1 -t)c 0,1 ]/[1 -t 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 ])x t = ([(1 -t)(1 -2t)c 0,1 + t(1 -t)c 1,1 -1 + t 2 ]/[1 -t 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 ])x t = (1 -t)([(1 -2t)c 0,1 + tc 1,1 -1 -t]/[1 -t 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 ])x t .So it follows thatEQUATIONNote that if we set c 0,1 = c 2 and c 1,1 = 1, we recover (Shi et al., 2023, Lemma 13 ) with σ = 2. Denote P ⋆ = proj ε, M (P). Combining ( 54) and ( 53) we get that P ⋆ is associated with (X t ) t∈[0,1] such that for any t ∈ [0, 1] we havedX t = {([(1 -2t)c 0,1 + tc 1,1 -1 -t]/[1 -t 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 ]) + ε}X t dt + √ 2dB t .Hence, we get thatX t = exp[ 1 2 G(t, c 0,1 , c 1,1 , ε)]X 0 + 2 t 0 exp[-G(s, c 0,1 , c 1,1 , ε)]ds exp[G(t, c 0,1 , c 1,1 , ε)] 1/2 Z,where Z ∼ N (0, Id) is independent from X 0 and for any t ∈ [0, 1], c 0,1 , c 1,1 , ε > 0 we haveG(t, c 0,1 , c 1,1 , ε) = 2 t 0 [(1 -2t)c 0,1 + tc 1,1 -1 -t]/[1 -t 2 + t 2 c 1,1 + 2t(1 -t)c 0,1 ]dt + 2εt.In addition, we defineg cov (c 0,1 , c 1,1 , ε) = exp[G(1, c 0,1 , c 1,1 , ε)], g var (c 0,1 , c 1,1 , ε) = exp[G(1, c 0,1 , c 1,1 , ε)] 1 + 2 1 0 exp[-G(t, c 0,1 , c 1,1 , ε)]dt .Hence, we have thatE[X 0 X ⊤ 1 ] = g cov (c 0,1 , c 1,1 , ε)Id, E[X 1 X ⊤ 1 ] = g var (c 0,1 , c 1,1 , ε)Id.Therefore, since for any n ∈ N, we have that P 2n+1 = proj ε, M (P 2n ) andP 2n+2 = P 2n+1 0,1 Q |0,1 , we define (c n 0,1 , c n 1,1 ) n∈N such that for any n ∈ N E P 2n [X 0 X ⊤ 1 ] = c n 0,1 Id, E P 2n [X 1 X ⊤ 1 ] = c n 1,1Id. Note that for any n ∈ N, we have thatE P 2n+1 [X 0 X ⊤ 1 ] = c n+1 0,1 Id, E P 2n+1 [X 1 X ⊤ 1 ] = c n+1 1,1 Id.Since P 0 = (π 0 ⊗ π 1 )Q |0,1 we get that c 0 0,1 = 0 and c 1,1 = 1. We have that for any n ∈ Nc n+1 0,1 = g cov (c n 0,1 , c n 1,1 , ε), c n+1 1,1 = g var (c n 1,1 , c n 1,1 , ε), which concludes the proof.Forward-backward. Next, we consider the following sequences of path measures (P n, ) n∈N and (P n, ) n∈N . We set P 0, = P 0, = (π 0 ⊗ π 1 )Q |0,1 . For any n ∈ N, we define P 2n+2, = P 2n+1, 0,1 Q |0,1 and P 2n+2, = P 2n+1, 0,1 Q |0,1 , i.e. P 2n+2, = proj R (P 2n+1, ) and P 2n+2, = proj R (P 2n+1, ). In addition, we define P 2n+1, = proj ε, M (P 2n, ) such that for any t ∈ [0, 1], P 2n+1, is associated with (X t ) t∈ [0, 1] wheredX t = {(E P 2n, 1|t [X 1 | X t ] -X t )/(1 -t) + εX t }dt + √ 2dB t , X 0 ∼ π 0 ,with ε ∈ R. Similarly, we define P 2n+1, = proj ε, M (P 2n, ) such that for any t ∈ [0, 1], P 2n+1, is associated with (Y 1-t ) t∈ [0, 1] wheredY t = {(E P 2n, 0|t [X 0 | Y t ] -Y t )/(1 -t) + εY t }dt + √ 2dB t , Y 0 ∼ π 1 .Proposition I.2 (Forward-Backward updates): For any n ∈ N, we have that P 2n+1, 0,1 = N (0, Σ n+1, Id) and P 2n+1, 0,1 = N (0, Σ n+1, Id) whereΣ n+1, = Id c n+1, 0,1 Id c n+1, 0,1 Id c n+1, 1,1 Id , Σ n+1, = c n+1, 0,0 Id c n+1, 0,1 Id c n+1, 0,1Id , and for any n ∈ N c n+1,1,1 = f var (c n, 0,0 , 1, c n, 0,1 , ε), c n+1, 0,1 = f cov (c n, 0,0 , 1, c n, 0,1 , ε), c n+1, 0,0 = f var (1, c n, 1,1 , c n, 0,1 , ε), c n+1, 0,1 = f cov (1, c n, 1,1 , c n, 0,1 , ε).The proof is similar to the one of Proposition I.1.Proof. Let P = (P 0,1 )Q |0,1 where P 0,1 is a Gaussian random variable with zero mean and covariance matrix Σ ∈ R 2d×2d such thatΣ = c 0,0 Id c 0,1 Id c 0,1 Id Id ,where Id is the d-dimensional identity matrix and c 0,1 , c 0,0 > 0. We denote P ⋆ = proj ε, M (P). We have that P 1|t is a Gaussian random variable with zero mean. We now compute its covariance matrix. First, we have thatE[X t X ⊤ 1 ] = (1 -t)E[X 0 X ⊤ 1 ] + tE[X 1 X ⊤ 1 ] = [(1 -t)c 0,1 + t]Id. We also have that E[X t X ⊤ t ] = (1 -t) 2 E[X 0 X ⊤ 0 ] + t(1 -t)(E[X 1 X ⊤ 0 ] + E[X 1 X ⊤ 0 ]) + t 2 E[X 1 X ⊤ 1 ] + 2t(1 -t)Id = [(1 -t) 2 c 0,0 + t 2 + 2t(1 -t)c 0,1 + 2t(1 -t)]Id = [2t -t 2 + (1 -t) 2 c 0,0 + 2t(1 -t)c 0,1 ]Id.Therefore, we get that for any t ∈ [0, 1] andx t ∈ R d E P 1|t [X 1 | X t = x t ] = ([(1 -t)c 0,1 + t]/[2t -t 2 + (1 -t) 2 c 0,0 + 2t(1 -t)c 0,1 ])x t .Hence, we have that for any t ∈ [0, 1] andx t ∈ R d E P 1|t [X 1 | X t = x t ] -x t = ([(1 -t)c 0,1 + t]/[2t -t 2 + (1 -t) 2 c 0,0 + 2t(1 -t)c 0,1 ] -1)x t = ([(1 -t)c 0,1 + t -2t + t 2 -(1 -t) 2 c 0,0 -2t(1 -t)c 0,1 ] /[2t -t 2 + (1 -t) 2 c 0,0 + 2t(1 -t)c 0,1 ])x t = ([(1 -t)(1 -2t)c 0,1 -(1 -t) 2 c 0,0 -t(1 -t)]/[2t -t 2 + (1 -t) 2 c 0,0 + 2t(1 -t)c 0,1 ])x t = (1 -t)([(1 -2t)c 0,1 -(1 -t)c 0,0 -t]/[2t -t 2 + (1 -t) 2 c 0,0 + 2t(1 -t)c 0,1 ])x t .Finally, we have that for any t ∈ [0, 1] andx t ∈ R d (E P 1|t [X 1 | X t = x t ] -x t )/(1 -t) = ([(1 -2t)c 0,1 -(1 -t)c 0,0 -t]/[2t -t 2 + (1 -t) 2 c 0,0 + 2t(1 -t)c 0,1 ])x t . (55)Note that if we set c 0,1 = c 2 and c 0,0 = 1, we recover (Shi et al., 2023, Lemma 13 ) with σ = 2. Denote P ⋆ = proj ε, M (P). Combining ( 55) and ( 53) we get that P ⋆ is associated with (X t ) t∈[0,1] such that for any t ∈ [0, 1] we havedX t = {([(1 -2t)c 0,1 -(1 -t)c 0,0 -t]/[2t -t 2 + (1 -t) 2 c 0,0 + 2t(1 -t)c 0,1 ]) + ε}X t dt + √ 2dB t .Hence, we get thatX t = exp[ 1 2 H(t, c 0,1 , c 0,0 , ε)]X 0 + (2 t 0 exp[-H(s, c 0,1 , c 0,0 , ε)]ds exp[H(t, c 0,1 , c 0,0 , ε)]) 1/2 Z,where Z ∼ N (0, Id) is independent from X 0 and for any t ∈ [0, 1], c 0,1 , c 1,1 , ε > 0 we haveH(t, c 0,1 , c 0,0 , ε) = 2 t 0 [(1-2t)c 0,1 -(1-t)c 0,0 -t]/[2t-t 2 +(1-t) 2 c 0,0 +2t(1-t)c 0,1 ]dt+2εt.In addition, we defineg cov (c 0,1 , c 0,0 , ε) = exp[ 1 2 H(1, c 0,1 , c 0,0 , ε)], g var (c 0,1 , c 0,0 , ε) = exp[H(1, c 0,1 , c 0,0 , ε)] 1 + 2 1 0 exp[-H(t, c 0,1 , c 0,0 , ε)]dt .Hence, we have thatEQUATIONRemember that P 0, = P 0, = (π 0 ⊗ π 1 )Q |0,1 . In addition, for any n ∈ N, we haveP 2n+2, = P 2n+1, 0,1 Q |0,1 and P 2n+2, = P 2n+1, 0,1 Q |0,1 , i.e. P 2n+2, = proj R (P 2n+1,) and P 2n+2, = proj R (P 2n+1, ). In addition, we also have P 2n+1, = proj ε, M (P 2n, ) and P 2n+1, = proj ε, M (P 2n, ). We also define (c n, 0,1 , c n, 1,1 ) n∈N such that for any n ∈ NE P 2n, [X 0 X ⊤ 1 ] = c n, 0,1 Id, E P 2n, [X 1 X ⊤ 1 ] = c n, 1,1 Id. Finally, we define (c n, 0,1 , c n, 1,1 ) n∈N such that for any n ∈ N E P 2n, [X 0 X ⊤ 1 ] = c n, 0,1 Id, E P 2n, [X 0 X ⊤ 0 ] = c n, 0,0 Id.n ∈ N c n+1, 0,1 = g cov (c n, 0,1 , c n, 0,0 , ε), c n+1, 1,1 = g var (c n, 0,1 , c n, 0,0 , ε), c n+1, 0,1 = g cov (c n, 0,1 , c n, 1,1 , ε), c n+1, 0,0 = g var (c n, 0,1 , c n, 1,1 , ε).In addition, we have that c n, 0,1 = c n, 0,1 = 0 and c 0, 1,1 = c 0, 0,0 = 1. This concludes the proof.Error accumulation. In Proposition I.1 and Proposition I.2, we derive the sequences corresponding to the evolution of the variance and the covariance throughout the DSBM iterations in forwardforward mode or forward-backward mode. In what follows, we showcase the behavior of these sequences for different values of ε > 0. We recall that ε corresponds to the error made in the Markov projection, i.e. proj M is replaced by proj ε, M in the forward-forward mode and proj M is replaced by proj ε, M and proj ε, M in the forward-backward mode. First, if we consider the perfect scenario, i.e. ε = 0, then we observe that both the forward-forward mode and the forward-backward mode satisfy that E P 2n [X 1 X ⊤ 1 ] = Id, see Figure 10 and Figure 11 . Additionally, we can show that in the perfect scenario, i.e. ε = 0, then both the forward-forward mode and the forward-backward mode satisfy that lim n→+∞ E P 2n [X 1 X ⊤ 0 ] = ( √ 2 -1)Id, see Figure 10 and Figure 11 . However, as ε increases the behavior between the forward-forward sequence and the forward-backward sequence significantly differs. More precisely, the error explodes as ε increases along the DSBM iteration for the forward-forward mode. On the contrary, in the forward-backward mode, the error remains bounded along the DSBM iterations, see Figure 10 and Figure 11 . In this section, we provide details on the scaling of the loss function we implement when training our online version of DSBM. We adapt the method of (Karras et al., 2022, Appendix B. 2) to the case of bridge matching. We only present our derivations in the case of the forward training of the online version of DSBM, i.e. ( 9). The preconditioning of the loss described in this setting can be readily extended to the forward-backward loss we consider in practice, i.e. the parametric version of (10).We consider the following objective function for any t ∈ [0, 1]EQUATIONWe also define for any t ∈ [0, 1] andx t ∈ R d , v θ t (x t ) = c o t nn θ t (c i t x t ) + c s t x t . Hence, c i t is an input scaling function, c ot is an output scaling function and c s t is a skip-connection function. During the training of the online version of DSBM, P will be given by P n , where P n = P v θn , where the sequence (θ n ) n∈N is given by ( 9). Here, we apply the principles of Karras et al. (2022) to the case where P = (π 0 ⊗ π 1 )Q |0,1 , i.e. at initialisation of the sequence. In what follows, we assume thatEQUATIONUsing (58), we have that for any t ∈ [0, 1]E Pt [∥X t ∥ 2 ] = (1 -t) 2 E π0 [∥X 0 ∥ 2 ] + t 2 E π1 [∥X 1 ∥ 2 ] + εt(1 -t)d = [(1 -t) 2 + t 2 + εt(1 -t)]d.We set c i t so that E[∥c i t X t ∥ 2 ] = d for every t ∈ [0, 1]. Hence, we have that for any t ∈ [0, 1] Next, we rewrite (57). For any t ∈ [0, 1] we have thatc i t = 1/ (1 -t) 2 + t 2 + εt(1 -t).ℓ t = λ t E P [∥c o t nn θ t (c i t X t ) + c s t X t -X1-Xt 1-t ∥ 2 ] = (c o t ) 2 λ t E P [∥nn θ t (c i t X t ) -[-c s t X t + X1-Xt 1-t ]/c o t ∥ 2 ] = (c o t ) 2 λ t E P [∥nn θ t (c i t X t ) -[ X1-(1+c s t (1-t))Xt 1-t ]/c o t ∥ 2 ]Hence, we get that for any t ∈ [0, 1], T t = [ X1-(1+c s t (1-t))Xt 1-t ]/c oE P1,t [∥X 1 -(1 + c s t (1 -t))X t ∥ 2 ] = (1 + c s t (1 -t)) 2 E Pt [∥X t ∥ 2 ] + E π1 [∥X 1 ∥ 2 ] -2(1 + c s t (1 -t))E P1,t [⟨X t , X 1 ⟩] = (1 + c s t (1 -t)) 2 E Pt [∥X t ∥ 2 ] + d -2(1 + c s t (1 -t))td Hence, we get that for any t ∈ [0, 1] (c o t ) 2 = ((1 + c s t (1 -t)) 2 E Pt [∥X t ∥ 2 ]/d + 1 -2(1 + c s t (1 -t))t)/(1 -t) 2 . We now minimise (c o t ) 2 with respect to (1 + c s t (1 -t)). We get that 1 + c s t (1 -t) = t/(E Pt [∥X t ∥ 2 ]/d). Hence, we get that c s t = t/[(1 -t)((1 -t) 2 + t 2 + εt(1 -t))] -1/(1 -t).With that choice, we get that for any t ∈ Karras et al. (2022) , the weighting function λ t is set so that the weight in front of the regression loss is equal to one for all times t ∈ [0, 1]. Hence, Karras et al. (2022) suggests to set λ t = 1/(c o t ) 2 . However, in practice we observe better results by letting λ t = 1. This means that the effective weight is given by 1/(c o t ) 2 . Therefore, for any t ∈ [0, 1] we have[0, 1] (c o t ) 2 = (1 -t 2 /((1 -t) 2 + t 2 + εt(1 -t)))/(1 -t) 2 . In(c i t ) 2 = (1 + (ε -2)t(1 -t)) -1 , c s t = ((ε -2)t -1)/(1 + (ε -2)t(1 -t)), (c o t ) 2 = (1 + t + (ε -2)t(1 -t))/(1 -t).In this section, we delve deeper into the specifics of each experiment, implementation details, and share additional results.We consider two ways of parameterising the vector fields: as in DSBM, we can use two separate neural networks to approximate the forward and backward vector fields, or we can use a single neural network that is conditioned on the direction. In the latter case, we do the conditioning in a similar fashion to how DDM's neural networks, U-Nets or MLPs, are conditioned on time embeddings. After all, if we work with continuous time variables t ∈ [0, 1], then the direction signal s ∈ {0, 1} can be thought of as a target time. Thus, we perform the same initial transformations on t and s, i.e. computing sinusoidal embeddings followed by a 2-layer MLP, and use the concatenated outputs in adaptive group normalisation layers (Dhariwal and Nichol, 2021; Hudson et al., 2023; Perez et al., 2018) .To optimise our networks, we use Adam (Kingma and Ba, 2015) with β = (0.9, 0.999), and we modify the gradients to keep their global norm below 1.0. We re-initialise the optimiser's state when the finetuning phase starts.All image samples in the paper are generated using EMA parameters as it has been known to increase the visual quality of resulting images (Song and Ermon, 2020) . Sampling is also the integral part of DSBM's finetuning stage, both iterative and online. Here, we have two options: sample with EMA or non-EMA parameters. The non-EMA sampling might be easier to implement, while EMA sampling results in a more stable training and slightly better quality, e.g. see AFHQ samples in Figure 23 and Figure 24 for comparison.For every model used in the paper, we provide hyperparameters in Table 3 . illustrated in Figure 13 , 2-Wasserstein distance can be very noisy even with 10K samples in the test set. To mitigate this variance, we averaged the 2-Wasserstein distance across five random sets of 10K samples per run, and then averaged these results across multiple runs. Despite these measures, we recommend a future redesign of these 2D experiments to facilitate more robust comparisons between methods. To parameterise the forward and backward drifts, we use a 2-layer MLP network with 256 hidden units. To process time variables, we compute sinusoidal time embeddings, followed by a 2-layer MLP with 256 hidden units and 50 output units. The resulting time embeddings are then concatenated with X t , so the drift networks receive 100-dimensional input vectors.For iterative DSBM finetuning, we perform 40K steps with varying number of outer iterations, i.e. when we switch between training the forward and the backward networks. Alternating every 5K steps, corresponds to 8 outer DSBM iteration. Similarly, changing the direction every 1K steps, leads to 40 outer iterations.We do not have a cache dataloader like in the original DSBM implementation2 , thus we generate training samples on the fly by sampling either from the forward or the backward network. For this simple task, we also do not use EMA.During training and evaluation, we use Euler-Maruyama method with 100 equidistant time steps between 0 and 1. The covariance is evaluated using 10K samples.Additional comparison α-DSBM vs OT bridge-matching. We consider the scalar Gaussian setting. We highlight the dependence of OT bridge matching (256/16/8) on the batch size, as the mini-batch OT coupling can be far from the true OT coupling as the batch size decreases. All experiments are run with similar compute and architecture/training hyperparameters, see Table 5 .Ablation of the hyperparameter α. Instead of letting α be determined by Adam and adaptive for α-DSBM, we explicitly set it by using Stochastic Gradient Descent (SGD) with learning parameter α; see Figure 14 . We also ran online α-DSBM with α = 10 -1 but the training diverges in this case. We use the same U-Net architecture with hyperparameters given in Table 3 .For DSBM finetuning, we perform 30 outer iterations, i.e. alternating between training the forward and the backward networks, while at each outer iteration a network is trained for 5000 steps. We do not have a cache dataloader and generate training samples on the fly by sampling either from the forward or the backward network with EMA parameters.During training and evaluation, we use Euler-Maruyama method with 30 equidistant time steps between 0 and 1. For evaluation, we compute FID based on the whole MNIST training set of 60000 examples and a set of 4000 samples that were initialised from each test image in the EMNIST dataset. MSD is computed between 4000 initial EMNIST test examples and their corresponding MNIST samples.In Figures 15-18, we provide forward and backward samples, i.e. EMNIST → MNIST and MNIST → EMNIST, from models that differ in parameterisation, finetuning methods, and sampling strategy. For all the models above, we used ε = 1. ⋆= argmin v L(v, P n ),P n+1 = argmin P∈R(Q) KL(P v n+1 ⋆ | P),since we have that P n = P v n ⋆ . Hence, our online procedure Algorithm 1, which corresponds to the discretisation of the flow of path measures (3) can be rewritten asP n+1/2 = P v n+1 ⋆ , v n+1 ⋆ = Gradientstep(L(v, P n ), P n+1 = argmin P∈R(Q) KL(P v n+1 ⋆ | P).Therefore, our proposed algorithm can be seen as an incremental version of the Maximisation-Maximisation algorithm associated with DSBM instead of an incremental version of the Expectation-Maximisation algorithm discussed in(Neal and Hinton, 1998).https://github.com/yuyang-shi/dsbm-pytorch 
 
 Appendix:  


Decision: This paper introduces a novel approach to solving Entropic Optimal Transport problems via the Schrödinger Bridge framework, demonstrating theoretical convergence and empirical improvements in image translation tasks. The reviewers generally found the work to be strong, with two reviewers (xpT6 and Ti9G) giving it an accept score of 7, showing clear support for the paper.

While Reviewer ovZP had concerns about the clarity of assumptions and comparisons with SOTA methods, these were sufficiently addressed during the discussion, leading to an increased score of 6. Reviewer KqVg also maintained a score of 6, acknowledging the incremental contributions but remaining positive overall.
Given the overall positive scores, particularly strong positive attitudes from two of the reviewers, I recommend accepting this paper.

Reviews: review: paper_summary: paper_summary: The paper proposes a novel algorithm for mass transport problems, aiming to compute maps that transport one distribution to another. This paper introduces the Schrödinger Bridge Flow algorithm, a dynamic entropy-regularized version of OT, eliminating the need to train multiple DDM-type models. The algorithm discretizes a flow of path measures, ensuring the Schrödinger Bridge as the only stationary point. The paper demonstrates the algorithm's effectiveness on various unpaired data translation tasks, showcasing its potential to solve Entropic Optimal Transport problems with improved computational efficiency and practical implementation compared to existing methodologies.
strength_weakness: strength_weakness: 1.Assumptions in Theorem 3.1: Theorem 3.1 only mentions "under mild assumptions" without specifying them in the main text. It is necessary to state these assumptions clearly in the main body of the paper and provide references to justify their rationality. While a more detailed explanation can be provided in the appendix, the main text should at least include a summary of the assumptions to ensure the reader understands the conditions under which the theorem holds.
2.Experimental Results and Comparisons: The experimental results on some datasets are not particularly impressive, even though the authors claim that their method approximates OT maps more accurately than existing methods. The results in Table 1 do not show significant progress or state-of-the-art (SOTA) performance on certain datasets. Additionally, it is confusing that the authors explicitly mentioned Rectified Flow, Flow Matching, and OT-Flow Matching in the literature review, but did not include comparisons with these methods in the experiments. For instance, Rectified Flow has conducted experiments on mass transport, and other related methods, such as Denoising Diffusion Bridge Models from ICLR 2024, have tackled similar tasks. A direct comparison of the proposed method's performance with these existing methods would provide a clearer picture of its advantages and limitations.
3.Discussion of Related Work: The paper could benefit from a more comprehensive discussion of related work, including diffusion models, flow matching, and the application of Wasserstein gradient flow-based models in mass transport tasks. This would help in better contextualizing the authors' innovations and contributions. By situating their work more clearly within the existing literature, the authors could highlight the unique aspects and potential advantages of their proposed method, enhancing the readers' understanding of its significance.
questions: questions: 1.See Weaknesses
2.I am curious about the complexity and cost of the training method proposed by the authors. Could the authors provide more detailed quantitative data, such as training time comparisons with related methods?
limitations: limitations: Although the authors proposed a novel model to address mass transport problems and provided some theoretical proofs, the method is simulation-based. I believe that the training cost is higher compared to simulation-free methods, and the model performance does not achieve state-of-the-art (SOTA) results. Additionally, there is a lack of direct comparison with existing SOTA methods.
paper_summary: paper_summary: This work introduces $\alpha$-DSBM, a new way of training DSBM-like models, which does not require a Markovian projection at each step and eliminates the need to train multiple models. The main advantage over previous DSBM-based approaches is that $\alpha$-DSBM only needs to train a single model with a single model, thus exhibiting a more stable training procedure.

The authors thoroughly contextualize their work within related work and provide detailed theoretical derivations to motivate their approach. Empirically $\alpha$-DSBM is validated by comparing to existing DSBM-based approaches across different unpaired image translation tasks as well as some toy data examples.
strength_weakness: strength_weakness: - While the experimental section of the paper thoroughly analyzes and compares DSBM and its different flavors, including the proposed method, a comparison with other competing methods could further support $\alpha$-DSBM through empirical evidence. Specifically, a comparison to bridge/flow matching with mini-batch OT similar to [1, 2], a highly optimized CycleGAN as in [3], and adversarial-based OT methods like [4] could be good candidates for further comparisons. The authors also mention these as competing methods.


[1] Alexander Tong and Kilian Fatras and Nikolay Malkin and Guillaume Huguet and Yanlei Zhang and Jarrid Rector-Brooks and Guy Wolf and Yoshua Bengio. "# Improving and generalizing flow-based generative models with minibatch optimal transport". IN TMLR 2024.

[2] Luca Eyring and Dominik Klein and Théo Uscidda and Giovanni Palla and Niki Kilbertus and Zeynep Akata and Fabian Theis. "Unbalancedness in Neural Monge Maps Improves Unpaired Domain Translation". In ICLR 2024.

[3] Dmitrii Torbunov and Yi Huang and Huan-Hsin Tseng and Haiwang Yu and Jin Huang and Shinjae Yoo and Meifeng Lin and Brett Viren and Yihui Ren. "UVCGAN v2: An Improved Cycle-Consistent GAN for Unpaired Image-to-Image Translation". In Arxiv 2023.

[4] Beomsu Kim and Gihyun Kwon and Kwanyoung Kim and Jong Chul Ye. "Unpaired Image-to-Image Translation via Neural Schrödinger Bridge". In ICLR 2024.
questions: questions: - How does $\alpha$-DSBM perform compared to bridge matching with mini-batch OT sampling empirically? While it is clear that these methods introduce significant errors because of the mini-batch approximation, it is unclear how this affects results empirically compared to $\alpha$-DSBM. I think this is the most important additional competing work (apart from DSBM), as these methods share their overall goal.
- Does it make sense to combine $\alpha$-DSBM with mini-batch OT sampling for the initial pretraining? How would this impact $\alpha$-DSBM?
limitations: limitations: - The main limitation of $\alpha$-DSBM is that new data needs to be generated with the current model during its fine-tuning, making the procedure not simulation-free and, thus, more expensive. This is sufficiently addressed in the paper.
paper_summary: paper_summary: The authors consider developing a new algorithm for a Schrödinger Bridge problem of translating between two probability distributions. Motivated by the fact that current Schrödinger Bridge (SB) approaches either use mini-batch optimal transport techniques or require training diffusion at every iteration, the authors propose discretizing the Schrödinger Bridge Flow (SB Flow) of path measures that converge to SB. They call this discretization $\alpha$-Iterative Markovian Fitting ($\alpha$-IMF) procedure. Using the fact that Markovian path measures characterizing the considered flow can be parametrized by diffusion vector field $v$, authors propose a non-parametric scheme for updating $v_n$ corresponding to $n$ discretization step of SB Flow and prove its correctness. To implement the procedure in practice, authors prove a parametric counterpart to the non-parametric update scheme that allows for optimization over a vector field $v_{\theta}$ parametrized as a neural network called $\alpha$ Diffusion Schrödinger Bridge Matching ($\alpha$-DSBM). Furthermore, the authors computationally validate the proposed algorithm on both vector and image state spaces. The authors compare their $\alpha$-DSBM algorithm to other DSBM versions [3] in the Gaussian setting. Image state space experiments were held on image-to-image translation tasks, i.e., MNIST $\leftrightarrow$ EMNIST, Wild $\leftrightarrow$ Cat on AFHQ 64/256, with quantitative evaluation and comparison with the original DSBM algorithm [3].
strength_weakness: strength_weakness: - The proposed $\alpha$-DSBM algorithm [3] is, in some sense, an advancement of the previously known DSBM algorithm [3]. It is expected to be more efficient because it removes the necessity to learn a distinct model at each iteration. However, the paper lacks a study of image quality boost given the same computational budget. Thus, it is not clear whether the authors achieved the goal of developing a more efficient algorithm.
- The convergence of parametric $\alpha$-DSBM is not proved. While it may be a hard theoretical problem to prove, at least a more extensive empirical evaluation of the setup with the known ground truth solution would be beneficial. The presented setup with a scalar covariance matrix is too simple. It would be a more solid argument to evaluate the Gaussian setup with a full covariance matrix [6] or use a mixture benchmark [7] and compare it with other neural network methods. 
- The FID metric is similar for DSBM and $\alpha$-DSBM for the AFHQ experiment in Table 1 considering measured standard deviation. Thus again highlighting the lack of study with the same computational budget to justify that $\alpha$-DSBM is more efficient.
The chosen AFHQ dataset may not be the best choice considering the small test set, which may introduce bias in measuring the FID metric. In turn, the visual comparison is hard since most images produced by DSBM and $\alpha$-DSBM and given in Figure 19 are similar, and it is hard to deduce which algorithm is better. A bigger dataset such as Celeba [5] may help to solve this issue. At least a bigger dataset would give a good estimate of FID metrics. 
- It seems that some hyperparameters for DSBM AFHQ experiments are missing. The number of pre-training iterations, grad updates per IMF iteration, and IMF iterations are not specified. What is the training time and computational budget? 
- There is no study on *cornerstone* hyperparameter $\alpha$. What is the value of $\alpha$ used in the experiments? Is it equal to the learning rate of the diffusion model parameters?
- The absence of code may cause trouble with reproducing the presented results.
questions: questions: - Since the AFHQ dataset is quite small, experiments on unpaired male-to-female translation on the Celeba [5] dataset could provide a more meaningful comparison (as it was in DSBM[3]).
- How do DSBM and $\alpha$-DSBM compare in terms of computational budget? Could authors provide a quantitative comparison?
- How does variation of $\alpha$ affect speed of convergence and deviation from the SB solution?
- Can authors provide experimental results of $\alpha$-DSBM on more complex Gaussian distributions or multimodal low dimensional data, i.e., gaussian mixtures? 
- Is it possible to extend the proposed approach to GSBM [4], where, as far as my understanding goes, authors change the reciprocal projection step compared to DSBM [3]?
- Can this method with small $\epsilon$ can be applied to the generation problem? This could drastically reduce the computational budget of such and be applied in modern text-to-image generative models [2] similar to the work [1]. It would be great to see such experimental results. 
- Is it possible to perform *parametric* updates for Sinkhorn flow in a similar to $\alpha$-DSBM way?

Overall, the paper looks very promising. I will consider increasing my score if the authors address some of the listed weaknesses and questions, especially related to the same computational budget comparison, and evaluate their method on one of the listed more complex setups with the known ground truth distributions.

[1] Liu, Xingchao, et al. "Instaflow: One step is enough for high-quality diffusion-based text-to-image generation." The Twelfth International Conference on Learning Representations. 2023.

[2] Esser, Patrick, et al. "Scaling rectified flow transformers for high-resolution image synthesis." Forty-first International Conference on Machine Learning. 2024.

[3] Shi, Yuyang, et al. "Diffusion Schrödinger bridge matching." Advances in Neural Information Processing Systems 36 (2024).

[4] Liu, Guan-Horng, et al. "Generalized Schrödinger Bridge Matching." The Twelfth International Conference on Learning Representations.

[5] Liu Z. et al. Deep learning face attributes in the wild //Proceedings of the IEEE international conference on computer vision. – 2015. – С. 3730-3738.

[6] Hicham Janati, Boris Muzellec, Gabriel Peyré, and Marco Cuturi. Entropic optimal transport between unbalanced gaussian measures has a closed form. Advances in neural information processing systems, 33:10468–10479, 2020.

[7] Gushchin N. et al. Building the bridge of schrödinger: A continuous entropic optimal transport benchmark //Advances in Neural Information Processing Systems. – 2023. – Т. 36. – С. 18932-18963.
limitations: limitations: The authors adequately addressed the limitations and potential negative societal impact of their work.
paper_summary: paper_summary: This paper proposes $\alpha$-IMF, which is an incremental version of IMF method in DSBM paper. Specifically, this work demonstrate the convergence properties of $\alpha$-IMF and implement it through an online learning method. Here, $\alpha$ is implicitly reflected. Moreover, the functional flow interpretation provides a basis for online learning, allowing for online learning as opposed to the iterative approach of DSBM. From an engineering perspective, the authors suggest parameterizing forward and backward control simultaneously by adding a forward/backward indicator (0 or 1) as an input variable to the neural network. The method shows feasible results on toy (Gaussian) data and image-to-image translation tasks.
strength_weakness: strength_weakness: - The experiment on practical data is insufficient and the performance improvement is incremental. The paper only presents real-world I2I on Cat <-> Wild. Moreover, the performance of bidirectional online learning improves only incrementally. Totally, it needs more evaluation on other kinds of dataset, e.g. Male <-> Female / Handbags <-> Shoes. 
- It would be beneficial to quantify the advantages of bidirectional learning or online learning over DSBM. For instance, comparing the actual wall-clock time or GPU memory usage would be useful. 
- The current algorithm reflects $\alpha$ implicitly, and so there is no analysis provided on $\alpha$. If there were any intuition or insights into the role of $\alpha$ (step size) or if an algorithm could be developed to explicitly control $\alpha$, it would offer more valuable insights.
questions: questions: - If we perform this algorithm with full gradient descent algorithm, is it possible to explicitly control $\alpha$?
- Can we develop exact algorithm that reflect $\alpha$ explicitly? If so, can we obtain some experimental results that compares $\alpha$?
limitations: limitations: The authors adequately addressed the limitations.

Comments: title: Reminder: Respond to Reviewer Questions by August 13, 11:59 PM AoE
comment: Thank you for your detailed reply and additional results. It addresses most of my concerns, and  I updated my score. I suggest the authors to include these discussions and the additional results in their revision.
We thank the reviewer for their additional comment to which we answer below. 

> Thank you for addressing many of the initial concerns in your rebuttal. However, I believe that if the main contribution is claimed to be a simpler implementation with only one neural network and faster convergence, it is crucial to provide more quantitative data to support these claims.

We will clarify this point in our introduction but we believe that our contribution is two-fold.  

* First, we introduce a new theoretical framework to analyze online schemes aimed at solving the Schrodinger Bridge (SB). To the best of our knowledge, this is the first time an online training algorithm is proposed to solve SB along with a theoretical framework. This is what motivates the $\alpha$-IMF scheme.
* Second, we introduce our practical implementation of $\alpha$-IMF, $\alpha$-DSBM. With the simplification given by $\alpha$-IMF this can be seen as a fine-tuning approach of any bridge/flow matching procedure. To further improve the scalability of our approach we only use one bidirectional network compared to existing implementations which use two networks.

In the rebuttal we conducted additional experiments on CelebA to validate our findings. Because of the lack of space we could not include evaluation metrics in our initial rebuttal but report them below. We compare DSBM and $\alpha$-DSBM in the setting where both are evaluated with a single network (in order to answer the reviewer comment about “compare your model with other single neural network setups”). Both DSBM and $\alpha$-DSBM are run with the same number of gradient evaluations. The hardware set up and hyperparameters are similar to our AFHQ experiment. We report both alignment metrics and quality metrics.

Regarding the L2, LPIPS evals (lower is better):
* DSBM: 0.159 (L2) / 0.451 (LPIPS)
* alpha-DSBM: 0.05 (L2) / 0.376 (LPIPS)

We also report the Inception Score (higher is better) evaluations for the base models, DSBM and alpha-DSBM:
* Base training: 2.29
* DSBM: 2.88
* alpha-DSBM: 3.13

We hope that this additional experiment resolves the reviewer’s concerns.

> Providing details on the reduction in neural network parameters to quantitatively demonstrate the simplicity of your approach.

We highlighted in the rebuttal that in our implementation we halve the number of parameters compared to DSBM. The results regarding the use of a bidirectional or a unidirectional network are reported in Table 1 in the original manuscript.

> Clarifying the convergence rates in comparison with baseline methods, particularly considering the fine-tuning duration mentioned in your rebuttal experiments (I noticed that in the rebuttal experiments, the authors' fine-tuning time was much longer than the training time. Does this contradict the authors' statement?). 

We want to highlight that the finetuning is done for significantly fewer training steps than the initial base training (more precisely in the case of CelebA for instance: 100k base training and 10k finetuning). The reason the finetuning takes  longer than the base training is because our method (as well as DSB(M) and related methods) is not _sampling-free_. This is a well-known limitation of this line of work and we do not claim to solve this problem in our paper. We acknowledge this limitation in the original discussion of our paper. We agree that deriving “sampling-free” or almost sampling-free improvements of our methodology is an interesting future work. 

> the theoretical assumptions are further clarified,

We have detailed the theoretical assumptions in our original answer to the reviewer. If any of the discussed assumptions remain unclear, we are willing to clarify them. We highlight that the detailed discussion of these theoretical assumptions will be included in the revised version of the paper as well as the discussion regarding the parameters reduction of the method and the additional comparison with DSBM and related methods.
Thanks for the answers and added experiments. It would be very interesting to see comparisons to [1,2,4] also in the image translation settings. I think this setting would provide a more meaningful empirical comparison. Also thanks for pointing to [1,2] for the convergence of DSBM with mini-batch OT.
Thank you for addressing many of the initial concerns in your rebuttal. However, I believe that if the main contribution is claimed to be a simpler implementation with only one neural network and faster convergence, it is crucial to provide more quantitative data to support these claims.

I recommend:

Including additional experiments that compare your model with other single neural network setups, highlighting computational efficiency and performance (In the rebuttal, the author provides some additional comparative experiments).
Providing details on the reduction in neural network parameters to quantitatively demonstrate the simplicity of your approach.
Clarifying the convergence rates in comparison with baseline methods, particularly considering the fine-tuning duration mentioned in your rebuttal experiments (I noticed that in the rebuttal experiments, the authors' fine-tuning time was much longer than the training time. Does this contradict the authors' statement?).
If these points are addressed in the revised version and the theoretical assumptions are further clarified, I would be inclined to increase my score to a 5-6. Thank you for your efforts to improve the manuscript.
Thank you for your reply and additional results. Since my concerns and questions were addressed, I updated my score.
We thank you for your review and appreciate your time reviewing our paper.
The end of the discussion period is close. We would be grateful if we could hear your feedback regarding our answers to the reviews. We are happy to address any remaining points during the remaining discussion period.
Thanks a lot in advance!
Thanks a lot for your comments which have improved the overall quality of the paper. 
Here are the additional metrics requested below (the hardware setup is identical to AFHQ):
* Base model training time: 16 hours
* Finetuning alpha-DSBM/Finetuning DSBM: 7 hours

The finetuning of DBSM corresponds to two IMF iterations (one forward and one backward).

Regarding the L2, LPIPS evals (lower is better):
* DSBM: 0.159 (L2) / 0.451 (LPIPS)
* alpha-DSBM: 0.05 (L2) / 0.376 (LPIPS)

We are also reporting the Inception Score (higher is better) evals for the base models, DSBM and alpha-DSBM:
* Base training: 2.29
* DSBM: 2.88
* alpha-DSBM: 3.13

Hence, alpha-DSBM improves over DSBM for the same amount of compute. We see that the main reason for the underperformance of DSBM is that for that amount of compute there is little deviation of the alignement score, i.e. the model is still close to the base model. We will incorporate these numbers in the final version of the paper. We hope that these numbers resolve the concerns of the reviewer.
Dear Authors,
\
\
A quick reminder that the discussion phase with reviewers is open until **August 13, 11:59 PM AoE**. Please take some time to engage with Reviewer Ti9G and see if you can address the raised concerns in response to your rebuttal. Your engagement is crucial for ensuring that the reviewers have all the information they need to make informed decisions.
\
\
Thank you for your continued efforts and contributions to the review process.
\
Best regards,  
Area Chair, NeurIPS 2024
Thank you for the response. I would like to keep my score to 6.
I thank the authors for their reply, which resolves most of my concerns and questions. To finalize, could you please provide the metrics for the quality (e.g., FID) and similarity (e.g., L2 cost or LPIPS) for the Celeba dataset experiment? Also, could you please provide estimates on the GPU time you used to train DSBM and alpha-DSBM?

Please incorporate Celaba, Guassian, text2img experiments, and other changes into the final revisions.