Diffusion & Adversarial Schrödinger Bridges via Iterative Proportional Markovian Fitting

## Abstract

Abstract The Iterative Markovian Fitting (IMF) procedure based on iterative reciprocal and Markovian projections has recently been proposed as a powerful method for solving the Schrödinger Bridge problem. However, it has been observed that for the practical implementation of this procedure, it is crucial to alternate between fitting a forward and backward time diffusion at each iteration. Such implementation is thought to be a practical heuristic, which is required to stabilize training and obtain good results in applications such as unpaired domain translation. In our work, we show that this heuristic closely connects with the pioneer approaches for the Schrödinger Bridge based on the Iterative Proportional Fitting (IPF) procedure. Namely, we find that the practical implementation of IMF is, in fact, a combination of IMF and IPF procedures, and we call this combination the Iterative Proportional Markovian Fitting (IPMF) procedure. We show both theoretically and practically that this combined IPMF procedure can converge under more general settings, thus, showing that the IPMF procedure opens a door towards developing a unified framework for solving Schrödinger Bridge problems.

## 1 Introduction

Diffusion models inspired by the Schrödinger Bridge (SB) theory, which connects stochastic processes with the optimal transport theory, have recently emerged as a powerful approach for numerous applications in biology (Tong et al., ; Bunne et al., ), chemistry (Somnath et al., ; ), computer vision (Liu et al., ; Shi et al., ) and speech processing (Chen et al., ). Most of these applications are dedicated either to supervised translation, e.g., image super-resolution and inpainting (Liu et al., ) or unpaired domain translation, such as image style-transfer (Shi et al., ) or single-cell data analysis (Tong et al., ).

In this paper, we focus specifically on unpaired domain translation, where SB-based algorithms are typically used since they enforce two key properties: the similarity between input and translated object (referred to as the optimality property) and that the input domain is translated to the target domain (referred to as the marginal matching property).

Early works (De Bortoli et al., ; Vargas et al., ; Chen et al., ) on using Schrödinger Bridge for unpaired domain translation employed the well-celebrated Iterative Proportional Fitting (IPF) procedure (Kullback, ), also known as the Sinkhorn algorithm (Cuturi & Doucet, ). The IPF procedure is initialized with a simple prior process that satisfies the optimality property. It then refines this process iteratively through optimality-preserving transformations until the marginal matching property is achieved. In each iteration, IPF decreases the forward KL-divergence $\text{KL}(q^{*}||q)$ between the current approximation $q$ and the ground-truth Schrödinger Bridge $q^{*}$. However, due to approximation errors, in practice, IPF may suffer from the ”prior forgetting”, where the marginal matching property is achieved, but the optimality is lost (Vargas et al., ; ).

Recently, the Iterative Markovian Fitting (IMF) procedure (Shi et al., ; Peluchetti, ; Gushchin et al., ) was introduced as the promising competitor to IPF. Contrary to IPF, the IMF starts from the stochastic process with the marginal matching property, while the optimality is achieved during the IMF iterations. Unlike the IPF procedure, the IMF procedure at each iteration decreases reverse KL-divergence $\text{KL}(q||q^{*})$ between the current approximation $q$ and the ground-truth Schrödigner Bridge $q^{*}$. This iterative procedure can also be seen as a generalization of rectified flows (Liu et al., ) to stochastic processes, which is used (Liu et al., ; Yan et al., ) for the modern foundational generative models such as Stable Diffusion 3 (Esser et al., ). In analogy to IPF, the IMF procedure may also accumulate errors but in approximating data distributions due to a nonperfect fit at each iteration, leading to losing marginal matching property.

In practice, the IMF is implemented as a bidirectional procedure alternating between learning forward and backward processes either by diffusion-based models in the Diffusion Schrödinger Bridge Mathcing (DSBM) algorithm (Shi et al., ) or GANs in Adversarial Schrödinger Bridge Matching (ASBM) algorithm (Gushchin et al., ). This heuristic of alternating between learning forward and backward processes helps to stabilize the IMF training and overcome the accumulation of errors and loss of marginal matching property. In this work, we explore the theoretical aspects of this heuristic approach and make the following key contributions.

Main contributions:
We show that the heuristic bidirectional IMF procedure used in practice closely relates to IPF. Namely, we discover that it, in fact, secretly utilizes IPF iterations. Due to this, we propose to call this bidirectional IMF procedure Iterative Proportional Markovian Fitting (IPMF) and conjecture that it not only converges under more general settings than was previously thought but opens a promising way of developing a unified framework for solving the SB problem.

- 1.
Theory I. We theoretically prove that the IPMF procedure converges for 1-dimensional Gaussian distributions \wasyparagraph3.2. Even this proof is non-trivial and involves significant complexity, as IPF and IMF minimize different divergences, leading to interference. Specifically, the IPF decreases forward KL divergence at each iteration but may increase the reverse KL divergence, while IMF decreases forward KL divergence and may increase the reverse KL divergence.
- 2.
Practice I. We empirically validate the convergence of IPMF in the multi-dimensional Gaussian setup, for which closed formulas can be derived in discrete time \wasyparagraph4.1.
- 3.
Practice II.
We empirically validate through a series of standard experiments, including toy 2D setups \wasyparagraph4.2, the Schrödinger Bridge benchmark \wasyparagraph4.3, colored MNIST images and real-world image dataset - Celeba \wasyparagraph4.4, that IPMF procedure converges.

These contributions show that the IPMF procedure has significant potential to combine many previously introduced SB methods, including IPF and IMF-based, with both discrete (Gushchin et al., ; De Bortoli et al., ) and continuous time (Shi et al., ; Peluchetti, ; Vargas et al., ), together with their online versions (De Bortoli et al., ; Peluchetti, ; Karimi et al., ). Moreover, the forward-backward framework of IPMF could enable rectified flows to avoid error accumulation, making them even more powerful in diffusion-based generative modeling.

Notations. We fix $N\geq 1$ intermediate time moments $0=t_{0}<t_{1}<\dots<t_{N}<t_{N+1}=1$ together with $t_{0}=0$ and $t_{N+1}=1$. We consider discrete stochastic processes with a finite second moment, entropy, and those time-moments as the elements of the set $\mathcal{P}_{2,ac}(\mathbb{R}^{D\times(N+2)})$ of probability distributions on $\mathbb{R}^{D\times(N+2)}$. For any such $q\in\mathcal{P}_{2,ac}(\mathbb{R}^{D\times(N+2)})$, we write $q(x_{0},x_{t_{1}},\dots,x_{t_{N+1}})$ to denote its density at a point $(x_{0},x_{t_{1}},\dots,x_{t_{N}},x_{1})\in\mathbb{R}^{D\times(N+2)}$. For convenience we also use the notation $x_{\text{in}}=(x_{t_{1}},\dots,x_{t_{N}})$ to denote the vector of all intermediate-time variables.

For considering the continuous version of Schrödinger Bridge we denote by $\mathcal{P}(C([0,1]),\mathbb{R}^{D})$ the set of continuous stochastic processes with time $t\in[0,1]$, i.e., the set of distributions on continuous trajectories $f:[0,1]\rightarrow\mathbb{R}^{D}$. We use $dW_{t}$ to denote the differential of the standard Wiener process. We denote by $p^{T}\in\mathcal{P}(\mathbb{R}^{D\times(N+2)})$ the discrete process which is the finite-dimensional projection of $T$ to time moments $0=t_{0}<t_{1}<\dots<t_{N}<t_{N+1}=1$. In what follows, KL is a short notation for the Kullback-Leibler divergence, while $H$ is a short notation for the differential entropy.

## 2 Background

This section recalls the Schrödinger Bridge (SB) problem (\wasyparagraph2.1). Then, we describe procedures for solving SB: Iterative Proportional Fitting (IPF, \wasyparagraph2.2) and Iterative Markovian Fitting (IMF, \wasyparagraph2.3). We discuss the practical implementation of the IMF procedure called Bidirectional IMF in \wasyparagraph2.4.

The Schrödinger Bridge problem can be considered both in discrete and continuous time, which are equivalent. The discrete-time and continuous-time versions of IPF and IMF procedures are also similar. In the main text, we operate only with the discrete-time setup to avoid a lot of repetitions, which will harm the flow of the paper. Also, the discrete-time version provides the explicit formulas for the projections used in IPF and IMF procedures, which makes it sufficiently easier to explain the main idea of our paper. We present the analogical facts about the continuous setup in Appendix A.

### 2.1 Schrödinger Bridge (SB) Problem

Schrödinger Bridge problem formulation. To formulate the Schrödinger Bridge problem we consider the Wiener process $W^{\epsilon}$ with the volatility $\epsilon>0$ which starts at some distribution $p_{0}$, i.e., the process given by the stochastic differential equation (SDE): $dx_{t}=\sqrt{\epsilon}dW_{t}$, $x_{0}\sim p_{0}$. We denote by $W^{\epsilon}_{|x_{0},x_{1}}$ the stochastic process $W^{\epsilon}$ conditioned on values $x_{0},x_{1}$ at times ${t=0,1}$, respectively. This process $W^{\epsilon}_{|x_{0},x}$ is called the Brownian Bridge (Ibe, , Chapter 9). The Schrödinger Bridge problem (Schrödinger, ) with the Wiener prior between distributions $p_{0}$ and $p_{1}$ in the discrete-time setting (De Bortoli et al., ) formulates as follows:

$$ $\min_{q\in\Pi(p_{0},p_{1})}\text{KL}(q(x_{0},x_{\text{in}},x_{1})||p^{W^{ \epsilon}}(x_{0},x_{\text{in}},x_{1})),$ (1) $$

Here $p^{W^{\epsilon}}(x_{0},x_{\text{in}},x_{1})$ is the time-discretization of $W^{\epsilon}$, which is given by $p^{W^{\epsilon}}(x_{0},x_{\text{in}},x_{1})=p_{0}(x_{0})\prod_{n=1}^{N+1}
\mathcal{N}(x_{t_{n}}|x_{t_{n-1}},\epsilon(t_{n}-t_{n-1})I_{D})$. In turn, $\Pi(p_{0},p_{1})\subset\mathcal{P}_{2,ac}(\mathbb{R}^{D\times(N+2)})$ is the subset of discrete-stochastic processes with marginals $q(x_{0})=p_{0}(x_{0})$ and $q(x_{1})=p_{1}(x_{1})$.

Static Schrödinger Bridge problem.
One may decompose the objective (1) as follows:

$$ $\displaystyle\text{KL}(q(x_{0},x_{\text{in}},x_{1})||p^{W^{\epsilon}}(x_{0},x_ {\text{in}},x_{1}))=$ $\displaystyle\text{KL}(q(x_{0},x_{1})||p^{W^{\epsilon}}(x_{0},x_{1}))+\int \text{KL}(q(x_{\text{in}}|x_{0},x_{1})||p^{W^{\epsilon}}(x_{\text{in}}|x_{0},x _{1}))q(x_{0},x_{1})dx_{0}dx_{1}.$ (2) $$

i.e., KL divergence between $q$ and $p^{W^{\epsilon}}$ is a sum of two terms: the 1st represents the similarity of the processes’ joint, marginal distributions at start and finish times $t=0,1$, while the 2nd term represents the average similarity of conditional distributions $q(x_{\text{in}}|x_{0},x_{1})$ and $p^{W^{\epsilon}}(x_{\text{in}}|x_{0},x_{1})$. Since conditional distributions $q(x_{\text{in}}|x_{0},x_{1})$ can be chosen independently of $q(x_{0},x_{1})$ we can consider $q(x_{\text{in}}|x_{0},x_{1})=p^{W^{\epsilon}}(x_{\text{in}}|x_{0},x_{1})$. In this case $\text{KL}(q(x_{\text{in}}|x_{0},x_{1})||p^{W^{\epsilon}}(x_{\text{in}}|x_{0},x
_{1}))=0$ for every $x_{0},x_{1}$ and it leads to the Static Schrödinger Bridge problem:

$$ $\displaystyle\min_{q\in\Pi(p_{0},p_{1})}\!\!\text{KL}\big{(}q(x_{0},x_{1})||p^ {W^{\epsilon}}(x_{0},x_{1})\big{)},$ (3) $$

In turn, the static SB objective can be expanded as (Gushchin et al., , Eq. 7):

$$ $\displaystyle\text{KL}(q(x_{0},x_{1})||p^{W^{\epsilon}}(x_{0},x_{1}))=\int \frac{||x_{1}-x_{0}||^{2}}{2\epsilon}dq(x_{0},x_{1})-H(q(x_{0},x_{1}))+C,$ (4) $$

which is up to an additive constant is equivalent to the objective of entropic optimal transport (EOT) problem with the quadratic cost (Cuturi, ; Peyré et al., ; Léonard, ; Genevay, ).

### 2.2 Iterative Proportional Fitting (IPF)

Several first works on Schrödinger Bridge (Vargas et al., ; De Bortoli et al., ) propose methods based on the IPF procedure (Kullback, ). The IPF-based algorithm is started by setting the process $q^{0}(x_{0},x_{\text{in}},x_{1})=p_{0}(x_{0})p^{W^{\epsilon}}(x_{\text{in}},x_
{1}|x_{0})$. Then, the algorithm alternates between two types of IPF projections $\text{proj}_{1}$ and $\text{proj}_{0}$ which are given by (De Bortoli et al., , Proposition 2):

$$ $\displaystyle q^{2k+1}=\text{proj}_{1}\big{(}\underbrace{q^{2k}(x_{1})\prod_{n =0}^{N}q^{2k}(x_{t_{n}}|x_{t_{n+1}})}_{q^{2k}(x_{1})q^{2k}(x_{0},x_{\text{in}} |x_{1})}\big{)}\stackrel{{\scriptstyle\text{def}}}{{=}}p_{1}(x_{1})\underbrace {\prod_{n=0}^{N}q^{2k}(x_{t_{n}}|x_{t_{n+1}})}_{q^{2k}(x_{0},x_{\text{in}}|x_{ 1})},$ (5) $\displaystyle q^{2k+2}=\text{proj}_{0}\big{(}\underbrace{q^{2k+1}(x_{0})\prod_ {n=1}^{N+1}q^{2k+1}(x_{t_{n}}|x_{t_{n-1}})}_{q^{2k+1}(x_{0})q^{2k+1}(x_{\text{ in}},x_{1}|x_{0})}\big{)}\stackrel{{\scriptstyle\text{def}}}{{=}}p_{0}(x_{0}) \underbrace{\prod_{n=1}^{N+1}q^{2k+1}(x_{t_{n}}|x_{t_{n-1}})}_{q^{2k+1}(x_{ \text{in}},x_{1}|x_{0})}.$ (6) $$

Thus, these projections replace marginal distributions $q(x_{1})$ and $q(x_{0})$ in $q(x_{0},x_{\text{in}},x_{1})$ by $p_{1}(x_{1})$ and $p_{0}(x_{0})$ respectively. This sequence $q^{k}$ monotonically decreases the forward KL-divergence $\text{KL}(q^{*}||q^{k})$ at each iteration and converges to the solution of the Schrödinger Bridge $q^{*}$. However, since the prior process $p^{W^{\epsilon}}$ is used only at the initialization, the nonperfect fit in practice at some iteration may lead to deviation from the SB solution. This problem is called ”prior forgetting” and was discussed in (Vargas et al., , Appendix E.3). The continuous analog of the IPF procedure is considered in (Vargas et al., ) and uses inversions of diffusion processes (see Appendix A.2).

### 2.3 Iterative Markovian Fitting (IMF)

The Iterative Matkovian Fitting (IMF) procedure (Peluchetti, ; Shi et al., ) was recently proposed as a strong competitor to the IPF, which does not suffer from the ”prior forgetting” problem of IPF. In turn, the discrete-time analog of IMF (D-IMF) has been recently proposed by (Gushchin et al., ) to accelerate the inference of the process learned by IMF. The procedure is initialized with any process $q^{0}\in\Pi(p_{0},p_{1})$. Then the procedure alternates between reciprocal $\text{proj}_{\mathcal{R}}$ and Markovian $\text{proj}_{\mathcal{M}}$ projections:

$$ $\displaystyle q^{2k+1}=\text{proj}_{\mathcal{R}}(q^{2k})\stackrel{{ \scriptstyle\text{def}}}{{=}}q^{2k}(x_{0},x_{1})p^{W^{\epsilon}}(x_{\text{in}} |x_{0},x_{1}),$ (7) $\displaystyle q^{2k+2}=\text{proj}_{\mathcal{M}}(q^{2k+1})\stackrel{{ \scriptstyle\text{def}}}{{=}}\underbrace{q^{2k+1}(x_{0})\prod_{n=1}^{N+1}q^{2k +1}(x_{t_{n}}|x_{t_{n-1}})}_{\text{forward representation}}=\underbrace{q^{2k+ 1}(x_{1})\prod_{n=0}^{N}q^{2k+1}(x_{t_{n}}|x_{t_{n+1}})}_{\text{backward representation}}$ (8) $$

Thus, the reciprocal projection $\text{proj}_{\mathcal{R}}(q)$ creates a new (in general, non-markovian) process by using the joint distribution $q(x_{0},x_{1})$ and $p^{W^{\epsilon}}(x_{\text{in}}|x_{0},x_{1})$. The latter is called the discrete Brownian Bridge. In turn, the Markovian projection $\text{proj}_{\mathcal{M}}(q)$ uses the set of transitional densities $\{q(x_{t_{n}}|x_{t_{n-1}})\}$ or $\{q(x_{t_{n}}|x_{t_{n+1}})\}$ to create a new markovian process starting from $q(x_{0})$ or $q(x_{1})$ respectively. Unlike the IPF procedure, this sequence $q^{k}$ monotonically decreases the reverse KL-divergence objective $\text{KL}(q^{k}||q^{*})$ at each iteration and converges to the solution of the Schrödinger Bridge $q^{*}$. The continuous time version of the IMF procedure is considered in (Shi et al., ; Peluchetti, ) and uses similar Markovian and reciprocal projections (Appendix A).

### 2.4 Bidirectional IMF

Since the result of the Markovian projection (8) can be represented both by forward and backward representation, in practice, neural networks $\{q_{\theta}(x_{t_{n}}|x_{t_{n-1}})\}$ (forward parametrization) or $\{q_{\phi}(x_{t_{n}}|x_{t_{n+1}})\}$ (backward parametrization) are used to learn the corresponding transitional densities. In turn, starting distributions are set to be $q_{\theta}(x_{0})=p_{0}(x_{0})$ for forward parametrization and $q_{\phi}(x_{1})=p_{1}(x_{1})$ for the backward parametrization. In practice, the alternation between forward and backward representations of Markovian processes is used in both implementations of continuous-time IMF by DSBM algorithm (Shi et al., , Algorithm 1) based on diffusion models and discrete-time IMF by ASBM algorithm (Gushchin et al., , Algorithm 1) based on the GANs. So, this bidirectional procedure can be described as follows:

$$ $\displaystyle q^{4k+1}=\underbrace{q^{4k}(x_{0},x_{1})p^{W^{\epsilon}}(x_{ \text{in}}|x_{0},x_{1})}_{\text{proj}_{\mathcal{R}}(q^{4k})},\quad q^{4k+2}= \underbrace{p(x_{1})\prod_{n=0}^{N}q_{\phi}^{4k+1}(x_{t_{n-1}}|x_{t_{n}})}_{ \text{backward parametrization}},$ (9) $\displaystyle q^{4k+3}=\underbrace{q^{4k+2}(x_{0},x_{1})p^{W^{\epsilon}}(x_{ \text{in}}|x_{0},x_{1})}_{\text{proj}_{\mathcal{R}}(q^{4k+1})},\quad q^{4k+4}= \underbrace{p(x_{0})\prod_{n=1}^{N+1}q_{\theta}^{4k+3}(x_{t_{n}}|x_{t_{n-1}})} _{\text{forward parametrization}}.$ (10) $$

Thus, only one marginal is perfectly fitted, e.g., $q_{\theta}(x_{0})=p_{0}(x_{0})$ in the case of forward representation, while the other marginal is only learned, e.g., $q_{\theta}(x_{1})=\int\underbrace{q_{\theta}(x_{0})}_{=p_{0}(x_{0})}\prod_{n=1
}^{N+1}q_{\theta}(x_{t_{n}}|x_{t_{n-1}})dx_{0}dx_{1}\cdots dx_{N}\approx p_{1}
(x_{1})$. It was observed that such approximation errors do not accumulate in this bidirectional version of IMF (Shi et al., ; Peluchetti, ; Gushchin et al., ), while the usage of only forward or backward parametrization accumulates errors and lead to divergence (De Bortoli et al., , Appendix I).

## 3 Iterative Proportional Markovian Fitting (IPMF)

In this section, we show that the heuristical procedure of bidirectional IMF \wasyparagraph2.4 is, in fact, the alternation implementation of both IPF and IMF projections and state that this heuristic in fact defines the new unified Iterative Proportional Markovian Fitting (IPMF) procedure \wasyparagraph3.1. Next, in section \wasyparagraph3.2, we provide the theoretical analysis of convergence of this IPMF procedure together with the determination of the convergence rate in the case of 1-dimensional Gaussian.

### 3.1 Bidirectional IMF is IPMF

Here, we analyze theoretically what the heuristical bidirectional IMF does. We recall that the IPF projections of proj${}_{0}(q)$ given by (6) and proj${}_{1}(q)$ given by (5) of the markovian process $q$ is in fact just change the starting distribution from $q(x_{0})$ to $p_{0}(x_{0})$ and $q(x_{1})$ to $p_{1}(x_{1})$.
Now we note that the process $q^{4k+2}$ in (9) is obtained by using a combination of Markovian projection projM given by (8) in forward parametrization and IPF projection proj1 given by (5):

$$ $\displaystyle q^{4k+2}=p(x_{1})\prod_{n=0}^{N}q^{4k+1}(x_{t_{n}}|x_{t_{n+1}})= \underbrace{\text{proj}_{1}\big{(}q^{4k+1}(x_{1})\prod_{n=0}^{N}q^{4k+1}(x_{n} |x_{n+1})\big{)}}_{\text{proj}_{1}(\text{proj}_{\mathcal{M}}(q^{4k+1}))}.$ $$

In turn, the process $q^{4k+4}$ in (10) is obtained by using a combination of Markovian projection projM given by (8) in backward parametrization and IPF projection proj0 given by (6):

$$ $\displaystyle q^{4k+3}=p(x_{0})\prod_{n=1}^{N+1}q^{4k+3}(x_{t_{n}}|x_{t_{n-1}} )=\underbrace{\text{proj}_{0}\big{(}q^{4k+3}(x_{0})\prod_{n=1}^{N+1}q^{4k+3}(x _{n}|x_{n-1})\big{)}}_{\text{proj}_{0}(\text{proj}_{\mathcal{M}}(q^{4k+3}))}.$ $$

Thus, we can represent the bidirectional procedure IMF given by (10) and (9) as follows:

$$ Iterative Proportional Markovian Fitting (Discrete time setting) $\displaystyle q^{4k+1}=\underbrace{q^{4k}(x_{0},x_{1})p^{W^{\epsilon}}(x_{ \text{in}}|x_{0},x_{1})}_{\text{proj}_{\mathcal{R}}(q^{4k})},\quad q^{4k+2}= \underbrace{p(x_{1})\prod_{n=0}^{N}q^{4k+1}(x_{t_{n-1}}|x_{t_{n}})}_{\text{ proj}_{1}(\text{proj}_{\mathcal{M}}(q^{4k+1}))},$ $\displaystyle q^{4k+3}=\underbrace{q^{4k+2}(x_{0},x_{1})p^{W^{\epsilon}}(x_{ \text{in}}|x_{0},x_{1})}_{\text{proj}_{\mathcal{R}}(q^{4k+1})},\quad q^{4k+4}= \underbrace{p(x_{0})\prod_{n=1}^{N+1}q^{4k+3}(x_{t_{n}}|x_{t_{n-1}})}_{\text{ proj}_{0}(\text{proj}_{\mathcal{M}}(q^{4k+3}))}.$ $$

Hence, the Bidirectional IMF procedure, in fact, alternates between the two projections of IMF ($\text{proj}_{\mathcal{MR}}$) during which the process ”became more optimal” (step towards optimality property) and two IPF projections (proj0 and proj1) during which the marginal fitting improves (step towards marginal matching property). Because of it, we have called this (bidirectional IMF) procedure, which starts from any starting process $q^{0}(x_{0},x_{\text{in}},x_{1})$ as Iterative Proportional Markovian Fitting (IPMF). We say that one IPMF step consists of these two projections of IMF ($\text{proj}_{\mathcal{MR}}$) and two projections of IPF. We hypothesize that this combined procedure should converge from any starting process $q^{0}(x_{0},x_{\text{in}},x_{1})$, unlike IPF and IMF procedures, which require a specific form of the starting process. We provide the analogical analysis for the continuous version in Appendix A.3. Next, we provide a theoretical analysis of the convergence of the IPMF procedure for 1D Gaussians.

### 3.2 Theoretical Analysis for 1D Gaussians

In this section, we analyze the case when $p_{0}=\mathcal{N}(\mu_{0},\sigma_{0}^{2})$ and $p_{1}=\mathcal{N}(\mu_{1},\sigma_{1}^{2})$ are 1D Gaussians and the initial process $q^{0}(x_{0},x_{in},x_{1})$ has Gaussian joint marginal $q^{0}(x_{0},x_{1})$ at times $t=0,1$. We prove that in this case $q^{4k}(x_{0},x_{1})$ converges to the solution $q^{*}(x_{0},x_{1})$ of the static SB problem (3).

To begin with, we do some preparations. We define a function $\Xi:[-1,1]\times\mathbb{R}_{+}\times\mathbb{R}_{+}\rightarrow\mathbb{R}$:

$$ $\Xi(\rho,\sigma,\sigma^{\prime})\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{ \rho}{\sigma\sigma^{\prime}(1-\rho^{2})}\in(-\infty,+\infty).$ $$

###### Lemma 3.1 (2D Gaussian plans as entropic optimal transport plans) .

Consider a 2-dimensional Gaussian distribution $q(x_{0},x_{1})$ with marginals $p=\mathcal{N}(\mu,\sigma^{2})$ and $p^{\prime}=\mathcal{N}(\mu^{\prime},(\sigma^{\prime})^{2})$ and correlation $\rho\in(-1,1)$ between its components:

$$ $q(x_{0},x_{1})=\mathcal{N}\left(\begin{pmatrix}\mu\\ \mu^{\prime}\end{pmatrix},\begin{pmatrix}\sigma^{2}&\rho\sigma\sigma^{\prime} \\ \rho\sigma\sigma^{\prime}&(\sigma^{\prime})^{2}\end{pmatrix}\right).$ $$

Let $\chi=\Xi(\rho,\sigma,\sigma^{\prime})$. Then $q$ is the unique minimizer of the following problem:

$$ $\min_{q^{\prime}\in\Pi(p,p^{\prime})}\bigg{\{}\int\frac{\chi}{2}\|x_{0}-x_{1} \|^{2}q^{\prime}(x_{0},x_{1})dx_{0}dx_{1}-H\big{(}q^{\prime}\big{)}\bigg{\}},$ (11) $$

where $H(q^{\prime})=-\int q^{\prime}(x_{0},x_{1})\log q^{\prime}(x_{0},x_{1})dx_{0}
dx_{1}$ is the differential entropy of a distribution.

Problem (11) is the optimal transport problem with the transport cost $\chi\cdot\|x_{0}-x_{1}\|^{2}/2$ and entropy regularization (with weight $1$), see (Cuturi, ; Genevay, ). Thus, our lemma states that any Gaussian distribution is a so-called entropic OT plan between its marginals for certain transport cost. In fact, now for every Gaussian distribution we can assign a value $\chi=\Xi(\rho,\sigma^{\prime},\sigma)$ explaining for which cost does this distribution solve the entropic transport problem. We will call this value the optimality coefficient. We emphasize that if the optimality coefficient is $\chi=\epsilon^{-1}$, then $q$ is the static SB (3) between its marginals for the prior $W^{\epsilon}$, recall (4). We note that our lemma for $\rho>0$ can be derived from (Janati et al., , \wasyparagraph3), (Bunne et al., , \wasyparagraph5); the non-trivial part is when $\rho\leq 0$. We provide a simple yet general proof for all $\rho\in(-1,1)$ using the entropic optimal transport benchmark constructor (Gushchin et al., , Thm. 3.2). Now we state our main result.

###### Theorem 3.2 (Quantitative convergence of IPMF for 1D Gaussians) .

Let $p_{0}=\mathcal{N}(\mu_{0},\sigma_{0}^{2})$ and $p_{1}=\mathcal{N}(\mu_{1},\sigma_{1}^{2})$ be 1D Gaussians. Assume that we run IPMF procedure in the continuous time or in discrete time with $N=1$ intermediate point, starting from some 2D Gaussian distribution(^1^11We assume that $q^{0}(x_{0})=p_{0}(x_{0})$, i.e., the initial process starts at $p_{0}$ at time $t=0$. This is reasonable as anyway after the first IPMF round the process will satisfy this property thanks to the IPF projections involved. We also consider $\rho_{0}\in(-1,1)$ only: if $\rho_{0}\in\{-1,1\}$, after the first IMF step it changes to $\in(-1,1)$.)

$$ $q^{0}(x_{0},x_{1})=\mathcal{N}\left(\begin{pmatrix}\mu_{0}\\ \nu\end{pmatrix},\begin{pmatrix}\sigma_{0}^{2}&\rho_{0}\sigma_{0}s_{0}\\ \rho_{0}\sigma_{0}s_{0}&s_{0}^{2}\end{pmatrix}\right)\text{ with }\rho_{0}\in( -1,1),s_{0}>0$ $$

and denote the joint distribution obtained after $k$ IPMF steps by

$$ $q^{4k}(x_{0},x_{1})=\mathcal{N}\left(\begin{pmatrix}\mu_{0}\\ \nu_{k}\end{pmatrix},\begin{pmatrix}\sigma_{0}^{2}&\rho_{k}\sigma_{0}s_{k}\\ \rho_{k}\sigma_{0}s_{k}&s_{k}^{2}\end{pmatrix}\right).$ $$

Denote $\chi_{k}\stackrel{{\scriptstyle\text{def}}}{{=}}\chi(\rho,\sigma_{0},s_{k})$. Then the following bounds hold true:

$$ $\displaystyle|s_{k}^{2}-\sigma_{1}^{2}|\leq\alpha^{2k}|s_{0}^{2}-\sigma_{1}^{2 }|,\qquad|\nu_{k}-\mu_{1}|\leq\alpha^{k}|\nu_{0}-\mu_{1}|,\qquad|\chi_{k}- \frac{1}{\epsilon}|\leq\beta^{2k}|\chi_{0}-\frac{1}{\epsilon}|.$ $$

where factors $\alpha$ and $\beta$ depend on IPMF type (discrete or continuous), initial parameters $s_{0},\nu_{0},\rho_{0}$, marginal distributions $p_{0},p_{1}$ and $\epsilon$. In particular, $\lim_{k\rightarrow\infty}\rho_{k}=\rho^{*}$, where $\rho^{*}$ is the correlation of the static SB solution $q^{*}$ between $p_{0},p_{1}$, namely, $\rho^{*}=(\sqrt{\sigma_{0}^{2}\sigma_{1}^{2}+\frac{\epsilon^{2}}{4}}-\frac{
\epsilon}{2})/(\sigma_{0}\sigma_{1})$.

Idea of our proof. To prove our theorem, we analyze formulas for IMF step in Gaussian case from (Peluchetti, ; Gushchin et al., ). We find that IPF steps keep the optimality coefficient $\chi_{k}$ while improving the marginal matching property of $q(x_{0},x_{1})$. In turn, IMF makes $\chi_{k}$ closer to $\frac{1}{\epsilon}$ while not affecting the marginals of $q^{4k}$. As a result, IPMF at each round improves both properties.

###### Lemma 3.1 (2D Gaussian plans as entropic optimal transport plans) .

###### Theorem 3.2 (Quantitative convergence of IPMF for 1D Gaussians) .

## 4 Experimental Illustrations

In this section, we empirically support that the IPMF procedure converges under a more general setting, specifically for any starting process, unlike IPF and IMF. We show in \wasyparagraph3.1 and Appendix A that the bidirectional IMF procedure and the introduced IPMF procedure differ only in the initial starting process. Since both practical implementations of continuous-time IMF (Shi et al., , Algorithm 1) and discrete-time IMF (Gushchin et al., , Algorithm 1) use the considered bidirectional version, we use practical algorithms introduced in these works, i.e., Diffusion Scrödinger Bridge Matching (DSBM) and Adversarial Schrödinger Bridge Matching (ASBM) respectively.

Experimental setups. We consider multivariate Gaussian distributions for which we have closed-form IPMF update formulas \wasyparagraph4.1, an illustrative 2D example, the Schrödinger Bridges Benchmark (Gushchin et al., ) and real life image data distributions, i.e., the colored MNIST dataset and the Celeba dataset (Liu et al., ). All technical details can be found in the Appendix D.

Starting processes. In our experiments, we focus on running the IPMF procedure from different initializations, which we call starting processes. We construct starting processes by considering different couplings $q^{0}(x_{0},x_{1})$ and using the Brownian Bridge process $W^{\epsilon}_{|x_{0},x_{1}}$. Thus, for each considered coupling $q^{0}(x_{0},x_{1})$ we construct starting process as $q^{0}(x_{0},x_{\text{in}},x_{1})=q^{0}(x_{0},x_{1})p^{W^{\epsilon}}(x_{\text{
in}}|x_{0},x_{1})$ for the discrete-time case and $T^{0}=\int W^{\epsilon}_{|x_{0},x_{1}}dq^{0}(x_{0},x_{1})$ for the continuous-time case (Appendix A). For all the experimental setups, we consider the starting processes induced by coupling $q^{0}(x_{0},x_{1})=p_{0}(x_{0})p_{1}(x_{1})$, which represent the IMF starting process (used in IMF procedure) and by $q^{0}(x_{0},x_{1})=p_{0}(x_{0})p^{W^{\epsilon}}(x_{1}|x_{0})$, which represent IPF starting process (used in the IPF procedure). We also consider a set of different couplings, which cannot be used either for starting processes of IMF or IPF procedure specifically for each setup to showcase that the IPMF procedure converges under more general assumptions. The results of DSBM and ASBM algorithms starting from different starting processes are denoted as (D/A)SBM-*name of coupling*, e.g., the results for the DSBM using the IMF starting process would be denoted as DSBM-IMF.

### 4.1 High Dimensional Gaussians

Figure: (a) Forward KL-divergence during IPMF steps.
Refer to caption: extracted/5898938/pics/high-dim-gauss/forward_gauss.png

We conduct experiments using analytical formulas for the Gaussian case for the discrete IMF from (Gushchin et al., , Theorem 3.8). We follow setup from (Gushchin et al., , Section 5.2) and consider Schrödinger Bridge problem with the dimensionality $D=16$ and $\epsilon=0.3$ for centered Gaussians ${p_{0}=\mathcal{N}(0,\Sigma_{0})}$ and ${p_{1}=\mathcal{N}(0,\Sigma_{1})}$. To construct $\Sigma_{0}$ and $\Sigma_{1}$, sample their eigenvectors from the uniform distribution on the unit sphere and sample their eigenvalues from the log uniform distribution on $[-\log 2,\log 2]$. We use the same $p_{0}$ and $p_{1}$ for all experiments.

We run the IPMF procedure for 100 IPMF steps (each IPMF step consists of 2 IPF projections and two Markovian-Reciprocal projections as stated in \wasyparagraph3.1). We present in Figure 1 the forward KL-divergence $\text{KL}(q^{4k}(x_{0},x_{1})||q^{*}(x_{0},x_{1}))$ and reversed KL-divergence $\text{KL}(q^{*}(x_{0},x_{1})||q^{4k}(x_{0},x_{1}))$ between $q^{4k}(x_{0},x_{1})$ from each IPMF step and the solution of static Schrödinger Bridge $q^{*}(x_{0},x_{1})$. We observe the fastest convergence for the IMF starting process, which is induced by the coupling $q^{0}(x_{0},x_{1})=p_{0}(x_{0})p_{1}(x_{1})$ following by the independent coupling $p_{0}\rightarrow p_{0}$ staring process $q^{0}(x_{0},x_{1})=p_{0}(x_{0})p_{0}(x_{1})$ and the worst speed of convergence for the IPF starting process with the coupling $q^{0}(x_{0},x_{1})=p_{0}(x_{0})p^{W^{\epsilon}}(x_{1}|x_{0})$. Interestingly, the independent $p_{0}\rightarrow p_{0}$ process, despite not being classic initialization (it was not used in IPF or IMF procedures), beats the IPF starting process and shows the same convergence speed as the IMF at the first $40$ IPMF steps.

### 4.2 Illustrative 2D example

Here, we consider the SB problem with $\epsilon=0.1$ with $p_{0}$ as the 2D Gaussian distribution and $p_{1}$ as the Swiss roll distribution. As previously mentioned, we train DSBM and ASBM algorithms using IMF and IPF starting processes. Additionaly we consider Independent $p_{0}\rightarrow p_{0}$ starting processes induced by the coupling $q^{0}(x_{0},x_{1})=p_{0}(x_{0})p_{0}(x_{1})$. We present starting processes and results in Figure 2. In all the cases, we observe convergence in the target distribution.

Figure: (a) ${x_{0}\sim p_{0}}$, ${x_{1}\sim p_{1}}$.
Refer to caption: extracted/5898938/pics/swiss_roll/swissroll_input_target.png

### 4.3 Evaluation on the SB Benchmark

We use the SB mixtures benchmark proposed by (Gushchin et al., ) with ground truth solution to the Scrödinger Bridge to test ASBM and DSBM with IMF, IPF and Independent $p_{0}\rightarrow p_{0}$ (i.e., induced by $q^{0}(x_{0},x_{1})=p_{0}(x_{0})p_{0}(x_{1})$ coupling), starting processes.

The benchmark provides continuous distribution pairs $p_{0}$, $p_{1}$ for dimensions $D\in\{2,16,64,128\}$ that have known SB solutions ($q^{*}$ for discrete setup and $T^{*}$ for continuous setup), for volatility $\epsilon\in\{0.1,1,10\}$. To evaluate the quality of the recovered SB solutions, we use the $\text{cB}\mathbb{W}_{2}^{2}\text{-UVP}$ metric as proposed by (Gushchin et al., ) and provide results in Table 1. In addition we study how all the approaches learn the target distribution in Appendix C.

As can be seen, DSBM and ASBM starting from all the processes at $\epsilon\in\{1,10\}$ yield quite similar results, but on the $\epsilon=0.1$ DSBM and ASBM with IPF and Independent $p_{0}\rightarrow p_{0}$ starting processes metric do experience a slight decrease.

**Table 1: Comparisons of $\text{cB}\mathbb{W}_{2}^{2}\text{-UVP}\downarrow$ (%) between the static SB solution $q^{*}(x_{0},x_{1})$ and the learned solution on the SB benchmark. The best metric is bolded. Results marked with $\dagger$ are taken from (Gushchin et al., ) and (Gushchin et al., ). The results of DSBM and ASBM algorithms starting from different starting processes are denoted as (D/A)SBM-*name of starting process***
|  |  | $\epsilon=0.1$ | $\epsilon=1$ | $\epsilon=10$ |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | Algorithm Type | $D\!=\!2$ | $D\!=\!16$ | $D\!=\!64$ | $D\!=\!128$ | $D\!=\!2$ | $D\!=\!16$ | $D\!=\!64$ | $D\!=\!128$ | $D\!=\!2$ | $D\!=\!16$ | $D\!=\!64$ | $D\!=\!128$ |
| Best algorithm on benchmark^† | Varies | $1.94$ | $13.67$ | $11.74$ | $11.4$ | $1.04$ | $9.08$ | $18.05$ | $15.23$ | $1.40$ | $1.27$ | $2.36$ | $\mathbf{1.31}$ |
| DSBM-IMF | IPMF | $1.21$ | $4.61$ | $9.81$ | $19.8$ | $0.68$ | $\mathbf{0.63}$ | $5.8$ | $29.5$ | $0.23$ | $5.45$ | $68.9$ | $362$ |
| DSBM-IPF | $2.55$ | $17.4$ | $15.85$ | $17.45$ | $0.29$ | $0.76$ | $4.05$ | $29.59$ | $0.35$ | $3.98$ | $83.2$ | $210$ |  |
| DSBM-Ind$(p_{0},p_{0})$ | $2.72$ | $11.7$ | $16.5$ | $17.02$ | $0.41$ | $0.92$ | $\mathbf{3.7}$ | $29$ | $0.16$ | $3.91$ | $101$ | $255$ |  |
| ASBM-IMF^† | $0.89$ | $8.2$ | $13.5$ | $53.7$ | $\mathbf{0.19}$ | $1.6$ | $5.8$ | $\mathbf{10.5}$ | $\mathbf{0.13}$ | $0.4$ | $1.9$ | $4.7$ |  |
| ASBM-IPF | $3.06$ | $14.37$ | $44.35$ | $32.5$ | $0.18$ | $1.68$ | $9.25$ | $20.47$ | $\mathbf{0.13}$ | $\mathbf{0.36}$ | $2.28$ | $4.97$ |  |
| ASBM-Ind$(p_{0},p_{0})$ | $3.99$ | $15.73$ | $39.3$ | $40.32$ | $0.18$ | $1.68$ | $6.16$ | $12.8$ | $\mathbf{0.13}$ | $0.38$ | $\mathbf{1.36}$ | $\mathbf{2.6}$ |  |
| $\text{SF}^{2}$M-Sink^† | $\mathbf{0.54}$ | $\mathbf{3.7}$ | $\mathbf{9.5}$ | $\mathbf{10.9}$ | $0.2$ | $1.1$ | $9$ | $23$ | $0.31$ | $4.9$ | $319$ | $819$ |  |

### 4.4 Unpaired Image-to-Image translation

To test our approach on real data, we consider two unpaired image-to-image translation setups: colorized 3 $\rightarrow$ colorized 2 digits from the MNIST dataset with 32$\times$32 resolution size and male$\rightarrow$female faces from the Celeba dataset with 64$\times$64 resolution size.

Figure: Figure 3: Inverted 7 starting process, i.e., reciprocal process with marginals $p_{0}$ and $p^{\text{inv}7}$, visualization.
Refer to caption: extracted/5898938/pics/cmnist/traj.png

Colored MNIST. We construct train and test sets by RGB colorization of MNIST digits from corresponding train and test sets of classes ”2” and ”3”. We train ASBM and DSBM algorithms starting from the IMF process. In addition, we test starting process induced by the independent coupling of the distribution of colored digits of class ”3” ($p_{0}$) and the distribution of colored digits of class ”7” with inverted RGB channels ($p^{\text{inv}7}(x_{1})$), we call this process Inverted 7, i.e., $q^{0}(x_{0},x_{1})=p_{0}(x_{0})p^{\text{inv}7}(x_{1})$. We visualize the Inverted 7 starting the process in Figure 3. Further technical details can be found in the Appendix D. We learn DSBM and ASBM on the train set of digits and visualize the translated test images in Figure 4.

We observe that both DSBM and ASBM algorithms starting from both IMF and Inverted 7 starting process fit the target distribution of colored MNIST digits of class ”2” and preserve the color of the input image during translation. This supports that the IPMF procedure converges to the same solution, which resembles the solution of the Schrödinger Bridge.

Figure: (a) ${x\sim p_{0}}$
Refer to caption: extracted/5898938/pics/cmnist/ipmf_cmnist_input.png

Figure: Figure 5: Convergence of models to target distribution. FID plotted as a function of IPMF iterations for all the presented setups.
Refer to caption: extracted/5898938/pics/celeba/IPMF_nfe_100_f_eps_1_Celeba64_IMF_FID.png

Celeba.
In this setup, we consider the variation of the IMF starting process called IMF-OT, where the starting process is induced by mini batch optimal transport coupling $q^{\text{OT}}(x_{0},x_{1})$ (Tong et al., ) given by the solution of the discrete optimal transport between two mini-batches, $\{x_{0,n}\}_{n=1}^{N}$ and $\{x_{1,n}\}_{n=1}^{N}$. As in the previous experiments, we also consider Independent $p_{0}\rightarrow p_{0}$ starting process. We use approximately the same number of parameters for DSBM and ASBM generator and use 10% of male and female images as the test set for evaluation, for other details see Appendix D. We provide qualitative results in Figure 6 and quantitative analysis through plotting FID as function of number of IPMF iterations in Figure 5.

We see from the qualitative results in Figure 6 that all trained models: 1) converge to the target distribution, 2) keep alignment between the features of the input images and generated images (e.g., the color of hair, background e.t.c.). However, it is worth noting that the samples generated by all models do differ, i.e., they generate different hair shapes and the presence or absence of glasses. We explain this as a different neural network optimization trajectory artifact. From the FID plot Figure5, we see that despite different starting processes and continuous or discrete time settings of DSBM and ASBM, all the models fit the target distribution equally well.

Figure: (a) ${x\sim p_{0}}$
Refer to caption: extracted/5898938/pics/celeba/ipmf_input_celeba.png

## 5 Discussion

Potential impact. The presented Iterative Proportional Markovian Fitting procedure shows a potential to overcome the error accumulation problem observed in distillation methods like rectified flows (Liu et al., ; ), which is used for the acceleration of the foundational image generation models, e.g. StableDiffusion 3 in (Esser et al., ). These distillation methods are based on the one-directional IMF procedure in the limit of Schrödinger Bridge hyperparameter $\epsilon=0$. However, the one-directional version is observed to accumulate errors and may even lead to the divergence (De Bortoli et al., , Appenix I). Furthermore, using the limiting case of $\epsilon=0$ makes it impossible to use the IPF procedure to restore marginals. The usage of a bidirectional version along with $\epsilon>0$ should both correct the marginals and make trajectories of diffusion more straight to accelerate the inference of diffusion models. We believe that considering such distillation techniques from the IPMF point of view may help to overcome the current limitations of these techniques.

Limitations.
While we show the proof of exponential convergence of the IPMF procedure for the 1D Gaussians in the continuous-time IMF and discrete-time IMF (with one inner point $N=1$), and present a wide set of experiments supporting this procedure, the proof of convergence of IPMF in the general case still remains a promising avenue for future work.

Reproducibility Statement. For all experiments presented, the full set of hyperparameters is shown either in Section 4 or in Appendix D. Derivations supporting theoretical claims Lemma 3.1 and Theorem 3.2 are included in the Appendix B.

## Acknowledgements

The work was supported by the Analytical center under the RF Government (subsidy agreement 000000D730321P5Q0002, Grant No. 70-2021-00145 02.11.2021).

## Appendix A Continuous-time Schrödinger Bridge Setup

### A.1 Schrödinger Bridge (SB) Problem in continuous-time

This section covers the continuous-time formulation of SB as its IPF and IMF procedures.
First, we introduce several new notations to better align the continuous version with the discrete-time version considered in the main text. Consider the Markovian process $T$ defined by the corresponding forward or backward (time-reversed) SDEs:

$$ $\displaystyle T:dx_{t}=v^{+}(x_{t},t)dt+\sqrt{\epsilon}dW_{t}^{+},\quad x_{0} \sim p_{0}(x_{0}),$ $\displaystyle T:dx_{t}=v^{-}(x_{t},t)dt+\sqrt{\epsilon}dW_{t}^{-},\quad x_{1} \sim p_{1}(x_{1}),$ $$

where we additionally denote by $W_{t}^{+}$ and $W_{t}^{-}$ the Wiener process in forward or backward time. We say $T_{|x_{0}}$ and $T_{|x_{1}}$ denotes the conditional process of $T$ fixing the marginals using delta functions, i.e. setting $p_{0}(x_{0})=\delta_{x_{0}}(x)$ and $p_{1}(x_{1})=\delta_{x_{1}}(x)$:

$$ $\displaystyle T_{|x_{0}}:dx_{t}=v^{+}(x_{t},t)dt+\sqrt{\epsilon}dW_{t}^{+}, \quad x_{0}\sim\delta_{x_{0}}(x),$ $\displaystyle T_{|x_{1}}:dx_{t}=v^{-}(x_{t},t)dt+\sqrt{\epsilon}dW_{t}^{-}, \quad x_{1}\sim\delta_{x_{1}}(x).$ $$

Moreover, we use $p(x_{0})T_{|x_{0}}$ to denote the stochastic process which starts by sampling ${x_{0}\sim p(x_{0})}$ and then moving this $x_{0}$ according the SDE given by $T_{|x_{0}}$, i.e., $p(x_{0})T_{|x_{0}}$ is short for the process $\int T_{|x_{0}}p(x_{0})dx_{0}$. Finally, we use the shortened notation of the process $T_{|0,1}(x_{0},x_{1})$ conditioned on its values at times 0 and 1, saying $p^{T}(x_{0},x_{1})T_{|0,1}(x_{0},x_{1})=\int T_{|0,1}(x_{0},x_{1})p^{T}(x_{0},
x_{1})dx_{0}dx_{1}$. This visually links the following equations with the discrete-time formulation.

Schrödinger Bridge problem. Considering the continuous case, the Schödinger Bridge problem is stated using continuous stochastic processes instead of one with predefined timesteps. Thus, the Schödinger Bridge problem finds the most likely in the sense of Kullback-Leibler divergence stochastic process $T$ with respect to prior Wiener process $W^{\epsilon}$, i.e.:

$$ $\min_{T\in\mathcal{F}(p_{0},p_{1})}\text{KL}(T||W^{\epsilon}),$ (12) $$

where $\mathcal{F}(p_{0},p_{1})\subset\mathcal{P}(C([0,1]),\mathbb{R}^{D})$ is the set of all stochastic processes pinned marginal distribution $p_{0}$ and $p_{1}$ at times 0 0 and $1$, respectively. The minimization problem (12) has a unique solution $T^{*}$ which can be represented as forward or backward diffusion (Léonard, ):

$$ $\displaystyle T^{*}:dx_{t}=v^{*+}(x_{t},t)dt+\sqrt{\epsilon}dW_{t}^{+},\quad x _{0}\sim p_{0}(x_{0}),$ $\displaystyle T^{*}:dx_{t}=v^{*-}(x_{t},t)dt+\sqrt{\epsilon}dW_{t}^{-},\quad x _{1}\sim p_{1}(x_{1}),$ $$

where $v^{*+}$ and $v^{*-}$ are the corresponding drift functions.

Static Schrödinger Bridge problem. As in discrete-time, Kullback-Leibler divergence in (12) could be decomposed as follows:

$$ $\text{KL}(T||W^{\epsilon})=\text{KL}(p^{T}(x_{0},x_{1})||p^{W^{\epsilon}}(x_{0 },x_{1}))+\int\text{KL}(T_{|x_{0},x_{1}}||W^{\epsilon}_{|x_{0},x_{1}})dp^{T}(x _{0},x_{1}).$ (13) $$

It has been proved (Léonard, ) that for the solution $T^{*}$ it’s conditional process is given by $T^{*}_{|x_{0},x_{1}}=W^{\epsilon}_{|x_{0},x_{1}}$. Thus, we can set $T_{|x_{0},x_{1}}=W^{\epsilon}_{|x_{0},x_{1}}$ zeroing the second term in (13) and minimize over processes with $T_{|x_{0},x_{1}}=W^{\epsilon}_{|x_{0},x_{1}}$. This leads to the equivalent Static formulation of the Schrodinger Bridge problem:

$$ $\min_{q\in\Pi(p_{0},p_{1})}\text{KL}(q(x_{0},x_{1})||p^{W^{\epsilon}}(x_{0},x_ {1})),$ (14) $$

where $\Pi(p_{0},p_{1})$ is the set of all joint distributions with marginals $p_{0}$ and $p_{1}$. Whether time is discrete or continuous, the decomposition of SB leads to the same static formulation, which, is closely related to Entropic OT as shown in (4).

### A.2 Iterative Proportional Fitting (IPF) for continuous-time

Following the main text, we describe the IPF procedure for continuous-time setup using stochastic processes. Likewise, IPF starts with setting $T^{0}=p_{0}(x_{0})W^{\epsilon}_{|x_{0}}$ and then it alternates between following projections:

$$ $\displaystyle T^{2k+1}=\text{proj}_{1}\left(p^{T^{2k}}(x_{1})T^{2k}_{|x_{1}} \right)\stackrel{{\scriptstyle\text{def}}}{{=}}p_{1}(x_{1})T^{2k}_{|x_{1}},$ (15) $\displaystyle T^{2k+2}=\text{proj}_{0}\left(p^{T^{2k+1}}(x_{0})T^{2k+1}_{|x_{0 }}\right)\stackrel{{\scriptstyle\text{def}}}{{=}}p_{0}(x_{0})T^{2k+1}_{|x_{0}}.$ (16) $$

As in the discrete-time case, these projections replace marginal distributions $p^{T}(x_{1})$ and $p^{T}(x_{0})$ in the processes $p^{T}(x_{1})T_{|x_{1}}$ and $p^{T}(x_{0})T_{|x_{0}}$ by $p_{1}(x_{1})$ and $p_{0}(x_{0})$ respectively.
Similarly to discrete-time formulation, the sequence of $T^{k}$ converges to the solution of the Schrödinger Bridge problem $T^{*}$ implicitly the reverse Kullback-Leibler divergence $\text{KL}(T^{k}||T^{*})$ between the current process $T^{k}$ and the solution to the SB problem $T*$. Additionally, it should be mentioned that projections are conducted by numerical approximation of forward and time-reversed conditional processes, $T_{|x_{0}}$ and $T_{|x_{1}}$, by learning their drifts via one of the methods: score matching (De Bortoli et al., ) or maximum likelihood estimation (Vargas et al., ).

### A.3 Iterative Markovian Fitting (IMF) for continuous-time

IMF introduces new projections that alternate between reciprocal and Markovian processes starting from any process $T^{0}$ pinned by $p_{0}$ and $p_{1}$, i.e. in $\mathcal{F}(p_{0},p_{1})$:

$$ $\displaystyle T^{2k+1}=\text{proj}_{\mathcal{R}}\left(T^{2k}\right)\stackrel{{ \scriptstyle\text{def}}}{{=}}p^{T^{2k}}(x_{0},x_{1})W^{\epsilon}_{|x_{0},x_{1}},$ (17) $\displaystyle T^{2k+2}=\text{proj}_{\mathcal{M}}\left(T^{2k+1}\right)\stackrel {{\scriptstyle\text{def}}}{{=}}\underbrace{p^{T^{2k+1}}(x_{0})T^{2k+1}_{M|x_{0 }}}_{\text{forward representation}}=\underbrace{p^{T^{2k+1}}(x_{1})T^{2k+1}_{M |x_{1}}}_{\text{backward representation}}.$ (18) $$

where we denote by $T_{M}$ the Markovian projections of the processes $T$, which can be represented as the forward or backward time diffusion as follows (Gushchin et al., , Section 2.1):

$$ $\displaystyle T_{M}:dx_{t}^{+}=v_{M}^{+}(x^{+}_{t},t)dt+\sqrt{\epsilon}dW_{t}^ {+},\quad x_{0}\sim p^{T}(x_{0}),\quad v_{M}^{+}(x^{+}_{t},t)=\int\frac{x_{1}- x^{+}_{t}}{1-t}p^{T}(x_{1}|x_{t})dx_{1},$ $\displaystyle T_{M}:dx_{t}^{-}=v_{M}^{-}(x^{-}_{t},t)dt+\sqrt{\epsilon}dW_{t}^ {-},\quad x_{1}\sim p^{T}(x_{1}),\quad v_{M}^{-}(x^{-}_{t},t)=\int\frac{x_{0}- x^{-}_{t}}{1-t}p^{T}(x_{0}|x_{t})dx_{0}.$ $$

This procedure converges to a unique solution, which is the Schrödinger bridge $T^{*}$ (Léonard, ). While reciprocal projection can be easily done by combining the joint distribution $p^{T}(x_{0},x_{1})$ of the process $T$ and Brownian bridge $W^{\epsilon}_{|x_{0},x_{1}}$, the Markovian projection is much more challenging and must be fitted via Bridge matching (Shi et al., ; ; Peluchetti, ).

Since the result of the Markovian projection can be represented (8) both by forward and backward representation, in practice, neural networks $v^{+}_{\theta}$ (forward parametrization) or $v^{-}_{\phi}$ (backward parametrization) are used to learn the corresponding drifts of the Markovian projections. In turn, starting distributions are set to be $p_{0}(x_{0})$ for forward parametrization and $p_{1}(x_{1})$ for the backward parametrization. So, this bidirectional procedure can be described as follows:

$$ $\displaystyle T^{4k+1}=\underbrace{p^{T^{4k}}(x_{0},x_{1})W^{\epsilon}_{|x_{0} ,x_{1}}}_{\text{proj}_{\mathcal{R}}(T^{4k})},\quad T^{4k+2}=\underbrace{p_{1}( x_{1}){T}^{4k+1}_{M|x_{1}}}_{\text{backward parametrization}},$ (19) $\displaystyle T^{4k+3}=\underbrace{p^{T^{4k+2}}(x_{0},x_{1})W^{\epsilon}_{|x_{ 0},x_{1}}}_{\text{proj}_{\mathcal{R}}(T^{4k+2})},\quad T^{4k+4}=\underbrace{p_ {0}(x_{0}){T}^{4k+3}_{M|x_{0}}}_{\text{forward parametrization}}.$ (20) $$

### A.4 Iterative Proportional Markovian Fitting (IPMF) for continuous-time

Here, we analyze the continuous version of the heuristical bidirectional IMF. First, we recall, that the IPF projections proj${}_{0}(T)$ and proj${}_{1}(T)$ given by (15) and (16) of the Markovian process $T$ is just change the starting distribution from $p^{T}(x_{0})$ to $p_{0}(x_{0})$ and $p^{T}(x_{1})$ to $p_{1}(x_{1})$.

Now we note that the process $T^{4k+2}$ in (19) is obtained by using a combination of Markovian projection projM given by (18) in forward parametrization and IPF projection proj1 given by (16):

$$ $\displaystyle T^{4k+2}=p_{1}(x_{1})T^{4k+1}_{M|x_{1}}=\underbrace{\text{proj}_ {1}\big{(}p^{T^{4k+1}}(x_{1})T^{4k+1}_{M|x_{1}}\big{)}}_{\text{proj}_{1}(\text {proj}_{\mathcal{M}}(T^{4k+1}))}.$ $$

In turn, the process $T^{4k+4}$ in (20) is obtained by using a combination of Markovian projection projM given by (18) in backward parametrization and IPF projection proj0 given by (15):

$$ $\displaystyle T^{4k+3}=p_{0}(x_{0})T^{4k+3}_{M|x_{0}}=\underbrace{\text{proj}_ {0}\big{(}p^{T^{4k+3}}(x_{0})T^{4k+3}_{M|x_{0}}\big{)}}_{\text{proj}_{0}(\text {proj}_{\mathcal{M}}(T^{4k+3}))}.$ $$

Combining these facts we can rewrite bidirectional IMF in the following manner:

$$ Iterative Proportional Markovian Fitting (Conitnious time setting) $\displaystyle T^{4k+1}=\underbrace{p^{T^{4k}}(x_{0},x_{1})W^{\epsilon}_{|x_{0} ,x_{1}}}_{\text{proj}_{\mathcal{R}}(T^{4k})},\quad T^{4k+2}=\underbrace{p_{1}( x_{1})T^{4k+1}_{M|x_{1}}}_{\text{proj}_{1}\left(\text{proj}_{\mathcal{M}}(T^{4 k+1})\right)}$ (21) $\displaystyle T^{4k+3}=\underbrace{p^{T^{4k+2}}(x_{0},x_{1})W^{\epsilon}_{|x_{ 0},x_{1}}}_{\text{proj}_{\mathcal{R}}(T^{4k+2})},\quad T^{4k+4}=\underbrace{p_ {0}(x_{0}){T}^{4k+3}_{M|x_{0}}}_{\text{proj}_{0}\left(\text{proj}_{\mathcal{M} }(T^{4k+3})\right)}.$ (22) $$

Thus, we obtain the analog of the discrete-time IPMF procedure, which concludes our description of the continuous setups.

## Appendix B Complete Theoretical Analysis for 1D Gaussians

In this section, we study behavior of IPMF with parameter $\varepsilon$ between 1D Gaussians $p_{0}=\mathcal{N}(\mu_{0},\sigma_{0}^{2})$ and $p_{1}=\mathcal{N}(\mu_{1},\sigma_{1}^{2})$. For convenience, we use notation $\varepsilon_{*}$ instead of $\varepsilon$. We prove that the process parameters with each step exponentially converge to desired values $\mu_{0},\mu_{1},\sigma_{0},\sigma_{1},\varepsilon_{*}$.

The section is organized as follows:

1) In Appendix B.1, we reveal the connection between arbitrary 2-dimensional Gaussian distribution and solution of entropic OT problem with specific transport cost, i.e., we prove our Lemma 3.1.

2) In Appendix B.2, we study the effect of IMF step on the current process. We show that after IMF, marginals remain the same, while the new correlation becomes close to the correlation of the static $\varepsilon_{*}$-EOT solution between marginals.

3) In Appendix B.3, we study the effect of IPF steps on the current process. We show that during these steps, the marginals become close to $p_{0}$ and $p_{1}$, while the optimality value does not change.

4) Finally, in Appendix B.4, we provide the complete proof of our main Theorem 3.2.

### B.1 2D Gaussian plans as entropic optimal transport plans

###### Proof of Lemma 3.1 .

The solution of $1$-entropic optimal transport with the transport cost $c(x_{0},x_{1})=\frac{\chi}{2}\|x_{0}-x_{1}\|^{2}$ is

$$ $q_{c}(x_{0},x_{1})=\operatorname*{arg\,min}_{q^{\prime}\in\Pi(p,p^{\prime})} \bigg{\{}\int\frac{\chi}{2}\|x_{0}-x_{1}\|^{2}q^{\prime}(x_{0},x_{1})dx_{0}dx_ {1}-H\big{(}q^{\prime}\big{)}\bigg{\}}.$ $$

We note that we can add any functions $f(x_{0})$ and $g(x_{1})$ depending only on $x_{0}$ or $x_{1}$, respectively, to the cost function, and the solution does not change:

$$ $\displaystyle\int\left[f(x_{0})+\frac{\chi}{2}\|x_{0}-x_{1}\|^{2}+g(x_{1}) \right]q^{\prime}(x_{0},x_{1})dx_{0}dx_{1}$ $\displaystyle=$ $\displaystyle\int f(x_{0})p(x_{0})dx_{0}+\int g(x_{1})p^{\prime}(x_{1})dx_{1}+ \int\frac{\chi}{2}\|x_{0}-x_{1}\|^{2}q^{\prime}(x_{0},x_{1})dx_{0}dx_{1}$ $\displaystyle\sim$ $\displaystyle\int\frac{\chi}{2}\|x_{0}-x_{1}\|^{2}q^{\prime}(x_{0},x_{1})dx_{0 }dx_{1}.$ $$

Thus, for any $\chi\in(-\infty,+\infty)$, we can rearrange the cost term $c(x_{0},x_{1})$ so that it becomes lower-bounded. For positive $\chi$, it already holds true. For negative $\chi$, we use

$$ $\displaystyle\tilde{c}(x_{0},x_{1})=-\chi\|x_{0}\|^{2}+\frac{\chi}{2}\|x_{0}-x _{1}\|^{2}-\chi\|x_{0}\|^{2}=-\frac{\chi}{2}\|x_{0}+x_{1}\|^{2}\geq 0,$ $$

where $\tilde{c}(x_{0},x_{1})$ is a lower bounded function.

Following (Gushchin et al., , Theorem 3.2), the conditional distribution $q_{c}(x_{1}|x_{0})$ with the lower bounded cost function $c$ can be expressed as:

$$ $\displaystyle q_{c}(x_{1}|x_{0})$ $\displaystyle\propto$ $\displaystyle\exp\left(-c(x_{0},x_{1})+f_{c}(x_{1})\right)=\exp\left(-\frac{ \chi}{2}\|x_{0}-x_{1}\|^{2}+f_{c}(x_{1})\right),$ $$

where function $f_{c}(x_{1})$ depends only on $x_{1}$. Moreover, we can simplify this distribution even more

$$ $\displaystyle q_{c}(x_{1}|x_{0})$ $\displaystyle\propto$ $\displaystyle Z_{x_{0}}Z_{x_{1}}\exp\left(\chi x_{0}x_{1}\right),$ (23) $$

where factors $Z_{x_{0}}$ and $Z_{x_{1}}$ depend only on $x_{0}$ and $x_{1}$, respectively.

Meanwhile, the conditional distribution of $q(x_{1}|x_{0})$ has a closed form, namely,

$$ $\displaystyle q(x_{1}|x_{0})$ $\displaystyle=$ $\displaystyle\mathcal{N}\left(\mu^{\prime}+\frac{\sigma^{\prime}}{\sigma}(x_{0 }-\mu),(1-\rho^{2})(\sigma^{\prime})^{2}\right)\propto\exp\left(\frac{-(x_{1}- \mu^{\prime}-\frac{\sigma^{\prime}\rho}{\sigma}(x_{0}-\mu))^{2}}{2(1-\rho^{2}) (\sigma^{\prime})^{2}}\right)$ (24) $\displaystyle=$ $\displaystyle Z_{x_{0}}Z_{x_{1}}\exp\left(\frac{\sigma^{\prime}\rho}{\sigma(1- \rho^{2})(\sigma^{\prime})^{2}}x_{0}x_{1}\right)=Z_{x_{0}}Z_{x_{1}}\exp\left( \frac{\rho}{\sigma\sigma^{\prime}(1-\rho^{2})}x_{0}x_{1}\right),$ $$

where factors $Z_{x_{0}}$ and $Z_{x_{1}}$ depend only on $x_{0}$ and $x_{1}$, respectively. Equating (23) and (24) terms which depend on $x_{0}$ and $x_{1}$ simultaneously, we obtain the required function

$$ $\Xi(\rho,\sigma,\sigma^{\prime})\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{ \rho}{\sigma\sigma^{\prime}(1-\rho^{2})}=\chi\in(-\infty,+\infty).$ (25) $$

The function $\Xi$ is monotonously increasing w.r.t. $\rho\in(-1,1)$ and, thus, invertible, i.e., there exists a function $P:(-\infty,+\infty)\times\mathbb{R}_{+}\times\mathbb{R}_{+}\to(-1,1)$ such that

$$ $P(\chi,\sigma,\sigma^{\prime})=\Xi^{-1}(\chi,\sigma,\sigma^{\prime})=\frac{ \sqrt{\chi^{2}\sigma^{2}(\sigma^{\prime})^{2}+1/4}-1/2}{\chi\sigma\sigma^{ \prime}}.$ (26) $$

The inverse function $P$ is calculated via solving quadratic equation w.r.t. $\rho.$
∎

###### Proof of Lemma 3.1 .

### B.2 IMF step analysis

In our paper, we consider both discrete and continuous IMF. By construction, IMF step does change marginals of the process it works with. However, for both continuous and discrete IMF, the new correlation converges to the correlation of the $\varepsilon_{*}$-EOT between marginals.

###### Lemma B.1 (Correlation improvement after (D)IMF step) .

Consider a 2-dimensional Gaussian distribution with marginals $p=\mathcal{N}(\mu,\sigma^{2})$ and $p^{\prime}=\mathcal{N}(\mu^{\prime},(\sigma^{\prime})^{2})$ and correlation $\rho\in(-1,1)$ between its components. After continuous IMF or DIMF with single time point $t$, we obtain correlation $\rho_{new}$. The distance between $\rho_{new}$ and EOT correlation $\rho_{*}=P(\nicefrac{{1}}{{\varepsilon_{*}}},\sigma,\sigma^{\prime})$ decreases as:

$$ $\displaystyle|\rho_{new}-\rho_{*}|$ $\displaystyle\leq$ $\displaystyle\gamma\cdot|\rho-\rho_{*}|,$ $$

where factor $\gamma$ for continuous and discrete IMF (with $N=1)$ is, respectively,

$$ $\displaystyle\gamma_{c}(\sigma,\sigma^{\prime})$ $\displaystyle=$ $\displaystyle\frac{\sqrt{\sigma^{2}(\sigma^{\prime})^{2}+\varepsilon_{*}^{2}/4 }-\varepsilon_{*}/2}{\sigma\sigma^{\prime}},$ (27) $\displaystyle\gamma_{d}(\sigma,\sigma^{\prime},t)$ $\displaystyle=$ $\displaystyle\frac{1}{1+\frac{t^{2}(1-t)^{2}\sigma^{2}(\sigma^{\prime})^{2}+t( 1-t)(t^{2}(\sigma^{\prime})^{2}+(1-t)^{2}\sigma^{2})\varepsilon+t^{2}(1-t)^{2} \varepsilon^{2}_{*}}{(1-t)^{2}((1-t)\sigma^{2}+t\sigma\sigma^{\prime})^{2}+t^{ 2}(t(\sigma^{\prime})^{2}+(1-t)\sigma\sigma^{\prime})^{2}+t(1-t)((1-t)\sigma+t \sigma^{\prime})^{2}\varepsilon_{*}}}.$ (28) $$

###### Proof.

Continuous case.
Following (Peluchetti, , Eq. 42), we have the formula for $\rho_{new}$:

$$ $\displaystyle\rho_{new}(\rho)$ $\displaystyle=$ $\displaystyle\exp\left\{-\varepsilon_{*}\frac{\tanh^{-1}\left(\frac{c_{1}}{c_{ 3}}\right)+\tanh^{-1}\left(\frac{c_{2}}{c_{3}}\right)}{c_{3}}\right\}>0,$ (29) $\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle\varepsilon_{*}+2(\sigma^{\prime})^{2}(\rho\sigma^{2}-(\sigma^{ \prime})^{2}),c_{3}=\sqrt{(\varepsilon_{*}+2(\rho+1)\sigma^{2}(\sigma^{\prime} )^{2})(\varepsilon_{*}+2(\rho-1)\sigma^{2}(\sigma^{\prime})^{2})},$ $\displaystyle c_{2}$ $\displaystyle=$ $\displaystyle\varepsilon_{*}+2\sigma^{2}(\rho(\sigma^{\prime})^{2}-\sigma^{2}).$ $$

The map $\rho_{new}(\rho)$ is contraction over $\rho\in[-1,1]$ with the contraction coefficient $\gamma_{c}(\sigma,\sigma^{\prime})\stackrel{{\scriptstyle\text{def}}}{{=}}
\frac{\sqrt{\sigma^{2}(\sigma^{\prime})^{2}+\varepsilon_{*}^{2}/4}-\varepsilon
_{*}/2}{\sigma\sigma^{\prime}}<1.$ The unique fixed point of such map is $\rho_{*}=P(\nicefrac{{1}}{{\varepsilon_{*}}},\sigma_{0},\sigma_{1})$, since IMF does not change $\varepsilon_{*}$-EOT solution. Hence, we derive a bound

$$ $|\rho_{new}(\rho)-\rho_{*}|=|\rho_{new}(\rho)-\rho_{new}(\rho_{*})|\leq\gamma_ {c}(\sigma,\sigma^{\prime})|\rho-\rho_{*}|.$ $$

Discrete case ($N=1$). In this case, we use notations from (Gushchin et al., ), namely, we denote covariance matrix $\begin{pmatrix}\Sigma_{0}&\Sigma\\
\Sigma&\Sigma_{1}\end{pmatrix}\stackrel{{\scriptstyle\text{def}}}{{=}}\begin{
pmatrix}\sigma^{2}&\rho\sigma\sigma^{\prime}\\
\rho\sigma\sigma^{\prime}&(\sigma^{\prime})^{2}\end{pmatrix}$.

The general formulas of DIMF step are given for time points $0=t_{0}<t_{1}<\dots<t_{N}<t_{N+1}=1$. Following (Gushchin et al., ), we have an explicit formula for reciprocal step. For any $1\leq k\leq N$, we have joint covariance between time moments

$$ $\displaystyle\Sigma_{t_{k},t_{k}}$ $\displaystyle=$ $\displaystyle(1-t_{k})^{2}\Sigma_{0}+2t_{k}(1-t_{k})\Sigma+t_{k}^{2}\Sigma_{1} +t_{k}(1-t_{k})\varepsilon_{*},$ $\displaystyle\Sigma_{t_{k+1},t_{k}}$ $\displaystyle=$ $\displaystyle(1-t_{k})(1-t_{k+1})\Sigma_{0}+[(1-t_{k})t_{k+1}+(1-t_{k+1})t_{k} )]\Sigma+t_{k}t_{k+1}\Sigma_{1}+t_{k}(1-t_{k+1})\varepsilon_{*},$ $\displaystyle\Sigma_{t_{1},0}$ $\displaystyle=$ $\displaystyle(1-t_{1})\Sigma_{0}+t_{1}\Sigma,$ $\displaystyle\Sigma_{1,t_{N}}$ $\displaystyle=$ $\displaystyle t_{N}\Sigma_{1}+(1-t_{N})\Sigma.$ $$

Matrices $\Sigma_{t_{k+1},t_{k}}$ and $\Sigma_{t_{k},t_{k}}$ depend on $\Sigma$. For markovian step, we write down an analytical formula for the new covariance $\Sigma_{new}$ between marginals:

$$ $\displaystyle f(\Sigma)\stackrel{{\scriptstyle\text{def}}}{{=}}\Sigma_{new}(\Sigma)$ $\displaystyle=$ $\displaystyle\prod_{k=0}^{N}\left(\Sigma_{t_{k+1},t_{k}}\cdot\Sigma_{t_{k},t_{ k}}^{-1}\right)\cdot\Sigma_{0},$ $$

The derivative of $f^{\prime}(\Sigma)$ is as follows:

$$ $\displaystyle f^{\prime}(\Sigma)$ $\displaystyle=$ $\displaystyle\left[\frac{(1-t_{N})}{\Sigma_{1,t_{N}}}-\frac{2t_{N}(1-t_{N})}{ \Sigma_{t_{N},t_{N}}}+\frac{t_{1}}{\Sigma_{t_{1},0}}\right]\cdot f(\Sigma)$ (30) $\displaystyle+$ $\displaystyle\sum_{k=1}^{N-1}\left[\frac{[(1-t_{k})t_{k+1}+t_{k}(1-t_{k+1})]}{ \Sigma_{t_{k+1},t_{k}}}-\frac{2t_{k}(1-t_{k})\Sigma}{\Sigma_{t_{k},t_{k}}} \right]\cdot f(\Sigma)$ $\displaystyle=$ $\displaystyle\sum_{k=0}^{N}\left[-\frac{t_{k+1}(1-t_{k+1})}{\Sigma_{t_{k+1},t_ {k+1}}}+\frac{t_{k+1}(1-t_{k})+t_{k}(1-t_{k+1})}{\Sigma_{t_{k+1},t_{k}}}-\frac {t_{k}(1-t_{k})}{\Sigma_{t_{k},t_{k}}}\right]\cdot f(\Sigma).$ $$

In the case of single point $t=t_{1}$ ($N=1$), we prove that the function $f(\Sigma)$ is a contraction map. The sufficient condition for the map to be contraction is to have derivative’s norm bounded by $\gamma_{d}<1$.

Firstly, we can write down the simplified formula $\rho_{new}(\rho)$ in our original notations:

$$ $\displaystyle\rho_{new}(\rho)=\frac{((1-t)\sigma+t\rho\sigma^{\prime})(t\sigma ^{\prime}+(1-t)\rho\sigma)}{(1-t)\sigma^{2}+2t(1-t)\rho\sigma\sigma^{\prime}+t ^{2}(\sigma^{\prime})^{2}+t(1-t)\varepsilon_{*}}.$ (31) $$

Next, we simplify derivative (30):

$$ $\displaystyle\Sigma_{0,t}$ $\displaystyle=$ $\displaystyle(1-t)\cdot\Sigma_{0}+t\cdot\Sigma,$ $\displaystyle\Sigma_{t,1}$ $\displaystyle=$ $\displaystyle t\cdot\Sigma_{1}+(1-t)\cdot\Sigma,$ $\displaystyle\Sigma_{t,t}$ $\displaystyle=$ $\displaystyle(1-t)^{2}\cdot\Sigma_{0}+2(1-t)t\cdot\Sigma+t^{2}\cdot\Sigma_{1}+ t(1-t)\varepsilon_{0}=(1-t)\cdot\Sigma_{0,t}+t\cdot\Sigma_{t,1}+t(1-t)\varepsilon,$ $\displaystyle f^{\prime}(\Sigma)$ $\displaystyle=$ $\displaystyle\frac{(1-t)\Sigma_{0,t}}{\Sigma_{t,t}}+\frac{t\Sigma_{t,1}}{ \Sigma_{t,t}}-2\cdot\frac{t\Sigma_{t,1}\cdot(1-t)\Sigma_{0,t}}{\Sigma_{t,t} \cdot\Sigma_{t,t}}.$ $$

We define new variables $A\stackrel{{\scriptstyle\text{def}}}{{=}}(1-t)\Sigma_{0,t},B\stackrel{{
\scriptstyle\text{def}}}{{=}}t\Sigma_{t,1}$ and $\tilde{\varepsilon}_{*}=t(1-t)\varepsilon_{*}$. We note that while $\Sigma\in[-\sqrt{\Sigma_{0}\Sigma_{1}},\sqrt{\Sigma_{0}\Sigma_{1}}]$ the value $A+B=(1-t)^{2}\cdot\Sigma_{0}+2(1-t)t\cdot\Sigma+t^{2}\cdot\Sigma_{1}\geq 0$.
Then, we restate $f^{\prime}$ as:

$$ $\displaystyle f^{\prime}$ $\displaystyle=$ $\displaystyle\frac{A}{A+B+\tilde{\varepsilon}_{*}}+\frac{B}{A+B+\tilde{ \varepsilon}_{*}}-\frac{2AB}{(A+B+\tilde{\varepsilon}_{*})^{2}}=\frac{(A+B)(A+ B+\tilde{\varepsilon}_{*})-2AB}{(A+B+\tilde{\varepsilon}_{*})^{2}}$ (32) $\displaystyle=$ $\displaystyle\frac{A^{2}+B^{2}+(A+B)\tilde{\varepsilon}_{*}}{(A+B+\tilde{ \varepsilon}_{*})^{2}}$ $\displaystyle=$ $\displaystyle\frac{A^{2}+B^{2}+(A+B)\tilde{\varepsilon}_{*}}{A^{2}+2AB+B^{2}+2 (A+B)\tilde{\varepsilon}_{*}+\tilde{\varepsilon}_{*}^{2}}=\frac{1}{1+\frac{2AB +(A+B)\tilde{\varepsilon}_{*}+\tilde{\varepsilon}_{*}^{2}}{A^{2}+B^{2}+(A+B) \tilde{\varepsilon}_{*}}}.$ (33) $$

We note that all terms in (32) are greater than 0 0:

$$ $0<f^{\prime}(\Sigma),\quad\Sigma\in[-\sqrt{\Sigma_{0}\Sigma_{1}},\sqrt{\Sigma_ {0}\Sigma_{1}}].$ (34) $$

In negative segment $\Sigma\in[-\sqrt{\Sigma_{0}\Sigma_{1}},0]$, the derivative $f^{\prime}$ is greater than in positive segment $[0,\sqrt{\Sigma_{0}\Sigma_{1}}]$, and edge value $f(-\sqrt{\Sigma_{0}\Sigma_{1}})>-\sqrt{\Sigma_{0}\Sigma_{1}}$. Thus, in negative segment, the convergence to the fixed point $\rho_{*}\sqrt{\Sigma_{0}\Sigma_{1}}>0$ is faster, than in positive segment.

For $\Sigma\in[0,\sqrt{\Sigma_{0}\Sigma_{1}}]$, we can bound the fraction in denominator of (33) by taking its numerator’s minimum at $\Sigma=0$ and its denominator’s maximum at $\Sigma=\sqrt{\Sigma_{0}\Sigma_{1}}$, i.e,

$$ $0<f^{\prime}\leq\gamma_{d}(\Sigma_{0},\Sigma_{1},t)<1,$ $$

$$ $\gamma_{d}(\Sigma_{0},\Sigma_{1},t)=\frac{1}{1+\frac{t^{2}(1-t)^{2}\Sigma_{0} \Sigma_{1}+t(1-t)(t^{2}\Sigma_{1}+(1-t)^{2}\Sigma_{0})\varepsilon+t^{2}(1-t)^{ 2}\varepsilon^{2}_{*}}{(1-t)^{2}((1-t)\Sigma_{0}+t\sqrt{\Sigma_{0}\Sigma_{1}}) ^{2}+t^{2}(t\Sigma_{1}+(1-t)\sqrt{\Sigma_{0}\Sigma_{1}})^{2}+t(1-t)((1-t)\sqrt {\Sigma_{0}}+t\sqrt{\Sigma_{1}})^{2}\varepsilon_{*}}}.$ $$

We note that $\gamma_{d}(\Sigma_{0},\Sigma_{1},t)$ is increasing function w.r.t. $\Sigma_{0},\Sigma_{1}.$

If we put into the function $f$ argument $\Sigma_{*}=\rho_{*}\sqrt{\Sigma_{0}\Sigma_{1}}$ corresponding to the $\varepsilon_{*}$-EOT correlation, DIMF does not change it. Hence, $\Sigma_{*}$ is the fixed point of $f(\Sigma)$, and we have

$$ $\displaystyle|\Sigma_{new}-\Sigma_{*}|=|f(\Sigma)-f(\Sigma_{*})|\leq\gamma_{d} (\Sigma_{0},\Sigma_{1},t)|\Sigma-\Sigma_{*}|.$ $$

Dividing both sides by $\sqrt{\Sigma_{0}\Sigma_{1}}$, we get

$$ $\displaystyle|\rho_{new}-\rho_{*}|\leq\gamma_{d}(\Sigma_{0},\Sigma_{1},t)|\rho -\rho_{*}|.$ $$

∎

###### Lemma B.2 ( $\chi$ improvement after (D)IMF step) .

Consider a 2-dimensional Gaussian distribution with marginals $p=\mathcal{N}(\mu,\sigma^{2})$ and $p^{\prime}=\mathcal{N}(\mu^{\prime},(\sigma^{\prime})^{2})$ and correlation $\rho\in(-1,1)$ between its components. After continuous IMF or DIMF with a single time point $t$, we obtain new correlation $\rho_{new}$, such that $|\rho_{new}-\rho_{*}|\leq\gamma|\rho-\rho_{*}|$ where $\rho_{*}=P(\nicefrac{{1}}{{\varepsilon_{*}}},\sigma,\sigma^{\prime})$ and $\gamma<1$ is from (27) for IMF and from (28) for DIMF. We have bound in terms of $\chi=\Xi(\rho,\sigma,\sigma^{\prime})$ and $\chi_{new}=\Xi(\rho_{new},\sigma,\sigma^{\prime})$:

$$ $\displaystyle|\chi_{new}-\nicefrac{{1}}{{\varepsilon_{*}}}|$ $\displaystyle\leq$ $\displaystyle l(\rho,\rho_{*},\gamma)\cdot|\chi-\nicefrac{{1}}{{\varepsilon_{* }}}|,$ (35) $\displaystyle l(\rho,\rho_{*},\gamma)$ $\displaystyle=$ $\displaystyle\left[1-(1-\gamma)\frac{(1-\max\{\rho_{*},|\rho|\}^{2})^{2}}{1+ \max\{\rho_{*},|\rho|\}^{2}}\right]<1.$ $$

###### Proof.

Monotone. The function $\rho_{new}$ from (29) for continuous IMF and from (31) for DIMF is monotonously increasing on $(-1,1)$. For continuous IMF, with the growth of $\rho$, $c_{3}$ grows faster than $c_{1}$ or $c_{2}$, hence, $\frac{c_{1}}{c_{3}},\frac{c_{2}}{c_{3}}$ and $\tanh^{-1}\left(\frac{c_{1}}{c_{3}}\right),\tanh^{-1}\left(\frac{c_{2}}{c_{3}}\right)$ decrease. Thus, the power in the exponent of $\rho_{new}$ and $\rho_{new}$ itself increase. For DIMF, the derivative of $\rho_{new}$ greater than zero based on (34).

The monotone means that the value $\rho_{new}$ always remains from the same side from $\rho$:

$$ $\begin{cases}\rho>\rho_{*}&\Longrightarrow\rho_{new}(\rho)>\rho_{*}=\rho_{new} (\rho_{*}),\\ \rho\leq\rho_{*}&\Longrightarrow\rho_{new}(\rho)\leq\rho_{*},\end{cases}$ (36) $$

The same inequalities hold true for $\chi=\Xi(\rho,\sigma,\sigma^{\prime}),\chi_{new}=\Xi(\rho_{new},\sigma,\sigma^
{\prime})$ and $\chi_{*}=\nicefrac{{1}}{{\varepsilon_{*}}}$ as well: if $\chi<\chi_{*}$, then $\chi_{new}<\chi_{*}$ and vice versa, since $\Xi(\rho,\sigma,\sigma^{\prime})$ is monotonously increasing w.r.t. $\rho$.

$\Xi$ Properties. In this proof, we omit arguments $\sigma,\sigma^{\prime}$ of $P(\chi,\sigma,\sigma^{\prime})$ and $\Xi(\rho,\sigma,\sigma^{\prime})$, because they do not change during IMF step. The second derivative of the function $\Xi(\rho)$ is

$$ $\frac{d^{2}\Xi}{d\rho^{2}}(\rho)=\frac{2\rho(3+\rho^{2})}{\sigma\sigma^{\prime }(1-\rho^{2})^{3}}.$ $$

Hence, we have $\frac{d^{2}\Xi}{d\rho^{2}}(\rho)\leq 0$ for $\rho\in(-1,0]$ and $\frac{d^{2}\Xi}{d\rho^{2}}(\rho)\geq 0$ for $\rho\in[0,1)$. It means that the function $\Xi(\rho)$ is concave on $(-1,0]$ and convex on $[0,1)$.

The function $\Xi(\rho)$ is monotonously increasing w.r.t. $\rho$, thus, decreasing of the radius $h\stackrel{{\scriptstyle\text{def}}}{{=}}|\rho-\rho_{*}|$ around $\rho_{*}$ causes the decreasing of $|\chi-\chi_{*}|$ around $\chi_{*}$. We consider two cases: $\chi>\chi_{*}$ and $\chi<\chi_{*}$.

Case $\chi>\chi_{*}$. We have $\rho=\rho_{*}+h,\chi=\Xi(\rho_{*}+h)=\Xi(\rho)$ and $\Xi(\rho_{*}+\gamma h)\geq\chi_{new}$. We compare the difference using convexity on $[0,1)$:

$$ $\displaystyle\chi-\chi_{new}$ $\displaystyle\geq$ $\displaystyle\Xi(\rho_{*}+h)-\Xi(\rho_{*}+\gamma h)\geq(\rho_{*}+h-(\rho_{*}-h \gamma))\cdot\frac{d\Xi}{d\rho}(\rho_{*}+\gamma h)$ $\displaystyle=$ $\displaystyle(1-\gamma)h\cdot\frac{d\Xi}{d\rho}(\rho_{*}+\gamma h).$ $$

Since the derivative of $\Xi$ is always positive, we continue the bound:

$$ $\displaystyle\Xi(\rho_{*}+h)-\Xi(\rho_{*}+\gamma h)\geq\min_{\rho^{\prime}\in[ \rho_{*},\rho_{*}+h]}\left|\frac{d\Xi}{d\rho}(\rho^{\prime})\right|(1-\gamma)| \rho-\rho_{*}|.$ $$

Next, we use Lipschitz property of $\Xi$, i.e.,

$$ $|\chi-\chi_{*}|=|\Xi(\rho)-\Xi(\rho_{*})|\leq\max_{\rho^{\prime}\in[\rho_{*}, \rho_{*}+h]}\left|\frac{d\Xi}{d\rho}(\rho)\right||\rho-\rho_{*}|,$ $$

and combine it with the previous bound

$$ $\displaystyle\chi-\chi_{new}\geq\Xi(\rho_{*}+h)-\Xi(\rho_{*}+\gamma h)$ $\displaystyle\geq$ $\displaystyle\frac{\min\limits_{\rho^{\prime}\in[\rho_{*},\rho]}\left|\frac{d \Xi}{d\rho}(\rho^{\prime})\right|}{\max\limits_{\rho^{\prime}\in[\rho_{*},\rho ]}|\frac{d\Xi}{d\rho}(\rho^{\prime})|}(1-\gamma)|\chi-\chi_{*}|.$ $$

Case $\chi<\chi_{*}$. We have $\rho=\rho_{*}-h,\chi=\Xi(\rho_{*}-h)=\Xi(\rho)$ and $\Xi(\rho_{*}-\gamma h)\leq\chi_{new}$. There are three subcases for $\chi,\chi_{new}$ positions around 0 0:

- 1.
For positions $\chi_{*}>\chi_{new}>\chi\geq 0$, we use convexity of $\Xi$ on $[0,1)$ and obtain
$\displaystyle\chi_{new}-\chi$
$\displaystyle\geq$
$\displaystyle\Xi(\rho_{*}-\gamma h)-\Xi(\rho_{*}-h)\geq(1-\gamma)h\cdot\frac{d
\Xi}{d\rho}(\rho_{*}-h)$
$\displaystyle\geq$
$\displaystyle\min_{\rho^{\prime}\in[\rho_{*}-h,\rho_{*}]}\left|\frac{d\Xi}{d
\rho}(\rho^{\prime})\right|(1-\gamma)|\rho-\rho_{*}|.$
- 2.
For positions $\chi_{*}>0\geq\chi_{new}>\chi$, we use concavity of $\Xi$ on $(-1,0]$ and obtain
$\displaystyle\chi_{new}-\chi$
$\displaystyle\geq$
$\displaystyle\Xi(\rho_{*}-\gamma h)-\Xi(\rho_{*}-h)\geq(1-\gamma)h\cdot\frac{d
\Xi}{d\rho}(\rho_{*}-\gamma h)$
$\displaystyle\geq$
$\displaystyle\min_{\rho^{\prime}\in[\rho_{*}-h,\rho_{*}]}\left|\frac{d\Xi}{d
\rho}(\rho^{\prime})\right|(1-\gamma)|\rho-\rho_{*}|.$
- 3.
For positions $\chi_{*}>\chi_{new}>0>\chi$, we use concavity of $\Xi$ on $(-1,0]$ and convexity of $\Xi$ on $[0,1)$ and obtain
$\displaystyle\chi_{new}-\chi$
$\displaystyle\geq$
$\displaystyle\Xi(\rho_{*}-\gamma h)-\Xi(\rho_{*}-h)=[\Xi(\rho_{*}-\gamma h)-
\Xi(0)]+[\Xi(0)-\Xi(\rho_{*}-h)]$
$\displaystyle\geq$
$\displaystyle(\rho_{*}-\gamma h)\cdot\frac{d\Xi}{d\rho}(0)+(h-\rho_{*})\cdot
\frac{d\Xi}{d\rho}(0)=(1-\gamma)h\cdot\frac{d\Xi}{d\rho}(0)$
$\displaystyle\geq$
$\displaystyle\min_{\rho^{\prime}\in[\rho_{*}-h,\rho_{*}]}\left|\frac{d\Xi}{d
\rho}(\rho^{\prime})\right|(1-\gamma)|\rho-\rho_{*}|.$

Overall, we make the bound

$$ $\displaystyle\chi_{new}-\chi$ $\displaystyle\geq$ $\displaystyle\min_{\rho^{\prime}\in[\rho_{*}-h,\rho_{*}]}\left|\frac{d\Xi}{d \rho}(\rho^{\prime})\right|(1-\gamma)|\rho-\rho_{*}|$ $\displaystyle\geq$ $\displaystyle\frac{\min\limits_{\rho^{\prime}\in[\rho,\rho_{*}]}|\frac{d\Xi}{d \rho}(\rho^{\prime})|}{\max\limits_{\rho^{\prime}\in[\rho,\rho_{*}]}|\frac{d \Xi}{d\rho}(\rho^{\prime})|}(1-\gamma)|\chi-\chi_{*}|.$ $$

For the function $\Xi(\rho)=\frac{\rho}{\sigma_{0}\sigma_{1}(1-\rho^{2})}$, the centrally symmetrical derivative is

$$ $\displaystyle\frac{d\Xi}{d\rho}(\rho)$ $\displaystyle=$ $\displaystyle\frac{1+\rho^{2}}{\sigma_{0}\sigma_{1}(1-\rho^{2})^{2}}.$ $$

The derivative $\frac{d\Xi}{d\rho}$ has its global minimum at $\rho=0$. It grows as $\rho\to\pm 1$, hence, the maximum value is achieved at points which are farthest from 0 0:

$$ $\displaystyle\max\limits_{\rho^{\prime}\in[\rho_{*},\rho]}\left|\frac{d\Xi}{d \rho}(\rho^{\prime})\right|$ $\displaystyle\leq$ $\displaystyle\frac{d\Xi}{d\rho}(\rho),$ $\displaystyle\max\limits_{\rho^{\prime}\in[\rho,\rho_{*}]}\left|\frac{d\Xi}{d \rho}(\rho^{\prime})\right|$ $\displaystyle\leq$ $\displaystyle\max\left\{\frac{d\Xi}{d\rho}(\rho_{*}),\frac{d\Xi}{d\rho}(|\rho| )\right\},$ $\displaystyle\min\limits_{\rho^{\prime}\in[-1,+1]}\left|\frac{d\Xi}{d\rho}( \rho^{\prime})\right|$ $\displaystyle\geq$ $\displaystyle\frac{1}{\sigma_{0}\sigma_{1}}.$ $$

Thus, we prove the bound

$$ $|\chi-\chi_{*}|-|\chi_{new}-\chi_{*}|=|\chi_{new}-\chi|\geq\frac{(1-\max\{\rho _{*},|\rho|\}^{2})^{2}}{1+\max\{\rho_{*},|\rho|\}^{2}}(1-\gamma)|\chi-\chi_{*}|.$ $$

$$ $|\chi_{new}-\chi_{*}|\leq\left[1-(1-\gamma)\frac{(1-\max\{\rho_{*},|\rho|\}^{2 })^{2}}{1+\max\{\rho_{*},|\rho|\}^{2}}\right]|\chi-\chi_{*}|.$ $$

∎

###### Lemma B.1 (Correlation improvement after (D)IMF step) .

###### Proof.

###### Lemma B.2 ( χ 𝜒 \chi italic_χ improvement after (D)IMF step) .

###### Proof.

### B.3 IPF step analysis

We use IPMF with parameter $\varepsilon_{*}$ between distributions $\mathcal{N}(\mu_{0},\sigma_{0}^{2})$ and $\mathcal{N}(\mu_{1},\sigma_{1}^{2})$. We start with the process $\mathcal{N}\left(\begin{pmatrix}\mu_{0}\\
\nu\end{pmatrix},\begin{pmatrix}\sigma_{0}^{2}&\rho\sigma_{0}s\\
\rho\sigma_{0}s&s^{2}\end{pmatrix}\right)$ with correlation $\rho\in(-1,1)$ .

One IPMF step consists of $4$ consecutive steps:

- 1.
IMF step refining current optimality value,
- 2.
IPF step changing final prior to $\mathcal{N}(\mu_{1},\sigma_{1}^{2})$,
- 3.
IMF step refining current optimality value,
- 4.
IPF step changing starting prior to $\mathcal{N}(\mu_{0},\sigma_{0}^{2})$.

We use the following notations for the covariance matrices changes:

$$ $\displaystyle\begin{pmatrix}\sigma_{0}^{2}&\rho\sigma_{0}s\\ \rho\sigma_{0}s&s^{2}\end{pmatrix}$ $\displaystyle\overset{IMF}{\Longrightarrow}$ $\displaystyle\begin{pmatrix}\sigma_{0}^{2}&\tilde{\rho}\sigma_{0}s\\ \tilde{\rho}\sigma_{0}s&s^{2}\end{pmatrix}\overset{IPF}{\Longrightarrow}\begin {pmatrix}(s^{\prime})^{2}&\rho^{\prime}s^{\prime}\sigma_{1}\\ \rho^{\prime}s^{\prime}\sigma_{1}&\sigma_{1}^{2}\end{pmatrix}$ $\displaystyle\overset{IMF}{\Longrightarrow}$ $\displaystyle\begin{pmatrix}(s^{\prime})^{2}&\hat{\rho}s^{\prime}\sigma_{1}\\ \hat{\rho}s^{\prime}\sigma_{1}&\sigma_{1}^{2}\end{pmatrix}\overset{IPF}{ \Longrightarrow}\begin{pmatrix}\sigma_{0}^{2}&\rho^{\prime\prime}\sigma_{0}s^{ \prime\prime}\\ \rho^{\prime\prime}\sigma_{0}s^{\prime\prime}&(s^{\prime\prime})^{2}\end{ pmatrix},$ $$

and for the means changes:

$$ $\begin{pmatrix}\mu_{0}\\ \nu\end{pmatrix}\overset{IMF}{\Longrightarrow}\begin{pmatrix}\mu_{0}\\ \nu\end{pmatrix}\overset{IPF}{\Longrightarrow}\begin{pmatrix}\nu^{\prime}\\ \mu_{1}\end{pmatrix}\overset{IMF}{\Longrightarrow}\begin{pmatrix}\nu^{\prime} \\ \mu_{1}\end{pmatrix}\overset{IPF}{\Longrightarrow}\begin{pmatrix}\mu_{0}\\ \nu^{\prime\prime}\end{pmatrix}.$ $$

###### Lemma B.3 (Improvement after IPF steps) .

Consider an initial $2$-dimensional Gaussian joint distribution
$\mathcal{N}\left(\begin{pmatrix}\mu_{0}\\
\nu\end{pmatrix},\begin{pmatrix}\sigma_{0}^{2}&\rho\sigma_{0}s\\
\rho\sigma_{0}s&s^{2}\end{pmatrix}\right),\rho\in(-1,1)$. We run IPMF step between distributions $\mathcal{N}(\mu_{0},\sigma_{0}^{2})$ and $\mathcal{N}(\mu_{1},\sigma_{1}^{2})$ and obtain new joint distribution $\mathcal{N}\left(\begin{pmatrix}\mu_{0}\\
\mu^{\prime\prime}\end{pmatrix},\begin{pmatrix}\sigma_{0}^{2}&\hat{\rho}\sigma
_{0}s^{\prime\prime}\\
\hat{\rho}\sigma_{0}s^{\prime\prime}&(s^{\prime\prime})^{2}\end{pmatrix}\right)$. Then, the distance between ground truth $\mu_{1},\sigma_{1}$ and new joint distribution parameters decreases as:

$$ $\displaystyle|(s^{\prime\prime})^{2}-\sigma_{1}^{2}|$ $\displaystyle\leq$ $\displaystyle\hat{\rho}^{2}|s^{2}-\sigma_{1}^{2}|,$ (37) $\displaystyle|\nu^{\prime\prime}-\mu_{1}|$ $\displaystyle\leq$ $\displaystyle|\hat{\rho}||\nu-\mu_{1}|,$ (38) $$

where $\hat{\rho}$ is the correlation after the second IMF step.

###### Proof.

During IPF step, we keep conditional distribution and change prior. For the first IPF, we keep the inner part $x_{0}|x_{1}=\mathcal{N}\left(\mu_{0}+\frac{\sigma_{0}}{\sigma}\rho(x_{1}-\mu),
(1-\rho^{2})\sigma_{0}^{2}\right)$. The same part should remain after IPF:

$$ $\mathcal{N}\left(\mu_{0}+\frac{\sigma_{0}}{s}\tilde{\rho}(x_{1}-\nu),(1-\tilde {\rho}^{2})\sigma_{0}^{2}\right)=\mathcal{N}\left(\nu^{\prime}+\frac{s^{\prime }}{\sigma_{1}}\rho^{\prime}(x_{1}-\mu_{1}),(1-(\rho^{\prime})^{2})(s^{\prime}) ^{2}\right),\quad\forall x_{1}\in\mathbb{R}.$ $$

This is equivalent to the system of equations:

$$ $\displaystyle\frac{\sigma_{0}}{s}\tilde{\rho}$ $\displaystyle=$ $\displaystyle\frac{s^{\prime}}{\sigma_{1}}\rho^{\prime},$ (39) $\displaystyle(1-\tilde{\rho}^{2})\sigma_{0}^{2}$ $\displaystyle=$ $\displaystyle(1-(\rho^{\prime})^{2})(s)^{\prime 2},$ (40) $\displaystyle\mu_{0}-\frac{\sigma_{0}}{s}\tilde{\rho}\nu$ $\displaystyle=$ $\displaystyle\nu^{\prime}-\frac{s^{\prime}}{\sigma_{1}}\rho^{\prime}\mu_{1}.$ $$

We simplify the system to

$$ $\displaystyle\rho^{\prime}$ $\displaystyle=$ $\displaystyle\frac{\sigma_{0}\sigma_{1}}{ss^{\prime}}\tilde{\rho},$ $\displaystyle(1-(\tilde{\rho})^{2})\frac{\sigma_{0}^{2}}{(s^{\prime})^{2}}$ $\displaystyle=$ $\displaystyle\left(1-(\tilde{\rho})^{2}\frac{\sigma_{1}^{2}}{s^{2}}\frac{ \sigma^{2}_{0}}{(s^{\prime})^{2}}\right),$ $\displaystyle 1-\frac{(s^{\prime})^{2}}{\sigma_{0}^{2}}$ $\displaystyle=$ $\displaystyle\tilde{\rho}^{2}\left(1-\frac{\sigma_{1}^{2}}{s^{2}}\right),$ (41) $\displaystyle\mu_{0}-\nu^{\prime}$ $\displaystyle=$ $\displaystyle\frac{\sigma_{0}}{s}\tilde{\rho}(\nu-\mu_{1}).$ $$

For the second IPF step, we similarly obtain the following transaction equations:

$$ $\displaystyle 1-\frac{(s^{\prime\prime})^{2}}{\sigma_{1}^{2}}$ $\displaystyle=$ $\displaystyle\hat{\rho}^{2}\left(1-\frac{\sigma_{0}^{2}}{(s^{\prime})^{2}} \right),$ $\displaystyle\mu_{1}-\nu^{\prime\prime}$ $\displaystyle=$ $\displaystyle\frac{\sigma_{1}}{s^{\prime}}\hat{\rho}(\nu^{\prime}-\mu_{0}).$ $$

Let us define $\kappa\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{s^{2}}{\sigma_{1}^{2}}$ and $\kappa^{\prime\prime}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{(s^{\prime
\prime})^{2}}{\sigma_{1}^{2}}$ and find the connection between them, namely,

$$ $\displaystyle\frac{(s^{\prime})^{2}}{\sigma_{0}^{2}}$ $\displaystyle=$ $\displaystyle 1-\tilde{\rho}^{2}\left(1-\frac{1}{\kappa}\right),$ $\displaystyle 1-s^{\prime\prime}$ $\displaystyle=$ $\displaystyle\hat{\rho}^{2}\left(1-\frac{\sigma_{0}^{2}}{(s^{\prime})^{2}} \right),$ $\displaystyle 1-\kappa^{\prime\prime}$ $\displaystyle=$ $\displaystyle\hat{\rho}^{2}\left(1-\frac{1}{1-\tilde{\rho}^{2}(1-\frac{1}{ \kappa})}\right)=\hat{\rho}^{2}\left(\kappa\frac{1-\tilde{\rho}^{2}(1-\frac{1} {\kappa})-1}{\kappa-\tilde{\rho}^{2}(\kappa-1)}\right)$ $\displaystyle=$ $\displaystyle-\hat{\rho}^{2}\left(\frac{\tilde{\rho}^{2}(\kappa-1)}{1+(1- \tilde{\rho}^{2})(\kappa-1)}\right),$ $\displaystyle\kappa^{\prime\prime}-1$ $\displaystyle=$ $\displaystyle\frac{\hat{\rho}^{2}\tilde{\rho}^{2}(\kappa-1)}{1+(1-\tilde{\rho} ^{2})(\kappa-1)}.$ $$

We define $r\stackrel{{\scriptstyle\text{def}}}{{=}}\kappa-1$ and $r^{\prime\prime}\stackrel{{\scriptstyle\text{def}}}{{=}}\kappa^{\prime\prime}-1$ and use them as the distance between current step parameters and true solution. Note, that values $r,r^{\prime\prime}$ lie within $(-1,+\infty)$, and we have transaction

$$ $r^{\prime\prime}=\frac{\hat{\rho}^{2}\tilde{\rho}^{2}r}{1+(1-\tilde{\rho}^{2}) r}.$ (42) $$

We also can estimate norm of the next $r^{\prime\prime}$ depending on previous one. We bound $0\leq\frac{\tilde{\rho}^{2}}{1+(1-\tilde{\rho}^{2})r}\leq\frac{\tilde{\rho}^{2
}}{1-(1-\tilde{\rho}^{2})}\leq 1$, since there are bounds $r>-1,1\geq\tilde{\rho}^{2}\geq 0$. Hence, we obtain

$$ $\displaystyle|r^{\prime\prime}|$ $\displaystyle\leq$ $\displaystyle\hat{\rho}^{2}|r|,$ $\displaystyle|(s^{\prime\prime})^{2}-\sigma_{1}^{2}|$ $\displaystyle\leq$ $\displaystyle\hat{\rho}^{2}|\sigma^{2}-\sigma_{1}^{2}|.$ $$

For the means, we have transaction

$$ $\displaystyle\mu_{1}-\nu^{\prime\prime}$ $\displaystyle=$ $\displaystyle\frac{\sigma_{0}}{s^{\prime}}\frac{\sigma_{1}}{s}\tilde{\rho}\hat {\rho}(\mu_{1}-\nu)=\frac{\tilde{\rho}\hat{\rho}(\mu_{1}-\nu)}{\sqrt{\kappa} \sqrt{1-\tilde{\rho}^{2}\left(1-\frac{1}{\kappa}\right)}}$ (43) $\displaystyle=$ $\displaystyle\hat{\rho}\operatorname{sign}(\tilde{\rho})(\mu_{1}-\nu)\sqrt{ \frac{\tilde{\rho}^{2}}{1+(1-\tilde{\rho}^{2})(\kappa-1)}}.$ $$

Similarly, we have a linear bound for the next norm difference:

$$ $|\nu^{\prime\prime}-\mu_{1}|\leq|\hat{\rho}||\nu-\mu_{1}|.$ $$

∎

###### Lemma B.4 (IPF step does not change optimality coefficient $\chi$ ) .

Consider an initial $2$-dimensional Gaussian joint distribution
$\mathcal{N}\left(\begin{pmatrix}\mu_{0}\\
\nu\end{pmatrix},\begin{pmatrix}\sigma_{0}^{2}&\rho\sigma_{0}s\\
\rho\sigma_{0}s&s^{2}\end{pmatrix}\right)$ with correlation $\rho\in(-1,1)$. We run IPF step between distributions $\mathcal{N}(\mu_{0},\sigma_{0}^{2})$ and $\mathcal{N}(\mu_{1},\sigma_{1}^{2})$ and obtain new joint distribution $\mathcal{N}\left(\begin{pmatrix}\nu^{\prime}\\
\mu_{1}\end{pmatrix},\begin{pmatrix}(s^{\prime})^{2}&\rho^{\prime}s^{\prime}
\sigma_{1}\\
\rho^{\prime}s^{\prime}\sigma_{1}&\sigma_{1}^{2}\end{pmatrix}\right)$. Then, IPF step does not change $\chi=\Xi(\rho,\sigma_{0},s)$, i.e.,

$$ $\Xi(\rho,\sigma_{0},s)=\Xi(\rho^{\prime},s^{\prime},\sigma_{1}).$ $$

###### Proof.

We use inversion $\rho=P(\chi,\sigma_{0},s)$. During IPF step, equations (39) and (40) hold true, i.e.,

$$ $\displaystyle(1-\rho^{2})\sigma_{0}^{2}$ $\displaystyle=$ $\displaystyle\left(1-\left(\frac{\sqrt{\chi^{2}\sigma_{0}^{2}s^{2}+1/4}-1/2}{ \chi\sigma_{0}s}\right)^{2}\right)\sigma_{0}^{2}$ $\displaystyle=$ $\displaystyle\left(\frac{\chi^{2}\sigma_{0}^{2}s^{2}-(\chi^{2}\sigma_{0}^{2}s^ {2}+1/4)-1/4+\sqrt{\chi^{2}\sigma_{0}^{2}s^{2}+1/4}}{\chi^{2}\sigma_{0}^{2}s^{ 2}}\right)\sigma_{0}^{2}$ $\displaystyle=$ $\displaystyle\frac{\left(\frac{\sqrt{\chi^{2}\sigma_{0}^{2}\sigma_{1}^{2}+1/4} -1/2}{\chi\sigma_{0}s}\right)\sigma_{0}}{s\chi}=\frac{\rho\sigma_{0}}{s\chi}.$ $$

Consequently, we have equality $\chi=\frac{\rho\sigma_{0}/s}{(1-\rho^{2})\sigma_{0}^{2}}$, and, since numerator and denominator remains the same due to (39) and (40), $\chi$ does not change.
∎

###### Lemma B.3 (Improvement after IPF steps) .

###### Proof.

###### Lemma B.4 (IPF step does not change optimality coefficient χ 𝜒 \chi italic_χ ) .

###### Proof.

### B.4 Proof of IPMF Convergence Theorem 3.2

###### Proof.

Notations. We denote the variance of the 0 0-th marginal after the $k$-th IPMF step as $s^{\prime}_{k}$. For the first one, we have formula $s^{\prime}_{0}=\sqrt{\sigma_{0}^{2}-\sigma_{0}^{2}\tilde{\rho}_{0}^{2}\left(1-
\frac{\sigma_{1}^{2}}{s^{2}}\right)}$, where $\tilde{\rho}_{0}$ is the correlation after the first IMF step. More explicitly, $\tilde{\rho}_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}\rho_{new}(\rho_{0})$, where $\rho_{new}$ is taken from (29) for continuous IMF and from (31) for DIMF. We denote optimality coefficients $\chi_{k}\stackrel{{\scriptstyle\text{def}}}{{=}}\Xi(\rho_{k},\sigma_{0},s_{k})$ and $\chi_{*}\stackrel{{\scriptstyle\text{def}}}{{=}}\nicefrac{{1}}{{\varepsilon_{*
}}}$.

Ranges. Due to update equations for $\chi$ (36) and for $s,\nu$ (42), the fact that IMF step keeps $s,\nu$, while IPF keeps $\chi$, the parameters $s_{k},\chi_{k}$ remain on the same side from $\sigma_{1},\frac{1}{\varepsilon_{*}}$, respectively. Namely, we have ranges for the variances $s_{k}\in[\sigma_{1}^{min},\sigma_{1}^{max}]\stackrel{{\scriptstyle\text{def}}}
{{=}}[\min\{\sigma_{1},s_{0}\},\max\{\sigma_{1},s_{0}\}],s^{\prime}_{k}\in[
\sigma_{0}^{min},\sigma_{0}^{max}]\stackrel{{\scriptstyle\text{def}}}{{=}}[
\min\{\sigma_{0},s^{\prime}_{0}\},\max\{\sigma_{0},s^{\prime}_{0}\}]$ and parameters $\chi_{k}\in[\chi^{min},\chi^{max}]\stackrel{{\scriptstyle\text{def}}}{{=}}[
\min\{\chi_{*},|\chi_{0}|\},\max\{\chi_{*},|\chi_{0}|\}].$

Update bounds. We use update bounds for $\chi$ (35) twice, for $s$ (37) and for $\nu$ (38), however, we need to limit above the coefficients $|P(\chi,\sigma,\sigma^{\prime})|$ and $l(P(\chi,\sigma,\sigma^{\prime}),P(\chi_{*},\sigma,\sigma^{\prime}),\gamma(
\sigma,\sigma^{\prime}))$ over the considered ranges of the parameters $\sigma\in[\sigma_{0}^{min},\sigma_{0}^{max}],\sigma^{\prime}\in[\sigma_{1}^{
min},\sigma_{1}^{max}]$ and $\chi\in[\chi^{min},\chi^{max}]$. The functions $P,l,\gamma$ are defined in (26), (35), (27) (or (28) with fixed $t$), respectively.

Since the function $|P(\chi,\sigma,\sigma^{\prime})|$ is increasing w.r.t. $\sigma,\sigma^{\prime}$ and $\chi$ symmetrically around 0 0, we take maximal values $\sigma_{0}^{max},\sigma_{1}^{max}$ and $\chi^{max}$. Similarly, the function $l(P(\chi,\sigma,\sigma^{\prime}),P(\chi_{*},\sigma,\sigma^{\prime}),\gamma(
\sigma,\sigma^{\prime}))$ is increasing w.r.t. all arguments symmetrically around 0 0. Hence, we maximize the function $|P|$ and the function $\gamma$, which is also increasing w.r.t. $\sigma$ and $\sigma^{\prime}$.

Final bounds. The final bound after $k$ step of IPMF are:

$$ $\displaystyle|s_{k}^{2}-\sigma_{1}^{2}|$ $\displaystyle\leq$ $\displaystyle\alpha^{2k}|s_{0}^{2}-\sigma_{1}^{2}|,$ $\displaystyle|\nu_{k}-\mu_{1}|$ $\displaystyle\leq$ $\displaystyle\alpha^{k}|\nu_{0}-\mu_{1}|,$ $\displaystyle|\chi_{k}-\nicefrac{{1}}{{\varepsilon_{*}}}|$ $\displaystyle\leq$ $\displaystyle\beta^{2k}|\chi_{0}-\nicefrac{{1}}{{\varepsilon_{*}}}|,$ $$

where $\alpha\stackrel{{\scriptstyle\text{def}}}{{=}}P\left(\chi^{max},\sigma_{0}^{
max},\sigma_{1}^{max}\right)$ and $\beta\stackrel{{\scriptstyle\text{def}}}{{=}}l(P(\chi^{max},\sigma_{0}^{max},
\sigma_{1}^{max}),P(\chi_{*},\sigma_{0}^{max},\sigma_{1}^{max}),\gamma(\sigma_
{0}^{max},\sigma_{1}^{max}))$ taking $l$ from (35), $\gamma$ from (27) for continuous IMF and from (28) with fixed $t$ for discrete IMF.∎

###### Proof.

## Appendix C Additional Experiments

### C.1 SB benchmark

We additionally study how well implementations of IPMF procedure starting from different starting processes map initial distribution $p_{0}$ into $p_{1}$ by measuring the metric $\text{B}\mathbb{W}_{2}^{2}\text{-UVP}$ also proposed by the authors of the benchmark (Gushchin et al., ). We present the results in Table 2. One can observe that DSBM initialized from different starting processes has quite close results and so is the case for ASBM experiments with $\epsilon\in\{1,10\}$, but with $\epsilon=0.1$ one can notice that ASBM starting from IPF and Independent $p_{0}\rightarrow p_{0}$ experience a decline in $\text{B}\mathbb{W}_{2}^{2}\text{-UVP}$ metric.

**Table 2: Comparisons of $\text{B}\mathbb{W}_{2}^{2}\text{-UVP}\downarrow$ (%) between the ground truth static SB solution $p^{T}(x_{0},x_{1})$ and the learned solution on the SB benchmark. The best metric over is bolded. Results marked with $\dagger$ are taken from (Gushchin et al., ) or (Gushchin et al., ).**
|  |  | $\epsilon=0.1$ | $\epsilon=1$ | $\epsilon=10$ |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | Algorithm Type | $D\!=\!2$ | $D\!=\!16$ | $D\!=\!64$ | $D\!=\!128$ | $D\!=\!2$ | $D\!=\!16$ | $D\!=\!64$ | $D\!=\!128$ | $D\!=\!2$ | $D\!=\!16$ | $D\!=\!64$ | $D\!=\!128$ |
| Best algorithm on benchmark^† | Varies | $0.016$ | $0.05$ | $0.25$ | $0.22$ | $0.005$ | $0.09$ | $0.56$ | $0.12$ | $0.01$ | $0.02$ | $0.15$ | $0.23$ |
| DSBM-IMF | IPMF | $0.1$ | $0.14$ | $0.44$ | $3.2$ | $0.13$ | $0.1$ | $0.91$ | $6.67$ | $0.1$ | $5.17$ | $66.7$ | $356$ |
| DSBM-IPF | $0.35$ | $0.6$ | $0.6$ | $1.62$ | $0.01$ | $0.18$ | $0.91$ | $6.64$ | $0.2$ | $3.78$ | $81$ | $206$ |  |
| DSBM-Ind $p_{0}\rightarrow p_{0}$ | $0.08$ | $0.38$ | $0.62$ | $1.72$ | $0.13$ | $0.18$ | $0.84$ | $7.45$ | $0.04$ | $3.72$ | $99.3$ | $251$ |  |
| $\text{SF}^{2}$M-Sink^† | $0.04$ | $0.18$ | $\mathbf{0.39}$ | $\mathbf{1.1}$ | $0.07$ | $0.3$ | $4.5$ | $17.7$ | $0.17$ | $4.7$ | $316$ | $812$ |  |
| ASBM-IMF^† | $\mathbf{0.016}$ | $\mathbf{0.1}$ | $0.85$ | $11.05$ | $\mathbf{0.02}$ | $0.34$ | $\mathbf{1.57}$ | $\mathbf{3.8}$ | $\mathbf{0.013}$ | $\mathbf{0.25}$ | $\mathbf{1.7}$ | $\mathbf{4.7}$ |  |
| ASBM-IPF | $0.05$ | $0.73$ | $32.05$ | $10.67$ | $0.02$ | $0.53$ | $4.19$ | $10.11$ | $0.002$ | $0.18$ | $2.2$ | $5.08$ |  |
| ASBM-Ind $p_{0}\rightarrow p_{0}$ | $0.36$ | $0.76$ | $16.33$ | $22.63$ | $0.07$ | $0.48$ | $1.93$ | $5.36$ | $0.04$ | $0.23$ | $1.04$ | $2.29$ |  |

### C.2 Celeba

In Figure7, we provide additional quantitative study of IPMF convergence in CMMD (Jayasumana et al., ) and Mean Squared Error between inputs and outputs of translation by plotting them as a function of IPMF iteration. Both metrics are calculated on Celeba male$\rightarrow$female (64 $\times$ 64) test set. We notice CMMD plot resembles one of FID, see Figure 5, and MSE plot is approximately the same for the DSBM and ASBM groups.

Figure: (a) CMMD
Refer to caption: extracted/5898938/pics/celeba/IPMF_nfe_100_f_eps_1_Celeba64_IMF_cmmd.png

## Appendix D Experimental Details

### D.1 General details

Authors of ASBM (Gushchin et al., ) kindly provided us the code for all the experiments. All the hyperparameters including neural networks architectures were chosen as close as possible to the ones used by the authors of ASBM in their experimental section. Particularly, as it is described in (Gushchin et al., , Appendix D), authors used DD-GAN () with Brownian Bridge posterior sampling instead of DDPM’s one and implementation from:

[https://github.com/NVlabs/denoising-diffusion-gan](https://github.com/NVlabs/denoising-diffusion-gan)

DSBM (Shi et al., ) implementation is taken from the official code repository:

[https://github.com/yuyang-shi/dsbm-pytorch](https://github.com/yuyang-shi/dsbm-pytorch)

Sampling on the inference stage is done by Euler Maryama SDE numerical solver (Kloeden, ) with indicated in Table 3 NFE.

Independent $p_{0}\rightarrow p_{0}$ starting process in all the experiments was implemented in mini batch manner, i.e., $\{x_{0,n}\}_{n=1}^{N}\sim p_{0}$ and $x_{1}$ batch $\{x_{1,n}\}_{n=1}^{N}\sim q^{0}(\cdot|\{x_{0,n}\}_{n=1}^{N})$ is generated by permutation of $\{x_{0,n}\}_{n=1}^{N}$ mini batch indices.

The Exponential Moving Average (EMA) has been used to enhance generator’s training stability of both ASBM and DSBM. The parameters of the EMA are provided in Table 3, in case the EMA decay is set to ”N/A” no averaging has been applied.

**Table 3: Hyperparameters of models from Celeba 5, SwissRoll 4.2, cMNIST 3 and SB Bench experiments 4.3. In ”Start process” column ”All” states for all the used options. ”N/A” state for either not used or not applicable corresponding option for the algorithm.**
| Model | Dataset | Start process | IPMF iters | IPMF-0 Grad Updates | IPMF-k Grad Updates |
| --- | --- | --- | --- | --- | --- |
| ASBM | Celeba | All | 20 | 200000 | 20000 |
| DSBM | Celeba | IMF | 20 | 100000 | 20000 |
| DSBM | Celeba | Ind $p_{0}\rightarrow p_{0}$ | 20 | 20000 | 20000 |
| ASBM | Swiss Roll | All | 20 | 400000 | 40000 |
| DSBM | Swiss Roll | All | 20 | 20000 | 20000 |
| ASBM | cMNIST | All | 20 | 75000 | 38000 |
| DSBM | cMNIST | All | 20 | 100000 | 20000 |
| ASBM | SB Bench | All | 20 | 133000 | 67000 |
| DSBM | SB Bench | All | 20 | 20000 | 20000 |

### D.2 Illustrative 2D examples

ASBM. For toy experiments the MLP with hidden layers $[256,256,256]$ has been chosen for both discriminator and generator. The generator takes vector of $(dim+1+2)$ length with data, latent variable and embedding (a simple lookup table torch.nn.Embedding) dimensions, respectively. The networks have torch.nn.LeakyReLU as activation layer with $0.2$ angle of negative slope. The optimization has been conducted using torch.optim.Adam with running averages coefficients $0.5$ and $0.9$. Additionally, the CosineAnnealingLR scheduler has been used only at pretraining iteration with minimal learning rate set to 1e-5 and no restarting. To stabilize GAN training R1 regularizer with coefficient $0.01$ (Mescheder et al., ) has been used.

DSBM. MLP with $[\text{dim}+12,128,128,128,128,128,\text{dim}]$ number of hidden neurons, torch.nn.SiLU activation functions, residual connections between 2nd/4th and 4th/6th layers and Sinusoidal Positional Embedding has been used.

### D.3 SB benchmark

Scrödinger Bridges/Entropic Optimal Transport Benchmark (Gushchin et al., ) and $\text{cB}\mathbb{W}_{2}^{2}\text{-UVP}$, $\text{B}\mathbb{W}_{2}^{2}\text{-UVP}$ metric implementation was taken from the official code repository:

[https://github.com/ngushchin/EntropicOTBenchmark](https://github.com/ngushchin/EntropicOTBenchmark)

Conditional plan metric $\text{cB}\mathbb{W}_{2}^{2}\text{-UVP}$ , see Table 1, was calculated over predefined test set and conditional expectation per each test set sample estimated via Monte Carlo integration with 1000 samples. Target distribution fitting metric, $\text{B}\mathbb{W}_{2}^{2}\text{-UVP}$, see Table 2, was estimated using Monte Carlo method and 10000 samples.

ASBM.
The same architecture and optimizer have been used as in toy experiments D.2, but without the scheduler.

DSBM. MLP with $[\text{dim}+12,128,128,128,128,128,\text{dim}]$ number of hidden neurons, torch.nn.SiLU activation functions, residual connections between 2nd/4th and 4th/6th layers and Sinusoidal Positional Embedding has been used.

### D.4 CMNIST

Working with MNIST dataset, we use regular train/test split with 60000 images and 10000 images correspondingly. We RGB color train and test digits of classes ”2” and ”3”. Each sample is resized to $32\times 32$ and normalized by $0.5$ mean and $0.5$ std.

ASBM.
The cMNIST setup mainly differs by the architecture used. The generator model is built upon the NCSN++ architecture (), following the approach in () and (Gushchin et al., ). We use 2 residual and attention blocks, 128 base channels, and $(1,2,2,2)$ feature multiplications per corresponding resolution level. The dimension of the latent vector has been set to 100. Following the best practices of time-dependent neural networks sinusoidal embeddings are employed to condition on the integer time steps, with a dimensionality equal to $2\times$ the number of initial channel, resulting in a 256-dimensional embedding. The discriminator adopts ResNet-like architecture with 4 resolution levels. The same optimizer with the same parameters as in toy D.2 and SB benchmark D.3 experiments have been used except ones that are presented in Table 3. No scheduler has been applied. Additionally, R1 regularization is applied to the discriminator with a coefficient of 0.02, in line with () and (Gushchin et al., ).

DSBM.
The model is based on the U-Net architecture (Ronneberger et al., ) with attention blocks, 2 residual blocks per level, 4 attention heads, 128 base channels, $(1,2,2,2)$ feature multiplications per resolution level. Training was held by Adam (Kingma & Ba, ) optimizer.

### D.5 Celeba

Test FID, see Figure 5 is calculated using [pytorch-fid package](https://github.com/mseitzer/pytorch-fid), test CMMD is calculated using [unofficial implementation](https://github.com/sayakpaul/cmmd-pytorch) in PyTorch. Working with CelabA dataset (Liu et al., ), we use all 84434 male and 118165 female samples ($90\%$ train, $10\%$ test of each class). Each sample is resized to $64\times 64$ and normalized by $0.5$ mean and $0.5$ std.

ASBM. As in cMNIST experiments D.4 the generator model is built upon the NCSN++ architecture() but with small parameter changes. The number of initial channels has been lowered to 64, but the number of resolution levels has been increased with the following changes in feature multiplication, which were set to $(1,1,2,2,4)$. The discriminator also has been upgraded by growing the number of resolution levels up to 6. No other changes were proposed.

DSBM. Following Colored MNIST translation experiment exactly the same neural network and optimizer was used.

## Appendix E Broader impact.

This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.

## Peer Reviews

confidence: 2
contribution: 2
initial_rating: 6
presentation: 3
questions: - In section 4.3, the results presented in Table 1 are difficult to interpret. What's the role of the "volatility" parameter $\epsilon$ and why it influences so much the performances of the different procedures?
soundness: 3
strengths: - The paper is very well written.
- The problem in consideration and the background are introduced very clearly, so that the idea of the paper is easy to comprehend. 
- The proof for 1D gaussian, although limited, could lead to generalisations that improve on the results that are available in the literature at the moment.
summary: In the context of Schrödinger Bridge (SB) problems, the paper analyses theoretically some heuristic algorithms typically used in practice in the literature. Specifically, the authors give a new interpretation to these heuristics, as a combination of well-known procedures used in the context of theoretical analysis of SB problems. The authors name the resulting procedure "Iterative Proportional Markov Fitting" (IPMF).

The main result of the paper is to show the converge of IPMF when both the distributions $p_0$ and $p_1$ are unidimensional gaussians.
Furthermore, the authors study heuristically the convergence of IPMF in more interesting settings, by performing numerical experiments on some typical datasets and benchmarks.
weaknesses: - Section 4, presenting the experiments, seems to be a bit disconnected from the previous sections. While at the beginning the authors say the aim of the section is to show the convergence of the different procedures in more generic settings, the focus shifts on the comparison between the performance of different procedures. Keeping the focus on the convergence of the procedures to the target could enhance the readability of the section.

[Minor] 
Typos:
- Line 172: Matkovian
- Line 374: staring
- Line 412: Additionaly
- Line 417: Scrödinger
- Line 528: Appenix

A careful reading of the paper is suggested to remove additional typos.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: please see above
soundness: 3
strengths: The paper is clear, and in the first pages it adopts a tutorial-like style. The presentation of the relevant background (IPF & IMF) is properly achieved and the proposal is delivered in context. The authors also prove a particular case of their proposal (Thm 3.2).

Experimentally, the paper includes empirical validation on synthetic (Gaussians of different dimensions) and image datasets (Celeba, MNIST). The experiments are complemented by both quantitative and qualitative performance evaluations.
summary: The paper targets the solution of the Schrodinger Bridge problem. It first revisits two standard methods to solve the problem: Iterative Markovian Fitting (IMF) and Iterative Proportional Fitting (IPF). Then, the authors identify a connection between both methods and the associated heuristics of their implementation. This leads them to propose the Iterative Proportional Markovian Fitting (IPMF). The authors present a theoretical formulation and analysis of IPMF, to then provide experimental validation over synthetic and real-world image datasets.
weaknesses: There are some inaccuracies in the text (Appenix, markovian, Matkovian, etc). 

Perhaps due to the lack of familiarity of this reviewer with the SOTA on SB: the main contribution of the paper is the identification that IPF & IMF can be understood as the same procedure. However, what are the practical implications of this? Is this a a purely conceptual statement, or does this fact help practical applications of SB? Perhaps the authors can be more specific about the practical advantages of their proposal

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. What's the technical difficulties for generalizing the proof to multi-variate Gaussians or other situations?
2. Can you provide more intuition on why IMPF can usually work for any specific starting process?
soundness: 3
strengths: This paper is well written and shows a novel interesting results on connecting the IMF and IPF with theoretical results for 1d gaussian cases and empirical results for higher dimensional datasets.
summary: This paper shows that the popular bidirectional IMF algorithm uses IPF implicitly and they provided a theoretical result proving that for 1d Gaussians, this algorithm IMPF can converges.
weaknesses: Although this paper shows an interesting theoretical results, it's only for 1d Gaussian cases. I feel like the contribution is relatively weak. If authors can say something more, i.e, for multi-variate Gaussian, this paper will become stronger.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: * Typo in line 59 “Mathcing” should be “Matching”

* Typo in line 374 – “staring” should be “starting”

* In the introduction the authors mention that the reverse Kullback-Leibler between the iterates of the IMF and the target EOT coupling decreases at each step of the IMF. However, I could not find such a result in [1, 2]. To the best of my knowledge, the theoretical study of the IMF is quite limited. 

* In Equation (9) and Equation (10) the forward and backward are done successively. Do you think it is possible to learn the backward and forward processes jointly? Is there any work investigating that direction? In that case, it would be interesting to understand the IPMF in that setting.

* It would be good to have some kind of diagram to illustrate the differences between IMF, IPF and IPMF as well as some explicit descriptions on when they are equal. Basically, the choice of the initial coupling matters a lot. If the initial coupling in the IPMF has the right marginals (and is in the reciprocal class) then we recover the IMF. If the initial coupling in the IPMF is Markov in the reciprocal class and has the right initial marginal then we recover the IPF. 

* What other starting coupling could be explored? This is related to my comment regarding the lack of applications in the “Weaknesses” section. 

* I want to point out that in Section 3.1. The rewriting of the IPMF is only meaningful if we start from a coupling that does not preserve the marginals, otherwise since the Markovian projection preserves the marginal there is no point applying the projection on the last marginal (or the first marginal). I think it is necessary to specify this. 

* It would be interesting to investigate what is the effect of approximating the projection in the convergence of the IMF, the IPF and the IPMF, similarly as in [3].

[1] Shi et al. – Diffusion Schrödinger Bridge Matching

[2] Peluchetti – Diffusion Bridge Mixture Transports, Schrödinger Bridge Problems and Generative Modeling

[3] Chen et al. – Provably Convergent Schr ̈odinger Bridge with Applications to Probabilistic
Time Series Imputation
soundness: 2
strengths: * I think the paper is nicely written. There is proper acknowledgement of prior work and the authors have carefully introduced complex notions such as Iterative Proportional Fitting, Iterative Markovian Fitting, Schrodinger Bridges and Entropic Optimal Transport. 

* The authors are also discussing an interesting (and often overlooked) point in the implementation of EOT based samplers which is the initial coupling. While I think the current paper still suffers from major issues, there is value in this study as it puts the emphasis on 1) how stable the IMF is with respect to the original coupling 2) the interaction between IMF and IPF
summary: In this paper, the authors propose a reinterpretation of the Iterative Markovian Fitting (IMF). Iterative Markovian Fitting was proposed in [1, 2]. This scheme performs Markovian and reciprocal projections in order to find the Entropic Optimal Transport (EOT). A much older scheme to find EOT is Iterative Proportional Fitting (IPF) [3, 4]. The authors propose to reinterpret a popular implementation of IMF with forward/backward projections to link this specific IMF implementation with the IPF algorithm. The procedure is called IPMF (Iterative Proportional Markovian Fitting). They show that the proposed procedure converges for one-dimensional Gaussian targets. They then illustrate experimentally the stability of the forward/backward implementation of IMF, i.e. IPMF, for several starting couplings including some that do not satisfy the requirements of IMF and IPF. The experiments are conducted on high dimensional Gaussian, on Colored MNIST and on CelebA. 

[1] Shi et al. – Diffusion Schrödinger Bridge Matching

[2] Peluchetti – Diffusion Bridge Mixture Transports, Schrödinger Bridge Problems and Generative Modeling

[3] Kullback – On information and sufficiency, Annals of Mathematics and Statistics

[4] Sinkhorn – A Relationship Between Arbitrary Positive Matrices and Doubly Stochastic Matrices
weaknesses: * I think the paper has limited novelty. Indeed, the link between the IMF procedure and the IPF was already identified in  [1]. In particular, in [1] the following coupling “We now relate Algorithm 1 to the classical IPF and practical algorithms such as DSB (De Bortoli et al., 2021). Instead of initializing DSBM with $\Pi_{0,T}^0$ given by a coupling between $\pi_0, \pi_1$ , if we initialize it by $\Pi_{0,T}^0 = \mathbb{Q}_{0,T}$ [...] then DSBM also recovers the IPF iterates used in DSB.” Hence the connection between IMF and IPF through the initial coupling was already acknowledged in [1], see Proposition 10 for an explicit statement. In that regard, the IPMF rewriting is not needed to show that the IMF procedure reverts to the IPF one when the coupling is initialized as the IPF one which is one of the main points of the paper. 

* Overall I think the paper has limited impact. If it should be read as a theory paper then the result is quite weak (one dimensional Gaussian random variables) and there is not so much insight that is gained from the proof of the theorem, i.e. it seems extremely unclear that the proof generalizes to other settings and what would be a path to reach such general convergence results. If the paper is to be seen as a methodology paper then I think that there should be a clear motivation to choose couplings that are different from the ones of the IMF (i.e. couplings with correct marginals) or the ones from the IPF (couplings with first right marginal and output of the Brownian motion or related stochastic process). So far, I do not see in the paper a clear use case of such couplings. Hence I find the paper poorly motivated.

* Finally, I find the experiments to be quite poor. The only settings that are run are high dimensional Gaussians, 2d examples and unpaired image translation with colored MNIST and Celeba. I found the quantitative study of the last two experiments to be quite poor. For instance, FID scores are presented in Figure 5 and MSE are presented in Appendix C.2 and in that case it seems that changing the coupling yields worse results regarding the quality of the obtained coupling. I think, the authors should investigate other couplings (I know that the minibatch OT one was investigated in the case of Celeba) in order to show that there is something to be gained from the IPMF perspective.

Meta-review: (a) Summary of the paper's claims and findings:

The paper introduces the Iterative Proportional Markovian Fitting (IPMF) procedure, combining aspects of the Iterative Markovian Fitting (IMF) and Iterative Proportional Fitting (IPF) algorithms. The authors theoretically demonstrate convergence of the IPMF procedure in the limited case of one-dimensional Gaussian distributions and present empirical results for high-dimensional Gaussian distributions, MNIST, and CelebA datasets. They argue that IPMF generalizes prior methods and offers practical flexibility for starting couplings, showing stability in broader settings.

(b) Strengths:

The paper is well-written, providing clear explanations of foundational methods (IMF, IPF) and their connections.
It addresses an important aspect of stability in Schrödinger Bridge implementations.
Experimental results are provided to support practical utility, including a new coupling derived using SDEdit.
The theoretical insight into the convergence for 1D Gaussian cases is rigorous and includes experimental support for broader scenarios.
(c) Weaknesses:

The theoretical results are limited to 1D Gaussian cases, offering little insight into generalization to higher dimensions or other distributions.
The novelty is weak, as prior work has already connected IMF and IPF. The proposed IPMF framework does not significantly extend the theoretical or practical frontier.
Experimental evaluations are narrow and lack comprehensive comparisons with state-of-the-art methods. While some new couplings are explored, the practical gains are marginal or unclear.
The methodology lacks a compelling motivation or application scenario where IPMF is preferable over existing techniques.

(d) Reasons for rejection:

The paper’s contribution is found to be incremental, with limited theoretical novelty and insufficient practical impact. The proposed framework reiterates known connections and does not provide a strong rationale for adopting IPMF in place of established methods. Despite the authors' responses and updates, the scope of the results and their applicability remain constrained.

Decision: Reject