Multi-Marginal Stochastic Flow Matching

• Multi-Marginal Stochastic Flow Matching (MMSFM) is a framework that extends score and flow matching to learn continuous-time dynamics from sparse, multi-time marginals with optimal transport regularization.
• The method leverages variational and probabilistic formulations, employing tools like conditional regression, analytic interpolants, and spline or neural techniques to ensure smooth, globally coherent trajectories.
• MMSFM has practical applications in fields such as single-cell biology, oceanography, and meteorology, offering robust interpolation and improved accuracy over pairwise methods.

Searching arXiv for papers on MMSFM and closely related multi-marginal flow matching frameworks. Multi-Marginal Stochastic Flow Matching (MMSFM) is a multi-marginal extension of simulation-free score and flow matching methods for learning continuous-time stochastic or deterministic dynamics from discrete snapshot distributions observed at several time points. In this setting, one is given marginals such as $\mu_{t_0},\mu_{t_1},\dots,\mu_{t_K}$ or $\rho_0,\dots,\rho_M$ and seeks a flow, ODE, or SDE whose pushforward laws interpolate all prescribed marginals while retaining a least-action or entropy-regularized optimal-transport structure. Across recent formulations, MMSFM is closely connected to multi-marginal Schrödinger bridges, measure-valued splines, conditional flow matching, and dynamic optimal transport with intermediate constraints (Chen et al., 2019, Lee et al., 6 Aug 2025, Kansal et al., 3 Jun 2026).

1. Problem class and motivating setting

MMSFM addresses the regime in which only sparse sample “snapshots” of an evolving distribution are available, often at irregularly spaced times, and no ground-truth trajectories are observed. The central objective is to infer a continuous-time stochastic flow that passes through each observed marginal at its corresponding time, such as $p_{t_i}=\rho_i$ or $\rho_{t_k}=\mu_{t_k}$, while respecting temporal ordering and avoiding reduction to a sequence of unrelated two-point transport problems (Lee et al., 6 Aug 2025).

A recurring formulation is to learn a time-dependent vector field $u_t(x)$ or $v_t(x)$ such that $\dot X_t=u_t(X_t)$ in the deterministic case, or

$dX_t = u_t(X_t)\,dt + \sigma\,dW_t$

in the stochastic case, with the resulting pushforward marginals matching all observed snapshots. Several papers explicitly contrast this with pairwise interpolation between adjacent snapshots. In the pairwise setting, one solves disjoint two-point problems and stitches them together; the cited literature associates this with loss of global optimality, temporal coherence, nonsmooth trajectories, or “kinks” at intermediate times (Theodoropoulos et al., 11 Jun 2025, Islam et al., 3 Oct 2025).

The setting is especially prominent for single-cell biology, oceanographic and meteorological data, air quality, image progression, and economics. A defining aspect of the recent MMSFM literature is that intermediate marginals are not treated merely as evaluation checkpoints; they are incorporated into the training target, the variational action, the conditional path, or the coupling itself (Theodoropoulos et al., 11 Jun 2025, Lee et al., 6 Aug 2025, Kansal et al., 3 Jun 2026).

2. Variational and probabilistic formulations

A foundational antecedent is the multi-marginal Schrödinger bridge problem in phase space. There, the reference process is the inertial Wiener prior

$dx_t=v_t\,dt,\qquad dv_t=dw_t,$

with law $\mathbb P_0$, and the goal is to find a path law $\rho_0,\dots,\rho_M$0 that matches prescribed positional marginals at multiple times while minimizing $\rho_0,\dots,\rho_M$1 on path space. By Girsanov’s theorem, any admissible law with finite KL divergence corresponds to adding an acceleration drift $\rho_0,\dots,\rho_M$2, and the optimization acquires a Benamou–Brenier form with Fokker–Planck constraints and positional-marginal constraints. In the zero-noise limit, this converges to a measure-valued spline action, and the stochastic problem becomes a Fisher-regularized version of spline interpolation in Wasserstein space (Chen et al., 2019).

Momentum Multi-Marginal Schrödinger Bridge Matching (3MSBM) makes this phase-space construction explicit by lifting the state to $\rho_0,\dots,\rho_M$3 and using

$\rho_0,\dots,\rho_M$4

The control objective is

$\rho_0,\dots,\rho_M$5

subject to the SDE and the positional constraints $\rho_0,\dots,\rho_M$6 for all $\rho_0,\dots,\rho_M$7. In this formulation, the bridge is conditioned on all marginals simultaneously, and the optimal acceleration is obtained through a value function with matrix and linear terms, yielding backward ODEs with jumps at each observation time (Theodoropoulos et al., 11 Jun 2025).

A distinct line formulates the multi-marginal problem through dynamic optimal transport with time-localized potential terms. The classical Benamou–Brenier action

$\rho_0,\dots,\rho_M$8

is augmented by

$\rho_0,\dots,\rho_M$9

where $p_{t_i}=\rho_i$0 may be $p_{t_i}=\rho_i$1, MMD$p_{t_i}=\rho_i$2, or KL. The full action softly steers the flow toward intermediate marginals, and in the limit $p_{t_i}=\rho_i$3, $p_{t_i}=\rho_i$4 with $p_{t_i}=\rho_i$5, it recovers the hard-constrained piecewise OT problem (Kansal et al., 3 Jun 2026).

In simulation-free stochastic formulations, the drift is often decomposed as

$p_{t_i}=\rho_i$6

with separate networks for the transport velocity and the score. This reparameterization allows flow matching and score matching to be trained jointly from analytically tractable conditional paths rather than by simulating the full SDE during training (Lee et al., 6 Aug 2025).

3. Conditional paths, interpolants, and measure-valued splines

The practical behavior of MMSFM depends heavily on the choice of conditional path used to generate regression targets. Recent work offers several constructions.

Formulation	Conditional path	Stated role
MMSFM for irregular time points (Lee et al., 6 Aug 2025)	Monotonic cubic Hermite spline over overlapping triplets with Gaussian bridge variance	Handles irregular snapshot timing
ALI-CFM (Kviman et al., 1 Oct 2025)	Neural interpolant $p_{t_i}=\rho_i$7	Matches intermediate marginals adversarially
SplineFlow / IMMFM (Rathod et al., 30 Jan 2026, Islam et al., 3 Oct 2025)	B-spline or piecewise-quadratic interpolant	Enforces smooth multi-marginal targets

In the irregular-time MMSFM formulation, the full time interval is partitioned into overlapping windows of $p_{t_i}=\rho_i$8 consecutive marginals, with $p_{t_i}=\rho_i$9 in the reported experiments. For each window, samples are aligned into matched tuples using a first-order-Markov approximation to the multi-marginal OT plan, and a Euclidean transport spline $\rho_{t_k}=\mu_{t_k}$0 is fit through the matched points. The paper uses monotonic cubic Hermite splines to avoid overshoots on irregular intervals. The conditional density is then Gaussian,

$\rho_{t_k}=\mu_{t_k}$1

with Brownian-bridge variance $\rho_{t_k}=\mu_{t_k}$2, $\rho_{t_k}=\mu_{t_k}$3, and analytic drift and score targets (Lee et al., 6 Aug 2025).

ALI-CFM replaces handcrafted interpolants with neurally parameterized curves. Intermediate marginals are imposed by a GAN-style loss at each observation time, and regularization toward linear, piecewise-linear, or low-curvature references is used to ensure uniqueness and smoothness. The paper states existence and almost-everywhere uniqueness for the linear-reference and piecewise-linear-reference regularizers under the given assumptions (Kviman et al., 1 Oct 2025).

SplineFlow uses B-spline interpolation across all anchors. The path is

$\rho_{t_k}=\mu_{t_k}$4

where $\rho_{t_k}=\mu_{t_k}$5 are Cox–de Boor basis functions of degree $\rho_{t_k}=\mu_{t_k}$6. The reported properties are local support, non-negativity, partition of unity, $\rho_{t_k}=\mu_{t_k}$7 smoothness, and approximation error $\rho_{t_k}=\mu_{t_k}$8 for smooth targets, compared with $\rho_{t_k}=\mu_{t_k}$9 for linear interpolation. The stated motivation is that B-splines do not exhibit Runge’s phenomenon and yield stable higher-order multi-marginal paths (Rathod et al., 30 Jan 2026).

IMMFM uses a piecewise-quadratic path over triples of consecutive observations. For $u_t(x)$0,

$u_t(x)$1

with $u_t(x)$2. This produces a continuously varying velocity and is designed to avoid the velocity jumps of naively chained linear paths (Islam et al., 3 Oct 2025).

A distinct construction appears in OTP-FM. There, the straight-line conditional FM target is corrected by potential-induced force terms derived from the variational action. The resulting trajectories are described as spatiotemporally flexible while remaining simulation-free during training (Kansal et al., 3 Jun 2026).

4. Learning objectives and algorithms

The dominant training paradigm is conditional regression onto analytically known conditional drifts or velocities. In the irregular-time MMSFM model, the final objective is

$u_t(x)$3

with stratified time sampling over each window and $u_t(x)$4. The learned SDE is

$u_t(x)$5

and training is simulation-free in the sense that full trajectories are not simulated during training (Lee et al., 6 Aug 2025).

3MSBM alternates a closed-form bridge step with a matching update. A minibatch of points is sampled from the current coupling, backward ODEs for $u_t(x)$6 and $u_t(x)$7 are solved segment by segment, a closed-form linear acceleration field $u_t(x)$8 is constructed, Gaussian phase-space trajectories are sampled analytically, and a parametric acceleration network $u_t(x)$9 is updated by $v_t(x)$0 regression. Because $v_t(x)$1 is matched to a bridge that exactly attains the pinned marginals, the paper states that the learned process remains consistent with $v_t(x)$2 throughout training (Theodoropoulos et al., 11 Jun 2025).

ALI-CFM is trained in two stages. First, the interpolant $v_t(x)$3 is learned with the adversarial interpolant objective by alternating discriminator ascent and interpolant descent. Second, a time-dependent network $v_t(x)$4 is trained by conditional flow matching on $v_t(x)$5. This separates marginal fitting at observation times from vector-field learning (Kviman et al., 1 Oct 2025).

OTP-FM extends conditional FM by adding OT-potential corrections to the target velocity and then learning a consistency model $v_t(x)$6 through improved MeanFlow self-distillation. The paper emphasizes a simulation-free algorithm, fixed-point updates for intermediate states, and few-step consistency inference that requires only two-step evaluations rather than ODE/SDE simulation at test time (Kansal et al., 3 Jun 2026).

A common algorithmic burden is the coupling stage. In the irregular-time MMSFM formulation, mini-batch OT is implemented with EMD, which has worst-case cost $v_t(x)$7 per mini-batch; in SplineFlow, the paper states that spline setup is minor in practice because $v_t(x)$8 and most cost lies in network evaluation (Lee et al., 6 Aug 2025, Rathod et al., 30 Jan 2026).

5. Empirical behavior and reported benchmarks

The reported evaluations span synthetic 2D systems, single-cell RNA-seq, spatial transcriptomics, ocean currents, air quality, image progression, cell tracking, and longitudinal neuroimaging. Common metrics include Wasserstein-2, Sliced-Wasserstein, Maximum Mean Discrepancy, Wasserstein-1, and Earth-Mover’s Distance (Theodoropoulos et al., 11 Jun 2025, Lee et al., 6 Aug 2025).

For 3MSBM, the reported benchmarks are Lotka–Volterra predator–prey, Gulf-of-Mexico ocean currents, Beijing air quality, and Embryoid Body single-cell sequencing. The paper reports that 3MSBM consistently attains the lowest interpolation error on held-out marginals, with an example of $v_t(x)$9 versus $\dot X_t=u_t(X_t)$0–$\dot X_t=u_t(X_t)$1 for pairwise methods, while preserving training-set marginals within $\dot X_t=u_t(X_t)$2. It also reports globally smooth trajectories, approximately $\dot X_t=u_t(X_t)$3 faster training than deep momentum mmSB (DMSB), and memory usage growing only linearly in $\dot X_t=u_t(X_t)$4 and $\dot X_t=u_t(X_t)$5 (Theodoropoulos et al., 11 Jun 2025).

For OTP-FM, the reported results include EB and CITE-seq single-cell datasets, Gulf of Mexico ocean currents, and Beijing air quality. The paper reports state-of-the-art held-out interpolation on EB 100D leave-two-out with $\dot X_t=u_t(X_t)$6 average MMD versus $\dot X_t=u_t(X_t)$7 for OT-CFM, training times of $\dot X_t=u_t(X_t)$8–$\dot X_t=u_t(X_t)$9 minutes on an L40S GPU, and interpolation improvements of $dX_t = u_t(X_t)\,dt + \sigma\,dW_t$0–$dX_t = u_t(X_t)\,dt + \sigma\,dW_t$1 over spline-based MMFM and 3MSBM on the ocean and air-quality benchmarks (Kansal et al., 3 Jun 2026).

For the irregular-time MMSFM model, the paper reports that Triplet ($dX_t = u_t(X_t)\,dt + \sigma\,dW_t$2) outperforms Pairwise ($dX_t = u_t(X_t)\,dt + \sigma\,dW_t$3) and MIOFlow on synthetic data, particularly under irregular timing; that on DynGen all methods struggle with true bifurcation but MMSFM still yields lower $dX_t = u_t(X_t)\,dt + \sigma\,dW_t$4; that on CITE-seq and Multiome, Triplet improves over Pairwise by $dX_t = u_t(X_t)\,dt + \sigma\,dW_t$5–$dX_t = u_t(X_t)\,dt + \sigma\,dW_t$6 in $dX_t = u_t(X_t)\,dt + \sigma\,dW_t$7 and MMD; and that on Imagenette, Triplet gives more stable losses and visually coherent intermediate images (Lee et al., 6 Aug 2025).

ALI-CFM reports that only ALI gives smooth interpolants and low-variance gradients on a synthetic “knot” task with up to $dX_t = u_t(X_t)\,dt + \sigma\,dW_t$8 marginals, that it yields smooth and accurate trajectories on a cell-tracking benchmark with $dX_t = u_t(X_t)\,dt + \sigma\,dW_t$9 video frames, and that on spatial transcriptomics it achieves the lowest EMD of $dx_t=v_t\,dt,\qquad dv_t=dw_t,$0 versus $dx_t=v_t\,dt,\qquad dv_t=dw_t,$1 for OT-CFM, $dx_t=v_t\,dt,\qquad dv_t=dw_t,$2 for OT-MMFM, and $dx_t=v_t\,dt,\qquad dv_t=dw_t,$3 for OT-MFM (Kviman et al., 1 Oct 2025).

IMMFM is evaluated on synthetic 2D benchmarks, Starmen video data, and ADNI, MS lesion, and GBM longitudinal imaging. The paper reports forecast MSE reductions of $dx_t=v_t\,dt,\qquad dv_t=dw_t,$4–$dx_t=v_t\,dt,\qquad dv_t=dw_t,$5 versus TFM and multi-marginal ODE baselines on synthetic data, PSNR improvements of $dx_t=v_t\,dt,\qquad dv_t=dw_t,$6–$dx_t=v_t\,dt,\qquad dv_t=dw_t,$7 dB and SSIM gains of $dx_t=v_t\,dt,\qquad dv_t=dw_t,$8–$dx_t=v_t\,dt,\qquad dv_t=dw_t,$9 on clinical imaging, and an increase in AD–CN classification accuracy from $\mathbb P_0$0 using only observed data to $\mathbb P_0$1 using IMMFM forecasts (Islam et al., 3 Oct 2025).

6. Relation to adjacent methods, limitations, and extensions

The literature draws a sharp distinction between global multi-marginal matching and pairwise stitching. One recurrent misconception is that multi-marginal flow matching is merely a post hoc concatenation of two-point flows. The cited papers state the opposite: global bridges, spline paths, or OT-potential actions are designed precisely to impose all marginals simultaneously and thereby preserve long-range temporal dependencies and temporal coherence (Theodoropoulos et al., 11 Jun 2025).

A second misconception is that all such methods are deterministic ODE models. In fact, some formulations are explicitly stochastic, using an SDE drift together with score matching and Brownian noise, while others learn a probability-flow ODE or a consistency model. This suggests that “MMSFM” is better understood through its constraint structure and training target than through a single fixed dynamical formalism (Lee et al., 6 Aug 2025, Kansal et al., 3 Jun 2026).

The reported limitations are also method-specific. In the irregular-time MMSFM model, mini-batch OT couplings can swap paths in branching flows, the first-order Markov approximation to the multi-marginal plan may fail when long-range couplings matter, and OT cost grows cubically in batch size (Lee et al., 6 Aug 2025). Earlier multi-marginal Schrödinger bridge solvers are described as suffering from trajectory caching or piecewise modeling, and OTP-FM is presented as an attempt to avoid these bottlenecks (Kansal et al., 3 Jun 2026). In bifurcating systems, the published evaluations explicitly note that all methods may struggle (Lee et al., 6 Aug 2025).

The extension space is broad. Reported directions include high-resolution video interpolation through phase-space lifts, control-aware generative models with nonzero deterministic drift, mixture-of-bridges for multimodal dynamics, approximate OT or entropic regularization for large mini-batches, and causal inference by comparing learned drifts under perturbations (Theodoropoulos et al., 11 Jun 2025, Lee et al., 6 Aug 2025). A plausible implication is that future work will continue to move between three poles already present in the literature: hard marginal constraints via Schrödinger bridges, soft intermediate constraints via OT potentials, and learned interpolants that trade analytic structure for flexibility.