Momentum Multi-Marginal Schrödinger Bridge Matching
Authors: Panagiotis Theodoropoulos
Augustinos D Saravanos
Evangelos Theodorou
Guan-Horng Liu
Venue: NeurIPS 2025 poster ()
Keywords: Diffusion models
Schrödinger bridge
Distribution matching
Trajectory Inference
Primary area: deep_learning
TL;DR: We introduce Momentum Multi-Marginal Schrödinger Bridge Matching (3MSBM), a novel matching framework that learns smooth trajectories for stochastic systems that satisfy multiple positional constraints.

Abstract: Understanding complex systems by inferring trajectories from sparse sample snapshots is a fundamental challenge in a wide range of domains, e.g., single-cell biology, meteorology, and economics. Despite advancements in Bridge and Flow matching frameworks, current methodologies rely on pairwise interpolation between adjacent snapshots. This hinders their ability to capture long-range temporal dependencies and potentially affects the coherence of the inferred trajectories. To address these issues, we introduce Momentum Multi-Marginal Schrödinger Bridge Matching (3MSBM), a novel matching framework that learns smooth measure-valued splines for stochastic systems that satisfy multiple positional constraints. This is achieved by lifting the dynamics to phase space and generalizing stochastic bridges to be conditioned on several points, forming a multi-marginal conditional stochastic optimal control problem. The underlying dynamics are then learned by minimizing a variational objective, having fixed the path induced by the multi-marginal conditional bridge. As a matching approach, 3MSBM learns transport maps that preserve intermediate marginals throughout training, significantly improving convergence and scalability. Extensive experimentation in a series of real-world applications validates the superior performance of 3MSBM compared to existing methods in capturing complex dynamics with temporal dependencies, opening new avenues for training matching frameworks in multi-marginal settings.

Momentum Multi-Marginal Schrödinger Bridge Matching
Panagiotis Theodoropoulos! Augustinos D. Saravanos', Evangelos A. Theodorou' '*, Guan-Horng Liu
Georgia Institute of Technology, 2FAIR at Meta, "Equal advising
Abstract
Understanding complex systems by inferring trajectories from sparse sample snap- shots is a fundamental challenge wide range of domains, e.g., single-cell biology, meteorology, and economics. Despite advancements in Bridge and Flow matching frameworks, current methodologies rely on pairwise interpolation be- tween adjacent snapshots. This hinders their ability to capture long-range temporal dependencies and potentially affects coherence the inferred trajectories. To address these issues, introduce Momentum Multi-Marginal Schrüdinger Bridge Matching (3MSBM), a novel matching framework that learns smooth measure-valle valued splines for stochastic systems that satisfy multiple positional con- straints. This is achieved by lifting the dynamics to phase space and generalizing stochastic bridges be conditioned conditional stochastic optimal control problem. The underlying dynamics are then learned by minimizing variational objective, having fixed the path induced by the multi-marginal conditional bridge. As a matching approach, 3MSBM learns trans- port maps that preserve intermediate marginals throughout training, significantly improving convergence and scalability. Extensive experimentation in a series real-world applications validates the superior performance of 3MSBM compared to existing methods in capturing complex dynamics with temporal dependencies, opening new avenues for training matching frameworks in multi-marginall settings.
several points, forming multi-marginal
Introduction
Transporting samples between probability distributions is fundamental problem in machine learning. Diffusion Models (DMs:[Ho al., 2020, Song et al., 2020]) constitute a prominent technique in generative modeling, which employ stochastic mappings through Stochastic Differential Equations (SDEs) to transport data samples to a tractable prior distribution, and then learn to reverse this process [Anderson, 1982, Vincent, 2011]. However, diffusion models present several limitations, e.g., the lack of optimality guarantees concerning the kinetic energy for their generated trajectories [Shi et al., 2023]. To address these shortcomings, principled approaches that stem from Optimal Transport [Villani et al., 2009] have emerged that aim to minimize the transportation energy of mapping samples between two marginals, To and T. In this vein, the Schrödinger Bridge (SB; [Schrödinger, 1931))-eqiivalen to Entropic Optimal Transport (EOT:(Cuturi, 2013, Peyré et al., 2019]>-has been of the most prominent approaches used in generative modeling [Vargas et al., 2021]. This popularity has been enabled by recent remarkable advancements of matching methods [Lipman al., 2022, Liu al., 2022a]. Crucially, these matching-bas frameworks circumvent the need to cache the full trajectories of forward and backward SDEs, mitigate the time-diccretization and "forgetting" issues encountered in earlier SB techniques, and maintain a feasible transport map throughout training. This renders them highly scalable and stable methods for training the SB [Shi al., 2023, Gushchin al., 2024a, Peluchetti, 2023, Liu et al., 2024, Rapakoulias al., 2024].
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Table 1: Comparison between our 3MSBM and ofate-of-the-an art multi-marriina algorithms in 1)
simulation- training, 2) smooth and coherent trajectories and 3) globally optimal coupling.
Simulation-Fire Training
Smooth Trajectories
Global Coupling
DMSB [Chen et al., 2023a]
MMFM [Rohbeck al., 2024]
SBIRR [Shen et al., 2024]
X
MMSFM [Lee al., 2025]
X
Smooth SB [Hong et al., 2025]
X
/
3MSBM (Ours)
/
However, many complex real-world scenarios provide sep- arate measurements at coarse time intervals [Chen al., 2023a]. In this case, solving several distinct SB problems between adjacent marginals and connecting the ensuing bridges leads to suboptimal trajectories that fail to account for temporal dependencies or model the dynamics without discontinuities [Lavenant et al., 2024]. To address this is- sue, a generalization of the SB problem is considered: the multi-marginal Schrödinger Bridge (mmSB) [Chen et al., 2019], in which the dynamics are augmented to the phase space. Incorporating the velocity and coupling it with the position leads to smooth trajectories traversing multiple time-indexed marginals in position space (Dockhorn et al., 2021, Chen et al., 2023a], illustrated in Figure 1. This approach enables to capture the underlying dynamics and better leverage the information-rich data of complex systems such as cell dynamics (Yeo al., 2021, Zhang et al., 2024a], meteorological evolution [Franzke et al., 2015], and economics [Kazakevičius et al., 2021]. This wide applicability of multi-marginal systems in real-world applications has spurred significant interest in developing algorithms to address these problems. However, existing methodologies either sacrifice optimality for scalability or vice versa.
3MSBM (Ours)
Vanilla SB
13
02
4s
Trajoctorics
Figure 1: Trajectory comparison be- tween by vanilla SB, and our 3MSBM.
On one end of the spectrum, many recent methods optimize locally between adjacent marginals, often exhibiting substantial scalability, yet failing to recover the global coupling. For instance, Shen al. [2024] propose simulation afreeeeeio--eeeeeeee iterative scheme that solves mmSB by performing trajectory inference with pairwise bridge optimization. Nevertheless, this pairwise formulation cannot enforce global trajectory consistency and thus depends informative priors, limiting practicality in real- world settings. Similarly, Multi-Marginal Flow Matching [MMFM; Rohbeck et al. [2024]] uses deterministic cubic-spline interpolation with precomputed piecewise couplings, impairing temporal coherence and dynamical fidelity. Furthermore, the stochastic counterpart of MMFM Martinaarginal Stochastic Flow Matching [(MMSMM) Lee al. [2025]] optimizes measured value splines over overlapping triplets of marginals, improving on MMFM's pairwise scheme; however, the lack of global coupling and the first-order dynamies still limit smooth, temporally coherent trajectories. On the other hand, methods that solve the mmSB without local approximations can recover the global coupling but generally scale poorly. Deep Momentum Murti-Marginal Schrödinger Bridge [DMSB; Chen et al. [2023a]] solves for the global mmSB coupling in phase space using Bregman iterations. However, it suffers from scalability limitations due to the need cache full SDE trajectories, leading to computational bottlenecks, error accumulation, and potential instability. Lastly, modeling the reference dynamics as smooth Gaussian paths achieved temporally coherent and smooth trajectories (Hong et al., 2025], though the belief propagation prohibits scaling in high dimensions. Table 1 summarizes the key attributes of these approaches alongside our proposed methodology.
In this work, we introduce Momentum Multi-Marginal Schrödinger Bridge Matching (3MSBM). novel scalable matching algorithm that solves the mmSB in the phase space. We lift the dynamics to phase space and minimize path acceleration, resulting in smooth low- low-curvature trajectories. This preserves trends often lost in first-order models crucial cruciall under sparse or irregular observations and captures non-linear transitions and inflection points for more realistic interpolation. We begin by deriving momentum Brownian bridges, yielding closed-fore form expressions for conditional bridges that traverse any arbitrary number of marginals. This eliminates the need for costly numerical integrations, improving efficiency and avoiding error accumulation. The resulting conditional path- path-optimal under stochastic optimal control and satisfying all positional marginal constraints is then held fixed while

---

we match parameterized drift to the prescribed trajectory, enabling simulation-free drift learning. Iterating these steps converges the globally optimal mmSB coupling. Algorithmically, our 3MSBM shares similar attributes with other matching frameworks [L.iu et al., 2024], as the model-induced marginals remain close to the ground truth throughout the optimization, unlike prior mmSB solvers that align with the targets only at convergence. Empirical results verify the efficiency and scalability of our framework, handling high-dimensional problems and outperforming state-of-the-art methods tailored to tackle multi-marg marginal settings, Our contributions are summarized as follows:
We propose 3MSBM, novel matching algorithm for learning smooth interpolation, while preserving multiple marginal constraints.
We present theoretical analysis extending the concept of Brownian Bridges to second-order systems with the capacity to handle arbitrarily many marginals.
Unlike prior mmSB methods, e-g. [Chen et al., 2023a], our method enjoys provable stable convergence and admits map that satisfies the marginal constraints throughout training.
Extensive experimentation demonstrates the enhanced scalability of 3MSBM, with respect to the dimensionality of the input data and the number of marginals.
2 Preliminaries
2.1 Schrödinger Bridge
The Schrödinger Bridge can be obtained through the following Stochastic Optimal Control (SOC) formulation, trying to find the unique non-linear stochastic process 2t E Rdd between marginals To. TT that minimizes the kinetic energy [Chen et al., 2016]
min L' Ep. [lluu.ll' l|] s.t. dz, = u,dt + odW., Mc-Pi
F0 TO
IT
TT
resulting in the stochastic equivalent [Gentil et al., 2017] of the fluid dynamic formulation in OT [Be- namou and Brenier, 2000]. Specifically, the optimal drift of Eq. (1) generates the optimal probability path Pt of the dynamic Schrödinger Bridge (dSB) between the marginals To. and TT. More recently, the advancement of matching algorithms [Gushchin et al., 2024b, De Bortoli et al., 2023] led to the development of highly efficient and scalable SB Matching algorithms. Representing the marginal bridges, Ps = Prio, d (Zo, Zr)
probability path pt as a mixture of endpoint-condition [Léonard, 2013], motivates a two-step alternating training scheme (Theodoropoulos al., 2024]. The first step entails fixing the coupling Tu. by drawing pairs samples (zo, ØT) and optimizing intermediate bridges between the drawn pairs. Subsequently, the parameterized drift uf is matched given the prescribed marginal path from the previous step, progressively refining the coupling induced by ul At convergence, these steps aim to construct a stochastic process whose coupling matches the static solution of the SB, i.e. Torr and optimally interpolates the coupling (æo, 21) Tō.r-
2.2 Momentum Multi-MMagiina Schrödinger Bridge
The Momentum Multi-Marginal Schrödinger Bridge (mmSB) extends the objective in Eq. (1) to traverse multiple marginals constraints. Additionally, the dynamics are lifted into second-order, incorporating the velocity, denoted with U, E Rd, along with the position Iz. Consequently, the I,, J. they de-
marginal distributions are also augmented nn :- n. (I. v). for n pend on both the position It, and the velocity et. We define the joint marginal probability path Pe(z,, v.) for E [), and the position and velocity marginals with (=()) = (c,, d)dv, and (.(w) (2; r)dr respectively. Application of the Girsanov theorem in the phase space yields the corresponding multi-marginal phase space SOC formulation [Chen et al., 2019]
0 at
/ Epp [Ilael1]jat
"
arg min
dm; Am,d + wedt + gdWe,
In gn =
Tn( v)duns -0 10.,,...,A NJ.
10
0 0
where m. = [xe, V E 122d
and g
=
3

---

3 Momentum Multi-Marginal Schrödinger Bridge Matching
We propose Momentum Multi-Marginal Schrödinger Bridge Matching (3MSBM), a novel matching framework, which incorporates the velocity into the dynamics to learn smooth measure- valued splines for stochastic systems satisfying positional marginal constraints over time. Following recent advances in matching frameworks, our algorithm solves Eq. (2), separating the problem into two components: 1) the optimization of the intermediate path, conditioned on the multi-marginal coupling, and 2) the optimization of the parameterized drift, refining the coupling.
3.1 Intermediate Path Optimization
Our analysis begins by formulating a multi-marginal conditional path that satisfies numerous con- straints at sparse intervals. Specifically, we first derive the optimal control expression for phase-space conditional bridge, conditioned at [m., for n E {0, 1,..., N}. Based on this expression for the optimal control, we then obtain a recursive formula for the conditional acceleration that interpo- lates through a set of prescribed positions on ~ qn (I We denote this set of fixed points as } = {To, ,1,..., .N) }, and the coupling q((z))).
Theorem 3.1 (SOC representation of Multi-Marginal Momentum Brownian Bridge (3MBB)). Con- sider the following momentum system interpolating among multiple marginals
ae min L. 2 Ilaall Σ
1
(m. (..,))(M,,, m.n)
(3)
s.t dmp = Amydt + uedt + gdW.
= ma
(4)
We define the value function as V.(ma) := mm7 P. + m. P. and first-order approximations, respectively. This formulation admits the following optimal control expression uz (m²) -gg" (m, + rT). For the multi-margina marginal bridge with {ma) fixed at {tn),
where Ps, are the second-
for n E 10. N), the dynamics of P, and rt obey the following backward ODES
P. APt + P.A -oio7
Pn = -P,7 + +R)- with PNI = =0
(5)
r1 Ar,,
Tn P. (P Rmm,)
where Pat = lim, Ps and rut == lim ru for / {s {s:8 (11, (3)...))..)
Our 3MBB presents natural extension of the well-established concept of the momentum Brownian Bridge [Chen and Georgiou, 2015]. For the derivation of the multi-marginal bridge, we apply the dynamic programming principle, recursively optimizing acceleration each segment while accounting for subsequent segments via the intermediate constraints, as illustrated in Figure 2. The terms Ps+. and ra+ capture the influence of the subsequent segment through the intermediate constraints PR, Tn. which serve as terminal conditions for the corresponding ODEs of the next segment. Importantly, from the terminal conditions in Eq. (5) it is implied that r: and hence the optimal control u/- woold depend on all subsequent pinned points after t {m,, t, 2 t}. as shown explicitly, in the acceleattin formulation derived in the next proposition.
0] o. At the limit when c - 0, the solution of 3MBB (Th. 3.1) admits
Proposition 3.2. Let R =
closed-form expression on every segment; in particular, for E [I I141:
N Σ j=n+1
atmelHEnt1 bot1 2 ()) = - +. CT(()
(6)
where {ān+1 : tut1 2 t} signifies the bridge is conditioned on the set of the ensuing points, X are static coefficients and CT CT CT (t) are time-varying coefficients specific to each segment. The proof. the definitions for these functions, and the Xj coefficient values are left for Appendix B.2.
The expression in Eq. (6) provides a recursive formula to compute the optimal conditional bridge for the segment E [tn, tn)1.. Notice that the linear combination Σj=n+1 dyæ captures the dependence of each segment on all next pinned points after t torrtt((eeee {mnt+1 t,+1 2 t).

---

Intermediate Murginal Costs
Pr. r2: Backward ODEs
0. Osts
4 1<
12 22 Sttt
Fa - (Pg R)
R
Ps
A - (P:. R)
R
Ps
(Pp. - R)
R
P
Pg R-
PoPPotrs' (Pgg)
T5
7. RR,,)
2.(22- Tz- m,)
-m,
1Rg
m
m2
m3
Figure 2: Visualization of the Dynamic Principle. Pt, r. Eq. (5) are solved backward, propagating the influence of future pinned points in preceding segments, through the intermediate constraints.
Remark 3.3. Importantly, the expression for the acceleration in Eq. (6) explicitly shows that the optimal bridge does not need to converge predefined velocities En the intermediate marginals.
As 0 in the intermediate state costs, the constraints on the joint variable m, shift to ensure the trajectory reaches the conditioned I time ln without explicitly prescribing any velocities the inter- mediate points, consistent with the principles of Bridge Matching.
2.0
0.5
1.
Time stog
0.5
1.0
Example- Bridge for 3 marginals To obtain the stochastic bridge of Eq. (3) between two marginals, i.e., N = 2, we consider the same value function approximation
-10
0.5
Pinned Points
in Th. 3.1, admitting the same
0.0
optimal control law ug (m+)
gg T P. P--1(m (m+ + rr). For simplicity, let us assume to = 0, 4 = 1, and l2 = 2. From Prop. 3.2, we obtain the following solution for conditional acceleration.
Figure 3: Bridge materializa- tions among pinned points.
a am T2)
(t- (t-2)z (i2 x²)
2-1 U4
3
t E [1,2)
(7)
a* (mm((11,
(t- (3t-777777) (88-7) (rx - 51) + (e-1)(31-)) 77th + (31-7)
30
12( 2)
(52 I1); E [0,
Notice that the segment t E [0. 1), the À coefficients of the linear combination > - Tj in Eq. (6)
found to be: λ 1, 12 B.2 presents more examples of conditional accelerations with more marginals. demonstrating the
1, whereas for t E [1,2] the sole coefficient is 0. Appendix
dependency of the functions (1), C2 (t), CT),CC(())) (t), along with the A; coefficient values on the number of the following marginals. Figure 3 depicts different materializations between 3 pinned points, illustrating the convergence of the bridges the conditioned points.
3.2 Bridge Matching for Momentum Systems
Subsequently, following the optimization of the 3MBB, we match the parameterized acceleration of given the prescribed conditional probability path p(m, H t, 2 (}), induced by the acceleration a" (m {Tn Ln 2 f)). Given that the 3MBB is solved for each set of points (En). we can marginalize to construct the marginal path Pt-
Proposition 3.4. Les us define the marginal path Pt
mixture of bridges pe(m:) = [ Pilla.) (ont,.)m,,le (m, : tn 2 (d))((())))) where Pel(a,) (mal :t, 2 t}) is the conditional proba- bility path associated with the solution of the 3MBB path in Eq 6. The parameterized acceleration that satisfies the FPE prescribed by the pt. is given by
Pt 2/ a.l(2.)) (mellan t, 2
(t, mi)
(8)
This suggests that the minimization of the variational gap to match at given pt given by
min Bg E Epsitens 1//
a((12d1)
Matching the parameterized drift given the prescribed path pt leads to a more refined coupling 2( (z,)), which will be used for the conditional path optimization the next iteration. The linearity of the system implies that we can efficiently sample m, = [x,, vel. Vtt [0, 77, from the conditional
5

---

Algorithm Momentum Multi-Marginal Schrödinger Bridge Matching (3MSBM)
1: Input: Marginals g(z)).)(e))....(((..)) a, K,T.
2: Initialize of, (((xx,)) = ((ro) & ((x1) 8 ... 8 ((xN), and vo ~ ~N(0
3: repeat
4: for to J
5: Calculate de](=,) using Eq. (6) for t from 0 to T
6: UN f- sdeint to. (dl(s,)) K)T)
7: Calculate (e]{2,) using Eq, (6) for from T to
8: Vo f sdeint(Ix. UN Od](=.) J,K,K, K)
9: end for
Update af, from Eq. (9) using (r)
11: Sample new Zo, T1....; 2N from af
12:
until converges
probability path (s (Wr, E.). the vector 14 and the covariance matrix Eg have analytic solutions (Särkkä and Solin, 2019]. More explicitly, we can construct m, through
X Ve
T Eo + E
L L
=
= pi
+
where co. 61 L)) and L, =
22
(10)
Dimctton: Forvani
Desctvm Backwand
where L, is computed following the Cholesky decomposition of the covariance matrix. The expressions for the mean vec- tor and the covariance matrix are in Ap- pendix C. We conclude this section on a theoretical analysis of our methodology converging to the unique multi- marginal Schrödinger Bridge solution by iteratively optimizing Eq. (3), and Eq. (9).
0.0
-0.5
-1.0
-1.5
13-20-15 10-05 05 10 15
Figure 4: Iterative propagation of dynamics using the optimal conditional acceleration to approximate the locity profile n,, (2,, at the intermediate marginals.
Convergence We establish the convergence of our alternating scheme - between interpolating path optimization and the coupling refinement the unique mmSB solution. Let denote with P the path measure associated with the learned dynamics, and with T' the path measure of the learned dynamics at the 7hh iteration of our algorithm.
Theorem 3.5. Under mild assumptions, our iterative scheme admits fixed point solution P Le.. KL(PP(P*)) - 0, and in particular, this fixed point coincides with the unique PramS
We provide an alternative claim the convergence proof [Shi et al., 2023]. We establish conver- gence to the global minimizer, based on optimal control principles and the monotonicity of the KL divergence, expressing Eq. (2) as KL divergence, due to the Girsanov Theorem [Chen et al., 2019].
3.3 Training Scheme
summary of our training procedure presented in Alg. The first step of our alternating matching algorithm is to compute the conditional multi- marginal bridge, fixing the parameterized coupling qº (f=,)) }) and sampling collections of points {Tn} ~ d ([z,)). Furthermore, the initial velocity distribution to ~ {u( = To( u)dz is needed, which is unknown in practice for most applications. To address this, following [Chen et al., 2023a], we initialize the velocity to NO T) and iteratively propagate the dynamics with the conditional acceleration in Eq. (6) for the given samples (lines 4-9 in Alg. 1). This approximates the true conditional distribution n, (v, via Langevin-like like dynamics. Notably, this process requires only few iterations (at most 10 iterations in our experiments), and does not involve backpropagation, thus adding negligible computational cost even in high dimensions. This iterative propagation yields the optimal conditional acceleration in Eq. (6), which induces the optimal conditional path. Solving the 3MBB for every set of {.,) enables to marginalize and construct the marginal path Pt- Subsequently, and match the parameterized acceleration (t, m2) minimizing the variational objective in Eq. (9). This optimization induces refined joint distribution o ((x,)) by propagating the dynamics through the augmented SDE in Eq. (2), which will be used in the first step of the next training iteration.
fix the marginal path

---

SBIRR
Smooth
SBIRR
Smooth SB
t3
to
t2
Time Marginal
22
Time Marginal
[3
MMFM
3MSBM((urs)
MMFM
3MSBM(Ourr)
t
t5
ts
t
t7 B
t,
t3
Trajectory comparison L.V among SBIRR, MMFM, Smooth SB, and our 3MSBM.
Figure 5:
Figure 6: Trajectory comparison on GoM among
SBIRR, MMFM, Smooth SB, and our 3MSBM.
Table 2:
Mean and Standard Deviation SWD distances over 5 seeds at the left-out marginals, and average SWD from the rest points belonging to the training set on the LV data for MIOFlow, SBIRR, MMFM, Smooth SB, and 3MSBM (Lower is better).
Method
SWD 1
SWD t3
SWD ts
SWD tr
Rest SWD
MIOFlow
1.53
1.49
1.22:0007
1.46::0.06
0,93:0004
SBIRR
0.17:0003
0.163
0.24
0.48
0.18
MMFM
0.21#0002
0.25
0.39
0.57=00.05
0.22
DMSB
0.64
0.67=0.06
0.98
0.63++0.04
0.43
Smooth SB
0.29
0.18
0.11
0.37
0.23
3MSBM
0.23+50.02
0.13
0.15
0.36
0.18
4 Experiments
We empirically validate the computational and performance benefits of our method. The simulation- free training scheme of our algorithm suggests that we avoid costly numerical simulations and approximation errors, consistent with benefits observed by prior matching methods [Shi et al., 2023, Liu et al., 2024]. For this reason, our algorithm capable of preserving high scalability while maintaining accuracy, demonstrated in the experiments below. We evaluate the performance of 3MSBM on synthetic and real-world trajectory inference tasks, such as Lotka- Volterra (Sec. 4.1), ocean current in the Gulf of Mexico (Sec. 4.2), single-cell sequencing (Sec. 4.4), and the Beijing air quality data (Sec. 4.3). We compare against state-of-the-ar methods explicitly designed to incorporate multi-marginal settings, such Deep Momentum Murti-Marginal Schrödinger Bridge (DMSB: Chen et al. [2023a]). Schrodinger Bridge with Iterative Reference Refinement (SBIRR; Shen et al. [2024]), smooth Schrodinger Bridges (smoothSB; Hong al. [2025]), and Multi-Marginal Flow Matching (MMFM; Rohbeck et al. [2024]), and against one additional Neural ODE-bases method: MIOFlow (Huguet et al., 2022], using the same metric in all datasets: the Sliced Wasserstein Distance (SWD). Additional results for these tasks with more metrics and baselines are provided in Appendix E.
4.1 Lotka- Volterra
We first consider a synthetic dataset generated by the Lotka-Volterra (LV) equations [Goel et al., 1971], which model predator-prey interactions through coupled nonlinear dynamics. The generated dataset consists of 9 marginals in total; the even-numbered indices are used to train the model (i.e., to, t2,,t, to. tg), and the remainder of the time points are used to assess the efficacy of our model to impute and infer the missing time points. In this experiment, we benchmarked 3MSBM against MIOFlow, SBIRR, MMFM, DMSB, and Smooth SB. Performance was evaluated using the SWD distance to the validation marginals to measure imputation accuracy, and the SWD distance to the training marginals assess how well each method preserved the observed data during generation. The results in Table 2 and Figure 5 show that 3MSBM outperforms the baseline models in inferring the marginals the missing points, yielding the lowest deviation from most left-out marginals, while also generating trajectories which preserve the training marginals more faithfully, as shown by the lower average SWD distance over the remaining points belonging to the training set. Additional results- with more metries and baselines-aue provided in Appendix E.2.
7

---

Table 3: Mean and SD over 5 seeds for SWD distances on the GoM for MIOFlow, SBIRR, MMFM, DMSB, Smooth SB, and 3MSBM (Lower is better).
Method
SWD 1
SWD 03
SWD l9
SWD 4
Rest SWD
MIOFlow
0.83 10.06
0.34 10.03
1.23 10.09
0.96 10.06
0.19 10.03
SBIRR
0.15 +0.03
0.11 +0002
0.11 #0.03
0.09 +0.04
0.13 #0.04
MMFM
0.23 +0.04
0.25 +0.08
+0.03
0.19 +0.04
0.14 +0.03
DMSB
0.22 +0.02
0.54 +0.04
0.39 +0.02
0.28 +0.04
0.09 +0.03
Smooth SB
0.17 +0.02
0.14 #0.04
0.10 #0.02
0.13 +0.02
0.08 10.01
3MSBM
0.14 +.0.02
0.14 +:002
0.08 +0.01
0.06 10.01
0.05 +0001
Table Mean and SD over 5 seeds for SWD distances on the Beijing air quality data for MMFM, DMSB, and 3MSBM (Lower is better).
Method
SWD f/
SWD ts
SWD iy
SWD t
Rest SWD
MMFM
17.51 +2.41
23.94 +1.9
32.56
2299
39.98
+3559
30.25 =2.17
DMSB
21.10
+2.65
21.92
1177
35.53
1388
35.75
14.13
33,42 2.42
3MSBM
17.70 +1.93
9.78 +1.58
22.23
+ 3.64
32.23
+3776
21.25 +1.63
Table 5: Mean and SID over 5 seeds for the MMD and SWD on Embryoid Body (EB) dataset for SBIRR, MMFM, DMSB, and 3MSBM (Lower is better).
Method MMD 1
SWD 1
MMD iy
SWD fy
Rest MMD
Rest SWD
SBIRR
0.70
0.80
0.73=:0.06
0.91
0.47=00.05
0.66:00.07 0.52+0.007
MMFM
0.37=0.02
0.59+
0.35+
0.76++0.04
0.22++0.02
DMSB
0.38
0..06
0.36
0.54
0.14
0.45
3MSBM(ours)
0.48 1..04
0.14
0.11 0.003
0.36
4.2
Gulf of Mexico
Subsequently, we evaluate the efficacy of our model to infer the missing time points in real-world multi- multinmarriial dataset. The dataset contains ocean-current snapshots the velocity field around a vortex in the Gulf of Mexico (GoM). It includes total of 9 marginals, with even- wenn-iddeedd time points (i.e., to, f2, fa, f6. t&) used for training, and the remaining are left out to evaluate the model's ability to impute and infer missing temporal states. For this experiment, we compared our 3MSBM against MIOFlow, SBIRR, MMFM, DMSB, and Smooth SB. The metrics used to evaluate performance were: SWD distance from the left-out points, and the mean of SWD distance from the training points, capturing how well the generated trajectories of each algorithm preserve the marginals that comprised the training set. Figure shows that 3MSBM generates smoother trajectories with more accurate recovery of the left-out marginals. In comparison, SBIRR produces noisier, kinked trajectories; MMFM struggles to capture the dynamics, leading larger deviations from the left-out marginals, while Smooth SB achieves the smoothest trajectories among the baselines. These observations are further confirmed by Table 3, where 3MSBM achieves the lowest SWD distances for most validation points and better preserves the training marginals. Additional results are provided in Appendix E.3.
4.3 Beijing Air Quality
To further study the capacity of 3MSBM to effectively infer missing values, we also tested it in the Beijing multi-site air quality data set [Chen, 2017]. This dataset consists of hourly air pollutant data from 12 air-quality monitoring sites across Beijing. We focus on PM2.5 data, an indicator monitoring the density of particles smaller than 2.5 micrometers, between January 2013 and January 2015, across 12 monitoring sites. We employed a slightly different setup than Rohbeck et al. [2024]. We focused on single monitoring site and aggregated the measurements collected within the same month. To introduce temporal separation between observations, we selected measurements from every other month, resulting in 13 temporal snapshots. For the imputation task, we omitted the data at tz2 ts. fs. and f1 while the remaining snapshots formed the training set. We benchmarked our 3MSBM method against MMFM with cubic splines and DMSB. Table 4 shows that 3MSBM achieved overall better imputation accuracy, yielding the smallest Sliced Wasserstein Distance (SWD) distances, while also better preserving the marginals consisting of the training snapshots compared to the baselines. Additional details and results more metries are left for Appendix E.4.
8

---

Ground Truth
SBIRR
DMSB
MMFM
3MS
to
PC2
13
Time Marginal
PCI
PCI SBIRR
PCI
Ground Truth
DMSB
MMFM
3MSBM(Ours)
PC4
t
PC3
PC3
PC3
PC3
PC3
Figure 7: Comparison of the evolution of the EB dataset on the 100-dimensional PCA feature space among SBIRR, MMFM, DMSB, and 3MSBM.
1-point Bridge
4-point Bridge
0.3
Time Marginal
07
022
10"
27
37522
0.0
Number of Pinned Points
Conditional bridges on GoM currents using LEFT: pinned point, RIGHT: 4 pinned points.
Figure 8:
Figure 9: Mean Wi2 distance of left-out marginals for vary- ing number of pinned points.
Figure 10: Exponential Decay of [Xjl coefficients for increas- ing number pinned points.
4.4 Single cell sequencing
Lastly, we demonstrate the efficacy 3MSBM to infer trajectories in nigh-dimensional spaces. In particular, we use the Embryoid Body (EB) stem cell differentiation tracking dataset, which tracks the cells through stages over 27-day period. Cell snapshots are collected five discrete day-intervals: to [0, t E [6,9], f2 E [12, 15], 13 E [18,21] 21], 4 [24, 27]. The training set consists of the even-indexed time-steps (i.e., to, 12, (4), while the rest are used as the validation set. We used the preprocessed dataset provided by [Tong et al., 2020, Moon et al., 2019], embedded in 100-dimensional principal component analysis (PCA) feature space. We compare 3MSBM with SBIRR, MMFM, and DMSB, evaluating performance using Sliced Wasserstein Distance (SWD)-u as in prior tasks and additionally Maximum Mean Discrepancy (MMD). As in previous experiments, assessed the quality of imputed marginals and the preservation of training marginals. As shown in Table 5, 3MSBM consistently outperforms existing methods across all metrics, achieving significantly more accurate imputation of missing time points and recovering population dynamics that closely match the ground truth, as illustrated in Figure/ 7. While DMSB also generated accurate PCA reconstructions (Figure 7), is noted that required approximately 2.5 times more training time than 3MSBM. A detailed comparison of the resource requirements differences between DMSB and 3MSBM given in Sec. 4.6. Further single-cell sequencing setup and expanded results, with additional baselines and more metrics, are deferred to App. E.5.
4.5 Ablation study number of pinned points
We ablate the number of pinned points used in Eq. (6). by modifying the linear combination as follows: Σ Xyā starting from K = n + 1, namely including only the nearest fixed point, and incrementally adding more up K N. Figure 14 demonstrates that incorporating multiple pinned points significantly improves performance on the GoM dataset compared to using only the next point, enabling better inference of the underlying dynamics, as illustrated Figure 8, albeit these benefits plateau beyond K = + 3. This phenomenon is further explained by the exponential decay of coefficients |Ayl for increasing number, depicted in Figure 10, thus rendering distant points negligible. Notice that for the bridge using the next pinned points, it holds |>l = 1>2 as found in Section 3.1. Consequently, practical implementations can adopt a truncated conditional policy, considering only the next k pinned points, thereby significantly improving efficiency without sacrificing accuracy.

---

Percernagee Training Time (%)
100
DMCB
Percentage of Allocated Merrory
100
L
40
40
20
Figure 11: Training time percentage (%) com- parison between 3MSBM and DMSB
Figure 12: Allocated memory percentage (%) comparison between 3MSBM and DMSB
Table Per-epoch training time (sec) increasing the marginals N, and the percent change [%)].
Table 7: Percentage of time spent in 3MBB and the Bridge Matching.
Num. of marginals
3
4
5
Dataset
Cond. Bridge
Bridge Matching
Time [s]
773
789
812
EB-100
3.3%
96.7%
Percent Increase [%]
2.1
5.1
GoM
7.4%
92.6%
4.6 Scalability of 3MSBM
We empirically demonstrate the scalability of 3MSBM. Notably, the matching- based training with the closed-form conditional bridge from 3MBB (Eq. obviatig numerical peegration-removve both computational overhead and approximation error.
Computational Resources Since 3MSBM and DMSB address the same problem, i.e the mmSB problem in Eq. 2. through different training scheme (matching-baseed and IPF-based respectively), we report the resources needed by each method. In particular, Figures 11 and 12 demonstrate that our 3MSBM faster wall-clock time on every dataset, while also requiring significantly less memory casily handling the high-dimenn dighedimensional single-cell sequencing task.
Ablation on the number of marginals Next, the number of marginals on EB-100 and GoM is ablated, i.c. increasing the total number of marginals N, while holding all other hyperparameters fixed (e.g., NEE, batch size, model size), and report per-epoch training time. In Table is shown that varying N has negligible impact on the per-epoch training time, indicating that 3MSBM scales well with the number of marginals. This insensitivity stems from our algorithmic design, as the analytic form of our conditional-bridd step is independent of N, and crucially.
Comparison between 3MBB and Matching Finally, we present breakdown of the time percent- age attributed in each of the two steps of our algorithm: i) the iterative propagation of the conditional dynamics (lines 4-9 Alg. 1) and ii) the Bridge Matching step (line 10 in Alg. 1) the EB-100 task. Table 7 demonstrates that the computational complexity introduced by the itcrative propagation of the conditional dynamics, remains negligible compared to the Bridge Matching step.
Conclusion and Limitations
In this work, we developed 3MSBM, novel matching algorithm that infers temporally coherent trajectories from multi-snapshot datasets, showing strong performance in high-dimensional settings and many marginals. Our work paves new ways for learning dynamic processes from sparse temporal observations, addressing a universal challenge across various disciplines. For instance, in single-cell biology, where can only have access to snapshots of data, our 3MSBM offers principled way to reconstruct unobserved trajectories, enabling insights into gene regulation, differentiation, and drug responses. While 3MSBM significantly improves scalability over existing multi-marginn methods, certain limitations remain. In particular, its effectiveness in densely sampled high-dimensional image spaces, such those encountered in video interpolation, remains unexplored. As future work, aim to extend the method to capture long-term dependencies large-scale image settings.
10

---

Acknowledgements
The authors would like to thank Tianrong Chen for the fruitful discussions and insightful comments that helped shape the direction of this work. This research is partially supported by the DARPA AIQ program through the DARPA CMO contract number HR00112520010. We would like to thank Dr. Pat Shafto, AIQ Program Manager, for useful technical discussions. Augustinos Saravanos acknowledges financial support by the A. Onassis Foundation Scholarship.
References
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-6851, 2020.
Yang Song, Jascha Sohl-Dickstein, Diederik Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:201 1. 133466, 2020.
Brian DO Anderson. Reverse-time diffusion equation models. Applications, 12(3):313-326 1982.
Stochastic Processes and their
Pascal Vincent. A connection between score matching and denoising autoencoders. Neural compula- lion, 23(7):1661-1674, 2011.
Yuyang Shi, Valentin De Bortoli, Andrew Campbell, and Arnaud Doucet. Diffusion schrödinger bridge matching, 2023.
Cédric Villani et al. Optimal transport: old and new, volume 338. Springer, 2009.
Erwin Schrödinger. Über die unkehrung der naturgesetze. Verlag der Akademie der Wissenschaften in Kommission bei Walter De Gruyter u..., 1931.
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. Advances in neural information processing systems, 26, 2013.
Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications data science. Foundations and Trends@ in Machine Learning, 116--6):355-677, 2019.
Francisco Vargas, Pierre Thodoroff, Austen Lamacraft, and Neil Lawrence. Solving schrüdinger bridges via maximum likelihood. Entropy, 23(9): 134, 2021.
Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXiv:22/0.02747, 2022.
Xingchao Liu, Lemeng Wu, Mao Ye, and Qiang Liu. Let us build bridges: Understanding and extending diffusion generative models, 2022a.
Nikita Gushchin, Sergei Kholkin, Evgeny Burnaev, and Alexander Korotin. Light and optimal schrödinger bridge matching. In Forty-first International Conference on Machine Learning, 2024a.
Stefano Peluchetti. Diffusion bridge mixture transports, schrödinger bridge problems and generative modeling. Journal of Machine Learning Research, 24(374):1-51, 2023.
Guan-Horng Liu, Yaron Lipman, Maximilian Nickel, Brian Karrer, Evangelos Theodorou, and Ricky TQ Chen. Generalized Schrödinger bridge matching. In International Conference Learning Representations, 2024.
George Rapakoulias, Ali Reza Pedram, Fengjiao Liu, Lingjiong Zhu, and Panagiotis Tsiotras. Go with the flow: Fast diffusion for gaussian mixture models. arXiv preprint arXiv:2412.09055, 2024.
Tianrong Chen, Guan-Horng Liu, Molei Tao, and Evangelos Theodorou. Deep momentum multi- marginal schrödinger bridge. Advances Neural Information Processing Systems, 36:57058- 57086, 2023a.
11

---

Hugo Lavenant, Stephen Zhang, Young Heon Kim, Geoffrey Schiebinger, et al. Toward math- ematical theory of trajectory inference. The Annals of Applied Probability, 34(1A):428-006 2024.
Young-Heon
Yongxin Chen, Giovanni Conforti, Tryphon T Georgiou, and Luigia Ripani. Multi-marginal schrödinger bridges. In International Conference on Geometric Science of Information, pages 725-732. Springer, 2019.
Tim Dockhorn, Arash Vahdat, and Karsten Kreis. Score-based generative modeling with critically- damped langevin diffusion. arXiv preprint arXiv: 2021.
Grace Hui Ting Yeo, Sachit Saksena, and David K Gifford. Generative modeling of single- cell time series with prescient enables prediction of cell trajectories with interventions. Nature communications, 12(1):3222, 2021.
Jiaqi Zhang, Erica Larschan, Jeremy Bigness, and Ritambhara Singh. scnode: generative model for temporal single cell transcriptomic data prediction. Bioinformatics, 40(Supplement 2):i146 ii154, 2024a.
Christian LE Franzke, Terence J O'Kane, Judith Berner, Paul D Williams, and Valerio Lucarini. Stochastic climate theory and modeling. Wiley Interdisciplinary Reviews: Climate Change, 6(1): 63-78, 2015.
Rytis Kazakevičius, Aleksejus Kononovicius, Bronislovas Kaulakys, and Vygintas Gontis. Under- standing the nature of the long-range memory phenomenon in socioeconomic systems. Entropy. 23(9):1125, 2021.
Yunyi Shen, Renato Berlinghieri, and Tamara Broderick. Multi-marginal schr\," odinger bridges with iterative reference refinement. arXiv preprint urXiv:2408.06227, 2024.
Martin Rohbeck, Edward De Brouwer, Charlotte Bunne, Jan-Christian Huetter, Anne Biton, Kelvin Y Chen, Aviv Regev, and Romain Lopez. Modeling complex system dynamics with flow matching across time and conditions. In The Thirteenth International Conference on Learning Representa- tions, 2024.
Justin Lee, Behnaz Moradijamei, and Heman Shakeri. multi-marginal stochastic flow matching for high-dimensional snapshot data irregular time points. arXiv preprint arXiv:2508 04351, 2025.
Wanli Hong, Yuliang Shi, and Jonathan Weed. Trajectory inference with smooth schr)." odinger bridges. arXiv preprint arXiv:2503.00S: 2025.
Yongxin Chen, Tryphon T Georgiou, and Michele Pavon. On the relation between optimal transport and schrödinger bridges: A stochastic control viewpoint. Journal of Optimization Theory and Applications, 169:671 1-61, 2016.
Ivan Gentil, Christian Léonard, and Luigia Ripani. About the analogy between optimal transport and minimal entropy. In Annales de Faculté des sciences de Toulouse: Mathématiques, volume 26, pages 569-600, 2017.
Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to the monge- kantorovich mass transfer problem. Numerische Mathematik, 84(3):375-385 3303,,33333033 2000.
Nikita Gushchin, Daniil Selikhanovych, Sergei Kholkin, Evgeny Burnaev, and Alexander Korotin. Adversarial schr\." odinger bridge matching. arXiv preprint 2rXiv:2405.14449 2024b.
Valentin De Bortoli, Guan-Horng Liu, Tianrong Chen, Evangelos A Theodorou, and Weilie Nie.
Augmented bridge matching. arXiv preprint arXiv:2311.06978, 2023.
Christian Léonard. survey of the schr\ " odinger problem and some of its connections with optimal transport. arXiv preprint arXiv: 15, 2013.
Panagiotis Theodoropoulos, Nikolaos Komianos, Vincent Pacelli, Guan-Horng Liu, and Evangelos Theodorou. Feedback schr)," odinger bridge matching. arXiv preprint Att.:2440.14053 2024.
12

---

Yongxin Chen and Tryphon Georgiou. Stochastic bridges of linear systems. IEEE Transactions on Automatic Control, (20):5566-331,01
Simo Särkkä and Amo Solin. Applied stochastic differential equations, volume 10. Cambridge University Press, 2019.
Guillaume Huguet, Daniel Sumner Magruder, Alexander Tong, Oluwadamilola Fasina, Manik Kuchroo, Guy Wolf, and Smita Krishnaswamy. Manifold interpolating optimal-transpor flows for trajectory inference. Advances in neural information processing systems, 35:29705-297188 2022.
Narendra S Goel, Samaresh C Maitra, and Elliott W Montroll. On the volterra and other nonlinear models of interacting populations. Reviews of modern physics, 43(2):231, 1971.
S Chen. Beijing multi-site air quality, uci machine learning repository, 2017.
Alexander Tong, Jessie Huang, Guy Wolf, David Van Dijk, and Smita Krishnaswamy. Trajectorynet: A dynamic optimal transport network for modeling cellular dynamics. In International conference on machine learning, pages
PMLR, 2020.
Kevin R Moon, David Van Dijk, Zheng Wang, Scott Gigante, Daniel B Burkhardt, William S Chen, Kristina Yim, Antonia van den Elzen, Matthew J Him, Ronald R Coifman, al. Visualizing structure and transitions in high-dimensional biological data. Nature biotechnology. 37(12): 1482-1492, 2019.
Valentin De Bortoli, James Thortoon, Jeremy Heng, and Arnaud Doucet. Diffusion schrödinger bridge with applications to score-based generative modeling. Advances in Neural Information Processing Systems, 34: 177099977779, 2021.
Tianrong Chen, Hoonn Liu, and Evangelos A Theodorou. Likelihood training of schr\ odinger bridge using hawwadd-bckwaadd sdes theory. arXiv preprint arXiv:21 1221,
Michele Pavon, Giulio Trigila, and Esteban G Tabak. The data-driven schrödinger bridge. Communi- cations on Pure and Applied Mathematicx, 1407):1545-1533, 2021.
Marcel Nutz. Introduction to entropic optimal transport. Lecture notes, Columbia University, 2021.
Promit Ghosal, Marcel Nutz, and Espen Bernton. Stability of entropic optimal transport and schrödinger bridges. Journal of Functional Analysis, 28969) 2022.
Flavien Léger. A gradient descent perspective on sinkhorn. Applied Mathematics & Optimization, 84 (2):1843-1855, 2021.
Kirill Neklyudov, Rob Brekelmans, Daniel Severo, and Alireza Makhzani. Action matching: Learning stochastic dynamics from samples. In International conference machine learning, pages 25858 25889. PMLR, 2023.
Linqi Zhou, Aaron Lou, Samar Khanna, and Stefano Ermon. Denoising diffusion bridge models. arXiv preprint axXiv:2309.66948 2023.
Michael S. Albergo, Nicholas M. Boffi, and Eric Vanden-Fijndenn Stochastic interpolants: A unifying framework for flows and diffusions, 2023.
Alexander Tong, Nikolay Malkin, Kilian Fatras, Lazar Atanackovic, Yanlei Zhang, Guillaume Huguet, Guy Wolf, and Yoshua Bengio. Simulation -fee schr\ odinger bridges via score and flow matching. arXiv preprint arXiv:2307.03672, 2023a.
Vignesh Ram Somnath, Matteo Pariset, Ya-Ping Hsieh, Maria Rodriguez Martinez, Andreas Krause, and Charlotte Bunne. Aligned diffusion schrödinger bridges. In Uncertainty Artificial Intelli- gence, pages 1985-19999 PMLR, 2023.
Guan-Horng Liu, Arash Vahdat, De-An Huang, Evangelos A Theodorou, Weili Nie, and Anima Anandkumar. SB Mogg-tt--mmage Schrödinger bridge. 2023.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
13

---

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. 00482, 2023b.
Kacper Kapusniak, Peter Potaptchik, Teodora Reu, Leo Zhang, Alexander Tong, Michael Bronstein, Joey Bose, and Francesco Di Giovanni. Metric flow matching for smooth interpolations on the data manifold. Advances in Neural Information Processing Systems, 37:13501 1-135042, 2024.
Quan Dao, Hao Phung, Binh Nguyen, and Anh Tran. Flow matching in latent space. arXiv preprint arXiv:2307.08690 2023.
Xingchao Liu, Chengyue Gong, and Qiang Liu. Flow straight and fast: Learning to generate and transfer data with rectified flow. arXiv preprint arXiv:2209.0300 2022b.
Amartya Banerjee, Harlin Lee, Nir Sharon, and Caroline MoosmåAller. Efficient trajectory inference in wasserstein space using consecutive averaging. arXiv preprint arXiv:2405. (rxXi:2405.19679 2024.
Zhenyi Zhang, Tiejun Li, and Peijic Zhou. Learning stochastic dynamics from snapshots through regularized unbalanced optimal transport. arXiv preprint arXiv:2410.00844, 2024b...
Tianrong Chen, Jiatao Gu, Laurent Dinh, Evangelos A Theodorou, Joshua Susskind, and Shuangfei Zhai. Generative modeling with phase stochastic bridges. arXiv preprint arXiv:231.007005, 2023b.
Evangelos Theodorou, Jonas Buchli, and Stefan Schaal. A generalized path integral control approach to reinforcement learning. The Journal of Machine Learning Research, 11:3137-3181, 2010.
14

---

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: Yes, the claims in the abstract and the contributions listed in the introduction reflect our paper's main contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in 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.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Yes, in Section clearly state the limitations of work.
Guidelines:
The answer NA means that the paper has no 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 separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model sell-specification, asymptotic approximations only holding locally). The authors should reflect 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 is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to 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 fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
15

---

Justification: Yes, Section B in the Appendix provides detailed and rigorous proofs for every theoretical result introduced in our paper.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proof's provided in appendix 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]
Yes we disclose all necessary information to reproduce our experimental results. In particular, we provide a detailed pseudocode of our methodology in Section 3, and in Section E in the Appendix we list all the hyperparameters for every conducted experiment.
Justification:
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 a dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, 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 often one good way to accomplish this, but reproducibility can also be provided via detailed 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, or other means that are appropriate to the research performed.
While NeurlPS 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 a 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 be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open- open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility 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, 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.
16

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Question:
Answer: [No]
Justification:
We have not released the code yet, but we plan to do so very soon.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https: /ttp:///nis..ooo cc/
public/ guide Soodeeudeeeis ss s onP olicy) for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is 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 NeurIPS code and data submission guidelines (httpa //nips. cc/ httpp//blcccccc/pppppppl/cccaclt cealli/cgiiie
Code Submis ss wwwisssssoooooooooo onooonnooiiooo illyyy for more details.
The authors should provide instructions on 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 much information 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: Section 4 in the main paper deseribes the experimental pipeline adopted for every experiment conducted in this paper to a level sufficient to reproduce our experimental results. The full details regarding hyperparameter selections are left for the Appendix E.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to 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: Yes, our experiments were conducted for 5 different seeds, giving us standard deviation intervals for each metrics.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
17

---

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 run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to 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 not show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: We describe the computational resources needed to conduct our experiments in the Appendix E.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate 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 as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn't 't make 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 https: /tt/rieurpppp wwtpcppplc/cccctccccccccccccccce ctt.:/ttiiccuiieeeee ines?
Answer: [Yes]
Justification:
We adhere to the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. if there special consid- eration 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: [Yes]
Justification: Yes, in Section 5, clearly state the broader impact of our work.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
18

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), faimess 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 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 an 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 to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks
Guidelines:
The answer NA means that the paper poses 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 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 challenging, and many papers do not require this, but encourage authors to take this into account and make best. faith effort.
12. Licenses for existing assets
Question: 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?
Answer: [Yes]
Justification: Datasets and code for the benchmark models used for the scope of this work have been properly cited and acknowledged.
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a 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.
19

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, papersw snnersw: swithcode com// datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information 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: [NA]
Justification: This paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing 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 details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information 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 to the NeurIPS Code of Ethics, 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 institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending 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.
20

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
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, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does involve LLMs as any important, original, or non-standard components.
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 our LLM policy (https: /tt////eeipppp/eeeeee http://uursscc///aaa.ieoeeeee wwt..eceeeneeeceeeeeeeeaeae 2025/ LLM) for what should or should not be described.
21

---

Contents
A Summary of Notation
22
Proofs
23
B.1 Proof of Theorem 3.1
23
B.2 Proof of Proposition 3.2
25
B.3 Proof of Proposition 3.4
28
B.4 Proof of Theorem 3.5
29
Gaussian Path
Extended Related Works
33
E Additional Details Experiments
34
E.1 General Information
34
E.2 Lotka-Volters
34
E.3 Gulf of Mexico
35
E.4 Beijing air quality
35
E.5 Single sequencing
36
E.6 Ablation study on
37
A Summary of Notation
In the Table below, we summarize the notation used throughout our work.
Table 8: Notation
Time coordinate
Ft
Position
Ut
Velocity
dg
Acceleration
Augmented Variable
Soft marginal constraint
gn (xx) (on)
nh Positional marginal th Velocity marginal nh Augmented marginal Augmented probability path Coefficients of impact of future pinned points of the nhth segment
n (n vn) vi) A(a)
(L))
Functions shaping the cond. bridge of the th segment Coupling over n marginals
{an) V.(mm)
Value Function
P
Value function second and first order approximations
p E
Gaussian path mean vector
Gaussian path covariance matrix
L
Cholesky decomposition the cov. matrix
22

---

B Proofs
B.1 Proof of Theorem 3.1
Theorem B.1 (SOC representation of Multi-Marginal Momentum Brownian Bridge (3MBB)). Con- sider the following momentum system interpolating among multiple marginals
[0] aj
arg min !'"
1 211ae1" dt 2(m,
"
(mn mnn T R(mm m,))
(11)
2
n=1
8.1 dm: = Am + + gdW, mg = mo-
(12)
We define the value function as Ve(m))=== bmm P. P mi Pr where Pt, rt are the second- and first-order approximations, respectively. This formulation admits the following optimal control expression ug (m²) -gg" P,- P,-(mm, + TT)). For the multi-marginal bridge with {,a)) J fixed at {tn).
for n E (0. N), the dynamics of Pe and It obey the following backward ODEs
P. = APi + P.A -
Pn = =P + +R) with PN+ = =0
(13)
à1- - Ari,
Tn (r (n
n Tut Rmm))
(14)
where Pa+ == lim Pt, and rut == lim rs. for t E & (t1, V tt2,
Proof.
We start our analysis by considering the second-order approximation of the value function:
1 (t,, 2 m" 2.m, + r.m.
(15)
where Qr and ra serve as second and first-order approximations. From the Bellman principle and application of the Ito's Lemma to the value function, we obtain the Caciioon--aab-----eeeeeeee (HJB) Equation:
V1 min - EE [leell 2dt]] + (Am, + =w) +
1 TT nom gg")
(16)
2
Solving for the optimal control uf, we obtain:
u Vmm = - ( + r()
(17)
Thus plugging it back to the HJB, we can rewrite it
V1 Vm VM + VT Am + (rrVVmeeee = 0
VI
(18)
2
Recall the definition for the value function
V(t, m. Qm+ r:Qm,
Substituting for the definition of the value function in the HJB yields the following PDE
2 mm m
Qgg -m. 000m::00000000mmm0mm0mm00mm000/0
Tr(V
0
(19)
Grouping the terms of the PDE quadratic in m,, we obtain the Riccati Equation for Q,
ATQ + QA Qgg
(20)
Then, grouping the linear terms yields
r +
- - 1:0 +++AQA=
(21)
Now, notice that substituting the Ricatti in Eq. 20 into Eq. 21, one obtains
r, = -r,
(22)
The solution of this ODE is
= (t,,s))..
(23)
23

---

where (t,, s) is the state transition function of the following dynamics dm, m,t) from t to and is defined as (tt, 8)
t-s 1
Finally, we define P(t) = Q(t)=1, and modify the Ricatti
in Eq. 20 follows
= ATQ Q''QQQ -QgT-
Q QQ = QQQQQQQ +Q +Q-'QQQQQQ 'QQQ Qqq
(24)
P = AP + PAT gg
yielding the Lyapunov equation
p A PAT gg
(25)
Therefore, have proved the desired ODEs for the first and second-order approximations Tı, Pt of the value function. Now, we have to determine the expressions for the terminal conditions for the ODEs in each segment.
Note that follows from the dynamic principle that these ODEs are backward, therefore ensuing segments will affect previous. These terminal conditions carry the information from the subsequent segments. More specifically, assume multi-marginal process with {mp+1} pinned at (ln+1) J. Now let us consider the two segments from both sides:
Segment for / (Ins In+1]
Segment + 1: fort E [tx+11 tn+2
To solve Pt, r; for Segment n. we have to account for the effect of Segment n + 1. To compute the value function at tn+1, accounting for the impact from Segment +
V Vn+++ + (ma+1 - m1+1) )) Ra41 (ma+1 mx+1)
(26)
where the first term encapsulates the impact from Segment n + 1. More explicitly, for the Segment + 1, at l = Ln+1, we define the value function
Vat++ = m n+++ Qu+1+ + Tn+1+ Qn+ mm42+
(27)
m 1+1+ Q M 1 + ((n + 1, 2)m mm+21
Therefore, for Segment the value function at the terminal time t,+1 is given by
Vn+1 = Vuti+ 1011=V+t++
= (m - mm++1) Bu+1 (m041 ma+1))
1
m (Qn+1t + R + + (1 Qn+1+
(28)
ma+1 Ru+1 M + const.
1 T
m ++1 2n+1 + Tn+1 1n+1Qn-1mn+t + const.
2
This suggests that the terminal constraints for the Ricatti equation Qt, and the reference dynamics vector re are given by
Qn+1 = Qu+1+ + Ru+1
(29)
Tn+1 == (#(t=+1, TMm - mn 1
(30)
Lastly, recall that the Lyapunov function is defined P. Q- thus the terminal constraints with respect to the Lyapunov function are given by
P (P--1+ + 11 Tn+1
P Pati (P-1 1+7+11
(31)
This implies that and hence the optimal control would depend on all preceding pinned points after (my. tú 2 t). Finally, notice that for the last segment t f.N-i,tN], holds that Px+ since there is no effect from any subsequent segment, hence Pty RØy- and rty simplifies to -max
D
24

---

B.2 Proof of Proposition 3.2
0
Proposition B.2. Let R =
④ c
the limit when e - the conditional acceleration Eq.
admits analytic form:
N
amel tn+1 2 = (+) - En+1) C2 (3()
Σ j=
XoZz. E tn+1)
(32)
where (In+1 : fn+1 2 1) signifies the bridge is conditioned on the set of the ensuing points, Xj are static coefficients and CP (t), (t), C3 are time-varying coefficients specific to each segment.
Proof. We start our analysis with the last segment N - 1, and then move to derive the formulation for an arbitrary preceding segment n.
Segment - 1: E [1N-
For the last segment, the terminal constraint t E 11 T) is given by Ps Ry' and hence TN simplifies -m.. Solving the backward differential equation Pt, for [tx-1. 7), with Pr PN Ry yields
0² (t (T3 (t-T) 2
(t- + (t-7)
P(()
(t - T)² (-TT)
(33)
-0²(t 1)) +
Additionally, solving for r; for E 1 71. with rr = -m yields:
Y2 (1, T)mr
(34)
where 4 (t.. .) is the transition matrix of the following dynamics dm, = A(t)m,dt, and is defined as Plugging Eq. 33, and 34 into Eq. 17 yields
1 t-s (1, 0
uy 3 TT- (ūT - xx) 3 7t
Vt [NN-1,
(35)
Ut
Therefore, regardless of the total number of marginals the C-functions for the last segments
always: C C(N--)(( = (-- C² ((N)))(( T-t Cg3g CN(N-1) (t) =t)
3
Segment nt I E [Ln, (n+1]
Now, move derive the conditional acceleration, for arbitrary segment n, with N i.e.
is not the last segment. Let us recall the optimal control formulation from Eq, (17)
ui (m²) = -gg P. (m, + (tt, (x+1)rn+1)
(36)
For convenience, let us define the following functions corresponding to the nh segment: E [tm, fat1)
2 (") (t) 3t 3t 33
(n) (t) 31 4 3f
[n] (l) = 61 - 6/ - 3
(m) (1) = 61 - -6/n+1 4
23
(37)
(n) (t)
4f (n) (t) =
25
41,+1
26
27 (t) = 6t 6t, +3
28
(m) (t) 6t 6t, 4
For this arbitrary segment n, we solve the backward ODE P., with the corresponding terminal, using an ODE solution software.
P AP, + P.A - -
P, = 1P-(11 R)-1, with PN- = 0
(38)
0
Given the structure of R =
with c + 0 this leads us to the following expression for P,, for
t E [tn, (n+1)
(1- (n- (a (m) (o)
((aaa)))(((a) ارتو (() 217)
2(a(t)) (())) (b))
2(a(t)) af)) (t)+ 2 (c)
P
(39)
2 (t-tn=1) 1 (a 1 (() (0)))) (0)((
20 (t=(a+)((a (n) ((1h (()+ (8)
2(a(n) (t) (m)
(atn) (t) + 3( (m) (t)
25

---

where (*),, g(m) are -egment-specffic coefficients that shape the conditional bridge for the given segment n. These are recursively computed using the coefficients of the subsequent segment + 1.
as follows:
(n)) a(nt1)= (10) (bn+1) +(((t1) 2/
(0) (ln+1)
(40)
g(m) 4(a(at1) + 3(m+1))
From the expression of Pt. we can obtain its inverse
6(a 6(al")) (1)+ (t) (t- ( (1-tn+18aa(a(m)))) (8)+
-3(a
(t- (--bet))a((al)aa) (0 (1) (t)
P. x²
(41)
3(a (m) (()
3(alm) (1) (1)
(((a+1))(a( (0) 2f (()) (0)
Hence, we can compute the terms: i) gg P.-. ii) gg Pe (L, (n+1) in the optimal control formula- tion follows:
gg P
(42)
3( ((") (c))
3(a (0 i3 =6
(t- (t-bx41) [ (() (8)
(()) (())
gg PT-
(43)
3(a 2 (t) (t)+ (1))
6(a
(a (1) (0) (0)
(t))
Therefore, we can rewrite the optimal control in Eq. (17) as
x Ve
" 1
(t)
C
(n)
(t)
+
( (t)
c (t)
Fn-1, 74 E [laybnti
(44)
where we have that
((((()) (m)
23 (t) g(m) (n) (t))
(t)
(t (n+12 * (1) (t) 3(m)22 (1) (t)
3(a(") (m) (t) + 2 (n) (t)
cf" (t) (t-tx+1) (a (1) 12)
26
(45)
g(m) 22 (t)
((a + 3(m)
0 (t) (a(m) (m) (t) ((n) 22 (t)
Now we proceed to the computation of the term In+1, for which it holds from Eq. (13), that
Tn+1 = Pn+1 (P-1 1P-1+7T41 In+1+ Rm,+1)
(46)
where the term rn+1- #(tn++1 tx+2 2 n2+ carries the impact from the future segments. Therefore, the first term recursively introduces the impact of further future pinned points through
Pn+1PT11 7t11 Pat1+ Pa+1P-11 741 ((t++1:f+42raa
(47)
Since Pat1 P-1 (1,441, (n4)) is also independent of time, we have that
Pat1 P-1 n+1+ a( 2242)
(48)
(m)
where X is some static coefficient, specific for the nh segment, i.c., different segments are character- ized by different coefficients.
The structure of this matrix in Eq. (48) suggests that r2 will be dependent only on the positional constraints, when multiplied with m, for j {n + 2,...,N) .N. Finally, given the linearity of the dynamics of rt, we can recursively add the impact of more pinned points
P-11 ((tn+1, f PPt1P-1 P ( (Pn+z( P Fat2+ RRmx+2)
(49)
26

---

This recursion leads to being expressed as a linear combination of those future pinned points, through the following expression
PatiPP P n+1+ In+1+
En E=j=tz A3 F.y
(50)
Subsequently, computation of P from Pu+1 = (Pa41+ + R) of R also leads to
along with the diagonal structure
-nnt1 (((a)
Patt R mm-1
(51)
In+1
where x(()) is also some static coefficient, specific to the nh segment. Consequently, this leads to the following expression for Fn+1
-In+1 (n4m) + Ejentz 1/72 2,
Ta+1
(52)
or equivalently
-nat1 (n)
Ta+1
En E Ej=n+1 A Bz
(53)
( This implies that r, and hence the optimal control depends on all preceding pinned points after {m :h+1 2 t}. as linear combination these points. found that the elements of the vector X [An+1, 2 Ax] depend only on the number of accounted pinned points Σj. and decay exponentially as this number increases, as illustrated in
where we defined (c()
Figure 10. other words, the values of x(")) decrease the further the corresponding Ty is located from the segment whose bridge we compute.
Remark B.3. It is highlighted that the sole dependency of the coefficients (()) (() is the number of future marginals.
B.2.1 Examples of Multi-Marginaa Bridges
At this point, we provide examples of multi-margine conditional bridges, elucidating that the structure of each segment is governed by the number of future marginals. For simplicity, let us denote with (n)), g(m), (n) the segment-specifie coefficients corresponding to the 1hth segment: t [tn, fn+1).
2-marginal Bridge The formulation for a two-marginal bridge coincides with the segment N 1 for multi-marginal bridge, when N = 2, and T More specifically, we have:
3
a (mil.5g T1) =
(t - 1)2 (5r - xx) 1 1-t Uz, Vt E 10,
(54)
3-marginal Bridge The formulation of the 3-marginal bridge is given by
3 (T 1)2
a' (mee,
(i2 xz)
T-- U2,
E (t,,T)
188t 18t, - 12
12t 12t 12
a (mme(Ti1,z))
(1 - (31. 3/, 4) (Ie (Z)= +
(1 1)(31 341 4) Vt
(3t 38q 4) (52 Z1).
t t [to, t1)
(55)
Notice that for 4 1, and L/ 2, we derive the same expression as in the Example Section Additionally, we see that the segment E (o, 41) is obtained by our generalized formula for a O. 1, and coincides with the expression of the N segment. Finally, verified that the last segment shares the same formulation with the same coefficients as the bridge of the 2-marginal case.
27

---

4-marginal Bridge The 4- 4-magiinal bridge further illustrates how the structure of each segment depends on the number of future marginals. In particular, following Remark B.3, the last two segments t E (12, 7) and E [ (2) share the same formulation as in the 3-marginal bridge, since they are conditioned on and 2-future marginals, respectively. To compute the first segment, E [4o, (1), we find that a (0) 4, 3(0) 4, and the (0) vector to be: A [-1.25, 1.5, -0.25) This results in the following bridge formulation:
3
3
a" 5z)
= (T 22
T --3 Ut+
t E [t2, T)
18t 18t2 12
a" (mie[i2, Fz)
(3r x22
(t. - 122 3/2 - 4)
12t 12t₂ 12
6
U1 +
(I3 F2),
t (2)
(t t2) - 312 4)
(3t 312- - 4)
36t + 36t1 21
24t 24t1 - 211
a' (melã1, T2,
=
(t 1 ) 16|- 7)
(It æ1) + (t: ( (1)(61 16 77)
Ur
12
+ (6t - 6t1 - 7)
(00.255g
+ 1.53. 1.25511
t E [tu t1)
(56)
5-marginal Bridge It is easy to see that the last three segments of the 5-marginal bridge follow the same structure as in the 4-marginal case. For example, based on Remark B.3, the coefficients for the third-to-lee lasst segment, E f2), are @1 = 31 4, and the vectors X2)) and A()), corresponding to the segments E [tz, ta) and E [t1. f2). respectively, match those of the third-to-last and second-to- last segments in the 4-marginal bridge. Finally, for the first segment, / E [0, 41,, we compute that: (0) 28, g(0) 32, and (0) [-1.667, 1.6, -0.4, 0.067]T. Substituting these coefficients into Eq. (45) yields the corresponding bridge formulation.
3 (T' t)
3
a" (m/) R1)
(I1 I²t)
Ut
te [t,,T) T)
T
188 18// -
a" a'(me Rs) =
(I) - Tz)
(1 - (2)) (31- 313- 4)
12t 1213 12
+ (L - 131 - 3/3 - 4)
Ut +
(34 Iz),
L ((2,la)
31 313 4
366 + 3612 - 21
241 - 2412 - 21
a" a'(me T2, 13. F4)
(0f 52) + (t tz2 t2) -61 7)
(t 12) (6t 6tz
12
+ (6t 61 7)
(-0.25%,
+ 1.533 - 1.25582)
E (2)
6 (45t 45t 26)
12 (15t 15t, 13)
a" amel 181, 23,
=
(t 1 (45t - 45t - 52)
(x) - 51) +
(t - t1) (45t - 45t1 52)) 116
90
45t 45 52 ( (-1.2677)
+ +1.3g - 0.43g + 0.0677,),
t Le (1:4)
(57)
B.3 Proof of Proposition 3.4
Proposition B.4. Let define the marginal path Pt mixture of bridges pe(me)
(Pt (mm : In 2 1) (ddg((x)))) where Pota, y (me) (ān : In 2 13) the conditional proba-
bility path associated with the solution of the 3MBB path in Eq 6. The parameterized acceleration
28

---

that satisfies the FPE prescribed by the Pt is given by
(t m+) J [a PP (an (a
(milkan - : in 2 #))dg((Zn))
(58)
Pt
This suggests that the minimization of the variational gap to match of given pe is given by
min E ) Eeucin 1 [aa](s.)) d
(59)
Proof. We want to show that the acceleration from Eq. preserves the prescribed path pt. The
momentum Fokker Plank Equation (FPE) is given by
1 (.p:(m,) = ty (m.) v.(a((())p((mm))) + 00² (m
(60)
We let the marginal be constructed mixture of conditional probability paths conditioned collection of pinned points {Ia Jne (= (ddnn)))ddd )ddddx dæn. Using this definition for the marginal path, obtains that
/ (((ta))0apmmmm @epe(me) = / ((
[VDe((m) 04 V. / ((( (ut dFr,)
Appe(m:)= = V. (V.pi(me)) = - -V. p
= a((()n)) vo (V. p((mullantt) dfan))
/ ((p
(61)
Hence remains be checked whether the following equality holds
az (p = q( ((((x))) (agp (maa (Tx (mm)) }) d(a,}
(62)
which suggests that the parameterized drift that minimizes the following minimization problem of (me) arg min He [Ilaa azja_ (m. (Tx
(63)
preserves the prescribed pt.
B.4 Proof of Theorem 3.5
Definition B.5 (Markovian Projection of path measure). The Markovian Projection of TP is defined as p.M. arg minyeM (T(PP....
Intuitively, the Markovian Projection seeks the path measure that minimizes the variational distance to P. In other words, it seeks the closest Markovian path measure the KI. sense.
Definition B.6 (Reciprocal Class and Projection). For multi-margine path measures, we say that is in the reciprocal class (Q) of Q E .M if
f ((nn)d(((2z))
PP =
namely, it shares the same bridges with Q. We define the reciprocal projection of TP as
projæ (P) = arg min TER(
KL(P|]
Similarly, the Reciprocal Projection yields the closest reciprocal path measure in the KL sense.
29

---

Lemma B.7. Let P be a Markov measure in the reciprocal class of E M such that J F,dd qn (. for n [O....,N}. Then, under some mild regularity assumptions on Q. 9n. it found that P is equal to the unique multi-marginat the Schrödinger Bridge TpmmS8
Proof. First let us assume that KL(PIQ) < oo, and that KL(gal f Qn < 0o for n E Assume Q E M, then by (Theorem 2.10, Theorem 2.12 Léonard [2013]). it follows that the solution of the dynamic SB IP must also be Markov measure. Finally, from the factorization of the KL, it holds that
(2n) (Zn) Q(xx))++ / /K..(P (T.((QQ) (t.
(64)
which implies that (T(( (z,)) (e S (T,(P) Q) with equality (when KT((P (Q) 20) if and
only if PI Piz.. Q )= Therefore, [P* is the (unique) solution mmSB if and only if disintegrates
above (Proposition 2.3 Léonard [2013]).
Lemma B.8. (Proposition in [Shi et al., 2023]] Let V M and T R(Q) and IFKL(PIV) co. and if (PIV) < oo we have
KL(PIV) = (L(P) (P)) + KL(proi
(65)
and if (P(T)) < oc, then
KL(PIT)
projR(e) (P)) + (L(proj (P))T)
(66)
Theorem B.9. Assume that the conditions of Lemma B.7, and B.8 hold. Then, our iterative scheme admits a fixed point solution IP*, i.e., KL(PP1PP + O, and in particular, this fixed point coincides with the unique pmmmS
Proof. We define the following path measures V projM( and projR(o)) (V). Assume the conditions for Lemma B.8 hold for IP, V and T. Then for any arbitrary fixed point V we can write
KL KL(VVVV) + KL(V
Σ KLLVVVV + KL(V 1=
(67)
Thus, it holds that KL KL(V'VV) os for every iteration i e Similarly, for P., and T, we obtain KT((P(IP') < xo and KT. oc for each E N for any arbitrary fixed IP' and TV.
Consequently, define the following function
V
KL V' KL(T PL
(68)
For two consecutive iterates i and i + 1, we have
:= KL(VVV)) (L(T)T) KL(PIIP')
gi++ KL(V++I + KL(T++ (T') + KL(Fi(P)
Using Lemma B.8, we can rewrite v as
Ve == KL(V' V') + KL(T (L | + KL(Vi+I KL + 1
= -KL(V*
KL
= -LLVV + KL(T'(T+t)) KL p + gi+1
Now, we take the sum of this telescoping series and obtain
KL(V + KL(T)T+!) + (L(
(69)
Note that To and The are finite (with To > (T0), since KL(PP(P' co, KL(VIV) < 00 and KL < so for every iteration i E N. Therefore, since we also have (T.((P)((+t) ( v and KL(T(T-+) > 0, get that
30

---

lim KL pi+1) ) 0,
lim KL(V*/V++) + 0,
lim KL(T 0.
Hence the iterates ps, V and T converge some fixed points Pd P V V*, and T + T
as + 0C.
From the factorization of the KL divergence, we have for consecutive projections of our algorithm:
KL(TOPPT)) (F.))
Bygo (nn)
KL
(a., (em)))
(z..
(70)
= KL (zx (Z..)))
since the bridges of T are the same with P. and Q, i.e. T(()
Pr Ment ((n) Then,
[(n))
KL(V KL (x.) P(an)) Byy
KL(V(()) |P
(an) (am)))
(on)
(71)
KL ( I (2n))
(ax)
since the coupling after the Markovian projection at iteration i remains the same for the reciprocal path measure at iteration + 1, namely y( = T6+1) Therefore, we can deduce {ax)
(an)
KL((VIP') (T(
(72)
We further assume that KL(TTIP*) oo, KL(VIPF* 0o. Therefore, the iterations of Eq. 70 and 71, yield KL(T (P) > KL(VIPP) 2 KL(T++I(P) for 2 0, implying con- vergence, since it is non-increasing and bounded below. Applying Lemma B.8, we obtain lim, (KL(TT) [1+* KL(V*(P*)
lim; KL(T'IV)) O. By definition of the lower semi- continuity of the KI. divergence, we have KI((V*(T*) S lim inf KL(V (T>>)), Additionally, by the definition of the KI. divergence, we have 0 < (L((VTTT). Finally, it also holds that
lim inf
( (
O. Combining all three claims have
< (L(V")T*) < inf (L [T*>* = O.
k-+00
Therefore, V T which also means that MMR(()) [Shi et. 2023]. Also, by construction, all the iterates of ui satisfy the positional marginal constraints PP,, On. hence also 1Pt de, Go Therefore, by Lemma B.7, IP* is the unique multi-marginal Schrödinger bridge pomSB
C Gaussian Path
The scalability of our framework is based on the capacity to efficiently sample from the conditional gaussian path induced by the solution of of the 3MBB path in Eq 6. Let us define the marginal path Pt as a mixture of bridges Pe (m,) = Ptifz., (m, t, 2 t} {))(((zz,)) conditional probability path associated with the solution of the 3MBB path in Eq The linearity of the system implies that we can efficiently sample m, [ vz]. VL E [0, T]. from the conditional probability path Potaa, = N(,,,, E,), as the mean vector μr and the covariance matrix
where Pa Pa( (me [Tn
in 2 (})
Σ, have analytic solutions [Särkkä and Solin, 2019].
Let us recall the optimal control formulation for the th segment.
xx Vz
0
I
Int1 1 E A5āz.
u = (t))
(aa)
(t)
+
((
ch (l)
Vt E [tn, tx+1] (73)
31

---

Thus, the augmented SDE for the corresponding is written
dm, -Iopo0
1
X Vi +
- (1) + (c (1)) E 1.m 1,
1) dt + gdW;
C²
(t)
(74)
where C (t), 20 (t), (t) are the segment specific functions, and X
(m) the coefficients for the future pinned points. To find the expressions the mean and the covariance, we follow [Särkkä and Solin, 2019], and consider the following ODEs for pt, and Eg respectively:
1
A = (()
(1)
(1)
+ ct) - ct (t) Lj=n+1 A/
(75)
0
Σ (t)
(t)
(t))
E2 (1) (t) _rJT gg
Mean ODEs
If we explicitly write the mean ODE system, we obtain the following two ODEs
Ax 4
j = -1 (1)pxz + C22 (t))o -C Ct" (t) + cl
Σ
(76)
j=++1
which corresponds the following second-order ODE
ju c.a² (t) (t) G. (t) c (f) + cl
Σ
(77)
j=++1
This ODE is then solved using an ODE solver software for the corresponding functions (()) of the respective segment n: t [t,, 1).
Covariance ODEs If we explicitly write the mean ODE system, we obtain the following two ODEs
Ege = 2Σgv
Σnv clah + Ch"h (t)2 + Evv
(78)
S.D -2 (1)E.. + 2c (4)E +
which corresponds to the following third-order ODE
2 5: 2 C CL" (1)E (1)"2cc (3=" ((m) (1) (1)) (0- (1) (t)
(1)
(79)
This equation, however, is hard to solve even using software packages. For this reason, we integrate the covariance ODEs using Euler integration, once at the beginning of our training. This procedure can be solved once and can be applied for any fixed coupling, during the matching, since the system of ODEs in Eq. (78) does not depend on any points m, but its sole dependence on time. This suggests that the computational overhead from simulating the covariance ODEs is negligible.
Tue [
Iea
Eev Duv
Given the expressions of μr
and N, =
one can obtain m, through
Dev
X, V1
L ED en
m.
+
L
L
(80)
L L 6p
I L
where the matrix L.
is computed following the Cholesky decomposition to the
covariance matrix, and
(0 Iza)
61
32

---

D Extended Related Works
Schrödinger Bridge Recently, in generative modeling there has been a surge of principled ap- proaches that stem from Optimal Transport [Villani et al., 2009]. The most prominent problem formulation has been the Schrödinger Bridge (SB; Schrödinger [1931]). In particular, SB gained sig- nificant popularity in the realm of generative modeling following advancements proposing a training scheme based on the Iterative Proportional Fitting (IPF), a continuous state space extension of the Sinkhorn algorithm to solve the dynamic SB problem [De Bortoli et al., 2021, Vargas et al., 2021]. Notably, SB generalizes standard diffusion models transporting data between arbitrary distributions To, Ti with fully nonlinear stochastic processes, seeking the unique path measure that minimizes the kinetic energy. The Schrödinger Bridge [Schrödinger, 1931] in the path measure sense is concerned with finding the optimal measure P that minimizes the following optimization problem
min (P(())Q) Po
no, Pi T
(81)
where Q is a Markovian reference measure. Hence the solution of the dynamic SB IP* is considered to be the closest path measure to Q. Another formulation of the dynamic SB crucially emerges by applying the Girsanov theorem framing the problem as Stochastic Optimal Control (SOC) Problem [Chen et al., 2016, 2021].
min /'' Mt.Pt
Ep, [| |dt s.t. Öp: Ōt
o² 2
= - (uppr) + Apt.
and
Po
To
P1
T1
(82)
Finally, note that the static SB is equivalent to the entropy regularized OT formulation [Pavon et al., 2021, Nutz, 2021, Cuturi, 2013].
min E ()o. Rdd
-0 - æ1 1| 1 (x0, x1) + <KL(T| @ x1)
(83)
This regularization term enabled efficient solution through the Sinkhorn algorithm and has presented numerous benefits, such as smoothness, and other statistical properties [Ghosal et al., 2022, Léger, 2021, Peyré et al., 2019].
Bridge Matching Peluchetti [2023] first proposed the Markovian projection to propose Bridge matching, while Liu et al. [2022a] employed it to leam representations in constrained domains. The Bridge matching objective offers a computationally efficient alternative but requires additional assumptions. To this front, Action Matching Neklyudov et al. [2023] presents a general matching method with the least assumptions, at the expense of being unfavorable to scalability. Additionally, recent advances have introduced more general frameworks for conditional generative modeling. Denoising Diffusion Bridge Models (DDBMs) extend traditional diffusion models to handle arbitrary source and target distributions by learning the score of a diffusion bridge, thereby unifying and generalizing methods such as score-based diffusion and flow matching [Zhou et al., 2023]. Similarly, the stochastic interpolant framework [Albergo et al., 2023] integrates flow- and diffusion- based approaches by defining continuous-timg stochastic processes that interpolate between distributions. These interpolants achieve exact bridging in finite time by introducing an auxiliary latent variable, offering flexible control over the interpolation path.
Recently, these matching frameworks have been employed solve the SB problem. DSBM ([Shi et al., 2023]) employs Iterative Markovian Fitting (IMF) to obtain the Schrodinger Bridge solution, while De Bortoli al. [2023] explores flow and bridge matching processes, proposing a modification to preserve coupling information, demonstrating efficiency in learning mixtures of image translation tasks. SF22 M [Tong et al., 2023a] provides a simulation-free objective for inferring stochastic dynamics, demonstrating efficiency in solving Schrödinger bridge problems. GSBM Liu et al. [2024] presents framework for solving distribution matching to account for task specific state costs. While these methods aim to identify the optimal coupling, [Somnath et al., 2023] and [Liu et al., 2023] propose Bridge Matching algorithms between priori coupled data, namely the pairing between clean and corrupted images or pairs of biological data from the static Schrödinger Bridge. Lastly, works that aim to improve the efficiency of these matching frameworks by introducing a light solver for implementing optimal matching using Gaussian mixture parameterization [Gushchin al., 2024a].
33

---

Flow Matching In parallel, there have been methodologies that employ deterministic dynamics. Lipman al. [2022] introduces the deterministic counterpart of Bridge Matching; Flow Matching (FM) for training Continuous Normalizing Flows (CNFs;Chenn et al. [2018) 8)) using fixed conditional probability paths. Further developments include Conditional Flow Matching (CHM) offers stable regression objective for training CNFs without requiring Gaussian source distributions or density evaluations [Tong et al., 2023b]. Metric Flow Matching (MFM), which learns approximate geodesics on data manifolds [Kapusniak et al., 2024], and Flow Matching in Latent Space, which improves computational efficiency for high-resolution image synthesis [Dao et al., 2023]. Finally, CFM retrieves exactly the first iteration of the Rectified Flow Liu al. [2022b]. which is an iterative approach for learning ODE models to transport between distributions.
Multi-Marginal Among the advancements of Flow matching models was the introduction of the Multi-Marginal Flow matching framework [Rohbeck al., 2024]. Similarly to our approach, they proposed a simulation-free training approach, leverages cubic spline-based flow interpolation and classifier-free guidance across time and conditions. TrajectoryNet [Tong et al., 2020] and MIOFlow (Huguet et al., 2022] combine Optimal Transport with Continuous Normalizing Flows [Chen et al., 2018] and manifold embeddings, respectively, to model non-linear continuous trajectories through multiple points. In the stochastic realm, recent works have proposed the mmSB training through extending Iterative Proportional Fitting a continuous extension the Sinkhorn algorithm to solve the dynamic SB problem [De Bortoli et al., 2021]- to phase space and adapting the Bregman iterations [Chen et al., 2023a]. Another approach alternates between learning piccewise SBs on the unobserved trajectories and refining the best guess for the dynamics within the specified reference class [Shen et al., 2024]. More recently, modeling the reference dynamics as a special class of smooth Gaussian paths was shown achieve more regular and interpretable trajectories [Hong et al., 2025]. Furthermore, the multi-marginal problem has been recently addressed by Deep Momentum Multi-Marginal Schrödinger Bridge [DMSB;Chen et al. [2023a]] proposed to solve the mmSB in phase space via adapting the Bregman iterations. More recently, an iterative method for solving the mmSB proposed learning piecewise SB dynamics within a preselected reference class [Shen et al., 2024]. In contrast, modeling the reference dynamics as smooth Gaussian paths was shown to achieve more temporally coherent and smooth trajectories [Hong et al., 2025], though the belief propagation prohibits scaling in high dimensions. Lastly, Wasserstein Lane-Riesenfeld (WLR) is a geometry- aware method to reconstruct smooth trajectories from point clouds. performs consecutive geodesic averaging in Wasserstein space, giving spline-like curves that can handle mass spliting/bifurcaations. achieving strong results on cell datasets [Banerjee et al., 2024].
E
Additional Details on Experiments
E.1
General Information
In this section, we revisit our experimental results to evaluate the performance of 3MSBM on variety of trajectory inference tasks, such as Lotka-Volterra, occan current in the Gulf of Mexico, single-cell sequencing. and the Beijing air quality data. We compared against statc-of-the-ar methods explicitly designed to incorporate multi-marginal settings, such Deep Momentum Multi Marginla Schrödinger Bridge (DMSB;Chen Refinement (SBIRR; Shen et al. [2024]), smooth Schrodinger Bridges (SmoothSB: Hong et al. [2025]), and Multi-Marginal Flow Matching (MMFM; Rohbeck et al. [2024]), along with two NeuralODE- methods: MIOFlow [Huguet et al., 2022] and DeepRUOT [Zhang et al., 2024b). We used the official implementations of all compared methods, with default hyperparameters. For all experiments with our 3MSBM, we employed the ResNet architectures from Chen et al. [2023b], Dockhorn et al. [2021]. We used the AdamW optimizer and applied Exponential Moving Averaging with a decay rate of 0.999. All results are averaged over 5 random seeds, with means and standard deviations reported in Section 4 and the tables below. Experiments were run on an RTX 4090 GPU with 24 GB of VRAM.
al. [2023a]), Schrodinger Bridge with Iterative Reference
F.2 Volteraa
We first consider synthetic dataset generated by the Lotka-Volterra (LV) equations [Goel et al., 1971], which model predator-pres interactions through coupled nonlinear dynamics. We used the dataset from [Shen et al., 2024] with the 5 training and 4 validation time points, with 50 observations
34

---

Table 9: Mean of SWD, MMD, W1, and W3 over seeds at left-out marginals the LV dataset for DeepRUOT, MIOFlow, SBIRR, MMFM, Smooth SB, and 3MSBM (Lower is better).
SWD
(b) MMD
Method
t3
ls
l7
Method
f3
ts
l7
DeepRUOT
0.30
0.16
0.22
0.44
DeepRUOT
3.02
0.20
0.39
0.48
SBIRR
0.17
0.18
0.24
0.48
SBIRR
0.44
0.53
0.46
0.50
MIOFlow
1.53
1.49
1.23
1.46
MIOFlow
7.09
6.52
6.29
5.28
Smooth SB
0.49
0.18
0.13
0.37
Smooth SB
3.23
1.21
0.32
0.42
DMSB
0.64
0.67
0.98
0.63
DMSB
3.46
3.43
1.16
1.74
MMFM
0.21
0.25
0.39
0.57
MMFM
0.52
0.60
0.63
0.77
3MSBM
0.29
0.13
0.11
0.37
3MSBM
4.33
0.72
0.40
0.40
(c) W1
(d) W2
Method
f1
f3
ls
t7
Method
t1
f3
ts
17
DeepRUOT
0.40
0.17
0.29
0.62
DeepRUOT
0.44
0.19
0.31
0.65
MIOFlow
2.53
2.12
1.76
1.53
MIOFlow
2.77
2.34
1.76
1.55
SBIRR
0.19
0.20
0.30
0.48
SBIRR
0.19
0.25
0.45
0.74
Smooth SB
0.29
0.27
0.22
0.68
Smooth SB
0.27
0.27
0.22
0.68
DMSB
0.99
0.74
0.60
0.98
DMSB
0.84
0.98
0.71
1.24
MMFM
0.24
0.43
0.57
0.75
MMFM
0.22
0.43
0.77
1.23
3MSBM
0.23
0.18
0.12
0.35
3MSBM
0.24
0.17
0.15
0.36
per time point. In particular, the generated dataset consists of 9 marginals in total; the even-numbered indices are used to train the model (i.e., to, f2, 14, fG, is), and the remainder of the time points are used assess the efficacy of our model impute and infer the missing time points. this experiment, we benchmarked 3MSBM against DeepRUOT, MIOFlow, DMSB, SBIRR, MMFM, and Smooth SB. Table 9 reports the mean performance over 5 seeds of each method with respect the SWD, MMD, W1, and W2 distances from the validation points. The hyperparameter selection for the LV with our method were: the diffusion coefficient was set to a = 0.3, and the learning rate was 10-1
E.3 Gulf of Mexico
Subsequently, we evaluate the efficacy of our model to infer the missing time points in real-world multi-marginal dataset. The dataset contains ocean-current snapshots of the velocity field around a vortex in the Gulf of Mexico (GoM). Similarly to the LV dataset, we used the big vortex dataset provided in [Shen et al., 2024], consisting of 300 samples across training and validation times. More explicitly, out of the total 9 marginals, the even-indes even-indexed time points (i.e., to, t2, tt, tG, ts) are used for training, and the remaining are left out evaluate the model's ability impute and infer missing temporal states. We compared our 3MSBM against DeepRUOT, MIOFlow, DMSB, SBIRR, MMFM, and Smooth SB for this experiment. Table 10a demonstrates the mean performance over seeds of each method with respect to the SWD, MMD, Wi, and W22 distances from the validation points. The hyperparameters used for the GoM experiment with our method were: batch size of 32 for the matching, the diffusion coefficient was set 0.3, and the learning rate was set equal to 10
E,4 Beijing air quality
We revisit our experiments using the Beijing multi-site air quality data set [Chen, 2017]. This dataset consists of hourly air pollutant data from 12 air-quality monitoring sites across Beijing. We focus on PM2.5 data, an indicator monitoring the density of particles smaller than 2.5 micrometers. from January 2013 to January 2015. We focused on single monitoring site and aggregated the measure- ments collected within the same month. To introduce temporal separation between observations, we selected measurements from every other month, resulting in 13 temporal snapshots. For imputation task, we omitted the data at /,, ls. (s. and /11, while the remaining snapshots formed the training set. Table la shows the mean performance over 5 seeds of each method in the SWD, MMD, Wi, and W2 distances from the validation time points, benchmarking 3MSBM method against
35

---

Table 10: Mean of SWD, MMD, Wi. and W/2 over 5 seeds at left-out marginals the GoM dataset for DeepRUOT, MIOFlow, SBIRR, MMFM, Smooth SB, and 3MSBM (Lower is better).
SWD
(b) MMD
Method
t3
17
Method
f3
ts
l7
DeepRUOT
0.21
0.33
0.23
0.21
DeepRUOT
2.35
1.75
1.30
1.10
MIOFlow
0.83
0.34
1.23
0.97
MIOFlow
6.52
4.13
6.88
6.11
SBIRR
0.15
0.11
0.11
0.09
SBIRR
4.65
0.82
1.15
0.35
Smooth SB
0.17
0.14
0.10
0.13
Smooth SB
4.98
4.02
0.84
0.28
DMSB
0.23
0.54
0.39
0.28
DMSB
4.50
4.40
4.57
3.37
MMFM
0.23
0.25
0.10
0.19
MMFM
0.91
1.20
0.49
0.51
3MSBM
0.14
0.14
0.08
0.06
3MSBM
3.99
2.92
0.87
0.52
(c) W1
(d) W2
Method
f
f3
15
t7
Method
f3
fs
17
DeepRUOT
0.29
0.29
0.20
0.25
DeepRUOT
0.32
0.43
0.36
0.33
MIOFlow
1.10
0.55
1.67
1.69
MIOFlow
1.11
0.48
1.68
1.70
SBIRR
0.28
0.15
0.11
0.15
SBIRR
0.24
0.13
0.21
0.17
Smooth SB
0.16
0.27
0.21
0.56
Smooth SB
0.22
0.27
0.21
0.16
DMSB
0.24
0.31
0.70
0.46
DMSB
0.28
0.30
0.72
0.46
MMFM
0.33
0.38
0.22
0.29
MMFM
0.33
0.32
0.19
0.31
3MSBM 0.17
0.21
0.09
0.12
3MSBM
0.20
0.18
0.07
0.09
Table 11: Mean of SWD, MMD, W1, and W22 over seeds at left-out marginals the Beijing Air Quality dataset for DeepRUOT, MIOFlow, MMFM, Smooth SB, and 3MSBM (Lower is better).
(a) SWD
(b) MMD
Method
f
f3
f5
f7
Method
f
f3
fs
f7
DeepRUOT
13.67
52.60
71.34
84.67
DeepRUOT
0.36
0.34
0.99
1.67
MIOFlow
46.64
79.06
76.06
60.87
MIOFlow
0.58
0.92
0.38
0.51
Smooth SB
28.61
28.81
35.79
32.90
Smooth SB
0.41
0.43
0.39
0.48
DMSB
21.10
21.92
35.53
35.75
DMSB
0.76
0.79
0.54
0.47
MMFM
17.51
23.94
32.56
39.98
MMFM
0.44
0.56
0.59
0.55
3MSBM
17.70
9.78
22.23
32.23
3MSBM
0.35
0.85
0.28
0.32
(c) Wi
(d) V/
Method
f
f3
f
Method
t
f3
f6
DeepRUOT
10.49
39.82
51.00
68.97
DeepRUOT
13.67
52.59
71.35
84.67
MIOFlow
31.79
56.35
57.89
45.36
MIOFlow
46.64
79.06
76.06
60.87
Smooth SB
23.73
23.88
24.70
28.61
Smooth SB
28.61
28.81
35.79
32.90
DMSB
58.79
32.70
40.22
42.06
DMSB
60.19
38.77
41.38
43.25
MMFM
28.08
26.40
37.73
51.12
MMFM
26.42
29.95
43.49
49.96
3MSBM
12.71
57.44
26.02
29.61
3MSBM
12.87
79.36
27.89
32.26
MMFM with cubic splines, DeepRUOT, MIOFlow, and Smooth SB. Note, for this experiment, did not benchmark against SBIRR, since we did not possess the corresponding informative prior measure. The hyperparameters used for the Beijing air quality experiment with method were: total number of samples of 1000 were used, with a batch size 64 for the matching, the diffusion coefficient was
0.2, and the learning rate was set 5 10-5
F.5 Single sequencing
Lastly, we revisit our experiments on the Embryoid Body (EB) stem cell differentiation dataset, which captures cell progression across stages over 27-day period. Following the setup in Section 4.4, used the preprocessed data from [Tong et al., 2020, Moon et al., 2019], embedded
36

---

in 100-dimensional PCA feature space. Cell snapshots were collected at five discrete intervals: to E [0.3]. i [6. 9]. tz E [12, 15). t3 [18. 21), t4 E [24, 27]. Below, we present the results of the comparison of our 3MSBM against DeepRUOT, MIOFlow, DMSB, SBIRR, MMFM. Table 10a demonstrates the mean performance over 5 seeds of each method with respect to the SWD, MMD, W1, and W2 distances from the validation points, i.e., at the snapshots 4,, and 13- Observing the results in Table 12 is evident that our 3MSBM consistently outperforms the state-of-the-ar algorithms in the high-dimensional EB-100 task across all metrics. The hyperparameters used for every EB experiment with our method were: a total number samples 1000 were used, with a batch size of 64 for the matching, the diffusion coefficient was set to 0.1, and the learning rate was set to 10
Table 12: Mean of SWD, MMD, W1, and W/ over seeds at left-out marginals on the EB-100 data for DeepRUOT, MIOFlow, SBIRR, MMFM, Smooth SB, and 3MSBM (Lower is better).
(a) SWD
(b) MMD
Method
t
f3
Method
f
f3
DeepRUOT
0.73
0.67
DeepRUOT
0.43
0.36
MIOFlow
0.84
0.94
MIOFlow
1.01
0.92
SBIRR
0.80
0.91
SBIRR
0.71
0.73
DMSB
0.58
0.54
DMSB
0.38
0.36
MMFM
0.59
0.76
MMFM
0.37
0.35
3MSBM
0.48
0.38
3MSBM
0.18
0.14
(c) Wi
(d) 14
Method
4
t3
Method
13
DeepRUOT
13.45
14.90
DeepRUOT
13.64
15.10
MIOFlow
13.20
13.57
MIOFlow
13.66
14.05
SBIRR
15.09
20.39
SBIRR
15.42
20.98
DMSB
14.08
15.22
DMSB
14.83
15.49
MMFM
13.61
14.64
MMFM
14.68
14.83
3MSBM
13.89
13.11
3MSBM
14.51
13.26
E.6
Ablation study on
In stochastic optimal control [Theodorou et al., 2010], the value a plays a crucial role in represent- ing the uncertainty from the environment the error in applying the control. As result, the optimal control policy can vary significantly with different degrees of noise. Figure 14 demonstrates the performance of our 3MSBM with respect varying noise in the EB and GoM datasets. We observe consistent performance in the training marginals across all tested values of a, whereas for the points in the validation set, increasing up point is deemed beneficial Sample trajectories in GoM in Figure 13 further verify this trend. Low noise values (e.g. = 0.05) cause the trajectories to be overly tight, whereas high noise (e.g. 1.0), the trajectories become overly diffuse. On the other hand, moderate noise values (e.g. a = 0.4) achieve a good balance, enabling well-spread trajectories matching the validation marginals.
it improves performance.
o 0.05
o = 0.4
o 1.0
Time Marginal
Figure 13: Comparison of the trajectories inferred on the Gulf Mexico current dataset for different
values of
37

---

EB
GoM
0.4
0.3
0.3
SWD
Validation Training
O0
0.2
Validation
0.2
Training
0.1
0.1
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
0
0
Figure 14: SWD from the marginals in the Validation and Training set for varying values of sigma on EB and GoM
38


Reviews:
[{"review_id": "32k4oyt0wJ", "reviewer": "Reviewer_dYDf", "summary": "This paper introduces Momentum Multi-Marginal Schrödinger Bridge Matching (3MSBM), a novel framework for inferring smooth, temporally coherent stochastic trajectories that simultaneously satisfy multiple positional constraints. By lifting dynamics into phase space (position + velocity), the authors extend classical Schrödinger Bridge theory to handle an arbitrary sequence of intermediate distributions.", "questions": "1. There are notable inconsistencies between expressions (2) and (3), (4). Specifically, (2) involves the term $u_tdt$, while (4) only contains \n$u_t$; in addition, (2) does not include 's.t.', and (4) lacks the definitions of $A$ and $g$.\n\n2. Is there any condition imposed on $u_t$? In particular, is it required that the first half of $u_t$ be zero?\n\n3. Essentially, is (2) a second-order equation? If so, why is it expressed in the current form?\n\n4. The ODE version, i.e., dynamic optimal transport, corresponds exactly to optimal transport. Does the method proposed in the paper admit a similar interpretation in terms of optimal transport?\n\n5. It is claimed that existing methods struggle to capture long-range temporal dependencies, which affects the coherence of the inferred trajectories. However, it remains unclear why formulating the problem using a second-order SDE addresses this specific limitation.", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Strengths\n1. A simulation-free, matching-based method that preserves all intermediate marginals throughout training, avoiding trajectory caching and time-discretization issues.\n\n2. Generalization of momentum Brownian bridges to the multi-marginal setting, with analytic expressions for the optimal control and conditional acceleration.\n\n3. Provable stable convergence and efficient implementation that scale to high dimensions and many marginals.\n\n4. Demonstrated superior performance over state-of-the-art multi-marginal methods (DMSB, SBIRR, Smooth SB, MMFM) on synthetic dynamics (Lotka–Volterra), real-world flows (Gulf of Mexico, Beijing air quality), and high-dimensional single-cell differentiation data.\n\nWeaknesses\n\nSeveral parts of the manuscript suffer from unclear wording and a high density of technical or grammatical errors, which hinder readability and understanding.", "was_revised": true, "final_justification": "The response has successfully resolved all the issues I had.", "rebuttal": "We would like to thank the reviewer for their positive comments, as well as their constructive criticism and further suggestions for improving our work. In the following, we address the points raised by the reviewer.\n\n\n### **1. Typo Corrections**\n\nWe appreciate the reviewer’s suggestions to improve the clarity and readability of our work. We have carefully reviewed the paper and made several corrections that will be incorporated into the revised version, as we can not update the manuscript during the rebuttal phase. More specifically, we have extensively addressed and corrected technical typos in Equations, such as: \n\n- Eq. (4) was missing a differential $dt$, hence the corrected $dm_t = Am_t dt+udt + \\mathbf{g}dW_t$.\n- Eq. (2) was corrected to include the ‘s.t.’ to explicitly denote the constraints of our optimization problem in Eq. (2).\n- The definitions of $A$ and $\\mathbf{g}$ in Eq. (4) follow from Eq. (2). However, we included them in Eq. (4) as well to improve readability.\n\nAdditionally, we thoroughly reviewed the manuscript and corrected grammatical errors. These include typographical errors (“phase state” → “phase space”), punctuation/parallel‑structure fixes, article misuse and noun phrase issues (e.g., “globally optimality coupling” → “globally optimal coupling”), and bracket mismatches in Algorithm 1. \n\nWe are confident these improvements have greatly enhanced readability and have addressed the clarity concerns raised.\n\n### **2. Notation clarifications**\n\nOur current notation follows prior multi-marginal and momentum works [1,3]. Our system is indeed second-order, consisting of the two following differential equations\n\n$$dx_t = v_tdt $$\n$$dv_t = a_tdt + \\sigma dW_t$$\n\nwhere $x_t$  is the state, $v_t$  is the velocity, $a_t$  is the acceleration, $W_t$  is the Brownian motion and $\\sigma$  is a diffusion coefficient for the noise being injected only through the velocity channel.  Thus, the two differential equations can be rewritten in a more compact form as \n\n$$\n\\begin{bmatrix}dx_t\\\\\\ dv_t\\end{bmatrix} = \\bigg(\\begin{bmatrix}0&1\\\\\\0&0\\end{bmatrix}\\begin{bmatrix}x_t\\\\\\ v_t\\end{bmatrix} + \\begin{bmatrix}0\\\\\\ a_t\\end{bmatrix} \\bigg)dt + \\begin{bmatrix}0&0\\\\\\0&\\sigma\\end{bmatrix}dW_t\n$$\nBy defining the augmented state $m_t =\\begin{bmatrix}x_t\\\\\\ v_t\\end{bmatrix}$, $A=\\begin{bmatrix}0&1\\\\\\0&0\\end{bmatrix}$, $u_t = \\begin{bmatrix}0\\\\\\a_t\\end{bmatrix}$, and $\\mathbf{g} = \\begin{bmatrix}0&0\\\\\\0&\\sigma\\end{bmatrix}$, we can write the augmented SDE above in a more concrete fashion, as follows:\n\n$$\ndm_t = (Am_t + u_t)dt + \\mathbf{g}dW_t\n$$\n\nTherefore, the control $u_t$  is not constrained in any manner, rather than it is by definition that the control is injected only through the velocity channel through the acceleration, and hence follows the aforementioned formulation.\n\nTo avoid confusion, we have revised Eq. (2) and (3) & (4) to enhance readability. In particular, we now explicitly include the differential equations as part of the constraints and have restructured the optimization objectives for better presentation. More explicitly, Eq. (2) is now expressed as:\n$$\n\\min_{a_t} \\int_0^1 \\mathbb{E_\\text{$p_t$}} \\Bigr[\\frac{1}{2}||a_t||^2\\Bigr ] dt,\n$$\n$$\n\\text{s.t.}\\quad dx_t =v_t dt, \\\\quad v_t= a_tdt+\\sigma dW_t, \\\\quad x_n\\sim q_{n} = \\int \\pi_n(x, v) dv_n, \\\\quad n= \\\\{ 0, 1, \\dots, N \\\\}\n$$\nand Eqs. (3) and (4) are rewritten as follows:\n$$\n\\min_{a_t} \\int_0^1  \\frac{1}{2}||a_t||^2 dt + \\sum_{n=1}^N (m_{n} - \\bar{m}_{n})^\\intercal R (m_n - \\bar m_n)\n$$\n$$\n\\text{s.t.}\\quad  dx_t =v_t dt , \\quad dv_t= a_tdt+\\sigma dW_t , \\quad m_0 = \\bar{m}_0 .\n$$\n\n### **3. Optimal Transport Interpretation**\n\nIn this work, we focus on stochastic dynamics to accurately model the inherent randomness present in many real-world systems, such as those in meteorology and biology. The multi-marginal Schrödinger Bridge (mmSB) problem [2] has been shown to be the stochastic counterpart of optimizing measured value splines in the Wasserstein space [4], which can be interpreted as the multi-marginal Optimal Transport (mmOT) problem [5]. Notably, our 3MSBM provably converges to the mmSB solution — i.e., the joint optimal coupling over the set of $N$  marginals. Therefore, our method provably solves the stochastic equivalent of the momentum formulation of the mmOT problem.\n\n### **4. Motivation for Second-Order Dynamics**\n\nTo capture coherent long-term dependencies in sparse and irregularly sampled time series, we incorporate second-order dynamics into our dynamical systems and optimize them to minimize the acceleration of the path. This results in minimum curvature paths leading to smoother and consistent trajectories. Crucially, the inclusion of velocity into the dynamics provides information to the model about past trends—through the acceleration—that is otherwise lost in first-order systems. In contrast, first-order systems are characterized by memoryless dynamics, which lack the capacity to preserve momentum trends, often resulting in interpolated trajectories with unnatural *kinks,* as illustrated in Figure 1 of the manuscript, or discontinuities—phenomena that are statistically and physically implausible. To mitigate, prior methods based on first-order models typically rely on highly informative priors that encode the system’s behavior. However, such priors are rarely available in real-world scenarios, making this assumption impractical. \n\nBy accounting for momentum, the model can represent non-linear transitions and inflection points in the underlying dynamics, enabling more realistic interpolation. This is especially critical when observations are limited to coarse or widely spaced timepoints, as it allows the system to preserve natural tendencies like speeding up or slowing down. Moreover, second-order dynamics allow us to more faithfully capture the behavior of many natural physical and biological systems, which inherently follow second-order laws—such as Newtonian or Hamiltonian dynamics—thereby making our model both more physically plausible and better aligned with real-world systems. \n\nIn our proposed methodology, leveraging second-order dynamics allows us to infer physically plausible trajectories across arbitrary numbers of time marginals. By jointly optimizing across all timepoints, our method retrieves globally optimal couplings and captures the dynamics in a statistically and physically consistent manner. This results in long-range temporal inference that is temporally coherent and accurately represents the underlying physics of the system. Our experiments in Section 5 demonstrate that this leads to significantly improved performance over state-of-the-art methods, both qualitatively and quantitatively. In the revision, we have modified the final two paragraphs of the introduction to include the benefits of the formulation of the problem using a second-order SDE as discussed above.\n\n---\n\n### References\n\n1. Deep momentum multi-marginal Schrödinger Bridge\n2. Stochastic bridges of linear systems\n3. Multi-marginal Schrödinger Bridges\n4. A smoothest curve approximation\n5. Measure-valued spline curves: an optimal transport viewpoint"}, {"review_id": "cC0yKaYrf6", "reviewer": "Reviewer_mTyy", "summary": "The paper proposes 3MSBM, a simulation-free matching algorithm for multi-marginal Schrödinger Bridge problems in phase space. By deriving recursive formulas for optimal acceleration, it achieves smooth, temporally coherent trajectories while preserving marginal constraints throughout training. Experiments across synthetic and real-world datasets show that 3MSBM outperforms prior methods in both accuracy and efficiency.", "questions": "----\n* If the algorithm were implemented using only forward dynamics, any suboptimality in the initial Markov projection (i.e., the objective in Equation (9)) would propagate forward, accumulating bias and ultimately leading to a biased estimate of the MMSB solution. It appears that the authors do adopt a forward–backward scheme as described in [1], which mitigates this issue. However, to the best of the reviewer's understanding, this important detail is not clearly stated in the main paper. Since many potential readers—especially those working in population dynamics—may not be familiar with the Schrödinger Bridge Matching (SBM) framework, the reviewer recommends that the authors explicitly describe this forward–backward structure in the main paper, particularly in Section 3.3.\n----\n\n* For Theorem 3.5, it appears that the notation $\\mathbb{P}$ refers to the path measure of the underlying stochastic process (i.e., the solution to the SDE). However, this notation is used without being explicitly introduced, and related concepts such as the path space are also not defined. Providing a brief clarification would improve readability, especially for readers less familiar with measure-theoretic formulations of stochastic processes.\n\n----\n\n* The 3MBB framework introduces the dynamic programming principle by simulating backward ODEs to compute the optimal control. However, when the number of marginals $N$ is large, this backward integration may accumulate numerical errors, leading to a biased estimate of the initial control $u_0$. This bias subsequently propagates through the forward dynamics, since the forward bridge is conditioned on a suboptimal initial segment. While the authors derive an analytic expression for the optimal control in Eq. (7) and Appendix B.2, the derivation is restricted to up to five marginals. A closed-form analytic solution for general N could eliminate such approximation errors and improve scalability. Given the structure shown in the appendix, is it possible to derive a general recurrence relation for the $N$-marginal case?\n----\n\n* As the reviewer understand, SBIRR refines a predefined prior process in both the Lotka–Volterra (LV) and Gulf of Mexico (GoM) experiments, where the chosen priors already provide a strong approximation of the underlying dynamics. However, in general settings, assuming access to such informative priors that already capture the system behavior is often unrealistic. In contrast, the proposed 3MSBM appears to model the dynamics effectively without relying on strong prior assumptions. Could the authors clarify whether 3MSBM also benefits from a similar modeling choice, or does it learn the dynamics purely from data?\n\n----\n\n* Prior approaches [2] rely on computationally expensive Bregman iterations to enforce multiple constraints, resulting in long training times. Compared to this, I wonder how much training time the proposed method actually saves. Since SBM methods are generally known to converge faster than IPF-type objectives, and DMSB is based on such an IPF framework, it would be helpful to clarify whether 3MSBM also provides a significant advantage in terms of training efficiency.\n\n---\n\n* (Minor comment) Equation (4) is missing the differential $\\mathrm{d}t$ for the control.\n\n----\n    [1] Diffusion Schrödinger Bridge Matching. \n    [2] Deep Momentum Multi-Marginal Schrödinger Bridge.", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "**Strengths**\n\n* 3MSBM is a new matching framework that extending Schrödinger Bridges (SB) under multiple marginal constraints in phase space. This allows the model to generate smooth and temporally coherent trajectories that satisfy multiple positional constraints across time.\n\n* This work proposes a matching-based algorithm by deriving recursive formulas for optimal acceleration tailored to multi-marginal settings. These formulas extend Brownian bridges to momentum systems and enable simulation-free training using analytic Gaussian paths, making it possible to generalize previous matching frameworks to the phase-space multi-marginal SB setting.\n\n* Experiments on diverse datasets, demonstrate consistent improvement over state-of-the-art methods. 3MSBM achieves lower transport errors and better marginal fidelity in both low and high-dimensional regimes.\n\n----\n\n**Weaknesses**\n\nSee questions below.", "was_revised": true, "final_justification": "Momentum-SB methods capture underlying physics directly from data via their momentum term, providing a compelling framework for data-driven physical modeling. The proposed 3-MSBM preserves this strength while markedly reducing training time and memory consumption compared with earlier momentum-SB approaches. Experiments confirm its effectiveness and scalability. Reviewers pointed out that several key concepts were missing from the main text and had to be inferred from the figures; the authors have committed to clarifying these points in the revision.", "rebuttal": "We truly appreciate the encouraging comments of the reviewer, as well as their useful remarks. In the following, we address the points raised by the reviewer.\n\n### **1. Comment on Forward-Backward scheme**\n\nWe thank the reviewer for their comment. This is true, our algorithm does involve a forward-backward scheme, which will be explicitly stated in Section 3.3 to avoid confusion.\n\n### **2. Comment on formally defining path measure $\\mathbb{P}$**\n\nWe thank the reviewer for identifying this omission. In the revision, we will formally define the path measure $\\mathbb{P}$ associated with the learned dynamics, and $\\mathbb{P}^i$ as the learned dynamics at the $i^\\text{th}$ iterate of our algorithm before Theorem 3.5.\n\n### **3. Recursive formula for conditional acceleration**\n\nWe thank the reviewer for this important and insightful comment. While we agree with the reviewer’s conjecture that naive backward simulation of the ODEs in Eq (5) to compute the optimal control can scale poorly as $N$ grows, we emphasize that this is not practically implemented by our 3MBB. Specifically, Eq. (5) was never simulated throughout 3MBB training. Instead, we are able to compute an analytic solution to Eq. (5) for arbitrary $N$, as proposed by the reviewer.\n\nThe recursive expression for arbitrary $N$ is formalized in Prop 3.2, where the analytic formulation ensures that both the computational cost and approximation error remain independent of the number of marginals, effectively addressing the scalability and stability concerns associated with numerical solvers. Below, we present the general recursive formula, broken down in two steps, for obtaining the conditional bridge for an arbitrary $n^{\\text{th}}$ segment with $n<N-1$, i.e., $n$ is not the last segment. \n\n1. Expressions for the time-dependent functions $C_1^n(t), C_2^n(t), C_3^n(t)$\n    \n    For convenience, let us define the following functions corresponding to the $n^\\text{th}$ segment: $t\\in[t_n, t_{n+1})$\n    \n    $$ z_1^{ (n)}(t) = 3t - 3t_{n+1} - 3  \\qquad z_2^{ (n)}(t) = 3t - 4 - 3t_{n+1} \\qquad z_3^{ (n)}(t) = 6t - 6t_{n+1} - 3 \\qquad z_4^{ (n)}(t) = 6t - 6t_{n+1} - 4 $$\n    \n    $$ z_5^{ (n)}(t) = 4t - 4t_{n+1} - 3 \\qquad z_6^{ (n)}(t) = 4t - 4 - 4t_{n+1} \\qquad  z_7^{ (n)}(t) = 6t - 6t_{n+1} + 3 \\qquad z_8^{ (n)}(t) = 6t - 6t_{n+1} + 4 $$\n    \n    It follows that for the arbitrary $n^{\\text{th}}$ segment, the corresponding time-dependent functions $C_1^n(t), C_2^n(t), C_3^n(t)$, are given by  \n    $$C_1^{\\small (n)}(t) = \\frac{-3(\\alpha^{(n)} z_3^{(n)}(t) + \\beta^{(n)} z_4^{(n)}(t))}{(t-t_{n+1})^2*(\\alpha^{(n)} z_1^{(n)}(t) + \\beta^{(n)} z_2^{(n)}(t))} $$\n    $$C_2^{(n)}(t) = \\frac{3(\\alpha^{(n)} z_5^{(n)}(t) + \\beta^{(n)} z_6^{(n)}(t))}{(t-t_{n+1})*(\\alpha^{(n)} z_1^{(n)}(t) + \\beta^{(n)} z_2^{(n)}(t))} $$\n    $$C_3^{(n)}(t) = \\frac{6(\\alpha^{(n)}+ \\beta^{(n)})}{(\\alpha^{(n)} z_1^{(n)}(t) + \\beta^{(n)} z_2^{(n)}(t))}$$\n    \n    where $\\alpha^{(n)}, \\beta^{(n)}$ are segment-specific coefficients that shape the conditional bridge for the given segment $n$. These are recursively computed using the coefficients of the subsequent segment $n+1$, as follows:\n    $$\\alpha^{(n)} = \\alpha^{(n+1)} z_1^{(n)}(t_{n+1}) + \\beta^{(n+1)} z_2^{(n)}(t_{n+1}) $$\n    $$\\beta^{(n)} = 4(\\alpha^{(n+1)} +\\beta^{(n+1)})$$\n    \n    with the coefficients of the second-to-last segment $N-1$ given by: $\\alpha^{(N-1)}=0, \\ \\beta^{(N-1)}=1$.\n    \n2. The $\\lambda^{(n)}$ coefficients\n    \n    Because of the recursive nature of the backward ODE of $r_t$, with the terminal conditions depending on all succeeding segments, we find that the coefficients $\\lambda^{(n)}$ converge to the following exponential decay with an absolute error $<10^{-6}$ for segments more than 6 marginals before the last one:\n    \n    $\\lambda_j^{(n)} = (-1)^{j}1.608e^{-1.317(j-1)}$ ,   for     $j = 2, \\dots, N$         and       $\\lambda_1 = -\\sum_{j=2}^N \\lambda_j^{(n)}$\n    \n    Conversely, if the number of future marginals is smaller than 5, the $\\lambda^{(n)}$ are given by the values shown in the examples in Appendix B.2.1.\n    \n\nWe had initially presented the bridges for up to 5 segments for simplicity, but we will enrich Appendix B.2.1 with the discussion above to better describe the procedure and the steps needed to derive the conditional bridge for an arbitrary $n^\\text{th}$ segment: $t\\in[t_n, t_{n+1})$.\n\n### **4. Discussion on SBIRR**\n\nWe emphasize that, unlike SBIRR, which assumes prior knowledge of the underlying system behavior, our proposed 3MSBM learns the dynamics and achieves highly accurate trajectory inference purely from the data. This is made possible by two key components of our framework:\n\n1. By incorporating the velocity into the state, our 3MSBM captures *momentum-like effects* that encode trends in the data—such as speeding up or slowing down—that are otherwise lost in first-order models. This is particularly important in settings with sparse or irregular observations, as it enables the model to infer smoother and more physically realistic trajectories through intermediate timepoints, including non-linear transitions and inflection points.\n2. Our conditional bridge formulation solves for the joint optimal coupling across all marginals, which allows the dynamics at any given time to incorporate information from *all future pinned points*. This provides a form of implicit \"look-ahead\" guidance, helping the model infer trajectories that are globally consistent and faithful to the underlying physical or biological system.\n\n### **5. Comparisons between 3MSBM and DMSB**\n\nOur 3MSBM as a matching framework provides enhanced scalability compared to prior IPF based methods solving the mmSB problem, such as Deep Momentum Multi-Marginal Schrödinger Bridge (DMSB;[1]). To empirically verify our claims:\n\n1. We compare the total training time of our algorithm to DMSB on all datasets and observe that 3MSBM is consistently and significantly faster in wall-clock time across all datasets.\n    \n    \n    | Time (hours:min) | EB-100 | GoM | LV | BAQ |\n    | --- | --- | --- | --- | --- |\n    | DMSB | 5:51 | 2:55 | 2:43 | 3:15 |\n    | 3MSBM | 2:18 | 0:34 | 0:21 | 2:03 |\n\n2. We ablate the batch size in the EB-100 dataset and report the peak memory usage compared to DMSB. We find that our method consistently requires less memory than the IPF-based DMSB, which encounters a memory error for larger batch sizes. This enables significantly larger batch sizes, even in high-dimensional settings, providing greater design flexibility without compromising performance.\n    \n    | Batch Size | DMSB | 3MSBM |\n    | --- | --- | --- |\n    | 256 | 6.1 GB | 1.8 GB |\n    | 512 | 12.3 GB | 3.3 GB |\n    | 1024 | -- | 6.7 GB |\n    | 2048 | -- | 13.8 GB |\n   \n   \n### **6. Typo Correction**\n\nWe thank the reviewer for identifying this typo. It has been corrected in the revised manuscript.\n\n---\n\n### References\n\n1.  Deep momentum multi-marginal Schrödinger Bridge"}, {"review_id": "p6soKX6m2c", "reviewer": "Reviewer_Pewn", "summary": "The authors present a method for Schrödinger bridge matching via a phase space formulation and for multiple marginals. It alternates between a closed-form momentum bridge and drift fitting, preserving all snapshot constraints. Authors claim convergence and faster training. Experiments on four datasets show lower transport errors than prior methods.", "questions": "- What is the capital R in Algorithm 1? Previously this was defined as a matrix.\n- How does training cost scale with the number of time points? Do training costs ultimately constrain the batch sizes one can use, how does this in turn affect accuracy?\n- The alternating scheme between interpolating paths and coupling refinements feels like it could be brittle for less smooth trajectories or multi-modal marginals. Do the authors have any sense of the limitations of their method in terms of these issues?\n- It would be great to see a table (or something like this) clearly showing the differentiate factors and or advantages of their method compared to the many other \"multi-marginal\" and/or \"momentum\" Schrödinger Bridge approaches.", "limitations": "Yes.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths**\n- The paper is well written.\n- The motivation of trying to capture long range dependence in trajectory inference type problems is interesting.\n- The numerical results feel well founded and generally point to a strong method.\n\n**Weaknessess**\n- There are many variants for \"multi-marginal\" and/or \"momentum\" Schrödinger Bridge approaches in the literature. The authors do not do a great job of clearly showing the clear differentiate factors and or advantages of their method (aside from improvements on some benchmark problems).\n- Scaling claims are somewhat unsubstantiated; how does the ability to training scale with number of samples, dimensional of samples etc? Where do other methods fail at scaling? What experiments demonstrate this? The authors leave high dimensional problems to future work.", "was_revised": true, "final_justification": "The author resolved all my major concerns. I think this paper is a reasonable additional to the conference, but maybe not an extremely strong one.", "rebuttal": "We thank the reviewer for their positive comments, constructive criticism, and helpful suggestions. Below, we address the raised points.\n\n### **1. 3MSBM vs prior works**\n\nIn the revision, we elaborated on Table 1 and modified the Introduction to better present the differentiating factors of our 3MSBM compared to prior momentum and multi-marginal methods. The table below compares our method against several state-of-the-art multi-marginal approaches tailored for trajectory inference, based on three key attributes:\n\n- Simulation-free training: indicating scalability\n- Globally optimal coupling: assessing whether the method retrieves the mmSB solution versus relying on local approximations.\n- Smooth and coherent trajectories: reflecting second-order dynamics\n\n|  | SBIRR | MMFM | MMSFM | DMSB | Smooth SB | 3MSBM |\n| --- | --- | --- | --- | --- | --- | --- |\n| Simulation Free Training | ☑ | ☑ | ☑ | ☒ | ☒ | ☑ |\n| Globally Optimal Coupling | ☒ | ☒ | ☒ | ☑ | ☑ | ☑ |\n| Smooth Coherent Trajectories | ☒ | ☑ | ☒ | ☑ | ☑ | ☑ |\n\nWe group the prior work into 2 classes.\n\n1. Several recent methods propose interpolating across multiple timepoints by locally optimizing between adjacent marginals. While these approaches are generally scalable, they do not recover the global coupling.\n    - SBIRR [1] is a simulation-free iterative scheme performing trajectory inference via pairwise bridge optimization. Nevertheless, this pairwise formulation fails to capture global trajectory consistency. As a result, SBIRR relies on informative priors to guide trajectory inference, making it sensitive to the choice of reference class and often impractical for real-world applications.\n    - MMFM [2] performs cubic spline-based interpolation using deterministic dynamics. While scalable, it relies on precomputed piecewise OT couplings and hence does not provide the joint optimal coupling. This impairs the temporal coherence of the trajectories and the accurate capture of the dynamics.\n    - MMSFM [3] proposes a stochastic flow matching that optimizes measured value splines in position space using overlapping windows of three consecutive marginals. While this improves upon MMFM's pairwise interpolation, the lack of global coupling and reliance on first-order dynamics hinder its ability to learn smooth, temporally coherent trajectories.\n2. An alternative class of methods solves the mmSB without local approximations. However, these methods generally present low to moderate scalability capacity.\n    - DMSB [4] solves the mmSB in phase space using Bregman iterations and successfully retrieves the mmSB coupling. However, it suffers from scalability limitations due to the need to cache full SDE trajectories, leading to computational bottlenecks, error accumulation, and potential instability.\n    - Smooth SB [5] replaces Brownian priors with smooth Gaussian paths in phase space, achieving more temporally coherent and smooth trajectories. However, the belief propagation prohibits scaling in high dimensions.\n\nOur proposed algorithm directly addresses these limitations. The differentiating attributes of our 3MSBM compared to prior methods can be summarized as follows: \n\n- We lift the dynamics to phase space and minimize path acceleration, resulting in smoother low-curvature trajectories. This allows the model to capture trajectory trends typically lost in first-order systems—especially important when observations are sparse or irregular. As a result, the system can represent non-linear transitions and inflection points, enabling more realistic interpolation.\n- We solve the multi-marginal Brownian Bridge (3MBB; Theorem 3.1) objective analytically, and derive closed-form solutions for the conditional bridges across any number of $N$ marginals. This eliminates the need for costly numerical integrations of Eq. (5), which would scale poorly and accumulate approximation errors. By also avoiding approximations (e.g., local optimization of adjacent marginals and truncated policies), our methodology converges to the mmSB coupling across $N$ marginals.\n- Our matching framework allows us to learn the parameterized drift given a prescribed path in a simulation-free manner. Moreover, the iterative solutions of the learned drift remain close to the ground truth marginals throughout training, as opposed to prior methods that approximate the true marginals merely after convergence.\n\n### **2.  Scalability of 3MBB**\n\nOur algorithm adopts a two-step architecture, each of which avoids costly numerical simulations.\n\n- Closed-form conditional bridge:\n    \n    Notably, the analytic solution to the 3MBB objective allows us to compute the conditional bridge in Eq. (6) without requiring numerical integration or discretization, thereby eliminating both computational overhead and approximation errors. Consequently, this step remains independent of the number of marginals, as the conditional bridges are derived in closed form.\n    \n- Bridge Matching:\n    \n    The second step regresses the parameterized drift on the path prescribed in the previous step, in a simulation-free fashion, unlike more computationally expensive prior IPF-based methods.\n    \n\nNext, we provide 4 additional experiments to validate our claims, which are detailed below.\n\n1. We compare the total training time of our algorithm to DMSB across all experiments and observe that our 3MSBM is consistently and significantly faster in wall-clock time across all datasets.\n    \n    | Time (hours:min) | EB-100 | GoM | LV | BAQ |\n    | --- | --- | --- | --- | --- |\n    | DMSB | 5:51 | 2:55 | 2:43 | 3:15 |\n    | 3MSBM | 2:18 | 0:34 | 0:21 | 2:03 |\n\n2. We ablate the number of marginals in the EB-100 and GoM datasets, and report the training time per-epoch. We observe that the number of marginals has a negligible impact on training time per epoch, indicating that our algorithm scales well with increasing number of marginals.\n\n    - **Ablation on EB-100**\n    | Num of marginals | Time per epoch |\n    | --- | --- |\n    | 3 | 773 sec |\n    | 4 | 789 sec |\n    | 5 | 812 sec |\n\n    - **Ablation on GoM**\n    | Num of marginals | Time per epoch |\n    | --- | --- |\n    | 5 | 263 sec |\n    | 7 | 270 sec |\n    | 9 | 290 sec |\n\n3. We present a breakdown of the time percentage attributed n each of the two steps of our algorithm: i) the iterative propagation of the conditional dynamics (*lines 4-9 in Alg. 1)* and ii) the Bridge Matching step *(line 10 in Alg. 1)* take in the EB-100 and GoM datasets. We observe that the computational complexity introduced by the first step of the algorithm, i.e., the iterative propagation of the conditional dynamics, remains negligible compared to the Bridge Matching step.\n    \n    \n    |  | Cond. Bridge | Bridge Matching |\n    | --- | --- | --- |\n    | EB-100 | 3.3 % | 96.7 % |\n    | GoM | 7.4 % | 92.6 % |\n\n4. We ablate the batch size in the EB-100 dataset and report the peak memory usage compared to DMSB. We find that our method consistently requires less memory than the IPF-based DMSB, which encounters a memory error for larger batch sizes. This enables significantly larger batch sizes, even in high-dimensional settings, providing greater design flexibility without compromising performance. \n    \n    \n    | Batch Size | DMSB | 3MSBM |\n    | --- | --- | --- |\n    | 256 | 6.1 GB | 1.8 GB |\n    | 512 | 12.3 GB | 3.3 GB |\n    | 1024 | -- | 6.7 GB |\n    | 2048 | -- | 13.8 GB |\n\n### **3. Robustness to un-smooth and multi-marginal data**\n\nWe emphasize that a key feature of our method is the utilization of conditional bridges that explicitly condition on all future pinned marginals, as shown in Eq. (6). This design allows the learned dynamics to effectively “look ahead,” providing a global temporal context that enhances robustness to both irregular transitions and multi-modal marginals. \n\nRegarding robustness to less smooth data, our method outperformed state-of-the-art baselines on the Beijing Air Quality dataset (Section 4.3), comprising a time series with abrupt fluctuations and high volatility. To further validate the robustness of our method under multi-modal scenarios, we additionally consider:\n\n1. A sequence of marginals starting from a univariate Gaussian $(t_0)$, transitioning to a 4-mode GMM $(t_1)$, followed by an 8-mode GMM $(t_2)$, and returning to a univariate Gaussian $(t_3)$.\n    \n    \n    |  | $\\mathcal{W}_2$   @   $t_0$ | $\\mathcal{W}_2$   @   $t_1$ | $\\mathcal{W}_2$   @   $t_2$ | $\\mathcal{W}_2$   @   $t_3$ |\n    | --- | --- | --- | --- | --- |\n    | DMSB | 0.05 | 0.06 | 0.4 | 0.09 |\n    | 3MSBM | 0.01 | 0.03 | 0.07 | 0.01 |\n\n2. A more complex geometry starting from a 9-mode GMM arranged in a rhomboid shape $(t_0)$, transitioning to five Swiss rolls, with four placed at the corners of a rectangle and one in the center $(t_1)$.\n    \n    \n    |  | $\\mathcal{W}_2$   @   $t_0$ | $\\mathcal{W}_2$   @   $t_1$ |\n    | --- | --- | --- |\n    | DMSB | 0.07 | 0.14 |\n    | 3MSBM | 0.03 | 0.05 |\n\nIn both settings, we compute the $\\mathcal{W}_2$ distance between generated samples and ground truth marginals, comparing our method to DMSB. Our approach consistently yields lower $\\mathcal{W}_2$ errors, demonstrating superior fidelity in capturing complex and multi-modal dynamics.\n\n### **4. Typo Correction**\n\nWe thank the reviewer for pointing out this typo. It has been corrected: $R$  refers to the cost matrices in the optimal control problem. In Algorithm 1, the number of iterations is denoted by 'Max_Iter' and the loop index by 'iter'.\n\n---\n\n### References\n\n1. Multi-marginal Schrödinger Bridges with iterative reference refinement.\n2. Modeling complex system dynamics with flow matching across time and conditions.\n3. Multi-marginal stochastic flow matching for high-dimensional snapshot data at irregular time points.\n4. Deep momentum multi-marginal Schrödinger Bridge\n5. Trajectory inference with smooth Schrödinger Bridges"}, {"review_id": "7NmTzoi1NM", "reviewer": "Reviewer_wHqP", "summary": "This paper proposes 3MBSM, a method for performing distributionally optimal interpolation based on multiple snapshot data to infer the continuous-time dynamics of a system. The authors address the Momentum Multi-Marginal Schrödinger Bridge problem (a type of Schrödinger Bridge formulated in phase space), which compensates for the shortcomings of existing methods in capturing long-term temporal dependencies and ensuring trajectory continuity by simultaneously considering the coupling among multiple marginal distributions and incorporating both the position and momentum information of the system. Initially, the paper derives an analytical solution for the conditional acceleration function \\(a(t)\\) using an optimal control approach, and then employs neural networks to approximate the unconditional \\(a(t)\\). The proposed method is evaluated on multiple datasets, and the experimental results demonstrate that its error metrics for inferring trajectories at missing time points are superior to those of mainstream methods.", "questions": "Based on the above weaknesses, I have the following questions:\n\n1. The method in the paper estimates the initial velocity distribution through forward and backward iterations, with each iteration requiring numerical integration of the SDE. Would multiple numerical integrations significantly increase the computational complexity of the algorithm and affect its scalability? Could this issue be discussed?\n\n2. The efficiency in sampling of the proposed method stems from the fact that if the conditional probability path is Gaussian, the mean and covariance of the Gaussian distribution have analytical solutions. This arises because \\(a(t)\\) depends linearly on \\(x(t)\\) and \\(v(t)\\) (in essence, this again comes from the second-order approximation of the value function). Is this linearization of the system too crude? Can it handle systems with more complicated dynamics and strong nonlinearities? Are there any potential solutions to this problem?\n\n3. Could more experiments be included? For example, on the current datasets, it would be valuable to compare the proposed method with all baselines using a unified evaluation metric; perform quantitative evaluations on the single-cell dataset; and compare the proposed method with Neural-ODE approaches such as Mioflow [1] and DeepRUOT [2].\n\n4. What is the rationale behind selecting a second-order dynamical system in the paper? Using a second-order dynamical system for interpolation to obtain smoother trajectories appears to be a trivial conclusion, similar to how higher-order polynomial interpolation yields smoother trajectories. Could the authors discuss what additional information about the system dynamics is captured by using a second-order dynamical system for interpolation?\n\n5. Could a comparison be made between the proposed method and the method in DMSB (since both papers tackle the same problem)? Specifically, how much computational overhead is reduced by the current method compared to DMSB? From an optimal control perspective, which method is superior? (For instance, is there a significant difference in the cost functions of the trajectories obtained by the two methods, and which one yields lower cost?)\n\nOverall, I provide a positive and optimistic evaluation of this work. If the authors could address the issues raised, with additional experiments or further discussion, I would be very willing to raise my score further.\n\n[1] Manifold Interpolating Optimal-Transport Flows for Trajectory Inference, Guillaume Huguet, D.S. Magruder, Alexander Tong et al., NeurIPS 2022.\n\n[2] Learning stochastic dynamics from snapshots through regularized unbalanced optimal transport, Zhenyi Zhang, Tiejun Li, Peijie Zhou, ICLR 2025.", "limitations": "While this study provides valuable insights, a key limitation lies in the limited consideration of unnormalized distributions in the limitations section.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:**\n\n1. The paper is generally well-written and presents extensive experiments on multiple datasets, overall offering a comprehensive study.\n\n2. The theoretical analysis in the paper is sufficiently thorough. By employing a second-order approximation of the value function and using optimal control methods, the paper derives the differential equation governing the value function and the corresponding optimal control law, as well as the analytical solution for \\(a(t)\\). Additionally, the paper proves that the proposed iterative framework has a unique fixed point.\n\n3. The idea presented in the paper is very interesting. In the field of trajectory inference, most current work is based on a first-order dynamical system framework; using a second-order dynamical system to describe the underlying dynamics in the data is a promising and worthy direction for exploration.\n\n\n**Weaknesses:**\n\n1. The paper employs a training strategy similar to Flow Matching: first, a conditional probability path is obtained by solving an optimal control problem, and then a neural network is used to regress the unconditional \\(a(t)\\). However, before each update of \\(a(t)\\) in the training algorithm, iterations are required to estimate the initial distribution of \\(v(t)\\), and each iteration involves numerical integration of an SDE. This limits the scalability of the algorithm.\n\n2. The paper effectively linearizes the system, making the control law \\(a(t)\\) linearly dependent on \\(x(t)\\) and \\(v(t)\\), which restricts the solution accuracy of the algorithm.\n\n3. The experimental evaluation in the paper is not sufficiently comprehensive. Although the paper plans to compare the proposed 3MBSM method with SBIRR, DMSB, SmoothSB, and MMFM, only a subset of these methods is compared on each dataset (for example, only MMFM and DMSB are compared on the Beijing Air Quality dataset). Additionally, the evaluation metrics vary across datasets (e.g., the Gulf of Mexico dataset uses the W1 and W2 distances, while other datasets employ SWD and MMD). In the high-dimensional single-cell dataset, only qualitative evaluations are provided, with no quantitative assessments.\n\n4. The motivation for using a second-order dynamical system in the paper is unclear; the justification of \"yielding smoother trajectories\" does not appear to be a decisive factor for adopting a second-order system. Furthermore, employing a second-order dynamical system introduces additional complexity, potentially making the fitting process more challenging.\n\n5. The paper does not evaluate the method from an optimal control perspective, for instance, by comparing the cost functions of the trajectories with those obtained by other methods.\n\n6. The papre dose not consider the unbalancedness of the dataset, which is very important in the single cell dynamics.", "was_revised": true, "final_justification": "The author's reply has shown me the great potential of this method in learning complex dynamics with stochastic effect, addressing all my concerns. Assuming all these things are incorporated in the revised version, I am very happy to raise my score to 5.", "rebuttal": "We would like to thank the reviewer for their positive remarks, as well as their useful comments and further suggestions for improving our work. In the following, we address each of the concerns and suggestions of the reviewer.\n\n### **1. Computational complexity of solving 3MBB**\n\nWe highlight that the step of iteratively propagating the conditional dynamics—to approximate the initial velocity— introduces minimal computational overhead. Below, we empirically validate this claim by presenting a breakdown of the time percentage attributed to each of the two steps of our algorithm: i) the iterative propagation of the conditional dynamics (*lines 4-9 in Alg. 1)* and ii) the Bridge Matching step *(line 10 in Alg. 1)* take in the EB-100 and GoM datasets. We see that the computational complexity introduced by the first step of the algorithm, i.e., the iterative propagation of the conditional dynamics, remains negligible compared to the Bridge Matching.\n\n| Dataset | Cond. Bridge | Bridge Matching |\n| --- | --- | --- |\n| EB-100 | 3.3 % | 96.7 % |\n|    GoM  | 7.4 % | 92.6 % |\n\nThis efficiency stems from two key factors: i) convergence to the initial velocity is typically achieved within only ~10 iterations, ii) most importantly, the backward ODEs admit analytic solutions, allowing us to compute the conditional acceleration in closed form. This eliminates the need to perform costly numerical integration to solve Eq. (5), which would cause computational bottlenecks and error accumulation.\n\n### **2. Discussion on linearization**\n\nWe would like to thank the reviewer for this important and insightful question. \n\nWe emphasize that the linear expression of the conditional $a_t$ in Eq. (6) is not an approximation but the exact solution to the optimization problem in Eqs. (3) & (4). This comes from the fact that its objective (Eq. (3)) is quadratic in the augmented state $m_t$, and that its dynamics (Eq. (4)) are control-affine with additive noises. In such cases, the optimization is of a standard Linear-Quadratic Regulator (LQR) form, where the value function is known to be quadratic [3]. That said, the quadratic form of value function in L142 and its backward ODEs in Eq 5 are exact—rather than approximated, and, consequently, the optimal conditional $a_t$ is linear in $(x_t,v_t)$. We stress this distinction, but also acknowledge that current writing may confuse readers. We’ll update this part in the revision, and we thank the reviewer for raising the comments.\n\nRegarding the expression of Eq. (3), which relaxes the original interpolation problem (i.e., solving Eq. (2) with $m_n$ pinned at $\\bar m_n$ for all $n \\in N$) to a soft quadratic penalty. We note that this relaxation follows prior frameworks [1,2]—which transforms the problem into a multi-marginal conditional stochastic optimal control problem—and the soft penalty can be made hard by taking the penalty matrix $R$ to infinity, as illustrated in Prop 3.2.\n\n### **3. Second-order dynamics**\n\nTo capture coherent long-term dependencies in sparse, irregularly sampled time series, we incorporate second-order dynamics into our multi-point inference framework. By minimizing path acceleration, the model achieves low-curvature, smooth and consistent trajectories. Velocity provides a memory, enabling the system to preserve trends such as acceleration or deceleration that are lost in first-order (memoryless) models. This enables the capture of nonlinear transitions and inflection points, leading to more realistic interpolations that align with the behavior of physical systems—many of which (e.g., Newtonian or Hamiltonian systems) are governed by second-order dynamics. This is particularly important when observations are coarse or widely spaced. In contrast, first-order models often produce unnatural kinks or discontinuities (see Fig. 1).\n\nOur approach overcomes these challenges by solving the momentum multi-marginal Schrödinger Bridge problem, jointly optimizing across all timepoints. This yields globally optimal couplings and enables physically plausible long-range inference, as evidenced by the superior performance in Section 4.\n\n### **4. Additional Experiments**\n\nWe thank the reviewer for suggesting a more comprehensive experimental comparison. Below, we present additional results on the same datasets used in the main paper (GoM, LV, EB-100, and Beijing Air Quality), following the original experimental setups. Alongside the trajectory inference baselines included in the paper (i.e, SBIRR [4], Smooth SB [5], DMSB [6], and MMFM [7]) we also evaluate against two NeuralODE-based methods: MIOFlow [8] and DeepRUOT [9]. We use the **same metric** in all datasets: the **Sliced Wasserstein Distance (SWD)**.\n\nNote that the Smooth SB could not scale to the 100-dimensional single-cell dataset (EB-100), and thus, results are not available for this setting. Our 3MSBM consistently achieves the lowest or near-lowest SWD across all datasets and timepoints, demonstrating its ability to generate smooth, temporally coherent interpolations even in the presence of complex, multi-modal dynamics. Especially, on the single-cell EB-100 dataset, 3MSBM shows strong imputation capabilities, accurately interpolating both left-out marginals and achieving the lowest SWD scores.\n- **SWD on GoM**\n| Method | $t_1$ | $t_3$ | $t_5$ | $t_7$ |\n| --- | --- | --- | --- | --- |\n| DeepRUOT | 0.21 | 0.33 | 0.23 | 0.21 |\n| MIOFlow | 0.83 | 0.34 | 1.23 | 0.97 |\n| SBIRR | 0.15 | **0.11** | 0.11 | 0.09 |\n| Smooth SB | 0.17 | 0.14 | 0.10 | 0.13 |\n| DMSB | 0.23 | 0.54 | 0.39 | 0.28 |\n| MMFM | 0.23 | 0.25 | 0.10 | 0.19 |\n| 3MSBM | **0.14** | 0.14 | **0.08** | **0.06** |\n\n- **SWD on LV**\n| Method | $t_1$ | $t_3$ | $t_5$ | $t_7$ |\n| --- | --- | --- | --- | --- |\n| DeepRUOT | 0.30 | 0.16 | 0.22 | 0.44 |\n| MIOFlow | 1.53 | 1.49 | 1.23 | 1.46 |\n| SBIRR | **0.17** | 0.18 | 0.24 | 0.48 |\n| Smooth SB | 0.49 | 0.18 | 0.13 | 0.37 |\n| DMSB | 0.64 | 0.67 | 0.98 | 0.63 |\n| MMFM | 0.21 | 0.25 | 0.39 | 0.57 |\n| 3MSBM | 0.29 | **0.13** | **0.11** | **0.37** |\n\n- **SWD on EB-100**\n| Method | $t_1$ | $t_3$ |\n| --- | --- | --- |\n| DeepRUOT | 0.73 | 0.67 |\n| MIOFlow | 0.84 | 0.94 |\n| SBIRR | 0.80 | 0.91 |\n| Smooth SB | -- | -- |\n| DMSB | 0.58 | 0.54 |\n| MMFM | 0.59 | 0.76 |\n| 3MSBM | **0.48** | **0.38** |\n\n- **SWD on BAQ**\n| Method | $t_2$ | $t_5$ | $t_8$ | $t_{11}$ |\n| --- | --- | --- | --- | --- |\n| DeepRUOT | **13.67** | 52.60 | 71.34 | 84.67 |\n| MIOFlow | 46.64 | 79.06 | 76.06 | 60.87 |\n| SBIRR | 47.84 | 60.40 | 55.03 | 59.12 |\n| Smooth SB | 28.61 | 28.81 | 35.79 | 32.90 |\n| DMSB | 51.10 | 37.92 | 40.53 | 35.75 |\n| MMFM | 17.51 | 23.94 | 32.56 | 39.98 |\n| 3MSBM | 17.70 | **9.78** | **22.23** | **32.23** |\n\n### **5. Comparisons between 3MSBM and DMSB**\n\nIndeed, our 3MSBM as a matching framework provides enhanced scalability compared to prior IPF methods solving the mmSB problem. To empirically verify our claims:\n\n1. We compare the total training time of our algorithm to DMSB on all datasets and observe that 3MSBM is consistently and significantly faster in wall-clock time across all datasets.\n  \n    | Time (hours:min) | EB-100 | GoM | LV | BAQ |\n    | --- | --- | --- | --- | --- |\n    | DMSB | 5:51 | 2:55 | 2:43 | 3:15 |\n    | 3MSBM | 2:18 | 0:34 | 0:21 | 2:03 |\n\n2. We ablate the batch size in the EB-100 dataset and report the peak memory usage compared to DMSB. We find that our method consistently requires less memory than the IPF-based DMSB, which encounters a memory error for larger batch sizes. This enables significantly larger batch sizes, even in high-dimensional settings, providing greater design flexibility without compromising performance.\n    \n    \n    | Batch Size | DMSB | 3MSBM |\n    | --- | --- | --- |\n    | 256 | 6.1 GB | 1.8 GB |\n    | 512 | 12.3 GB | 3.3 GB |\n    | 1024 | -- | 6.7 GB |\n    | 2048 | -- | 13.8 GB |\n\nAdditionally, we evaluate our method from an optimal control perspective by computing the objective $\\int_0^1 \\|a_t\\|^2dt$ . We compare the performance of our 3MSBM with DMSB in terms of this control cost, which directly reflects the mmSB objective. Across all experiments, 3MSBM consistently achieves lower control cost, indicating greater optimality, while satisfying the marginal constraints more accurately, as shown in the experiments above.\n\n| $\\int_0^1 \\|a_t\\|^2dt$ | EB-100 | GoM | LV | BAQ |\n| --- | --- | --- | --- | --- |\n| DMSB | 2347.32 | 31.86 | 105.52 | 204.7 |\n| 3MSBM | 1160.96 | 18.45 | 89.75 | 158.9 |\n\n### **6. Unbalancedness**\n\nWhile unbalanced dynamics were not within the scope of our current work, we agree with the reviewer on their importance for modeling single-cell systems, and we thank them for this insightful suggestion. We plan to investigate this direction by extending our mmSB objective through the lens of unbalanced OT, in the spirit of the DeepRUOT framework. Specifically, we believe it is possible to incorporate a regularization term in the mmSB objective accounting for cell birth and death rates, analogous to the works in [1, 2]. That said, due to this generalization, it is likely that the analytic solution derived in Proposition 3.2 might not be applicable in this instance. Nevertheless, we agree with the reviewer that this is an intriguing direction for future work.\n\n---\n\n### References\n\n1. Generalized Schrödinger bridge matching\n2. Feedback Schrödinger Bridge matching\n3. Stochastic bridges of linear systems\n4. Multi-marginal Schrödinger Bridges with iterative reference refinement\n5. Trajectory inference with smooth Schrödinger Bridges\n6. Deep momentum multi-marginal Schrödinger Bridge\n7. Modeling complex system dynamics with flow matching across time and conditions.\n8. Manifold Interpolating Optimal-Transport Flows for Trajectory Inference\n9. Learning stochastic dynamics from snapshots through regularized unbalanced optimal transport."}]

Author rebuttal:

Meta-review: This paper introduces momentum multi marginal Schrödinger bridge matching, a method for inferring smooth and continuous-time trajectories of stochastic systems from multiple, sparse snapshots. The core innovation is to address this as a multi marginal problem in phase space. This method learns a single, globally coherent trajectory rather than relying on pairwise interpolation. 

Reviewers unanimously vote to accept this work after rebuttal. The strong rebuttal by the authors addressed all initial reviewer concerns and led two reveiwers to increase their score. The main strengths of this work are its theoretical soundness identifying and solving a key gap in the existing literature. Paired with strong empirical evaluation.

Initial weaknesses were noted around clarity and initial experiments, both fully resolved during the rebuttal. For these reasons I recommend acceptance at this time.
Decision: Accept (poster)