Robust Generalized Schr\"{o}dinger Bridge via Sparse Variational Gaussian Processes
Authors: Minyoung Kim
Venue: ICLR 2026 Poster (2026)
Keywords: Gaussian processes
Distribution matching
Diffusion models
Bayesian statistics
Primary area: generative models
TL;DR: We propose a generalized Schr\"{o}dinger Bridge matching algorithm that is robust to uncertain/noisy stage costs, by formulating a Gaussian process inference problem on the marginal probability paths.

Abstract: The famous Schr\"{o}dinger bridge (SB) has gained renewed attention in the generative machine learning field these days for its successful applications in various areas including unsupervised image-to-image translation and particle crowd modeling. Recently, a promising algorithm dubbed GSBM was proposed to solve the generalized SB (GSB) problem, an extension of SB to deal with additional path constraints. Therein the SB is formulated as a minimal kinetic energy conditional flow matching problem, and an additional task-specific stage cost is introduced as the conditional stochastic optimal control (CondSOC) problem. The GSB is a new emerging problem with considerable room for research contributions, and we introduce a novel Gaussian process pinned marginal path posterior inference as a meaningful contribution in this area. Our main motivation is that the stage cost in GSBM, typically representing task-specific obstacles in the particle paths and other congestion penalties, can be potentially noisy and uncertain. Whereas the current GSBM approach regards this stage cost as a noise-free deterministic quantity in the CondSOC optimization, we instead model it as a stochastic quantity. Specifically, we impose a Gaussian process (GP) prior on the pinned marginal path, view the CondSOC objective as a (noisy) likelihood function, and infer the posterior path via sparse variational free-energy GP approximate inference. The main benefit is more flexible marginal path modeling that takes into account the uncertainty in the stage cost such as more realistic noisy observations. On some image-to-image translation and crowd navigation problems under noisy scenarios, we show that our proposed GP-based method yields more robust solutions than the original GSBM.

Published as a conference paper at ICLR 2026
ROBUST GENERALIZED SCHRÖDINGER BRIDGE VIA SPARSE VARIATIONAL GAUSSIAN PROCESSES
Minyoung Kim
Samsung AI Center Cambridge, UK mi kimm1 @mma coom
ABSTRACT
The famous Schrödinger bridge (SB) has gained renewed attention in the gener- ative machine learning field these days for its successful applications in various areas including unsupervised image-to-image translation and particle crowd mod- eling. Recently, a promising algorithm dubbed GSBM was proposed to solve the generalized SB (GSB) problem, an extension of SB to deal with additional path constraints. Therein the SB is formulated as a minimal kinetic energy conditional flow matching problem, and an additional task-specific stage cost is introduced as the conditional stochastic optimal control (CondSOC) problem. The GSB a emerging problem with considerable room for research contributions, and we introduce a novel Gaussian process pinned marginal path posterior inference as a meaningful contribution in this area. Our main motivation is that the stage cost in GSBM, typically representing task-specific obstacles in the particle paths and other congestion penalties, can be potentially noisy and uncertain. Whereas the current GSBM approach regards this stage cost as noise-free deterministic quantity in the CondSOC optimization, we instead model it as a stochastic quantity. Specifically, we impose Gaussian process (GP) prior on the pinned marginal path, view the CondSOC objective as (noisy) likelihood function, and infer the posterior path via sparse variational free-energy GP approximate inference. The main benefit is more flexible marginal path modeling that takes into account the uncertainty the stage cost such as more realistic noisy observations. On some image-to-image translation and crowd navigation problems under noisy scenarios, we show that our proposed GP- GP-based method yields more robust solutions than the original GSBM.
INTRODUCTION
The bridge matching problem is to find a stochastic/ordinary differential equation (SDE/ODE) that bridges two given distributions (Lipman et al.] 2023 Albergo & Vanden-Lijndenn 2023 Neklyudov
et al. 2023 Villani
2009 Peyré & Cuturi (017). ). Its application areas are enormous, which also
subsumes the recent generative modeling (e.g.. conditional or unconditional image generation) (Sohl
Dickstein et al. 2015 Chen et al. 2018 Ho et al.] 2020 Song et al. (021). The famous Schrödinger
bridge (SB) problem (Fortet] 1940) Kullback] T 1968; Rüschendorf 1995 De Bortoli et al., 202
Vargas et al. 2021 Shi et al. 2023), a widely known problem for century, is a special instance of the bridge problem. The SB problem has recently received renewed attention in the generative machine learning and related fields for its applications to unsupervised image-to-image translation (Kim)
2024), image inpainting (Wang et al. 2021). , fluid modeling (Shi et al., 2023), molecular simulations Noé et al. 2020), robot navigation (Liu et al. 2018). and more.
Among these recent trends and large body of works, the generalization of SB (GSB), incorporating additional costs on marginal paths, has started drawing attention as an extension of the SB problem. One of the most promising methods called GSBM (Liu et al., 2024) was proposed, in which the minimal kinetic energy conditional flow matching problem (Lipman et al.) 2023 Tong et al., 2023), and an additional task-specific stage cost is introduced in the conditional stochastic optimal control (CondSOC) problem. Incorporation of task task-specific stage costs can give better guidance towards desired probability paths, or penalize prohibitive paths, which is desirable for realistic crowd navigation problems and latent-variable generative modeling problems.
SB is formulated
GSBM opens up intriguing research possibilities, among which in this paper we are particularly interested in the CondSOC path optimization part. Although GSBM originally adopted deterministic

---

Published as a conference paper at ICLR 2026
spline-based optimization for CondSOC, which can be effective in many existing domains, there several rooms to improve. One is more flexible modeling. However, simply increasing polynomial degrees (e.g., cubic splines) is often prone to fail in practice mainly due to numerical issues, while path integral resampling (Kappen) 2005) for non-Gaussian path modeling is computationally expensive with some technical issues discussed in Sec. [
Another option for improvement is to take into account uncertainty in path modeling in order to lead to more robust GSB solutions.
The Bayesian Gaussian process (GP) approach can nicely fulfill these two desiderata. To this end, we view the CondSOC objective as the likelihood function, and impose a GP prior path measure to the mean and stdev paths. With sparse variational free energy formulation adapted from (Titsias) 2009). we can infer the posterior path measure effectively. Although this can incur additional computational complexity in the CondSOC optimization compared to GSBM's spline optimization, we are able to achieve significantly more flexible and robust solutions than GSBM by properly handling uncertainty in potentially noisy stage cost observation. Also, our variational GP inducing variables merely take up roughly twice as many in GSBM's. Tested on several crowd navigation and image-to-image translation problems in noisy observation scenarios, we show that our GP-based approach is more robust noise than GSBM by capturing uncertainty that may reside in the GSB problems.
BACKGROUND
While here briefly provide background material required to understand our approach, for more comprehensive and detailed background, we refer the readers who are less familiar with this topic to our unified bridge algorithm (UBA) framework in Appendix A (also in stand-alone article (Kim) 2025)). It not only gives introduction to the (conditional) flow matching (CFM) and Schrödinger bridge (SB) problems, but also offers a unified SDE framework that subsumes both CFM and SB into one.
Bridge problems. The SB is special instance of the more general problem class called the bridge problem which aims to find an SDE (or sometimes an ODE) that bridges two given distributions. More specifically, given two distributions no(-) and m (-), the goal is to find an SDE, specifically the drift function u, (x), with specified diffusion coefficient a (a =0 for ODE),
dx, ue (x.) + adW, Zo To~))
(1)
with the constraint that the state at time f== conforms to the other target Ti, that is, I1 (.). We can also find the reverse-time SDE for the bridge the other way around, i.e., starting from Ti and landing at To. Once solved, the solution to the bridge problem can give us the ability to sample from one of o{ given the samples from the other. For instance, in typical generative modeling, mo is usually a tractable density like Gaussian, while πı is a target distribution that we want to sample from. In the bridge problem, however, To can also be an arbitrary distribution beyond tractable densities.
Schrödinger bridge (SB) has an additional constraint that the probability path measure should be the one that is closest to given reference path measure (Appendix (A.2.2). Denoting the path measure of (1) by Pu and the reference measure by pref, the SB can be defined as the following optimization:
min (L(PI||pr))) s.t. (x)) no (2n), (1) (x)) n (21)
Without loss of generality, we use the Brownian SDE, dx, = odW, as a reference measure pref Several algorithms have been proposed to solve the SB: the traditional iterative filtering (IPF) algorithm, its extensions with modern deep neural nets, and mini-batch versions of (entropic) optimal transport couplings in CFMs. See details and related references in Appendix A...33 and A.2.1 The computational overheads and known drawbacks of these methods were addressed in the recent iterative projection method (IMF) also known as the DSBM algorithm (Shi et al. 2023) (Appendix A.2.4).
Generalized Schrödinger bridge (GSB) problem. The standard SB has been extended by certain state costs introduced to express application-specific preferences or penalties for state paths { in SB (Chen et al.) 2015 Liu et al., 2022, Chen [023 Liu al.) 2024). For instance, the state cost can encode particle interaction costs (Gaitonde et al. 2021). quantum potentials (Philippidis et al. 1979). or geometric costs (Liu et al., 2024). The GSB problem is what mainly deal with in this paper, and we particularly focus on the recent GSB matching (GSBM) algorithm (Liu et al. 2024) for its superiority to its predecessors including (Liu et al. 2022 in terms of the quality of solutions.
2

---

Published as a conference paper at ICLR 2026
For smooth exposition, we describe the GSBM algorithm within our UBA framework (Appendix A or (Kim) (. particular by showing how it can be extended from the minimal kinetic UBA form of the IMF/DSBM algorithm (Shi et al. 2023) described in Alg. The following is concise summary.
For the given state (or stage) cost (z),, repeat the following two steps (A and B) for sufficient number times to find a neural-net based SDE, dz, = va(t, dt + odW;
(a) x )l1 | va(ct) dt, (2o, zi)xz)))((0,,1 (3) (a((oo))
(Step-A) arg min (1urTcc
where (r(x x1) = NN; Me, 12 is the pinned marginal at for coupling (4o, #1) sampled from:
initially (when the previous iterate gold is not available) Q20 21) pugold (1)) driven by dr, = Uigoia (1, se)dt + adW, 30 ~ T
(
and a (x4 ro, v1) determined from (x( = j4 + at - pt), 0t o== -t/r//
by the following formula (f+ time derivative):
(5)
Once (e( is optimized for each coupling (Io. x1). we update the neural net by:
(Step-B)
0 +- arg min E,
(2 (2)) - no(t, 22))l2
(6)
The rationale behind this stochastic optimal control (SOC) formulation is detailed in Appendix A.3.3 The appearance of the stage cost (x( in 3 is the only difference between DSBM and GSBM (i.c.. it can be shown that the above procedure recovers DSBM exactly if 0). Following (Liu et al. 2024) we call the optimization (3) the CondSOC problem. In CondSOC, the stage cost V is added to the kinetic energy term where the stage cost. is task-specific, and expresses our problem-specific preferences or penalties (e.g., obstacles in crowd navigation or desired particle paths). In the GSBM algorithm, the mean (4 and stdev "tt functions the Gaussiar (x(22zz 11), i.e., the optimization variables of CondSOC, are parametrized as splines at some knot points.
GAUSSIAN PROCESS MODELING FOR PINNED MARGINAL PATHS
We denote the CondSOC objective in GSBM as J(P.; V.) where P. (P #1)}e is the pinned marginal path, and V. := [V((-)} is the task task-specific stage cost function. That is,
J J.P:V.)) =
(VV))=/En((liee (r((ottooo [ x1)1 11) Vi(x.)) dt
(7)
where au (e(2z/20. 20. 21) is determined from P. by (5. So. GSBM aims to find P. by solving: arg minp, [[P.; V.). Note that in GSBM the pinned marginal path P. is the deterministic op- timization variables of the CondSOC, and they aim to find a point estimate. We instead treat P. as random variate, and place some prior distribution Pprior (P.). For instance, preference of the SB's linear pinned marginal path over other less smooth ones through this prior. regularizer enforce our prior preference. With
impose
Then " lllog log priior (P.) )" can be introduced some balancing coefficient (2 0), what we call the regularized-pan GSBM can be written as:
arg min J(P.: V.)) T log P Pprior (P.) arg max P Pprior (P.) (Pt: (8) P.
Prior
Likelibood
where the equivalence comes from exponentiation. This can be seen as a maximum-a-posterioo (MAP) solution. But we go one step further to consider the posterior distribution,
Prost (Pa)
Pprtoor (P.) (xp(-N(P::Vo)/7
Although in (iu et al. 2024 they also proposed non Gaussian path modeling using the so-called path integral resampling, there are several technical issues, including: i) the related on (x 220 æ1) may not admit the non-Gaussian Ps P= 21) as its marginal distribution, and ii) the importance sampling resampling procedure is computationally very overwhelming practice. So, we adhere the Gaussian path recipe described here.
3

---

Published as a conference paper at ICLR 2026
We specifically consider Gaussian process (GP) prior for the mean and stdev functions (N. freh and % := (1))). For computational tractability, we impose independent prior modeling with dimension-wise factorization for p4- To ensure positivity of % we use the parameterization V a /t(1 -t) log(1 e in which we impose GP prior on the unconstrained 9. Specifically,
d
Pprior (P.) = GP( GP((.) = II gP(r): mid Myd) mp(
(10)
j=1
where mir4) and k( are mean and covariance functions of GP for j and %. We use superseript j for the j-th dimension (e.g., nee is the j-th function of the d-dimensional vector function (o)). A main benefit of this Bayesian path modeling is that by treating the potentially noisy CondSOC objective (e.g., resulting from noisy stage cost V) as a (stochastic) likelihood function, we can have a more robust path solution as a posterior process than simply point estimate that fully trusts the objective J. Another benefit is that GP can lead to more natural and flexible path modeling.
Since we deal with (Io, F1 )-pinned marginal paths, we need to ensure 14g = zo. 141 = 21, To = 7n1 O for consistency with conditioning. For 7/ we already meet these boundary conditions due to our construction. For pe, we deal with the conditional process by conditioning mg = and μ¹ æf for each j, which is also GP. More specifically, the conditional GP prior for pd can be expressed as:
Ad (mu = zd uî = =1) ~ GP(M (!!!)) where
(11)
Mr mtt + M,,,,,1) (k
(0.1))(0.))) " my ot -m²l
(12)
1 ke.ed - 1o.1. (k (o.1), (f0,1,,
(13)
where kA, for index sets A, B is defined to be the (| x [B)) matrix whose (i entry is KA,
In all our experiments, the mean function for 14 is chosen as the linear interpolation, m = (1 (1-t)a0 t)xo
1 and the mean function of m is set to constant log(e 1), which in turn leads to o /((t) for yt. This choice makes the prior mean coincide with the solution of DSBM in the original SB problem. For the covariance functions, we use the squared exponential kernel function (See ablation study on different kernel choices in Sec. 5.4). where the related kernel hyperparameters, denoted by 17, are chosen by empirical Bayes (as discussed in the variational inference section below). For simplicity, we assume the same kernel function and hyperparameters for all dimensions of μtc.
3.1 SPARSE VARIATIONAL (FREE-ENERGY) GP POSTERIOR INFERENCE
Due to the intractability of the posterior P post (P.) defined in
we adopt the sparse variational
GP approximate inference (Titsias 2009 Deztouli & Bonilla a 2015 Matthews et al. 2016 Bauer
al. 2016) that is based on the variational free energy principle. First off, we consider the fully factorized GP as tractable variational process family, specifically Q(P.) = Q(u.) Q(%+)
I1 2(x)) 2. Both Q(jd) and @(f+) are chosen to be Gaussian processes parametrized by the inducing- point processes. For exposition we focus on 2(rr) here since derivations for 2(%.) are similar. For notational simplicity we will often drop superscript dependency on μ (and also j).
We choose n inducing input (time) points Z = (I1...., (x) with 0 <<t,- <... <,,, < 1. These pscudo inputs can be seen as representative time points for the posterior process in that knowing the values of μ4 at f Z has decisive effects inferring function values at the other input points. Further insights underlying the principle can be found in the nice survey Although the inducing input points can be learned as well as model and variational parameters in the BLBO learning, we instead fix them as equal-spaced points in [0. for simplicity (Sec. for ablation on n). Following the variational free-energy principle (Titsias 2009), we define 2(1². as:
Setthonero-Candeer
& Rasmussen 2005).
((v=)= / Q(uz) Q(uz)Q((o)))) ((zz)lz))dzz
(14)
where #z= . [4,1 IT. known the inducing variables, is a n-dimensional vector of function values at the inducing inputs2 Z. The key idea is that we model only Q(Hz) as a n-variate learnable Gaussian distribution while having Q(valuz) equal to P Pprior (Po (uz), the conditional process
4

---

Published as a conference paper at ICLR 2026
induced from the prior GP, which is also GP. This not only helps avoid modeling difficult conditional process Q(peluz), but also the integration in (14) becomes tractable leading to closed- form GP Q(e.). More specifically, define:
2(12) (uiz (H),)H)) for j l....,,
(15)
where CHi,, Sh. are n-dimensional variational parameter vectors by assuming diagonal covariance. Note that when compared to GSBM's spline parametrization, assuming the number of spline knot points is the same as here, they require (n) parameters to represent knot function and derivative values for each dimension j. we have the same complexity in terms of parametrization. And the conditional prior P Perior (Mo (uz) ) is also GP with mean and covariance functions written as:
Elpp]] M₁ Tx,2 2, (uz Mz), Cov(Uus (4) I.,,t L,, T.2. l.z.t
(16)
where M and L are the mean and covariance functions of the boundary-condittonee prior GP
(II
Due to the Gaussianity of both terms, the integration (14) can be solved analytically, yielding a Gaussian process ((.)) whose mean and covariance functions are:
Eelud] M, + L1 L, Mz),
(17)
Covo m4) L.,: - L, S L,, L
(18)
The variational posterior of the stdev function 2(5.). following the same derivations, becomes Gaussian process as well, with the following mean and covariances:
Eelid = m2 + k7 (12 2 (C+ m
(19)
Cove 7x) = k,, k²,zz z. (k ) 84 (x),2)) 1 2
12,z .2) )
(20)
where Z' is the inducing input times for n. and (O S are the variational parameters for the induc- ing variables fiz.. Overall, == ((C"J,su };; C SS constitute all our variational parameters.
ELBO learning. The ELBO objective can be derived from the negaaiivv of the KL divergence between ((P.) and (P.) (Appendix for derivations). We minimize the negative ELBÓ,
min Eorr. (P(.: (o)::.)) + (L(Q(P2,2))1PP)) prsor (P2.z.) A.J.T
(21)
with respect to the variational parameters A as well as the model parameters 17 (prior) and (likeli- hood). The variational-pria KL divergence, the second term in (T). is confined to n-dim Gaussians, and easy compute in closed forms (Appendix for for details). The first term of (1 can be computed
and optimized by reparametrized Monte-Cane sampling using (17 (Kingma & Welling 2014). For the reparametrized samples Pa. that is, ()): we apply (5) compute J(P.VV V.) where
the time derivatives (e.g., jis) are computed as (J++ - /1//Al A/ utilizing the nearby covariance structure (say, Cov (Pe: (c-A:)) and At + O. This ensures maintaining computational graphs a sole function of the model parameters (n.T), is
for backpropagation. Since the lower bounded by the ELBO where the gap exactly KL(Q(P.)|( PPppstt (Pa)). minimizing (2 with respect to A with fixed (1), -) guarantees to reduce the gap, while minimizing it with respect to (17, T) with fixed A potentially improves the (log- eeviencee. The latter optimization essentially amounts to performing empirical Bayes (evidence maximization), a principled way to do model selection. Our algorithm, dubbed GP-GSBM, summarized in Alg.
4 RELATED WORK
The core related works on Schrödinger bridge and generalized SB were discussed in Sec. and Appendix A Since our work is on application of Gaussian processes (Rasmussen & Williams 2006 to the generalized Schrödinger bridge problem, here we briefly highlight two recent papers that applied GP to the bridge/flow problem, and discuss how they are different from ours.
Stream-level flow-matching with GP (Wei & Ma 2025). This work is based on Tong et al. (2023 CFM framework whereas in our GP-GSBM the probability path GP is applied to the GSB problem.
Since uo is GP 17-178 18), f. is also GP whose mean and cov functions involve time derivatives of 17 18 But, we used approximation (Pe+ae / At with small At 0.01, which we found numerically more stable.
5

---

Published as a conference paper at ICLR 2026
Algorithm 1 Our GP-GSBM algorithm for robust generalized Schrödinger bridge matching.
Input:
The end- point distributions wo and ı (i.e., samples from them) and the state cost Vi()..
Repeat until convergence or a sufficient number times:
0. Collect hatch pairs (1) sampled from the current estimate QTa, x1) ) in
x1)
1. Solve the ELBO optimization in (21 to get variational parameters A and model parameters (7, 7).
2. Prepare batch of (20, 4,,,, .r) and (arx//20, 21) as follows:
1) Sample ~ [0, 1] uniformly at random.
2) Sample (pt, ~ ((P)). using (0) with the learned parameters A and yy for Q.
3) Sample 2; ~ (2 [4, 78 n)
4) Compute (22
in (5) using approximate time derivatives.
Update the neural net of the SDE by running the following for Nanet iterations:
0 enenn (1 (20, (20, ve (1,2
(22)
They have conditional GP derived by conditioning the prior GP on several fixed observed stream points. So this the key difference. Recall that we take the sparse variational GP posterior inference while turning the CondSOC objective into a likelihood function. A limitation of their approach is that there are some task-specific costs on the marginal paths (e.g.. obstacles specified by the stage cost V), there is way to penalize the GP path that passes through the obstacles systematically. That is, their GP framework may not be applicable to GSB problems unless one manually marks and provides the stream points circumvent the obstacles, which is cumbersome and costly.
Flow matching with GP prior (Kollovich et al. 2025). The setup and goal of this work are highly different from ours. They introduce GP as an additional prior structure on the distribution mo for generation of the time-series data Ø1_ Hence 2p is the same dimensional time-series. Then they use either Langevin dynamics (Durmus Éric Moulines, 2017) or posterior sampling DDhriwall Nichol| 2021] Kollovich al.) (202)) in denoising diffusion models (Sohl-Dickstein et al. 22115 Ho ct 20020; Song et al.] 22021) with the mini-batch-OO -a------CFF CFM Iearned from Tn- Hence this work is highly different from approach of imposing a GP prior on the pinned marginal paths.
EXPERIMENTS
We test the performance of our GP-GSBM on several crowd navigation and image generation/transle tion problems. We mainly compare ours with GSBM (Liu al. 2024) and DSBM (Shi et al. the former deals with the stage cost deterministically while the latter ignores the stage cost. The de- tails of experimental setups including hyperparameters and neural net architectures are in Appendix
(023):
C
Due to the lack of space, some parts of our experimental results are moved to Appendix D
5.1 CROWD NAVIGATION PROBLEMS
We follow the problem setups and experimental settings from (Liu et al. 2022 2024) where we focus on the geometric surface LiDAR problem (Sec. 5.1.1 and 2D obstacle problems (Sec. 5.1.2).
5.1.1 LIDAR GEOMETRIC SURFACE STAGE COSTS
We tackle the crowd navigation problem where the path cost is incurred through a complex geometric surface. Similar to the setup in (Liu et al., (024)), we deal with the surfaces observed through the LiDAR 3D seans of the Mt. Rainier (Legg & Anderson 2013). where the surfaces form a manifold. We define the closeness to the manifold and the height as the path cost. Formally,
Vx(x) = (1) x]|² + ex)(LL
(23)
where (x) is the projection of the 3D point to the manifold M, which is done by tangent plane projection estimated by z's k-nearest neighbor points (k = 20) on M. The height of the projected point is the z-coordinate value [=)
This stage cost function encourages the navigation path to stay close to the surfaces and avoid high altitudes. Fig- shows the sample paths from the competing models after training. We see that DSBM (leftmost), as expected, yields straight linear paths due to its ignorance of the stage cost, resulting in high CondSOC loss summarized in Table 1 Our GP-GSBM successfully discovers the low

---

Published as a conference paper at ICLR 2026
(Naise-free scenario)
(Noisy scenario)
DSBM
GP-GSBM (Ours) Stream-level GP
GSBM
GP-GSBM (Ours)
Figure I: LiDAR path samples for different models in 3D (top rows) and 2D (bottom) views. Starting from the samples To shown as red points on the left, we generate samples for лi as blue points on the right. Path samples visualized in 5 uniform time points from t Ł= 0 (red) t= (blue).
Table 1: LiDAR results, CondSOC objective values for the learned models. Figures in parentheses the Wasserstein distances between the true target क and the learned ones.
DSBM
GSBM
STREAM-LEVEL GP 7012.6 + 61.9 (0.15)
GP-GSBM (OURS)
NOISE-FREE OBS.
7747.0 + 76.4 (0.04)
6199.3 = 47.3 (0.04)
5925.0 + 65. (0.03)
NOISY 12686.9 + 150.9 (0.04)
8506. + 63.6 (0.04)
12679.1 + 579.9 (0.16)
8300.0 + 67.6 (0.04)
altitude regions bypassing the saddle points as shown in the second column. GSBM (Liu et al. 2024 also able to find out the desired two viable pathways similarly as ours (figure omitted here due to the similarity ours). However, our GP-GSBM achieves even lower CondSOC cost perhaps due to more ficxible pinned marginal path modeling than GSBM's deterministic spline optimization. Table summarizes the CondSOC objective values for the learned models, averaged over 10 random runs.
Noisy observation setup. Since the projection to the manifold and reading the projected point can be noisy process in practice, mimic this noisy process by injecting random noise the projected points. The results in Table show that our GP-GSBM attains the lowest loss among the competing methods, signifying its benefit of robust path estimation through the Gaussian process inference under this noisy likelihood scenario. The visualized path samples in Fig. (rightmost column) also shows that the desired two viable pathways are still well discovered by model.
5.1.2 STAGE COSTS FROM OBSTACLES AND MEAN-FIELD INTERACTIONS
We next consider the crowd navigation problems with stage costs driven by obstacles and mean-field interactions. Similar to previous works (Liu et al. 2022 2024), the stage cost is defined as:
Vx(x) Lobs (x) + Aent log (x Xegst (1 || x' 2'(13)
(24)
where (x) is the path marginal distributions obtained from the pinned marginals pr (z/ x1) and the coupling Q1 The obstacle cost Toots defined the distance between the state æ and the problem-specifie obstacle. We may have the log p. (xx) term as the entropy cost encouraging diversified paths, and the congestion term penalize densely packed particles along the paths. We test the competing models on two problem sets, St unnel that has two ellipsoidal obstacles (Fig. and CMM with three Gaussian blobs as obstacles (Fig. 5)
The learned μr. of the pinned marginals are shown in Fig. and Fig. 5 while Table 2 reports the final CondSOC loss values averaged over 10 random runs. DSBM always finds a straight line, and incurs high CondSOC loss due to its ignorance of obstacles, i.e., the stage cost V. For Stunnel (Fig- 4.. in the deterministic obstacle scenario, GSBM identifies an optimal path that circumvents the obstacles sharply while retaining minimal kinetic energy. Our GP-GSBM also finds similar solution as GSBM's, and the posterior stdevs (red-shaded areas) are small in this case due to the deterministic nature of the problem. The model selection empirical Bayes finds small (=s 0.1), implying that the model that puts more emphasis on the likelihood than the prior. This is reasonable and promising considering that the likelihood in this scenario is quite trustworthy.
Unlike in the original setup where the obstacles were always present, and the bridge matching algorithms can learn the stage cost function deterninistically, devise uncertain setup by
7

---

Published as a conference paper at ICLR 2026
Table 2: Stunnel and GMM results. CondSOC objective values for the learned models. Figures in parentheses are the Wasserstein distances between the true target Ti and the learned ones.
DSBM
GSBM
GP-GSBM (OuRS)
DETERMINISTIC
18628.8 (0.02)
492.94 = 1.13 (0.03)
488.78 + L.07 (0.02)
STUNNEL
UNCERTAIN
9549.2 (0.03)
502.20 + 1.76 (0.04)
452.30 + 1.78 (0.03)
DETERMINISTIC
19824.2 (2.91)
97.4 + 0.67 (2.31)
85.3 1 0.54 (1.81)
GMM
UNCERTAIN
13232.4 (2.96)
101.6 + 1.02 (2.83)
89.2 + 0.62 (2.19)
regarding the presence of obstacles random process. Specifically we randomly turn on and off the obstacles with probability p p.0.5. This has the effect of turning the CondSoc objective function into an uncertain quantity, and we aim to judge which models are the most robust such noisy CondSOC objective function. We anticipate that our GP-GSBM will be more robust to GSBM (Liu et al. 2024) since we treat the CondSOC loss a likelihood function, that is, probabilistic model, rather than GSBM's deterministic treatment. As shown in Table 2 GSBM incurs even higher CondSOC loss than the deterministic case. GSBM is baffled in this case because it only optimizes the cost without taking into account the prior preference of straight line interpolation between To and J1, which is preferred for the obstacle absent cases. On the other hand, our GP-GSBM learns higher prior weight (R= 1.0) from model selection, which encourages the marginal path to stay closer to the straight line interpolation (prior), preferable for obstacle absent cases. Furthermore, the uncertain stage cost is reflected in the stdev shading in our GP posterior where we observe even higher posterior uncertainty shown as larger red-shaded areas in Fig. | (the rightmost column).
5.2 UNPAIRED CAAGE-TO-IMAN TRANSLATION
In this section we test the performance of our GP-GSBM on unsupervised image-to-image translation problem. Specifically, we deal with images of cats (regarded as πо) and dogs (as πΙ) from the AFHQ dataset (Chol et al.; 2020). Although SB's s optimal transport objective guides a bridge matching algorithm to learn coupling between samples from सo and कi that is optimal in the ambient pixel- space L.2 distance sense, this can often incur artifacts in the bridging paths. In other words, more intrinsic (latent) structure of the image manifolds cannot be taken into account.
To this end, following the latent space guidance proposed in (Liu et al. (2024), we incorporate the latent spherical-linear interpolated reconstruction error as the stage cost. More specifically,
Vi (zx) || dec(2e)ll1:
24 slerp(t, enc(xo), enc(x))
(25)
where enc(x) and x = dec(2) are the encoder and decoder of the pre-trained VAE model (Kingma Welling 2014). ) As claimed in (Liu et al. (0024), this stage cost helps preserving the underlying latent structures through semantically meaningful pinned marginal P. (2x (r, ). yielding faster training convergence and better couplings than models without the V term such as DSBM (Shi et al.,
2023). In addition to this observation, we further conjecture that the cost V itself is inherently noisy for several reasons, most notably that the learned VAE model may not be perfect in representing the true image manifold, and also the SLERP interpolation is only a proxy for the optimal latent paths. This is where our GP-GSBM can be especially beneficial for more robust path learning and bridging by handling such uncertainty in the stage cost in a principled Bayesian manner.
We trained our GP-GSBM with images of size (64 x 64) from cats and dogs, roughly 5K images from each. We first visualize how the prior and posterior pinned marginals (xr x1) differ for dog (t= 1) to cat 0)= translation. After training the model, we pair the coupling of (81.20) by running the trained (reverse-time) SDE model starting from 21 ~ π1, which is shown in Fig. 2(aa) (more and magnified images in Fig. 6). From this coupling, we show the mean images from Pi x1) for the prior and posterior (before and after the variational GP inference) in the middle. Visually clearly, the posterior exhibits semantically more meaningful path samples than the prior with less artifacts since the posterior takes into account the latent- preeett-reeeving stage cost, forcing the samples to stay closer to the image manifold. In Fig. 2 b) (also Fig. 7. the SDE generation progress is shown with intermediate time steps. We also visualize the generated cat images in Fig. ( (also Fig. 8). As a quantitative comparison, the FID scores of the generated cat images reported in Table Our GP-GSBM attains the lowest FID score, better than (deterministic) GSBM by capturing potential uncertainty/inaccuraay that may reside in the stage cost observation.
8

---

Published as a conference paper at ICLR 2026
I-0
Prior
Post
Prior
Post
Prior
Post
Prior
Post
(a) Prior and posterior of | x1)
(b) Generation by SDE
(c) Generated images
Figure 2: Results of GP-GSBM on AFHQ dog (t = 1, rightmost columns in (a), (b)) + cat (t
leftmost columns in (a), (b)) translation. See Appendix D.I for more and enlarged images.
Table 3: AFHQ dog + cat image generation FID scores. Figures for DSBM (Shi al. 2023 and
GSBM (Liu et al. 2024) are excerpted from (Liu et al. 2024).
DSBM
GSBM
STREAM-LEVEL GP
GP-GSBM (OURS)
14.16
12.39
18.77
10.21
5.3 COMPARISON TO STREAM-LEVEL GP
As discussed in Sec. 4 our approach has several key differences from the conditional GP path modeling for CFM models in the recent work (Wei & Ma 2025). Most notably, we regard GSBM's CondSOC objective as likelihood function for GP posterior inference, whereas they used the linear interpolation velocity ay (x Zo; Ø1) = Σ. as a conditional GP prior, thus difficult to incorporate task-specific stage costs. Here we perform some empirical tests highlight the differences. We compare ours with the stream-level GP method on the LiDAAR crowd navigation (Fig. [ Table and AFHQ image-to-image translation (Table 3. We adapted the stream-level GP within our UBA framework for these bridge matching problems. Overall the stream-level GP incurs high costs with lower distribution matching performance due the ignorance of the stage costs. For additional results on non-GSB problem, refer to Appendix D.2 for CIFAR-10 10 image generation experiments.
5.4 ABLATION STUDY AND RUNNING TIME (COMPUTATIONAL COMPLEXITY)
Sensitivity the number of induc- ing points (n) and the kernel func- tion. We used the default n = 15 (Stunnel) and = 30 (LiDAR). To see how sensitive our GP-GSBM is
Stunnel
UDAR
mintstic Uncertain
Conssooi loss
700
ComtiiCK tuae
8000
Noise free Noiy
600
2000
500
to this hyperparameter, we vary in Fig. 3 which shows that the per- formance is not very sensitive to n unless it is too small. We also test kernel types other than the default squared exponential kernel. As shown in The polynomial kernel (degree selected from grid of to 10) slightly lags behind it.
6000
20
10
200
-0 40
50
Number ontindc polnns (m)
Number of Incucting points (n)
Figure 3: Sensitivity to the number of inducing points (n).
Running time. In Appendix B.2 we analyzed asymptotic complexity of GP-GSBM compared to GSBM. Due to the GP modeling, our algorithm has extra cost of O(n²d) (also (n kernel inversion cost when we do model selection). Even though n (the number of inducing points) is constant and usually no greater than 30 in practice, this may be computationally demanding especially for high- dimensional states. However, we suggested some workaround to reduce this cost in Appendix B.2 In this case, on LiDAR we got comparable time, GP-GSBM: 1.70 second per ELBO iteration, and
Table 4: CondSOC objectives for different kernels on Stunnel (left) and LiDAR (right).
STUNNEL
SQUARED EXP.
POLYNOMIAL
LIDAR
SQUARED EXP.
POLYNOMIAL
DETERMINISTIC
485.78 = 1.07
556.50 + 1.47
NOISE-FRRF NoisY
5925.0 + 65.4
6604.8 1: 55.4
UNCERTAIN
452.30 = 1.78
435.37 + 9.31
$300.0 + 67.6
8625.2 :: 65.8

---

Published as a conference paper at ICLR 2026
GSBM: 1.61 second per CondSOC iteration, when run single RTX-4090 GPU. We conduct more theoretical and empirical study to reduce further computational overhead as ongoing/future research.
CONCLUSION AND LIMITATION
We have proposed novel Gaussian process approach to marginal path modeling for the generalized Schrödinger bridge problem. By imposing a GP prior on the pinned marginal path and viewing the CondSOC objective can lead to more flexible and robust solutions than the existing methods on problems under noisy observations and uncertainty. On several crowd navigation and amage-to-imane translation problems we have empirically demonstrated these benefits. A limitation of our current approach, however, that it incurs higher computational complexity arising from GP posterior inference involving kernel inversion. Although we suggested some heuristic workarounds, further study on finding more principled solutions to reduce time complexity is needed, which leave as future research.
noisy likelihood function, the inferred posterior path in our model
REFERENCES
Michael Albergo and Eric Vanden-Fijnden. Building Normalizing Flows with Stochastic Inter- polants, 2023. International Conference on Learning Representations.
Matthias Bauer, Mark van der Wilk, and Carl Edward Rasmussen. Understanding probabilistic sparse Gaussian process approximations, 2016. In Advances in Neural Information Processing Systems.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations, 2018. In Advances in Neural Information Processing Systems.
Yongxin Chen. Density control interacting agent systems. IEEE Transactions on Automatic Control, 69(1))24 260, 2023.
Yongxin Chen, Tryphon Georgiou, and Michele Pavon. Optimal steering of inertial particles diffusing anisotropically with losses, 2015. American Control Conference (ACC).
Yunjey Choi, Youngjung Uh, Jaejun Yoo, and Jung-Woo Ha. Stargan v2: Diverse image synthesis for multiple domains, 2020. IEEE Conference on Computer Vision and Pattern Recognition (CVPR).
V. De Bortoli, J. Thornton, J. Heng, and A. Doucet. Diffusion Schrödinger Bridge with Applications to Score-based Generative Modeling, 2021. In Advances in Neural Information Processing Systems.
A. Dezfouli and V. Bonilla. Scalable inference for Gaussian process models with black-box likelihoods, 2015. In Advances in Neural Information Processing Systems.
Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans image synthesis, 2021. In Advances in Neural Information Processing Systems.
Alain Durmus and éric Moulines. Nonasymptotic convergence analysis for the unadjusted langevin algorithm. The Annals of Applied Probability, 27(3):155-1187,,
2017. 2017.
R. Fortet. Résolution d'un système d' 'ééquations de M. Schrödinger. Journal de Mathématiques Pures el Appliqués, 1:83-105, 1940.
Jason Gaitonde, Jon Kleinberg, and Éva Tardos. Polarization in Geometric Opinion Dynamics, 2021.
ACM Conference on Economics and Computation.
Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. In Advances in Neural Information Processing Systems, 2018.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising Diffusion Probabilistic Models, 2020. In Advances in Neural Information Processing Systems.
Hilbert Kappen. Path integrals and symmetry breaking for optimal control theory. Journal of Statistical Mechanics: Theory and Experiment, 11(P11011)1 1-:1244, 2005.
10

---

Published as a conference paper at ICLR 2026
Beomsu Kim, Gihyun Kwon, Kwanyoung Kim, and Jong Chul Ye. Unpaired Image-to-Iman Transla- tion via Neural Schrödinger Bridge, 2024, International Conference on Learning Representations.
Minyoung Kim. A Unified Framework for Diffusion Bridge Problems: Flow Matching and Schrödinger Matching into One. In arXiv preprint, 2025. URR/ttts:://are http://aaxxx...r/abb...o/ owg/ 15533217566
Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization, 2015. International Conference on Learning Representations.
Diederik P. Kingma and Max Welling. Auto-encoding variational Bayes, 2014. In Proceedings of the Second International Conference Learning Representations, ICLR.
Marcel Kollovich, Abdul Fatir Ansari, Michael ScinlkeeSSheeider, Jasper Zschiegner, Hao Wang. and Yuyang Bernie Wang. Predict, refine, synthesize: Self-guiding diffusion models for probabilistic time-series forecasting, 2023. In Advances in Neural Information Processing Systems.
Marcel Kollovieh, Marten Lienen, David Lüdke, Leo Schwinn, and Stephan Günnemann. Flow Match- ing with Gaussian Process Priors for Probabilistic Time Series Forecasting, 2025, International Conference on Learning Representations.
S. Kullback. Probability densities with given marginals. The Annals of Mathematical Statistics, 39 (4):1236-1243, 1968.
Nicholas Legg and Scott Anderson. Southwest Flank of Mt.Rainier, WA. 2013. URL https
pyppenttporraplyy.og//eettttttte 0wtp:/.3....2011...... 269010
Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. 2023. International Conference on Learning Representations.
Guan-Horng Liu, Tianrong Chen, Oswin So, and Evangelos A Theodorou. Deep Generalized Schrödinger Bridge, 2022. In Advances in Neural Information Processing Systems.
Guan-Horng Liu, Yaron Lipman, Maximilian Nickel, Brian Karrer, Evangelos A. Theodorou, and Ricky T. Q Chen. Generalized Schrödinger Bridge Matching, 2024, International Conference on Learning Representations.
Xingchao Liu, Chengyue Gong, and Qiang Liu. Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow, 2023. International Conference on Learning Representations.
Zhiyu Liu, Bo Wu, and Hai Lin. A mean field game approach to swarming robots control, 201 8.
American Control Conference (ACC).
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization, 2019. International Conference Learning Representations.
A. Matthews, Hensman, R. E. Turner, and Z. Ghahramani. On sparse variational methods and Kullback-Leibler divergence between stochastic processes, 2016. International Conference on Artificial Intelligence and Statistics.
Kirill Neklyudov, Rob Brekelmans, Daniel Severo, and Alireza Makhzani. Action Matching: Learning Stochastic Dynamics from Samples, 2023. International Conference on Machine Learning.
Frank Noé, Alexandre Tkatchenko, Klaus-Robert Müller, and Cecilia Clementi. Machine learning for molecular simulation. Annual Review of Physical Chemistry, 71:361-390, 2020.
Gabriel Peyré and Marco Cuturi. Computational optimal transport, 2017. Center for Research in Economics and Statistics Working Papers.
C. Philippidis, C. Dewdney, and B. J. Hiley. Quantum interference and the quantum potential. Nuovo Cimento B, 52(1):15-28, 1979.
Joaquin Quiñonero Cinbonerr-Candela and Carl Edward Rasmussen. A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research, 6:1939-1959, 2005.
11

---

Published as a conference paper at ICLR 2026
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2006.
L. C. G. Rogers and D. Williams. Diffiusions, Markov processes, and Martingales. Vol. 2. Cambridge University Press, Cambridge, 2000.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-Net: Convolutional networks for biomedical image segmentation, 2015. International Conference on Medical Image Computing and Computer- assisted Intervention.
L. Rüschendorf. Convergence of the iterative proportional fitting procedure. The Annals of Statistics, 23(4): 1604)1160-1174, 1995.
Yuyang Shi, Valentin De Bortoli, Andrew Campbell, and Arnaud Doucet. Diffusion Schrödinger Bridge Matching, 2023. In Advances in Neural Information Processing Systems.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics, 2015. International Conference on Machine Learning.
Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising Diffusion Implicit Models, 2021.
International Conference on Learning Representations.
M. K. Titsias. Variational learning of inducing variables in sparse Gaussian processes, 2009, In Proceedings of the Twelfth International Conference on Artificial Intelligence and Statisties.
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. In arXiv preprint, 2023. URL hatptttppss... /ttpps://arzaa /ttp///aavvv//////a//eee/eaaaa./ org/bbs/ 2302. 00482]
B. Tzen and M. Raginsky. Theoretical guarantees for sampling and inference in generative models with latent diffusions, 2019. In Conference on Learning Theory.
F. Vargas, P. Thodoroff, A. Lamacraft, and N. Lawrence. Solving Schrödinger Bridges via Maximum Likelihood. Entropy, (3(9):1134, 2021.
Cédric Villani. Optimal Transport: Old and New. Springer Berlin, Heidelberg, 2009.
Gefei Wang, Yuling Jiao, Qian Xu, Yang Wang, and Can Yang. Deep Generative Learning via Schrödinger Bridge, 2021. International Conference Machine Learning.
Ganchao Wei and Li Ma. Stream-level Flow Matching with Gaussian Processes, 2025. International Conference on Machine Learning.
12

---

Published as a conference paper at ICLR 2026
Appendix
Table of Contents
A Unified Framework for Diffusion Bridge Problems (UBA) (Sec. A
Technical Details of Gaussian Process GSBM (GP-GSBM) (Sec.
Experimental Details (Sec.
Additional Experimental Results (Sec.
D
Revision: Thought Experiment and Background Gaussian Processes (Sec.
A A UNIFIED FRAMEWORK FOR DIFFUSION BRIDGE PROBLEMS: FLOW MATCHING AND SCHRÓDINGER MATCHING INTO ONE
The bridge problem is to find a stochastic differential equation (SDE), or sometimes an ordinary differential equation (ODE), that bridges two given distributions. The application areas of the bridge problem are enormous, among which the recent generative modeling (e.g., conditional or unconditional image generation) the most popular. Also the famous Schrödinger bridge problem, widely known problem for century, is special instance of the bridge problem. Two most popular algorithms to tackle the bridge problems in the deep learning era are: (conditional) flow matching and iterative fitling al algorithms, where the former confined to ODE solutions, and the latter specifically for the Schrödinger bridge problem. The main contribution this section is in two folds: i) We provide concise reviews of these algorithms with technical details to some extent: ii) We propose a novel unified perspective and framework that subsumes these seemingly unrelated algorithms instantiate the
(and their variants) into one. In particular, we show that unified framework Flow Matching (FM) algorithm, the (mini-batch) optimal transport FM algorithm, the (mini batch) Schrödinger bridge FM algorithm, and the deep Schrödinger bridge matching (DSBM) algorithm as its special cases. We believe that this unified framework will be useful for viewing the bridge problems in a more general and flexible perspective, and in turn can help researchers and practitioners to develop new bridge algorithms in their fields.
A.1 (DIFFUSION) BRIDGE PROBLEMS
The diffusion bridge problem, or simply the bridge problem, can defined as follows.
Bridge problem. Given two distributions ( and (-) in R4, find SDE, more specifically, find the drift function uz(z) where u R(0, 1] x Rd + Rd with specified diffusion coefficient a. dæ, ui (xt) + adW,, zo To~))o))
(26)
that yields in m1(+). can aim to find a reverse-time SDE (or both). That is, find (x(r) in
(-). Here (W:): h is the Wiener process the Brownian motion. Alternatively one
غ dz (xt)dtt + o a Wi, x))m((+
(27)
that yields Fo m.)+).)
Note that if we specify o a= =0, 0, then our goal is to find an ODE that bridges the two distributions To(-) and Ti (-). Once solved, the solution to the bridge problem can give us the ability to sample from one of the (T(01) given the samples from the other, simply by integrating the learned SDE. The application areas the bridge problem are enormous, among which the generative modeling (e.g.. conditional or unconditional image generation) is the most popular. For instance, in typical generative modeling, To is usually a tractable density like Gaussian, while m is a target distribution that we sample from. In the bridge problem, however, mo also be an arbitrary distribution beyond tractable densities like Gaussians, and we do not make any particular assumption on Ay and na as long have samples from the two distributions.
Two instances of the bridge problem. There are two interesting special instances of the bridge problem: the Schrödinger bridge problem and the ODE bridge problem.
13

---

Published as a conference paper at ICLR 2026
ODE bridge problem. We strictly restrict ourselves ODEs, i.e., o
Schrödinger bridge problem. With a > 0, there is an additional constraint that the path measure of the SDE 26). denoted by p", is closest to given reference SDE path measure pref, That is,
min (L(P))) s.t. (rro) no (zo)) Pf m (21)
(28)
where Pe of a path measure P indicates the marginal distribution at time t. In this case to have finite KL divergence, a of pd has to be set equal to Fref of pref, i.e. Oref, (from the Girsanov theorem).
A.2 FLOW MATCHING AND SCHRÓDINGER BRIDGE MATCHING ALGORIHTM
Among several existing algorithms that aim to solve the bridge problems, here we focus on two recent
matching algorithms: (conditional) flow matching (Lipman et al. 2023] Tong et al. 2023 Albergo
Vanden-Eijnn Eijnden 2023 Liu et al. 2023) and Schrödinger bridge matching algorithms (De Bortoli
et al. 2021 Vargas et al. | Shi et al., 2023). These algorithms were developed independently:
the former aimed to solve the ODE bridge problem while the latter the Schrödinger bridge problem.
In this section we review the algorithms focusing mainly on the key ideas with some technical details,
but being mathematically less rigorous for better readability.
In Sec. A.3
we propose a novel unified framework that subsumes these two seemingly unrelated algorithms and their variants into one.
A.2.1 (CONDITIONAL) FLOW MATCHING FOR ODE BRIDGE PROBLEMS
The Flow Matching (FM) ( Lipman et al. 2023 or its extension Conditional Flow Matching (CFM) (Tong et al., 2023) is one promising way to solve the ODE bridge problem. The key idea of the FM is quite intuitive. We first design some marginal distribution path (Zt() with the boundary conditions Po = To. Pi = TT, We then derive the ODE dæ, = (2 dt that yields (xx)) its marginal distributions. The drift (xx)) is approximated by neural network (t, xx) with parameters by solving:
min (x ue (t, xx)l
(29)
Once solved, we can generate samples from Ti (or To) approximately by simulating: dx, (t, 2z dt) To ~ To (resp., that designing the marginal path (xx(): satisfying the boundary condition often difficult. And is this issue that motivated the CFM.
T1). However, one of the main limitations of this strategy is
Conditional Flow Matching (CFM). To make the path design easier, we introduce some latent random variable to condition It. Although CFM derivations hold regardless of the choice of it typically chosen as the terminal random variates z = (Io, I1), and we will follow this practice and notation. Specifically, in CFM we design the so-called pinned marginal path (P(xx (0o, 21) he and the coupling distribution Q( .1) subject to the condition (o = To(Io). Pi (11) = Ti (11) where (I())) is defined as
P((a))= / Pa( on an) Q on )d((ro 21)
(30)
Then we derive the ODE dx, un (xx (zo 21 dtt that yields (2 x1)h as its marginal distribu- tions for each (zo, æ1), which admits closed form if [P(z) x1)} are Gaussians (Lipman et al. 2023 We then approximate Elu (r. x1 ) [ar), the conditional expectation derived from the joint P( Zo, 21 Q 21), by neural network vo (t, xx) by solving the following optimization:
min E [u, To. 21) - vo (t,t,
(31)
where the expectation is taken with respect to the joint P( Jao X1 Q((u,, .11 and uniform t. Sur- prisingly, it can be shown (Tong et al. 2023) that the gradient of the objective in 31 coincides with that in (29 for (xx)) defined
ue(x)=
E (ul x1 ) Pa |20, #1)]
(32)
P((x))
14

---

Published as a conference paper at ICLR 2026
hence sharing the same training dynamics 29 Therefore the optimal (t,t, of 31 is good estimate for (zz). Since the ODE dz, = uy ( admits (2t)) as its marginal distributions, so does dr, (o, 2,) approximately.
Many existing flow matching variants including FM (Lipman et al. 2023). Stochastic Interpola- tion (Albergo & Vanden- |ijdenn 2023), and Rectified FM (Liu et al. 2023) be viewed special instances of this CFM framework. For instance, these models can be realized by having a straight line (linear interpolation) pinned marginal path (33) with vanishing variance while one boundary (e.g..
TO) is fixed as standard normal (0,
Limitations of CFM. A reasonable choice for the coupling distribution Q z1) is the Optimal Transport (OT) or the entropic OT between To and πT. The pinned marginal Pt x1) can be chosen Gaussian with the linear interpolation between Zo and ø its mean, more specifically,
x1)
(r x1) [ttz + (1 3
(33)
for some scheduled variances S2. When the combination of the entropic OT Q x1) and S = Fref /eefv(111/ is used, it can shown that the marginals {Pr(ar)) coincide with the marginals of the Schrödinger bridge with the Brownian motion reference pref : dx, = d.W However, the main limitations of CFM (Tong et al. 2023) are: CFM solutions confined to ODEs, hence unable to find the optimal SDE solution to general bridge problems including the Schrödinger bridge problem; ii) CFM itself does not provide recipe about how to solve the entropic OT problem exactly what is called SB-CFM proposed in (Tong et al. (2023) only approximates it with the Sinkhorn-Knopp solution for minibatch data, which is usually substantially different from the population entropic OT solution, that is, the solution to the Schrödinger bridge problem (See the SB static view in Sec. A.2.2).
A.2.2 SCHRÓDINGER BRIDGE PROBLEM
The Schrödinger bridge problem can be defined as (28) where we assume a zero-drift Brownian SDE with diffusion coefficient Fref for the reference path measure pref. That is,
pref dx, = drej Io mo-))(()
(34)
We denote by PsB the Schrödinger bridge path measure, i.e., the solution to 28 In the literature, there are two well-known views for PsB: the static view and the optimal control view.
The static view has a direct link to the entropic optimal transport (EOT) solution, more specifically PSB p hejo..) = PEOT x1) pref (( (
(35)
where PEOT (1) is the EOT joint distribution solution with the negative entropy regularizing coefficient 20 More formally,
PFOT (wo, PROT(HO, x1) = arg min Ep(2a.2)) - x1112 (P((o,, xI)) p(
(36)
s.t. P(xo)= = =o(oo) P(x1) = TT (I1))
(37)
where H. indicates the Shannon entropy. We call pref (11 the pinned reference process, which admits a closed-form Gaussian expression for the specific choice (34). Although the product form (35), i.e., the product of the boundary joint distribution and the pinned path measure, does not in general become Markovian (e.g., Itō SDE representable), the Schrödinger bridge a well-known exception where there exists unique SDE that yields psB as its path measure.
Alternatively, it is not difficult to derive an optimal control formulation for the Schrödinger bridge problem. Specifically, PSB can be described by the SDE that has the minimum kinetic energy among those that satisfy the bridge constraint. Letting
pd da. = ve( )dt + GrefdW,, xo - mo(-)
(38)
we have PSB_ p where o² is the minimizer of the following problem:
(x:))) min Epe S 20
s.t. Pf((1) 7 (21)
(39)
Next we summarize two recent algorithms that solve the Schrödinger bridge problem exactly (at least in theory): Iterative Proportional Filtering (IPF) and Iterative Markovian Fitting (IMF).
15

---

Published as a conference paper at ICLR 2026
A.2.3 ITERATIVE PROPORTIONAL FILTERING (IPF)
IPF aims to solve the Schrödinger bridge problem (28) by alternating the forward and reverse half bridge (HB) problems until convergence. More specifically, with initial pa = -rre,, we solve the followings for 1......
(Reverse HB) p2n-1= = arg min (L(PIIP22-2) s.t. Pi(x1) = = (21)
(40)
P
(Forward HB) p2n arg min (L(P||P22~)) s.t. Po(xu) = To (20)
(41)
It can be shown that lim pm + PSB (oortet 1940 Kullback 1968 Rüschendorf 1995). It is
not difficult to show that the optimal solution of (40) 41) can be attained by time-reversing the
SDE of the previous iteration. This fact was exploited recently in (De Bortoli et al. 2021 Vargas al. 2021) yield neuwoll based IPF algorithms where the score V log P" ( that appears in time reversal is estimated either by regression estimation (De Bortoli et al. 2021 or maximum
likelihood estimation (Vargas et al., (202T) However, the main drawback these IPF algorithms is that they are iimulation-based methods, thus very expensive to train.
A.2.4 ITERATIVE MARKOVIAN FITTING (IMF)
Recently in (Shi et al. 2023), the concept of path measure projection was introduced, specifically the Markovian and reciprocal projections that preserve the boundary marginals of the path measure. This idea was developed into novel matching algorithm called the iterative Markovian fitting (IMF) that alternates applying the two projections starting from the initial path measure. Not only it shown to converge to PSB but the algorithm is computationally more efficient than IPF without relying on simulation-banl learning. A practical version of the algorithm dubbed Deep Schrödinger Bridge Matching (DSBM).
We begin with discussing the two projections.
Reciprocal projection. They define the reciprocal class of path measures to be the set of path measures that admit pref (-|rr, x1) as their pinned conditional path measures. That is, the reciprocal class R is defined as:
R = (P: P(ao, P( pres (-Iro,
(42)
The reciprocal projection of path measure P. denoted by 1l (P). is defined the path measure in the reciprocal class that is closest to P in the KL divergence sense. Formally,
1R (P) = arg min (L(P|[R) = 2( 2 (-|z0, 21) RER
(43)
where the latter equality can be casily derived from the KL decomposition property. So it basically
says that the reciprocal projection of P is simply done by replacing P(. by that of p
while keeping the coupling (x0,
Markovian projection. They also define the Markovian class the set of any SDE-representable path measures with diffusion coefficient Oref. That is,
= (P: dr. gr (x, )dt + drejdW, for any vector field g)
(44)
The Markovian projection of a path measure P, denoted by IlM(P), is defined similarly as the path measure in the Markovian class that is closest P in the KL divergence sense,
[P = arg min KL(PI|M) MEM
(45)
In Shi et al. 2023) (Proposition 2 therein), it was shown that (PM(P) can be expressed succinetly for reciprocal path measures P. Specifically, for P R. we have IIM(P) = p=" where p= is described by the SDE: dr, = oz (zt)dt
Bre,dW,, Zy P(xg) where
of (sx) = Ep(nin m.)) lofes Va, log pref (21 = Ep( laxs)
T1 - D4
(46)
1-
where the latter equality comes immediately from the closed-form pref (2:1))) = -N oee (1
T). Also it shown that the marginals are preserved after the projection, that is, P ()) = Pr(.)
16

---

Published as a conference paper at ICLR 2026
for all t E [0, 1]. This means that applying any number of Markovian (and also reciprocal) projections to path measure P always preserves the boundary marginals Po(2r) and Pi (z1). And this is one of the key theoretical underpinnings of their algorithms called IMF and its practical version DSBM (details below) to solve the Schrödinger bridge problem.
IMF and DSBM algorithms. Conceptually the IMF algorithm can be seen as a successive alternat- ing application of the Markovian and reciprocal projections, starting from any initial path measure pe that satisfies pa (xzo) = no (2)) and pd (x1) - т (xi) (e.g., pD mo(20)) J (x)))) (re)) (- (-(ro, 21)
is a typical choice). That is, for n = 1,2
p2n-1= = (
p2n I (I
(47)
Not only do all {P*) Jn20 meet the boundary conditions (i.e., Pe To. PT = it can be also shown that they keep getting closer to pSB, and converge to PSB (i.e., KL(PP+1||| PSBB < [L(P")||SB) and 1im pm pSB) (Shi al. 2023) (Proposition and Theorem 8 therein). The reciprocal projection is straightforward as it only requires sampling from the pinned process (re/() 21) that is done by running dr, = 4=7 di + JrejdW, with (#o, Xi) taken from the previous path measure. However, the Markovian projection involves the difficult P(z1 (x2) in (46) from the previous path measure P. To circumvent P(r (aze), they used the regression theorem by introducing a neural network (t, x) to approximate ef (x) and optimizing the following:
Jn20
L'" (/ 21) oies Vx, log pref (z) (r]
arg min
(48)
where now the cached samples (z1, 24) from the previous path measure P can be used to solve (48). Although theoretically ug- (t, (z) with ideally rich neural network capacity and perfect optimization, in practice due to the neural network approximation error, the boundary condition is not satisfied, i.e., pP--1 j Ti. Hence to mitigate the issue, they proposed IMF's practical version, called the Diffusion Schrodinger Bridge Matching (DSBM) algorithm (Shi et al. 2023). The idea to do Markovian projections with both forward and reverse-time SDEs in an alternating fashion where the former starts from T and the latter from TT, which was shown to mitigate the boundary condition issue.
A.3 A UNIFIED FRAMEWORK FOR DIFFUSION BRIDGE MATCHING PROBLEMS
Our proposed unified framework is described in Alg. It can be seen as an extension of the CFM algorithm (Togg et al. 2023) where the difference is that we consider the SDE bridge instead of the ODE bridge (i.e.. the diffusion term in step 2). But this difference is crucial, as will be shown, allowing to resolve the limitations of the CFM discussed in Sec. A It also makes the framework general enough to subsume the IMF/DSBM algorithm for the Schrödinger bridge problem and various ODE bridge algorithms as special cases. We also emphasize that even though this small change of adding the diffusion term in step 2 may look minor, its theoretical consequence, specifically our theoretical result in Theorem A. has rarely been studied in the literature by far.
We call the unified framework Unified Bridge Algorithm (UBA for short). Note that UBA described in Alg. can deal with both ODE and SDE bridge problems, and if the diffusion coefficient vanishes, it reduces to CFM for ODE bridge. Similarly as CFM, under the assumption of rich enough neural network functional capacity and perfect optimization solutions, our framework guarantees to solve the bridge problem. More formally, we have the following theorem.
Theorem A.I (Our Unified Bridge Algorithm (UBA) solves the bridge problem). If the neural network uo (t. x) functional space is rich enough to approximate any function arbitrarily closely, and if the optimization in step 3 can be solved perfectly, then each iteration of going through steps 1-3 in Alg. 2 ensures that dx, = volt, Zz)dt + odW,, xo mo(-) (after the optimization in step 3) admits {PA (. of (49) as its marginal distributions.
The proof can be found in Sec. A.4 The theorem says that after each iteration of going through steps 1-3, it is always guaranteed that the current SDE admits (Pr defined in step as marginal distributions. Since Po(-) = =-(-) and (.(.) < (-).), the bridge problem is solved. Depending on the design choice, one can have just one iteration to solve the bridge problem. Under certain choices,
however, it might be necessary run the iterations many times find the desired bridge solutions
17

---

Published as a conference paper at ICLR 2026
Algorithm 2 Our Unified Bridge Algorithm (UBA) for bridge problems.
Input:
The end- end-point distributions To and (i.e., samples from them). Repeat until convergence or a sufficient number times:
Choose pinned marginal path { (t (ar x1)): and coupling distribution Q( 21) such that
P-(-)= = +-)+) and Pi(-) = (-)-) where
(r(x,) f (i
(49)
2. Choose σ 2 0, and find un (n| 21) such that the SDE
dx, (a )dl + adW,
(50)
admits {Pi (x) xi) } as its marginals. (Note: many possible choices for a and un (2/2o,, 21))
3. Solve the following optimization problem with respect the neural network es (l)
min E, (P.teoor))) Ie (2, z1) (t x+)l
(51)
Return: The learned SDE dz, = (t, x: dt odW, as the bridge problem solution.
(e.g., mini-batch OT-CFM Tong et al. 2023) and the IMF/DSBM Schrödinger bridge matching
algorithm (Shi et al. 2023 illustrate in Sec. I and Sec. A.3.3 respectively).
In the subsequent sections, we illustrate how several popular ODE bridge and Schrödinger bridge algorithms can be instantiated as special cases of our UBA framework.
A.3.1 A SPECIAL CASE: (MINI-BATCH) OPTIMAL TRANSPORT CFM TTNG ET AL,
2023
Within our general Unified Bridge Algorithm (UBA) framework (Alg. 2). we select (2 Fo, Q .1) and u,(x (I) ai) as follows. First in step 1,
11),
Pe x1) = (2x;;; (1 - t)xo + t mmin
(52)
Q x1) pm00 (I), := Σ Σ
8(z0 =0))(m1 = xr)ppgo mOT
(53)
BE JEB,
where mmin + and ((=ihe Bu (e))) is the mini-batch data, and {ppor Jay is a (l (Bi))
mini-batch OT solution matrix learned with the mini-batch as training data. That is,
1
p mor = arg min Σ Pullas - all S.L. 10
Σ JEB:
Σ Paj B
(54)
B
B1
It is worth mentioning that Onain + 0 is required have boundary consistency for Po ( x1)
t -= 0 and = 1. Note also that (r ) is always fixed over iterations while Q x1) varies iterations depending on the mini-batch data sampled. Note that in our UBA framework, each iteration allows for different choices of P. (z) Zo. æ1) and Q 21).
In step 2, we choose a and define u, (2x (2o, 21) be constant (independent on straight line vector from To to æ i.c.,
U2 (x; Mo, F1 ) = Jj To which can be shown to make the ODE dx, = ue( [2o x1)dt admit Pe(xe x1) as its marginal distributions (Tong et al. 2023
(55)
The above choices precisely yield the mnini-batch) optimal transport CFM (OT-CFM) introduced in (Tong al. 2023).
A.3.2 A SPECIAL CASE: (MINI-BATCH) SCHRÖDINGER BRIDGE CFM (TONG ET AL. 2023
Within our general Unified Bridge Algorithm (UBA) framework (Alg 2). we select Px ( #1). Q(To, 21) and ux(r) (ro 21) as follows. First step 1,
(z .1) pr (2 (ate xI) Nar;; (1 t)zo + tx 1 (11)) f)T)
(56)
18

---

Published as a conference paper at ICLR 2026
Q(wo, -mmpo (ro, = Σ E 8(m0 zz))(x1 =1)p FOT (57)
iEBo
where (aab)ieBo {aT))en) ) is the mini-batch data, and (p") )is is a BBl x |Bil) is the mini- batch entropic OT solution matrix with the negative entropy regularizing coefficient 20 re (e.g., from the Sinkhorn-Knope algorithm) learned with the mini-batch as training data. Although the samples from the coupling distribution Q( 2: in (57) over the iterations in Alg. conform to the data distributions To and 1 marginally, the mini-batch entropic OT solution is generally substantially different from the population entropic OT solution (the optimal solution of the Schrödinger Bridge).
In step 2, we choose o = O, and define uz zo. x1) to be:
1- 2/
ut (3+ (ro mi) = 2f(1 t) (It (L201 + (1 - -t)xro) + 51 - 320
(58)
which can be shown to make the ODE dx.- = (e (ro 21 )dt admit Pe (2o, v1) as its marginal distributions (Tong et al. 2023).
The above choices precisely yield the (mini-batch) Schrödinger Bridge CFM (SB-CFM) introduced in (Tong al.. 2023). Although the marginals of SB-CFM match those of the Schrödinger bridge solution, the entire path measure not since it only finds an ODE bridge solution.
A.3.3 A SPECIAL CASE: DEEP SCHRÖDINGER BRIDGE MATCHING (DSBM) (SHI ET AL.. 2023
As discussed in A.2.4 the IMF/DSBM algorithm is based on the IMF principle where starting from P( hx((o.)) (20)) (a1)
Pref ((ra (oe( ro, 21), repeatedly and alternatively applying the projections P < (P) and P + (P) leads to convergence to the Schrödinger bridge solution. How does this algorithm
in the framework of our Unified Bridge Algorithm (UBA) in Alg. 2 We will see that specific choice of (. (ro x1). ). (zz. 21) (in step 1) and un (x .11 (in step 2)
precisely leads to the IMF algorithm. We describe the algorithm in Alg. 3
In step I, the pinned path marginals P. (2 are set to be equal to pr (x//20 x1) which can be written analytically
Gaussian (59). The coupling ( x1) is defined to be the coupling distribution
pes (ro, 21 ) that is induced from the SDE in the previous iteration (step 3). p : dx, = vo vr, )dl adW, z. In the first iteration where no is available yet, we set (Zzo, x1) == (z())) (z1). ). We need to check if the boundary condition for (49) is satisfied. This will be done shortly in the following paragraph. In step 2, we fix a == Fref= and set u4 (x (ro, .11) == o² 2, log pref (an (2+)) = 1=z
In step we update 0 by solving the optimization (62). the same as 51. with the chosen Pr i.). Q(1o, 211) and ue (a #1).
Now we see how this choice leads to the IMF algorithm precisely. First, due to Doob's h- transform (Rogers & Williams 2000), the SDE dz, = u, (xx )dt adW, with the choice (61) admits ( (xr(aro, 31 )r as its marginals for any (Jn,, ..). Next, the step 3, if optimized perfectly and ideally with zero neural net approximation error, is equivalent to II (II (P9%)s)) )) where Potd is the optimized / in the previous iteration This can be easily understood by looking at the Markovian projection (P) written in the optimization form (88): The expectation is taken with respect to P(x,, 21) that matches Q( 21 P. (ara [20, 21) in 62). and which is exactly the reciprocal projec- tion of pia since Q _1) = pr (æo, x1) by construction. Lastly, we can verify that the choice in step 1 ensures the boundary conditions Po(-)
=o(-), P-) (-). This is because Q(ra, .1) always satisfies (()))) = To (2)) and Q(x1) = To (=(=)): initially Q #1) = To ( m (41) obvi- ously, and later as Q(20, x1) pue (zo. 21) results from the Markovian projection (from step
in the previous iteration). We recall from Sec. A.2.4 that the Markovian projection preserves the boundary conditions.
IMF/DSBM algorithm as UBA in minimal kinetic energy forms. We can reformulate the IMF algorithm within our UBA framework using the minimal kinetic form as described in Alg. fact it can be shown that Alg. 3 and Alg. 4 Fare indeed equivalent, as stated in Theorem A.5 in Sec. A.4.2
'Initially when there is no previous Oats available, the step 3 is equivalent to (Pint)) where pinit no (x)) (1x))) pref (-I2ro, 21) which is already in the reciprocal class R.
19

---

Published as a conference paper at ICLR 2026
Algorithm 3 IMF/DSBM algorithm (Shi et al., 2023) special instance of our UBA.
Input:
The end- noi-point distributions To and (i.e., samples from them).
Repeat until convergence or a sufficient number times:
Choose pinned marginal path { Pt (ar2zo, x1)): and coupling distribution (71) as follows:
(a 21) == p. (2a 01) = (ax (1 t) + tæı, tee - t)1)
(59)
@=o, = (T (21)
initially (if 0 is not available)
(60)
pro (2o. (I) mi)
otherwise
Choose = Oref, and set uz zu, 21) as:
me (=s o Va, log pref (2)11))
31 2t
(61)
1-6
3. Solve the following optimization problem with respect to the neural network es (l,
min E, (..0f (ul (Zx x1) (t(t xx)]
(62)
Return:
The learned SDE dz; = (t, xe + odW, as the bridge problem solution.
Algorithm 4 IMF/DSBM algorithm (Shi et al. 2023) minimal kinetic energy form in our UBA.
Input:
The end- point distributions To and m (i.e., samples from them).
Repeat until convergence or a sufficient number times:
1. Choose Pz 20, z1) = (2; pe, 7FT) where (He, (M) are solutions to the following optimization:
(msrr)) arg min f'' Ep(celso.a))
[ (1)o0 dt
where
(63)
(1) = dpc dt + ac(x pe). at 7 (
1 dy
otet 2%
(64)
dt
Choose a coupling distribution Q( 21) as follows:
To( (21) (ro, 11)
initially (if 0 is not available)
Q(To, =
(65)
otherwise
Choose a == Oref, set u. (2lko. .i) = as (alwo, .1).
3. Solve the following optimization problem with respect to the neural network vy x)::
min ( blue (are [w0, x1) us(t,
(66)
Return: The learned SDE dz, vo(t, dt) + odW, the bridge problem solution.
Our proof in Sec. A.4.2 relies on some results from the stochastic optimal control theory (Tzen & Raginsky 2019). Then what is the benefit of having this stochastic optimal control formulation for the IMF algorithm? Compared to Alg it has more flexibility allowing us to extend or re-purpose the bridge matching algorithm for different goals. For instance, the Generalized Schrödinger Bridge Matching (GSBM) (Liu et al.) 2024) adopted a f'ormulation similar to Alg. 4 in which they introduced the stage cost function that is minimized together with the control norm term. The final solution SDE would not be the Schrödinger bridge solution, but can be seen as generalized solution that takes into account mobblem-specife stage costs. Hence the algorithmic framework Alg. 4is especially beneficial for developing new problem setups and novel bridge algorithms.
A.4 THEOREMS AND PROOFS FOR UNIFIED BRIDGE ALGORITHM (UBA)
A.4.1
PROOF OF THFOREM A.1
To prove Theorem A.T
show the following three lemmas in turn:
20

---

Published as a conference paper at ICLR 2026
1. First, in Lemma
show that after step of Alg. 2 is done, the SDE dr, uy(2)) + odW, Zo ~ mo(-) admits {P) as its marginal distributions where ux(2) P(iz B (2) zo, 21 (x)
A.2
2. Under the assumptions made in the theorem, that is, i) the neural network (a(t, x)'s func- tional space is rich enough approximate any function arbitrarily closely; ii) the step
3 of Alg. 2 is solved perfectly, we will show that the solution to step 3 is (o(t,z)
E/u, (2: To Fi ||,,= This straightforwardly comes from the regression theorem, but we will elaborate it in greater detail in Lemma A.3 below.. We then show that this conditional expectation equals (x(r) defined in 668) i.e., no(1, (1,, = ux(x). This will complete the proof, and assert that dx, ug- (t, dt + TreydWi, 2o 3o its marginals.
(- admits {P. (2;)) } as
3. Practically, the training dynamics of the gradient descent for step 3 of Alg. can be shown to be identical to that of minimizing E (t. (t,, although the proof is very similar to the result (Tong et al. 2023). Hence, in practice, even without the assumptions of the ideal rich neural network functional capacity and perfect optimization, we can continue to reduce the error between uo (t, x) and (z(x) in the course of gradient descent for step 3.
_x(r)]1 This is done in Lemma A.4
Lemma A.2. Suppose (Pr(zIro, æ1)ht be the marginal distributions of the SDE dx. (x a1)dt adW, for given xy and xı. In other words, step of Alg. is done. For
Pa(a)) = / Pi =1 x1) 2n )d( 21)
(67)
1
uc(x) = (x(x) E ( (ul Pt æ1)],
(68)
the SDE dx, dzs))dt + odW, Io To(-) has marginal distributions {P(()))
Proof. For the given To and Iı, we apply the Fokker-Planck equation to the SDE dx, (2 d1) odW, with its marginals (P((I)ro,
=
0 (( v1) - div (P.((o(, x1 (e (a) 21))
a² AP((c)
(69)
at
2
where div is the divergence operator and A is the Laplace operator. Now derive the Pokker-Planck equation for the target SDE as follows:
① P((x))
a 0t
/
(70)
at
f =
(71)
0t
= / (- (iiv(PP((c/20
+ d²
A 21) Q(z0 xn (d)(
(72)
2
a²
A /
= div EQ [r (2/ 21 P.(x ]]
2
(73)
o² AP(())
div (((x)) Px(x)
Eq
Ue (o (1)
(74)
o²
(x)
2 (P)())
(75)
This establishes the Fokker-Planne equation for the SDE da, = u,(xx) dt + adW,, xo ~ no(-), to
which {P) is the solution. This completes the proof of Lemma A.2
Lemma A.3. Under the assumptions of ideal rich neural network capacity and perfect optimization made in the theorem, the solution step of Alg. 2 that is.
0r arg min E 1x, v1 ) Pt (o) Il (1.. (ol/wo x1) uo vo(l, x+)
(76)
satisfies vo- (1, EE (1/4
ux( where ux(x) is defined in 68.
21

---

Published as a conference paper at ICLR 2026
Proof. We prove the second equality first. The expectation (2u (2y [2y] is taken with respect to the distribution R(xo, 21 (z)) defined to be proportional to (z ( I1 (1Q
/ R(zu, 41 ( (2) To. 21 2 21)
E uu (z, (o. 21 [=1]
(77)
Po Po 21 (1) 21)
P x1 x1) (1)(z, -
(78)
Pa (r no, an ( 01)
ue (2, )d)o 21)
(79)
P((2))
1
ʃ
Pe( T1
2n ( (x) 21)(zo,
(80)
P((,,)
1
BQ( (1sla 21 (2 bo 21]
(81)
Pz(2))
(xx)
(82)
We now prove the first equality. Although this straightforwardyy comes from the regression theorem, but here we will elaborate it in greater detail. Due to the assumptions, the optimization (76) can be written in functional form as:
v" = arg min E, v(.--)
(1 x1) vft, x1)1
(83)
where its optimizer v² (t) equals the optimizer vg- (t, x) of (6 In the functional optimiza- tion (83, the objective completely decomposed over i, and we can equivalently minimize
E PP ( [ x1) vtt, 22)1 for each t. We take the functional gradient with
respect to v(t. -). For ease of exposition, we will use simpler notation where we minimize Ep(v.))lffyy, g 12 with respect to the function g(-). Hence there direct correspondence:
Ies 2 to (J0, .1). fly 2) (1(21/20, x1), and g to v. We see that at optimum.
Qg(y) 2p(y, <)((v) fly, 2))dz
(84)
leading to gy) = [y 2) lv]. This regression theorem implies v (t, = Elue(Iz/ x1 [=a] This completes the proof Lemma
Lemma A.4. VoEe,te)lloo(t, for each I, where Pt (x( and ux(x)) are defined as 66 and (88) respectively.
wr(x) 1a |/ u.(x
Proof.
Pa) (1))) P (t, 20o(1, x) me(xr)]
(85)
V (l 2V.E ((dxx)Tu(()))
(86)
The first term can be written as:
(oBp,c)lloo(tt=))1 Vo / [too
(87)
= Vo / /leo(t, P/ (0 2dd(,, x1,
(88)
Vol (
(89)
The second term can be derived as:
Vo vo (1, Tu(()) Vo (n((,))))nnx)))(((())()))
(90)
= Vo v (u
(91)
Ve / [o(t,=)Tu((z200.
uo(2 21) 21)0011100000)
(92)
V VolEp xt) (x) (2| m1)
(93)
22

---

Published as a conference paper at ICLR 2026
Combining (9) and 13). and noting that || æ1)ll is independent of 8, complete the proof Lemma A.4
A.4.2
EQUIVALENCE BETWEEN ALG. AND ALG.
Theorem A.5. Under the same assumptions as those Theorem A. Alg. 3 and Alg. lead to the same SDE solution, which is the Schrödinger bridge matching.
Proof. The proof goes as follows: i) We first show that the optimization problem (63 64 in step 1 in Alg. can be constraint of Gaussian marginals { (1 [o r1)); 5 ii) We then relax it to an unconstrained version, instance of stochastic optimal control problem with fixed initial; iii) The latter is then shown to admit Gaussian pinned marginal solutions following the theory developed in (Tzen & Raginsky 2019). thus proving that Gaussian constraining does not essentially restrict the problem. It turns out that the optimal pinned marginals and the optimal control have exactly the linear interpolation forms in (59) and (61). respectively, which completes the proof.
constrained minimal kinetic energy optimal control problem with the
and view the unconstrained problem
In step 2, since we set up (I (r] . x1) ar (x(2)/0 21) where is the solution to 663 64). we will use the notation u place of a throughout the proof. First, we show that the SDE dr. Ue ue (xx [ro, x1 dd
odW, with initial state Ip at t=0 and u satisfying (64). admits (Fr x1 - (xx;; 7/ T) as marginal distributions. Note that we must have (Ag = Io. μ1 = 21, To - 0, and Y1 - due to the conditioning (pinned process). This fact is in fact an extension of the similar one for ODE cases TTngg et al. 2023). We will do the proof here for SDE cases. We will establish the Fokker-Plame Planck equation for the SDE, and derive:
aP, (x 21)
PlogN(az Me: I!
(94)
Ōt
2é d (x Mc)
Mc)TTL | Mel 13
= =
(95)
7
where 14 and % are the time derivatives. We also derive the divergence and Laplacian as follows:
div ( (ar x1 (x)ro,)
||x - s 12((1
(3( (-
Lé dl +
(x (x-) Th)T(, p4.
2%
d
(96)
Ftt
77
98
2
|| - 10 d
AP, (x = Tr (V
(x( x1)
(97)
Vt.
From 95
96
and 997. we can establish the following equality, and it proves the fact.
aP: (z ro, 21)
div ( mm) + o²
(2( 21)
(98)
Ot
From the above fact,
re-state the step of Alg. follows:
(Step re-stated) Choose Pt (t (xl 21) as the marginals of the SDE, dx, = ( x1)dt adW, with initial state To 1=0 where u. the solution the constrained optimization:
min (.(( Ju
(1
} (1/ m1)1 dt s.t. ( (P((z)aa, are Gaussians
(99)
Instead of solving (99) directly, we try to deal with its unconstrained version, i.e., without the Gaussian marginal constraint. To this end we utilize the theory of stochastic optimal control with the fixed initial state (Tzen & Raginsky 2019), which we adapted for our purpose below in Lemma A.6
In Lemma A.6
adopt the Dirac's s delta function g(-) ō., for the terminal cost to ensure that the SDE dr, ( jro, 41 )dl + odW, with initial state do lands at æ1 as the final state. Then the unconstrained version of (99). which is perfectly framed as optimal control problem in Lemma A.6
has the optimal solution written as:
u (2/ 21) V²V, log E pref (8z, |xa 2=| o² o² log pref (ril(1 x)
(100)
23

---

Published as a conference paper at ICLR 2026
which coincides with (6 in Alg. B Now, due to Doob's h-transform (Rogers & Williams 2000).
the SDE dx. We ( 20 21 )dt + adW, with the choice 61 (100) admits (P. (x) x0, (1)
) as
its marginals. In other words, Pt (r x1) = P (ar (ro, x ), which Gaussian, meaning that the constrained optimization (99) and its unconstrained version essentially solve the same problem.
Noting that (59 and 61 Alg. are equivalent to (63) and in Alg. we conclude that Alg.
and Alg. 4 are equivalent.
Now we describe Lemma A.6
Lemma A.6 (Stochastic optimal control with fixed initial state; Adapted from Theorem 2.1 in Tzen & Raginsky 2019)). Let pd be the path measure of the SDE: dx; (r) + odW, starting from the fixed initial state rv. For the stochastic optimal control problem with the immediate cost 202 [,,((zz) (2+) at time and the terminal cost log 1/g(zz) ) at final time t = 1 for any function g. the
cost-to-go function defined as:
f1(x) = Eps J.
202 [bb(( -log))))llloaaaaa
(x1) 2t =
(101)
has the optimal control (i.e., the optimal drift bx(x))
by (x) arg min (z) o² log (g(x))) x]]
log
(102)
where pref is the Brownian path measure with diffusion coefficient a.
Proof. We utilize the (simplified) Feynman-Kaa Kac formula, saying that the PDE,
(he(z) me (x( V (x + ( 0, (-( x-)
(103)
Öt
has solution hạ(2) = E[g [x) = æ] where the expectation is taken with respect to the SDE,
dæ, = (x,)) a.W.
Now we plug in 14 = 0, q g. and let (x) == log (xr). Note that hx(r) is always positive since is positive, and henee u2 (x) well defined. Then by some algebra, we see the following PDE:
@v((x) 1 Tr (
o²
2 [[))(()))
vI(+) logg(.)
(104)
Ot 2
has solution By (z) = log E p [g ((1) = æt. Note that we can write " ||Vvv(<>)]' as the following variational form,
| |V -min 2
min (xv( +
(105)
20
where the minimum is attained at 6* = - (z))). So (x(x) = - -log pres ( (x11))= solution to:
= x] is the
Ovz(x) 1 Øt TT (T(("))01(1))
|B|12
min (z((()+
202 e1(.) logg(.)
(106)
Note that (T6) is the Hamilton-Jacob-Belllne lem with the immediate cost 2/ (xx and the terminal cost log 1/g(21). In fact, vx( log E Eprer [g(x1)) 222, zz| is the (optimal) value function, and the optimal control, i.e., the solution to 105) function form, is: by (x) = -
equation for the stochastic optimal control prob-
V log Eprer (g(xn) Jar, = x:.
B TECHNICAL DETAILS OF GAUSSIAN PROCESS GSBM
B. ELBO DERIVATIONS
First, the posterior path measure can be written as:
P Ppost (P.) P Pprion (P.) exp(-J((P.. V.)/T)
(107)
E(n,T)
24

---

Published as a conference paper at ICLR 2026
where E1), r) is the evidence that depends only on the model parameters 7 (prior kernel hyperpa- rameters) and T (the likelihood hyperparameter), and defined as a normalizer:
En,, r) Ep prrtor (Pa) ( V.)/T)]
(108)
We start with the KL divergence between the variational posterior and the true posterior measures.
dQ(P.)
(L(O((.)) post (P.) = (E((P))
log
(109)
posttt (P.)
log E(1, r) + Re()))) [PP: Vo)/T] + (L( priir (Po)
(110)
Now let Y and Y' be arbitrary time index sets in [0, 1]. The last term in 110 can be written as:
KL KLpYy Salyy,He .Z yy, (2.))1111 P (NY, (UY,,z, 1z, ;))))))))
(111)
Q(zz) Perior ( Q(fzz) Pprror
Be
log P Pprior (Hz) Ppstorttyy 12) Perior ) Pprtor(fsT
(112)
log (uz2) Q(īz.)
Be
(113)
Pprior (uzz) Pprior (7z.)
KL((QUz)) prior (pZ)) + (L(((((( (Prrioo((
(114)
(m)) (HN (2 KL 2.)) (115)
where 1122) comes from our construction of the variational density in [4] with ((aluz)) set
equal to prior (p. (uz) (and similarly for 7). Eq. (155 is due our Gaussian modeling of the inducing variables in (15). Note that this derivation starting from index sets Y and Y' results in the quantity (155) that does not depend on them. That is, (15) can be safely plugged in for the
last term of 11 which we denote by (L((((22.)) (Pz.z)). Due to the non-negativity of
(L((((.)) post (P.)), we have:
E(n 2 - -EEurr) U(P.:V))/+ + KL
prior (P2,2)))
(116)
where the right hand side is the evidence lower bound objective (ELBO). Lastly, the variational-prior KL divergence Gaussians, and can be written in closed forms follows:
Pprroo (Pz.z:) )), is confined to n-dim
(L(2(PZ,z))1PP))) prior (Pz,z:))
(117)
Tr
+ +R²))
log [T )-)
(118)
2
- + (C
- -m -n' + log 162, 2.1 log 1
(119)
where we assumed the same kernel functions over dimensions j of pc. denoted by L == 1/2,1
Vj. n'n = 12' and defined as:
d
S"
Su. R'
CO.. - , T" =
log S (120)
d
1=1
j=1
B.2 COMPLEXITY ANALYSIS
In this section analyze the computational complexity of our GP-GSBM and GSBM. Since the SDE training part in both methods are identical, we mainly compare the ELBO learning of our GP-GSBM with the corresponding CondSOC optimization in GSBM. As described in Sec. 3 the key steps in ELBO learning consist of: i) sampling (He, () from the variational posterior Q. ii) computing Ore (xte (2n which involves time derivatives (e.g., j4), iii) evaluate the likelihood. ((..; V.), and
25

---

Published as a conference paper at ICLR 2026
Table 5: Complexity analysis for ELBO learning in our GP-GSBM and CondSOC optimizaiinn in GSBM (Liu et al. (2024).
GSBM (LIU ET AL.
2024
GP-GSBM (OURS)
SAMPLE/ESTIMATE (Mt, Ye) COMPUTE a (21//zo.. z1)
O(NTnd)
O(NTnd + NNT")) O(Td)
((Td)
EVALUATE L.KEL.HOOD/L.OOS EVALUATE KL((((2))(Prior)
O(T (d + Tv))
- (v)) O(Ñn² + +Nn'd)
iv) evaluate the variational-pria -rior KL term. The former three steps decide the first term of (1 while the latter determines the second term. The first three steps have corresponding counterparts in GSBM: i) estimating (/4,7) 7/)) from the spline model, ii) computing or(3); (1)) and iii) evaluating the loss in (3. It does not have the KL. computation.
First, the (Me, "i)) sampling step in our GP-GSBM. Following 17 20 and recalling that used diagonal 5 and the same kernel function across all dimensions j. sampling (/4, 7r) takes O(NTnd NTn") time where T' is the batch size, N is the number of ELBO learning iterations, and Ñ = 1 if no model selection for the kernel hyperparameters 17 is required, and N = N otherwise. This is because if the kernel function fixed, the kernel matrix inversion can be amortized over iterations by caching inverse kernel matrices. Using the same notation in GSBM, where n denotes the number of spline knot points, the spline estimation of (1e, Tt.,, corresponding to the sampling our case, takes O ONNTnd).
Computing ar (x |z0, x1) in 5 takes (TT)) for both our approximate time derivatives and GSBM's auto-differentiation. Evaluating (P':; V.) in 7 takes O (d + Tv)) where Ty is the cost of computing (z GSBM incurs the same cost for the loss evaluation. The variational-prio KL term, only appearing in our GP-GSBM, takes O(N Nndd) complexity following 117 120).
We summarize these results in Table 5 Comparing with GSBM, have additional cost, mainly (n²)) kernel inversion cost when model selection is performed for the kernel hyperparameters and O(n2d) in the KL divergence computation. This may be costly for large d even though we have constant n (usually no greater than 30). One work around is to use diagonalized kernels in (188 119). which reduces the complexity to linear in n. In practice, we observed no significant difference from using full kernel matrices. In this case, GP-GSBM incurs the same asymptotic complexity as GSBM when we do not use model selection during ELBO learning.
C EXPERIMENTAL DETAILS
We mostly follow the experimental settings, default hyperparameters, and network architectures from GSBM (Liu et al. 2024) and their official code base") DSBM (Shi et al. (2023) is performed by turning off the CondSOC optimization GSBM, which initializes the marginal path as linear interpolation of DSBM. Our GP-GSBM is also implemented based on their code base where we similarly take alternating forward and backward SDE training and integration for coupling (1) samples. For the kernel implementation, we adopted the GPyTorct Iibrary (Gardner et al. 2018).
Problem setups. For the stage cost for crowd navigation in 24). we use Xohs 1500, Aent 0, Auget 50 for Stunnel, and Aohe = 1500, Aent = 5, Angst = 0 for GMM. For the noisy scenario in LiDAR, the manifold projection altered by (I + 0.5c) in multiplicative way where e ~ (0, I).
Hyperparameters. The numbers of inducing points in our GP-GSBM, n for μ4, and n' for Tt correspond to the numbers of spline knot points in GSBM. And we match them for both models for fair comparison. Specifically we chose: Stunnel and GMM (n 15. n' 30), LiDAR (n 30, n' = 30), and CIFAR- 10 and AFHQ (n = 8, n' = 8). For all experiments, our variational inference ELBO learning, corresponding to the CondSOC optimization in GSBM, takes the number of iterations: 1000 (Stunnel), 2000 (GMM), 200 (LiDAR), and 100 (CIFAR-10 10 and AFHQ). Following (Liu et al.) 2024). we use the Adam optimizer (Kingma & Ba 2015) for AFHQ (and CIFAR-10), SGD with momentum for LiDAR, and SGD for the rest. For the SDE noise level, o = 1.0 is used for crowd
zed sccr-schrodi nger b bridge--hhiiddmmtcciiiii
26

---

Published as a conference paper at ICLR 2026
(Deterministic obstacle scenario)
(Uncertain obstacle scenario)
M₁
H
Q(44
H4
@(1+)
7
Y
(Gr,)
12-
X
X
DSBM
GSBM
GP-GSBM (Ours)
GSBM
GP-GSBM (Ours)
Figure 4: Stunnel results. The two ellipsoidal obstacles are shown as gray regions. For each panel, show Pt. Tt (top and middle rows) of the learned pinned marginals (r x1) = NW4,, 7f I), and samples from the pinned marginal at 5 uniform time points (bottom row). The two boundary distributions are isotropic Gaussians where To is shown blue points, and æi as purple. For DSBM (Shi et al. 2023) and GSBM (Liu et al. (2024), the learned μl,, "t are shown. For our GP- GSBM the learned variational posteriors ((1+)) and Q(Y)) are depicted as means and 2x pointwise stdevs (95% confidence intervals), the latter as red-shaded regions.
navigation tasks, and = 0.5 for AFHQ. For CIFAR-10 (a non-GSB task), small known to be preferred as alluded in (Shi et al. 2023), and we use 0.001. For SDE integration, we take 1000 time steps. For any other hyperparameters that are not discussed here, we take the default values in the official code base of (Liu al. 2024).
Model selection. Due to the extra cost of the gradient-based kernel hyperparameter (7) search, we use grid-based empirical Bayes. More specifically, for the squared exponential kernel that used in our experiments, the width parameter is chosen from the candidate set {0.05, 0.1, 0.5, 1.0), where select the that yields the highest ELBO score. This allows for expedited GP inference since the kernel inversion can be done only once at the beginning and amortized over variational inference iterations. For the likelihood hyperparameter
gradient-based empirical Bayes.
Network architectures. We describe the architectures for the SDE network Up (t. x).. For the crowd navigation problems, adopt the fully connected network, the same as (Liu al. [022], that has 4 to 5 residual blocks with sinusoidal time embeddings. For the AFHQ task, we use the U- Net (Ronneberger et al. 2015) architecture with implementation from (Dhariwal & Nichol 2021)_ For the CIFAR-IO experiment, we follow the same U-Net architecture that implemented in CFM (Tong
ct al. 2023), which uses 4 heads and 16 attention resolutions with two residual blocks. All networks except for the CIFAR-10 case are trained from scratch, and optimized with Adam (Loshchilov &
Hutter 2019). For CIFAR-10, we warm-started training with the trained CFM checkpoint similarly
as (Shi et al. 2023 For inference, we use the exponential moving average with rate 0.9999.
Computing platforms. All experiments for Stunnel and GMM were conducted on a single RTX-2080- Ti GPU while two RTX-2080-Ti GPUs were used for CIFAR- 10. For LiDAR, single RTX-4090 GPU was used. For AFHQ, we train our GP-GSBM on 32 V100 GPUs (4 nodes and 8 GPUs per node) where a similar setup used in GSBM training in (Liu et al. 2024).
D ADDITIONAL EXPERIMENTAL RESULTS
D.1 MORE FIGURES FOR AFHQ IMAGE-TO-IMAGE TRANSLATION (DOG CAT)
After training the model, we pair the coupling of (æ1, 2D) ) by rumning the trained (reverse-time SDE model starting from #1 mi, which is shown in Fig. 6 where #1 (data) samples the rightmost.
27

---

Published as a conference paper at ICLR 2026
(Deterministic obstacle scenario)
(Uncertain obstacle scenario)
M4
Q(4)
1
2(4))
Y
7
((%)
DSBM
GSBM
GP-GSBM (Ours)
GSBM
GP-GSBM (Ours)
Figure 5: GMM results. The three blob obstacles are shown as gray regions. For each panel, we show 14,, "t (top and middle rows) of the learned pinned marginals Pt (2t ( N NN4,,, 7 T), and samples from the pinned marginal at 5 uniform time points (bottom row). The two boundary distributions are mixtures of isotropic Gaussians where Tg shown blue points (4 Gaussian blobs in the center), and и as purple (8 outer Gaussian blobs). For DSBM (Shi et al. 2023) and GSBM (Liu al. 2024). the learned (4:, Yt are shown. For our GP-GSBM, the learned variational posteriors Q and ((r)) are depicted as means and 2x pointwise stdevs (95% confidence intervals), the latter as red-shaded regions.
Table 6: (CIFAR-10) FID scores on this non-GSB image generation problem. The results from CFM and Stream-level GP are excerpted from (Wei & Ma 2025).
CFM
STREAM-LEVEL GP
GP-GSBM (OuRS)
I-CFM
OT-CFM
I-CFM
OT-CFM
3.75 1 0.006
3.74 1 0.009
3.62 0.008
3.75 0.009
3.65 0.008
columns, and zo (generated) samples on the leftmost. From this coupling. we show the mean images from Pt (x: x1) for the prior and posterior (before and after the variational GP inference) in the middle. Visually clearly, the posterior exhibits semantically more meaningful path samples than the prior with less artifacts since the posterior takes into account the latent- preettpreserving stage cost, forcing the samples to stay closer to the image manifold. In Fig. 7 the SDE generation progress is shown with intermediate time steps. We also visualize the generated cat images in Fig.
D.2 10 RESULTS (NoN-GSB COMPARISON STREAM-LEVEL GP)
We follow the typical CIFAR-10 experimental setup from Tong et al. 2023 Wei & Ma 2025) where To A00 I) and कi is the CIFAR-I0 data. First we that this experimental setting may not be an ideal example for SB matching problem as indicated by the fact that CFM with independent coupling (I-CFM) is comparable or sometimes better than mini-batch OT ba-ed coupling (OT-CFM).
For the FID computation, we use 50K samples generated and averaged over 20 random runs. The results are summarized in Table ( where means and standard errors are reported. Comparing our GP-GSBM, stream-level GP, and CFM, most methods perform equally well on this (non-GSB) noise- to-image generation problem. Some images generated by our model are shown in Fig. (diffusion from To to TT) and Fig. (final generated samples).
28

---

Published as a conference paper at ICLR 2026
t=0
t=1
Prior
Post
Prior
Post
Prior
Post
Prior
Post
Prior
Post
Prior
Post
Prior
Post
Figure 6: GP-GSBM on AFHQ dog (t=11 1, rightmost columns) - cat (t= 0, leftmost columns). Mean images from (a 21) for the prior and posterior (before and after the variational GP inference).
29

---

Published as a conference paper at ICLR 2026
t=0
t=1
Figure 7: GP-GSBM on AFHQ dog (t=1, 1, rightmost columns) - cat (t=0, leftmost columns). The (reverse-time) SDE generation progress.
30

---

Published as a conference paper at ICLR 2026
Figure 8: GP-GSBM on AFHQ dog + cat. Generated cat images starting from dog data samples.
31

---

Published as a conference paper at ICLR 2026
Figure 9: (CIFAR-10) SDE generation snapshots by our GP-GSBM. From the leftmost columns
(samples from To (O, I)) the generated images Ti on the rightmost columns.
Figure 10: (CIFAR-10) Generated images from (0,
from./((0, I) by our GP-GSBM.
32

---

Published as a conference paper at ICLR 2026
E REVISION
This section contains added materials and new results from the reviews/rebuttal. We thank all the anonymous reviewers for their insightful and constructive comments and suggestions.
E.1 ToY (THOUGHT) EXPERIMENT: WHY GP-GSBM Is ROBUST COMPARED TO GSBM
In this section we provide the key intuition on our (Bayesian) GP path modeling for the generalized Schrödinger bridge problem, illustrating why our GP-GSBM more robust to noise than the (deterministic) GSBM (Liu et al. 2024) algorithm.
To focus on the core insight, we consider highly simplified scenario, in particular, we make the following simplifications: i) the CondSOC cost (or the negative log-likelihood score in our Bayesian model) in ( queried from a black box oracle, meaning that it returns (stochastic) cost for a path query, ii) the posterior inference in our GP-GSBM is done exactly without resorting to approximation strategies like the sparse variational GP. For the former, although the CondSOC cost is composed of the kinetic cost and the problem-specific cost (z), we treat it single number retrieved from the oracle for simplicity.
In the 2D state space, we have two end points xy = (-1,0) and x1 (+1, 0) sampled from the
coupling distribution Q( Let's consider N = 1001 paths connecting the two end points whose trajectories are given by parabolic functions Y a(X² - 1), that is, [z,]2- a(æz,]]
for
a i/500, i E +1.... +-500} where I2 = ([ra]1, [x+]2). They are visualized in Fig. 11 We assume that only these N paths incur finite (< 00) cost (to be defined shortly), whereas any other path connecting To and 1 incurs +oo CondSOC cost.
Among these N 1001 paths, let's say all paths incur CondSOC cost +1 except for the horizontal line path Y = 0 (i.e., a = O, shown as red in Fig. [1] that incurs - e where = 0.001 most of the time. Specifically, the CondSOC cost for the horizontal line path stochastic/noisy defined as:
Jy= f1-6
with prob. 0.999 with prob. 0.001
(121)
100
Whereas the other 1000 paths incur cost +1 with probability 1, the path Y = 0 gives us slightly lower cost 1 - most of the time, but can have the rare event of incurring extremely high cost 10
Now we are going to train the GSBM model and our GP-GSBM model, but let's assume that this rare event of the horizontal line path Y = O was not observed in the training data due to its rare event probability (0.001). We also note that due to the finite number of distinct paths with finite (< oo) cost, the optimal variance of the pinned marginal paths p( [2 21) for any sensible model must be 0 (i.e., " 0), and hence we only focus on the mean μ4c
First, GSBM is deterministic point estimator, and thus it will learn the horizontal line Y 0 path as its optimal path since incurs slightly lower cost than any other 1000 paths, and we didn't observe the rare event in the training data. We can compute GSBM's (population) expected cost test time (denoted by risk), which is
RiskosBM S (1 - c) 0.999 10 0.001 s 1001
(122)
Apart from GSBM's deterministic solution that only finds point estimate (a single path) and exclude all the other possibility, the posterior in our Bayesian model offers us distribution over set of plausible paths under the given observations (observed costs). When we consider very flat GP prior (i.e., nearly uniform prior over those N paths), and we use the temperature T = 1 for the likelihood model, the posterior will be proportional to ed, more specifically.
1+
if path =Y
NN
P Ppoet (path) =
path = "Y =Y=a for a E { 500....; + 500
(N
b
otherwise
(123)
33

---

Published as a conference paper at ICLR 2026
Xg
X1 (+1,0)
(-1,0)
Figure 11: Toy (thought) experiment setup. N 1001 paths connecting zo = (-1,0) and x1 (+1,0) are shown, which incur finite (< oo) CondSOC costs. The horizontal line path Y 0 (shown red) incurs stochastic noisy cost Jy defined in (T2I). The rest 1000 paths (blue) incur fixed cost J = 1. The GSBM (deterministic) algorithm finds the red Y = 0 path as a point estimate, which yields the risk s 1001. Our (Bayesian) GP-GSBM learns the near-uniform posterior distribution over those 1001 paths as in 123 achieving the risk R2 2.
The risk of GP-GSBM then computed as: RisKGP-GSBM
Epparr (parchpaahh) [Enoise [Jpari]]]
(124)
(N 1) -- (N - 1)e te 11- ((1 (()-c) 00.999 + 10 0.001) + (N 1)e +e-l+e 1tc
c -+
1
This thought experiment, albeit quite simplified, provides the key intuition on why our Bayesian (GP) formulation for the generalized Schrödinger bridge problem is more robust to noise than (deterministic) GSBM (Liu al. 2024).
E.2 BACKGROUND ON GAUSSIAN PROCESSES (GP)
The GP models are special Bayesian models where we impose prior distribution over function, and the distribution has the property that for any imput index set, the function values on the set follow Gaussian. Specifically, input point & E X, is a latent variable. The Gaussian process is characterized by the mean function m X IR and the positive definite covariance (or kernel) function k X x X R. defined as:
consider the latent function
X - IR where fæ). the function value at
E[f(z) = m(x), Cov(f(z), = k(x,x')
(125)
for any x, E X. For any input set X = (-1,,-.,, the latent variables ( jointly follow the multivariate Gaussian with mean (m(2)),.... (=()))).....,(;,))) )) and the covariance ma- trix K(X) whose (i,j) entry is k(z,, xz). The GP prior distribution P(f(-)) is often denoted by GP(F():()) -),, and used this notation in our main paper.
(*x,)
Once we impose a GP prior on latent function f()), and once we have likelihood model P(DIS(-)) that defines the likelihood of the observation (or observed data) D with respect the function f(-). we can follow the Bayes rule to infer the posterior distribution over the functions. That is,
P(f(-)D) & GP(F()) P(D)f())
(126)
Note that there is no restriction on how the likelihood model P(DIf(-) is defined. Each observed instance in DD can be dependent on some or all of the function values [f(x))=EA JaEA for some A C X. For instance, in the traditional GP regression (Rasmussen & Williams (006)) where we model each input/output pair y f(x) + for some noise the likelihood of each observed instance (z, y) D depends on f(x) at the point æ alone. However, for the CondSOC based likelihood model (7) and 99) in our GP-GSBM, and the likelihood at each time f depends on not just (j(tt,, (t)), but also those in their infinitesimal neighborhood (p(t+ + At). 7tt At)) with At 0. This is because of the term a, defined in (5).
have two types of latent functions μs (or μa(-)) and % (or 7()),,
Except for a few special cases (c.g., GP regression with Gaussian noise), most likelihood models in practice including the one in our GP-GSBM, lead to non-GP posterior P))())D) For tractable
34

---

Published as a conference paper at ICLR 2026
solutions, one has resort to approximate GP inference algorithms (Rasmussen & Williams, 2006). In this paper, we adopted the sparse variational free-energy approximate inference algorithm which is
known to take various advantages (Titsias 2009 Dezfouli & Bonilla 2015 Matthews al. 2016 Bauer et al. 2016).
E.3
ENLARGED VERSION OF FIG.
In Fig.
12
we show an enlarged version of the images in Fig. |
for better readability.
(Noise-free scenario)
-2
DSBM
GP-GSBM (Ours)
Stream-level GP
(Noisy scenario)
GSBM
GP-GSBM (Ours)
Figure 12: LiDAR path samples for different models 3D (top rows) and 2D (bottom) views.
Starting from the samples of ma (red points the left), generate samples for blue points the right. Path samples are visualized in uniform time points from / Ł=0 0 (red) to f t- 1 (blue).
35


Reviews:
[{"review_id": "DNfiRznZn8", "reviewer": "Reviewer_gihK", "summary": "This manuscript tackles the generalized Schrodinger Bridge (GSB) problem that incorporates an extra stage cost function on top of the original SB problem. Given that the previous efficient GSB solver only offers a point estimate of the path measure, the authors propose to treat the GSB as a posterior inference problem. Specifically, they treat the GSB objective as a stochastic likelihood function and leverage a Gaussian path prior, and aim to estimate the posterior distribution. To address the intractable posterior, the authors propose to leverage fully factorized variational inference, where the variational posteriors are tractable GPs. Empirical evaluation demonstrates that the proposed solution has competitive performance with the GSBM baseline and superior performance on the scenarios where there is uncertainty or under a proxy cost function.", "questions": "See weaknesses.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The writing is clear, with good contextualization of related works.\n* The motivation of this paper is good.\n* The proposed technique is sound, with well-executed experiments.", "weaknesses": "Overall, I have no major concerns about this manuscript. \n\n### Minor\n* The implications or limitations of some design choices should be highlighted. For example, will the independent prior (regularization) be suitable for all applications? The same question applies to the homogeneous kernel function across dims (though this is common) and the form of the variational posterior. However, adding complexity to these designs may further increase computational cost.\n* The figure's resolution can be improved. The current form is somewhat hard to read when the paper is printed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "hZ6QWwnXb1", "reviewer": "Reviewer_Jgjn", "summary": "This study proposes a variational solution that applies a Gaussian Process (GP) to efficiently solve the existing Generalized Schrödinger Bridge (GSB) problem. To define the Gaussian process, Gaussian distributions are first established based on the data from two domains, $x_0$ and $x_1$. A GP prior is then constructed between them, and learning is conducted by optimizing the Evidence Lower Bound (ELBO) over intermediate timesteps. Using this approach, the proposed method achieves superior performance compared to existing approaches.", "questions": "See Weaknesses Section.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- Simplification through a variational approach: The proposed method applies a Gaussian Process (GP) to the GSB problem, allowing convenient formulation of intermediate distributions through a simple prior. This variational perspective simplifies the overall problem structure while maintaining flexibility.\n\n- Visualization of results: The paper thoroughly validates the proposed methodology through diverse visual analyses presented in both the main manuscript and the appendix. These visualizations effectively demonstrate the distinctions and advantages of the proposed approach.", "weaknesses": "### Major Weaknesses\n- Time Complexity: As also mentioned by the authors, the proposed method has significantly higher time complexity compared to the original GSBM. According to Table 5 in the Appendix (Page 26), the method additionally requires $O(n^3)$ time. A strategy to further reduce the sampling time is needed. Considering the current computational cost, the performance improvement over GSBM does not appear dramatic.\n\n- “Sparse” Variational Gaussian Process: The proposed method constructs a flow by interpolating simple single Gaussian kernels defined on two domains, $x_0$ and $x_1$. While this simplifies the problem, it does not fully capture the global distribution, making the selection of data pairs highly influential to overall performance.\n\n- Uncertain Performance Difference: The proposed approach modifies the GSBM structure only by formulating $V_t$ as a GP problem. Based on the reported results, including those in the Appendix, the improvement over GSBM appears only incremental. It would strengthen the paper to include concrete experimental cases where GP-GSBM successfully solves problems that GSBM fails to address.\n\n### Minor Weaknesses\n\n- Choice of Primary Area: The proposed method utilizes a Gaussian Process (GP) as one possible solution to the GSB problem. In my opinion, this work may fit more appropriately under the area of generative modeling rather than GP itself.\n\n- Formatting in the Related Work Section: It might be better to remove quotation marks from the paragraph titles in the Related Work section for cleaner presentation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "pVjKjtRRie", "reviewer": "Reviewer_ofu8", "summary": "This paper extends the Schrödinger Bridge (SB) framework by introducing a robust generalized formulation that unifies diffusion-based generative modeling and SB matching under a common theoretical lens. The authors develop a new objective function incorporating robustness to model mismatch and derive corresponding training dynamics. The work provides both theoretical guarantees and empirical validation on complex datasets, aiming to improve the stability and expressiveness of SB-based generative models. But to be honest, this paper is too mathematically dense and is hard for me to follow.", "questions": "Could the authors provide a more intuitive explanation or diagram illustrating how the generalized robust SB formulation differs from the classical one, and how robustness manifests in practice?", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper tackles an ambitious and mathematically sophisticated problem, presenting a unifying and robust extension to the Schrödinger Bridge framework. The theoretical development appears rigorous, and the results—both in proofs and experiments—suggest meaningful advances in understanding and improving bridge-based generative modeling. The connection drawn between generalized SB matching and diffusion-based models is conceptually strong and of potential significance to the community.", "weaknesses": "The paper is mathematically dense and can be challenging to follow, even for readers with backgrounds about SB (not GSB). Some of the derivations could benefit from more intuition or interpretive commentary. It is also difficult, without deep verification of each step, to fully assess the soundness of the theoretical proofs.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "MUFqbjqiGj", "reviewer": "Reviewer_PMLn", "summary": "This paper proposes a GP prior GSBM algorithm for robust generalization. This allows a sparse variational free energy formulation that ensures effective posterior path measures using GP-based mean and covariance derivation on pinned marginals. The noticeable benefits are robustness in noisy scenarios and generality for various SDE settings. The authors presented various experiments demonstrating its effectiveness.", "questions": ".", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The proposed method and ELBO optimization method can be applied in various circumstance.\n* The overall manuscript is clearly written.\n* The authors put considerable effort into experiments, and they successfully found the method's strength in robustness.", "weaknesses": "* Although it is referred to as generalization, I think the GP prior is another constraint or additional premise on SB that would only favor some portion of SDE problems in theory. I think a GP-based formulation can be overly restrictive for very high-dimensional modalities with additional geometric constraints that have branching paths as solutions.\n* The precise notion of robustness is unclear, and why the GP prior recovers the \"original\" solution amidst noise is black-boxed.\n* There is a room for further theoretical discussion on the GP prior which might benefit the clarity manuscript.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)