Learning non-equilibrium diffusions with Schrödinger bridges: from exactly solvable to simulation-free
Authors: Stephen Y. Zhang
Michael Stumpf
Venue: NeurIPS 2025 poster ()
Keywords: Schrodinger bridge
nonequilibrium dynamics
flow matching
Primary area: machine_learning_for_sciences
TL;DR: We study the Schrodinger bridge problem with a general Ornstein-Uhlenbeck reference dynamics that allows for modelling of non-equilibrium dynamics.

Abstract: We consider the Schrödinger bridge problem which, given ensemble measurements of the initial and final configurations of a stochastic dynamical system and some prior knowledge on the dynamics, aims to reconstruct the "most likely" evolution of the system compatible with the data. Most existing literature assume Brownian reference dynamics, and are implicitly limited to modelling systems driven by the gradient of a potential energy. We depart from this regime and consider reference processes described by a multivariate Ornstein-Uhlenbeck process with generic drift matrix $\mathbf{A} \in \mathbb{R}^{d \times d}$. When $\mathbf{A}$ is asymmetric, this corresponds to a non-equilibrium system in which non-gradient forces are at play: this is important for applications to biological systems, which naturally exist out-of-equilibrium. In the case of Gaussian marginals, we derive explicit expressions that characterise exactly the solution of both the static and dynamic Schrödinger bridge. For general marginals, we propose mvOU-OTFM, a simulation-free algorithm based on flow and score matching for learning an approximation to the Schrödinger bridge. In application to a range of problems based on synthetic and real single cell data, we demonstrate that mvOU-OTFM achieves higher accuracy compared to competing methods, whilst being significantly faster to train.

Learning non-equilibrium diffusions with Schrödinger bridges: from exactly solvable to simulation-free
Stephen Zhang" University of Melbourne
Michael P. H. Stumpf University of Melbourne
Abstract
We consider the Schrödinger bridge problem which, given ensemble measurements of the initial and final configurations of a stochastic dynamical system and some prior knowledge on the dynamics, aims to reconstruct the "most likely" evolution of the system compatible with the data. Most existing literature assume Brownian reference dynamics, and are implicitly limited to modelling systems driven by the gradient of potential energy. We depart from this regime and consider reference processes described by multivariate matrix E R When is asymmetrie, this corresponds to a non-equilibrium play: this is important for applications
process with generic drift
system in which non gradient forces to biological systems, which naturally exist out-of-cquillbrium. In the case derive explicit expressions that characterise exactly the
non-grodien
Gaussian marginals, solution of both the static and dynamic Schrödinger bridge. For general marginals, we propose MVOU-OT -MFMM, simulation-fres algorithm based flow and score matching for learning an approximation to the Schrödinger bridge. In application a range of problems based on synthetic and real single cell data, we demonstrate that MVOU-OTFI achieves higher accuracy compared to competing methods, whilst being significantly faster train.
1
Introduction
We are interested in reconstruction stochastic dynamics of individuals from static population snapshots. This is a central problem with applications arising across the natural and social sciences. whenever longitudinal tracking of individuals over time is either impossible or impractical [56] 35 48
43 34 M.. In simple terms, consider system of indistinguishable particles a;, RRd undergoing some unobserved temporal dynamics. The practitioner observes the system to have distribution 2g ~ Po initial time 0 and later be x P1 final time t 1. The question is then: can we, under some suitable assumptions, reconstruct the continuous-tin behaviour of the system for the unobserved time interval 0 < 1?
The Schrödinger bridge problem (SBP), by now a centrepiece of the theoretical literature on this topic, places this task theoretical footing in terms of mathematical formulation of this problem large deviations principle on the path space (28) and a (stochastic) least action principle intimately linked to optimal transportation theory [2]. Given reference process that encodes prior knowledge on the dynamics, the SBP can be understood as identifying the "most likely" evolution of the system that is compatible with the snapshot observations.
The SBP and related topics have enjoyed a great deal of recent interest from both applied and theoretical perspectives. Many applications arise from biological modelling of cell dynamics 42]. The majority existing work assumes, either explicitly or implicitly, that the system of interest is
potential driven, that is, driven by the gradient of potential energy [26 13 52
5 51 57 6]. In fact,
*Corespondence syzesyz. id.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Schrödinger Bridge Problem
Q-bridges via control farmalism
01
min
(IP
7: Po. P.1-PO.1 Multivanate Orrstein-Uhierteeee Reference
P
dPo (z, Q
DO
dX,- = AX, dt dB:
Disintegration identity
X Lu X1 x dx. (AX,
dt or dB.
Frow network
matching 8,)(()) V Ip. (2) (2)))) L a
PI
Score network
vr(x))
(.))))
2,(2) ))) RS
Po
The Schrödinger bridge problem with multivariate Omstein-UUlenbbeck reference process (T) can be solved via generalised entropic transport problem and characterisation of the Q-bridges. For non -ou-Gaasiiin endpoints, score and flow matching provide route to building neural approxima- tions without simulation.
Figure I:
most studies do not consider any prior drift, corresponding to the setting where the reference process purely an isotropic Brownian motion. few studies S 45 15] allow for scalar Brmstein-Uhlenber (OU) processes, motivated by applications in the generative modelling domain. Since scalar OU dynamics have unidimensional drift which can always be written as the gradient of a scalar function, all of these are also potential driven systems.
Changing the choice for the reference process provides route dealing with systems which potential driven, i.e. driven by a drift that is con-conservvttve vector field. Abundant motivations for modelling these dynamics arises from biological systems and other kinds of active matter; these exist 19 17 exhibiting inreversible dynamics at -qu-equiii
naturally far from equilibrium [60 37
25 allow for reference processes described by SDE with a generic drift term, they rely computationally expensive simulation procedures such as numerical integration and suffer from accuracy issues in high dimensions. Here we explore complementary approach and concentrate our attention on a family of linear reference dynamics as a middle ground between physical relevance and analytical tractability. Specifically, we consider reference processes arising in a class of linear drift-diffusion dynamics described by the SDE
steady state [27]. While few studies [53
dX; A(X; m)) dt + dB..
(1)
These processes are also known as multivariate Ornstein- Ornstein-Uhlenbbeck Uhlenbeck (mvOU) processes, often used as model of non-equilibrium systems the physics literature (27 20 [8]. Indeed, for asymmetric drift matrix A the system (I) has a drift that is no longer the gradient of a potential. As an example in the case of isotropic diffusion part, a non-equilbbrum steady state exists where the dynamics are irreversible and exhibit nonzero probability currents while the population density is unchanged [8 20 While allowing a broader range of reference dynamics, the statistics of (I) remain analytically tractable in particular it is Markovian Gaussian process with explicit formulae for its mean, covariance, and transition kernel. This provides avenues to efficient solution of the SBP that sidesteps the need for numerical integration, and as we will see, even analytical expressions in the Gaussian case. particular, the family of SBP problems that we model is strictly larger than the standard (Brownian) SBP, which we recover as special case when A, m =
oI, when A is asymmetric and has all cigenvalues with negative real
The contributions of our paper are twofold. Leveraging results on mvOU processes, we characterise for Gaussian endpoints the Gaussian Schrödinger Bridge (GSB) with reference dynamics described by (1). To handle general marginals, we develop a simulation free training procedure based on score and flow matching [52] to solve the SBP. We conclude by demonstrating our results in a range of synthetic and real data examples. We find that our approach solves the generalised SBP faster and more accurately than comparable simulation-bart based algorithms, the cost of assuming the form [ of the reference process. Our findings highlight the tradeoff between (i) analytically tractable but more restrictive models and (ii) more expressive paramaraettric models, which are much less tractable to learn both in terms of computational cost and statistical complexity.

---

2 Background and related work
2.1 Schrödinger bridges
In what follows we provide concise summary of the SBP. For details we refer readers to 28 for in-depth discussion. We work in X' = R and denote by P(X) ) the space of probability distributions on X. Let C[ 1]. X) denote the space of continuous paths wi [0, 1] - X valued X, and informally we will refer to path measures as probability measures supported on C[[o, I). Viewing stochastic processes as random variables valued in C([0, 1], XY), path measures prescribe their law. Let Q be general path measure describing reference Markov stochastic process. For prescribed initial and final marginals po, P1 (()), the Schrödinger bridge problem can be written
min PEP(C((O.1).2)) :o Po (L(P))),
(SBP-dyn)
where the minimum is taken over all candidate processes absolutely continuous with respect to Q that are compatible with the observed data at t = 0, 1. This dynamic form of the SBP, while elegant, unwieldy for practical purposes owing the formulation on path space. well known result (2 connects the dynamical SBP with its static counterpart:
min Poi Po)eII(pc.P1)
(T. QQ1).
(SBP-static)
Furthermore, the static and dynamic SBP are interchangeable via the disintegration identity regarding the optimal law P":
x TP* (-) - ( (()))..)))))
where Qay denotes the law of the Q-bridge conditioned at (0, and (1,z1). In other words, solution of (SBP-dyn) amounts to solving (SBP-static) for Fūi followed by construction IP* as per (2) by taking mixtures of Q-bridges.
The problem (SBP-static) can be reformulated as an entropy-regularise optimal transportation problem [2] where the effective cost matrix is the log-transine -og-transttion kernel under Q This admits efficient solution via the Sinkhorn-Knopp algorithm [14] the discrete case. This provides a practical roadmap to constructing dynamical Schrödinger bridges by first solving (SBP-static), then using construction of Q-bridges with build a solution of (SBP-dyn).
2.2 Probability flows and (score, flow) -matching
Probability flows Consider a generic drift-diffusion process in d dimensions with drift (x()) and diffusivity o,, described by an Itô diffusion whose marginal densities evolve following the corresponding Fokker-Planck equation (FPE):
dx,, v. (X,) dt + o, dB,,
Ōapi(x) - (p)()))(()) + V (D)V,.()),
(3)
where D. = T is the diffusivity matrix. Equivalently, the FPE can be rewritten in the form of a continuity equation involving a probability flow field u, [30 3 2 52 31]:
(xp:(x) = (x))))
ux(x) = v. - D.V log px)))
(4)
By recognising the form of the continuity equation in it is apparent that the family of marginal dis- tributions {p: (m)) generated by the dynamics specified in (3) are also generated by the probability flow (PF)-DD:
X, =u(()) Xo Po
(5)
The central quantity that allows us convert between the SDE (3) and the PF-ODE (4) the gradient of the log- den-iiyy V_ log (x =: sz)),, also known as the score function. Under mild regularity conditions, knowledge the score s(x) allows for sampling from p(x²) (z) via Langevin dynamics: the
SDE d,X, {(x))d dW, has stationary distribution p(x²) 46). While typically one needs to
resort to approximations learn & [55] 46 47, in setting of the Brownian (and as we will see, the
mvOU) bridge, one has closed form expressions for the score and the objective.

---

Conditional flow matching Stated in its original form, let t (z(x) be a family of marginals on E [0, 1] satisfying the continuity equation (p),(z) -V- (x)) where u, is a time-dependent vector field. Suppose further that (x( admits representation as mixture
p((x)) dq( (2(=))))=(z)2))
where q(=) is some latent variable and (+)= are called the conditional probability paths. We introduce conditional flow fields u that generate the conditional probability paths, i.e. (epe)((r) = - (Ptta (2) (1)((()22)) Then, the insight presented in [30 Theorem 1J is that, in fact,
u a:) [(x))))
U = (x|22)
(<
can be easily verified by checking the continuity equation. When pr (x())) is not tractable but 7(=), (x) and (xr) are, as is the case for the Schrödinger bridge (2) conditional flow matching is useful: Lipman et al. (30 Theorem 2] prove that minimising
(+) E Ex-pis (g,,(())) (x))2
(6)
equivalent to regression on the true (marginal) vector field u4.
Simulation-free Schrödinger bridges Previous work 152 39 has exploited the connection
between the dynamic and static SBP problems to create solutions to (SBP-dyn) without the need to numerically simulate SDEs. In particular, [52] proposed to utilise score matching and flow matching simultaneously to learn the probability flow and score of the dynamical SBP with Brownian reference process. Crucially, they exploit availability of closed form expressions for the Brownian bridge and minimise the objective
LIsr]* as (0,sp)
E
(2,((z) (z) (=)( + \e1l5p.((z)
- (t(( (2)11 (7)
188 A.E P
where π P01 is the optimal coupling solving (SaP-static). In the above (x(())))1)) 8 St((xo.a1) are, respectively, the flow and score of the Brownian bridge conditioned on (aro x1 ). for which there
readily accessible closed form expressions [22 Eq. 8]: u,(x)) = 2111(] 1-21 (x xx)) + (æ)- xo) and
(x) = (o²t(1 t))
x) with o to + (1 t)x.. Motivating the use of the conditional
objective 7 in practice, the authors prove [52] Theorem 3.2] that minimising
7 is equivalent to
regressing against the SB flow and score. Since 1% can be obtained by solving an static entropic
optimal transport problem.
2.3 Related work
The generalisation of dynamical optimal transport and the Schrödinger bridge to linear reference dynamics has been studied [T 9 10 from the viewpoint of stochastic optimal control. However, as was pointed out by [5], these studies have primarily focused such as existence and uniqueness. In particular, the result for the Gaussian case [9] is in terms of a system of coupled matrix differential equations and does not lend itself to straightforward computation. More generally,
theoretical aspects of the problem,
SDEs corresponding to a continuous time iterative proportional fitting scheme [53 44 have been proposed for general reference processes.
For the Gaussian case, the availability closed form solutions is by now classical [3 49 and more recent work provides analytical characterisations of Gaussian retraiy-regulrrriedd transport 5 22 32 5 provides explicit formulae that characterise the bridge marginals as well as the force (bridge control) for Brownian reference well as class of scalar OU references. To the authors' knowledge, however, all explicit formulae for Gaussian Schrödinger bridges assume either Brownian or scalar OU dynamics. The application of flow matching [30 3] as a simulation-fis -free technique to learn Schrödinger bridges was proposed in [52). but the authors restrict consideration to a Brownian reference process. Recently [33] apply flow matching for stochastic linear control systems, however they are concerned with interpolating distributions and do not study the SBP.
Finally, concurrent work 38 proposes a related flow matching scheme for approximately solving the Q-SBP where the reference process Q is described by a general, potentially non-linear diffusion. This comes however at the cost of using learned neural approximations for bridges and log-transition densities of Q. Our work treats the complementary case of linear reference dynamics, in which we may avail ourselves of analytical formulae for these quantities.

---

3 Schrödinger bridges for non-equilibrium systems
3.1 Multivariate
bridges
As is evident from 2 the dynamical formulation of the Q-SBP amounts to solution of the static Q-SBP (SBP-static) together with characterisation of the Q-bridges, i.e. the reference process Q conditioned on initial and terminal endpoints. For + described by linear SDE of the form (.) explicit formulae for the Q-bridges are available and we state these in the form of the following theorem.
Theorem 1 (SDE characterisation of mvOU bridge, adapted from results in 8 Consider the dimensional mvOU process (T). Conditioning
20) and (T,TT) TT), the bridges of this process Y, X1I( X = Zo. XT TT) 0 & L are generated by the SDE
dY: (A(Y, m) + c(((xom))) dt d.B, Cal(mo.m))) - (Y, kt).
(8)
where A, = 50 TT- c A ds and k, = - - (m) m.
We remark that this characterisation can be found implicitly in the results of 8 however, their results were conceived for more general setting of time- varying coefficients and hence is stated in terms of the state transition matrix, which does not in general admit an explicit formula. In our case, computation of the control term A, relies upon unidimensional integral and crucially this is independent of the endpoints, meaning it incurs a one-off computational cost. Although the derivation of Theorem follows the same procedure S . provide details in the Appendix since the work of Chen and Georgiou [8) is written for a control audience. Using the SDE characterisation (8) of the Q-bridge and the fact that it remains a Gaussian process, we obtain explicit formulae for the score and flow fields as well as mean and covariance functions of the Q-bridge. One can furthermore check that the Brownian bridge formulae are recovered as special case.
key quantity
Theorem 2 (Score and flow for multivariate conditioned on (0, vo). (T, TT), denote by Pe(( the density at time O < 1 < T. Then, score function and probability flow of the bridge are respectively
bridge). For the mvOU bridge
Sil(wa.rr) (x) = Vx log Pillan, (x) = 2(((co. (Mel(20,0e)) -x)
(9)
U U1 (T) (x) A(x m) Cc((20,)eT) f00 (1((00)0))
(10)
where Cel(mng m is the bridge control from Theorem 7 and (Ho)( E ((( are the mean and covariance of Pil(so.wr7]
Σ ((),))))) = 4; D (T-) $ T (T- =: 2:
(11)
(1) = fi 2o D (T -)
T (xT PT 20
(12)
In the above, μº denotes the mean at time t of unconditioned mvOU process started from x), and the function $, Sy (t-s))Aar (ttss) ds independent of the endpoints (Io, ZT). Consequently, E(( depends only on I,
write Ω. Eq((-)) to make this explicit.
While all our theoretical results are stated in terms of a generic diffusion matrix a, we remark that in practice can always assume to be diagonal. This is formalised in the following informal lemma, for which we provide a formal statement and proof in Appendix [..]
Lemma 1. Up to a orthogonal change of coordinates, any linear SDE with generic drift and diffusion matrix, is equal in its law to another linear SDE with transformed drift and diagonal diffusion matrix.
3.2 The Gaussian case
We derive explicit formulae for Gaussian Schrödinger bridges with general linear reference dynamics. Our results generalise those in (5) to mvOU reference processes, and we verify in the appendix that we recover several results of (5] Theorem 3, Table 1] as special cases. Both derivations use the disintegration property (2) together with analytical expressions for the optimal coupling and bridges. Theorem 3 (Characterisation of mvOU-GSB). Consider mvOU reference process Q described by and Gaussian initial and terminal marginals at times I. = O and l. = T respectively:
Po Naa, A), pT (b,
5

---

In what follows. we write E, for the covariance of the unconditioned process mass. Define the transformed means and covariances
started at a point
1/2 TA (a m) + m)), a = ΣT
A ET e
1/2 TA
AeT
Ae
TA
1/2
ΣT b
=T 1/2 BET
and let C be the cross-covariunce term of the Gaussian entropic optimal transport plan 27 Theorem 1] between Po (aa, A) and Pr = (Nb B) with unit diffusivity, i.e.
c A 1/2 A BA 12 +
1111)
1/2
A 1/2
I.
(13)
2
Setting D. T (T T for brevity, define
2, (c- (T- - r,
B, 1,ET
1/2
c== -1
(I-e-
m
(14)
Then the Q-GSB is a Markovian Gaussian process for 0 & t & T with mean and covariance
V E1 Y 21 B.b + Ct.
(15)
att Vly, Y, = Qet + A2, 2,B + BBT 21 + BBBBB.
(16)
Furthermore, the Q-B (SBP-dyn) is described by a SDE
dY, (ir S. E (Y. v1)) dt + dB,
Y ~ N(a A).
(17)
where S. (er))))))) (1) A,A2 + 21CB3, + BBC 21 B,BB T
Importantly, show in the Appendix that recover from the results of Theorem two prior results of S - specifically, we obtain results for a Brownian and scalar OU reference process as special cases. Furthermore, in [5] the authors remark that the quantity S, determining the GSB drift, is symmetric although S is asymmetric in general. The interpretation of this result is that the GSB is driven by a time-varying potential and is thus of gradient type. By contrast, observe empirically that the mvOU-GSB drift is asymmetric when the underlying reference process Q is not of gradient
type.
3.3 The Non-Gaussian case score and flow matching
When the endpoints (pu, (r) are non- non-Gausse Gaussian, direct analytical characterisation of the solution to (SBP-dyn) is out of reach. Here we show that combining the exact characterisation of the Q- g-biriges (Section 3.1) with score and flow matching (Section 2.2) yields a simulation-free estimator of the Schrödinger bridge when Q is mvOU process (I. Specifically, we propose to first exploit the static problem (SBP-static) which can be efficiently solved via Sinkhorn-Knopp iterations [4 with an analytically tractable cost function. Next, we use the property (2) together with the characterisation of (o-bridges presented in Theorem 2 The approach of [52), which considered a Brownian reference, thus arises as a special case for A- = 0, m 0, ol.
Proposition 1. When Q is mvOU process (SPP- statich shares the same minimiser as min 2 [[(T 14 Σx ( 1T) dr (xFo, aT) H(xlpp X pr)
(18)
REIT(PO.TT)
We remind that æ + μ is affine and computation of Er relies on a one-off time integration that does not depend on the endpoints (Section A.T. Equipped with the solution to [8] we parameterise the unknown probability flow and score function of (SBP-dyn) by u, (x) and s (x)
where (0, sp) are trainable parameters. We seek to minimise the loss
L(0, sp) wang.em)-e/.orrj
E
[1] "[][ >x](sf(2) 81 (Fo. (x((0))))))))]]] (2)
(19)
In the practical case when (Po, Pr) are discrete, sampling from t, (Tp, 4r) U[0,T] & π is trivial once (18) is solved. Sampling from Pt (m).arr) and evaluating (*)(=))))) (z)(2))) amounts to invoking the results of Theorem 2 Again, computations of these quantities involve one- solution of ID matrix-valued integral: this can be solved to desired accuracy and then cached and queried so it poses negligible computational cost, and this is what we do practice. As in the Brownian case, the following result establishes the connection of the conditional score and flow matching loss (19) to (SBP-dyn), and a generalisation and restatement of results from 52 Theorem 3.2, Proposition
3.4].

---

Algorithm mvOU-OTTM: score and flow matching for mvOU-Schrödingee Bridges
Input: Samples {ar {2r from source and target distribution at times f 0, T, mvOU reference parameters (A, m, D f batch size B.
Initialise: Probability flow field ui (x), score field of (w)).
Po N Σ Ō N Ei= 60/
Form empirical marginals
sinkhorn(C, Po, reg 1.0)
Solve 18
with log-kernel cost
while not converged do (ar, ) - sample(x), [t:) sample(U[0, 7))
for & j & B do
2 - (H4, (w,, xz)) eg (()))))
Sample from mvOU-bridge (Thm.
E - L (z) ue , a ) + 11,1off (z) ,t,
(Flow, score mcotcingatching loss
end for
I = B Σド= Pj
() Step Vol, V.L.
end while
Theorem 4 (mvOU-OTFM solves .f > 0 for all (t, (99) shares the same
gradients as the unconditional loss: L(0 EEEU[, Ez-p. [llu (2) (z))( + >e (z) {.(=)][]] Furthermore, if the coupling æ solves (SBP-static) as given in Proposition | then for
(uf, s achieving the global minimum of (9)) the solution to (SBP-dyn) is given by the SDE
dX,- = (X. D (Xi) df dBB.
We stress here that the linearity assumption is only imposed on the reference process Q in (I. and not on the SBP solution P. In the general non-Gaussian case the mvOU-SBP dynamics will still be nonlinear, much in the same way as for the Brownian SBP.
Further, non-linear reference dynamics can be connected to linear dynamics of the form by linearisation about fixed point of the drift. That is, for smooth reference drift fax with mo fixed point, fxx) = O, one has f(x)= = (oz f) (xo) vo) + O(I zoll²). This naturally motivates the study of mvOU processes with A (zoo and Bo.
4 Results
Gaussian marginals: benchmarking accuracy As a first visual illustration of Theorem 3 we show in Fig. the family marginals solving (SBP-dyn) for the same pair of Gaussian marginals in dimension d 10, where the reference process Q is taken to be (a)(i) a mvOU process with high-dimensional rotational drift and (b)(i) a standard Brownian motion. While both provide valid interpolations of Po P1. it is immediately clear that the dynamics are completely different. To provide some further insight, we show in (i) and (b) the time-dependent vector field generating each bridge. According Theorem 3 both processes are time-dependent linear SDEs of the form (77). For the mvOU reference the asymmetric nature of the drift matrix is evident, while for the Brownian reference the symmetry of the drift matrix corresponds potential driven dynamics [5].
In the Gaussian case, our explicit formulae for the OU-GSB allow us empirically quantify the accuracy of our neural solver. Since existing computational methods for solving the Q-Schrödinger bridge for non-gradient Q are limited, we compare MVOU-OTFM against Iterative Proportional Maximum Likelihood (IPML) (53) and Neural Lagrangian Schrödinger Bridge (NLSB) 25 The former is based on a forward-backwarrd SDE characterisation of the Q-SB and alternates between simulation and drift estimation using Gaussian processes, while the latter uses Neural SDE frame- work. For d ranging from 2 to 100, we solve (SBP-dyn) with each method and report Table [ the average marginal eиror, measured in the Bures- Wasserstein metric 49] and the average vector field error in L In all cases, we find that MYOU-OTFM achieves the highest accuracy among all solvers considered. Among the others, NLSB is generally the most accurate followed by IPML, while BM-OTFM performs the poorest since the reference process is misspecified.
In addition to being more accurate, MVQU-OTFM is extremely fast to train owing to being simulation free, taking approximately 1-2 minutes to train on CPU for d 50, regressing the score and flow networks against (19). NLSB took 15+ minutes to train GPU requiring backpropagation through
7

---

SDE solves, while IPFP requires iterative SDE simulation. On the other hand, our approach directly exploits the exact solution to Q-bridges rather than relying on numerical approximations.
GSB marginals Pt mVOU-CS8
GSB vector field Ut
1=000
BM-G5B
0.00
0.50
0.75
1.00
Figure 2: Gaussian Schrödinger Bridges. (a) (1) Marginals of the Gaussian Schrödinger bridge (d = 10) with mvOU reference (MVOU-GSB) (ii) Time dependent vector field generating the bridge. (b) Same as (a) but for Brownian reference.
Marginal (BW))
Force (L²)
MVOU-OTFM BM-OTFM IPMI. (-)
[PML.(-) (-)
MYOU-OTFM BM-OTFM
TPML.
2
0,19:0011
8,40++.77
1.55=1.53
5.65+1.11
1.21
3.56=0.25
12.23++0.77
10.31
7,59=0.33
0.23
9.06100666
6.3415.64
3.2410998
1.1610.26
12.2710.12
7.7220.15 10.08:0.75
10
0.59 1.036
8.95
3.00
3.00.10.3
1.36100.33
3.82 0.13
12.33.10.13
10.9611.11
8.1510.25
4-222
8.81+1119
6.97422.11
5.034-0.00
1.27
12.54440.16
12.42+0.00
11.0242211
50
00.36
11,74=0.37
9.03-0021
8.32+0.63
6.39++0.13
6.25-0.17
13.26++0.15
14.14=1.00
14,40+00.2
100
6.84 0.778
15.14::0.95
16.19:1187
14.38
17.40
10.45
15.53
15.30
16.49:00.07
Table I: Marginal error measured in Bures Wasserstein metric and force error measured in L" norm
for the Schrödinger bridge learned between two Gaussian measures in varying dimension d.
mHOU-OTFM
61)
SDE
(ii0
PF-DDE
(b) (1)
SM-0TFM
(ii)
(ii)
0
22
2
100 00
60
0.5
0.5
0.0
25
25
25
Figure 3: Non-Gaussian marginals. (a)() Sampled trajectories from the mvOU-SB learned between non-Gaussian marginals in d = 10, shown in (t, zo) coordinates with score field shown in background. (a)(i, iii) Sampled SDE (stochastic) and PF-ODE (deterministic) trajectories shown in (æo, 21) coordinates, and reference drift shown in background. (b) Same as (a) but for Brownian reference.
Non- Gan-Gausiian marginals We demonstrate an application of MVOU-OTFM for general non- Gaussian marginals in Fig. where the reference dynamics are the same those in Fig. Again, the contrast between the mvOU and Brownian reference dynamics is clear, and we show the distinction between sampling of trajectories via the SDE 3 and PF-ODE formulations
Repressilator dynamics with iterated reference As application to a system of biophysical importance, we consider the repressilator [16], a model of synthetic gene circuit exhibiting oscil- latory, and hence con-cquili,rium, dynamics in cells. This system composed of three genes in loop, in which each gene represses the activity of the next (Fig. (a)). This model, implemented as a SDE [62], was also studied by [M]. which introduced a method, "Schrödinger Bridge with Iterative Reference Refinement" (SBIRR) building upon IPFP (53) to solve series of Schrödinger bridge problems whilst simultaneously learning an improved reference Q within a parametric family. In particular, they propose to alternate between solving the QBBP using IPFP as subroutine, and fitting of a parametric global drift describing Q. The same approach can be used with mvOU-OTFM, which we present in Algorithm 2 in the Appendix: we propose to apply Algorithm [ as a subroutine, and solve a regularised linear regression problem to update the mvOU reference parameters ( (A, m) at each outer step.
8

---

Repressiiator (B)()
(mvOU
(BM
[d)
Averaged leave one out error EYD
mvOU reference
Brownian reference
Learned
Figure 4: Repressilator dynamics. (a) Repressilator population snapshots and trajectories. (b) (i) Ground truth vector field. (ii, iii) Inferred multi-marginal SB vector field with (fitted mvOU, Brownian) references. (c) Ground truth linearisation of system and drift A learned by mvOU-OTFM. (d) Leave-one-ore error by iteration.
(e) Illustration of leave-one-our interpolation between two
example timepoints P:-1 (blue), Ps+1 (green) with learned mvOU reference vs. Brownian reference.
Emor metric
Iterate
Itente
Ttemte 2
Iterate 3
Iterate 4
SBIRR (mvOU)
SBIRR (MIP)
EMD
3,38 + 1.52
2.22 = 1.12
1.59 + 0.66
1.49 + 0.64
1.40 = 0.57
2.10 0,74
1.67 0.95
Energy distance
1.86 + 1.06
1.29 = 0.86
1.03 + 0.65
0.95 = 0.58
0.89 0.55
1.39 + 0.82
1.10 0.86
Table 2: Repressilator leave-one-ous interpolation error for mvOU-OTFM and SBIRR.
We use a similar setup 44 and sample snapshots for T 10 instants from the repressilator system. We show in Fig. 4 the sampled data along with example trajectories. We run Algorithm 2 for iterations and we compare in Fig. (b) the ground truth vector field to the SB vector field, both the fitted mvOU reference (at the final iterate of Alg. 2 and for a Brownian reference (as the first iterate of Alg. or equivalently the output of [22 As was found by 44 588, iterated fitting of global autonomous vector field for the reference process leverages reconstruction of dynamics better adapted the underlying system.
timepoint information and allows
We reason that the "best" mvOU process to describe the dynamics should resemble the linearisation of the repressilator system about its fixed point. Comparing the Jacobian of the repressilator system to the fitted drift matrix A in Fig. cc, shows a clear resemblance in particular the cyclic pattern of activation and inhibition is recovered. For quantitative assessment of performance, we show in Fig. dd) and Table 2 the averaged leave-one-out error for marginal interpolation task: for each remaining timepoints with timepoint t, held out. The trained
run on the model is evaluated on predicting t; from 1;-1. For both the earth mover's ' distance (EMD) and energy distance [4] we find that prediction accuracy improves with each additional iteration of Alg.
1 T, Alg.
We compare to SBIRR 44 with two choices of reference process: (i) mvOU and (ii) a general reference family parameterised by feedforward neural network. The former scenario is thus comparable to mvOU-OTFM, whilst the latter considers a wider class of reference processes. We find that both methods offer an improvement on the null reference (iteration 0 of Alg. 22.. and as expected SBIRR with a neural reference family outperforms the mvOU reference family. However, mvOU-OTFM achieves a higher accuracy than either of these methods which suggests that, while SBIRR is able to handle general reference processes and is thus more flexible, mvOU-OTEM trains more accurately. Furthermore, InOU-OTFM significantly faster, training in 1-2 minutes on CPU. By comparison, training of SBIRR requires at least several minutes with GPU acceleration.
Dynamics-resolved single cell data We next apply our framework to single-cell RNA sequencing (scRNA-see seq) dataset 4 of the cell cycle, where partial experimental quantification of the vector field is possible via metabolic labelling, i.e. for each observed cell state æ; we also have velocity estimate é; [40]. We consider N 2793 cells in d 30 PCA dimensions (Fig. 5 (la)))). For each cell

---

(a)()
Yeyomoo
Viel
You (mvOU reference) (IV) Ves (Brownien reference)
(D)
Interpolation time
Xcde
F 20
Average OFTer
23
60
102
Figure 5: Cell cycle scRNA-seq. (a) Streamlines of (i) transcriptomic vector field calculated from metabolic labelling data (ii) Fitted mvOU reference process, and Schrödinger bridge drift osn learned by OTFM with (ii) mvOU reference and (iv) Brownian reference. (b) Marginal interpolation error between first and last snapshot as a function of the reference velocity scale parameter, 7.
an experimental reading of its progression along the cell cycle is available via fluorescence, and cells are binned into T = 5 snapshots using this coordinate. From the paired state-velocity data {æ,, 6,] we fit parameters (A, m)) of mvOU process (T) by ridge regression (Fig. (ii). A limitation of "RNA velocity" " is that estimation of vector field direction is often significantly more reliable than magnitude [59]. We introduce therefore additional "scale" parameter y 0 by which to scale the mvOU drift A We reason that this amounts to a choice of matching the disparate and a priori unknown timescales of the snapshots and velocity data.
In order to select an appropriate scale, for 0 & Y & 100 we use the mvOU-GSB to interpolate between the first and last snapshots (see Fig. for illustration) and calculate the Bures- Wasserstein distance between the interpolant and each intermediate marginal (Fig. (b). From this it is clear that choosing y too small or too large leads to mvOU-GSB that does not fit well the data, while broad range of approximately 30 & 7 & 70 yields good behaviour. We fix 7 = 50 and use mvOU-OTFM to solve muulii-argiinal SBP with reference parameters (yA, m)) From Fig. 5 (iii) we see that the expected cyclic behaviour is recovered. On the other hand, using a Brownian reference (cquivalently, setting 7 0) fails to do so.
Finally, to further demonstrate robustness of our approach to the dimension d. we earry out the interpolation analysis for d = 50, 100 and provide numerical results in Table These settings of moderately high dimensionality are typical for single cell studies, where PCA projection to 20-50 dimensions is routinely used prior to any downstream analysis. In all cases, we are able to run mvOU-OTFM within minutes on CPU.
Discussion
Motivated by applications to non-equilibrium dynamics, study the Schrödinger bridge problem where the reference process is taken to be in a family of multivariate CmmstinnUUleebbee (mvOU) processes. These processes are able to capture some key behaviours of non-equilibrium systems whilst still being analytically tractable. Leveraging the tractability, we derive an exact solution of the bridge between Gaussian measures for mvOU reference. We extend our approach non-Gaussian measures using a score and flow matching approach that avoids the need for costly numerical simulation. We showcase the improvement in both accuracy and speed of our approach in variety of settings, both low and high-dimensional scenarios.
Limitations A key limitation of our method is that our analytical expressions involve matrix inversions and powers of positive definite matrices which scale as O(d3) with the dimension d. The same limitation was identified in (5) as also applying to their approach. In practice, we find that application to problems with d & 100 work well. Scaling our methods to even higher dimensional settings left to future work, and we expect that tools from the Gaussian process literature will be helpful for this. Regarding flow matching, the link between Alg. and (BPP-dyn) requires the solution to (SBP-static) to be computed exactly. In all our experiments, we achieve this by computing couplings using all available samples. In large-scale settings where the number of data points N is very large or possibly streaming, a common approach is to use minibatch couplings [50]. In such settings, the use of minibatches results in learned flows being biased. An in-depth discussion of this provided in 24), which advocates for the infrequent computation of OT couplings on large batch sizes rather than the commonly used per-iteration minibatch OT.
10

---

References
[1] Guy Abel and Nikola Sander. Quantifying global international migration flows. Science, 343(6178):1520-1222 2014.
[2] Michael S Albergo. Nicholas M Boffi, and Eric Vanden Vanden-Eijnden. Stochastic interpolants: A unifying framework for flows and diffusions. arXiv preprint arXiv:2303.08797, 20263
[3] Michael S Albergo and Eric Vanden- Dijnden. Building normalizing flows with stochastic interpolants. arXiv preprint arXiv:2209. 15571, 2022.
[4] Nico Battich, Joep Beumer, Buys De Barbanson, Lenno Krenning, Chloé S Baron, Marvin E Tanenbaum, Hans Clevers, and Alexander Van Oudenaarden. Sequencing metabolically labeled transcripts in single cells reveals mrna turnover strategies. Science, 367(6482211151-11555 2020.
[5] Charlotte Bunne, Ya-Ping Hsich, Marco Cuturi, and Andreas Krause. The schrödinger bridge between gaussian measures has closed form. In International Conference Artificial Intelligence and Statistics, pages 5802-5833. PMLR, 2023.
[6] Charlotte Bunne, Laetitia Papaxanthos, Andreas Krause, and Marco Cuturi. Proximal opti- mal transport modeling of population dynamics. In International Conference on Artificial Intelligence and Statistics, pages 6511 6511-6528. PMLR, 2022.
[7] Charlotte Bunne, Geoffrey Schiebinger, Andreas Krause, Aviv Regev, and Marco Cuturi. Optimal transport for single-cell and spatial omics. Nature Reviews Methods Primers, 4(1):58, 2024.
[8] Yongxin Chen and Tryphon Georgiou. Stochastic bridges of linear systems. IEEE Transactions on Automatic Control, 61(2):526-531, 2015.
[9] Yongxin Chen, Tryphon Georgiou, and Michele Pavon. Optimal steering of linear stochastic system to a final probability distribution, part i. IEEE Transactions on Automatic Control, 6185):115--116) 2015.
[10] Yongxin Chen, Tryphon T. Georgiou, and Michele Pavon. Optimal stecring of linear stochastic system to a final probability distribution, part ii. IEEE Transactions on Automatic Control, 61(5):1170-11888 2016.
[11] Yongxin Chen, Tryphon T Georgiou, and Michele Pavon. Optimal transport over a linear dynamical system. IEEE Transactions on Automatic Control, 6255)22137-2 2016.
[12] Yongxin Chen, Tryphon T Georgiou, and Michele Pavon. Stochastic control liaisons: Richard sinkhorn meets gaspard monge on a schrodinger bridge. Siam Review, 63(2):249-315 2021.
[13] Lénaïc Chizat, Stephen Zhang, Matthieu Heitz, and Geoffrey Schiebinger. Trajectory inference via mean-field langevin in path space. Advances in Neural Information Processing Systems, 35:16731-1677 2022.
[14] Mareo Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.
[15] Valentin De Bortoli, James Thomton, Jeremy Heng, and Arnaud Doucet. Diffusion schrödinger bridge with applications to basee generative modeling. Advances in Neural Information Processing Systems, 34: 17695-17709, 2021.
[16] Michael B Elowitz and Stanislas Leibler. A synthetic oscillatory network of transcriptional regulators. Nature, 40366767) 2000.
[17] Xiaona Fang, Karsten Kruse, Ting Lu, and Jin Wang. Nonequilibrium physics in biology. Reviews of Modern Physics, 91(4):045004, 2019.
[18] Mathieu Gilson, Enzo Tagliazucchi, and Rodrigo Cofré. Entropy production of multivariate omstein-uhlenbeck processes correlates with consciousness levels in the human brain. Physical Review E, 107(2))024121, 2023.

---

[19] Federico S Gnesotto, Federica Mura, Jannes Gladrow, and Chase P Broedersz. Broken detailed balance and non-equilibrium dynamics in living systems: a review. Reports on Progress in Physics, (166)0666 2018.
[20] Claude Godrèche and of menstein-uuhenbee processes. Journal of Physies A: Mathematical and Theoretical, 52(3):035002, 2018.
Luck. Characterising the nonequilibrium stationary states
[21] Hicham Janati. Advances in Optimal transport and applications to neuroscience. PhD thesis, Institut Polytechnique de Paris, 2021.
[22] Hicham Janati, Boris Muzellec, Gabriel Peyré, and Marco Cuturi. Entropic optimal transport between unbalanced gaussian measures has a closed form. Advances in neural information processing systems, 13:10468-10479 2020.
[23] Donald Kirk. Optimal control theory: an introduction. Courier Corporation, 2004.
[24] Michal Klein, Alireza Mousavi- Hosseini, Stephen Zhang, and Marco Cuturi. On fitting flow models with large sinkhorn couplings. arXiv preprint arXiv:2506.05522 255666 2025.
[25] Takeshi Koshizuka and Issei Sato. Neural lagrangian schr odinger bridge: Diffusion modeling for population dynamics. arXiv preprint wxXiv:220.00485: 2022.
[26] Hugo Lavenant, Stephen Zhang, Young- Young-Heon Kim, and Geoffrey Schiebinger. Towards mathematical theory of trajectory inference. arXiv preprint arXiv:2102.09204, 2021.
[27] Melvin Lax. Fluctuations from the nonequilibrium steady state. Reviews of modern physics, 32(1):25, 1960.
[28] Christian Léonard. A survey of the schrodinger problem and some of its connections with optimal transport. arXiv preprint arXiv: 1308.03000000 02000.0215.000 2013.
[29] Anders Lindquist. On feedback control of linear stochastic systems. SIAM Journal on Control, 11(2):323-343, 1973.
[30] Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXiv:2210.02747, 2022.
[31] Suryanarayana Maddu, Victor Chardès, Michael Shelley, et al. Inferring biological processes with intrinsic noise from cross-sectional data. arXiv preprint arXhv:2410.07501, 2024.
[32] Anton Mallasto, Augusto Gerolin, and Hà Quang Minh. Entropyy-regularized 2- 2-ssssrsneeinn distance between gaussian measures. Information Geometry, 3(1):289-323,203)
[33] Yuhang Mei, Mohammad Al-Jarrah, Amirhossein Taghvaei, and Yongxin Chen. Flow matching for stochastic linear control systems. arXiv preprini arXiv:2412.00617 2024.
[34] Christopher E Miles, Scott A McKinley, Fangyuan Ding, and Richard B Lehoueq. Inferring stochastic rates from heterogeneous snapshots of particle positions. Bulletin of mathematical biology, 86(6):74, 2024.
[35] Boris Muzellec, Richard Nock, Giorgio Patrini, and Frank Nielsen. Tsallis regularized opti- mal transport and ecological inference. In Proceedings of the AAAI conference on artificial intelligence, volume 31, 2017.
[36] Ingram Olkin and Friedrich Pukelsheim. The distance between two random vectors with given dispersion matrices. Linear Algebra and its Applications, 48:257-263, 1982.
[37] Juan MR Parrondo, Christian Van den Broeck, and Ryoichi Kawai. Entropy production and the arrow of time. New Journal of Physics, 11(7):073008, 2009.
[38] Katarina Petrović, Lazar Atanackovic, Kacper Kapusniak, Michael M Bronstein, Joey Bose, and Alexander Tong. Curly flow matching for learning non-gradient field dynamics. In ICLR 2025 Workshop Machine Learning for Genomics Explorations, 2025.
12

---

[39] Aram- Alexande Pooladian and Jonathan Niles-Weed. Plug-in estimation of schr odinger bridges. arXiv preprint arXiv:2408. 1686, 2024.
[40] Xiaojie Qiu, Yan Zhang, Jorge D Martin-Rufino, Chen Weng, Shayan Hosseinzadeh, Dian Yang, Angela N Pogson, Marco Y Hein, Kyung Hoi Joseph Min, Li Wang, al. Mapping transcriptomic vector fields of single cells. Cell, 185(4):600-11, 2022.
[41] Maria L Rizzo and Gábor J Székely. Energy distance. wiley interdisciplinary reviews: Compu- tational statistics, 8(1):27-38, 2016.
[42] Geoffrey Schiebinger. Reconstructing developmental landscapes and trajectories from single- cell data. Current Opinion in Systems Biology, 27:100351, 2021.
[43] Geoffrey Schiebinger, Jian Shu, Marcin Tabaka, Brian Cleary, Vidya Subramanian, Aryeh Solomon, Joshua Gould, Siyan Liu, Stacie Lin, Peter Berube, et al. Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming. Cell, 176(4):928-943, 2019.
[44] Yunyi Shen, Renato Berlinghieri, and Tamara Broderick. Multi-marginal schr)." odinger bridges with iterative reference refinement. arXiv preprint arXiv:2408. 06277, 2024.
[45] Yuyang Shi, Valentin De Bortoli, Andrew Campbell, and Arnau Doucet. Diffusion schr odinger bridge matching. arXiv preprint arXiv:2303.16852. 2023.
[46] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
[47] Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced score matching: A scalable approach to density and score estimation. In Uncertainty artificial intelligence, pages 574-584. PMLR, 2020.
[48] Andrew M Stuart and Marie-Theree Wolfram. Inverse optimal transport. SIAM Journal Applied Mathematics,
619, 2020.
[49] Asuka Takatsu. Wasserstein geometry of gaussian measures. 2011.
[50] Alexander Tong, Kilian Fatras, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector- Brooks, Guy Wolf, and Yoshua Bengio. Improving and generalizing flow-based generative models with minibatch optimal transport. arXiv preprint arXiv:2302.0044 2023.
[51] Alexander Tong, Jessie Huang, Guy Wolf, David Van Dijk, and Smita Krishnaswamy. Trajecto- rynet: A dynamic optimal transport network for modeling cellular dynamics. In International conference on machine learning, pages 9526-9536. PMLR, 2020.
[52] Alexander Tong, Nikolay Malkin, Kilian Fatras, Lazar Atanackovic, Yanlei Zhang, Guillaume Huguet, Guy Wolf, and Yoshua Bengio. Simulation schrl," odinger bridges via score and flow matching. arXiv preprint NiX.ii22377.032222. 2023.
[53] Francisco Vargas, Pierre Thodoroff, Austen Lamacraft, and Neil Lawrence. Solving schrödinger bridges via maximum likelihood. Entropy, 23(9):1144) 2021.
[54] Pat Vatiwutipong and Nattakorn Phewchean. Alternative way to derive the distribution of the multivariate ornstein -uhhinnbuhebb process. Advances in Difference Equations, 2019:1-7, 2019.
[55] Pascal Vincent. A connection between score matching and denoising autoencoders. Neural computation, 23(7):1661-1674 2011.
[56] Caleb Weinreb, Samuel Wolock, Betsabeh K Tusi, Merav Socolovsky, and Allon M Klein. Fundamental limits on dynamic inference from single-cell snapshots. Proceedings of the National Academy of Sciences, I 115(10):E2467-E2477 2018
[57] Grace Hui Ting Yeo, Sachit D Saksena, and David K Gifford. Generative modeling of single-cell time series with prescient enables prediction of cell trajectories with interventions. Nature communicationss 12(1):3222, 2021.
13

---

[58] Stephen Y Zhang, Joint trajectory and network inference via reference fitting, In Machine Learning in Computational Biology, pages 72-85. PMLR, 2024.
[59] Shijie C Zheng, Genevieve Stein-O'Brien, Leandros Boukas, Loyal A Goff, and Kasper D Hansen. Pumping the brakes ma velocity by understanding and interpreting ma velocity estimates. Genome biology, 24( (4)22)246, 2023.
[60] Ligang Zhu, Songlin Yang. Kun Zhang, Hong Wang, Xiaona Fang, and Jin Wang. Uncov- ering underlying physical principles and driving forees of cell differentiation and reprogram- ming from single-cell transcriptomics. Proceedings of the National Academy of Sciences, 121(34):e20154012,, 2024.
14

---

A Theoretical results
Quantity
Description
Introduced
A E Lddx
m Rd
Multivariate Ornstein- Uhlenbeck drift parameters
Eq.
:00
Multivariate Ornstein- Uhlenbeck diffusivity
Eq.
P
Schrüdinger bridge path measure Reference process path measure Unconditional mean & cov. of mvOU process started at Conditional density of mvOU bridge
Section
2.1
Section
142 E
Section
Pe (e
Section
mvOU bridge control between (wo, 2r)
Section
(#o.wy): (m) Ω
Conditional mean and covariance mvOU bridge
Section 3.1
Conditional covariance process of mvOU bridge
Scction 3.1
U (Gaa mnl) mm)
Conditional probability flow and score field of mvOU bridge
Section
(a,
(bb B)
Initial and terminal mvOU-GSB marginals Transformed mvOU-GSB marginals
Section
3.2
((a, A), A.,
Section 3.2
C
Cross-covarianes of entropic transport plan Key quantities for the mvOU-GSB mvOU-GSB mean process
Section 3.2
21, Br,
Section
Section
E,, E,
mvOU-GSB variance and covariance process SDE drift matrix of mvOU-GSB
Section
Scction
Table 3: Glossary of some key notations and quantities used in the statements of our theoretical results.
A.1 Some calculations multivariate
processes
For convenience, we first collect some results about multivariate Ornstein-Uhlenbeck processes of the form ( a detailed discussion of mvOU processes can be found in e.g. (54). For a lime-invariant process with drift matrix and diffusion o,, we have that
X X xo) -N(H, Ex),
where
E, = I" (t ((-T)) oo ((t-T)A dr.
=
14 8o
(20)
X6, ~ N(4, En) and writing Xoe = Xo + X₁, some tedious but straightforward applications of conditional expectations and the tower property reveal that
En
EgeLAT
Xor N(por Eod), Mor po
+ e
no.
Eos
LA Σo
Σ3 + elA E0e!AA
(21)
More generally, for a process dx, = (AX; + b)dt + odB, one has:
iA p1 e xo + (e ()(A -b), Σ, = S ( (t-s)a ToT (t-s) ds.
(22)
The expressions for the covariance of the mvOU process can be obtained from the fact that the covariance evolves following Lyapunov equation. For the case of constant coefficients, we state the following result.
Lyapunov equation solution
It is easy to verify that G. AGG G.A + Q: has solution
G ((=))A (1-S)A ds.
A.2 Bridges of multivariate Uhlenbeck processes
This problem has been studied by Chen al. in X however the material therein is geared towards a control audience and considers more general case where all coefficients are time-dependent, We will re-derive the results that we will need for processes of the form (T). Additionally, while in practice we typically use ol, state some results in the general setting of non-isotropic noise.
15

---

Derivation of the dynamical OU-bridge (Theorem Consider mvOU process
dx2 (AX, b) dl + d.,,
(23)
Now form the controlled version pinned (0, .o) and (7, ær), where (x()) is additional force arising from the conditioning of the process at an endpoint:
dx; (AX; + b (x();)) dt + dB,
(24)
The main result that we need is that (x(x) is itself a linear time-dependent field whose coefficients are independent of (ro, TT). In the above we consider the simplest case of classical optimal control, where the control is directly added to the system drift and there is no unobserved system.
Consider the system started (0, Ta) and pinned to (T æT). For now, we make all statements for time-dependent coefficients and the diffusion is non-isotropic. We then form the Lagrangian:
min E [[
{[[
[[""
|ulll'dt sup X ærl?
sup min E X
uu [ul/ + X| xr
Fixing some the resulting problem is known as linear caatarr-uaaaticc--c the control literature:
quadratic-Gaussian Gaussian problem
29
T +
X 2 2| -2r11
min EE
2 uul
dt
(25)
2
for which can write HJB equation:
VV())) + min A[Vi(( u((x)) + + [1l(((())1
(26)
subject to the terminal boundary condition Vi(()) = - =,||| and the operators A, A are defined
AA[ = (A,æ b. + ue vn.) _ v. (o.o,
(27)
A. (A,x b. V. 1
/ v (0.0, Vxf),
(28)
these are the generators of the controlled 24 and uncontrolled (23 SDEs respectively. Substituting all these in, we find that the HJBE is
(vV((=) + A.VV((x)) min (u, V Vi)
2 |luull
(29)
follows from the first order condition on the "inner" problem that u V.Vi
(=v)())) + A[V](=) [VV((2))1 IV
(30)
Use an ansatz that the value function quadratic:
1
Vi(x) 2 M, + c. de V.V((z) M,x c..
(31)
Then:
A.VK] x TT (A, (A]MM) + (A) e, + M,b,, x) (b,, c.) [tr(o,o, Mt). Now note that the quadratic form only depends on the symmetric part of the matrix, i.e.
(32)
x: A, Ma. 3 (A) M, + M,A,) xc
The HIBE is therefore
2" M. + x
(AA M, M,A,)æ + (AT ct M, x +.b trr(o,o, Mi)
(33)
(Mæ + (e) (M,
2
16

---

So
M, A, M, + M,A, M. 4 + A cz + M,b, M,C.-
(34)
(35)
with the boundary condition Mr AI and = -XTT.
Let A, M, so that A, -M, M,M, A,M, Substituting, we find that
A, = A,A, + A,A,-1 -
(36)
Actually we can rewrite the value function in a different way, which is more convenient for us:
Va(x) 1 ke) M1( ke) const
in which case c, = - k, and so we have ky ar. The corresponding ODE for k;
k A,ks + b
(37)
Let us now specialise to the case time-invariant coefficients, which allows us to write down explicit expressions for the solutions:
k Ak: ke (xT + A A A - Rewriting the SDE drift to have the form Ax + = A(x - m) we have that -Am =
(38)
A b -m, when A nonsingular. Then:
k = (T--)A (xT m) +
(39)
Now we deal with the quadratic term A,, Let G. Ar so that G - Ar--. Then G. satisfies
o,G, = G,A - AG, +I, Go = 0.
(40)
So
G 1- ds.
(41)
Substituting back, we find that
T- - ds.
A e - ds
(42)
The bridge control is therefore
Uel(2o.a7) -A, k:))
(43)
For the general case of non-constant coefficients, we express them as solutions of ODEs with terminal boundary conditions.
k A,k; bt: kT ET;
(44)
A A.A, + A,A, - I, =
(45)
In practice, both of these equations can be approximately solved by numerical integration in time. Finally we remark here that the diffusion component does not play role in the control, which is classical.
As studied by X for more general scenario of time-varying processes, the mVOU process and its conditioned versions are Gaussian processes. Under these assumptions, (X,, X,) has joint mean
"Statie" Ornstein-Uhlenber bridge statisties (Theorem 2
μ4g + p4.
(46)
and covariance
(t ()-
(47)
c((-s)A
17

---

where L is the covariance of the process started from point mass at time t:
i e ( oo (t-=)AA ds,
e = Pi-. + ele-e) B (t-8)AA
(48)
By taking the Schur complement, we have that X, X, is distributed with variance
V[x, [i]] D e
(t- (1-a)A A
D ( ((=s))
(49)
and mean
X X ,4, + e
(X, Me).
(50)
More generally, we can consider the three-point correlation (X,, Xs X.). HA + ss ① jr and the covariance is
Then the mean is
4 (+-r)a
$ (t((t-))A
(51)
Using the Schur complement again, we have that (X,, r. X x3) X has covariance
V[X, X. X,] = d e (*-r)A 4 e (t (t- & (e
(52)
For the case of time-varying coefficients, need to introduce the state transition matrix Lta which describes the deterministic aspects of evolution between $ < under (A.). That is, for dynamics d A, has z2 V t. and L (t, (t r) (r, for t > > In this case there is not an explicit expression for Φ,. Instead, let (Φ,), be the unconditional variance evolutions for such process started at 0, obtained by solving
29
Φ.
4. A, + oio, 4p= 0.
(53)
Then the more general result from g for the covariance of (X,, X., X,) is
+
$ TTT
$ 1
V
4
(54)
V
V +
From this result, the same Schur complement computation yields expressions for the conditional covariances and mean:
V[X, r, Xx|X,]
4, M T
D V 4; M 4s:
(55)
V[X,
e XX]
4. + V D
V t0
(56)
E X X,] M4. + V.T.O
Φ, 'X,, pt)
(57)
Formal statement and proof of Lemma
Lemma (Formal statement of Lemma T.. dx, AX,d + ad be a linear SDE with generic drift A Rddd and diffusion
Razdd There exists orthogonal matrix U Radd such that Y, = U X, obeys a linear SDE with transformed drift matrix UT AU and a diagonal diffusion matrix A (A), where (AI,..., >
Proof. Let D L.0oT. This is positive semidefinite and therefore can its spectral decomposition can be written D JUATU where U is orthogonal and A diag(A) A) 2 0. Consider the SDEs
dx. AX,dt + odB,
(a.)
and
dx,- = AX,dl + UAUTdB.
(b.)
Let f()) be a twice differentiable test function. SDE (a.) has generator (z)) = (Q,f(x)) Ax)
T f(z)((z))) Computing the generator for SDE (b.), (cf)(z) = (0,f(x), Ax) + ttr(UAU Ozz fe))(U)) Note that AU (z))(UU ur(o Our (())r)) Hence both SDEs share the same generator, and so equal in law.
Since orthogonal transformations of Brownian increments are also Brownian increments, we have that U dB, = dB,. Substituting into (b.), we have dx,, = AX,dt + U UdB, Then, d(UT X.) = TTAX AX,dt + AdB, - dY: = U AUYdd + AdB,
18

---

A.3 Simulation-free Schrödinger bridges with linear reference dynamics
Having characterised the solution of the SBP for general reference processes and that have derived the score and flow for the Ornstein-Uhlenbeck bridge, the path clear towards a simulation- free scheme for learning Schrödinger bridges where the reference dynamics are given by a linear SDE.
For solution of the static SBP problem 18 (equivalently, SBP-static),, the transition kernel of the reference dynamics is given by
dQor x1) dæ & dæ
dQo dæ
Ze exp
ax 11:00 E (x, 1400)
(zo)
dQu:
dQo
log
2 ar
1 80)
T "x (a) 14()) + log Zt - log
(aro)
(58)
dx & dæ
dx
where 14: denotes the mean time / conditional on (0, mg). Notably, the last two terms depend only zo- It is classical result (and very easy to show) that these kinds of terms do not affect the minimiser of (T8) and so they are immaterial. The cost function use thus effectively
1 ClEp, xt) (())) 2 (x (,
(x. E
(x;
(59)
Once a solution ा of the static SBP is on hand, we want to utilise the stochastic regression objective on the conditional score and flow. Sampling (To. æx) from π, recall that ((m..m)) Q (#a so we want sample from the Q-bridge using [II). [2]:
1 ~~ U1O, 1].
(x0.a1)
NNHe(( Ee((ar.as:))
(60)
Equations (9) and 12 give us the score and flow æz.
A.4
Characterisation of the Q-GSB (Theorem
We will construct the Q-GSB utilising its characterisation Our approach is similar to that of S we first obtain explicit formulae for the static Q-GSB (SBP-static) and then build towards the dynamical Q-GSB (SBP-dyn) by using the characterisation of Q-bridges.
Static Q-GSB The standard Gaussian EOT problem has well known solution 32 Here we use the notations of [21] and define
OT (a, (a, B) == min rell(a.3) 2 x 10-2 y1l2 d v) o² H( & 3),
(61)
where a > 0 the regularisation level and A(aa A) and B (b,, B). The solution the Gaussian case is given by [1 Theorem 1]
61 in
1 _ lla bl2 1
1 B (A, B)
OT O (a,
(62)
2
1 (=)+)))
_ ( B) tr(A) tr(B) (t.(()) + o² 22 det
(63)
1/2 BA + g² + A 1/2
1/2
σ²
C 1/2
2 - I.
(64)
We are, of course, interested in slightly different problem, namely, for Po Naa, A) and pr (b, B), we seek to solve minzell(po.rr) KL(T (o. Expanding the definition of KL and
19

---

simplifying, we get
dm dQoor
SBP-statich
min xelT(po.Pr)
dr log
(65)
-s dt log min xe(l(
dQ2OT
+ H( @ + H(Poldz)) + H drr)
dzzg @ dær
(66)
S dr log
dQo.T
min =E
+ H(
(67)
dæo @ dær
-f dr( UT) log
(dQrooFzTZ0o) H(T PP)
min (P))
(68)
der
where by = we denote equality of the objective up to additive constants which do not affect the minimiser We note that the last line is exactly (61) with a, = Pos PT and and modified cost. Further, recognise dQrroo(<> or Thus.
(2)l /ro///r as the Q-transition density, and Qr10 = N( 2o Er).
1 rell(po.pr) 2 where we remind that mo is the mean at time conditional on starting at (0, zo). and Σ is the covariance at time started at point mass, i.c. Eg- 0.
[ (IT HT
SBP-static)
min
(Z En (0r -TT 20) du Zv. ②T) + H T Po X pr) (69)
In particular, have E, +. using the notation of Theorem Recall that in the case of constant coefficients, we have the following explicit expressions:
H 2b e e'A( m) + m,
(70)
$ (t-s))ard oo (t( ds.
(71)
Define the following changes of variable (Io, xT) - (T.. TT), injective whenever Ey is positive definite:
To To(wo) Σn 1/2 azo PT
Er Tr(xT)
ΣT A/T
# ग and note that the relative entropy is Let Po = (To) and PT (Tr) and ग = (To x Ty
invariant this change coordinates. So the problem 69 is equivalent to
1
min TE 2
Tr - Toll du FT
H(T & Pr).
(72)
This exactly 8 (Po: PT) as per (6] with Pa
- means and covariances are
M( A) and Pr Nb, B). The transformed
a = Σ
1/2 TA (a - m) + m)
A = ET
1/2TA AE TA EnT 1/2
b Σ 1/2
B = ΣТ BET
We remark that for time-dependent coefficients approach still applies, however computation of the transformation coefficients will be terms of the general state transition matrix instead of matrix exponentials. In what follows, focus results for constant coefficients for simplicity and their practical relevance.
The solution to 72 is therefore
c
(wo, 2T) ~ N
a
B
where C is the cross-covariance the transport plan, given by (64)
20

---

Dynamic Q-GSB: marginals and covariance Let (z FT) denote the Q-bridge pinned at (T). Then by Theorem 2
(ao 2T) = =N(H E(((eo N (M) Ste).
Invoking the inverse mappings , ,T and substituting the expression for Ma (x from 2 we get
e (7-1)A
ri
1/2
कτ +
[I (I-e-(T-1)
(73)
21
2.
20 + BB +
(74)
where r, (Te(T-1)4 T For clarity, we state first general formula:
no [M]]
Eoo Σio
Eoi
)
Lemma 2. Let (Xo, X1) ~ N
Let
Σ11
Yxo.x, N(AX, + BX C2S.
Then, [| A4 + Bug + c and V[Y] 52+ A + BEnB' + AE BE AT.
Application to (74 gives us the following formula for the variance of the bridge at time I:
V v[,,] Ω AM B,BB, A,B. B.C 21 =: E.
(75)
Substituting and simplifying, we get the following expressions for each of the terms:
T - - (1) ) A Ac (I - (-
(76)
B, B,B3, T,Br
(77)
CBB (c- p E 1/2 C T
(78)
Similarly, the mean can be computed as
E [[X,] 2.. + BB Ct
(79)
(e
c (a- m) + r. m) +
Dt-
(80)
To calculate the covariance of the SB process, we use again the disintegration property from (52) for the Q-bridges and () < s <<st< Ti
We have
v[x,, X, Xr] (1=
1 0 AT Φ ( (T-)) =: n.
(81)
Now,
Vxx,, X = B[(X, - - -v)
(82)
E[X,X, 1 V.u.
(83)
E| X X,X, Xr v.v.
(84)
and
E|X,X,1 Xg. X = V(X,, X, X,. X MoliXo. Hr) in which the first (variance) term doesn 't actually depend (Xn, XT). So,
(85)
E[X,X,1] N.t Ex xx) [Ms][ Xx) (x(o
(86)
In fact let's switch to the "mapped" coordinates X Xx for the endpoints. We abuse notation and omit the inverse map in what follows. Then, from previously,
M2((x,xx)) = X, + B,Xr + Ct.
(87)
Expanding, collecting and cancelling terms, we have that
(X,,,xx) (| X,
V.
(88)
21
X |2. +
22
X + B, X3 + B.V
18
(89)
2 A 22 B.C 21 BBBB.
(90)
Putting everything together, we get
V[,,, S2t A2,
C.B + B.C 21 B.BB. =: =:
911)
21

---

Dynamic Q-GSB: SDE representation Expanding the covariance (TD. we find that
We proceed via the generator route also used by [5].
V[ X X+n] = Bi,thh = Sl,c+ 421 2hh + 2I. 1hh BBC 21th + B,BB, 1th (92)
E + hi (c. A2, + 2,, + BBC a, + 2,,B
(93)
E + h. + ...
(94)
where we have set the key quantity
T. S. = (e.: (r,r) (1) + 21 A., +
2 + B.C 21, + B,B
T
(95)
Now, (Xe+h, X²) are jointly Gaussian with mean and covariance
Vi+h V:
E 1ith
+ E
(96)
So
[[X++1|Xi = 1ith - E 3 (X, vi)
(97)
= =X, + h h + S. =. ' - vi)
(98)
v[X++|Xx1] = - ET
-
(99)
= -hEE, - (S. s!
(100)
Let x - u())) be twice differentiable test function. Then:
E [u((x+n)) = x]
(101)
= E (u) + (dau( (x,((()))))))))) x.)))(X++a: - X.)+
(102)
2 (X x.) T622 (X))) - X1) X, = x
= =u(x) + h(C,,(x)))
(iz, + S. و (a v,)) +
1
(103)
_ 2
(Xi+h X,) (o (X, (X+th - X1) X;
- =((x)+ + ( h(e)u(x))) (i: + S. (a v1)) + h
2 1r(l (E, (S. + S,1) +
(104)
Subtracting u(æ)) and taking the limit h - 0, it is clear from the definition of the generator that the SDE drift is
i S. _. (x-1x))
(105)
This takes the same form as the equation found in [. however our formulae for S, allow us to apply it to any linear reference SDE, not necessarily ones with scalar drift. Additionally, we empirically verify that for asymmetric A the matrix S, is asymmetric. This contrasts with the symmetric nature of the drift for the gradient-type setting.
Now we want to work out what each of these terms are practice. Effectively we need to compute 2, B, and (2,.St,)))))
21 = (Ae- -('--4 - T,
(106)
B, 1/
(107)
r, (c ( TATA A,A CT Φ-
(108)
(2,5),))(1) = $, {1{AT-
(T (T-))A
(T(e(T-1) (-AD, + e oo
e'A ooTT
(109)
22

---

Remark on the case of time-varying coefficients In this case we can still assume without loss of generality that we still work on the time interval [0, by shifting the time coordinate if necessary. Let (4, s) be the state transition matrix associated with A, and e be the solution at time / to (53). Then we still have Ey & for the transition kernel covariance the cost, i.e. the covariance started from point mass. For the mean of the reference process started from Zo. we have the generalised formula
(r())) (),00) S: ((, ds,
(110)
and we remark that its inverse po² always exists since & is never singular:
1 (y) (1, -1 ( [ v (((,,)).m. dss)
(111)
Then
0 E/11/ pTT
(112)
A
(
(113)
b Σ 1/2b.
(114)
Σ 1/2 BET 1/2
(115)
Recall that the general expression for the Qrbiiege conditioned on (0, æo). (T, æT) is
[l((=mm11) = (zo)) + o V (T, t)
- (x)))
(116)
Letting Ty MT (E) To) x0)) and øT Σ/ Tr and substituting the expression for pt
we have
MUIT (zo)), = (v 2 10 (T, (T+ ) Σ 1 2 To
+ p v (T, t)
(117)
+ I"" (1, ((1,))Am,,l. A.m,
T
Set r, = +, (T, T7' then
2, (I, [.) E/
(118)
B, 1,E// 1/2
(119)
- J. ((,,s)
c
(120)
Case: Brownian motion Let's check this against the results of Bunne and Hsieh [5] who consider a class of reference processes
dY, (c.Y2 + az) dt + gs dB,.
This corresponds to a special case of ours, where the drift is scalar-valued. Setting St = w. 0z =
0, their results give the marginal parameters of the GSB connecting Nu,o En) and M(m, E1). In what follows, using their results and notations of [ Table 1], we have that r; t, F2 = 1 - t, ('), wit, Pt= t. Then:
Ht = Fi44 + T1 + ((t) - r(((),
(121)
(1 t) tpri,
(122)
Σ, -, [3 + +rT C.) + x(t,t)(1 pc)1
(123)
= (1 1)² 1 + (11 (),,) + C..) + w22(1 I)I.
(124)
Let us compute the same quantities using our formulae. We set A 0, m = and T 1. Then,
11 xو = wI. Then,
A= A
B = - BB
23

---

Also from the underbraced expressions, we have
21, w(1 I)I,
2, wtI
Assembling our results, we get
E 1 (I r2) a + r,b
= (1 - 1-)aaa- + tb
v| Ω2 2 + B, 2.UB C,C
S2 + + (1-11) A 1 - t) C
Now C is the transport plan obtained from the scaled input Gaussians, of variance A = w//wwBBBB B/W² (in the case of Brownian reference process, there is no linear mapping via the flow map). For the scaled measures, EOT is calculated with E o² 1. It stands to reason that once we scale everything back, we should get C C/w². To verify this, recall that for the untransformed problem
arg min 1 1 x-- dT H( arg min xEll(a,
--
1 2
| wll dm + w H(T @ 3).
Then for the RHS problem, we have that
WA
w
1/2 BA
1/2
I) 1/122
1/2 13/2A-11:
I
2
1/2
+
1
1/1/2
1/2
I)
C²
Also it is easy to verify that Qr Σ ((g.00) = w2((1 I)I.
v[x,] - 22 t)I ()- fBB + 1(1 - () + CT),
(125)
where is the EOT plan between the unscaled measures A). B3 with E w², This is exactly the same result as the one derived in [3].
Case: Centered scalar OU process A AI, m From 15 Table 1], we have that 7. sinh(At)/sinh(M) (A), F. = sinh(At) coth(At) - sinh(At) coth(A), Ç = 0. The mean is then
pr T:Mo T.
(126)
With ((t,t)) = wee At sinh(Xt)/X A and P - A(1-t) sinh sinh(A), the variance is
ES E. + + + (C. + C..)) + x(t, t) - P.)
(127)
where x²2 = w² /in For convenience, we will check things term by term and using Mathe- matica.
Check first the mean. We have:
22r
1 c1-c-2A6
-22t
P. e-=
I
(128)
21
2-2
Then, plugging into our formula,
E[X.J - (ell-oo)
at r b.
(129)
Now the variance. We have that
( (i>sin((A)) I,
Σ, w²
(130)
-
A
A = aA,
(131)
w sinh((A)
B We à sinh())
B BBBB
(132)
2, (ca-D)A - T E|!²
(133)
B. = 1/2
(134)
24

---

From here it is straightforward to verify that
21,A2 A
B,BB, =riB.
(135)
Now note that
c= = A 1/2
1/2
BA
1/2
11)
1
1/2
1/2
I
1
(136)
2
(ag)
1/2
(a.3) 1//2
= aB
A
BA
1)
A-
(137)
2
VapC 1/va3
(138)
where C Cl/va denotes the EOT covariance for entropic regularisation level 1/vaB. Then:
C.B. + BBC 21 = (e~01-1) 1 r,) 3. 1/var + C cl/va))
(139)
It is also easy to check that
w (An
vaB
X
Finally,
S2 = Ee (=)om)) = 21 -
2A(1
0 x(t, t) P.)
(140)
We have verified that all the terms in the expression for the variance agree.
A.5
Proof of Theorem
Proof. We follow the arguments of [52] with application to the generalised Schrödinger bridge not much changes for Q as the reference and we reproduce detail the arguments as follows. For the equality of gradients, it suffices to show this for the flow matching component of the loss, since the score matching component can be handled using the exact same arguments conditional flow between (Fo, TT) and u, be the SB flow defined by
[22 Let Ue((2sg.rr) be
Pe|[mng (x)
uz(x) = (x
(x) (
Pr(x) = E ((g Pt((mo) (x).
Let ui be the neural flow approximation with parameter 0. Assuming that Pr (x > 0 for all (t, x) we then have:
Volk (eg E. Be (m)( (=](=))))r))))]1 Vol uu (x))) (141)
(o En Exe my [|(x(())) [| (xx) 2( (z), Ue (e|( (x)) (142)
-oR (z)) (ze) [x)] + [=)]2 2 (u,(2)) un (2:) Pelin =y) (u, (z), ux())) (u) (= (x)))
(143)
(144)
where in the last line terms not depending on are zero and we use the fact that E Be(( Fn fxx) Now, using the relation between us and u, (T)
E Ex)m(() (u, uz(x:)) ddp( (u, uc()))
(145)
ddp((a)) ( (()())) f di((wo, aT) Pr pr (x(x)
(146)
-f dm ØT) S dæ (x) (.(xx) Pel(aoo (x)
Uel(220.00)) (x)
pr(ar)
(147)
= -f d d xT) (pilkro ( (uf(a)) (z), (x((m)))rr (x))
(148)
(o n Ex E 2r) (x)) ( 2e
(149)
Peikon
25

---

So we conclude that the two gradients are equal. Clearly, the unconditional loss can be rewritten over candidate flow and score fields (ú, 8) as
8) li u 11 (ddpo())) + || (dpp
where denote by [| (2dd [o, 11 x Rd - Rd This loss is zero iff i u and á = s (dpr x da)msmmost everywhere.
St dt (1)
for test function h,
Let P denote the law of the Schrödinger bridge per dyBP and write Pa(x) mean its marginal at time I. Then in ( and for Q mvOU process (I). identifying:
dlP OT = where x is prescribed in Proposition 1 and
Q 2 Pel( (æ) where Pr((20.0r) is defined in Theorem 2
substituting all these into one has
pr(x)) = d xy DT))
Consequently, there exists ux(x) such that (ap)(a) = -v (p)()))(())). Writing (x V log Pa (z) to be the score, recognising terms from the probability flow ODE, follows that that the SDE
dx, = (u,( + Ds, odB.
(150)
generates the marginals (x) of the Schrödinger bridge. Since characterised as a mixture of Q-bridges, it follows that X, defined by the SDE (150) generates the Markovianisation of P (see
Appendix of [52). Moreover, is the unique process that both Markov and mixture of
Q-bridges 228 45
A.6 Additional background
Hamiltnn-Jacoi-Bellaan equation Consider the deterministic control problem on [0, 77:
Vo(zo) min
C(xr, we) di D(TT).
(151)
Jo
subject a dynamics é; F( ue). V((x) is the value function for the problem starting at (t, æ) up to the final time T', The function ør + D(xT) specifics a cost on the final value. The corresponding Hamilton-Jacobi Macobi-Belllan (HJB) equation 23 Section 3.11] is
a Vz() + min (V.Vi(z), F(r, ux(x))) + C(a, (x))) = 0,
(152)
subject to the final boundary condition Vr(x) = D(x). A heuristic derivation is as follows. Consid- ering time interval (t, + öt) and a path Te, it's clear that
++t
Clr,, u,)ds + V Vetoa(zeta))
Va(a.) min
(153)
Using Taylor expansion and the constraint 2; = F un) we have to leading order that
V++81(Fe+ V((x;)) + (d.V((x.) <V, Varre, F(x,, mz))) 81 + O
(154)
Substituting back and also approximating the integral, we find that
Vi(x,) min Clz,, ut) V((x,) (&.Vi(x.) + F.Vi(zs),FFr, uc))) öt
(155)
Cancelling terms, rearranging and taking limit öt 0, get the desired result.
26

---

Stochastic saaniltonnJacob--elllnen case X, driven by an SDE of the form
Now we consider the stochastic variant of the HJB, in which
dX,= F(X,, we) di + o2 dB..
(156)
The value function thus in expectation:
Vo(z)) min I
u)) dt d(
(157)
J0
Carrying out the same expansion as before, we have
(x(2 min E
(x),, + Vetor Vetor(Xe+00)
(158)
Note that in the above, (x X = xz) is a random variable and hence so is Vi+ (XX+6t)))) which is why it appears in the expectation. Using Itô's formula to expand this, we get
Vi(X,) + avv(X))) v (a:o) V.V.())))) de (V2V((X)).
Vi(X,) + (VV(()+) V (a,0, V.V(())))) + (V2V(()), (VVV(X)) ui)) dt Ve+sel (ee+tt)
2
2
+ (Vee((X)) o,dBi).
Plugging this in, cancelling terms, and noting that the final term has zero expectation, we find that 0 2 (X,,) + min (AA[V)(X,,, me) + C(X,,4))
(159)
where A is the generator of the SDE governing Xx, i.c.
Aff u) (F(x, u). v.f(=))+ + V (o1o, Vzf).
(160)
The HJB equation for stochastic control problem is therefore
0 = de V(r) + min (< uxxx) + C(r, uz()))
(161)
iee this is the same as the HJB for the deterministic case, except with a diffusive term arising from the stochasticity.
B Experiment details
For (mvOU, BM)-OTFM and IPFP, all computations carried out by CPU (8x Intel Xeon Gold 6254). NL.B and SBIRR computations were accelerated using a single NVIDIA L40S GPU. Code is available at httph: / htt///gttttttt higtt........ www wwwp:////////ete/////////// _BPP
B.1 Gaussian benchmarking
We consider settings of dimension d {2, 5, 10, 25, 50) with data consisting of an initial and terminal Gaussian distribution and a mvOU reference process constructed as follows. We sample d x 2 submatrix U from a random x orthogonal matrix. Then for the reference process, A is constructed as the d x d matrix
A U
0
0
UT.
225
Similarly, we let m = [1 -1j
For the marginals, we take N(uo, Eo) and N( En) where
140 U (-5 -0.5]
Σo U
0.1 0.005
UT + O.II.
0.005
0.1
1.1
-2
#1 U [0.5
2.5]
Σ =
U + 0.II.
-2
1.1
For each d, the initial and final marginals are approximated by = 128 samples { N(H4, Eo), {of AN ~ NW,, Σ,). We fix o = and compute the exact marginal parameters Ex)eeo using the formulas of Theorem and reference parameters (A, m). Additionally, compute the SDE drift (sp æ) of the mvOU-GSB using (17)
27

---

(mvOU, BM)-OTFM We apply mvOU-OTFM as per Algorithm 1 to learn a conditional flow matching approximation the Schrödinger bridge. We choose parameterise the probability flow and score fields using two feed-forward neural networks uo(t, x) NNo + 1, (t, a)) and sp(t, m) (d 1, (t, ..), each with [64, 64, 64] hidden dimensions and ReLU activations. We use batch size of 64 and learning rate 10 for 2,500 iterations using the AdamW optimiser. For BM-OTFM, the same training procedure was used except the reference process was taken to be Brownian motion with unit diffusivity.
IPFP We use the IPFP implementation provided by the authors of [53]. specifically using the variant of their algorithm based on Gaussian process approximations to the Schrödinger bridge drift [33] Algorithm 2]. We provide the mvOU reference process as the prior drift function, i.e. a - A(x m) and employ an exponential kernel for the Gaussian process vector field approximation. We run the IPFP algorithm for 10 iterations, which is double that used in the original publication. Since IPFP explicitly constructs forward and reverse processes, at each time 0 & L & 1, IPFP outputs two estimates of the Schrödinger bridge marginal, one from each process. We consider both outputs and distinguish between them using the (-,+)) employ only the forward drift.
-)) symbols in Table 1 For the vector field estimate we
NL.B We use the NLSB implementation provided by the authors of [25] Following Section 4.6]. for mvOU reference process dx,, ft dt + a d,,,, the generalised SB problem (SBP-dyn) can be rewritten as a stochastic control problem
[
min E the
uu 1 dt
subject dX, (fr u.) dt + a dB..
On the other hand, the Lagrangian SB problem [25 Definition 3.1] is
[Ta L(t, x. (x)))d subject to d.,, Ur dl + dB,.
min E Mts
This shows that the appropriate Lagrangian to use is L(t, 22 vx)) = - We apply
this in the mvOU setting by taking f(()) = A(x m). The NLSB approach backpropagate through solution of the SDE to directly learn a neural approximation of the SB drift v. We train NLSB the same hyperparameters as used in the original publication, using the Adam optimiser with learning rate 10- for total of 2,500 epochs.
Metrics We measure both the approximation error of Gaussian Schrödinger bridge marginals as well as the error in vector field estimation. For the marginal reconstruction, for each method and for each 0 t & 1 we sample points integrated forward in time and compute the mean f4, and variance E, from samples. We then compare the ground truth marginals (1u; Σ,) computed from exact formulas using the Bures Wasserstein metric (62). For the vector field, for each method and for each time 0 & / & 1 we sample M 1024 points from the ground truth SB marginal N N(1; Σ,) and empirically estimate | - usBllN((,..)) from samples. All experiments were repeated over 5 independent runs and summary statistics for metries are shown in Table [
B.2 Gaussian mixture example
We consider initial and terminal marginals sampled from Gaussian mixtures. At time 0, we sample from
N (U -0.5] ...IUU)). N ( 0.5] 0.0625UUU)
in 1:1 ratio, and at time t = 1 we sample in the same fashion from
N 2 2.5] 0.010U)), N (U12.5 2.5] 0.25UU7).
Here, the matrix U and reference process parameters (A, m) are the same as used for the previous Gaussian example for d = 10. All other training details are the same as for the Gaussian example.
28

---

B.3 Repressilator example
Based on the system studied in
we simulate stochastic trajectories from the system
dx
B
dt d
d/
1
k))
701
10
dr
B
dl dB( 922
k
n
3
1
(162)
df
1 (I1/ki)²
7
1
dr3
B
dt dB(
0.1
dt
1+ (I2/k:)"
with initial condition (V[[1, 1, 2]. 0.01.) and simulated for the time interval é [0, 10] using the Euler-Maruyama discretisation. Snapshots were sampled at T = 10 time points evenly spaced on [o, 10], each comprising of 100 samples.
We employ Algorithm 2 to this data, running for iterations starting from an initial Brownian reference process A m O. For the inner loop running mvOU-OTFM (Algorithm D parameterise the probability flow and score as in the Gaussian example, with feed-forward networks of hidden dimension [64, 64, 64] and ReLU activations. We run mvOU-OTFM for 1, 1,000 iterations using the AdamW optimiser, a batch size of 64 and a learning rate of 10-2 At each step of the outer loop. we employ ridge regression fit the updated mvOU reference parameters (A, m)) We do this using the standard RidgeCV method implemented in the scikit-lean package, which automatically selects the regularisation parameter.
We also carry out hold-one-out runs where for 2 & i & 9 (i.e. all time-points except for the very first and last), Algorithm ii applied to T - 1 snapshots with the snapshot at t, held out. Once the mvOU reference parameters are learned, forward integration marginal at We report the reconstruction error in terms of the earth-mover distance (EMD) and energy distance [41] Table shows results averaged over held-out timepoints, and full results (split by timepoint) are shown in Table
the learned mvOU-SB is used to predict the
Since the system marginals are unimodal we reason that they can be reasonably well approximated by Gaussians. For each 2 < i & 9 we fit multivariate Gaussians to the snapshots f;-1, t, f and use the results of Theorem 3 with the fitted mvOU reference output by Algorithm 2 to solve the mvOU-GSB between Pt. Pt.... This is illustrated in Figure 4 for i 4, and we show the full results for all timepoints in Figure 0 comparison with the standard Brownian GSB.
We apply SBIRR 481 using the implementation provided with the original publication, which notably includes an improved implementation of IPFP (53 that utilises GPU acceleration. We provide the reference vector field A(x m. and seek to learn reference process in of two families: (1) mvOU processes, i.e. we consider the family of reference drifts Á(x m) where A, m are to be fit, or (2) general drifts, i.c. we parameterise the drift using a feed-forward neural network with hidden dimensions [64, 64, 64]. For each choice of reference family,
SBIRR
Algorithm 1 of for outer iterations and 10 inner IPFP iterations, as was also done in the original paper.
B.4 Cell cycle scRNA-seq
The metabolic labelled cell cycle dataset of 4 is obtained and preprocessed following the tutorial available with the Dynamo 40) package. This gives dataset of N 2, 793 cells, embedded in 30 PCA dimensions. In addition to transcriptional state (o.). Dynamo uses metabolic labelling data to predict the transcriptional velocity [6 for each cell.
To fit the reference process parameters (A, m)) method in scikit-learn. We train mvOU-OTFM using [ with = 0.3, parameterising the probability flow and score as previously using feed-forward networks of hidden dimensions [64, 64, 64] and train with a batch size of 64, learning rate of 10-2 for a total of 1,000 iterations.
again ridge regression via the RidgeCV
29

---

mv3U reference
Brewnian reference
mVGU reference
Brownian referance
A-2
(a) i =
(b) =
Brownion reference
-1
-2
-1
(c) i 3
(d) 4
mvOU reforenee
Brewnian reforence
mvVU raference
Brownion reference
R13
Rvs
(c) 5
i=
Brownian
-1
i=
Figure 6: Repressilator mvOU-GSB interpolation. Using the learned mvOU reference process, we interpolate between P;-1 (blue) and p:-1 (green). Middle timepoint is shown in red.
Omsttenn Uhitebee reference
Wiener reference
5.3
5.0
2.5
-25
25
-50
5.0
10
Figure 7: Cell cycle mvOU-GSB interpolation. Using the learned mvOU reference process and scale factor 7 50, we interpolate between the first snapshot P1 (blue) and last snapshot pr (green). All computations are done in d = 30 and shown in leading PCs.
30

---

Scale factor
0.0
10.0
20.0
30.0
40.0
50.0
60.0
70.0
80.0
90.0
100.0
0.25
40.18
30.8
22.55
16.34
12.77
11.07
10.26
10.11
10.69
11.93
13.65
d
50
0.5
28.59
18.47
11.53
9.62
12.36
15.32
12.52
10.64
11.02
13.96
0.75
15.31
10.31
7.47
7.06
8.18
8.9
8.5
7.9
8.09
9.42
11.92
0.25
44.49
36.82
29.97
24.19
19.99
17.48
16.13
15.29
14.78
14.7
15.18
100
0.5
32.48
23.75
17.12
13.55
13.79
16.84
19.78
20.21
18.25
15.7
14.22
0.75
19.33
14.95
12
10.74
11.19
12.68
13.98
14.36
14.03
13.69
13.82
Table 4: Single cell interpolation results for d = 50, 100.
Algorithm 2 Iterated reference fitting with mvOU-OTFM
Input:
Samples {ar ]==1
N from multiple snapshots times [t,)]=1- initial mvOU reference parameters
(A, m), diffusivity D _00
Initialise:
Probability flow field un (ax). score field of (æ)) FitReferencee
Define:
V) :- minAA. IV A(X m)12 \|A Form empirical marginals
+
Pj N - 8 9. $
while converged do
(w $6 < fitOTFM (mm (T)....))))) 0. - (ul, + Daf, (a. 1 & N,, EEjET
Fit flow and score with reference parameters
SDE drift
A. m FitReference
1
Update reference parameters
end while
Leave-one-ont marginal interpolation
Error metric
Iterate 0
Iterate
Iterate 2
Iterate
Inerate 4
SBIRR (mvOU)
SBIRR (MLP)
3.59
0.17
3.14
0.12
2.28
0.15
2.08
0.11
2,02
: 0.11
2.24 : 0.28
2.71
0.41
5.20
0.47
2.59
1 0.29
1.62
0.28
1.27
0.25
1.13
0.16
3.13 1 0.62
1.29
0.92
3
3.23
0.24
1.42 0.18
1.10
0.14
0.86
0.08
0,83
0.10
2.67
+ 0.85
1.39
0.55
4
1.48
0.20
0.52
1 0.05
0.47
0.05
0.47
1 0.06
0.48
1 0.06
1.38 1 0.46
0.94 1 0.28
EMD
0.65
1.29
0.71
6
6.18
0.41
3.42
2.18 118-00 0.28
1.9 0.38
1.75 + 0.17
2.40 + 0.32
1.96 1.41
2.56
+ 0.25
3.09
3.09
2.13 0.46
2.23 0.55
1.93
1.45 -
0.51 0.13
2.29
0.26
2.12
1.22 0.33
181
181
0.16
1.93 0.6
2.14 0.40
3.00
1 0.09
2.15 0.D7
2,24 0.10
210 D.D7
2.05 1 0.08
0.18
I.56 0.26
2
3.53
0.21
0.19
1.66 0.20
1.38 0.21
1.26 +
2.61 0.33
2.05 0.61
3
2.29
+ 0.17
0.11
1.00-.01 0.10
0.80 0.06
0.78 + 0.07
2.10 + 0.52
1.25 0.43
Energy
0.06
0.04
5
1.10 0.18
0.27
0.45 155-0. 0.13
0.55 0.04
0.51 0.07
0.64 0.10
0.66 0.37
6
2.49 0.12
1.12
1.13
1.01 0.22
0.91 0.11
1.43 0.17
0.90 0.64
7
0.97 0.13
1,39 0.65
1.02
0.1
1.5 0.25
0.25
0.65 1 0.09
0.17 0.08
8
0.54
+ 0.11
0.57 4 0.03
0.56 0.11
0.55 0.10
0.55 + 0.04
0.48 0.08
0.36 0.06
Table 5: Full results for repressilator example
31

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The claims made in the abstract and introduction accurately reflect this paper's contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes].
Justification: We discuss the limitations of our work, specifically on the linearity of the reference process and the computational cost of our algorithm.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiatiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
32

---

Answer: [Yes]
Justification:
We provide all assumptions and complete proofs for all theoretical results
stated.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulae, and proofs in the paper should be numbered and cross referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: In our paper we provide all information needed to reproduce our main experi- mental results.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
33

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We provide data and code with instructions to reproduce our experimental results in our supplemental material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips... cc/
Poliiyllllyyyoyylo for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
//nips.ccppblle httpc///uuidessoeeeuubbiiiis
(CodeSubiisss wwtauubmmsssooni ssiiiioiiccy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: These details are all available in supplemental material.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:
[Yes]
Justification:
We provide error bars for the main experimental results of our paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for
example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
34

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
show in tables or
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: These details are available in our supplemental material.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer:
[Yes]
Justification:
All research in this paper conforms with the NeurIPS code of ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- cration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: Our work is foundational and we do not foresee any direct path to any societal implications arising from our work.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
35

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All assets used in our research are appropriately credited.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in package should be provided. For popular datasets, paperswithcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
36

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer:
[Yes]
Justification:
We provide the full code and data in our supplementary materials.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Not applicable.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Not applicable.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
37

---

16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: No LLM usage towards this paper.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to L.LM policy (https: for what should or should not be described.
38


Reviews:
[{"review_id": "Awc7KlZktp", "reviewer": "Reviewer_xnYH", "summary": "The paper considers a special case of Schroedinger's bridge where marginals are Gaussian and the dynamics are an OU process instead of Brownian motion. The method is developed and demonstrated on some example data sets.", "questions": "I appreciate this is standard in the SBP literature, but I'm always a bit confused as to why it's useful to look at just over 2000 genes, with their expression compressed to 30 PCs. I suspect it's more likely about getting scaling on the data set to work. It would be great to see a larger example that feels a little more plausible.\n\nI didn't quite understand why a larger data set wasn't run at the end. I appreciate the comparisons are with methods that are computationally expensive. But can't you go much larger with this approach? If not why not? Would be good to clarify. I wold have gone higher for significance if we'd been able to see such an example", "limitations": "yes", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The paper does a good job of connecting to the physics literature, and mapping the ideas onto language that is easily understood by machine learning researchers.\n\nIn a relatively compressed space the paper does a good job of covering the related work and how the proposed algorithm fits into the wider landscape.\n\nAgain, given the challenges of a technical exposition, the paper does a good job of giving confidence in the method in the main text while providing additional resources in appendices and references.\n\nThe experiments show some illustrative examples and also a larger example on gene expression data.", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time to assess our work and for your positive evaluation and suggestions!\n\n> The paper considers a special case of Schroedinger's bridge where marginals are Gaussian and the dynamics are an OU process instead of Brownian motion. The method is developed and demonstrated on some example data sets.\n\nWe would like to clarify that the **classical** Schrodinger’s bridge (with Brownian reference) is in fact a **special case** of the mvOU-SBP we consider (setting zero drift and isotropic diffusion), rather than the other way around. \nWe would also like to clarify that the case of Gaussian marginals is **only one part of our contribution**, for which we provide thorough analytical results. For the non-analytically tractable (but much more application-relevant) case of **non-Gaussian** initial and final measurements, we introduce an algorithm based on **score and flow matching** to solve the mvOU-SBP.\n\n> I appreciate this is standard in the SBP literature, but I'm always a bit confused as to why it's useful to look at just over 2000 genes, with their expression compressed to 30 PCs. I suspect it's more likely about getting scaling on the data set to work. It would be great to see a larger example that feels a little more plausible.\n\nThank you for raising the point of dimensionality reduction and pre-processing in the single cell data example. Indeed, as you point out, related SBP works like [2] and, more broadly, applications on single-cell RNA sequencing data use these kinds of preprocessing pipelines as the norm. The reviewer points out two choices: \n\n - (i.) **Filtering** to ~2000 highly variable genes\n - (ii.) **PCA** to 30 components.\n\nWe provide additional results scaling up to 100 PCs (see below, Table 1), and explain the necessity of these preprocessing steps in what follows. \n\nIn our case, these preprocessing steps were the same ones used for preparing the dataset in the original publication, see [4, STAR Methods]. In the present work, our main focus is on **theoretical and methodological development** of the SBP with multivariate OU reference, so we did not re-process the data in the initial submission. \n\n**(i. Filtering)**: Many genes are non-informative for the biological process of interest (in this case, the cell cycle) and consist only of technical artifacts or zeros. Filtering by variance is therefore almost **universally used as best practice**, see e.g. “Feature selection” section of [1]: typically \"1,000 and 5,000 HVGs are selected for downstream analysis\". \n\n**(ii. PCA)**: This is again a **standard pre-processing practice** to select a “natural” set of coordinates for downstream analysis and helps with the extremely noisy nature of single cell data. Additionally, the theory of SBP has an underlying Gaussian assumption (e.g. that the reference process is driven by a white noise). PCA performs a kind of whitening transformation to the data, making it more favourable for modelling.\n\n> I didn't quite understand why a larger data set wasn't run at the end. I appreciate the comparisons are with methods that are computationally expensive. But can't you go much larger with this approach? If not why not? Would be good to clarify. I wold have gone higher for significance if we'd been able to see such an example\n\nTo provide some context, for applications of Schrodinger bridges to single cell data\n\n* [2]: 20-30 dim. PCA\n* [3, 5]: 5 dim. PCA\n\nOne limitation which also applies to any method leveraging Gaussian OT [2, 6] is that linear solves involving PSD matrices of dimension $d \\times d$ are required, with a intrinsic computational complexity of $O(d^3)$. We realize that this limitation could be more clearly stated, and we will do so in the revised paper. We point out that there may be ways to alleviate this, e.g. leveraging low-rank approximations of covariance matrices, and this is a potential direction for future work. \n\n**Motivated by reviewer `xnYH`'s feedback, we reprocessed the cell cycle data and scaled up the number of PCs to 50 and 100 dimensions.** \n\n### Table 1 - cell cycle interpolation (equivalent to Figure. 5b)\nMarginal interpolation error between first and last snapshot as a function of the reference velocity scale parameter, $\\gamma$.\n#### (a) 50 PCs\n|t|0.0|10.0|20.0|30.0|40.0|50.0|60.0|70.0|80.0|90.0|100.0|\n|----:|------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|--------:|\n|0|0|0|0|0|0|0|0|0|0|0|0|\n|1|40.18|30.8|22.55|16.34|12.77|11.07|10.26|10.11|10.69|11.93|13.65|\n|2|28.59|18.47|11.53|9.62|12.36|15.32|15|12.52|10.64|11.02|13.96|\n|3|15.31|10.31|7.47|7.06|8.18|8.9|8.5|7.9|8.09|9.42|11.92|\n|4|0|0|0|0|0|0|0|0|0|0|0.01|\n\n#### (b) 100 PCs\nt|0.0|10.0|20.0|30.0|40.0|50.0|60.0|70.0|80.0|90.0|100.0|\n|----:|------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|--------:|\n|0|0|0|0|0|0|0|0|0|0|0|0|\n|1|44.49|36.82|29.97|24.19|19.99|17.48|16.13|15.29|14.78|14.7|15.18|\n|2|32.48|23.75|17.12|13.55|13.79|16.84|19.78|20.21|18.25|15.7|14.22|\n|3|19.33|14.95|12|10.74|11.19|12.68|13.98|14.36|14.03|13.69|13.82|\n|4|0|0|0|0|0|0|0|0|0|0|0|\n\nWe apply also mvOU-OTFM (Alg.1) for different number of PCs, as done in Figure 5(a) of our paper. Since this year there is no possibility to provide images in the rebuttal, we are unable to provide the visualisation of results for $d = 50, 100$. \nHowever, regarding the point of scalability, we report the following runtimes for mvOU-OTFM **on CPU**:\n\n### Runtimes for mvOU-OTFM (Alg. 1) for cell cycle example\n* d = 30, 1m13s\n* d = 50, 1m49s\n* d = 100, 5m20s\n\nWe would like to point out that these runtimes are **significantly faster**, even **without using GPU**, compared to methods relying on numerical integration such as NLSB and SBIRR. \n\nTo further demonstrate scaling our method, we extend the Gaussian results (Table 1 of our paper) from 50 to 100 dimensions. \n\n### Table 2 - Gaussian SB results, $d = 100$\n#### (a) Marginal error (Bures-Wasserstein dist.)\n|   dim | OU-GSB              | BM-GSB           | IPML(←)          | IPML(→)          | NLSB             |\n|------:|:--------------------|:-----------------|:-----------------|:-----------------|:-----------------|\n|   100 | **6.84 $\\pm$ 0.78** | 15.14 $\\pm$ 0.95 | 16.19 $\\pm$ 1.87 | 14.38 $\\pm$ 0.38 | 17.40 $\\pm$ 0.13 |\n#### (b) Force error (L2)\n\n|   dim | OU-GSB               | BM-GSB           | IPML             | NLSB             |\n|------:|:---------------------|:-----------------|:-----------------|:-----------------|\n|   100 | **10.45 $\\pm$ 0.18** | 15.53 $\\pm$ 0.14 | 15.30 $\\pm$ 0.39 | 16.49 $\\pm$ 0.07 |\n\n\n\n[1] Luecken MD, Theis FJ. Current best practices in single‐cell RNA‐seq analysis: a tutorial. Molecular systems biology. 2019 Jun;15(6):e8746.\n\n[2] Bunne C, Hsieh YP, Cuturi M, Krause A. The schrödinger bridge between gaussian measures has a closed form. InInternational Conference on Artificial Intelligence and Statistics 2023 Apr 11 (pp. 5802-5833). PMLR.\n\n[3] Shen Y, Berlinghieri R, Broderick T. Multi-marginal schr\\\" odinger bridges with iterative reference refinement. arXiv preprint arXiv:2408.06277. 2024 Aug 12.\n\n[4] Qiu X, Zhang Y, Martin-Rufino JD, Weng C, Hosseinzadeh S, Yang D, Pogson AN, Hein MY, Min KH, Wang L, Grody EI. Mapping transcriptomic vector fields of single cells. Cell. 2022 Feb 17;185(4):690-711.\n\n[5] Vargas F, Thodoroff P, Lamacraft A, Lawrence N. Solving schrödinger bridges via maximum likelihood. Entropy. 2021 Aug 31;23(9):1134.\n\n[6] Thornton J, Cuturi M. Rethinking initialization of the sinkhorn algorithm. InInternational Conference on Artificial Intelligence and Statistics 2023 Apr 11 (pp. 8682-8698). PMLR."}, {"review_id": "qqYvmLX0q1", "reviewer": "Reviewer_kfao", "summary": "The paper studies the Schrödinger bridge problem with a multivariate Ornstein-Uhlenbeck process as a reference measure. By allowing for a potentially asymmetric drift matrix, this model can be used to study systems exhibiting non-gradient out-of-equilibrium dynamics which are relevant, as stated by the authors, in the contexts of biological and social systems. The paper develops both exact solutions for Gaussian marginals and a simulation-free learning method for general marginals. The approach is based on flow and score matching to avoid expensive numerical simulations. This results in a scalable process for recovering stochastic dynamics from initial and final conditions which is compared to previous methods in different settings.", "questions": "I think the case for the paper would be stronger if the authors included a more thorough discussion on the dynamical meaning of assuming a linear drift and state-independent diffusion coefficients. As it is now it just seems to be for theoretical convenience as it allows for the explicit closed-form expressions that their methods are based on. But there is no methodological discussion on when could these assumptions be reasonable and when not and if, in the future, any of these could be removed or softened.", "limitations": "This is discussed on the previous section.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper is well written and structured and discusses a topic which has generated attention lately. The main theoretical contributions (Theorems 2 and 3) are clean and interesting. These are complemented with comprehensive simulations that show a clear improvement with respect to previous methods. All this makes the paper a good contribution to the venue. However, it appears to me that the main drawback of the paper is that it does not include a strong discussion of the limitations of the strong dynamical assumptions of the model.", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your time taken to evaluate our submission and for the positive assessment of our work. \n\n> it appears to me that the main drawback of the paper is that it does not include a strong discussion of the limitations of the strong dynamical assumptions of the model.\n\n**We agree that a more in-depth discussion of the limitations of the multivariate Ornstein-Uhlenbeck (mvOU) reference process is warranted, and we will add in the revised paper a discussion of this.**\n\nBy considering a multivariate Ornstein-Uhlenbeck reference described by the SDE $dX_t = A (X_t - m) dt + \\sigma dB_t$, our framework can describe a family of SBP processes that are **strictly larger than the Brownian SBP**. In other words, we can recover the standard SBP problem (such as the ones studied in [1, 2]) by taking $A, m = 0$. \n\nWe want to emphasize that the linearity assumption applies **only** to the **reference** dynamics $\\mathbb{Q}$, and not the SB dynamics $\\mathbb{P}$. By modelling the SB flow field using neural networks, the mvOU-OTFM algorithm (Alg. 1) can deal with **general non-linear dynamics** of $\\mathbb{P}$ in the case of general, non-Gaussian marginals. \n\nIndeed, as the reviewer points out, our framework is still not fully general since it does not allow for reference dynamics with non-linear drift, such as those that can be modelled by NLSB [3] and IPFP [4]. Our work **highlights the tradeoff** between (i) analytically tractable but more restrictive models and (ii) more expressive non-parametric models, which are much less tractable to learn, both in terms of computational cost and statistical complexity. \n\n> I think the case for the paper would be stronger if the authors included a more thorough discussion on the dynamical meaning of assuming a linear drift and state-independent diffusion coefficients. As it is now it just seems to be for theoretical convenience as it allows for the explicit closed-form expressions that their methods are based on. But there is no methodological discussion on when could these assumptions be reasonable and when not and if, in the future, any of these could be removed or softened.\n\n**Thank you for pointing this out and we agree that further discussion of this is desirable. We will add discussion of all of these points in the revised paper. Also see above for our response about the assumption of a linear drift.** \n\n* For **application to non-linear systems**, the linear reference dynamics can be **justified by local linearization** of the nonlinear system about a fixed point, assuming that the observed data arise in the vicinity of the fixed point. This is a **canonical analysis technique** for studying nonlinear dynamics in the physics literature [5].  \n \n* While we hinted at this link to linearization about a fixed point (line 263), we will make this connection more explicit: **linearization of a non-linear stochastic dynamical system about a stable point yields a multivariate OU process** where the drift matrix is the Jacobian of the nonlinear system: for a nonlinear SDE $dX_t = f(X_t) dt + \\sigma dB_t$ with a fixed point at $f(0) = 0$, one can consider a linearized dynamics $d\\tilde{X}_t = A \\tilde{X}_t dt + \\sigma dB_t$ where $A = (\\partial_x f)(0)$.\n\n* We would like to emphasize that the benefit of tractability is not only theoretical but also **methodological/computational** – we attribute availability of closed form expressions to the speed and accuracy of our method in our experiments (see Tables 1 and 2 in initial submission), relative to competing methods which rely upon numerical integration for more expressive but expensive solvers.\n\nRegarding scope for assumptions to be relaxed, we note that concurrent work [6] proposes CurlyFM, which trains a neural approximation to the reference process instead of our exact computations. While this may be a direction along which limitations could be relaxed, we leave a detailed comparison to future work.\n\n> In the appendices, many formula fall outside of the margins of the documents. Some reformatting will be needed to void this.\n\nThank you, this has been rectified.\n\n[1] Tong A, Malkin N, Fatras K, Atanackovic L, Zhang Y, Huguet G, Wolf G, Bengio Y. Simulation-free schr\\\" odinger bridges via score and flow matching. arXiv preprint arXiv:2307.03672. 2023 Jul 7.\n\n[2] Bunne C, Hsieh YP, Cuturi M, Krause A. The schrödinger bridge between gaussian measures has a closed form. InInternational Conference on Artificial Intelligence and Statistics 2023 Apr 11 (pp. 5802-5833). PMLR.\n\n[3] Koshizuka T, Sato I. Neural lagrangian schr\\\" odinger bridge: Diffusion modeling for population dynamics. arXiv preprint arXiv:2204.04853. 2022 Apr 11.\n\n[4] Vargas F, Thodoroff P, Lamacraft A, Lawrence N. Solving schrödinger bridges via maximum likelihood. Entropy. 2021 Aug 31;23(9):1134.\n\n[5] Godrèche C, Luck JM. Characterising the nonequilibrium stationary states of Ornstein–Uhlenbeck processes. Journal of Physics A: Mathematical and Theoretical. 2018 Dec 18;52(3):035002.\n\n[6] Petrović K, Atanackovic L, Kapusniak K, Bronstein MM, Bose J, Tong A. Curly flow matching for learning non-gradient field dynamics. InICLR 2025 Workshop on Machine Learning for Genomics Explorations 2025."}, {"review_id": "1qR4OX1m8T", "reviewer": "Reviewer_xh1L", "summary": "This paper studies Schrodinger bridge problem with the reference processes being a multivariate OU process with generic drift. The authors derived a explicit form of the SDE which the bridges should satisfy, as well as a score function and probability flow. With Gaussian marginals, the process can be solved explicitly. In the non-gaussian case, the authors proposed mv-OU-OTFM and show it solves the dynamic. Numerical experiments are done to confirm their claims.", "questions": "See above.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "1. The authors give some introduction in Schrodinger bridges, flow matching, probability flows, etc., which might help readers from a broader background.\n\n2. Although Theorem 1 is adapted from an existing work [1], stating it in this context is still helpful. Similar situation for theorem 2. But I think the authors can be more clear about these results were proved somewhere else, for example, by adding the citation following the name of the theorems.\n\n3. Theorem 3 gives an explicit expression of the Gaussian bridge. The proof is build on the one in [2]. While it generalizes results in [2], one could argue the proof technique is not novel from a mathematic perspective. Note that S should be in upper case in the reference.\n\n4. The authors proposed mvOU-OTFM as the score and flow matching algorithm for OU Schrodinger bridges. The proof idea in Theorem 4 is adapted from [3]. \n\n5. Authors done thorough experiments showing mvOU-OTFM outperforms IPML, NLSB, etc. Impressively, with 4 iteration, mvOU-OTFM achieves better results than SBIRR considering it runs much faster. In page 8, I think by algorithm 2, the authors mean algorithm 1? Table 2 is presented but never referred to. \n\n\n[1] Chen et al. Stochastic bridges of linear systems. 2015.\n\n[2] Bunne et al. The Schrödinger bridge between gaussian measures has a closed form. 2023\n\n[3] Tong et al. Simulation-free Schrodinger bridges via score and flow matching. 2023.", "was_revised": true, "final_justification": "My main concern was the mathematical contribution can be considered weak, but the authors explained to me how their work is different from the existing ones, and I am convinced that this paper has enough theoretical justification with good experimental results. Hence I am updating my rating.", "rebuttal": "Thank you for the time taken to review our work and for your suggestions.\n\n> Although Theorem 1 is adapted from an existing work [1], stating it in this context is still helpful. Similar situation for theorem 2. But I think the authors can be more clear about these results were proved somewhere else, for example, by adding the citation following the name of the theorems.\n\nThank you for these comments and we agree that the motivation for stating these results in this way is to **ensure cohesion and flow** in the text.\nWe agree that the results of Theorem 1 can be read implicitly in [1] and **we will add “adapted from Chen et al.” in the theorem title** in the revised text. \nOn the other hand, the results of Theorem 2 as stated are **not readily available** from any papers that we are aware of. As such, we do not put a citation in the theorem statement. We will add some discussion on the closest related results in the literature, and point out that the result is built upon Theorem 1.\n\n> Theorem 3 gives an explicit expression of the Gaussian bridge. The proof is build on the one in [2]. While it generalizes results in [2], one could argue the proof technique is not novel from a mathematic perspective. \n\nWe agree that parts of the proof technique for Thm. 3 overlap with the details of [2]. Indeed, both the approaches of [2] and ours proceed via the **well-known and classical characterisation of the SBP as a reciprocal process** (see e.g. Eq. 2 of our submission, also [4, Eq. 1.1]) and then utilising the characterisation of $\\mathbb{Q}$-bridges. **In our view, this is the “canonical” path towards building the SBP where it is analytically tractable.**\nWe realise that the proof approach is a helpful discussion and context to provide in the main text, and we will add a few sentences about this and the relation to the approach of [2] in the revised paper. \n\nThe step of the proof that draws from [2] is the derivation of the drift of the SBP SDE via **calculating the SDE’s generator**, and we acknowledge this in the proof (see Supplementary Material, _\"Dynamic $\\mathbb{Q}$-GSB : SDE representation\"_). \nWe argue that our derivation **streamlines** the argument compared to [2], making the derivations easier to follow while being more general. In particular, [2] utilises \"heavy machinery\" such as the \"central identity of QFT\" (see [2], Preliminaries for the Proof of Theorem 3) to compute the generator. In our analysis, we find that this is in fact unnecessary and standard arguments suffice. \n\n> Note that S should be in upper case in the reference.\n\nThank you, indeed the name should be capitalised. We will rectify this in the revised paper.\n\n> Authors done thorough experiments showing mvOU-OTFM outperforms IPML, NLSB, etc.\n\n**We point out that we have now added additional experiments, scaling up to $d = 100$ for the Gaussian example and $d = 50, 100$ for the cell cycle example (see Tables 1, 2 of our response to `xnYH`)**\n\n> In page 8, I think by algorithm 2, the authors mean algorithm 1?\n\nAlgorithm 2 refers to the \"iterated reference fitting\" approach that uses Alg. 1 as an inner subroutine. Due to issues with space, the pseudocode for Alg. 2 is contained in the appendix (see Supplementary Material). Thank you for bringing it to our attention that this may be confusing -- in the revised paper, we will either include it in the main text (using the extra 10th page for camera ready), or put an explicit reference to appendix. \n\n>   Table 2 is presented but never referred to.\n\nThank you for pointing this out, this has been rectified. \n\n[1] Chen et al. Stochastic bridges of linear systems. 2015.\n\n[2] Bunne et al. The Schrödinger bridge between gaussian measures has a closed form. 2023\n\n[3] Tong et al. Simulation-free Schrodinger bridges via score and flow matching. 2023.\n\n[4] Léonard C. A survey of the schr\\\" odinger problem and some of its connections with optimal transport. arXiv preprint arXiv:1308.0215. 2013 Aug 1."}, {"review_id": "Wd7d3Qf02F", "reviewer": "Reviewer_uH9r", "summary": "The study frames the Schrödinger bridge (SB) problem around a multivariate Ornstein–Uhlenbeck reference process, keeping the math tractable while covering non-equilibrium systems. It provides a closed-form solution for Gaussian marginals (mvOU-GSB) and, for arbitrary marginals, introduces MVOU-OTFM, a simulation-free method that learns flow and score via conditional flow matching rooted in exact mvOU-bridge formulas. Across several benchmarks, the proposed approach trains quickly and achieves the lowest error rates among all baselines.", "questions": "* The reference process ought to be treated as a prior that is chosen before any data are seen. In SBIRR the prior dynamics remain fixed and are not directly inferred from the observations. Such an experiment in single cell-data, If we instead fit the OU parameters directly to the data, how is that fundamentally different from simply augmenting SBP-Dyn with additional parameters?\n   \n* The experiments are limited to low-dimensional settings with small sample sizes, such as Repressilator or cell-dynamics data. Can the proposed algorithm still produce an optimal transport plan for larger-scale datasets?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n* The paper delivers the first closed-form Schrödinger bridge solution for Gaussian marginals under a fully multivariate, possibly asymmetric OU reference, broadening the scope of Gaussian SB theory.\n\n* Leveraging these formulas, MVOU-OTFM eliminates SDE simulations by directly regressing flow and score fields, yielding order-of-magnitude speed-up.\n\n* From synthetic mixtures to 30-dimensional single-cell RNA-seq, MVOU-OTFM consistently achieves the lowest error.\n\nWeaknesses:\n\n* The paper’s claim that extending a scalar OU SDE to a multivariate one is a major contribution is not persuasive. Multivariate OU processes and their closed-form solutions are already well known. A stronger advance would be to use a full non-diagonal diffusion matrix and to show a simulation-free bridge algorithm. Instead, the authors keep a diagonal diffusion matrix in their experiments, so the theorem feels too limited to count as a significant contribution on its own.\n\n* The central aim of any Schrödinger Bridge algorithm is to learn dynamics that generate the optimal coupling, since that coupling is notoriously hard to identify for real data. Section 3.3 claims that a single Sinkhorn solve achieves this goal, yet the claim is fragile. Sinkhorn tackles an entropically regularised transport problem, its cost and memory scale poorly with dimension and sample size, and the regularisation introduces a persistent bias unless the regularisation strength tends to zero. Under these conditions the resulting dynamics cannot be considered optimal SB solution. Established SB methods improve accuracy through iterative proportional fitting (IPF) [1,2] or iterative Markovian fitting (IMF) [3,4]. Please add either empirical results or a theoretical argument that compares the proposed single Sinkhorn step to these iterative procedures. Otherwise the use of “SB: exactly solvable” in the title feels misleading.\n\n----\n    [1] Diffusion schrödinger bridge with applications to score-based generative modeling.\n    [2] Likelihood Training of Schrödinger Bridge using FBSDEs Theory.\n    [3] Diffusion Scrödinger Bridge Matching.\n    [4] Diffusion Bridge Mixture Transports, Schrödinger Bridge Problems and Generative Modeling", "was_revised": true, "final_justification": "Although the mvOU process is well known and its derivations already exist, I acknowledge that a closed-form matching objective for the dynamic-SB problem with an mvOU prior has not been thoroughly treated in the literature. This work applies mvOU to the SB problem and demonstrates its effectiveness on multi-marginal data. The main limitation is the reliance on a pre-computed OT coupling, which can bias the coupling in high-dimensional settings. The manuscript should explicitly discuss this approximation error.", "rebuttal": "Thank you for your time to read our paper and we appreciate your detailed feedback.\n\n> ... extending a scalar OU SDE to a multivariate one ... is not persuasive. Multivariate OU processes ... are already well known.\n\nWe **agree** with the reviewer that multivariate OU (mvOU) processes are well-known – their analytical tractability is **exactly what we exploit** to develop our theoretical and computational results.\n\nAs the reviewer acknowledges, our characterisation of the the Schrodinger bridge with mvOU reference are **novel** and we consistently achieves the **lowest error** with much less computation compared to approaches that are more general at the expense of relying on simulation.\nWe attribute these improvements to the **use of closed form expressions for the flow and score**. \n\nWe point out that the gap between scalar OU SDEs (as in generative modelling literature) and their generic multivariate counterparts is **non-trivial** – in fact, a scalar OU process is _always_ gradient-driven. Thus, while sufficient for generative modelling tasks, scalar OU dynamics are **not appropriate for modelling non-equilibrium systems in physical systems**. \n\nInstead, mvOU processes allow for a **generic drift matrix**. The motivation for mvOU dynamics arises naturally -- for a nonlinear system (such as the oscillator in our Figure 4), a **canonical analysis technique** [1, 2] in the physics literature is to study its linearization -- this is a **linear SDE** where the **drift matrix is exactly the Jacobian of $f(x)$**. \n\nConcretely, an SDE $dX_t = f(X_t) dt + \\sigma dB_t$ with a fixed point at $f(0) = 0$, can be linearized: $d\\tilde{X}_t = A \\tilde{X}_t dt + \\sigma dB_t$ and $A = (\\partial_x f)(0)$. \n\nIn the case of gradient dynamics $f(x) = \\nabla_x U(x)$, note that $A = (\\partial_{xx} U)(0)$ is **symmetric**. The relevance of considering an **asymmetric drift matrix** is to model non-gradient systems. \n\n> A stronger advance would be to use a full non-diagonal diffusion matrix and to show ... algorithm\n\nThank you for the suggestion. We acknowledge that our examples use diagonal diffusion. As we explain below, we provide theoretical grounds that considering diagonal diffusion and generic drift is sufficient. \n\n- **(i).** **All** of our theoretical results (Thms. 1-4) and our algorithms (Algs. 1-2) are already for a **generic diffusion matrix** $\\bm{\\sigma}$ (see e.g. line 164, 213)\n- **(ii).** **Any** Linear SDE with a non-diagonal diffusion matrix is **equivalent** to one with a diagonal diffusion matrix under a change of basis.\n\nFor **(i)**, while we typeset vector quantities in bold, we will make this convention explicit.\n\nFor **(ii)**, note that **any diffusion matrix** $D = \\tfrac{1}{2} \\sigma \\sigma^\\top$ can be diagonalized by a change of basis, since it is PSD. \n\nTherefore we can **generically work with the diagonal diffusion setting**, since given any non-diagonal diffusion we can choose an orthogonal basis of $\\mathbb{R}^d$ under which it becomes diagonal. **On the other hand**, asymmetry in the drift component **cannot** be eliminated by the same trick. \n\nWe formalise this as follows:\n\n**[Lemma.] Up to a change of basis, any linear SDE with generic drift and diffusion matrix is equal in its law to a linear SDE with diagonal diffusion.**\n\nLet $d X_t = A X_t dt + \\sigma d B_t$ be a linear SDE with generic drift $A \\in \\mathbb{R}^{d \\times d}$ and diffusion $\\sigma \\in \\mathbb{R}^{d \\times d}$. There exists an orthogonal basis $U \\in \\mathbb{R}^{d \\times d}$ such that $Y_t = U^\\top X_t$ obeys a linear SDE with drift matrix $U^\\top A U$ and a diagonal diffusion matrix $\\Lambda = \\mathrm{diag}(\\lambda)$, $\\lambda = (\\lambda_1, \\ldots, \\lambda_d) \\ge 0$\n\n_Proof_ Let $D = \\frac{1}{2} \\sigma \\sigma^\\top$, this is PSD and can be written $D = \\frac{1}{2} U \\Lambda^2 U^\\top$ where $U$ is orthogonal and $\\Lambda = \\mathrm{diag}(\\lambda) \\succeq 0$. Consider the SDEs\n(a.) $$ d X_t = A X_t dt + \\sigma d B_t$$ and \n(b.) $$ d X_t = A X_t dt + U \\Lambda U^\\top d B_t.$$\n(a.) has generator $(\\mathcal{L} f)(x) = \\langle \\partial_x f(x), Ax \\rangle + \\frac{1}{2} tr(\\sigma^\\top \\partial_{xx} f(x) \\sigma)$.\n\nComputing the generator for (b): $(\\tilde{\\mathcal{L}} f)(x) = \\langle \\partial_x f(x), Ax \\rangle + \\frac{1}{2} tr( U \\Lambda U^\\top \\partial_{xx} f(x) U \\Lambda U^\\top)$.\nNote that $tr( U \\Lambda U^\\top \\partial_{xx} f U \\Lambda U^\\top) = tr(\\partial_{xx} f U \\Lambda^2 U^\\top) = tr(\\sigma^\\top \\partial_{xx} f \\sigma)$. \nHence both SDEs share the same generator, and so are equal in law. \n\nSince orthogonal transformations of Brownian increments are also Brownian increments, we have that $U^\\top dB_t = d\\tilde{B}_t$. Substituting into (b.), we have $d X_t = A X_t dt + U \\Lambda d \\tilde{B}_t$. Then,\n$$ d(U^\\top X_t) = U^\\top A X_t dt + \\Lambda d \\tilde{B}_t \\Longrightarrow d Y_t = U^\\top A U Y_t dt + \\Lambda d \\tilde{B}_t.$$\n$□$\n\n> Section 3.3 claims that a single Sinkhorn solve achieves this goal, yet the claim is fragile...\n\nWe aim to solve the SBP with **discrete** initial and final data, in the form of **empirical distributions**. This is different from the setting of estimating the **continuum SBP from samples** where the initial and final _densities_ must be estimated as well as the law of the bridge (e.g. as in [3]). The **SBP between empirical distributions**, the problem we consider, is well-defined as-is (see e.g. Section 6 of [4]), and that is the viewpoint we take.\n\nPut differently, for source and target $(\\mu, \\nu)$ supported on $\\{ x_1, \\ldots, x_M \\}$ and $\\{ y_1, \\ldots, y_N \\}$ respectively, the $\\mathbb{Q}$-SB takes the form $\\mathbb{P} = \\sum_{i, j} \\pi_{ij} \\mathbb{Q}^{(x_i, y_j)}$, where $\\{ \\pi : \\sum_{ij} \\pi_{ij} = 1 \\}$ is the solution to an **entropic OT problem**.\n\nThat this is the case can be read from Equations 1.1-1.2 of [4]; alternatively Proposition D.1 of [5] or Prop 2.1 of [6]. This is exactly what we discuss in Section 2.1 of our paper, and the SB between empirical distributions is **precisely** what we aim to learn.\n\nRegarding regularisation, we point out that it is a standard fact in SBP literature that the entropic term in the \"static\" EOT problem (see Eq.(SBP-static) in our paper, also Prop. 2.1 of [6]) is **tied to the noise level of the system**, and therefore **must not be taken to zero**. This is in contrast to Monge map estimation, in which any $\\varepsilon > 0$ introduces a bias.\n\nWe acknowledge that these points deserve an extended discussion, and we will add this to the revised paper.\n\n> Established SB methods improve accuracy through (IPF) [1,2] or iterative Markovian fitting (IMF) [3,4]\n \nWe gently point out that **we currently compare to IPFP**: [7, 8] are based on IPF with various drift families. We show our approach is more accurate and faster.\n\nIn sum, the use of a single Sinkhorn step for learning the SB between two empirical distributions is justified both from a theoretical point of view (see our answer to the question above and [4, 5, 6]), as well as empirical observations. For the Gaussian mixture benchmarks (Table 1 of our paper), \n\n**Note that we have now increased the number of dimensions from 50 to 100 (see Table 2 in our response to `xnYH`)**\n\n> Otherwise the use of “SB: exactly solvable” in the title feels misleading.\n\nTo clarify, the term “exactly solvable” in our title refers to our **analytical results for the Gaussian case**, owing to the existence of matrix-analytic expressions for quantities involving Gaussian marginals and mvOU dynamics.\n\n**In this setting, our results provide an exact solution to the SBP**. We see that ‘exactly solvable’ may confuse readers, and we will make the abstract and introduction more explicit to avoid risking ambiguity. \n\n> The reference process ought to be treated as a prior that is chosen before any data are seen. In SBIRR the prior dynamics remain fixed and are not directly inferred from the observations. Such an experiment in single cell-data, If we instead fit the OU parameters directly to the data, how is that fundamentally different from simply augmenting SBP-Dyn with additional parameters?\n\nTo clarify, in SBIRR prior dynamics are **not fixed but updated at each step** [8, Alg. 1], thus using information across observations to update the prior\n\nIn our cell cycle example, OU params are fit to RNA velocity data [9] -- we construct the prior using information that is **external** to that used for fitting SB. \n\n> The experiments are limited to low-dimensional settings with small sample sizes... Can the proposed algorithm still produce an optimal transport plan for larger-scale datasets?\n\n**Scaling in $N$**: [Alg. 1] samples **batches** at each iteration so it can be scaled to **arbitrarily large datasets**. Theoretically one should solve EOT over the full dataset, but minibatch OT is widely used to scale to large, even streaming data [6]\n\n**Scaling in $d$**: **In Table 1,2 of our response to `xnYH`**, we report new results for single cell $d = 50, 100$ and the synthetic Gaussian for $d = 100$. Considering most single cell studies use $d < 100$ PCs, our method can scale to the $d$ used in practice\n\n[1] Godrèche C, Luck JM. Characterising the nonequilibrium stationary states of Ornstein–Uhlenbeck processes\n\n[2] Lax M. Fluctuations from the nonequilibrium steady state\n\n[3] Pavon M, Trigila G, Tabak EG. The data‐driven schrödinger bridge\n\n[4] Léonard C. A survey of the schrodinger problem and some of its connections with optimal transport.\n\n[5] Lavenant H, Zhang S, Kim YH, Schiebinger G. Towards a mathematical theory of trajectory inference\n\n[6] Tong A, Malkin N, Fatras K, et al Simulation-free schrodinger bridges via score and flow matching\n\n[7] Vargas F, Thodoroff P, Lamacraft A, Lawrence N. Solving schrödinger bridges via maximum likelihood\n\n[8] Shen Y, Berlinghieri R, Broderick T Multi-marginal schrodinger bridges with iterative reference refinement\n\n[9] Qiu X, Zhang Y, et al Mapping transcriptomic vector fields of single cells"}]

Author rebuttal:

Meta-review: The paper generalizes the classical Schr\"odinger's bridge to a case where the reference path measure is induced by an OU process with asymmetric linear coefficient matrix. Analytical solution is worked out for a special case where both marginals are Gaussian, which is not surprising, but a simulation-free algorithm is also proposed for general marginals. The method is then demonstrated on relatively low-dimensional data sets. Both reviewers and I found the idea interesting, and we appreciated that details are worked out. For the benefit of the authors, however, I'd like to point out that asymmetric A means the process is irreversible, but not non-equilibrium (OU is still geometric ergodic in your context where A is negative definite, which means it converges to an equilibrium distribution). I'd also avoid saying asymmetric A means non-conservative forces (as in the abstract), because dissipative A is never conservative,while on the other hand A=[0 1; -1 0] for example gives perfectly conservative dynamics. However, the correction/removal of these statements does not really affect the machine learning component of this paper, and based on reviewers' consensus, I'm happy to recommend acceptance. The authors are encouraged to take all the discussions into consideration when preparing for a revised/final version.
Decision: Accept (poster)