PDPO: Parametric Density Path Optimization
Authors: Sebastian Gutierrez Hernandez
Peng Chen
Hao-Min Zhou
Venue: NeurIPS 2025 poster ()
Keywords: Optimal Transport
Generalized Schrodinger Bridges
Lagrantian Costs
Primary area: deep_learning
TL;DR: In this work we focus on action minimizing optimal transport, where a density path minimizes the cost of transportration between two prescribed densities. The cost is determined by a kinetic energy and a potential energy term.

Abstract: We introduce Parametric Density Path Optimization (PDPO), a novel method for computing action-minimizing paths between probability densities. The core idea is to represent the target probability path as the pushforward of a reference density through a parametric map, transforming the original infinite-dimensional optimization over densities to a finite-dimensional one over the parameters of the map. We derive a static formulation of the dynamic problem of action minimization and propose cubic spline interpolation of the path in parameter space to solve the static problem. Theoretically, we establish an error bound of the action under proper assumptions on the regularity of the parameter path. Empirically, we find that using 3–5 control points of the spline interpolation suffices to accurately resolve both multimodal and high-dimensional problems. We demonstrate that
PDPO can flexibly accommodate a wide range of potential terms, including those modeling obstacles, mean-field interactions, stochastic control, and higher-order dynamics. Our method outperforms existing state-of-the-art approaches in benchmark tasks, demonstrating superior computational efficiency and solution quality.

PDPO: Parametric Density Path Optimization
Sebastian Gutierrez Hernandez'
Peng Chen?
Haomin Zhou'
School of Mathematics, Georgia Institute of Technology
School of Computational Science and Engineering, Georgia Institute of Technology
{shern3, pchen402 hwmzoou@@atechh.co
Abstract
We introduce Parametric Density Path Optimization (PDPO), novel method for computing action-minimiza paths between probability densities. The core idea to represent the target probability path as the pushforward of a reference density through parametric map, transforming the original infmitt-dimeensionll optimization over densities to a finite-dimensional one over the parameters of the map. We derive static formulation of the dynamic problem of action minimization and propose cubic spline interpolation of the path in parameter space to solve the static problem. Theoretically, assumptions on the regularity of the parameter path. Empirically, we find that using 3-5 control points of the spline interpolation suffices to accurately resolve both multimodal and high-dimensional problems. We demonstrate that PDPO can flexibly accommodate wide range of potential terms, including those modeling obstacles, mean-field interactions, stochastic control, and higher-order dynamics. Our method outperforms existing state-of-the-ard approaches in benchmark tasks, demonstrating superior computational efficiency and solution quality. Source code htttps://github. http://www....oo/ooooa/a//t/////
establish error bound of the action under proper
Introduction
Optimal transport (OT) theory has become powerful and widely used framework for quantifying discrepancies and constructing interpolations between probability distributions, with rapidly expand- ing applications in machine learning, statistics, and computational science (27]. It offers a principled way to compare distributions by computing the minimal cost required to transform one distribution into another.
While the computation of optimal transport (OT) plans and geodesies between probability distributions is now relatively well understood [31 22 3
34 25 13, many real-world applications require solving more complex, constrained transport problems- specin-iiiiii.llll....... identifying optimal paths between distributions subject to additional requirements. For instance, one may need to guide a swarm of robots from one configuration another while avoiding obstacles, or interpolate between data distributions in way that preserves semantic structure [20] 29 17].
These tasks can be formulated action-minimizinn problems governed by the "principle of least action": Given initial and terminal densities Po, E P(Rd), where P(Rd) denotes the space of probability densities over Rd, the objective is to find time-dependent density path p(t. -) E P(R)) and a continuous velocity field v: [0, 1] Rd Rd that minimize the action functional:
inf '' f K(p,v) dx dt p.v Jgd
(1)
subject to 0,p+ (pv) = 0, p(0,-) Po. p(1,-) =
(2)
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Here, (p, u) denotes the transportation energy, and (()) represents potential term that captures interactions among particles with the environment. When (p. v) 3 olllp p and F(p) this reduces to the classical Wasserstein geodesic problem. Directly solving such metion-mininizing problems poses substantial mathematical and computational challenges. low dimensions (e.g.. 2 or 3), the associated PDE system derived from the first-order optimality conditions- can be tackled using classical numerical methods 5 6 [4]. However, these approaches do not scale well with dimension and become computationally infeasible in high-dimensional settings.
Recent advances in machine learning have greatly expanded the range sigh-dimensional problems that can be tackled effectively. A leading example is the Generalized Schrödinger Bridge Method (GSBM) [20]. originally developed for stochastic optimal control (SOC) problems [24]. GSBM learns forward and backward vector fields by modeling conditional densities and velocities, drawing inspiration from stochastic interpolants [2]. and approximates Gaussian path statisties via spline-based optimization. In parallel, [29) introduced an algorithm for problems with linear energy potentials, leveraging Kantorovich duality and amortized inference to efficiently compute c-transforms and transportation costs. For Mean-Field Games, APAC-Net [17] casts the primal-dual formulation convex-concave saddle point problem, trained using a GAN-style adversarial framework.
Our method builds upon prior work on parametric probability distributions [21] 33 [5], extending these ideas to address boundary-valued action- mimion--iiiizzinggg density problems. Our formulation
an extension of the static OT formulation to the parameter space, whereas in [16], they extend the dynamic formulation. In the remainder of this section, we outline our methodology, with full technical details deferred to Section 3
Figure la offers overview of our approach: curves in parameter space 0(t) @ (left) induce density paths P. E (IRR)) via pushforward maps Pr (Ta())# # À from fixed reference density À (top), while samples of these densities are obtained via direct evaluation (=())=) with À (right).
While spline-based methods have previously been applied particle-level density transport (2 and Gaussian settings [20]. key innovation of our work is the use of cubic Hermite splines in the parameter space of neural network. We fix a set of control points {0,,) 1 with i, = i/(K + and optimize over them to discover eettonnminimizzig density paths, which presents the first spline-based parametric approach in a learned map setting. Figures Tb [b] Tc illustrate our method on an obstacle avoidance task, showing the pushforward of À at the optimized control points and the resulting continuous path of probability densities.
The optimization of the K + 2 control points proceeds in the following steps: Step 1: Initialize pair of boundary parameters, Og and 0,, that accurately approximate the target boundary densities
DuleA
22
Dilet
(b) Control points
Pazester
Sumple
(a) The space of probability densities P(Bd) (top), the admissible parameter set θ C Rd (left), and the sampling space Idd with an obstacle (right)
6.50
0.5
(c) Density path
Figure Visualization of our framework. (a) Illustrates the three main components: the space of probability densities ((Rd) (top), the admissible parameter space θ C Rd (left), and the sampling space I containing an obstacle (right). (b) and (c) depict the pushforwards of the optimized control points (To, ) for O, ...,, and the continuous spline trajectory (To(e) ) X, respectively. Time is color-coded, with red indicating it = 0 and yellow indicating =

---

Po and P1. Step 2: With the boundary parameters (6o, 01) fixed, optimize the interior control points (0t,)TT to minimize the action associated with the coupling (To,, To, 1. Step 3: Optimize the boundary parameters (fo, f1) in order to improve the quality of the coupling while keeping the interior control points fixed. Repeat steps 2 and 3 iteratively until convergence.
Contributions. Our main contributions are: (1) We introduce Parametric Density Path Optimization (PDPO), novel framework that transforms dimensional problem by parameterizing pushforward maps with cubic Hermite splines. Our approach achieves accurate density paths with as few as 3-5 control points. (2) We generalize the parametric framework to support a wide class of oetion-minimiin control (via Fisher information), pcceleration-regulee within a unified optimization scheme based on learned transport maps. (3) We empirically validate PDPO across several challenging benchmarks -iniuurkk-iicuudinng obstacle avoidance, cntropy-constrainn transport, and -inensninnii opinion dynamics achieving up to 7% lower action, up 10x improved boundary accuracy, and 40-80% faster runtimes compared to state-of-the-arr baselines.
density path optimization into a finite-
problems. including stochastic optimal dynamics, and mean-field interactions,
Background
This section presents the mathematical foundations of optimal transport, its dynamic formulation, and parametric methods that enable scalable and efficient application in machine learning contexts.
2.1 Optimal transport and Wasserstein distances
In what follows, we only consider distributions μ that are Lebesgue dominated, thus there is function p Fd - Rzo such that for any measurable set A, μ(A) = JA p(x)dz. Since are restricting ourselves to this case, we use the terms distribution and density interchangeably for convenience.
The Wasserstein-2 distance between two probability distributions Po, P1 E P((RR) is defined as:
inf / ||
(2 PI)
xEII(po.PL) R
where 1l(po, P1) denotes the set all joint distributions with marginals po and pt 33 When is absolutely continuous with respect to the Lebesgue measure, the Benamou-Brenier formulation 4 provides dynamic reinterpretation of this distance:
W2 pi) innf inf ('' | P.v
1 [llo(t, x1) p(,, (((t))d))
(4)
subject to 0.p V (pv) = 0, p(0,-) -0,-)=Po. Po:
p(1,-) = Pt.
(5)
This dynamical formulation interprets the Wasserstein distance as the minimum kinetic energy required continuously transport the mass of po into Pt over the time interval [0, 1].
2.2 Generalized action-minimizing problems in Wasserstein space
We consider the following generalized action-minimizing problem:
F(p(t)) dt inf 1 |lu ((),,
1 ((t, p(t,)))
(6)
subject to 3,p V (pv) = 0, p(0.-) -) = po:
p(1,.) = Pl,
where the functional F(p) incorporates additional constraints or potentials, defined as:
F(p) = Kin / ((x))(x)da + Kn / U(p(z))ds + K2
W(x - -))((r))())ddda
Here, the three terms correspond respectively to: external potentials (e.g., obstacles or environmental forces), internal energy (e.g., entropy or Fisher information), and interaction energy (e-g., repulsion or attraction between mass clements). The constants Ko, K1, K2 E 14, model the strength of each term.

---

In our setup, the Fisher Information potential, given by
04 /
FI(p)
| log pll p(x)dx,
(8)
plays a crucial role. As shown in prior work [1] 23 and detailed in Appendix A that the following deterministic control problem:
is well established
! [ inf Jo g 2 lltt,, (1),,, z)dr + F(p(t)) + FI(p(t)) 0:0
subject
0.p+ (pe) 0, p(0.-) po. p(l,-) Pl.
is equivalent to the stochastic optimal control (SOC) problem:
inf [. J P:
1 2 llu 2)1 x) F(p(t)) dt
2 Ap. p(0,-) po. p(l,-)
subject to 3,p + (pv)
= PL
This equivalence allows us to recast SOC problems as deterministic citiinnmiiimzzation tasks augmented with Fisher Information regularization term. While the deterministic action-ninimization formulation offers a principled framework for modeling action-minimizing dynamics, directly solving the resulting nnfnite-dimensionall optimization problem over probability densities is computationally prohibitive in high dimensions. To address this, we next introduce parametric representations that reduce the problem a finite-dimensional setting, enabling tractable and scalable computation.
2.3 Parametric pushforward representations of probability densities
To render action-minimizing problems computationally tractable, we focus on parameterized families of probability distributions, thereby reducing the original infinite-dimensionall optimization problem to finite-dimensional one over the parameter space. This strategy builds on recent developments in parameterized Wasserstein dynamics for solving initial-value Hamiltonian systems and gradient flows [5].
The work most closely related to our approach is that of [8]. which formulates the Schrödinger Bridge (SB) problem between Gaussian distributions as an action- action- minimization task. Leveraging the fact that the SB solution between Gaussians remains Gaussian, they parameterize the evolving density by its mean and covariance, and derive Euler-Lagrange equations for these parameters by projecting the density-level optimality conditions onto the parametric space.
We now introduce the core elements of parametric pushforward formulation, starting with the definition of some key concepts. For further background, we refer the reader to [6]
Definition 1. We consider a parameter space θ C IR. of dimension D. We call T IRd x x Bd a parametric mapping, for every E O, the function T(,, A), or To(.) used interchangeably. is Lebesgue measurable. The corresponding parametric density space is defined as
Po: (po (To) e) P(R.)).
where (To) # denotes the pushforward of a fixed reference density X under the map Te, see Fig- ure top). Moreover, given a curve in the parameter space 0(t) E O indicated by E [0, 1]. the corresponding pushforward density path is defined as Pa(t) == (To(t))
# (
see Figure Ta top), and
the sample trajectories are given by (z) (Ttt)22) X, see Figure Tal right)
3 Action Minimizing Problems via Parametric Pushforwards
3.1 Projection of the action-minimizing problem to the parametric setting
To motivate our framework, we begin by revisiting an alternative formulation of the classical Wasserstein-2 distance defined 33. Given x, y IRd, define the set of admissible paths
[(x, w) = (n( C²((0, 1), Rd) 7(0) x, x11) = w}. Then the classical formula
|| y/ = inf f( 1I dt)

---

allows us to write the Wasserstein-2 distance as
J (t) dt dx v). xell(po.P)) Jdxxd (m))
inf [
(Po
inf
(9)
In other words, one first chooses a coupling æ between the endpoints and then pays, for each pair (2r, y), the minimal kinetic energy needed drive particle from æ to y- Our framework adopts exactly this viewpoint: we cast the dynamic optimization in (6) as a static coupling problem whose cost still evaluated via dynamic (least-action) criterion.
For generalized minimiinini problems, the action of a path is
[6'
d di )
A(r) Egnipo
dt +
)(.
(10)
We make the following restrictions guarantee the existence of the density in all the trajectories of coupling, Define the set
Tpo.pn ( a diffeomorphic, T#( p,, det( > 0}.
Since the Po and P1 have densities, the Monge map 7 is diffeomorphism that satisfies det (DT) 4 and the set pp # 0. If E To. then it is smoothly isotopic to the identity id, see Lemma Chapter 6 77. By the definition of diffeomorphism being smoothly isotopic, for maps T E To the
set (Tr) consisting of the flow maps C¹ ( 1] x R. ).), such that ()(x)) == ((t, x)) describes the position at time t of a particle starting at æ, ie.. (x(x)) = x, 7e(-) diffeomorphic for each
t [0, 1], and (10: 71) Po is non-empty. We also define the set of couplings
i = t E U0 p1) ग = (id, Tx) (po) for some map Ta E TTo.pp1)
accumulated potential over the evolving density path. The cost associated with coupling E II then defined as
<(T) = inf A(7). (ET(a)
Our primary objective is the following action-minimizing problem:
c(TT) inf inf A(7). =ef((p.P1))
inf
(11)
=EN((...P)) 7E (m)
To the best of our knowledge, this formulation is novel. While related ideas appear in [29 and (32] those works are limited to linear potentials or "lifted" Lagrangians. In contrast, our formulation includes internal and interaction energy terms, which significantly alter the structure of the problem. In particular, the cost of each particle trajectory depends not only on its velocity and position, but also on the global density path P:.
The following theorem, whose proof can be found in Appendix B establishes the equivalence between the dynamic and static formulations given in Equations (6) and respectively.
Theorem 1. The dynamic formulation in Equation is equivalent to the static formulation in Equation III.
With these definitions in place, now formulate the action-minimizing problem from parametric perspective. Let the parameter space θ the parametric map T, and the reference density X be fixed. Define the of admissible boundary parameters as
+ = (00, 01) E R RD.
(To.) p=;
Given parameter curve with boundary condition fo--1 O0-11 := {f0--1 E C[0 A(0) fo. (11 = 0,), we define its action (by slight abuse of notation) as
A( = EE-a [/')
d dt
(owo
at F((T(())) A)dt.
(12)
5

---

where (o() denotes the transported trajectory at t of a particle 2 drawn from the reference
distribution with density X. The corresponding parametric retion-minimizattin problem is given by:
inf
inf 00-1E00-1
A(00-+1)
(13)
(0c.0)
This formulation parallels the coupling-band problem in Equation but operates directly in the parameter space of the transport maps.
While [16] explored the Lagrangian formulation of Equation (6 in parameter space from theoretical perspective, our approach in Equation (33) introduces an alternative formulation based on coupling costs. This perspective not only offers conceptual clarity but also lends itself naturally to numerical implementation.
Motivated by (20), employ cubic Hermite spline approximation of the parameter curve 6 in
solving the optimization problem (13). The following theorem provides an estimate for the error introduced in the objective action. The full set of assumptions and proof are provided in Appendix C
Theorem 2. Let 00-+1 be a parameter curve with i continuous derivatives, and ler 0---1 be its cubic Hermite spline approximation using K uniformly spaced control (collocation) points. Then the error in the action satisfies A[00-=1] - A[0o-n1] = =Ohh=11 where K. min(1, 4).
While action iisiionmiiimizzing curves can computed analytically or via numerical integration for simple parametric families (c.g., affine maps) and linear potentials (see [33]). such formulations are inherently limited in expressivity. To overcome this, we adopt Neural ODEs as flexible and expressive parametric maps for modeling complex density dynamics.
3.2 Pushforward map via Neural ODE
Neural Ordinary Differential Equations (Neural ODEs) 9 define continuous-time transformations via parameterized differential equations:
(T) vo(T, (TT, (00 2)
(14)
dr
where vo [0, Rdd Rd is neural network parameterized by is typically standard normal distribution, and 7 denotes the ODE time variable. The resulting transport map is defined as:
T= = = L/'' co(T, ((T,)b(T)))))) dr.
(15)
We adopt Neural ODEs as parametric maps due to their ability to support efficient computation of both the entropy log(p) and the score V log(p) through auxiliary ODE systems (35). Details of these complementary entropy and score equations are provided in Appendix
3.3 Optimization framework
The static formulation 13) defines bilevel optimization problem: the boundary conditions ensure feasibility by satisfying (To,) pi, i 0, 1. The curve f0-11 minimizes the action and deter- mines the coupling cost. From the perspective of the spline control points (O,, J+ K++ the boundary parameters are On, 0,, while the interior path is defined by {0c, [K This motivates our optimization strategy:
Step 1: Initialization. First, initialize the boundary parameters 8n and 01 such that (Te.)s À are good approximations of p; for O, 1. We employ Flow Matching (FM) for its efficiency. We in general consider this step outside of the main algorithm. In Table we include in the last column the FM training times. A parameter can be used for any functional, so once a boundary parameter has been obtained, it can be reused across multiple problems, see Appendix F.1.2 Second, initialize the control points {0c. 1==1 as the linear interpolation of fg and 0j. Then run few iterations of path optimization with F(p) = 0. In Appendix F.I.I and Appendix F2 we show the computational advantages each part of the initialization step.
Step Path optimization. Perform gradient-based path optimization update the interior control points {0e, === to minimize the action (12), while keeping the boundary parameters fixed, see details
K

---

in Algorithm 1 The action functional T2 is evaluated using the trapezoidal rule with N time steps and M Monte Carlo samples. The gradient is computed using automatic differentiation. See Algorithm|
Step 3: Boundary/coupling optimization. Optimize the boundary parameters (0o, 02) while keeping the interior spline control points [0,,,,, fixed, see details in Algorithm 2 Specifically, we formulate
this as a penalized optimization problem:
Σ Lo.L((TT) pi) + A(00-11- 1=0
EEZA Po.
min E
(16)
Here L is dissimilarity metric between distributions, for which we use the FM loss function. The weights a;, = 0, I, balance boundary accuracy and coupling optimality, The action A(0 (12 is
added as penalty term, which is estimated at the current control points, see Algorithm 2
We repeat the optimizations of step and step alternatively until convergence, producing both an optimal coupling and a minimal-action path between Po and P1. see Algorithm [] Observe as long as the action lower bounded, step will converge to a local minimum. The dissimilarity assumed to be nonnegative, thus step also converges to local minimum.
4 Experiments
We conduct comprehensive evaluation of our method across a variety benchmark scenarios to demonstrate its accuracy, efficiency, flexibility, and scalability. In all experiments, we compare against GSBM (20). the current state-of-the-are which has been shown to outperform earlier approaches. Where appropriate. we also include comparisons with Neural Lagrangian Optimal Transport (NLOT) important to note that NLOT is limited to linear potentials, while APAC-Net is specifically designed for mean-field games (MFG) and employs soft terminal costs instead of enforcing hard distributional constraints.
29 and APAC-Net [17].
In
E.3
we include an ablation study on the number of control points.
4.1
Implementation details
We implement our method in PyTorch 28 and run all experiments on an AMD 7543 CPU + NVIDIA RTX A6000 GPU. The Wasserstein-2 distance is approximated using the POT library [12) with 3,000 samples. For estimating the action, use Monte Carlo integration with 3,000 samples and compute the time integral using the trapezoidal rule with a step size of At = 1/50. Each experiment is repeated across three random trials with different seeds. Additional implementation and experimental details are provided in Appendix F In Table see the boundary conditions pu and p1, and Table the for the PDPO algorithm.
4.2 Obstacle avoidance with mean-field interactions
We demonstrate the effectiveness of our method through three challenging obstacle avoidance scenarios involving mean-field interactions. See Appendix F.4 for their mathematical definitions.
S-Curve with Congestion (SCC) from [17 Particles navigate around two obstacles arranged in an "S" "S"-shaped configuration while interacting with one another. PDPO achieves superior performance compared to existing methods (GSBM and APAC-Net), yielding lower action values, more accurate boundary approximation, and faster runtimes (see Table
Notably, the interpolated density path generated by PDPO (Figure 2b exhibits clear nonlinearity between control points. As shown in Figure 2a a linear interpolation of the density between control points 0r, (red) and 0r, (orange) would intersect the obstacles. This highlights how PDPO's spline- based parameterization effectively exploits the geometry of the parameter space to avoid infeasible paths. Note that both GSBM and APAC-Net violate the obstacle avoidance constraint in this scenario.
V-Neck with Entropy and Fisher Information (VNEFT) from 20) This stochastic optimal control (SOC) problem involves particles navigating narrow channel while minimizing entropy, see Figure [Tb PDPO outperforms GSBM by achieving lower action values, comparable boundary accuracy, and significantly faster runtime (see Table I]. When GSBM is restricted to the same runtime as
7

---

Table 1: Comparison for obstacle avoidance with mean- meen-fielld interactions. The quantities (2(p:)
i = 0, 1 denote the empirical are the Wasserstein-2 distances at times computed using the POT library
12 The FM time the sum of the FM training times for both boundaries.
Problem (method)
A(pe]
Su E.p, V(X)dtt (W2(Po), W2(p:))
Time
FM Time
SCC (PDPO)
39.90 + 0.16
1.12 + 0.02
(0.028, 0.014)
11m
134s
SCC (GSBM)
40.02 .75
4.05 = 1.48
(0.21, 0.25)
21m
SCC (Apac-net)
39.95 0.25
1.31
(0, 0.,03)
12m
VNEFI (PDPO)
148.91 1.71
38.290 = 0.43
(0.082, 0.078)
70m
77s
VNEFI (GSBM)
156.80 4.31
74.28 1.23
(0.086, 0.062)
115m
GMM (PDPO)
79.46 + 0.62
0,41 0.35
(0, 0.359)
8m
3m10s
GMM (GSBM)
88.77 3.31
14.59 3.26
(0.806, 0.258)
123m
GMM (NLOT)
82.06 + 2.80
0.006 + 0.006
(0, 0.629)
117m
Opteeeed MBA355 Ke 1p0
12
S.
400
160
2
16
0.25
0.33
0.75
100
0.5
0.25
3.55
8.00
025
0.55
0.75
1.0
0.5
0.35
0.5
075
100
(a) (To.,)
X
(b) PDPO (To,) X
(c) GSBM
(d) APAC-Net
Figure 2: S-curve-C: C: (a) Pushforward densities at control points, (b) Pushforward densities along the interpolated trajectory, (c) Density path produced by GSBM, (d) Density path produced by APAC-Net
PDPO (details in Appendix F.4 its action value remains similar, but its boundary approximation deteriorates by an order of magnitude.
Gaussian Mixture obstacle (GMM) from 20 In this benchmark, particles must move between multimodal source and target distributions while avoiding a Gaussian mixture-shaped obstacle, see Figure 13 for the density path in the appendix. This example demonstrates PDPO's flexibilityy in selecting different reference distributions. By setting the reference distribution = Po. there no approximation eтror at time t 0. In Table [ we compare against GSBM and NLOT. Our approach achieves slight improvement of 0. 0.04% in the action, but the computational time was reduced by approximately 85% from nearly 2 hours to less than 10 minutes. While GSBM achieves marginally better boundary representation for Pl, our method offers comparable optimization quality with substantially faster computation.
Note that for GSBM, the reported (2(P)) values are computed by evolving samples drawn from P1 via the backward SDE to obtain the empirical approximation of
4.3 Generalized momentum minimization
Building on the momentum Schrödinger Bridge (SB) framework by 10 we extend action- minimizing problems to incorporate acceleration alongside additional potential terms. This leads to novel formulation that minimizes the following objective:
A(7) En-ro [
mas di + F(er.)dl.
Figure illlstratee PDPO solutions for the obstacle avoidance with mean-fiell feend interaction class of problems. The resulting paths are noticeably smoother, compare the top and bottom of Figure
To our knowledge, no prior method has addressed these generalized settings within a unified momen- tum formulation, underscoring PDPO's flexibility beyond classical Wasserstein geodesics.

---

Optim bed deraty path
Derenton
Ommman
162
2.1
5.00
6a
68
Time 02
Timm 11
(a)
(b)
(c)
Optimtze demtyy 24m
107
retinteee sendy outh
nized draaty a3mh
100
6.4
5.0
Drmention
5-8
2.
L4
Omeraon
5.04
4.32
180
-25
-50
1.5
211459
O kn
1.44
0003
0.7
100
9:43
00:3
1:0
D.2
03.3
100
6.0
(d)
(c)
Figure 3: Solutions for generalized momentum-minimizeee problems (a-c). Solutions for kinetic energy minimizing problems (d-f).
4.4 High-dimensions opinion depolarization
We examine the opinion depolarization problem following 33) and [19]. Here, p(t,-) RR Rzo
models the opinion density of population, where each coordinate represents an agent's opinion on topic. The particle-level opinion dynamics evolve as
dxx(t) fpolarize pzt)); -)) dl polaaiie (z(();;((1,))))
with
Fpotatize (z: p(t, Ey-p, y. a(x, = {
if sign (t) sign (y, Ex))
otherwise
The function a(x, y: {) evaluates opinion alignment at random information & sampled independently of p(t, -). These dynamics cause unpolarized initial condition p(0,-) (e.g.. Gaussian) to polarize into opinion clusters. We seek velocity a(L, x) preventing polarization by enforcing an unpolarized terminal condition p( (see Appendix F.5).
The action is
d di (ooo(())1
+ [''' F((To() a (A))
Apolarize (0o-.1) ;= BENA
FigureAshows PDPO achieves quality comparable to [20]. as evidenced by 2D PCA projections
and directional similarity histograms terminal time. These histograms show cosine similarity distributions between opinion vector pairs, measuring alignment independent of magnitude. PDPO's key advantage computational efficiency: completing in 23 minutes versus over 5 hours for GSBM. We attribute this success to: (I) the spline path 0,-, readily reaches the non-polarized target at t from initialization, and (2) connects boundary conditions independently of folaaiie- Note that training boundary conditions 6g and 0 required 29m35s and mm12s, respectively.
5 Computational complexity and limitations
The primary computational cost our algorithm arises from pushforward evaluations. When action integral is approximated using N time points, our method requires evaluating N Neural ODEs (NODEs). Each NODE involves integrating velocity field using 10 time steps. For expectation estimates, use M Monte Carlo samples per NODE.

---

20
-10
-10
-10
-20
20
20
20
20
-20
Dissiilanty plot
Dissimlianty piot
Dissiilianty plot
2000
1000
1000
200
800
400
400
400
200
200
Figure 4: Opinion depolarization in 1000 dimensions. Top: 2D PCA projections of the distributions. Bottom: directional similarity histograms. Left: target unimodal distribution. Middle: PDPO solution.
Right: DeepGSB solution.
When the functional F includes entropy and Fisher Information terms, the cost increases substantially. Computing the entropy requires solving additional ID ODE, while the Fisher Information involves a separate d-dimensional ODE, where d is the problem dimension. In total, this results in O(N x M x d)) ODE evaluations for problems involving both terms.
In Appendix F.6
we showcase the scalability of our method with a 50-dimensional Schrödinger Bridge example, representing the fully computed.
case with Fisher Information we have success-
Currently, our implementation is limited to NODEs with MLP architectures, see Appendix which limits applications to image-based transport problems, unlike GSBM [20]
Conclusions
We introduced Parametric Density Path Optimization (PDPO), general framework for solving oetion-minimizn problems over probability densities by leveraging parameterized pushforward maps. This approach transforms the original -iinnte-dimmnsionall optimization over density paths into a tractable, finite-dimensional problem in parameter space, using cubic splines with a small number of control points.
central conceptual distinction of PDPO is that its outcome a single trained neural network, but collection of K + trained models (To,, JKK whose parameters lie on spline-defined trajectory. Together, these models define a continuous transformation in density space through the pushforward of the interpolated parameters. From this perspective, PDPO learns dynamical path in the parameter space, a time-dependent family of mappings that approximate the evolution of the probability mass.
Experimental results show that PDPO matches or surpasses state-of-ther the-art methods in accuracy, while requiring substantially less computation time.
7 Broader Impact
The study of action-minimizing problems in density space might have applications in the study of population dynamics.
Acknowledgments and Disclosure of Funding
Chen acknowledges partial support from NSF #2325631, #2245111. This research is partially supported by NSF grant DMS-230746. The authors would also like to acknowledge Benjamin Burns for the helpful discussions and the rendering of Figure and Shu Liu for the helpful discussions.

---

References
[1] Yves Achdou and Italo Capuzzo-Dolcetta. Mean Field Games: Numerical Methods. SIAM Journal on Numerical Analysis, 48(3):113--116,, January 2010. ISSN 0036-1429, 7170- doi: 10.1137/090758477. URL http:://ppbss: htt.ttp://eegeeeee
siam. org/ doi/ 1t10//00.1177700754477777
[2] Michael S. Albergo, Nicholas M. Boffi, and Eric Vanden- Stochastic Interpooaas A Unifying Framework for Flows and Diffusions, March 2023. URL htttp//arxiv..or moond-mat..
2303. 08797 MXXKv:2203.08797
[3] Brandon Amos, Lei Xu, and J. Zico Kolter. Input Convex Neural Networks. In Proceedings of the 34th International Conference on Machine Learning, pages 146-155. PMLR, July 2017. mlr. press/ v70/ amos17b. html
URL httpt:ttttppppeeeeeeeoeeeeeeeeee
ISSN: 2640-3498.
[4] Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to the Monge-Kantoroon mass transfer problem. Numerische Mathematik, 84(3):375-393 January 2000. ISSN 0945-3245. doi: 10. 10. (0000000001000000000000 URL littps://dot ort.:/...../000..0..........
5002110050002
[5] Jean-David Benamou and Guillaume Carlier. Augmented Lagrangian Methods for Transport Op- timization, Mean Field Games and Degenerate Elliptic Equations. Journal of Optimization The- ory and Applications, 167(1);1 October 2015. ISSN 0022-3239.
doi: 10. 1007/ (00557011--075-- URL http::/linkk springer com/ 10. 10100 072-
[6] Joe Benton, George Deligiannidis, and Arnaud Doucet. Error bounds for flow matching methods. Transactions on Machine Learning Research, 2024. ISSN 2835-8856. URL https
[7] E. H. Brown and John W. Milnor. Topology from the Differentiable Viewpoint. The American Mathematical Monthly, 74(4):461, April 1967. ISSN 00029890. doi: 10.2307/2314613. URL
https://www. jstor org/ sttssssaaleeeeeeee 21111111101111111rrroorr
[8] Charlotte Bunne, Ya-Ping Hsieh, Marco Cuturi, and Andreas Krause, The Schrödinger Bridge between Gaussian Measures has a Closed Form. In Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, pages 5802-5833. PMLR, April 2023. URL
attps://proceedingg: mlr https//v0//buunneee htmll
html
ISSN: 2640-3498.
[9] Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural Ordinary Differential Equations.
[10] Tianrong Chen, Guan-Horg Liu, Molei Tao, and Evangelos Theodorou. Deep Momentum Multi-Marginal Schrödinger Bridge. November 2023. URL https://openreview net/
forum?id-ykvvvogeee
[11] Yongxin Chen, Tryphon T. Georgiou, and Michele Pavon. On the Relation Between Opti- mal Transport and Schrödinger Bridges: A Stochastic Control Viewpoint. Journal of Op-
timization Theory and Applications, 169(2):671-691, May 2016. ISSN 1573-2878. doi: 00.1007//10957--055008033 URL https://tttt...... wtt.:oi..rgg1...1 1007 / http://11..100/000000010------00 -109977-0550000--2
[12) Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gau- theron, Nathalie T. H. Gayraud, Hicham Janati, Alain Rakotomamonjy, levgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer. POT: Python Optimal Transport. Journal of Machine Learning Research, 22(78):1-8, 2021. ISSN 1533-7928. URL http: //jmlr org/ wwwpg//appra//222. 20111110-411111 html
[13] Han Huang, Jiajia Yu, Jie Chen, and Rongjie Lai. Bridging mean-field games and nor- malizing flows with trajectory regularization. Journal of Computational Physics, 487: 112155, August 2023. ISSN 00211-9911. doi: 10. 10.1016//.jep.p 20... URL hhtps www.:sccencediree.. com ott.apace....ctttaaooa.aaooaeaea
[14] Matt Jacobs, Flavien Léger, Wuchen Li, and Stanley Osher. Solving Large-Scale Optimization Problems with a Convergence Rate Independent of Grid Size. SIAM Journal on Numerical Anal- ysis, 57(3): 100-1123, January 2019. ISSN 0036-1429, 1095-7170. doi: 10. 10.1137/18M1 18640X.
URL /tattps://epuute siam. org/ wwt.://......11111111

---

[15] Yijie Jin, Shu Liu, Hao Wu, Xiaojing Ye, and Haomin Zhou. Parameterized Wasserstein gradient flow. Journal of Computational Physies, sciencedirect. com http://ccceccccccarrrriiiooooo. 10. jep.244..... p:1/$0021999124000022
March 2025. ISSN 0021-9991, doi:
URL httpt::///// /wwp:
[16] Wuchen Li and Stanley Osher. Constrained dynamical optimal transport and its Lagrangian for- mulation, September 2018. BRLphttp:/aaxxi hwrggabb/18007 00937 arXiv:1807.0099 [math].
[17] Alex Tong Lin, Samy Wu Fung, Wuchen Li, Levon Nurbekyan, and Stanley J. Osher. Alternating the population and control neural networks solve high-dimensional stochastic mean-field games. Proceedings of the National Academy of Sciences, 1 118(3 xs2224731122 August
2021. doi: 111003/pnas.20247111
URL BtLhttpp: //wwww pnws. pnas wtg//oi///..oo/oo 10t/://.....0//ppp Publisher: Proceedings of the National Academy of Sciences.
2024713118
[18] Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow Matching for Generative Modeling. September 2022. URL https://openrreview.cc net/
forum/1d=Pqy/RDccJJge
[19] Guan- Liu, Tianrong Chen, Oswin So, and Evangelos Theodorou. Deep General- ized Schrödinger Bridge. October 2022. URL https://openrerievv.ne wttpr/iwwwwe//foumm/.e.m fp33Nsh005
[20] Guan-Horng Liu, Yaron Lipman, Maximilian Nickel, Brian Karrer, Evangelos A. Theodorou, and Ricky T. Q. Chen. Generalized Schrödinger Bridge Matching. arXiv, April 2024. URL
nttp://arxive wwtpg////////o/eoooooo 02233 arXiv:2310.02233.
[21] Shu Liu, Wuchen Li, Hongyuan Zha, and Haomin Zhou. Neural Parametric Fokker-Planck Equation. SIAM Journal on Numerical Analysis, 6085):1385-111 1148885--44,,0-11 June 2022. ISSN 0036-1429, doi: 10. 10.1133/20 1449260 344986. URL https: //epsbs siam. org 113/2/01441// Publisher: Society for Industrial and Applied Mathematics.
[22] Shu Liu, Stanley Osher, and Wuchen Li. A Natural Primal-Dual Hybrid Gradient Method for Adversarial Neural Network Training on Solving Partial Differential Equations, December 2024.
URL httt:://aaxiiii http://wwwxvrrrr///////// Setpri/worrg/bb/221..553 06278 arXiv:241 120224111.6627 1.06278
[math].
[23] Flavien Léger and Wuchen Li. Hopf-Cole transformation via generalized Schrk"odinger bridge arXiv:1901.09053
problem, January 2019. URL [math].
[24] Toshio Mikami. Stochastic optimal transport revisited. SN Partial Differential Equations and Applications, 2(1):5, February 2021. ISSN 2662-2963, 2662-2971. doi: 10.1007/ -42985--000000559-3 URI. http://arxiv org/ abs//00/.. 2003. 11811 arXiv:2003.1181 [math].
[25] Derek Onken, Samy Wu Fung, Xingjian Li, and Lars Ruthotto. OT-Flow: Fast and Accurate Con- tinuous Normalizing Flows via Optimal Transport. Proceedings of the AAAI Conference on Ar- tificial Intelligence, 35(10):9223-95) 9232, May 2021. ISSN 2374-3468. doi: 10. 100116/iiaai.3355 17113. URL https///::///je 10.
aaai. org/ index. phtp/AAt//rrtccle/viieeee
Number:
[26] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang. Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An Imperative Style. High- Performerfom Deep Learning Library, December 2019. URL http:://aaxxxi org/ abs/ 1912.01703 1rXiv:1912.011 [cs].
[27] Gabriel Peyré and Marco Cuturi. Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning, 11(5-6):355 607, 2019. ISSN 1935-8237. doi: 10.1561/220000073. http:tt:://aee http:///xdd.iiii/o//ooo////....t org 111000001111222000000000

---

[28] Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, Jinkyoo Park, and Stefano Ermon, Torchdyn: Implicit models and neural numerical methods pytorch.
[29]
Aram-Alexandre Pooladian, Carles Domingo-Enrich, Ricky T. Q. Chen, and Brandon Amos. Neural Optimal Transport with Lagrangian Costs. June 2024. URL httts:///ppenree
pet//forun?1d=xx4pa22sssy
[30] Simon Schweighofer, David Garcia, and Frank Schweitzer. An agent-based model of multi- dimensional opinion dynamics and opinion alignment. Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(9):093139, September 2020. ISSN 1054-1500.
doi: 10.0665000007233.
URL / wtt.:////...111100000000000000
[31] Alexander Tong, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector-Brooks, Kilian Fatras, Guy Wolf, and Yoshua Bengio. Improving and generalizing bosed generative models with minibatch optimal transport, October 2023. URLthttp://rree /ttp://///xiv..//// www.ii.og//a........ooooo wrgg/bs/////0//
00482 arXiv:2302.00482 [cs].
[32] Cédric Villani. Optimal Transport, volume 338 of Grundlehren der mathematischen Wis- senschaften. Springer, Berlin, Heidelberg, 2009, ISBN 978-3-540-71049-3 978-3-540- 71050-9, doi: 10. 111100/998-8--50-7005----9 URL http:://link.sprinngg compeer.commm cwm/10. 1007/
778-3-540-71050--
[33] Hao Wu, Shu Liu, Xiaojing Ye, and Haomin Zhou, Parameterized Wasserstein Hamiltonian Flow. SIAM Journal Numerical Analysis, 63(1):360 395, February 2025. ISSN 0036- 1429, 1095-7170. doi: 10.1137/33M15928IXX URL httph::t///eeee siam org/ 1137/
23M159281X
[34] Yujia Xie, Minshuo Chen, Haoming Jiang, Tuo Zhao, and Hongyuan Zha. On Scalable and Efficient Computation of Large Scale Optimal Transport. April 2019. URL https
[35] Mo Zhou, Stanley Osher, and Wuchen Li. Score-based Neural Ordinary Differential Equations for Computing Mean Field Control Problems, January 2025. URL battp://arxta org/ /b// 2409. 16471 arXiv:2409.1647 [math].

---

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

---

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

---

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

---

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

---

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

---

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

---

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

---

0, + V (p) 0,p + (pu) = o² Ap, A Deterministic and SOC problems
We now show the equivalence of
inf ( / (1le(t,)))) (() F(o(t)) + FI(((()) P:
(17)
subject to
Oep V (po)= =0, p(0,-) po, p(1,-) pl,
and
J(¹ See 1 (1), (((,))/ F(e(!))
inf (lt.r).(t..)))
(18)
a² subject Oep (pn) Ap, p(o,-) = po: p(1,-) = p1-
Proof. In what follows, we use the notation 6, for the L gradient.
Let S be a Lagrange multiplier for (77. then KKT conditions for 17 are given by
VS P. a.c.
0,5 7 / Fp)F(p) F,FI(D)
0,p V (pv) O.
Define & =S+ o2/22og/ and : VΦ, p a.e. Then the pair (p, u) satisfies
" VΦ,
a.e
2
which is in turn the KKT conditions of the optimization problem
B Equivalence of dynamic and static frameworks
Here we show the equivalence between
inf ('' / |lu(t, ll((, =)) + F(p(t)) dt 8.0
subject to Oep + v (po) = 0, p(0.-) = Pa. p(1,.) Pr.
and
c(TT) inf inf A(7). xefl(p..1))
inf
=Ef)(p..P)) (
Proof.
Observe that local minimizer of 66) (p, v) defines a C C'-diffeemorphic curve 7a,v as the flow
= x [' v(s, 7(s, xs),
and the coupling p = (a, x). The coupling up.u and path Np are an admissible pair for (II) From the definition of Np..: follows that ()).((t) (Tp. # (po) p, and
S. [ (((1,x))PP(() ((t, xddxdt = A(Tan):
then S. I 2 |o(t, x)) dd + F(P) dt A(tinw) + { F(py,)dt.

---

Thus 66) upper bounds I1 Now, given local minimizer of [11 (T, 7), the density Pn,n(t,)) (2(()) (Po) and the velocity field
d dt (t, x(0,z) Pu.
((t,
define admissible pair of (6. As before,
A(7) S. F(py)dt J/¹ / [l (r + F(Pr,y)dt.
Thus, 11 upper bounds
C Wasserstein error bounds for interpolation
Theorem Assume that T Lipschitz in its second variable with constant C, the sense, Zz Hd and any 0, 0 RDD
((12,))T((z))) $ C 0||
then
W2(po. Pb) - C Ø||
(19)
Proof.
From the static definition of the Wasserstein distance Equation (3) have
W2(po. P8) = (/ 0) - T(2, -((,
C||@ 011
We recognize that T being Lipschitz in the above sense is strong, We leave for future research, relaxing it to locally Lipschitz, a more feasible condition.
Let t0-11 and to-s1 denote two curves in the parameter space. Following the hypothesis on Theorem it follows that Vt [0, 1]
Corollary 1.
W2(Pore Pāt))) I C((((t) - (()
Corollary 2. Let 00-+1 be curve in parameter space with continuous derivatives up to order 2 1, and let d--+1 be the piecewise cubic Hermite spline interpolation of 0(t) at K uniformly spaced points {t. = K+1
K+1 Define K = min(l, 4), then
100 - 0o--1l S Mr max (10())))) {E(0,1]
where h R4T and Ma is constant depending Furthermore for all E [0, 1), the Wasserstein distance between the density paths time is given by
With g(x), the K derivative of the path O(t).
Wa(Petn: Pa()) I CM2 max EE(0,1) []0(() h*"
where C is the Lipschitz constant from Theorem
Theorem 4. Let 00--1 be curve in parameter space with continuous derivatives up to order 1 2 2, and let o be the cubic Hermite spline interpolation ((() at K uniformly spaced points {4: k= Define == min(l,4).
K+1
Assume that:
1. F Lipschitz with respect the Wasserstein-2 distance with constant My.
2. [,,)((2,0]] < Mor constant Ma for all and 6
3. 0gT(z, 0) is Lipschitz continuous in 0 for all z with Lipschitz constant Ca.-

---

4. (10(+)]l < Mer.for some constant Ma for all E [0, 1].
Then the error in the action functional bounded by:
[A[0o-1] Al0o-+1]1 S CA max
((0())
EE(O.1)
where h K-+1 and CA is defined below at the end of the proof.
Proof. We split the proof into two parts, one part a bound for the kinetic energy term, and the other on a bound for the potential energy term. The bound for the potential energy follows from the Lipschitz continuity of F and Corollary
L (P(oe)) F(Pāts))ltt -
f."
My WP Po(t) )tt S MM M ht max [1(()(())
$e|0,1]
The bound for the kinetic energy can be derived by assumptions 2-4 as -
B d
d ti ( (+) dt
di
<
EEnA
(x))
- dt ores) (2) /|/
dt
E=~A [(11 d (2( (
d d = at Ro(o(()
+ dt 28
S Mar + Mar (..r1)))1((EE.ll1007):0000000000 -
= M + MM /
S MaaT(Mer + Maa) / Bans [11
+ (((6)T)(,,0()))((()) (100))00000))0000000000000000110
MM,T(Mer + Ma²) / Ex~> (1001/00000000 0'(1ll
+ +I((O,T)(2,000) - (2,T))
Ma + M& /'' -1) 0'())) + Co. - ((()))10((1)) (1)1) de
< MM + Ma)
Ca. Mch" Mar) max [[0())(E)].
(ElU,1)
combination of the above two bounds concludes with
CA max Max Myhh, MarT (M + Mg. + Ca, MgMMh))
We recognize that hypothesis on Theorem 4 which requires to be Lipschitz with respect to the W2 metric, is strong assumption for general functionals F. When F consists of linear potentials and bilinear interaction potentials whose corresponding ((x) and W (x--y are Lipschitz continuous with constants My and Mw. the functional F(p) = SEs ((z)p(z)dz + S W(x p)()z)(())dzdd)

---

is also Lipschitz. To show this, let v. P(RR) and Ty-v denote the Monge map between them, then
F(u) F(u)]= =
[Ex=x[VXX]] + Ex~uEyy (Xx
Y)I)
- Ex-al V(X) - V(I +
- -
2 (My + +2Mw))Wa((.))
Observe that the last inequality follows from using the Lipschitz condition on V and W, and the following application of Hölder inequality,
Ex=pIXX-T--o (X)Il -exxxlIX-XTii (xx = (Ex~ullX - (x)))(
D NODES
See 9 for the standard reference on NODEs. For the NODE
dT)(T, = vo(T, ()T,))), (000 2) = X,
dr
initialized at p(0, 2) NN 1), let p(T,,) = (T, ) then
d log(p(T, x)) (==((T,))) V
(20)
dr
and d(V, log p(T, x)) dr
<=)((;,)) - V(V vo ((T,))))
(Veo(T, (T, =) =)))" logglo, ()
(21)
with initial conditions log(p(0,,)) 2)) log(pv(>)), where PN((2) is the density of standard normal distribution, and Vz log (0, = -z respectively. Although is not necessary to use a Gaussian density for the initial condition, is common in the literature, these quantities might not be accessible for other densities.
The result for 20 is standard and we do not include the proof here. Although Equation 21 has been proved before, see e.g., (35], we provide a proof for completeness.
Proof. For this proof, assume p(T, x)) has continuous second-order derivatives with respect and æ. In what follows, we use the short hand notation & to indicate (b(T,2) 2), and v instead of vg.
d V. log(p) = 0 Va log(p) v log(p)
dr
V -
0 V log(p)
-V (pv) + V log(p)
= VV T V. v) + v2 log(p)
= =-Vee" v. vr.
The first equality follows from the dynamical system de from the continuity equation O.p -V (pv), and the fourth equality from p -'up V log(
,े v(T, (b, 2)), the third equality follows
Algritthss/mmleemntttion details
We two separate optimizers: one for the boundary conditions (0u, 01) (Algorithm 2 and another for the control points {0,, 1 (Algorithm [}. We use Adam optimizers [26] in both cases. For the learning rate scheduler, used StepLR cosine, with specifics provided for each experiment.
key implementation challenge is preserving PyTorch's computational graph. Specifically, the 'oodch.load' operation does not retain the graph for muddient-based optimization. This is particularly

---

problematic in Algorithm where the interpolated points f(t,) used in the trapezoidal integration must track the computational graph with the spline control points enable proper gradient flow during optimization. When 'torch.load' is used to assign models, this computational graph is lost, disconnecting the density path from the control point optimization.
To address this issue, we re-implemented the MLP architecture from scratch, passing weights and biases as explicit function arguments rather than storing them as internal model parameters, see Algorithm 4 While effective, this approach introduces limitations to architectural flexibility: testing new network architectures requires either hard-coding them developing alternative methods to preserve computational graphs when loading pre- pre-triind models. We leave overcoming this limitation future research.
Algorithm Path optimization
Require:
Parametric points 10 + potential function F, N 2 X, number of optimization steps Q1.
1: Initialize spline in parameter space with (0r K+1
2: for number of iterations Q1 do
3: Obtain {O(t)))} by evaluating the current spline at time steps bj NV
4: Evaluate the N + pushforwards (Te(t,) (2:)) N.M
==0,==1
5: if Entropy or Fisher Information then
6:
Obtain entropy Equation (20) score from
7: end if
8: Evaluate the Equation 12) by integrating time by trapezoidal rule and evaluating the
expectations via Monte Carlo.
9: Update weights (0c; } via gradient descent.
10: end for
return {0.
Algorithm 2 Coupling optimization
Require: Sampling function from Po. sampling function from P1, potential function F, control points {0.] J==n K+1 weights aig, ar, s X, number of optimization steps Q2-
1: for number of iterations Q2 do
2: Sample M points from { { Po, and { P1.
3: Evaluate -E
E L,T, (zx).
4: Obtain {0(t;)}} by evaluating the current spline at time steps tj = N-
5: Evaluate the N pushforwards (To(t,) (2x) N.M
(=x)
6:
if' Entropy or Fisher Information then
7:
Obtain entropy Equation (00) score from 21
8:
end if
9: Evaluate the Equation 12 by integrating in time by trapezoidal rule and evaluating the expectations via Monte Carlo.
10: Update do 0 by minimizing 16 via gradient descent.
end for
12:
return 0g, 01
E.1
Hyperparameter Sensitivity and Tuning
Neural ODE architecture:
The neural ODE model defines the parameter space and must be sufficiently expressive to capture the boundary distributions. Key considerations:
Architecture sizing: Use larger networks than minimally required for boundary fitting. For example, while [2,64.4] architecture can learn our V-neck boundaries accurately, we needed [2, 128,4] to capture the solution complexity when potential terms significantly alter the density shape.

---

Algorithm 3 PDPO
Require:
Sampling function for Pa. sampling function for pt. reference density À potential function F, number of control points K, weights coo, di initialized 0o, 01, number of total iterations Q, number of path optimization steps Q1, number of coupling optimization steps Q2, number of geodesic-warmup steps Q23-
1: Initialize b,, using an equispaced linear interpolation of 0o and 0
2: Run Q3 steps of geodesic warmup.
3: for number of iterations Q do
4: Sample .M points ( X.
5: Update {0., } by minimizing the action of the points { MM using Algorithm !
6: Update (0n, 01) using the points { in Algorithm 2
7: end for
8: return {0., K
Algorithm 4 Parameterized MLP Forward Pass
Require: Architecture arch = [d, wo I.] where d is input/output dimension, w is hidden width, is number of layers
Require: Input e [Raxd, parameter vector 0 E R'
Require: Time-varying flag time_ varriin E {true, false]
1: din < if time_varree varying d + else d
Adjust input dimension D Current position in
2: idx +-
3: h - I
4: // Imput Layer
5: Extract W1 E Ruxda from 0lidx idx + w din]
6: idx idx + dmn
7: Extract bj E R= from 0fidx idx + w]
8: idx idx +
9: h + oWW1 +b)
b a is activation function
// Hidden Layers
11: for L = 2 to
do
Extract Wi E Rxx from Olidx idx w2]
13: idx +- idx + w²
Extract by R" from 0iidx idx + w]
15: idx + idx + w
h h be)
17: end for
18: // Output Layer
19: Extract Wr. E Rdxw from Olidx idx + d w]
20: idx idx + d. w
21: Extract by E Id from 0fidx idx + d]
h Wy.h by
No activation on output
23: return
If the potential energy terms expected to create new modes or dramati- cally change boundary shapes, increase the hidden dimension.
Rule of thumb:
Neural ODE Integration: Two factors affect forward simulation quality:
Solver choice: We use midpoint integration, as it is standard in flow-based generative models.
Integration steps: 10 steps suffice for most problems. Use least the number of steps required by your boundary NODEs for accurate sampling from the target distributions.
Action Approximation (Time Steps N): This an important parameter that controls the temporal discretization accuracy.

---

Table 2: Action convergence with respect time discretization steps N.
N
Action
Kinetic Energy
Potential Energy
10
66.06
37.38
28.68
20
42.09
37.77
12.82
30
30.31
29,46
0.788
40
30.20
29.72
0.735
50
30.10
29.37
0.725
Recommendation: Start with N=20 minimum and increase until the action converges (typically N=30-50). The dramatic drop from N=10 N=20 demonstrates why sufficient discretization is crucial. This parameter is related to the complexity of the problem: longer distances and more complex energy landscapes require finer temporal discretization.
Monte Carlo Samples M: The number of samples M affects the quality of the action estimate. Even for high-dimenn dige-dimmnsional problems (e.g., 50D with M=1000), we maintain solution quality shown in the Appendix.
Computation of log(p) and V log(p): Following the scaling approach in the CNF literature [4), use the Hutchinson trace estimator to make the ODE for log(p) computationally tractable. To the best of knowledge, there are unbiased estimator techniques for V log(p), and consider developing such estimators an interesting direction for future work.
F Additional Numerical Results
In this subsection, we report all the experimental details. In we report the hyperparameters for the algorithms and the boundary conditions for each problem. The notation [x.y.]] Table defines the architecture of the networks, x is the input dimension, is the number of neurons per
layer, and is the number of layers. We assumed the value the constants ag = 01 in Algorithm which we report
in Table
Here we report that 20 is under CC BY-NC licence, has no license.
is under CC BY-NC 4.0 License, and
F.1 Pre-training strategy
F11.1 Zero initialized boundary parameters
Our framework can use pre-trained parameters 0g and 0, for the boundary conditions Po and P1 respectively. pre-train parameter, we use Flow Matching (FM) [T8] learn velocity field to sample from dataset. The choice of training framework for the boundary parameters motivates to adopt the FM loss function as our dissimilarity metric provides an upper bound on the W2 distance between the target density and the pushforward density. Therefore, minimizing the FM loss function guarantees the feasibility criteria. The simulation- free training scheme of FM offers a computationally efficcint method for pretraining boundary conditions and guaranteeing the feasibility of the boundaries throughout the density path optimization algorithm. In Figure 5a we show the pushforward of the spline obtained from initializing all the control points 0. In Figure 5b we show the resulting optimized path using the PDPO strategy. The running time for this case is 24 32s, leading to action value of 40.31. Since the collocation points are initialized to O, we did not use a geodesic-warmup Compare this solution with the one reported Section 4.2
demonstrated by [6], the FM loss function
F.1.2 Reusability of boundary parameters
crucial feature of PDPO its boundary parameter reusability: the boundary condition models Ag and 01 can be trained independently of action functional's I's interior dynamics. This independence enables flexible reuse of the boundary parameters.
Specifically, when solving multiple action minimization problems that share boundary constraints, PDPO offers the following flexibility: after training library of generative models with shared architecture (0 pm I .p....,..... }p2} for distributions {p1, .2.....PP. P, any action minimization problem

---

Initial density path
x10
Optimized density path
x10.
7.2
2
2
6.4
7
1
1
Potential Energ
0
4.8
0
4.0
-1
3.2
2.4
-2
2
-2
0
2
1.6
Dimension
Diieension
0.00
0.50
0.75
1.00
1
0.50
1.00
0.8
0.0
(a) Zero initialized boundaries
(b) Optimized path for zero initialized bds Figure No pertaining for Og and Oj.
Table 3: Comparison for the SCC problem. The quantities W2:p),, i = 0, 1, denote the empirical Wasserstein-2 distances at times = 0 and computed using the POT library [12].
Problem (method)
A(pel
Ep, [V()) dt (Wz(Po), W2(P1))
Time
SCC (PDPO)
432.23 440.02 + 0.75
3.42
(0.12, 0.11) (0.59, 1.76)
36 12s
SCC (GSBM)
3.12
53 m 13
with boundary conditions in this set can reuse the parameters of the pretrained models, regardless of the direction of transport the specific action functional. For example, if Problem transports from Pi to P. and Problem 2 transports from Pj Pi, both problems can use the same pretrained pair (0p, Ony) by simply swapping their roles as initial and terminal conditions.
This modularity significantly reduces computational cost when solving families of related problems. For instance, in the opinion depolarization example, once Og and 0 are trained (29m35s and 31m12s: respectively), they can be reused for alternative action functionals, different polarization dynamics, or even reversed boundary conditions, provided the marginal distributions remain unchanged.
We treat the pretraining phase for the boundary models as a one-time investment.
To demonstrate this reusability feature, we reuse the pretrained boundary-condttion models from VNEFI example in a different setting, under modified action functional and with reversed flow direction. In the VNEFI problem, the boundary conditions are
Po = -== ..57)
P1 P1===N {H))
Denote the corresponding pretrained models by 0₂ and 0, (training times: 0.38 and 0.39 s, respec- tively).
We then apply these same boundary models to solve the SCC problem from 20 reversing the flow direction so that Po = / and P1 = μ.. Figure ba shows the pushforward of the control points, Figure 6heth interpolated trajectory, and Figure 6c the corresponding GSBM solution. Quantitative results are provided in Table 3
As before, PDPO achieves lower action value and shorter training time than GSBM. It is important to emphasize that the boundary parameters are reused from the VNEFI example, so the cost of FM training is amortized across problems. Both methods yield similar constraint-violation levels, with GSBM showing a slightly smaller residual.
F22 Geodesic warmup
In the geodesic warmup, we optimize the linearly initialized control points using Algorithm T with F(p) = 0. This approach effectively reduces the computational cost of our algorithm. As we show in

---

Optimized (
Optimized censity path
K10
Optimized denstyy path
1.188
Dimension 2
Dimeraton
0.924
0.924
Dimension
0.792
0.660
0.660
0.645
0.528
0.528
0.516
6. 0.3968
0.396
0.3872
Pimension
0.204
Dimension
0.264
Dimension
0.250
0.00
0.25
0.50
0.75
1.00
0.132
0.00
0.25
0.50
0.75
1.00
0.132
0.00
0.25
0.50
0.75
1.00
0.129
0.000
0.000
0.000
(a) (Tor, A
(b) PDPO (Te)
(c) GSBM
Figure Comparison for SCC with new boundary conditions.
Figure 7 using geodesic warmup drastically improves the speed of convergence and allows PDPO reach lower local minimum in the s-curve example.
To reduce the computational cost of the geodesic warmup, we consider at most N 15 points for the trapezoidal rule. The solution computed without geodesic warmup took 5m 32ss, whereas the solution with geodesic warmup took 5m 57s. Thus, the computational overhead of the geodesic warmup is negligible.
In Figure 8 we show a comparison of solution without geodesic warmup (Figures Sa and 8b) and with geodesic warmup (Figures 8 and 8d).
Action Comparison
100
Witheat geodesic-wermug With grocesic warmpp
Action Value
(0-
-30
3.0
300
Iteration
Figure 7: Comparison of action for solutions with and without geodesic warmup.
iraua deedyy 2a11
Drwalnn
122
16
0.00
235
0.5
0.75
1.0
6.00
625
3.55
0.5
1.00
0.5
0.25
2.50
235
120
0.00
0.3
0.33
0.5
(a) Linear interpolation pa(h)
Solution geodesic warmup
without
(c) Geodesic warmup
(d) Solution with geodesic warmup
rameter space
Figure 8: Comparison solutions with and without geodesic warmup.
F.3
Ablation Control Points
In Figure we can see the comparison of the difference in the solution obtained by PDPO when varying the number of control points K. As we can when the number of control points increases, the quality of the solution also increases. Observe that the transition between the obstacles is smoother the number of control points increases.

---

Table 4: Experimental set-up
GMM
V-neck
S-tunnel
Opinion
d
2
2
1000
3
3
3
N
30
60
30
20
20
M
1000
1000
1000
1000
5000
Architecture
[2,256,4]
[2,128,4]
[2,64,4]
[2,128,4]
(18,0,128,4)
Epochs
15
15
18
10
10
Coupling opt steps
20
20
20
20
20
Path opt. steps
30
20
20
20
Geodesic warmup steps
100
100
100
200
200
100000
100000
100000
100000
10000
(No, Ki, K2)
(50,0.0)
(3000.50.0)
(100, .5) --2
(0,50,0)
(0,50,0)
Mean Po
e
1Giix
0.
0
0
Mean P1
8ém
.,...,3
0
0
4
0.5
0.25
Covariance of Po
0.51
0.11
diag(
0.25 )
diag(
:
0..55
Covariance of P1
0.51
0.011
31
31
((T.22AA135)) xie?
Optimized (1T)))))5))
x10²
x1o?
12
72
7.2
Dimerston2
Dimencion
4.8
48
4.0
4.0
312
32
32
2.42
2.42
-3
-1
-2
0
Amennson
Dimension
Time
0.0
Time 00
Time 03
0.0
(a) Optimized control points K (b) Optimized control points
(c) Optimized control points K
Optimized
x1e
Optimized denatty path,
x10²
Optimized deratty path
x10²
72
72
72
Dimrrsion
4
48
4.0
4.
4
32
32
325
2.4
2.4
Cimensicn
Cimension
Dimension
0.00
0.25
0.50
0.75
1.00
0.
0.00
0.25
0.50
0.75
1.00
08
0.00
0.25
0.50
0.75
1.00
08
0.0
0.0
000
(d) Optimized density path K
(e) Optimized density path
(0) Optimized density path
Figure 9: Comparison of solutions with increasing number of control points K.
F.4
Obstacle Avoidance with Mean Field Interactions
Here provide the definitions for the entropy E(p) and congestion potential C(p).
E(p) = [ log(p)pdz and C(p)= Izaxa || y11 p(
2

---

100
80
0.042
0.045
90
3.000
0.035
0.030
AVix.1
50
02.025
00
0.015
50
20
0.00
0.010
0.088
Epoch
10
Epoch
15
Epoch
10
Figure 10: Quantities of interest with uncertainty estimates along the training process, S-curve.
Remark: In Table 1 we reported the W2 distance of the boundaries for GSBM. The samples the terminal time = are generated by the forward solver, while the samples at the initial time t = produced by the backward solver.
S-curve The definition of the problem and source code were taken from 17.. We refer to it for definition of the obstacle. APAC-Net was more sensitive to the coefficient Ko and K22; there we used
Ri = 5 and K2 = 1.. The action reported Table] T for APAC-Net was obtained by evaluating its code with the values of Ko, K2 reported in In Figure [0]ee we show the W2 distance with the boundaries, action, kinetic energy, and potential function along the training process. In these plots, an epoch is the completion of coupling path optimization iterations. In these plots, we can clearly see that the feasibility condition is met from epoch 0 thanks to our pre-training strategy.
solution using
Both schedulers in this experiment are StepLR. The learning rate for the coupling optimizer is 10 with step size of 10, and y 0.9. The learning rate for the path optimizer 5 x 10-4, step size of 10 and 0.1.
V-Neck The definition and source code were taken from 20
When GSBM's training time is constrained to match PDPO's s training time, its action value is the same as before, but the boundary approximation noticeably worse, W2(po, Po) 0.12 + 0.002, W2(P1, Pi) 0.107 + 0.002 Ih 15m 24s.
report the quantities of interest along the training epochs. See 11 pushforward of the three control points, and the comparison of the PDPO and GSBM solutions.
In Figure 12
the
Both schedulers in this experiment are StepLR. The learning rate for the coupling optimizer is 10-1 with step size of 5, and 7 = 0.1 The learning rate for the path optimizer is 10 step size of and 7 = 0.25
Optimized (01))100)) K. 3xig!
Optinized density path
10
Optimized density path
x10
1.35
5.0
7.5
1,800
5.0
1.4
Dimendion
1.20
2.5
Dimension
1.2
0.90
1.350
0.0
1.05
0
0.75
1.325
-2.5
0 0.900
-2.5
6.8
Fetental
0.60
-5.0
0.675
-5.0
R.e
0.152
10
0.450
-10
Dimendion
10
Dimension
0.00
0.25
0.50
0.75
L.0
0.00
0.25
0.50
0.75
1.0
02
0.00
0.25
0.50
0.75
1.00
0.15
Time (0
Time (t)
lime
(a) Pushforward with control points
(b) PDPO
(c) GSBM
Figure 11: V-neck-E-FI a) Pushforward with control points. b) Pushforward with interpolated curve. Solution by GSBM.
GMM The definition and source code were taken from 20 See Figure 13 for the comparison the solutions. In this example, we use X Po- To guarantee this is satisfied, define dg == 0. the zero

---

0130
300
230
275
111 ( jar
0.20
250
10.15
-50
0.110
230
-100
0.105
0.12
175
0.11
150
-150
Spoch
bpoch
Epech
Epach
Spach
Figure 12: Quantities of interest with uncertainty estimates along the training process, vneck
Octiniied (151000
Octiriedd ecndry
Octiniedd eceoyy
Oetoraed Ecedyy
6.48
4.32
100
5Mg
-20
2.45
-10
2.33
2.0
2.42
144
1.44
1.1
0:97
100
6x00
5:54
6.7
1.0
6x0
6:00
0.10
0.00
0.20
009
(a) Pushforward with control points
(b) PDPO
(c) GSBM
(d) NLOT
Figure 13: Comparison for GMM example
vector. In Figure 14 we report the quantities of interest along the training epochs. See 13 to the pushforward of the five control points, and the comparison of the PDPO and GSBM solutions.
The schedulers in this experiment are cosine for the coupling optimization and StepLR for the path optimization. The learning rate for the coupling optimizer scheduler is To 5, Tauh 2, 7mmn 0.001, step size of 3, and 7 = 0.9
10-6, the setup for the cosine 10-6, The learning rate for the path optimization is
F,5 Opinion Depolarization
In this problem, opinions X(t) E R 1000 evolve according to polarizing dynamic:
dx(() dt
fpolarize Cacuriieex((f (I(t):Pe))
0.400
200
120
3.00062
92-
1.00
100
0.00063
0.385
1.00
8B-
40-
84-
0.00577
0.380
Epoch
15
13
Figure 14: Quantities of interest with uncertainty estimates along the training process, GMM

---

with
if sign((a, Ei) = sign( (v. Et)
Spolarize (ri Pt, Ex) = Byno, (a(a, y. 5. d(x,y,)))
otherwise
When this dynamic evolves without intervention, opinions naturally segregate into groups with diametrically opposed views. However, the desired outcome is unimodal distribution.
To solve this problem, we follow (20 and incorporate the polarizing dynamics as base drift or prior See | for reference in OT problems with prior velocity fields. Specifically, the action optimize in the parameter space is defined by
[ 1
d dt Taco (2)1 + F((Ta()) ) (A))
Apolatzo
( (f0-1) = BSNd
|I pootaiiie (2); (2); (Ta())
The potential energy term a congestion cost that encourages particles to maintain a distance from each other. We follow the experimental setup in [20] and [19] Because of the dimension, we can only simulate deterministic dynamics, as specified in 30), Figure [5 we compare the optimized density trajectory and directional similarity histogram. These plots show that both methods obtain nonpolarized trajectory.
50D Scrhödinger Bridge
The SB between Gaussians has closed-form solution g In Table 5 we evaluate PDPO's accuracy by computing two different path discrepancies: fo W (Pecn) ((t) and fu W (Pece), ((t)) Here, p(t) represents the theoretical SB solution between the original Gaussian boundaries, p(t) is the SB solution between the approximated boundaries (To, ) and (To, ) 2, and poct) is the density path recovered by PDPO. The first integral measures how closely PDPO approximates the true SB solution, while the second evaluates how well PDPO solves the counday--approxingg
problem.
While the tabular results quantify the global accuracy of our method, Figure [6] provides a more intuitive visualization through 2D projection of sample trajectories. the figure, blue points represent samples from our PDPO solution Pe(t)= green points show the theoretical solution p(t), and red points display the accurately approximates individual particle trajectories throughout the transport process, confirming the effectiveness of our approach even in high- nigh-dimensional settings.
matched solution p(t). This visualization demonstrates that PDPO
Table 5: SB between Gaussians
d
So W (Pet)) P(t))dt
fo W3 (Prc)) d(
50D
0.778 + 0.0018
0.789 + 0.001

---

Density Visualization and Directional Similarity Histogram 0.333
t
0.000
0.567
1.000
10
10
10
-5
-5
-5
-5
100
10
10
-10
10
-10
10
1000
1000
1000
1000
600
600
600
600
400
400
400
400
200
200
200
200
8.00
0.25
0.50
0.75
1.00
8.00
0.25
0.50
0.75
1.00
8.00
0.25
0.50
0.75
1.00
800
0.25
0.50
0.75
1.00
Hight Agres
Neatral
Hight Disogree
(a) PDPO
Density Visualization and Directional Similarity Histogram 0.667
0.000
0.333
t 1.000
10
10
10
10
-5
-5
-5
-5
-10
-10
-10
-10
-10
10
-00
10
10
1000
1000
1000
1000
600
600
500
400
400
400
400
200
200
200
200
8.00
1.00
8.00
1.00
8.0
8.00
0.75
Iighly Agree
Neuthal
Hight/ Uneqrre
(b) GSBM
Figure 15: Comparison of optimized opinion polarization dynamics.

---

0.00, MSE 0.0114
t= 0.24, MSE 0.0074
[= 0.41, MSE 0.0240
This parh approo bd The pors
Asprcorsciss pats
Tue poth meen
Yur pars appnns mein
Asprosmaze par moin
Figure 16: Particle comparison two random directions.


Reviews:
[{"review_id": "TYcAVc4Sjn", "reviewer": "Reviewer_CE31", "summary": "This paper introduces a method for the optimisation of paths followed by probability distributions under trajectory constraints. The proposed approach transforms the original Wasserstein-2 optimal transport problem into an equivalent dynamic reinterpretation and formulates a constrained least-action optimisation objective. This dynamic formulation more easily allows for the incorporation of additional trajectory constraints. Probability distributions and their paths and their paths are represented by parametric models. Cubic splines parameterise the paths, with endpoint densities given by flow matching models to approximate the target distributions and intermediate points approximated via neural ODEs. Theoretical results regarding the equivalence of the formulations and the approximation error are derived. Experiments on synthetic benchmarks evaluate the algorithms against the most relevant baselines in the area.", "questions": "* The final algorithms have quite a few components and moving parts (flow matching models, neural ODEs, cubic spline parameters, Monte Carlo estimators for entropy, Fisher Information, etc.). How are hyper-parameters tuned? The appendix presents an ablation study for the number of control points, but the paper would benefit from a discussion on how sensitive the method would be to the choice of models involved and the machinery involved in the numerical integration and optimisation steps.\n\n* With the spline-based framework, the method seems flexible enough to incorporate more than 2 fixed endpoints. Have the authors considered problems where intermediate points also need to satisfy hard constraints?\n\n* Can you revise and explain the definition of $\\tilde{\\Pi}$?", "limitations": "Limitations are quite briefly mentioned in the conclusion, without much detail or discussion.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\n\n* The paper tackles a hard problem with an interesting solution that connects different fields.\n\n* Theoretical guarantees are provided.\n\n* The discussion of experimental results seems to be solid, insightful and to go into enough depth for a conference paper.\n\n* Significant computational gains with competitive performance are demonstrated.\n\n\n### Weaknesses\n\n* The paper is quite mathematically heavy at times with notation that turns out a bit loosely defined and a few potential mistakes.\n\n    - In the introduction and subsequent restatements, the action functional equation (1) has a minimisation over two function-valued variables $\\rho$ and $v$, but they are indexed (under the $\\inf$) as scalars representing their point evaluations $\\rho(t, x)$ and $v(t, x)$, which is confusing.\n    - How can a joint probability density $\\pi$ be equated to a point tuple as $\\pi = (x, T(_\\pi(x))$? These are two objects of different nature. The same apparent mistake occurs first between lines 134 and 135 and it]s then repeated throughout the text. The set $\\tilde{\\Pi}$ they define then appears in Eq. 11, which is part of the result in Theorem 1, making the latter seem invalid.\n\n* \"Opinion\" dynamics is involved in the experiments, but the term is never explained or defined in the main text.\n\n* Quite a few of the references to the appendix are invalid in the main text, though fixed in the supplementary material.", "was_revised": true, "final_justification": "I've increased my review score given that the rebuttals have addressed most of my concerns and other relevant points I found raised in other reviews. There are a few details to clarify in the revised paper, which should be taken seriously, given the feedback from the reviews, and I'd like encourage the authors to do so.", "rebuttal": "You're absolutely right about the notation problems. Let us clarify:\n\n## **Equation (1) notation:**\n\n We'll change $\\rho(t,x), v(t,x)$ to simply $\\rho, v$ under the infimum to be mathematically correct.\n\n## **Definition of $\\tilde{\\Pi}$**\n\nWe agree that the notation $\\pi(x,y) = (x,T_{\\pi}(x))$ is wrong. The reviewer is correct, and the definition of $\\tilde{\\Pi}$ needs updating, although **Theorem 1 keeps its validity** as we demonstrate below.\n\nTo define the set of couplings properly, we first define the projections $p_x$ and $p_y$ as the projections of $\\mathbb{R}^d \\times \\mathbb{R}^d$ onto the first and second $\\mathbb{R}^d$ components, respectively. Given $z \\in \\mathbb{R}^d \\times \\mathbb{R}^d$ with $z = (x,y)$ where $x,y \\in \\mathbb{R}^d$, we have $p_x(z) = x$ and $p_y(z) = y$.\n\nGiven probability measures $\\mu_0, \\mu_1$, the set of couplings between $\\mu_0$ and $\\mu_1$ is:\n$$\\Pi(\\mu_0,\\mu_1) :=\\\\{\\pi \\in \\mathcal{P}(\\mathbb{R}^d \\times \\mathbb{R}^d): (p_x)\\\\#(\\pi) = \\mu_0, (p_y)\\\\#(\\pi) = \\mu_1\\\\}$$\n\nIn our paper, we focus on couplings of the form $\\pi = (id,T)\\\\#(\\mu_0)$, where $T: \\mathbb{R}^d \\to \\mathbb{R}^d$ is a measurable map such that $T\\\\#(\\mu_0) = \\mu_1$, and $id(\\cdot)$ is the identity map: $id(x) = x$.\n\nWe acknowledge that in the paper, we made an abuse of notation using measures and densities interchangeably. Here we provide the rigorous treatment: we use $\\mu$ for measures and $\\rho$ for probability densities.\n\n**Corrected definition of $\\tilde{\\Pi}$:**\n$$\\tilde{\\Pi} := \\{\\pi \\in \\Pi(\\mu_0,\\mu_1): \\pi = (id,T_{\\pi})\\\\#(\\mu_0) \\text{ for some injective measurable map } T_{\\pi}: \\mathbb{R}^d \\to \\mathbb{R}^d\\}$$\n\n**Why this is well-defined:** The hypothesis that $\\mu_0$ is absolutely continuous with respect to the Lebesgue measure allows us to use Brenier's theorem, which guarantees that the coupling $\\pi^{\\*}$ achieving optimal transportation cost lies on a graph. That is, $\\pi^{\\*} = (id,T^{\\*})\\\\#(\\mu_0)$ where $T^{\\*}$ is the Monge map. Since we also assume $\\mu_1$ is absolutely continuous, this map is bijective, ensuring $\\tilde{\\Pi} \\neq \\emptyset$.\n\n**Proof that Theorem 1 remains valid:**\n\nThe optimization problems:\n$$\\inf_{\\rho,v} \\int_0^1 \\int_{\\mathbb{R}^d} \\frac{1}{2}\\|v(t,x)\\|^2 \\rho(t,x) dx + F(\\rho(t)) dt$$\nsubject to $\\partial_t \\rho + \\nabla \\cdot (\\rho v) = 0, \\rho(0,\\cdot) = \\rho_0, \\rho(1,\\cdot) = \\rho_1$ (dynamic problem)\n\nand \n$$\\inf_{\\pi \\in \\tilde{\\Pi}(\\rho_0,\\rho_1)} c(\\pi) = \\inf_{\\pi \\in \\tilde{\\Pi}(\\rho_0,\\rho_1)} \\inf_{\\gamma \\in \\Gamma(\\pi)} \\mathcal{A}(\\gamma)$$\n(static problem) are equivalent.\n\n**Direction 1 (Dynamic → Static):** A local minimizer $(\\rho,v)$ of the dynamic problem defines the $C^1$-diffeomorphic curve:\n$$\\gamma_{\\rho,v}(t,x) = x + \\int_0^t v(s,\\gamma(s,x)) ds$$\nand coupling $\\pi_{\\rho,v} = (id,\\gamma(1,\\cdot))\\\\#(\\mu_0)$. \n\nFrom the definition of $\\gamma_{\\rho,v}$, the distribution $\\mu_{\\gamma_{\\rho,v}(t)} = (\\gamma_{\\rho,v}(t,\\cdot))\\\\#(\\mu_0)$ has density given by the change of variables formula.\nand $\\rho_{\\gamma_{\\rho,v}(t)}(x) = \\rho(t,x)$, thus:\n$$\\int_0^1 \\int_{\\mathbb{R}^d} \\frac{1}{2}\\|v(t,x)\\|^2 \\rho(t,x) dx + F(\\rho(t)) dt = \\mathcal{A}(\\gamma_{\\rho,v})$$\n\nTherefore, the dynamic problem upper bounds the static problem.\n\n**Direction 2 (Static → Dynamic):** Given a local minimizer $(\\pi,\\gamma)$ of the static problem, let $\\mu_{\\pi,\\gamma} = (\\gamma_t(\\cdot))\\\\#(\\mu_0)$ with corresponding density $\\rho_{\\pi,\\gamma}(t,\\cdot)$. Define the velocity field:\n$$v_{\\pi,\\gamma}(t,\\gamma(t,x)) = \\frac{d}{dt}\\gamma(t,x)$$\n\nThe pair $(\\rho_{\\pi,\\gamma}, v_{\\pi,\\gamma})$ satisfies the conditions of the dynamic problem. As before:\n$$\\mathcal{A}(\\gamma) = \\int_0^1 \\int_{\\mathbb{R}^d} \\frac{1}{2}\\|v_{\\pi,\\gamma}(t,x)\\|^2 \\rho_{\\pi,\\gamma}(t,x) dx + F(\\rho_{\\pi,\\gamma}(t)) dt$$\n\nThus, the static problem upper bounds the dynamic problem.\n\nSince both problems upper-bound each other, they are equivalent\n\n## **Extensions and Future Work**\n\n**Intermediate constraints:** We have thought about it, and plan to work on a PDPO extension for these type of problems. In that direction, we preview that a modification on how the coupling works in Theorem 1 will have to be made. This will also have an effect on how the boundary parameters (and fixed intermediate points) need to be updated. \n\nOther feasible directions for extending our work are:\n- \"Overdamped Brownian Motion in a Force Field\" see section 8 in [1]\n- \"Density Control With Multiple Species\", see section 4 [2] \n- Optimal transport with nonlinear diffusion equations, see [3]\n\n\n## **Hyperparameter Sensitivity and Tuning**\n\n**Flow Matching Architecture:**\nThe FM model defines the parameter space and must be sufficiently expressive to capture the boundary distributions. Key considerations:\n- **Architecture sizing**: Use larger networks than minimally required for boundary fitting. For example, while [2,64,4] can learn our V-neck boundaries accurately, we needed [2,128,4] to capture the solution complexity when potential terms significantly alter the density shape.\n- **Rule of thumb**: If the potential energy terms are expected to create new modes or dramatically change boundary shapes, increase the hidden dimension.\n\n\n**Neural ODE Integration:**\nTwo factors affect *forward quality:\n- **Solver choice**: We use midpoint integration, as it is standar in flow based generative models.\n- **Integration steps**: 10 steps suffice for most problems. Use the maximum of the number of steps that your boundary NODEs require for accurate sampling from target datasets.\n\n**Action Approximation (Time Steps N):**\nThis an important parameter. \n\n| N | Action | Kinetic Energy | Potential Energy |\n|---|--------|----------------|------------------|\n| 10 | 66.06 | 37.38 | 28.68 |\n| 20 | 42.09 | 37.77 | 12.82 |\n| 30 | 30.31 | 29.46 | 0.788 |\n| 40 | 30.20 | 29.72 | 0.735 |\n| 50 | 30.10 | 29.37 | 0.725 |\n\n**Recommendation**: Start with N=20 minimum, increase until action converges (typically N=30-50). The dramatic drop from N=10 to N=20 shows why sufficient discretization is crucial. This parameter is related to the complexity of the problem, the longer the distance and the more complex the landascape,\n\n\n**Monte Carlo samples M**: If the number of samples affects the quality of the estimate. Although even for the 50D case with 100,0 we were able to maintain the quality of the solution as shown in the Appendix\n\n**Computation of $\\log(\\rho)$ and $\\nabla\\log(\\rho)$** The scaling issue of the computation of $\\log(\\rho)$ in the CNF paper [4], the Hutchinson trace estimator is used to make the ODE for $\\log(\\rho)$ computaionally tractable. To the best of our knowledge, there are no estamator techinques for $\\nabla\\log(\\rho)$ and we consider this to be an interesting problem for future development. \n\n\n## **Reference Formatting**\n\nWe'll update citations to use proper author names: \"as shown by Chen et al. [10]\" instead of \"as in [10]\".\n\n## References\n[1] On the Relation Between Optimal Transport and Schrödinger Bridges: A Stochastic Control Viewpoint Yongxin Chen, et all.\n[2]Density Control of Interacting Agent Systems, Yongxin Chen, IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 69, NO. 1, JANUARY 2024.\n\n[3] Nonlinear mobility continuity equations and generalized  displacement convexity. J. A. Carrillo, S. Lisini, G. Savare , D. Slepcev. Journal of Functional Analysis Volume 258, Issue 4, 15 February 2010, Pages 1273-1309.\n\n[4] FFJORD: FREE-FORM CONTINUOUS DYNAMICS FOR SCALABLE REVERSIBLE GENERATIVE MODELS. Will Grathwohl et al.  ICLR 2019"}, {"review_id": "Tawn3aJwlV", "reviewer": "Reviewer_5bmt", "summary": "This paper introduces Parametric Density Path Optimization (PDPO), a novel method for calculating the optimal, action-minimizing paths between different probability densities. The core innovation is to reframe this complex, infinite-dimensional problem into a more manageable, finite-dimensional one. It achieves this by representing the probability path as the transformation (or \"pushforward\") of a simple, fixed reference density through a parametric map, which is modeled by a Neural ODE. The path is then optimized by interpolating a small number of control points (typically 3-5) in the parameter space using cubic Hermite splines, which defines the continuous evolution of the density.\n\nPDPO demonstrates superior performance compared to existing state-of-the-art methods across several challenging benchmark tasks, including obstacle avoidance, entropy-constrained transport, and high-dimensional opinion dynamics. The method is not only more computationally efficient, achieving 40-80% faster runtimes, but it also produces higher-quality solutions, with 7-39% lower action and up to 10 times better boundary accuracy. A key strength of PDPO is its flexibility, as it can be applied to a wide range of problems involving external potentials, mean-field interactions, and higher-order dynamics, such as acceleration minimization, within a single, unified framework.", "questions": "### **1. On the Source of Computational Efficiency**\n\nRegarding the computational efficiency of PDPO compared to GSBM, could you clarify the contribution of the **Flow Matching (FM) initialization**? The paper notes that PDPO achieves significantly faster runtimes. Is it possible that this speed-up is largely attributable to the FM initialization step? To help disentangle the gains from the core PDPO framework versus the initialization strategy, it would be insightful to see a comparison where the GSBM baseline is also initialized using a similar FM-based approach.\n\n***\n\n### **2. Clarification on Constraint Violation in Figure 2**\n\nIn the obstacle avoidance experiment shown in Figure 2, the paper claims that both GSBM (2c) and APAC-Net (2d) violate the obstacle constraints, while PDPO effectively avoids them. However, on visual inspection, the optimized PDPO path (2b) also appears to have some samples crossing into the obstacle regions. Could you please provide a **quantitative measure of the constraint violation** (e.g., the percentage of the mass inside the obstacle) for all three methods? This would offer a more objective comparison of their performance on this task.\n\n***\n\n### **3. Request for Concrete Implementation Details**\n\nTo improve the paper's reproducibility and provide a clearer understanding of the methodology, could you offer more concrete implementation details or a working example for two key aspects of the framework?\n\n* **Parameterization of the Pushforward Map**: While the paper states that the map $T_{\\theta}$ is a Neural ODE, a more detailed description of the network architecture ($v_{\\theta}$) and how its parameters $\\theta$ are structured would be very helpful.\n* **Flow Matching Initialization**: A step-by-step walkthrough of how Flow Matching is employed to initialize the boundary parameters $\\theta_0$ and $\\theta_1$ would greatly clarify the practical application of this crucial first step.", "limitations": "### Scalability in High-Dimensional Stochastic Optimal Control (SOC)\n\nA primary limitation is the method's computational cost when dealing with problems that involve Fisher Information, which is used to solve Stochastic Optimal Control (SOC) problems.\n\n* The cost to compute the action functional increases substantially with the problem's dimension, `d`.\n* Calculating the Fisher Information requires solving a separate `d`-dimensional Ordinary Differential Equation (ODE).\n* The authors state that a 50-dimensional Schrödinger Bridge was the \"highest-dimensional case with Fisher Information we have successfully computed,\" highlighting a practical cap on its scalability for this class of problems.\n\n***\n\n### Architectural Restrictions\n\nThe current implementation of PDPO has architectural constraints that limit its range of application.\n\n* The framework uses Neural ODEs (NODEs) built with Multilayer Perceptron (MLP) architectures.\n* This reliance on MLPs restricts the method's direct use for certain data types, particularly images, where specialized architectures like U-Nets are more effective.\n* This is a notable limitation compared to methods like GSBM, which can be applied to image-based transport problems.\n\n***\n\n### Methodological and Evaluative Limitations\n\nThe methodology and evaluation presented in the paper have some implicit limitations.\n\n* **Dependence on Initialization:** The three-step optimization process relies on a strong initial approximation of the boundary distributions using Flow Matching (FM). While effective, this also means the final solution's quality is dependent on this prerequisite step.\n* **Potential for Unfair Comparison:** The performance gains reported over baselines like GSBM could be partially attributed to the specialized FM initialization. The paper does not provide an ablation study comparing performance when baselines are also given the benefit of a similar initialization, which makes it difficult to isolate the gains from the core PDPO algorithm alone.\n* **Convergence to Local Minima:** The optimization algorithm, being gradient-based, is guaranteed to converge to a local minimum, not necessarily a global one. The final solution is therefore sensitive to the initial path, again underscoring the importance of the initialization procedure.\n* **Qualitative Evaluation:** In some experiments, such as obstacle avoidance, the claims of superiority are based on visual inspection of plots. The paper does not provide quantitative metrics for constraint violation (e.g., the amount of probability mass inside an obstacle), relying on a qualitative assessment that can be subjective.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\n\n* **Novel and Theoretically-Grounded Formulation:** The paper introduces Parametric Density Path Optimization (PDPO), a novel method that reframes a complex infinite-dimensional optimization problem into a more tractable finite-dimensional one. The core strength lies in its static formulation of the action-minimizing problem, which the authors claim is novel for handling internal and interaction energy terms, unlike previous work. This new formulation is supported by theoretical guarantees, including a theorem proving its equivalence to the classical dynamic formulation (Theorem 1) and an error bound for the cubic spline approximation (Theorem 2).\n\n* **Superior Performance and Efficiency:** The paper provides strong empirical evidence that PDPO outperforms state-of-the-art methods in benchmark tasks.\n    * **Accuracy:** It achieves 7-39% lower action costs and up to 10x better boundary accuracy compared to baselines like GSBM and NLOT.\n    * **Speed:** It demonstrates significantly faster runtimes, ranging from 40-80% on some tasks to over 10x faster in high-dimensional cases (e.g., nearly 2 hours for GSBM vs. under 10 minutes for PDPO in the GMM experiment, and over 5 hours vs. 23 minutes for the opinion dynamics problem).\n    * **Simplicity:** A key advantage is its ability to resolve complex, high-dimensional problems accurately using only 3-5 control points for the spline interpolation.\n\n* **Flexibility and Broad Applicability:** PDPO is designed as a general framework that can flexibly accommodate a wide variety of problems within a unified optimization scheme. It is shown to be effective for tasks involving:\n    * External potentials and obstacles.\n    * Stochastic Optimal Control (SOC) through a Fisher Information term.\n    * Mean-field interactions and congestion.\n    * Higher-order dynamics, such as acceleration-regularized problems, which the authors claim PDPO is the first method capable of solving in this context.\n\n### Weaknesses\n\n* **Potentially Unfair Comparison Due to Initialization:** The proposed optimization strategy begins with an efficient initialization of the boundary parameters using Flow Matching (FM). While effective, this raises a question about the fairness of the performance comparison. The significant speed-up over GSBM might be partially due to this specific initialization strategy rather than solely the core PDPO framework. The paper does not provide an ablation study or results for a GSBM baseline enhanced with a similar FM initialization, which would be necessary to isolate the true source of the efficiency gains.\n\n* **Reliance on Qualitative Visual Evidence for Key Claims:** In the obstacle avoidance experiment (Figure 2), the paper claims that competing methods violate the constraints while PDPO does not. However, visual inspection of the PDPO result (Figure 2b) suggests that some samples may still lie within the obstacle. The paper does not provide a quantitative measure of constraint violation (e.g., the percentage of probability mass inside the obstacle) to substantiate its claim of superiority, relying instead on visual interpretation which can be ambiguous.\n\n* **Acknowledged Scalability and Architectural Limitations:** The authors are transparent about certain limitations of the current method.\n    * **Scalability in SOC:** The computational cost for problems involving entropy and Fisher Information increases substantially with the problem's dimension. The authors note that the highest-dimensional SOC problem they successfully computed was 50-dimensional, implying that scalability in this domain remains a challenge.\n    * **Model Architecture:** The current implementation uses Neural ODEs with MLP architectures. This restricts its direct application to other data modalities, such as images, where methods like GSBM can use more suitable architectures (e.g., U-Nets). The authors explicitly state this as a limitation for image-based transport problems.", "was_revised": true, "final_justification": "I understand the effort the authors made, but I believe this does not affect my final rating, and I decided to leave the score unchanged.", "rebuttal": "Thanks for the detailed review and recognition of PDPO's novelty and flexibility!\n\n## **Q1: On the Source of Computational Efficiency**\n\nGreat question! Let us clarify the contribution of Flow Matching (FM) initialization and why our speedups are genuine.\n\n**PDPO algorithm does not rely on FM**, although we strategically selected it. Let's review what exactly is needed to setup PDPO, why FM was selected, and how the initialization works.\n\nPDPO solves the nested minimization problem:\n$$\\inf_{(\\theta_0,\\theta_1)\\in \\Theta_{0}^1}\\inf_{\\theta_{0\\to1} \\in \\Theta_{0\\to1}}\\mathcal{A}(\\theta_{0\\to1})$$\n\nThe outer optimization requires parameters in the set:\n$$\\Theta_0^1 =  \\\\{(\\theta_0,\\theta_1)\\in \\mathbb{R}^D\\times \\mathbb{R}^D: (T_{\\theta_i})\\\\#,\\lambda = \\rho_i, i = 0,1\\\\}$$\n\nThis is essentially a **generative modeling question**: given reference density $\\lambda$, target dataset $\\rho_0$, and mapping $T_{\\theta}$, how do we train $T_{\\theta}$ to generate samples of $\\rho_0$? \n\nWe don't need FM specifically - we need **any generative model** that:\n1. Uses invertible mappings (NODEs) for computing $\\log(\\rho)$ and $\\nabla_x\\log(\\rho)$\n2. Has an efficient loss function to evaluate $(T_{\\theta_i})\\\\#\\lambda \\approx \\rho_i$\n\nExamples include Flow Matching, Stochastic Interpolants, Diffusion, CNFs, etc. We chose FM because it has a **simple and fast loss function** to evaluate during boundary optimization steps.\n\n**Key insight:** Once trained, the parameter $\\theta_0$ representing density $\\rho_0$ can be **reused across different action minimization problems** with different functionals $F$ and different terminal conditions $\\rho_1$ **without retraining**.\n\n**Intialization strategy:** The reviewr asks if its possible to a initialize GSBM with a similar FM-based approach. The answer is yes. \n\nGSBM uses two NODEs, a forward and bacwkard that define vector fields that transport $\\rho_0\\to\\rho_1$ and $\\rho_1\\to\\rho_0,$ respectively. To use a FM like initalization. We propose to initialize the forward and backward velocity fields as pre-trained NODEs. We use the GSBM method/algorithm to train them with a zero potential function and $\\sigma = 0.1.$ This reduces to train using Stochastic Interpolants. The initialization strategy was tested on the S-curve problem but there were no noticeable effects. We pretrained the network for approximately 40s each and verified that they sampled accurately the target density. The the training time and action value remain almost the same as before. \n\nTo clarify the computational cost of the FM initialization, we report the FM training details. In the table GMM1 refers to $\\rho_1$ in the GMM example. In that case $\\lambda = \\rho_0$ so no training is needed\n\n| $\\rho$       | T      |\n| -------------------------------------------------------------------------------------------------- | ------ |\n| $\\mathcal{N}(\\begin{bmatrix} -2 \\\\ -2 \\end{bmatrix},0.1I)$         | 36s    |\n| $\\mathcal{N}(\\begin{bmatrix} 2 \\\\ 2 \\end{bmatrix},0.01I)$          | 38s    |\n| $\\mathcal{N}(\\begin{bmatrix} -11 \\\\ -1 \\end{bmatrix},0.5I)$        | 41     |\n| $\\mathcal{N}(\\begin{bmatrix} 11 \\\\ 1 \\end{bmatrix},0.5I)$      | 36     |\n| GMM1   | 3m10s  |\n| $\\mathcal{N}(\\vec{0}_2,\\text{diag}\\begin{bmatrix} 0.5 \\\\ 0.25 \\end{bmatrix})$     | 41s    |\n| $\\mathcal{N}(\\vec{0}_2,3I)$    | 41s    |\n| $\\mathcal{N}(\\vec{0}_{1000},\\text{diag}\\begin{bmatrix} 4 \\\\ 0.25 \\\\ \\vdots \\\\ 0.25 \\end{bmatrix})$ | 29m35s |\n| $\\mathcal{N}(\\vec{0}_{1000},3I)$     | 31m12s |\n\n By examining the number of parameters, we noticed that PDPO has fewer parameters in most cases. To test whether parameter count explains PDPO's speed advantage, we solve S-Curve ($\\sigma = 0$) and V-Neck ($\\sigma = 1$) examples using identical MLP architectures [2+1,128,4] for both methods. Both problems use $\\rho_0 = \\mathcal{N}(\\begin{bmatrix} -11 \\\\ -11 \\end{bmatrix},0.5I)$ and $\\rho_1 = \\mathcal{N}(\\begin{bmatrix} 11 \\\\ 11 \\end{bmatrix},0.5I)$ with the same pretrained boundary parameters $\\theta_0, \\theta_1$ (training times: 0.33s and 0.35s respectively).\n\n| Problem | Method | Time | $\\mathcal{A}$ | $W_2^2(t=0)$ | $W_2^2(t=1)$ |\n\n|---------|--------|------|-------------|--------------|--------------|\n\n| S-Curve | PDPO | 12m34s | 471.86 | 0.099 | 0.115 |\n\n| S-Curve | GSBM | 25m35s | 485.69 | 1.816 | 1.835 |\n\n| V-Neck | PDPO | 1h10m | 148.91 | 0.085 | 0.081 |\n\n| V-Neck | GSBM | 1h16m | 153.80 | 0.192 | 0.283 |\n\nWith identical architectures, PDPO still outperforms GSBM in both action value and boundary approximation quality, demonstrating that performance differences are not due to parameter count alone. This also illustrates boundary parameter reusability across different problems.\n\n\n## **Q2: Clarification on Constraint Violation in Figure 2**\n\nYou're absolutely right - we need quantitative measures! Here's the breakdown of $\\mathbb{E}_{X\\sim\\rho(t)}[V(X)]$ where $V$ is **only the obstacle potential**:\n| Method   |t = 0.1 |t = 0.2  | t = 0.3  | t = 0.4  |t = 0.5  |t = 0.6 | t = 0.7  | t = 0.8  | t = 0.9 | $\\int_0^1\\mathbb{E}_{X\\sim\\rho(t)}[V(X)]dt$  |\n| -------- | ------- | ------- | ------- | ------- | ------- | ------- | ------- | ------- | ------- | ----- |\n| PDPO     | 0.03 |0.72 | 3.21 | 0.0 | 0.0     | 0.0     | 4.31     | 0.0     | 0.0     | 0.99   |\n| GSBM     | 3.56 | 6.37 | 2.37 | 0.34$ | 0.4$ | 3.46 | 0.14      |    3.76     |    1.24     |  2.57     |\n| APAC-Net |   0.77  | 4.75    |    0.52     |    0.02     |    0.63     |    0.91     |     7.27    |    0.82     |    0.46     |   1.48    |\n\nThis **quantitatively** shows that PDPO finds the path with less constraint violation.\n\n## **Q3: Request for Concrete Implementation Details**\n\n**Parametrization of the Pushforward Map:**\n\nYou've identified a fundamental paradigm mismatch that's worth clarifying. PDPO treats neural network parameters as function arguments in a functional programming style, while PyTorch follows an object-oriented paradigm where parameters are object attributes.\n\nThe Core Challenge: When we interpolate spline control points $\\theta(\\tau_j)$ and need to evaluate the pushforward $(T_{\\theta(\\tau_j})\\\\#\\lambda$, PyTorch's '''torch.load()''' operation treats $\\theta(\\tau_j)$ as static model states rather than dynamic function arguments. This breaks the computational graph connecting $\\theta(\\tau_j)$ back to the optimizable control points.\n\nCurrent Solution: We've implemented parametric MLPs that accept parameters as explicit function arguments:\n\n```python\ndef forward(self, x, theta):\n    return parametric_mlp_forward(x, theta)\n```\nWe acknowledge this currently restricts us to architectures we can reimplement from scratch. This is a genuine limitation for practitioners wanting to use existing complex architectures like U-Nets for image-based transport. \n\nWhy JAX/Flax Resolves This: JAX's functional programming paradigm naturally aligns with our approach. In Flax, parameters are function arguments by design. This eliminates the computational graph issues entirely and allows direct use of any Flax architecture without implementation. We're actively developing the JAX/Flax implementation and expect to demonstrate arbitrary architecture support (including U-Net equivalents) in follow-up work.\n\n\n**Architecture details:**\n\n| Problem | Architecture | Meaning |\n|---------|-------------|---------|\n| V-neck | [2,128,4] | 2D input, 128 hidden units, 4 layers |\n| GMM | [2,256,4] | 2D input, 256 hidden units, 4 layers |\n| Opinion | [1000,128,4] | 1000D input, 128 hidden units, 4 layers |\n\n**Flow Matching Initialization step-by-step:**\n\nFirst, select the architecture for the NODE. For example, an MLP with 4 layers each with 128 neurons.\n\nThen, for each boundary condition:\n1. Generate training dataset from target distribution\n2. Train NODE using FM loss with a random coupling.\n3. Store state dictionary: `torch.save(model.state_dict())`\n\nWe recommend training boundary conditions $\\theta_0$ and $\\theta_1$ **in parallel** for efficiency.\n\nTo initialize the PDPO algorithm:\n1. Load trained networks: `torch.load('path_to_state_dictionary')`\n2. Convert to tensor using the `state_dictionary_to_tensor` function for each network, obtaining the tensors $\\theta_0$ and $\\theta_1$\n3. Initialize control points: $\\theta_{t_i} = (1-i/K)\\theta_0 + (i/(K+1)\\theta_1.$\n4. Using the torch tensors $\\{\\theta_{t_i}\\}_{i = 0}^{K+1}$ use Algorithm 1 using parametric_mlps as described before. \n\n\n## References\n[1] GSBM"}, {"review_id": "KOoymrJ4HK", "reviewer": "Reviewer_Nr1t", "summary": "The proposed method, called Parametric Density Path Optimization (PDPO), utilizes an action minimization framework incorporating additional constraints such as external potential energy and interaction and internal energy to find an action-minimizing path between probability densities. This generalizes the classical approach used in the case of Wasserstein distances which can also be seen via the action-minimizing viewpoint.\n\nThe authors work with a parameterization which allows for a finite dimensional representation making the action minimization problem computationally tractable. The authors furthermore present a theorem establishing the equivalence between a  static and dynamic representation viewpoint of their constrained action-minimization approach (Theorem 1).\n\nA cubic Hermite spline approximation of parameter space curves together with a neural ODE approach is then used to solve the optimization problem practically and an error estimate for the spline-based approximation is provided (Theorem 2).\n\nExperiments are presented that include two-dimensional settings of obstacle avoidance problems, showing improved action cost and significantly reduced computation time vs baseline methods.\nThe authors also extend their approach to incorporate an acceleration term and show 2d examples where PDPO can find solutions also in this setting. Additionally, a 1000-dimensional case of opinion depolarization problem is considered where the method again has improved computational time as compared to the baseline GSBM approach.", "questions": "- The current 2D maps you considered still have a relatively simple spatial structure. Do you think your method will generalize to complex maps directly, or do you see challenges in convergence as problem complexity is increased?", "limitations": "- It would be interesting to add some discussions of potential failure cases of the proposed method as part of the limitations. An experiment where the approach could not be successfully applied yet or where other challenges of interest arise may be helpful to guide readers towards open problems within this area.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- overall, the paper is well-written and clearly describes the chosen problem and solution approach\n\n- significant improvements in computational complexity and improvements or matching performance in minimization of the objective function over the baseline approaches appear empirically supported by the presented experiments.\n\n- the work tackles a rather general action minimization path finding problem between probability densities of relevance to several application domains in the NeurIPS community.\n\n- the paper presents both the authors' theoretical insights and formulation (e.g. Theorem 1 and 2) as well as a practical algorithm and implementation that combines a neural ODE formulation in combination with the proposed parameterize to arrive at a practically feasible solution.\n\nWeaknesses:\n- The chosen obstacle avoidance scenarios are interesting, but may at present still be quite far from the required real-world complexity to be useful in end-applications such as in robotics, where much larger maps with significant detail would likely have to be considered to be of practical relevance. \n- To assess the practical significance of the proposed approach it may be helpful to showcase an experiment with much larger problem complexity as well.", "was_revised": true, "final_justification": "The work makes a valuable contribution in my view. Following the rebuttal comments by the authors, I maintain a \"borderline accept\" recommendation primarily due to the lack of an additional more realistic/higher complexity experimental evaluation scenario.", "rebuttal": "We thank the reviewer for their positive assessment of PDPO's theoretical contributions, computational efficiency, and practical relevance. Below, we address their constructive questions:\n\n# Q1, challenges for complex maps.\nFor more complex environments, we recommend tuning the number of control points K and the number of discretization steps for computing the action, N. We usually try to keep these parameters as low as possible, as they increase the computational cost of our method. \n\n#Open problems and limitations\n\n- Intermediate constraints or Multi-Marginal Constraints:  Extending PDPO to enforce constraints at intermediate timepoints (beyond the two-endpoint setting) introduces new coupling structures in the variational formulation. We are currently exploring how to modify Theorem 1 accordingly, which would also affect how intermediate control points are updated.\n- \"Overdamped Brownian Motion in a Force Field\" Section 8 in [Chen, 2024] studies an energy-minimizing formulation of overdamped Langevin dynamics in force fields.\n- \"Density Control With Multiple Species\", see section 4 [2]. Coupled transport with cross-species interactions.\n- Optimal transport with nonlinear diffusion equations. We want to make PDPO  compatible with nonlinear mobility equations such as those arising in nonlinear Fokker–Planck-type models or porous medium equations, e.g., in [3]\n- The Fisher Information term appears in action functionals corresponding to entropy-regularized or stochastic control problems. In our current implementation, the main bottleneck lies in evaluating the right-hand side of Equation (20), which involves the divergence of the score function. In contrast to $\\log(\\rho)$, efficient Hutchinson-like estimators for $\\nabla(\\nabla\\cdot v)$  are less developed. Exploring scalable approximations to this term—especially using trace estimators or score-based models—would help push the dimensional limits of PDPO for stochastic optimal control.\n\n\n[1]Density Control of Interacting Agent Systems, Yongxin Chen, IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 69, NO. 1, JANUARY 2024.\n\n[2] Nonlinear mobility continuity equations and generalized  displacement convexity. J. A. Carrillo, S. Lisini, G. Savare , D. Slepcev. Journal of Functional Analysis Volume 258, Issue 4, 15 February 2010, Pages 1273-1309.\n\n[3] FFJORD: FREE-FORM CONTINUOUS DYNAMICS FOR SCALABLE REVERSIBLE GENERATIVE MODELS. Will Grathwohl et al.  ICLR 2019"}, {"review_id": "4W9kPnZigV", "reviewer": "Reviewer_Q47s", "summary": "This paper introduces Parametric Density Path Optimization (PDPO), a framework that casts the infinite‑dimensional least‑action problem over density curves into a finite‑dimensional spline optimization in the parameter space of a push‑forward map. The authors select neural ODEs for the map, discretize the parameter curve with cubic Hermite splines, and prove that using only 3 – 5 control points suffices to keep the action error $O(h^{\\kappa-1})$.  Optimization proceeds by alternating a “path step’’ on interior control points with a “boundary step’’ on the endpoints (Algorithms 1‑3).  Experiments on obstacle‑avoidance transport, generalized momentum minimization, acceleration‑regularized dynamics, and a 1000‑d opinion‑depolarisation benchmark show 7–39 % lower action and up to 10× faster wall‑time than GSBM, NLOT and APAC‑Net in the authors’ setup.", "questions": "- Is it possible to  provide theoretical guarantees for the convergence of the alternating optimization procedure? Under what conditions does the bilevel optimization converge to a global or local optimum? The current analysis is insufficient for a rigorous optimization method.\n- How sensitive is the method to the initialization of boundary parameters? Can you provide a systematic analysis of when the Flow Matching initialization will succeed, and what happens when it fails to find good boundary parameters?\n- Computational Complexity: The $O(N \\times M \\times (1 + d))$ scaling for Fisher information terms seems prohibitive for truly high-dimensional problems. Can you provide a more detailed computational analysis and comparison with baselines using the same computational resources and implementation quality?\n- Can you provide more rigorous statistical analysis including confidence intervals, significance tests, and fair computational comparisons? The current evaluation with only three seeds and potentially unfair runtime comparisons undermines the claimed performance improvements. Please also answer the point raised on taking initialization into runtime as I raised above.\n- What are the precise assumptions required for the two theorems , and how realistic are they in practice? When do the error bounds become meaningful, and how do they relate to the quality of the final solution?", "limitations": "Yes but not enough.", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "- The paper addresses an important problem in optimal transport and provides a computationally efficient alternative to existing methods that struggle with high-dimensional settings.\n- I think the paper’s strongest point is the clean static‑coupling reformulation (Eq. 11) and its parametric version (Eq. 13) that lets one exploit off‑the‑shelf autodiff once everything lives in the parameter space. The spline error bound is convenient and, to my knowledge, new for neural push‑forward families.  \n- Empirically PDPO is much faster than GSBM on the same GPU while matching or improving action and endpoint fidelity, especially in the 1 000‑d test where runtime drops from five hours to 23 minutes.\n\nQuality‑wise, however, several issues remain:\n \n- First, the dynamic/static equivalence (Theorem 1) merely states without the assumptions on the regularity of $\\pi$. Moreover, the non‑emptiness of $\\Gamma(\\pi)$ looks restrictive but is not discussed in the main text, so I am not sure how often the theory covers the practical neural‑ODE setting. There are also no convergence guarantees for the proposed bilevel optimization procedure, which alternates between updating interior control points and boundary parameters, which I think is critical.\n\n- Second, the experimental evaluation, while showing improvements, has notable limitations. Many experiments only compare against GSBM, and the baseline comparisons may not be entirely fair due to different implementations and computational frameworks. The statistical analysis is minimal - using only three random seeds.\n\n- Third, from a technical standpoint, the method requires careful initialization of boundary parameters, and the authors acknowledge this is done \"outside of the main algorithm.\" This suggests the approach may be sensitive to initialization and require problem-specific tuning. The scalability analysis reveals concerning computational complexity of $O(N \\times M \\times (1 + d))$ for problems involving Fisher information terms, which severely limits applicability to high-dimensional problems. \n\n- Fourth, I find this statement in Section 3.3 troublesome: **\"We in general consider this step outside of the main algorithm. A parameter can be used for any functional, so once a boundary parameter has been obtained, it can be reused across multiple problems.\"** I think this makes the runtime comparison completely unfair.\n- Finally, the paper also lacks important details about convergence criteria and stopping conditions for the alternating optimization procedure. The claim that \"as long as the action is lower bounded, step 2 will converge to a local minimum\" is insufficiently rigorous without proper analysis of the optimization landscape.\n\n- Clarity-wise, the authors have done a decent job. Just one point that I am not a big fan of putting the algorithm boxes in the Appendix, especially with the split of supplementary material for NeuRIPS submission.", "was_revised": true, "final_justification": "I think that the rebuttals (include ones to other reviewers) reveal substantial gaps in the original submission. The missing theoretical machinery, and implementation details about parametric MLPs represent fundamental content needed for understanding your work. Furthermore, while the 1000D opinion problem provides some evidence of scalability, the method remains limited to 2D toy problems and 50D Gaussians, with acknowledged architectural restrictions to manually-implemented MLPs. Notably, GSBM demonstrates image transport capabilities (AFHQ dataset, Section 4.2) that PDPO cannot handle due to your architectural limitations. This selective comparison on problems favorable to the method's constraints undermines claims of general applicability.\n\nGiven these clarifications, I raise my score from 2 to 3 (Borderline Reject). I think the core ideas have merit but require substantial revision before publication.", "rebuttal": "We appreciate the detailed review, but several fundamental misunderstandings about our method and factual errors about our experimental setup require correction.\n\n**Error bound claim**: The reviewer states we \"prove that using only 3–5 control points suffices to keep the action error $O(h^{k-1})$.\" This is **not** what we claim. Our theoretical analysis (Theorem 4, Appendix C) shows that **if** the path $\\theta(t)$ is $C^l$ with $l \\geq 2$, then the action error scales as $O(h^{\\kappa-1})$, where $\\kappa = \\min(l,4)$. Separately, we empirically observe that **3–5 control points suffice in practice** to capture high-quality paths across diverse problems. These are distinct claims - one theoretical, one empirical.\n\n**Dynamic/static equivalence assumptions**: Theorem 1 **does not** require regularity assumptions on $\\pi$. **As shown** in lines 137-139, we guarantee $\\tilde{\\Pi} \\neq \\emptyset$, it contains optimal transportation coupling $\\pi^{\\*}$ as the McCann interpolation is a diffeomorphic flow. Since **the traditional Kantorovich problem has no regularity requirements** on $\\pi$, neither does our static formulation.\n\n**Non-emptiness of $\\Gamma(\\pi)$**:  **We addressed this** for the Kantorovich coupling $\\pi^{\\*}$ in lines 137-139. We add a detailed discussion for a general $\\pi = (id,T_\\pi)\\\\#(\\rho_0).$ If we assume $T_\\pi$ is a diffeomorphism, then $\\Gamma(\\pi) \\neq \\emptyset$ iff the identity map $id$ and $T_\\pi$ are smoothly isotopic[1]. Two diffeomorphisms $f_0, f_1$ are smoothly isotopic if there exists a smooth homotopy $\\\\{f_t\\\\}_{0 \\leq t \\leq 1}$ where each $f_t$ is a diffeomorphism [1]. By Lemma 2, Chapter 6 [2], a map $T$ is smoothly isotopic to $id$ if $\\det(DT) > 0$ everywhere.\n\nWe define:\n$$\\mathcal{T}(\\rho_0,\\rho_1) = \\\\{T:\\mathbb{R}^d \\to \\mathbb{R}^d : T \\text{ diffeomorphic}, T\\\\#(\\rho_0) = \\rho_1, \\det(DT) > 0\\\\}$$\n\nSince the Monge map $T^{\\*}$ satisfies $\\det(DT^{\\*}) > 0$, we have $\\mathcal{T}(\\rho_0,\\rho_1) \\neq \\emptyset$. We modify $\\tilde{\\Pi}$ to:\n$$\\tilde{\\Pi} = \\\\{\\pi \\in \\Pi(\\rho_0,\\rho_1): \\pi = (id,T_\\pi)\\\\#(\\rho_0), T_\\pi \\in \\mathcal{T}(\\rho_0,\\rho_1)\\\\}$$\n\nThe proof of Theorem 1 remains unchanged with this new definition.\n\nWe thank the reviewer for raising this point. Providing a rigorous discussion of the non-emptiness of $\\Gamma(\\pi)$ allowed us to make our theoretical framework more precise, and we appreciate the opportunity to clarify this subtle but important assumption.\n\n***Neural ODE coverage**: The reviewer appears to misunderstand how NODEs are used in our framework. **Clarification** NODEs are used to define the parametric maps $T_\\theta$. Flows/paths are approximated by the curves $T_{\\theta(t)}(z)$ for $z \\sim \\lambda$, $t \\in [0,1]$.\n\n**Convergence guarantees**:  The reviewer demands convergence analysis for our bilevel optimization. We have two separate problems:\n1. Density space (11): For linear coercive actions, [4] Chapter 7 shows unique transport plans (Corollary 7.23)\n2. Parameter space (13): For linear coercive actions, we expect global minimizers given sufficient parametric expressiveness.\nConvergence analysis for non-linear potentials represents a fundamental open challenge affecting all neural approaches in optimal transport, not a limitation specific to our method.\nNumerically, our convergence tables show consistent action reduction:\n| Bilevel optimization cycles | Action value $\\mathcal{A}$ |\n| -------------------- | ------------- |\n| 0   | 18562.34 |\n| 7   |     476.81 |\n| 15   |  473.19  |\n| 23    |    471.40  |\n| 31    |   471.71    |\n We currently use maximum iterations as the stopping criteria. Common alternatives in bilevel optimization include gradient norm thresholds or loss change thresholds. Testing these could improve convergence detection and efficiency. In Figures 9-13, **we demonstrated** consistent action convergence across experiments with uncertainty bounds. \n\n**Baseline comparisons**: We compare NLOT and APAC-Net **when** applicable.  NLOT is **limited** to linear potentials while APAC-Net uses soft terminal costs instead of hard distributional constraints. GSBM is the **only** baseline capable of handling all settings we consider.  GSBM is the **SOTA**, the only method with equivalent flexibility to ours. The only other comparable method is DeepSBM [4], which GSBM already outperforms.\n\n**Implementation and computational fairness**: As stated in the paper we use **official implementations** with **recommended parameters** on identical hardware (AMD 7543 CPU + NVIDIA RTX A6000 GPU). To explore the computational fairness question, we compare the architectural differences:\n| Method | GMM | S-Curve | V-Neck | Opinion 2D | Opinion 1000D |\n|--------|-----|---------|--------|------------|---------------|\n| GSBM | 83,202 | 83,202 | 83,202 | 1,152,770 | 1,664,744 |\n| PDPO | 133,122 | 8,706 | 33,794 | 33,794 | 290,280 |\n\nPDPO achieves superior performance with simpler and often smaller networks. \n\nTo test whether parameter count explains PDPO's speed advantage, we solve S-Curve ($\\sigma = 0$) (with the GSBM setup) and V-Neck ($\\sigma = 1$) examples using **identical MLP  architectures** [2+1,128,4] for both methods. Both problems use $\\rho_0 = \\mathcal{N}(\\begin{bmatrix} -11 \\\\ -11 \\end{bmatrix},0.5I)$ and $\\rho_1 = \\mathcal{N}(\\begin{bmatrix} 11 \\\\ 11 \\end{bmatrix},0.5I).$ To show that our initialization scheme is agnostic to the functional, we use the same pretrained boundary parameters $\\theta_0, \\theta_1$ in both problems (training times: 0.33s and 0.35s, respectively). \n\n| Problem | Method | Time | $\\mathcal{A}$ | $W_2^2(t=0)$ | $W_2^2(t=1)$ |\n|---------|--------|------|-------------|--------------|--------------|\n| S-Curve | PDPO | 12m34s | 471.86 | 0.099 | 0.115 |\n| S-Curve | GSBM | 25m35s | 485.69 | 1.816 | 1.835 |\n| V-Neck | PDPO | 1h10m | 148.91 | 0.085 | 0.081 |\n| V-Neck | GSBM | 1h16m | 153.80 | 0.192 | 0.283 |\n\nWith identical architectures, PDPO still outperforms GSBM in both action value and boundary approximation quality, demonstrating that performance differences are not due to parameter count alone. This also illustrates boundary parameter reusability across different problems. Addressing the reviewer's fourth point. \n\n**Flow Matching initialization costs**: The complete FM training times are:\n| Distribution | Epochs | Time |\n|-------------|--------|------|\n| 2D Gaussians | 1k | ~40s each |\n| GMM | 2k | 3m10s |\n| 1000D Opinion | 4k | ~30m each |\nIncluding FM training time, PDPO still outperforms baselines significantly.\n\n**Boundary parameter sensitivity**: Once trained, boundary parameters **are robust and reusable across problems with different potentials**, as shown in this rebuttal. We treat initialization as \"outside\" because it's **general** and **amortized**. The same $\\theta_0$ can be reused for multiple optimization problems, as demonstrated in the table showing the performance of the models with the same architecture. We emphasize this is a feature, not a limitation. Furthermore, Appendix F.1 demonstrates successful optimization starting from $\\vec{0}$ parameters, achieving action value 40.31 in 24m32s without geodesic warmup. This shows our method's robustness to initialization.\n\n## Computational Complexity & Limitations\n\n**Fisher Information scaling**: The $O(N \\times M \\times (1+d))$ complexity for Fisher Information is fundamental to computing Fisher Information in any method, not unique to PDPO. As future work, we plan to explore Hutchinson-type estimators, see [5] for $\\nabla \\log(\\rho)$ to address this limitation.\n\n**Statistical significance**: Three seeds are a standard practice in this field (GSBM, DeepSBM). Additionally, our convergence curves with uncertainty estimates across optimization iterations provide strong evidence.\n\n\n**Flow Matching success conditions**: Flow Matching is SOTA in generative modeling [6], so we expect to find good boundary parameters with appropriate tuning. Our boundary optimization combines FM loss with action minimization, maintaining generation quality while improving coupling optimality. Furthermore, the FM component of our algorithm can be replaced with other models such as Stochastic Interpolations [7].\n\n**Theorem assumptions**: **All assumptions are explicitly** stated in Appendix C, Theorem 4. We detail smoothness requirements, Lipschitz continuity, bounded derivatives, and their relationship to action error bounds. We also provide corollaries for spline interpolation and Wasserstein error bounds.\n\n[1] Introduction to Differential Topology Joel W. Robbin and Dietman A. Salamon. Lecture course “Differential Geometry II” held by the\nsecond author at ETH Z¨urich in the spring semester of 2018.\n[2] From the Differentiable Viewpoint John W. Milnor. The University Press of Virginia, Charlottesville (1981).\n[3] Villani\n[4] Deep Generalized Schrödinger Bridge. Guan et al. 36th Conference on Neural Information Processing Systems \n[5]FFJORD: FREE-FORM CONTINUOUS DYNAMICS FOR SCALABLE REVERSIBLE GENERATIVE MODELS. Will Grathwohl et al. ICLR 2019\n[6] Flow Matching Guide and Code. Yarlon Lipman et al. Nips 2024\n[7]Stochastic Interpolants: A Unifying Framework for Flows and Diffusions CoRR 2023"}]

Author rebuttal:

Meta-review: This paper presents a new way to find the most efficient paths between probability densities, under the theoretical framework of optimal transport. It converts an infinite-dimensional optimization into a finite dim over parameters, using cubic splines to map those paths. The approach handles complex scenarios like obstacles, mean-field effects, and stochastic controls, and it performs better than current top methods on benchmarks.
The proposed algorithm PDPO demonstrates clear improvements in runtime compared to existing methods. 
As pointed out by a reviewer, the only baseline the authors claimed that can be compared with is GSBM. The usability of proposed method could be more convincing with comprehensive comparisons.
Three reviewers chose accept or borderline accept. One reviewer suggested borderline reject. After reading the paper and discussions, the AC tend to accept.
Decision: Accept (poster)