Learning of Population Dynamics: Inverse Optimization Meets JKO Scheme

## Abstract

Abstract Learning population dynamics involves recovering the underlying process that governs particle evolution, given evolutionary snapshots of samples at discrete time points. Recent methods frame this as an energy minimization problem in probability space and leverage the celebrated JKO scheme for efficient time discretization. In this work, we introduce iJKOnet , an approach that combines the JKO framework with inverse optimization techniques to learn population dynamics. Our method relies on a conventional end-to-end adversarial training procedure and does not require restrictive architectural choices, e.g., input-convex neural networks. We establish theoretical guarantees for our methodology and demonstrate improved performance over prior JKO-based methods.

## 1 Introduction

Figure: Figure 1: iJKOnet working scheme: our method minimizes the gap between the optimal values of parametric ($\theta$) JKO functional and suboptimal values obtained at ground truth population measures.
Refer to caption: images/method_overview.png

Modeling population dynamics is a fundamental challenge in many scientific domains, including biology (schiebinger2019optimal; moon2019visualizing), ecology (ayala1973competition), meteorology (fisher2009data; sigrist2015stochastic; verma2024climode; price2025probabilistic), transportation flows in urban networks (medina2022urban), and epidemiology (wang2021literature; kosma2023neural), among others. The task is to infer the underlying stochastic dynamics of a system – typically modeled by stochastic differential equations (SDEs) – from observed marginal distributions at discrete time points. While this problem has been studied extensively in settings where individual trajectories are available (krishnan2017structured; li2020scalable; brogat2024learning), such data are often unavailable in practice. In many real-world scenarios, we only observe population-level data, where it is infeasible to continuously track individual entities. Instead, we are forced to rely on temporally separated and mutually independent snapshots of the population.

In single-cell genomics (macosko2015highly), for example, measuring the state of individual cells typically involves destructive sampling – a process where cells are destroyed during measurement, preventing any further observation of their future behavior. As a result, the data consist only of isolated profiles of cellular populations taken at discrete time points. Reconstructing the continuous developmental trajectories of cells from such fragmented data poses a challenge. A similar situation arises in financial markets (gontis2010long; yang2023neural), where analysts often have access only to marginal distributions of asset prices at specific times, making it necessary to infer the underlying dynamics that govern these distributions. In crowd dynamics (maury2010macroscopic; maury2011handling; wan2023scalable), modeling the temporal evolution of population densities is crucial for understanding and managing pedestrian flows. In this case, data on individual trajectories are also typically unavailable, and only aggregate distributions at various time points are observed.

A promising approach to address the absence of particle trajectories is the Jordan–Kinderlehrer–Otto (JKO) scheme (jordan1998variational), which models the evolution of a particle system as a sequence of distributions that gradually approach the minimum of a total energy functional while remaining close to the previous distributions. However, implementing the JKO scheme involves solving an optimization problem over the space of probability measures, which is computationally demanding. The first attempt to leverage JKO for learning population dynamics was introduced in (bunne2022proximal, JKOnet). While innovative, this approach relies on a complex learning objective and is limited to potential energy functionals – meaning it cannot capture stochasticity in the dynamics. A recent work, (terpin2024learning, JKOnet^∗), proposes replacing the JKO optimization step with its first-order optimality conditions. This relaxation allows modeling more general energy functionals and reduces computational complexity, but it does not support end-to-end training. Instead, it requires precomputing optimal transport couplings (cuturi2013sinkhorn) between subsequent time snapshots, which limits its scalability and generalization.

Contributions. In this work, we present a novel approach for recovering the system dynamics based on observed population-level data. Our key contributions are as follows:

- 1.
Methodology: In \wasyparagraph[3.2](https://arxiv.org/html/2506.01502v2#S3.SS2), we cast the problem of reconstructing energy functionals within the JKO framework as an inverse optimization task. This perspective leads to a novel min-max optimization objective for population dynamics recovery.
- 2.
Algorithm: We equip our population dynamics methodology with a conventional end-to-end adversarial learning procedure (\wasyparagraph[3.3](https://arxiv.org/html/2506.01502v2#S3.SS3)). Importantly, our practical scheme does not pose restrictions on the architectures of the utilized neural networks, which contributes to the scalability.
- 3.
Theoretical Guarantees: Under suitable assumptions, we show that our method can accurately recover the underlying energy functional governing the observed population dynamics (\wasyparagraph [3.4](https://arxiv.org/html/2506.01502v2#S3.SS4)).

In \wasyparagraph[5](https://arxiv.org/html/2506.01502v2#S5) and Appendix [B](https://arxiv.org/html/2506.01502v2#A2), we evaluate our approach on a range of synthetic and real-world datasets, including single-cell genomics. The results show that our method demonstrates improved performance over previous JKO-based approaches for learning population dynamics.

Notation. Let ${\mathcal{X}}$ be a compact subset of $\mathbb{R}^{D}$ equipped with the Euclidean norm $\|\cdot\|_{2}$. Let ${\mathcal{P}}({\mathcal{X}})$ denote the set of probability measures on ${\mathcal{X}}$, and let ${\mathcal{P}}_{ac}({\mathcal{X}})$ denote its subset of probability measures absolutely continuous with respect to the Lebesgue measure. For $\rho\in{\mathcal{P}}_{ac}({\mathcal{X}})$, we use $\rho$ to denote both the measure and its density function with respect to the Lebesgue measure. For measures $\mu,\nu\in{\mathcal{P}}({\mathcal{X}})$, we denote the set of couplings (transportation plans) between them by $\Pi(\mu,\nu)$. For a measure $\rho\in{\mathcal{P}}({\mathcal{X}})$ and measurable map $T:{\mathcal{X}}\to{\mathcal{X}}$, we denote by $T\sharp\rho$ the associated push-forward measure, and $\text{id}_{\mathcal{X}}:{\mathcal{X}}\to{\mathcal{X}}$ is the identity mapping. $\nabla\cdot F=\sum_{d=1}^{D}\frac{\partial F_{d}}{\partial x_{d}}$ denotes the divergence operator for a continuously differentiable vector field $F:{\mathcal{X}}\to{\mathcal{X}}$.

## 2 Background

To describe the evolution of population measures, we use the theoretical framework of Wasserstein gradient flows (WGFs). In what follows, we first introduce the preliminary concept of Optimal Transport between probability measures (villani2008optimal; santambrogio2015optimal). Then we provide sufficient theoretical background on WGFs; see (santambrogio2017euclidean; FigalliGlaudo2023) for an overview and (ambrosio2008gradient) for a comprehensive study of WGFs theory. Finally, we get the reader acquainted with the JKO scheme, which is the cornerstone of our approach.

Optimal Transport. The (squared) Wasserstein-2 distance $d_{{\mathbb{W}}_{2}}$ between two probability measures $\mu,\nu\in{\mathcal{P}}({\mathcal{X}})$ is defined as the solution to the Kantorovich problem (kantorovich1942translocation):

$$ $d_{{\mathbb{W}}_{2}}^{2}(\mu,\nu)\stackrel{{\scriptstyle\text{def}}}{{=}}\min_{\pi\in\Pi(\mu,\nu)}\int_{{\mathcal{X}}\times{\mathcal{X}}}\|x-y\|^{2}_{2}\,\differential\pi(x,y),\vskip-1.42262pt$ (1) $$

where the distribution $\pi^{\ast}$ delivering the minimum to ([1](https://arxiv.org/html/2506.01502v2#S2.E1)) is called the optimal coupling (or plan). If $\mu\in{\mathcal{P}}_{ac}({\mathcal{X}})$, then the Brenier’s theorem (brenier1991polar) establishes equivalence between the Kantorovich formulation ([1](https://arxiv.org/html/2506.01502v2#S2.E1)) and Monge’s problem (monge1781memoire):

$$ $d_{{\mathbb{W}}_{2}}^{2}(\mu,\nu)=\min_{T:T\sharp\mu=\nu}\int_{{\mathcal{X}}}\|x-T(x)\|^{2}_{2}\,\differential\mu(x),\vskip-1.42262pt$ (2) $$

where the optimal transport map $T^{\ast}$ is known as the Monge map. Moreover, there exists a unique (up to an additive constant) convex potential $\psi^{\ast}:{\mathcal{X}}\to\mathbb{R}$ such that $T^{\ast}=\nabla\psi^{\ast}$ and $(\nabla\psi^{\ast})\sharp\mu=\nu$ (mccann1995existence). In this setting, the optimal coupling in ([1](https://arxiv.org/html/2506.01502v2#S2.E1)) is given by $\pi^{\ast}=[\text{id}_{{\mathcal{X}}},\nabla\psi^{\ast}]\sharp\mu$.

Wasserstein Gradient Flows. For clarity, we first consider the concept of gradient flows in Euclidean space $\mathbb{R}^{D}$ and then move to the space of probability measures ${\mathcal{P}}({\mathcal{X}})$. For Euclidean space $\mathbb{R}^{D}$, a gradient flow is an absolutely continuous curve $x(t)$ starting at $x_{0}\in\mathbb{R}^{D}$ that minimizes an energy functional $J:\mathbb{R}^{D}\to\mathbb{R}$ ”as fast as possible”. To find such a curve, one needs to solve an ordinary differential equation (ODE) (teschl2012ordinary) of the form $x^{\prime}(t)=-\nabla J(x(t))$ with initial condition $x(0)=x_{0}$. The same idea is applicable to the space of probability measures ${\mathcal{P}}({\mathcal{X}})$. If we equip ${\mathcal{P}}({\mathcal{X}})$ with the Wasserstein-2 distance $d_{{\mathbb{W}}_{2}}$, we obtain the Wasserstein space ${\mathbb{W}}_{2}({\mathcal{X}})=({\mathcal{P}}({\mathcal{X}}),d_{{\mathbb{W}}_{2}})$ – a complete and separable metric space (bogachev2012monge). In this case, for an energy functional ${\mathcal{J}}:{\mathcal{P}}({\mathcal{X}})\to\mathbb{R}$, the gradient flow in ${\mathbb{W}}_{2}({\mathcal{X}})$, called the Wasserstein gradient flow (WGF), is an absolutely continuous curve $\rho_{t}:\mathbb{R}_{+}\to{\mathcal{P}}({\mathcal{X}})$ starting at $\rho_{0}$ that follows the steepest descent direction of ${\mathcal{J}}$, i.e., solves

$$ $\partial_{t}\rho_{t}=-\nabla_{{\mathbb{W}}_{2}}{\mathcal{J}}(\rho_{t}),\quad\rho_{(t=0)}=\rho_{0},\vskip-1.42262pt$ (3) $$

where $\nabla_{{\mathbb{W}}_{2}}{\mathcal{J}}(\rho_{t})$ denotes the Wasserstein gradient in ${\mathbb{W}}_{2}$ given by $\nabla_{{\mathbb{W}}_{2}}{\mathcal{J}}(\rho)=-\nabla\cdot(\rho\nabla\frac{\delta{\mathcal{J}}}{\delta\rho})$ with $\frac{\delta{\mathcal{J}}}{\delta\rho}$ denoting the first variation (chewi2025) of energy ${\mathcal{J}}$. Thus, equation ([3](https://arxiv.org/html/2506.01502v2#S2.E3)) can be rewritten in the form of the continuity equation, expressing mass conservation under the velocity field $v_{t}$:

$$ $\partial_{t}\rho_{t}+\nabla\cdot\left(\rho_{t}v_{t}\right)=0,\;v_{t}=-\nabla\frac{\delta{\mathcal{J}}}{\delta\rho}(\rho_{t}).\vskip-1.42262pt$ (4) $$

There exists an intriguing connection between WGFs and certain partial differential equations (PDEs) (evans2022partial). In particular, different energy functionals in the Wasserstein space correspond to distinct PDEs (santambrogio2015optimal), some of which are associated with diffusion processes appearing in practice (gomez2024beginner; bailo2024aggregation).

Examples of PDEs as WGFs. Consider the free energy functional (carrillo2003kinetic, Eq. (1.3)):

$$ ${\mathcal{J}}_{\text{FE}}(\rho)=\underbrace{\int_{\mathcal{X}}V(x)\,\differential\rho(x)}_{{\mathcal{V}}(\rho)}+\underbrace{\int_{{\mathcal{X}}\times{\mathcal{X}}}W(x-y)\,\differential\rho(x)\differential\rho(y)}_{{\mathcal{W}}(\rho)}+\underbrace{\int_{\mathcal{X}}U(\rho(x))\,\differential x}_{{\mathcal{U}}(\rho)},\vskip-1.42262pt$ (5) $$

where ${\mathcal{V}}$, ${\mathcal{W}}$, and ${\mathcal{U}}$ correspond to the system’s potential, interaction, and internal energies, respectively. The PDE corresponding to the WGF driven by this functional is known as the aggregation-diffusion equation (ADE). It describes the evolution of density $\rho_{t}$ under a corresponding velocity field $v_{t}$ ([4](https://arxiv.org/html/2506.01502v2#S2.E4)). This evolution reflects a balance of three effects: drift – driven by the potential $V(x)$, modeling an external field (e.g., gravity, electric potential); interaction – governed by the symmetric interaction kernel $W(x-y)$, accounting for non-local effects (e.g., particle interactions, long-range forces); and diffusion – modeled by the internal energy ${\mathcal{U}}(\rho)$, representing the energy associated with the local state of the system (e.g., thermal energy, chemical energy). A weak solution to the equation exists if the no-flux conditions: $\nabla V\cdot\vec{n}=0$ on $\partial{\mathcal{X}}$ holds, where $\vec{n}$ is the outward normal vector on the boundary, and if $W\in{\mathcal{C}}({\mathcal{X}}\times{\mathcal{X}}$), i.e. $W$ is a continuous function on ${\mathcal{X}}\times{\mathcal{X}}$. These conditions ensure that no mass crosses the boundary of ${\mathcal{X}}$; mass can only be redistributed within ${\mathcal{X}}$.
For a detailed discussion of the existence and classification of solutions, see (gomez2024beginner).

Such energy functionals are ubiquitous in real-world applications, including physics (carrillo2021phase), biology (keller1971model; potts2024distinguishing; potts2024aggregation), economics (fiaschi2025spatial), machine learning (suzuki2023uniformintime; chizat2024infinite; nitanda2024improved), and nonlinear optimization (bailo2024cbx), to name a few. For further discussion and references, see the recent surveys (carrillo2019aggregation; gomez2024beginner; bailo2024aggregation).

To demonstrate the flexibility of the energy formulation in ([5](https://arxiv.org/html/2506.01502v2#S2.E5)), we examine a few representative cases. When the system’s energy is purely internal and given by the negative differential entropy, ${\mathcal{U}}(\rho)=-{\mathcal{H}}(\rho)\stackrel{{\scriptstyle\text{def}}}{{=}}\int\rho(x)\log\rho(x)\differential x$, continuity equation ([4](https://arxiv.org/html/2506.01502v2#S2.E4)) reduces to the classical heat equation: $\partial_{t}\rho=\nabla^{2}\rho$ (vazquez2017asymptotic). Alternatively, when the energy is ${\mathcal{J}}_{\text{FP}}(\rho)={\mathcal{V}}(\rho)-\beta{\mathcal{H}}(\rho)$, the resulting PDE is the linear Fokker-Planck equation with diffusion coefficient $\beta$:

$$ $\partial_{t}\rho_{t}=\nabla\cdot\left(\nabla V(x)\,\rho_{t}\right)+\beta\nabla^{2}\rho_{t},\vskip-1.42262pt$ (6) $$

A fruitful branch of theoretical and practical research stems from the connection between the Fokker-Planck PDE and stochastic differential equations (SDEs) (risken1996fokker; bogachev2022fokker), the latter describes the stochastic evolution of particles. In particular, equation ([6](https://arxiv.org/html/2506.01502v2#S2.E6)) is equivalent (weinan2021applied) to the following Itô SDE:

$$ $\differential X_{t}=-\nabla V(X_{t})\,\differential t+\sqrt{2\beta}\,\differential W_{t},\vskip-1.42262pt$ (7) $$

where $X=\{X_{t}\}_{t\geq 0}$ is a $\mathbb{R}^{D}$-valued stochastic process and $W=\{W_{t}\}_{t\geq 0}$ is a standard Wiener process (sarkka2019applied). In other words, if $X_{t}\sim\rho_{t}$ evolves according to ([7](https://arxiv.org/html/2506.01502v2#S2.E7)), then the density $\rho_{t}$ evolves according to the Fokker-Planck equation ([6](https://arxiv.org/html/2506.01502v2#S2.E6)) in the space of probability measures.

WGFs provide a compelling framework for modeling PDEs and their associated SDEs across various domains, giving rise to a range of methods for approximating solutions to WGFs. These include deep learning-based approaches (mokrov2021large; alvarez2021optimizing; pmlr-v202-altekruger23a), particle-based methods such as Stein Variational Gradient Descent (SVGD) (liu2016stein) and its extensions (das2023provably; tankala2025beyond), as well as classical discretization techniques in Wasserstein space (carrillo2022primal, see \wasyparagraph1.2). Many of these advances build on the celebrated Jordan-Kinderlehrer-Otto (JKO) scheme (jordan1998variational), which we review next.

The JKO Scheme. A classical method for solving ODEs in Euclidean space is the implicit Euler scheme. For an (Euclidean) gradient flow $x^{\prime}(t)=-\nabla J(x(t))$ and $\tau>0$, this scheme approximately solves this ODE by iteratively applying the proximal operator (parikh2014proximal, \wasyparagraph1.1): ${x_{k+1}=\operatorname{prox}_{\tau J}(x_{k})=\operatorname*{arg\,min}_{x\in\mathbb{R}^{D}}\!\big\{J(x)+\frac{1}{2\tau}\|x-x_{k}\|_{2}^{2}\big\}}$. Compared to the explicit Euler scheme, this approach offers improved numerical stability (butcher2016numerical). Jordan, Kinderlehrer, and Otto (jordan1998variational) extended this idea to the space of probability measures, introducing a variational time discretization of the Fokker–Planck equation ([6](https://arxiv.org/html/2506.01502v2#S2.E6)), now known as the JKO scheme:

$$ $\rho^{\tau}_{k+1}=\operatorname*{arg\,min}_{\rho\in{\mathcal{P}}({\mathcal{X}})}\left\{{\mathcal{J}}(\rho)+\dfrac{1}{2\tau}d_{{\mathbb{W}}_{2}}^{2}(\rho,\rho^{\tau}_{k})\right\}=\operatorname{JKO}_{\tau{\mathcal{J}}}(\rho^{\tau}_{k}),\quad\rho^{\tau}_{0}=\rho_{0},\vskip-1.42262pt$ (8) $$

where $\tau>0$ is the time step. As $\tau\to 0$, the sequence ${\rho^{\tau}_{k}}{,k\in{\mathbb{N}}}$ converges to the continuous solution $\rho_{t}$ of the Fokker–Planck equation ([6](https://arxiv.org/html/2506.01502v2#S2.E6)). Later, the JKO scheme’s convergence was generalized (ambrosio2008gradient, Thm. 4.0.4) for the free energy functional ${\mathcal{J}}_{\text{FE}}$ ([5](https://arxiv.org/html/2506.01502v2#S2.E5)). Note that since $d_{{\mathbb{W}}_{2}}^{2}(\rho,\rho^{\tau}_{k})\geq 0$, the energy is non-increasing along the sequence: ${\mathcal{J}}(\rho^{\tau}_{k+1})\leq{\mathcal{J}}(\rho^{\tau}_{k})$. This monotonicity property is often utilized in JKO-based algorithms (salim2020wasserstein) in Wasserstein space.

## 3 iJKOnet Method

We begin by formally stating the problem addressed in our work (\wasyparagraph[3.1](https://arxiv.org/html/2506.01502v2#S3.SS1)), then we develop our methodology (\wasyparagraph[3.2](https://arxiv.org/html/2506.01502v2#S3.SS2)). Finally, we discuss the practical (\wasyparagraph[3.3](https://arxiv.org/html/2506.01502v2#S3.SS3)) and theoretical aspects (\wasyparagraph[3.4](https://arxiv.org/html/2506.01502v2#S3.SS4)) of our approach.

### 3.1 Problem statement

In our work, we address the problem of recovering the underlying energy functional ${\mathcal{J}}^{\ast}$ that governs the evolution of a density $\rho_{t}\in{\mathcal{P}}_{ac}({\mathcal{X}})$ in ([4](https://arxiv.org/html/2506.01502v2#S2.E4)) based on marginal population measures (bunne2022proximal; lavenant2024toward). Specifically, we are given independent samples from marginals $\{\rho_{k}\}_{k=0}^{K}$ of the evolving distribution $\rho_{t}$ at corresponding time points $t_{0}<t_{1}<\dots<t_{K}$. Importantly, each distribution $\rho_{k}$ may be represented by a different number of samples. As discussed in \wasyparagraph[2](https://arxiv.org/html/2506.01502v2#S2), the JKO scheme approaches the continuous-time dynamics as the step sizes $\Delta t_{k}=t_{k+1}-t_{k}$ are small enough. This motivates our modeling assumption that the ground truth sequence of measures $\{\rho_{k}\}_{k=0}^{K}$ follows the scheme $\rho_{k+1}=\operatorname{JKO}_{\Delta t_{k}{\mathcal{J}}^{\ast}}(\rho_{k})$.

Although the time intervals $\Delta t_{k}$ between observations may vary, corresponding to non-uniform step sizes in the JKO scheme, we assume equal spacing $\Delta t_{k}=\tau$ in the remaining text for simplicity.

### 3.2 Methodology and loss derivation

In this section, we assume $\rho_{k}\in{\mathcal{P}}_{ac}({\mathcal{X}})$ and minimization is always taken over $\rho\in{\mathcal{P}}_{ac}({\mathcal{X}})$ in order to utilize Brenier’s theorem. The key idea of our method builds on inverse optimization techniques (chan2025inverse) and is demonstrated in Figure [1](https://arxiv.org/html/2506.01502v2#S1.F1). Thanks to assumption $\rho_{k+1}=\operatorname{JKO}_{\tau{\mathcal{J}}^{\ast}}(\rho_{k})$, we can derive an inequality that becomes an equality if a candidate functional ${\mathcal{J}}$ matches the ground truth functional ${\mathcal{J}}^{\ast}$, since $\rho_{k+1}$ is obtained via the JKO step associated with ${\mathcal{J}}^{\ast}$:

$$ $\min_{\rho}\left\{{\mathcal{J}}(\rho)+\frac{1}{2\tau}d_{{\mathbb{W}}_{2}}^{2}(\rho_{k},\rho)\right\}\leq{\mathcal{J}}(\rho_{k+1})+\frac{1}{2\tau}d_{{\mathbb{W}}_{2}}^{2}(\rho_{k},\rho_{k+1}).$ (9) $$

Moving the right-hand side to the left yields an expression that is always upper-bounded by zero, regardless of the choice of ${\mathcal{J}}$. Maximizing the resulting gap with respect to ${\mathcal{J}}$ encourages the candidate functional to approximate the true functional ${\mathcal{J}}^{\ast}$. This yields the following objective:

$$ $\displaystyle\max_{{\mathcal{J}}}\sum_{k=0}^{K-1}\Big[$ $\displaystyle\min_{\rho}\Big\{{\mathcal{J}}(\rho)+\frac{1}{2\tau}d_{{\mathbb{W}}_{2}}^{2}(\rho_{k},\rho)\Big\}-{\mathcal{J}}(\rho_{k+1})-\underbrace{\frac{1}{2\tau}d_{{\mathbb{W}}_{2}}^{2}(\rho_{k},\rho_{k+1})}_{{\text{independent of ${\mathcal{J}}$, $\rho$}}}\Big]=$ $\displaystyle\max_{{\mathcal{J}}}\sum_{k=0}^{K-1}$ $\displaystyle\min_{\rho}\Big[{\mathcal{J}}(\rho)-{\mathcal{J}}(\rho_{k+1})+\frac{1}{2\tau}d_{{\mathbb{W}}_{2}}^{2}(\rho_{k},\rho)\Big]+\text{Const}.\vskip-2.84526pt$ (10) $$

Note that the distributions $\rho$ in ([10](https://arxiv.org/html/2506.01502v2#S3.E10)) can equivalently be expressed via transport maps $\rho=T^{k}\sharp\rho_{k}$, where $T^{k}:{\mathcal{X}}\to{\mathcal{X}}$. Furthermore, thanks to ([2](https://arxiv.org/html/2506.01502v2#S2.E2)), we have $d_{{\mathbb{W}}_{2}}^{2}(\rho_{k},\rho)\leq\int_{{\mathcal{X}}}\|x-T^{k}(x)\|^{2}d\rho_{k}(x)$. Bringing these two observations together, we end up with our final loss objective equivalent to ([10](https://arxiv.org/html/2506.01502v2#S3.E10)):

$$ $\!\!\max_{{\mathcal{J}}}\min_{T^{k}}{\mathcal{L}}({\mathcal{J}},T^{k})\!\stackrel{{\scriptstyle\text{def}}}{{=}}\!\max_{{\mathcal{J}}}\min_{T^{k}}\!\sum_{k=0}^{K-1}\Big[{\mathcal{J}}(T^{k}\sharp\rho_{k})-{\mathcal{J}}(\rho_{k+1})+\frac{1}{2\tau}\int_{\mathcal{X}}\|x-T^{k}(x)\|^{2}_{2}\rho_{k}(x)\,\differential x\Big].$ (11) $$

The key idea is that the inner minimization over maps $T^{k}$ corresponds to performing a JKO step for a given energy functional ${\mathcal{J}}$; that is, each optimal map $T^{k}_{{\mathcal{J}}}$ pushes $\rho_{k}$ to the JKO-updated distribution $\hat{\rho}_{k+1}=\operatorname{JKO}_{\tau{\mathcal{J}}}(\rho_{k})$. The outer maximization then drives $\mathcal{J}$ toward the true functional $\mathcal{J}^{\ast}$, causing the pushforward distributions $\hat{\rho}_{k+1}=T^{k}_{{\mathcal{J}}}\sharp\rho_{k}$ to approach $\rho_{k+1}$, see illustration in Figure [1](https://arxiv.org/html/2506.01502v2#S1.F1).

### 3.3 Practical aspects: method parametrization and learning procedure

We denote by $\theta\in\Theta$ and $\varphi\in\Phi$ the parameters of sufficiently expressive function classes used to approximate the candidate functional ${\mathcal{J}}$ and the transport maps $T^{k}$, respectively. Specifically, ${\mathcal{J}}_{\theta}$ and $T_{\varphi}^{k}$ are their parameterized counterparts.

Mapping Parameterization. Building upon (benamou2016discretization), prior works (mokrov2021large; alvarez2021optimizing; bunne2022proximal) typically parameterize the transport maps as $T_{\varphi}^{k}=\nabla\psi_{\varphi}^{k}$, where $\psi_{\varphi}^{k}$ are modeled using input-convex neural networks (ICNNs) (amos2017input). However, ICNNs suffer from poor scalability in high-dimensional settings (korotin2021neural). In contrast, since our objective ([11](https://arxiv.org/html/2506.01502v2#S3.E11)) imposes no convexity constraints, we parameterize $T_{\varphi}^{k}$ directly using standard architectures like MLPs or ResNets (he2016deep). This relaxation simplifies optimization and yields improved empirical stability.

Energy functional parametrization. Following (terpin2024learning), we parameterize the energy functional ${\mathcal{J}}$ as a specific instance of the free energy formulation in ([5](https://arxiv.org/html/2506.01502v2#S2.E5)), where the internal energy term is chosen to be scaled negative differential entropy:

$$ ${\mathcal{J}}_{\theta}(\rho)=\underbrace{\int_{\mathcal{X}}V_{\theta_{1}}(x)\,\differential\rho(x)}_{{\mathcal{V}}_{\theta_{1}}(\rho)}+\underbrace{\int_{{\mathcal{X}}\times{\mathcal{X}}}W_{\theta_{2}}(x-y)\,\differential\rho(x)\differential\rho(y)}_{{\mathcal{W}}_{\theta_{2}}(\rho)}+\underbrace{\theta_{3}\int_{{\mathcal{X}}}\log\rho(x)\,\differential\rho(x)}_{{\mathcal{U}}_{\theta_{3}}(\rho)=-\theta_{3}{\mathcal{H}}(\rho)}.$ (12) $$

Here $\theta=\{\theta_{1},\theta_{2},\theta_{3}\}$ represents the set of learnable parameters: $\theta_{1}$ and $\theta_{2}$ are parameters of neural networks that define the potential $V_{\theta_{1}}:{\mathcal{X}}\to\mathbb{R}$ and the interaction kernel $W_{\theta_{2}}:{\mathcal{X}}\to\mathbb{R}$, respectively; $\theta_{3}\in{\mathbb{R}}$ is a learnable scalar diffusion coefficient, see ([7](https://arxiv.org/html/2506.01502v2#S2.E7)).

Entropy evaluation. All the terms in the objective function ([11](https://arxiv.org/html/2506.01502v2#S3.E11)) can be estimated using Monte Carlo integration (metropolis1949monte), except for the internal energy term ${\mathcal{U}}_{\theta_{3}}$. To handle this remaining component, we follow (mokrov2021large, Theorem 1), (alvarez2021optimizing, Eq. (11)) and apply the change-of-variables formula, yielding:

$$ ${\mathcal{U}}_{\theta_{3}}(T_{\varphi}^{k}\sharp\rho_{k})={\mathcal{U}}_{\theta_{3}}(\rho_{k})-\theta_{3}\int_{{\mathcal{X}}}\log|\det\nabla_{x}T_{\varphi}^{k}(x)|\,\differential\rho_{k}(x),$ (13) $$

where $\nabla_{x}T_{\varphi}^{k}$ denotes the Jacobian of the invertible transport map $T_{\varphi}^{k}$. We note that, during training, the computation of $\log\det$ – even for non-invertible mappings $T_{\varphi}^{k}$ – can be effectively handled by modern programming tools; see Appendix [C.3](https://arxiv.org/html/2506.01502v2#A3.SS3) for details. To efficiently compute the gradient of the $\log\det$ term in ([13](https://arxiv.org/html/2506.01502v2#S3.E13)) with respect to $\varphi$, one can use the Hutchinson trace estimator (hutchinson1989stochastic), as proposed by (finlay2020train; huang2021convex). However, in our experimental setting (\wasyparagraph[5](https://arxiv.org/html/2506.01502v2#S5)), we find that computing the full Jacobian is tractable and sufficient. For estimating the negative differential entropy ${\mathcal{H}}(\rho_{k})$ in ([13](https://arxiv.org/html/2506.01502v2#S3.E13)), we employ the Kozachenko–Leonenko nearest-neighbor estimator (kozachenko1987sample; berrett2019efficient). Notably, entropy values for the population measures $\{\rho_{k}\}_{k=0}^{K}$ can be precomputed prior to training.

Learning Objective. To facilitate gradient-based optimization, we employ Monte Carlo estimation to approximate the loss in ([11](https://arxiv.org/html/2506.01502v2#S3.E11)). At each time step $t_{k}$ (for $k=0,\dots,K$), we draw a batch of $B_{k}$ samples $X_{k}=\{x^{1}_{k},\dots,x^{B_{k}}_{k}\}\sim\rho_{k}$ and optimize the following empirical loss objective:

$$ $\!\!\!\!\widehat{{\mathcal{L}}}(\theta,\varphi)\!\!=\!\!\sum_{k=0}^{K-1}\!\left[\widehat{{\mathcal{J}}}_{\theta}(T_{\varphi}^{k}(X_{k}))-\widehat{{\mathcal{J}}}_{\theta}(X_{k+1})+\frac{1}{2\tau}\sum_{j=1}^{B_{k}}\|x^{j}_{k}-T_{\varphi}^{k}(x^{j}_{k})\|^{2}_{2}\right]\!,$ (14) $$

where $T_{\varphi}^{k}(X_{k})$ denotes the batched pushforward used to compute the empirical estimates $\widehat{{\mathcal{J}}}_{\theta}$ of the functional ${\mathcal{J}}_{\theta}$. We optimize the objective in ([14](https://arxiv.org/html/2506.01502v2#S3.E14)) with respect to $\theta$ and $\varphi$ using a standard gradient descent–ascent scheme; implementation details are provided in \wasyparagraph[5](https://arxiv.org/html/2506.01502v2#S5) and Appendix [C](https://arxiv.org/html/2506.01502v2#A3), with the complete training procedure summarized in Algorithm [1](https://arxiv.org/html/2506.01502v2#algorithm1) therein.

### 3.4 Theoretical aspects

The central idea behind our loss formulation (\wasyparagraph[3.2](https://arxiv.org/html/2506.01502v2#S3.SS2)) is to minimize the inverse gap in ([9](https://arxiv.org/html/2506.01502v2#S3.E9)), which measures the discrepancy between the optimal value (left-hand side) and the expected value (right-hand side) of the JKO step for a candidate functional ${\mathcal{J}}$. While this objective is intuitively justified, it remains to be shown whether the minimizer ${\mathcal{J}}_{\text{min}}$ of ([11](https://arxiv.org/html/2506.01502v2#S3.E11)) truly approximates the ground-truth energy functional ${\mathcal{J}}^{\ast}$. Our Theorem [3.1](https://arxiv.org/html/2506.01502v2#S3.Thmtheorem1) addresses this question, demonstrating that, under suitable assumptions, optimizing ([11](https://arxiv.org/html/2506.01502v2#S3.E11)) indeed recovers ${\mathcal{J}}^{\ast}$ up to an additive constant that does not affect the dynamics governed by ([4](https://arxiv.org/html/2506.01502v2#S2.E4)). Full details and proofs are provided in Appendix [D](https://arxiv.org/html/2506.01502v2#A4).

We state our theorem for $K=1$. When considering $K>1$, the statement of the theorem holds true independently for each timestep $k=0,1,\dots,K$, and the results extend straightforwardly to this more general setting. We further assume that both the ground-truth functional ${\mathcal{J}}^{\ast}$ and the candidate functional ${\mathcal{J}}$ are purely of potential energy form, i.e., ${\mathcal{J}}^{\ast}(\rho)={\mathcal{J}}^{\ast}_{\text{PE}}(\rho)=\int_{{\mathcal{X}}}V^{\ast}(x)\,\differential\rho(x)$; ${\mathcal{J}}(\rho)={\mathcal{J}}_{\text{PE}}(\rho)=\int_{{\mathcal{X}}}V(x)\,\differential\rho(x)$. For notational simplicity, we denote the loss as ${\mathcal{L}}(V,T)$ in place of ${\mathcal{L}}({\mathcal{J}}_{\text{PE}},T)$. For technical purposes, we introduce the modified version of a potential as $V_{q}\coloneqq\tau V+\tfrac{1}{2}\|\cdot\|_{2}^{2}:{\mathcal{X}}\to\mathbb{R}$; subscript $q$ stands for “quadratic“.

###### Theorem 3.1 (Quality bounds for recovered potential energy) .

Let $\varepsilon(V)\stackrel{{\scriptstyle\text{def}}}{{=}}{\mathcal{L}}(V^{\ast},T_{V^{\ast}})-{\mathcal{L}}(V,T_{V})$ be the gap between the optimal and optimized value of inverse JKO loss ([11](https://arxiv.org/html/2506.01502v2#S3.E11)) with internal $\min_{T}$ problem solved exactly, i.e., $T_{V}\stackrel{{\scriptstyle\text{def}}}{{=}}\min_{T}{\mathcal{L}}(V,T)$. Let ${\mathcal{X}}$ be a convex set; (modified) potentials $V_{q}$ be strictly convex and $\tfrac{1}{\beta}$-smooth (see the definitions in Appendix [D](https://arxiv.org/html/2506.01502v2#A4)). Then there exists $C\!=\!C(\tau,\beta)$ such that following inequality holds:

$$ $\int_{{\mathcal{X}}}\big\|\nabla V^{\ast}(y)-\nabla V(y)\big\|^{2}\differential\rho_{1}(y)\leq C\varepsilon(V).$ (15) $$

Notably, equation ([15](https://arxiv.org/html/2506.01502v2#S3.E15)) compares gradients of the recovered and ground truth potentials, hiding the appearance of a redundant additive constant. Importantly, the assumptions of the theorem are not particularly restrictive in practice. Specifically, the smoothness of potentials can be ensured by employing smooth activation functions such as CELU (barron2017continuously), SiLU (hendrycks2016gaussian), SoftPlus, and others. Strict convexity, in turn, can be enforced through architectural design choices (amos2017input). Furthermore, the strict convexity of the (modified) potentials $V_{q}$ can often be assumed when the step size $\tau$ is sufficiently small. In our experiments (\wasyparagraph[5](https://arxiv.org/html/2506.01502v2#S5)), we parameterize $V$ using standard MLPs and observe that this approach performs adequately.

To the best of our knowledge, this work is the first to provide a quality analysis (Theorem [3.1](https://arxiv.org/html/2506.01502v2#S3.Thmtheorem1)) of JKO-based solvers for population dynamics. At present, our analysis focuses exclusively on the potential energy component. Extending this framework to incorporate interaction and internal energies presents an interesting direction for future research.

###### Theorem 3.1 (Quality bounds for recovered potential energy) .

## 4 Related Works

This section reviews research directions most relevant to our work. For an extended discussion of related works, see Appendix [A.1](https://arxiv.org/html/2506.01502v2#A1.SS1). First, we discuss existing methods for modeling the dynamics of $\rho_{k}$ given a known energy functional ${\mathcal{J}}$ (\wasyparagraph[4.1](https://arxiv.org/html/2506.01502v2#S4.SS1)). Next, we focus on approaches closely related to our setting – specifically, learning population dynamics from observed data (\wasyparagraph[4.2](https://arxiv.org/html/2506.01502v2#S4.SS2)) via the JKO scheme.

### 4.1 Solving aggregation-diffusion equation with given energy

This line of research focuses on methods for solving aggregation-diffusion equations derived from WGFs with given energy functionals of the form ([5](https://arxiv.org/html/2506.01502v2#S2.E5)). A variety of numerical approaches exist, including mesh-based schemes (cances2023finite; jungel2024convergent), particle methods (craig2016blob; campos2018convergence; carrillo2019blob), and JKO-based schemes (carrillo2022primal). For a comprehensive overview, see (bailo2024aggregation, \wasyparagraph5). Recently, the research field was empowered by deep learning techniques. In particular, (alvarez2021optimizing; mokrov2021large; fan2022variational; park2023deep; lee2024deep) employ gradient-based optimization and neural-network parametrization to solve WGFs with given energy ${\mathcal{J}}$ as to special cases of ([5](https://arxiv.org/html/2506.01502v2#S2.E5)).

### 4.2 Learning energy for aggregation-diffusion equation via JKO scheme

This line of research focuses on methods for learning population dynamics (\wasyparagraph[3.1](https://arxiv.org/html/2506.01502v2#S3.SS1)) using the JKO scheme, as in JKOnet (bunne2022proximal) and JKOnet^∗ (terpin2024learning). We discuss these methods in detail below, as they are the most closely related to our approach. Additional methods for learning population dynamics are reviewed in Appendix [A.1.1](https://arxiv.org/html/2506.01502v2#A1.SS1.SSS1).

JKOnet (bunne2022proximal) formulates the task of population dynamics recovery as a bi-level optimization problem aimed at minimizing the discrepancy between observed distributions $\rho_{k}$ and model predictions $\hat{\rho}_{k}$:

$$ $\displaystyle{\mathcal{L}}_{\texttt{JKOnet}}(\theta,\varphi)$ $\displaystyle=\sum_{k=0}^{K-1}d_{{\mathbb{W}}_{2}}^{2}(\hat{\rho}_{k},\rho_{k}),\quad\text{s.t. }\hat{\rho}_{0}=\rho_{0},\quad\hat{\rho}_{k+1}=\nabla\psi_{k}^{\ast}\sharp\hat{\rho}_{k},$ (16) $\displaystyle\psi_{k}^{\ast}$ $\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname*{arg\,min}_{\varphi:\psi_{\varphi}\in\texttt{CVX}}\begin{aligned} {\mathcal{J}}_{\theta}(\nabla\psi_{\varphi}\sharp\hat{\rho}_{k})+\frac{1}{2\tau}\int_{{\mathcal{X}}}\norm{x-\nabla\psi_{\varphi}(x)}^{2}\,\differential\hat{\rho}_{k}(x),\end{aligned}$ $$

where CVX denotes the set of continuously differentiable convex functions from ${\mathcal{X}}$ to $\mathbb{R}$. The transport maps $\psi_{\varphi}$ are parametrized either using ICNNs (bunne2022supervised) or through the ‘Monge gap’ approach (uscidda2023monge). However, the energy functional is limited to the potential energy term ${\mathcal{V}}_{\theta}(\rho)$, excluding interaction and internal components. This framework has two key limitations: (i) solving the bi-level optimization problem in ([16](https://arxiv.org/html/2506.01502v2#S4.E16)) is computationally challenging and requires specific techniques like unrolling optimizer’s steps, and (ii) extending the method to more realistic settings with richer energy structures is questionable due to computational complexity.

JKOnet^∗ (terpin2024learning) addresses the limitations of the original JKOnet by replacing the full JKO optimization problem ([8](https://arxiv.org/html/2506.01502v2#S2.E8)) with its first-order optimality conditions. This reformulation allows the use of more expressive energy functionals ([12](https://arxiv.org/html/2506.01502v2#S3.E12)). The method minimizes the following objective, utilizing a precomputed optimal transport plans $\pi_{k}$ between $\rho_{k}$ and $\rho_{k+1}$:

$$ $\displaystyle{\mathcal{L}}_{\texttt{JKOnet${}^{\ast}$}}(\theta)=\sum_{k=0}^{K-1}\int_{{\mathcal{X}}\times{\mathcal{X}}}\bigg\|$ $\displaystyle\nabla V_{\theta_{1}}(x_{k+1})+\int_{{\mathcal{X}}}\nabla U_{\theta_{2}}(x_{k+1}-y_{k+1})\,\differential\rho_{k+1}(y_{k+1})$ (17) $\displaystyle+$ $\displaystyle\theta_{3}\frac{\nabla\rho_{k+1}(x_{k+1})}{\rho_{k+1}(x_{k+1})}+\frac{1}{\tau}(x_{k+1}-x_{k})\bigg\|^{2}\,\differential\pi_{k}(x_{k},x_{k+1}).$ $$

The authors emphasize several advantages over JKOnet, including reduced computational costs and support for more general energy functionals. However, their method requires an additional optimization round to precompute optimal transport plans $\pi_{k}$, introducing extra sources of inaccuracy and rendering the approach non-end-to-end. Moreover, the authors of (terpin2024learning) employ discrete OT solvers to compute $\pi_{k}$, which may fail to accurately represent OT maps between the underlying distributions, particularly in high-dimensional settings (deb2021rates).

iJKOnet (Ours) integrates the strengths of both JKOnet and JKOnet^∗. From JKOnet, it inherits parameterized models for the energy functional ${\mathcal{J}}_{\theta}$ and the transport map $T_{\varphi}$, avoiding reliance on precomputed discrete optimal transport plans $\pi_{k}$, which can hinder flexible mappings between $\rho_{k}$ and $\rho_{k+1}$. From JKOnet^∗, it adopts a more expressive formulation of the energy functional ([12](https://arxiv.org/html/2506.01502v2#S3.E12)) – capable of capturing complex structures, with potential extensions to porous medium energies (alvarez2021optimizing), although such extensions are beyond the scope of the current work.

## 5 Experiments

Our implementation is primarily based on the publicly available code(^1^11[https://github.com/antonioterpin/jkonet-star](https://github.com/antonioterpin/jkonet-star)) from JKOnet^∗ (terpin2024learning), written in JAX (jax2018github), and will be made publicly available. Since our approach builds on the JKO framework, in \wasyparagraph[5.1](https://arxiv.org/html/2506.01502v2#S5.SS1), we first reproduce the experimental setup from JKOnet^∗ for learning a known potential energy. Although we found the original codebase well-structured and accessible, we identified several inconsistencies in the data generation process, which we discuss in detail there. Finally, in \wasyparagraph[5.2](https://arxiv.org/html/2506.01502v2#S5.SS2), we compare our method against popular non-JKO baselines on single-cell RNA-seq datasets. Extended comparisons can be found in Appendix [B](https://arxiv.org/html/2506.01502v2#A2).

JKO-based Models. As demonstrated in Appendix [B.1.2](https://arxiv.org/html/2506.01502v2#A2.SS1.SSS2) and [B.1.3](https://arxiv.org/html/2506.01502v2#A2.SS1.SSS3), learning the interaction and internal energy components directly from samples is challenging – likely because their estimation requires computing integrals over functions that must themselves be estimated. Consequently, throughout the experiments in this section, we use only the potential energy parametrization, i.e., we impose the inductive bias $\theta_{2}=\theta_{3}=0$. Following the terminology of (terpin2024learning), we refer to this model as iJKOnetV. For a fair comparison, we also implemented a time-varying potential energy parametrization, as described in Appendix B of the original paper, and denote this variant by iJKOnett,V. However, as discussed in Appendix [A.2](https://arxiv.org/html/2506.01502v2#A1.SS2), this parametrization does not correspond to a single, consistent energy functional governing the dynamics. Instead, it actually solves a sequence of $K$ independent optimal transport problems between consecutive snapshots (bunne2023learning). Moreover, because JKOnet (bunne2022proximal) is computationally expensive, we limit our comparisons to experiments with moderate dimensionality, with the corresponding results reported in Appendix [B.2](https://arxiv.org/html/2506.01502v2#A2.SS2). Technical details of model training are provided in Appendix [C](https://arxiv.org/html/2506.01502v2#A3).

Metrics. We evaluate the next-step distribution $\hat{\rho}_{k+1}$ from $\rho_{k}$ and compare it to the ground-truth distribution $\rho_{k+1}$ using the following metrics: (a) Earth Mover’s Distance (EMD) (or $d_{{\mathbb{W}}_{1}}$) (rubner1998metric) and $d_{{\mathbb{W}}_{2}}$, (b) Bures-Wasserstein Unexplained Variance Percentage (B$d_{{\mathbb{W}}_{2}}^{2}$-UVP) (korotin2021continuous), (c) $\mathcal{L}_{2}$-based Unexplained Variance Percentage (${\mathcal{L}}_{2}$-UVP) (korotin2021wasserstein), which measures the discrepancy between the ground-truth functional $F^{\ast}$ and its reconstruction $\widehat{F}$, and (d) Maximum Mean Discrepancy (MMD) (gretton2012kernel). Formal definitions and implementation details for all metrics are provided in Appendix [C.1](https://arxiv.org/html/2506.01502v2#A3.SS1).

### 5.1 Learning potential energy

As discussed in \wasyparagraph[1](https://arxiv.org/html/2506.01502v2#S1), the population dynamics setting assumes that particle trajectories are not tracked across time – i.e., each time snapshot should consist of independently sampled particles without temporal correlation. This implies that data at each time step should be regenerated from scratch. However, we identified an inconsistency in the original codebase[1](https://arxiv.org/html/2506.01502v2#footnote1): particle trajectories are preserved across time steps, resulting in temporally correlated samples. We refer to this as the paired setup, where each particle $x_{k}$ is directly linked to $x_{k+1}$ along a trajectory. In contrast, the intended, temporally uncorrelated setting is referred to as the unpaired setup. Table [3](https://arxiv.org/html/2506.01502v2#A2.T3) in Appendix [B.1.1](https://arxiv.org/html/2506.01502v2#A2.SS1.SSS1) shows that switching between setups has a substantial impact on performance. A detailed discussion of this phenomenon is provided in Appendix [A.3](https://arxiv.org/html/2506.01502v2#A1.SS3).

Following (terpin2024learning, \wasyparagraph4.1), we begin by examining how our method learns potentials in the 2D paired setup, which allows for a direct visual comparison with $\texttt{JKOnet}^{\ast}$. These experiments are conducted on the synthetic dataset from (terpin2024learning, Appendix B), with results presented in Figure [3](https://arxiv.org/html/2506.01502v2#A2.F3) and further details provided in Appendix [B.1.1](https://arxiv.org/html/2506.01502v2#A2.SS1.SSS1). We then repeat the same procedure in the 2D unpaired setup, analyzing how performance scales with the number of samples (2K and 10K). The corresponding results are shown in Figure [2](https://arxiv.org/html/2506.01502v2#S5.F2). They indicate that our method outperforms $\texttt{JKOnet}^{\ast}$ on nearly all potentials, and that for certain cases, increasing the number of samples does not improve performance, highlighting the difficulty of the unpaired setup.

Refer to caption: x1.png

### 5.2 Learning population dynamics for single-cell data

Following (chen2023deep) and (shen2025multimarginal), in this section, we apply our method to modeling the dynamics of single-cell RNA data, with extended comparisons provided in Appendix [B.2](https://arxiv.org/html/2506.01502v2#A2.SS2).

Dataset. We use the Embryoid Body (EB) dataset (moon2019visualizing) and follow the preprocessing pipeline described in (tong2020trajectorynet). The EB dataset comprises a collection of 5 timepoints sampled over a period of 30 days.

Non-JKO Models. We evaluate our method against classic algorithms for trajectory inference such as TrajectoryNet (tong2020trajectorynet), MIOFLOW (huguet2022manifold), DMSB (chen2023deep), NLSB (koshizuka2023neural) and the recently proposed MMSB (shen2025multimarginal). For these methods, the reported metrics are taken from the papers referenced in each table caption.

**Table 1: 5D experiment. $d_{{\mathbb{W}}_{2}}$ distance ($\downarrow$) comparison across $t_{2}$ and $t_{4}$. Results for baseline methods: (shen2025multimarginal).**
| Method | $\mathbf{t_{2}}$ | $\mathbf{t_{4}}$ |
| --- | --- | --- |
| Vanilla-SB | $1.49\pm 0.063$ | $1.55\pm 0.034$ |
| DMSB | $1.13\pm 0.082$ | $1.45\pm 0.16$ |
| TrajectoryNet | $2.03\pm 0.04$ | $1.93\pm 0.08$ |
| MMSB | $1.27\pm 0.028$ | $1.57\pm 0.048$ |
| $\text{JKOnet}^{\ast}_{V}$ | $1.145\pm 0.033$ | $2.529\pm 0.014$ |
| $\text{iJKOnet}_{V}$ (Ours) | $\mathbf{1.082\pm 0.011}$ | $\mathbf{1.147\pm 0.001}$ |
| $\text{JKOnet}^{\ast}_{t,V}$ | $4.414\pm 1.499$ | $2.771\pm 0.197$ |
| $\text{iJKOnet}_{t,V}$ (Ours) | $\mathbf{0.983\pm 0.037}$ | $\mathbf{0.849\pm 0.021}$ |

Leave-two-out in 5D. Following (shen2025multimarginal), we perform experiments using a leave-two-out setup. Since the EB dataset contains five timesteps, we remove the second ($\mathbf{t_{2}}$) and fourth ($\mathbf{t_{4}}$) timesteps from training data, and then evaluate how well our method can reconstruct the data from the remaining $\mathbf{t_{1}}$ and $\mathbf{t_{3}}$ timesteps. The results are demonstrated in Table [2](https://arxiv.org/html/2506.01502v2#S5.T2). Vanilla-SB corresponds to the IPF (GP) method (vargas2021solving); further details are provided in Appendix B.6 of (shen2025multimarginal). Our method outperforms previous approaches in terms of reconstruction metrics.

Leave-one-out in 100D. Following (chen2023deep), we conduct experiments using a leave-one-out setup. One of the timesteps $\mathbf{t_{2}}$, $\mathbf{t_{3}}$, or $\mathbf{t_{4}}$ is omitted, and we evaluate the method’s ability to reconstruct the distribution of the left-out timestep. The results are shown in Table [1](https://arxiv.org/html/2506.01502v2#S5.T1). Our method attains competitive performance with DMSB in the all-time setting (w/o LO) while relying on a simpler, simulation-free optimization routine that requires no trajectory caching, unlike DMSB.

**Table 2: 100D experiment. MMD distance ($\downarrow$). Comparison of methods across different leave-one-out splits. The results are averaged across 3 runs. Results for baseline methods: (chen2023deep).**
| Method | $\text{LO-}\mathbf{t_{2}}$ | $\text{LO-}\mathbf{t_{3}}$ | $\text{LO-}\mathbf{t_{4}}$ | w/o LO |
| --- | --- | --- | --- | --- |
| NLSB (koshizuka2023neural) | $0.38$ | $0.37$ | $0.37$ | $0.66$ |
| MIOFLOW (huguet2022manifold) | $0.23$ | $0.90$ | $0.23$ | $0.23$ |
| DMSB (chen2023deep) | $\mathbf{0.042\pm 0.020}$ | $\mathbf{0.033\pm 0.003}$ | $\mathbf{0.040\pm 0.020}$ | $\mathbf{0.032\pm 0.003}$ |
| $\text{JKOnet}^{\ast}_{V}$ (terpin2024learning) | $0.220\pm 0.025$ | $0.293\pm 0.018$ | $0.235\pm 0.006$ | $0.229\pm 0.052$ |
| $\text{iJKOnet}_{V}$ (Ours) | $\mathit{0.137\pm 0.001}$ | $\mathit{0.123\pm 0.001}$ | $\mathit{0.097\pm 0.002}$ | $\mathit{0.085\pm 0.024}$ |
| $\text{JKOnet}^{\ast}_{t,V}$ (terpin2024learning) | $0.575\pm 0.119$ | $0.619\pm 0.157$ | $0.456\pm 0.056$ | $0.477\pm 0.098$ |
| $\text{iJKOnet}_{t,V}$ (Ours) | $0.848\pm 0.043$ | $0.370\pm 0.233$ | $0.055\pm 0.007$ | $0.124\pm 0.243$ |

## 6 Discussion

Contributions. We introduce the novel iJKOnet method, which recovers energy functionals that govern the evolution of a system. Our work bridges the gap between inverse optimization theory and computational gradient flows. The proposed method outperforms previous JKO-based approaches in all comparisons and achieves comparable results to non-JKO methods on real-world tasks. We also establish theoretical guarantees for recovering dynamics driven solely by potential energy.

Limitations. Our method does not recover other types of internal energy functionals, such as those corresponding to the porous-medium equation (alvarez2021optimizing), nor does it handle time-dependent interaction energies, although theoretical results for such cases exist (ferreira2018gradient). Furthermore, it does not account for birth–death dynamics as in recent trajectory inference methods (zhang2025learning). The reliance on entropy estimation may also lead to performance degradation in higher-dimensional settings, and we find that learning interaction energies remains particularly challenging. As a result, jointly optimizing all energy parameters ($\theta_{1},\theta_{2},\theta_{3}$) can cause instability and convergence to inaccurate potential estimates. Finally, from a theoretical standpoint, our analysis currently provides guarantees only for potential energy (Theorem ([3.1](https://arxiv.org/html/2506.01502v2#S3.Thmtheorem1))); extending it to other energy types is an avenue for future work.

LLM Usage. Large Language Models (LLMs) were used exclusively to assist with sentence rephrasing and improving text clarity. All scientific content, results, and interpretations in this paper were produced solely by the authors.

## Appendix A Extended Discussion

We first review related work in \wasyparagraph[A.1](https://arxiv.org/html/2506.01502v2#A1.SS1) and then discuss the validity of time-varying potentials proposed by terpin2024learning in \wasyparagraph[A.2](https://arxiv.org/html/2506.01502v2#A1.SS2) in contrast to the stationary setting considered in this paper. Finally, in \wasyparagraph[A.3](https://arxiv.org/html/2506.01502v2#A1.SS3), we discuss potential reasons for the performance gap between the paired and unpaired setups introduced in \wasyparagraph[5.1](https://arxiv.org/html/2506.01502v2#S5.SS1).

### A.1 Related works

In this section, we first provide an extended discussion of related work on learning population dynamics (also known as trajectory inference) in \wasyparagraph[A.1.1](https://arxiv.org/html/2506.01502v2#A1.SS1.SSS1), followed by an overview of applications of the JKO scheme (jordan1998variational) in deep learning in \wasyparagraph[A.1.2](https://arxiv.org/html/2506.01502v2#A1.SS1.SSS2).

#### A.1.1 Learning population dynamics

In contrast to our approach, which models population dynamics as a Wasserstein gradient flow (WGF) (ambrosio2008gradient) of the free-energy functional ([5](https://arxiv.org/html/2506.01502v2#S2.E5)) (gomez2024beginner), several other fruitful approaches have been proposed in the literature. We review these alternatives below.

RNNs and Neural ODEs. One of the first works on learning population dynamics is (hashimoto2016learning), which used recurrent neural networks (RNNs) to learn the SDE ([7](https://arxiv.org/html/2506.01502v2#S2.E7)) as a WGF of a potential energy functional. Later, (chen2018neural, CNF) introduced neural ODEs as the continuous-time limit of RNNs for learning ODE-based dynamics, and (li2020scalable) subsequently extended this approach to handle SDEs. Building on these ideas, (erbe2023transcriptomic, RNAForecaster) applied neural ODEs to predict future transcriptomic states of single cells (battich2020sequencing).

Static OT. Another line of work models population dynamics by directly learning $K$ independent transport maps from time-snapshot data, using the static optimal transport formulation ([1](https://arxiv.org/html/2506.01502v2#S2.E1)) or ([2](https://arxiv.org/html/2506.01502v2#S2.E2)), without attempting to recover the unifying energy functionals guiding the system evolution. For example, (schiebinger2019optimal, Waddington-OT) employs discrete transport maps, whereas (bunne2023learning, CellOT) learns neural transport maps, as in (fan2023neural, Eq.(6)), (rout2022generative, Eq.(14)), and (korotin2023neural, Eq.(15)), see Appendix [A.2](https://arxiv.org/html/2506.01502v2#A1.SS2) for more details.

Dynamic OT. (tong2020trajectorynet, TrajectoryNet) models population dynamics using the dynamic optimal transport formulation of (benamou2000computational), parameterized by a regularized CNF. The work (huguet2022manifold, MIOFlow) later extended this approach by incorporating a geodesic autoencoder to better capture manifold structure. More recently, (wan2023scalable) proposed computational techniques that make dynamic OT scalable to higher-dimensional settings.

Flow Matching. The works in this paragraph extend Flow Matching (FM) techniques (lipman2023flow; albergo2023building; liu2023flow) in several directions. Conditional Flow Matching (tong2024improving, CFM) introduce simulation-free objectives for learning deterministic flows. Building on the celebrated work of (otto2001geometry), the Wasserstein space can be nominally (i.e., heuristically) viewed as a Riemannian manifold. Motivated by this observation, Wasserstein Flow Matching (haviv2025wasserstein, WFM) applies the idea of Riemannian Flow Matching (chen2024flow, RFM) to the Wasserstein space, i.e., performing FM over a distribution of distributions. Further developments, such as Meta Flow Matching (atanackovic2025meta, Meta FM), embed the population of samples with a Graph Neural Network (GNN) (liu2025graph) and use these embeddings as conditioning inputs for the learning vector field. Finally, Multi-Marginal Flow Matching (rohbeck2025modeling, MMFM) constructs a flow using smooth spline-based interpolation across time points and conditions, and regresses it with a neural network under the classifier-free guided Flow Matching framework (zheng2023guided). For a broader overview of FM for learning population dynamics, see the recent surveys (morehead2025go; li2025flow).

Action Matching. (neklyudov2023action, AM) introduced Action Matching (AM), a method that optimizes an action-gap objective. In contrast to Flow Matching, AM learns a vector field expressed as the gradient of an action, $\nabla s^{\ast}_{t}$, which uniquely defines the velocity field that transports particles optimally in the sense of optimal transport.
Building on this idea, (berman2024parametric) conditioned the action $s_{t,\mu}$ on physical parameters $\mu$, used (berman24bcolora, CoLoRa) for efficient action parametrization, and demonstrated its utility for surrogate modeling of classical numerical solvers, enabling fast prediction of system behavior across different physics parameter settings. Most recently, (neklyudov2024a, WLF) proposed Wasserstein Lagrangian Flows, a unifying framework that minimizes Lagrangian action functionals over the space of probability densities rather than the ground space, thereby encapsulating AM as a special case.

Schrödinger bridges. Two seminal works, (vargas2021solving, IPF (GP)) and (de2021diffusion, IPF (NN)), proposed solving the Schrödinger bridge problem using the Iterative Proportional Fitting (IPF) algorithm (gramer2000probability) in application to generative modeling, which alternates between forward and backward processes. The first approach employed Gaussian processes as parametrization, whereas the second relied on score-based neural networks. Later, bridge matching methods was developed (peluchetti2023non; peluchetti2023diffusion; liu2023learning), (shi2023diffusion, DSBM). tong2024simulation introduced simulation-free objectives for both deterministic and stochastic flows, further reducing computational overhead. Additional contributions include generalized setups for Schrödinger bridges (koshizuka2023neural, NLSB) and (liu2024generalized, GSBM), multi-marginal generalizations of the problem (chen2023deep, DMSB), (shen2025multimarginal, MMSB), (hong2025trajectory), and momentum-accelerated formulations (theodoropoulos2025momentum, 3MSBM).

Unbalanced OT. Several works jointly model marginal transitions and growth dynamics by minimizing the action in the Wasserstein-Fisher-Rao (WFR) metric (chizat2018interpolating), i.e., by solving the dynamical unbalanced optimal transport problem (chizat2018unbalanced). (tong2023learning) proposed BEMIOflow, an extension of (huguet2022manifold, MIOflow), that incorporates a neural network to predict cell growth and death rates continuously in time, thereby augmenting optimal transport with population-size dynamics. Later, (chen2022most) formulated an unbalanced Schrödinger Bridge problem, which was applied in practice by (zhang2025learning, DeepRUOT), following the developments in (sha2024reconstructing, TIGON). DeepRUOT was further simplified into a single-network formulation, Var-RUOT (sun2025variational), by exploiting optimality conditions; extended to include interaction modeling in (zhang2025modeling, CytoBridge); and adapted to use two independent networks to separately model distributional drift and mass growth (wang2025joint, VGFM), (zhang2025inferring). In a different direction, (klein2024genot, GENOT) generalized OT to simultaneously handle stochasticity, entropy regularization, and unbalanced mass transport, enabling applications to cross-modal and heterogeneous data. For a broader survey of trajectory inference methods and their applications to biology data, see (zhang2025deciphering; zhang2025integrating).

The concurrent work (andrade2025learning, iJKO) also connects inverse optimization with the JKO scheme ([8](https://arxiv.org/html/2506.01502v2#S2.E8)) but in the unbalanced setting. However, their focus is on sample complexity, the method addresses only potential energy, and experiments are restricted to low-dimensional cases.

Riemannian perspective. (scarvelis2023riemannian) proposed learning a metric tensor $A(x)$ that minimizes the average 1-Wasserstein distance on the learned manifold between pairs of consecutive population snapshots. Furthermore, (kapusniak2024metric, MFM) introduced Metric Flow Matching, which learns interpolants
$x_{t,\eta}=(1-t)x_{0}+tx_{1}+t(1-t)\varphi_{t,\eta}(x_{0},x_{1})$,
where $\eta$ are the parameters of a neural network $\varphi_{t,\eta}$ providing a nonlinear “correction” to straight-line interpolants (albergo2023building). The resulting velocity field minimizes a data-dependent Riemannian metric, assigning lower transport cost to regions with higher data density.

Splines in Wasserstein Space. In the seminal work of (schiebinger2019optimal), a piecewise linear OT interpolation method was proposed to infer cell trajectories. Subsequent works introduced higher-order piecewise polynomials (e.g., cubic splines) in Wasserstein space: (chen2018measure; benamou2019second) formulated a global cubic spline minimization problem, which was later extended by (chewi2021fast) to use Euclidean interpolation algorithm in $\mathbb{R}^{D}$ after a finding optimal Monge map ([2](https://arxiv.org/html/2506.01502v2#S2.E2)) between samples from consecutive measures, with further refinements by (clancy2022wasserstein; justiniano2024approximation). More recently, (banerjee2025efficient, WLR) proposed the Wasserstein Lane–Riesenfeld method, leveraging the classical subdivision algorithm of (lane1980theoretical). In parallel, (dyn2025subdivision) developed subdivision schemes for general metric spaces, including the Wasserstein space. (baccou2024subdivision; kawano2025ct) combines subdivision schemes with OT.

#### A.1.2 JKO scheme meets deep learning

Generative Modeling. (vidal2023taming, JKO-Flow) and (xu2023normalizing, JKO-iFlow) reinterpret the JKO scheme ([8](https://arxiv.org/html/2506.01502v2#S2.E8)) through the framework of CNFs (chen2018measure), applying it to 2D synthetic data and image generation tasks, respectively. Subsequently, (cheng2024convergence) provided theoretical convergence guarantees for this approach. (choi2024scalable, S-JKO) further accelerated WGF-based generative modeling by leveraging a semi-dual unbalanced OT formulation. They constructed the WGF of an $f$-divergence $D_{f}$ and used an (song2021scorebased, NCSN++) backbone – similar to (zhu2024neural, NSGF), who designed a generative model as the WGF with respect to the Sinkhorn divergence (peyre2019computational).

Variational Inference. (lambert2022variational) and (diao2023forward) study variational inference through the lens of Bures–Wasserstein gradient flows. (cheng2023particle, GWG) introduce a generalized minimizing movement scheme on the space ${\mathcal{P}}_{c_{h}}(\mathbb{R}^{D})$, where the transport cost is defined as $c_{h}(x,y)=g\!\left(\frac{x-y}{h}\right)h$ with $g$ a continuously differentiable Young function, and apply this framework to particle-based variational inference.

Datasets/Weights Learning. (bonet2025flowing) propose an explicit scheme that is computationally more efficient in practice than the implicit JKO scheme on the space ${\mathcal{P}}_{2}({\mathcal{P}}_{2}(\mathbb{R}^{D}))$, and apply this approach to modeling flow datasets viewed as random measures. (saragih2025flows) train a meta-model, based on JKO as well as other flow- and diffusion-based approaches, that generates dataset-conditioned classifiers by producing the neural network weights in a single forward pass.

### A.2 Time-varying potentials discussion

In this section, we examine the validity of using time-varying potentials, as proposed in (terpin2024learning, \wasyparagraph4.4). We demonstrate that, within this framework, the problem of reconstructing energy functionals (bunne2022proximal) through the JKO time-discretization scheme reduces to learning $K$ independent neural optimal transport maps (fan2023neural, Eq. (6)); (rout2022generative, Eq. (14)) that connect the data snapshots $\rho_{k}$ and $\rho_{k+1}$. To illustrate this, consider the loss in ([11](https://arxiv.org/html/2506.01502v2#S3.E11)) and assume that, instead of a single functional ${\mathcal{J}}$, we now assign a separate only potential energy functional ${\mathcal{J}}^{k}_{\text{PE}}=\int_{\mathcal{X}}V^{k}(x)\,\differential\rho_{k}(x)$ for each time step $k$:

$$ $\!\!\max_{{\mathcal{J}}^{k}_{\text{PE}}}\min_{T^{k}}{\mathcal{L}}({\mathcal{J}}^{k}_{\text{PE}},T^{k})\!\stackrel{{\scriptstyle\text{def}}}{{=}}\!\max_{{\mathcal{J}}^{k}_{\text{PE}}}\min_{T^{k}}\!\sum_{k=0}^{K-1}\Big[{\mathcal{J}}^{k}_{\text{PE}}(T^{k}\sharp\rho_{k})-{\mathcal{J}}^{k}_{\text{PE}}(\rho_{k+1})+\frac{1}{2\tau}\int_{\mathcal{X}}\!\|x-T^{k}(x)\|^{2}_{2}\rho_{k}(x)\differential x\Big].$ $$

It then becomes clear that the loss optimization decomposes into $K$ independent terms ${\mathcal{L}}^{k}({\mathcal{J}}^{k}_{\text{PE}},T^{k})$:

$$ ${\mathcal{L}}^{k}({\mathcal{J}}^{k}_{\text{PE}},T^{k})\!\stackrel{{\scriptstyle\text{def}}}{{=}}\!{\mathcal{J}}^{k}_{\text{PE}}(T^{k}\sharp\rho_{k})-{\mathcal{J}}^{k}_{\text{PE}}(\rho_{k+1})+\frac{1}{2\tau}\int_{\mathcal{X}}\|x-T^{k}(x)\|^{2}_{2}\rho_{k}(x)\,\differential x\to\max_{{\mathcal{J}}^{k}_{\text{PE}}}\min_{T^{k}}.$ (18) $$

By replacing $V^{k}$ with the dual potential $f_{\eta}$ and $T^{k}$ with $T_{\theta}$, we recover Equation (6) from (fan2023neural), which aligns with the objective proposed in (bunne2023learning, Eq. (9)) for modeling single-cell dynamics. Thus, (terpin2024learning) introduces a new strategy for minimizing this loss rather than directly addressing the recovery of ground-truth energy functionals (bunne2022proximal).

Fortunately, the theoretical framework developed by ferreira2018gradient and later extended by plazotta2016high establishes the validity of this approach even for time-varying free energy functionals ([5](https://arxiv.org/html/2506.01502v2#S2.E5)), not only for the potential energy case, in the limit as the time discretization tends to zero. However, this connection is not acknowledged in (terpin2024learning).

### A.3 Paired vs. unpaired setups

In this section, we discuss the substantial performance gap observed between the paired and unpaired setups, as illustrated in Table [3](https://arxiv.org/html/2506.01502v2#A2.T3). We suggest that this discrepancy stems from the intrinsic connection between the considered solvers, iJKOnet and $\texttt{JKOnet}^{\ast}$, and the underlying OT problem ([2](https://arxiv.org/html/2506.01502v2#S2.E2)).

Both solvers are designed to recover the optimal transport mappings $T^{k,\ast}$ between consecutive distributions $\rho_{k}$ and $\rho_{k+1}$, satisfying the pushforward relation $T^{k,\ast}\sharp\rho_{k}=\rho_{k+1}$. In $\texttt{JKOnet}^{\ast}$, these mappings are obtained explicitly by precomputing the OT maps between $\rho_{k}$ and $\rho_{k+1}$, whereas in our method, they are learned implicitly via neural networks $T^{k}_{\varphi}$.

In practice, the continuous distributions $\rho_{k}$ are replaced by their empirical counterparts $\hat{\rho}_{k}$, constructed from training samples $\{x_{k}^{1},x_{k}^{2},\dots,x_{k}^{N}\}=X_{k}\sim\rho_{k}$. This is precisely where the distinction between the paired and unpaired setups becomes critical:

- •
Paired. When samples are paired, they follow the ground-truth transport mappings, i.e.,
$T^{0,\ast}(x_{0}^{i})=x_{1}^{i},\quad T^{1,\ast}(x_{1}^{i})=x_{2}^{i},\quad\dots,\quad T^{k,\ast}(x_{k}^{i})=x_{k+1}^{i},\quad\dots$
In this case, estimating $T^{k,\ast}$ becomes substantially simpler and can be viewed as a supervised regression problem over the given sample pairs.
- •
Unpaired. In contrast, when the samples $X_{k}$, $k\in\{0,\dots,K\}$, are mutually independent, the quality of the estimated OT maps degenerate significantly. This degradation stems from the high sample complexity of Wasserstein-2 distances and maps (known to be relatively poor, see hutter2021minimax), as well as from the increased variance of the estimated transport mappings.

A more in-depth theoretical investigation of this phenomenon remains an interesting direction for future research.

## Appendix B Extended Comparisons

In this section, we first present extended synthetic experiments in \wasyparagraph[B.1](https://arxiv.org/html/2506.01502v2#A2.SS1), followed by results on real-world single-cell data in \wasyparagraph[B.2](https://arxiv.org/html/2506.01502v2#A2.SS2).

### B.1 Synthetic comparisons

In this section, we further investigate our method’s ability to learn potential energy (\wasyparagraph[B.1.1](https://arxiv.org/html/2506.01502v2#A2.SS1.SSS1)), interaction energy (\wasyparagraph[B.1.2](https://arxiv.org/html/2506.01502v2#A2.SS1.SSS2)), and internal energy (\wasyparagraph[B.1.3](https://arxiv.org/html/2506.01502v2#A2.SS1.SSS3)).

Since our approach builds on the JKO framework for modeling population dynamics, we adopt the experimental setup from JKOnet^∗ and compare our method against both JKOnet bunne2022proximal and JKOnet^∗ terpin2024learning on corresponding synthetic evaluation tasks.

Figure: Figure 3: Level curves of the true (green) and estimated (blue) potentials for the paired setup, following (terpin2024learning, Appendix F). These results can be directly compared with those in (terpin2024learning, Figure 6). Note that for the flat potential, the value range is near zero, as expected for the ground-truth potential.
Refer to caption: x5.png

#### B.1.1 Learning potential energy

Paired Setup. Following terpin2024learning, we evaluate our method in the paired setup (see details in \wasyparagraph[5.1](https://arxiv.org/html/2506.01502v2#S5.SS1)), which enables a direct visual comparison with $\texttt{JKOnet}^{\ast}$ (terpin2024learning). The experiments are conducted on the synthetic dataset from terpin2024learning, using a step size of $\tau=0.01$, $K=5$ time steps, and $N=400$ samples per step with 40% of test samples. Figure [3](https://arxiv.org/html/2506.01502v2#A2.F3) presents the ground-truth potentials $V(x)$ (green), defined in terpin2024learning, alongside the reconstructed potentials $V_{\theta}(x)$ (blue) learned by our method. As shown, our approach achieves performance comparable to $\texttt{JKOnet}^{\ast}$ in the paired setup.

Unpaired Setup. We repeat the experiments on a corrected version of the synthetic dataset, where samples are not temporally correlated. We select the six most challenging potentials in terms of convergence. Table [3](https://arxiv.org/html/2506.01502v2#A2.T3) reports the ratio of the final metrics between the unpaired and paired setups. As shown, switching to the unpaired setup significantly affects the method’s performance.

**Table 3: Comparison of paired and unpaired setups. Each value indicates the ratio of the final metric in the unpaired setup to that in the paired setup. Most ratios significantly exceed 1, highlighting the increased difficulty of the unpaired setting.**
|  |  | $\texttt{iJKOnet}_{V}$ (Ours) | $\texttt{JKOnet}_{V}^{\ast}$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Alias | Potential Name | B$d_{{\mathbb{W}}_{2}}^{2}$-UVP | EMD | $\mathcal{L}_{2}$-UVP | B$d_{{\mathbb{W}}_{2}}^{2}$-UVP | EMD | $\mathcal{L}_{2}$-UVP |
| FL | Flowers | 8066 | 90 | 297 | 135325 | 305 | 5809 |
| FN | Friedman | 556 | 5 | 2 | 263 | 6 | 3 |
| IG | Ishigami | 2259 | 38 | 36 | 4276 | 64 | 210 |
| WS | Watershed | 18135 | 220 | 1509 | 512720 | 1014 | 102592 |
| WP | Wavy plateau | 2 | 2 | 7 | 1 | 2 | 29 |
| ZR | Zigzag ridge | 158 | 23 | 13 | 155 | 25 | 31 |

#### B.1.2 Learning interaction energy

In this section, we evaluate how accurately our method and $\texttt{JKOnet}^{\ast}$ (terpin2024learning) can approximate the interaction energy in the unpaired setup. We exclude JKOnet from comparison, as its design does not support learning interaction energies.

We identified an inconsistency in the theoretical formulation presented in the $\texttt{JKOnet}^{\ast}$ paper (terpin2024learning). Specifically, as discussed in \wasyparagraph[2](https://arxiv.org/html/2506.01502v2#S2) (see also (santambrogio2015optimal, \wasyparagraph7.2)), the ground-truth interaction functional must be symmetric; that is, in ([5](https://arxiv.org/html/2506.01502v2#S2.E5)), the interaction kernel should satisfy $W(z)=W(-z)$. Incorporating this correction, we conducted an experiment to assess the ability of both methods to recover the interaction energy in a 2D unpaired setting.

Experimental Setup. The ground-truth interaction kernels are defined as $W(z)=\frac{1}{2}\big(W_{b}(z)+W_{b}(-z)\big)$, where $W_{b}$ denotes the base functionals listed in Table [3](https://arxiv.org/html/2506.01502v2#A2.T3). Both methods were restricted to use only the interaction energy component ${\mathcal{W}}_{\theta_{2}}$ in the energy parameterization ([12](https://arxiv.org/html/2506.01502v2#S3.E12)). The corresponding variants are denoted as $\texttt{iJKOnet}_{W}$ and $\texttt{JKOnet}^{\ast}_{W}$. Similar to \wasyparagraph[B.1.1](https://arxiv.org/html/2506.01502v2#A2.SS1.SSS1), we use a step size of $\tau=0.01$, $K=5$ time steps, and two sample size settings: a small-scale setup with $2\text{K}$ total samples ($N=400$ per step, with 40% reserved for testing) and a large-scale setup with $10\text{K}$ total samples ($N=2000$ per step, with 40% reserved for testing). Qualitative results are shown in Figure [4](https://arxiv.org/html/2506.01502v2#A2.F4), and quantitative ${\mathcal{L}}_{2}$-UVP metrics (refer to \wasyparagraph[C.1](https://arxiv.org/html/2506.01502v2#A3.SS1) for details) are reported in Table [4](https://arxiv.org/html/2506.01502v2#A2.T4).

Results. Neither method is able to recover the interaction energy accurately, likely due to biases introduced by the batched estimation of the interaction term ${\mathcal{W}}_{\theta_{2}}$ in ([12](https://arxiv.org/html/2506.01502v2#S3.E12)), which involves squaring the input. This approximation appears to cause both models to converge toward batch-specific minima.

**Table 4: Quantitative ${\mathcal{L}}_{2}$-UVP results for the unpaired setup in 2D interaction energy learning. Increasing the number of samples improves performance, though the effect is relatively limited.**
|  | Sample size: 10K | Sample size: 2K |  |  |
| --- | --- | --- | --- | --- |
| Potential | $\texttt{iJKOnet}_{W}$ (Ours) | $\texttt{JKOnet}^{\ast}_{W}$ | $\texttt{iJKOnet}_{W}$ (Ours) | $\texttt{JKOnet}^{\ast}_{W}$ |
| Flowers | 0.0009 | 0.0069 | 0.0032 | 0.0101 |
| Friedman | 0.0024 | 0.0093 | 0.0087 | 0.0010 |
| Ishigami | 0.0008 | 0.0066 | 0.0046 | 0.0087 |
| Watershed | 0.0009 | 0.0081 | 0.0043 | 0.0086 |
| Zigzag Ridge | 0.0017 | 0.0034 | 0.0060 | 0.0104 |

Figure: (a) (a) Ground-truth ${\mathcal{W}}^{\ast}$
Refer to caption: images/exp4/flowers_gt_level_curves_potential.png

#### B.1.3 Learning internal energy

In this section, we assess how accurately our method and $\texttt{JKOnet}^{\ast}$ (terpin2024learning) can estimate the diffusion coefficient $\theta_{3}$ in ([12](https://arxiv.org/html/2506.01502v2#S3.E12)) corresponding to the internal energy term under the unpaired setup. We exclude JKOnet from comparison, as its design does not support learning internal energy.

Experimental Setup. Following the protocol of terpin2024learning, we estimate the diffusion coefficient $\theta_{3}$ in 2D and 20D settings using both iJKOnet and $\texttt{JKOnet}^{\ast}$, i.e., the full parameterization in ([12](https://arxiv.org/html/2506.01502v2#S3.E12)). The same functional from Appendix [C.7](https://arxiv.org/html/2506.01502v2#A3.SS7) is used for both potential and interaction energies. The ground-truth diffusion levels are set to $\beta^{\ast}\in\{0.0,0.1,0.2\}$.

Results. Figure [5](https://arxiv.org/html/2506.01502v2#A2.F5) shows that iJKOnet fails to recover the ground-truth $\beta^{\ast}$ values, while $\texttt{JKOnet}^{\ast}$ provides a closer approximation. However, all predicted $\theta_{3}$ values tend to converge toward 0, deviating from the true levels $\beta^{\ast}\in\{0.0,0.1,0.2\}$. This indicates that accurately learning the internal energy remains challenging for both models.

Figure: Figure 5: Estimation error of the diffusion coefficient $\theta_{3}$ relative to the ground-truth values $\beta^{\ast}$. Blue bars correspond to iJKOnet and orange bars to $\texttt{JKOnet}^{\ast}$. The Y-axis indicates the absolute deviation between the estimated $\theta_{3}$ and true $\beta^{\ast}$ values.
Refer to caption: images/beta.png

### B.2 Learning single-cell dynamics

In this section, we present extended comparisons on real-world single-cell datasets, starting with 5D results and then considering 50D and 100D.

Datasets. Following \wasyparagraph[5.2](https://arxiv.org/html/2506.01502v2#S5.SS2), we consider the Embryoid Body (EB) dataset (moon2019visualizing) and additionally the Multiome (Multi) dataset (burkhardt2022multimodal). For both datasets, we apply the preprocessing pipeline of (tong2020trajectorynet). The Multi dataset contains single-cell measurements across four time points (days 2, 3, 4, and 7), which we use for high-dimensional experiments.

5D. Following (neklyudov2024computational, Table 2) and (terpin2024learning, \wasyparagraph4.4), we conducted experiments on the EB dataset in 5D. At each time step, we train for 1000 epochs on 60% of the data and compute the $d_{{\mathbb{W}}_{1}}$ distance between the observed distribution $\rho_{k}$ (remaining 40% of the data) and the one-step-ahead prediction $\hat{\rho}_{k}$. The results are shown in Table [5](https://arxiv.org/html/2506.01502v2#A2.T5). Our method outperforms several previous approaches when using a time-varying potential parametrization.

**Table 5: The results for the EB dataset in 5D are obtained by training on all time steps, evaluated for next-time-step prediction using the $d_{{\mathbb{W}}_{1}}$ metric. Results for non-JKO methods are taken from (neklyudov2024computational, Table 2):**
| Model | $t_{1}$ | $t_{2}$ | $t_{3}$ | $t_{4}$ | Mean |
| --- | --- | --- | --- | --- | --- |
| Neural SDE (li2020scalable) | 0.69 | 0.91 | 0.85 | 0.81 | 0.82 |
| TrajectoryNet (tong2020trajectorynet) | 0.73 | 1.06 | 0.90 | 1.01 | 0.93 |
| SB-FBSDE (chen2022likelihood) | 0.56 | 0.80 | 1.00 | 1.00 | 0.84 |
| NLSB (koshizuka2023neural) | 0.68 | 0.84 | 0.81 | 0.79 | 0.78 |
| OT-CFM (tong2024improving) | 0.78 | 0.76 | 0.77 | 0.75 | 0.77 |
| WLF-OT (neklyudov2024computational) | 0.65 | 0.78 | 0.76 | 0.75 | 0.74 |
| WLF-SB (neklyudov2024computational) | 0.63 | 0.79 | 0.77 | 0.74 | 0.73 |
| JKOnet (bunne2022proximal) | 1.53 | 1.27 | 1.13 | 1.41 | 1.34 |
| JKOnet${}^{\ast}_{V}$ terpin2024learning | 0.99 | 1.11 | 1.06 | 1.30 | 1.12 |
| iJKOnetV (Ours) | 0.92 | 1.11 | 0.95 | 1.21 | 1.05 |
| JKOnet${}^{\ast}_{t,V}$ (terpin2024learning) | 0.69 | 0.77 | 0.69 | 0.78 | 0.73 |
| iJKOnett,V (Ours) | 0.51 | 0.58 | 0.57 | 0.64 | 0.58 |

50D and 100D data. Following neklyudov2024computational, we conducted experiments on the Multi dataset in 50D and 100D using a leave-one-out setup, averaged over three runs. Models were trained on marginals from timepoint partitions [2,4,7] and [2,3,7], and evaluated on the corresponding left-out marginals at timepoints [3] and [4]. Training was performed for 5000 epochs.

The results are presented in Table [6](https://arxiv.org/html/2506.01502v2#A2.T6). Our method achieves comparable performance for the left-out [3] marginal. However, performance on the [4] marginal is poor, leading to a lower overall average. This is likely due to the longer time intervals between the learned marginals, which violate our assumption that consecutive marginals are obtained from JKO steps (see \wasyparagraph[3.1](https://arxiv.org/html/2506.01502v2#S3.SS1)), thereby reducing the accuracy of predicting the [4] marginal.

**Table 6: The results for Multi dataset for direct comparison with from (neklyudov2024computational, Table 1), averaged for 3 runs for leave-one-out setup, $\mathcal{W}_{2}$ metric.**
| Dimension |  | 50 |  |  | 100 |  |
| --- | --- | --- | --- | --- | --- | --- |
| Metric | BW-UVP | EMD | MMD | BW-UVP | EMD | MMD |
| $\text{JKOnet}_{V}^{\ast}$ | 128.729 $\pm$ 78.307 | 68.406 $\pm$ 6.055 | 0.0003 | 92.512 $\pm$ 38.075 | 72.639 $\pm$ 3.161 | 0.0002 |
| $\text{iJKOnet}_{V}$ (Ours) | 38.318 $\pm$ 18.325 | 50.560 $\pm$ 7.187 | 0.0003 | 38.991 $\pm$ 16.918 | 59.216 $\pm$ 6.869 | 0.0002 |
| $\text{JKOnet}_{t,V}^{\ast}$ | 89.264 $\pm$ 38.372 | 68.600 $\pm$ 9.846 | 0.0003 | 93.429 $\pm$ 34.332 | 78.674 $\pm$ 8.236 | 0.0002 |
| $\text{JKOnet}_{t,V}$ (Ours) | 36.156 $\pm$ 16.515 | 50.026 $\pm$ 7.727 | 0.0003 | 38.387 $\pm$ 18.525 | 59.318 $\pm$ 7.365 | 0.0002 |

## Appendix C Training and Model Details

### C.1 Metric computation details

In this section, we provide details on the evaluation metrics used in our experiments.

${\mathcal{L}}_{2}$-UVP. When the ground-truth functional $F^{\ast}$ (e.g., $V$ or $W$) is available, we assess the discrepancy between it and its reconstruction $\widehat{F}$ using the backward $\mathcal{L}_{2}$-based Unexplained Variance Percentage (${\mathcal{L}}_{2}$-UVP) metric introduced by korotin2021wasserstein, defined as follows:

$$ $\mathcal{L}_{2}\text{-UVP}(F^{\ast},\widehat{F})=100\cdot\frac{\tau^{2}\|\nabla\widehat{F}-\nabla F^{\ast}\|_{\rho_{k+1}}^{2}}{\mathrm{Var}(\rho_{k})}\%,$ (19) $$

where the norm $\|\cdot\|_{\rho_{k+1}}$ is computed with respect to the ground-truth distribution $\rho_{k+1}$, and $\tau$ denotes the time step size. Values close to 0% indicate that $\nabla\widehat{F}$ closely approximates $\nabla F^{\ast}$.

B$d_{{\mathbb{W}}_{2}}^{2}$-UVP. The Bures-Wasserstein UVP introduced by (korotin2021continuous) is defined as

$$ $\mathrm{B}d_{{\mathbb{W}}_{2}}^{2}\text{-UVP}(\rho_{k},\hat{\rho}_{k})\stackrel{{\scriptstyle\text{def}}}{{=}}100\cdot\frac{\mathrm{B}d_{{\mathbb{W}}_{2}}^{2}(\rho_{k},\hat{\rho}_{k})}{\frac{1}{2}\mathrm{Var}(\rho_{k})}\%,$ (20) $$

where the Bures-Wasserstein distance is given by

$$ $\mathrm{B}d_{{\mathbb{W}}_{2}}^{2}(\mathbb{P},\mathbb{Q})\stackrel{{\scriptstyle\text{def}}}{{=}}d_{{\mathbb{W}}_{2}}^{2}\left(\mathcal{N}(\mu_{\mathbb{P}},\Sigma_{\mathbb{P}}),\mathcal{N}(\mu_{\mathbb{Q}},\Sigma_{\mathbb{Q}})\right),$ (21) $$

with $\mu_{\mathbb{P}}$ and $\Sigma_{\mathbb{P}}$ denoting the mean and covariance of distribution $\mathbb{P}$, respectively.

EMD or $d_{{\mathbb{W}}_{1}}$. The Earth Mover’s Distance (EMD) (rubner1998metric) is defined as

$$ $\mathrm{EMD}(\rho_{k},\hat{\rho}_{k})\equiv d_{{\mathbb{W}}_{1}}(\rho_{k},\hat{\rho}_{k})\stackrel{{\scriptstyle\text{def}}}{{=}}\min_{\pi\in\Pi(\rho_{k},\hat{\rho}_{k})}\int_{{\mathcal{X}}\times{\mathcal{X}}}\|x-y\|\,\differential\pi(x,y),$ (22) $$

where $\rho_{k}$ and $\hat{\rho}_{k}$ denote the ground truth and predicted distributions at time step $k$.

MMD. Maximum Mean Discrepancy (MMD) (gretton2012kernel) between two distributions $\rho_{k}$ and $\hat{\rho}_{k}$ is defined as

$$ $\mathrm{MMD}^{2}(\rho_{k},\hat{\rho}_{k})=\left\|\mathbb{E}_{x\sim\rho_{k}}[h(x)]-\mathbb{E}_{\hat{x}\sim\hat{\rho}_{k}}[h(\hat{x})]\right\|_{\mathcal{H}}^{2},$ (23) $$

where $h$ is the feature map associated with the reproducing kernel $\mathbf{H}$, and $\|\cdot\|_{\mathcal{H}}$ denotes the norm in the corresponding RKHS (ghojogh2021reproducing). Here, $\rho_{k}$ and $\hat{\rho}_{k}$ are defined as in the EMD metric.

### C.2 Training details

Energy. To ensure stability, we accumulate gradients across all time steps $k=0,\dots,K$, using mini-batch estimates $\widehat{{\mathcal{J}}}_{\theta}$ of the energy function ${\mathcal{J}}_{\theta}$ (see Eq. ([12](https://arxiv.org/html/2506.01502v2#S3.E12))). For each time step $t_{k}$, we sample a mini-batch $X_{k}=\{x^{1}_{k},\dots,x^{B}_{k}\}\sim\rho_{k}$, with a fixed batch size $B$. The energy is then estimated as:

$$ ${\mathcal{J}}_{\theta}(\rho_{k})\approx\widehat{{\mathcal{J}}}_{\theta}(X_{k})=\frac{1}{B}\sum_{i=0}^{B}\left[V_{\theta_{1}}(x_{k})+\frac{1}{B}\sum_{j=0}^{B}W_{\theta_{2}}(x_{i}-x_{j})\right]-\theta_{3}\widehat{{\mathcal{H}}}(X_{k}),$ (24) $$

where $\widehat{{\mathcal{H}}}(X_{k})$ denotes the estimated entropy of $\rho_{k}$ discussed in the following section.

Map. In practice, we use two strategies to parameterize each transport map $T^{k}_{\varphi}$: for large-scale tasks, we assign a separate network to each time step; for moderate-dimensional tasks, we encode the time index $k$ as an additional input, i.e., $T_{\varphi}^{k}(x)=T_{\varphi}(x,k)$, which helps prevent overfitting. This approach performs better than using a network without any time embedding.

Training. We perform several optimization steps for the transport maps $T^{k}_{\varphi}$, parameterized by $\varphi$, while updating the energy function ${\mathcal{J}}_{\theta}$, parameterized by $\theta$, only once. This approach follows standard practice in min-max optimization. We experimented with common GAN stabilization methods (gulrajani2017improved) (e.g., gradient penalties), but found that a simple setup was most effective in practice. The overall procedure is outlined in the pseudocode of Algorithm [1](https://arxiv.org/html/2506.01502v2#algorithm1).

Figure: Algorithm 1 iJKOnet Training

Scalability. Our method scales well with the number of time points. Leveraging JAX (jax2018github) and Optax (deepmind2020jax), we avoid memory issues even for large batches and many time steps, since we do not backpropagate through time as in JKOnet (bunne2022proximal). Each step is processed independently, and the loss is computed by averaging outputs from separate networks, avoiding large computation graphs. Training time on GPU is comparable to $\texttt{JKOnet}^{\ast}$, but our approach eliminates the costly precomputation of OT couplings. In high-dimensional settings (50D–100D), $\texttt{JKOnet}^{\ast}$ exceeds GPU memory limits and requires CPU-based precomputations, making it substantially slower.

### C.3 Entropy estimation details

Prior to training, we estimated $\widehat{{\mathcal{H}}}(\rho_{k})$ for each $k=0,\dots,K$ using the Kozachenko–Leonenko nearest-neighbor estimator (kozachenko1987sample; berrett2019efficient) with 5 nearest neighbors. We used the publicly available implementation from (butakov2024mutual)(^2^22[https://github.com/VanessB/mutinfo](https://github.com/VanessB/mutinfo)). As discussed in \wasyparagraph[3.3](https://arxiv.org/html/2506.01502v2#S3.SS3), we use the following equation for estimating $\widehat{{\mathcal{H}}}(T_{\varphi}^{k}\sharp\rho_{k})$:

$$ $\widehat{{\mathcal{H}}}(T_{\varphi}^{k}\sharp\rho_{k})=\widehat{{\mathcal{H}}}(\rho_{k})-\int_{{\mathcal{X}}}\log|\det\nabla_{x}T_{\varphi}^{k}(x)|\,\differential\rho_{k}(x).$ (25) $$

Since our map is not required to be invertible during training, as is the case for gradient-based ICNN (amos2017input) parameterizations (bunne2022proximal), we compute the sign and the natural logarithm of the absolute value of the Jacobian determinant, both of which are supported by modern computational frameworks (jax2018github; paszke2019pytorch).

### C.4 Optimizer

We optimize two loss functions with respect to ${\mathcal{J}}_{\theta}$ and $T^{k}_{\varphi}$ using the Adam optimizer (kingma2014adam). For ${\mathcal{J}}_{\theta}$, we use parameters $\beta_{1}=0.9$, $\beta_{2}=0.999$, $\varepsilon=1\times 10^{-8}$, and a constant learning rate of $5\times 10^{-4}$, with gradient clipping applied using a maximum global norm of $10$. For $T^{k}_{\varphi}$, we use parameters $\beta_{1}=0.5$, $\beta_{2}=0.9$, and $\varepsilon=1\times 10^{-8}$, with a learning rate of $1\times 10^{-3}$. Training is performed with mini-batches of size $500$.

### C.5 Network architecture

The neural networks for the potential $V_{\theta_{1}}$ and interaction $W_{\theta_{2}}$ energies are multi-layer perceptrons (MLPs) with two hidden layers of size 64, using softplus activation functions, and a one-dimensional output layer. The neural network for the optimal transport maps $T^{k}_{\varphi}$ is also an MLP with two hidden layers of size 64. The time step $k$ is concatenated with the input $x^{i}_{k}$, and the network uses selu activations, with an output layer that matches the input dimension of $x^{i}_{k}$.

### C.6 Hardware

Experiments were conducted on a CentOS Linux 7 (Core) system with NVIDIA A100 GPUs. Most of the computation time was spent on metric evaluation. Leveraging JAX’s just-in-time compilation, our method completed $100$ training epochs in approximately one minute.

### C.7 Functionals

For easier representation of our experimental results, we use abbreviations for each potential, each abbreviation being shown in parentheses and the definitions can be found in (terpin2024learning, Appendix F): Wavy Plateau (WP), Double Exponential (DE), Rotational (RO), ReLU (RE), Flat (FT), Friedman (FN), Watershed (WS), Ishigami (IG), Flowers (FL), Bohachevsky (BC), Sphere (SP), Styblinski-Tang (ST), Oakley–Ohagan (OO), Zigzag Ridge (ZR), and Holder Table (HT).

## Appendix D Proofs

To begin with, we recall the notion of strong smoothness, see (beck2017first, §5.1). A function $V:{\mathcal{X}}\rightarrow\mathbb{R}$ is called $\frac{1}{\beta}$ - strongly smooth if it is continuously differentiable on ${\mathcal{X}}$ and $\forall x,y\in{\mathcal{X}}$ it holds:

$$ $\beta\|\nabla V(x)-\nabla V(y)\|_{2}\leq\|x-y\|_{2}.$ $$

Now we proceed to some auxiliary results.

###### Lemma D.1 (Solution to the JKO problem with potential energy is unique) .

Let ${\mathcal{J}}^{\ast}(\rho)=\int_{{\mathcal{X}}}V^{*}(x)\differential\rho(x)$, $\rho_{0}\in{\mathcal{P}}_{2,ac}({\mathcal{X}})$ and $\rho_{1}=\operatorname{JKO}_{\tau{\mathcal{J}}^{\ast}}(\rho_{0})$. Then $T_{V^{\ast}}=\operatorname*{arg\,min}_{T:{\mathcal{X}}\rightarrow{\mathcal{X}}}{\mathcal{L}}(V^{\ast},T)$ is the (unique) optimal transport map between $\rho_{0}$ and $\rho_{1}$. In particular, $T_{V^{\ast}}\sharp\rho_{0}=\rho_{1}$.

###### Proof.

Consider the functional $\rho\mapsto d_{{\mathbb{W}}_{2}}^{2}(\rho_{0},\rho)$, where $\rho\in{\mathcal{P}}_{2}({\mathcal{X}})$. In what follows, we prove that this functional is strictly convex.

Consider $\rho_{a},\rho_{b}\in{\mathcal{P}}_{2}({\mathcal{X}})$; $\alpha>0$. Let $\rho=\alpha\rho_{a}+(1-\alpha)\rho_{b}$. Note that $\rho_{0}$ is absolutely continuous. By the Brenier’s theorem, it holds:

$$ $\displaystyle d_{{\mathbb{W}}_{2}}^{2}(\rho_{0},\rho_{a})=\int\|x-T_{a}(x)\|_{2}^{2}d\rho_{0}(x)\,\,;\quad d_{{\mathbb{W}}_{2}}^{2}(\rho_{0},\rho_{b})=\int\|x-T_{b}(x)\|_{2}^{2}d\rho_{0}(x)$ $$

for the corresponding (unique) OT maps $T_{a}:T_{a}\sharp\rho_{0}=\rho_{a}$; $T_{b}:T_{b}\sharp\rho_{0}=\rho_{b}$.

Now consider the weighted sum:

$$ $\displaystyle\alpha d_{{\mathbb{W}}_{2}}^{2}(\rho_{0},\rho_{a})+(1-\alpha)d_{{\mathbb{W}}_{2}}^{2}(\rho_{0},\rho_{b})$ $\displaystyle=\int_{{\mathcal{X}}}\big(\alpha\|x-T_{a}(x)\|_{2}^{2}+(1-\alpha)\|x-T_{b}(x)\|_{2}^{2}\big)d\rho_{0}(x)$ $\displaystyle=\int_{{\mathcal{X}}}\|x-y\|_{2}^{2}d\pi_{ab}(x,y),$ $$

where plan $\pi_{ab}\in\Pi(\rho_{0},\rho)$ has the following conditionals:

$$ $\displaystyle\pi(y|x)=\alpha\delta\big(y=T_{a}(x)\big)+(1-\alpha)\delta\big(y=T_{b}(x)\big).$ $$

Since the OT map between $\rho_{0}$ and $\rho$ is unique (and deterministic), we have

$$ $\displaystyle d_{{\mathbb{W}}_{2}}^{2}(\rho_{0},\rho)=\min_{T\sharp\rho_{0}=\rho}\int_{{\mathcal{X}}}\|x-T(x)\|_{2}^{2}d\rho(x)<\int_{{\mathcal{X}}}\|x-y\|_{2}^{2}d\pi_{ab}(x,y),$ $$

which yields the strict convexity of $\rho\mapsto d_{{\mathbb{W}}_{2}}^{2}(\rho_{0},\rho)$.

The main statement of the Lemma follows from the strict convexity of ${\rho\mapsto\int_{{\mathcal{X}}}V^{\ast}(x)d\rho(x)+d_{{\mathbb{W}}_{2}}^{2}(\rho_{0},\rho)}$ and Brenier’s theorem.
∎

Now we are ready to prove our main quality bound theorem.

###### Proof.

(Theorem [3.1](https://arxiv.org/html/2506.01502v2#S3.Thmtheorem1)).

Prelimiary Note: The facts from convex analysis (properties of convex functions and their conjugates) which we use below could be found in (kornilov2024optimal, Lemma 1).

At first, we simplify the expression for the JKO loss:

$$ $\displaystyle{\mathcal{L}}(V,T_{V})$ $\displaystyle=\int V(T_{V}(x))d\rho_{0}(x)-\int V(y)d\rho_{1}(y)+\frac{1}{2\tau}\int\|x-T_{V}(x)\|_{2}^{2}d\rho_{0}(x)$ $\displaystyle=\frac{1}{\tau}\int V_{q}(T_{V}(x))d\rho_{0}(x)-\frac{1}{\tau}\int V_{q}(y)d\rho_{1}(y)-\frac{1}{\tau}\int\langle x,T_{V}(x)\rangle d\rho_{0}(x)$ (26) $\displaystyle\hskip 17.07164pt+\underbrace{\int\frac{\|y\|_{2}^{2}}{2\tau}d\rho_{1}(y)+\int\frac{\|x\|_{2}^{2}}{2\tau}d\rho_{0}(x)}_{\stackrel{{\scriptstyle\text{def}}}{{=}}C(\rho_{0},\rho_{1})},$ $$

where $V_{q}(y)=\tau V+\frac{1}{2}\|\cdot\|_{2}^{2}.$

From ([26](https://arxiv.org/html/2506.01502v2#A4.E26)) we note that $T_{V}=\operatorname*{arg\,min}_{T}\int\big(V_{q}(T_{V}(x))-\langle x,T_{V}(x)\rangle\big)d\rho_{0}(x)$. From the convex analysis, it follows that $T_{V}(x)=\nabla\overline{V_{q}}(x)$ delivers the minimum; $\overline{V_{q}}$ is the (Fenchel) conjugate of $V_{q}$. Note that $T_{V}$ is measurable since it is continuous almost surely (rockafellar-1970a, Thm.25.5). Substituting $\nabla\overline{V_{q}}$ into ([26](https://arxiv.org/html/2506.01502v2#A4.E26)) on par with Fenchel-Young equality yields:

$$ $\displaystyle{\mathcal{L}}(V,T_{V})$ $\displaystyle=-\frac{1}{\tau}\int\overline{V}_{q}(x)d\rho_{0}(x)-\frac{1}{\tau}\int V_{q}(y)d\rho_{1}(y)+C(\rho_{0},\rho_{1}).$ (27) $$

Note that $T_{V^{\ast}}\sharp\rho_{0}=\rho_{1}$ (Lemma [D.1](https://arxiv.org/html/2506.01502v2#A4.Thmtheorem1)). Therefore, eq. ([26](https://arxiv.org/html/2506.01502v2#A4.E26)) for ${\mathcal{L}}(V^{\ast},T_{V^{\ast}})$ could be simplified:

$$ $\displaystyle{\mathcal{L}}(V^{\ast},T_{V^{\ast}})=$ $\displaystyle\frac{1}{\tau}\underbrace{\int V_{q}^{\ast}(T_{V}(x))d\rho_{0}(x)}_{=\int V_{q}^{\ast}(y)d\rho_{1}(y)}-\frac{1}{\tau}\int V_{q}^{\ast}(y)d\rho_{1}(y)-\frac{1}{\tau}\int\langle x,T_{V^{\ast}}(x)\rangle d\rho_{0}(x)+C(\rho_{0},\rho_{1})=$ $\displaystyle-\frac{1}{\tau}\int\langle x,T_{V^{\ast}}(x)\rangle d\rho_{0}(x)+C(\rho_{0},\rho_{1}).$ (28) $$

Now we analyze the gap $\varepsilon(V)$ between optimal and optimized JKO losses. Leveraging ([27](https://arxiv.org/html/2506.01502v2#A4.E27)) and ([28](https://arxiv.org/html/2506.01502v2#A4.E28)) yields:

$$ $\displaystyle\tau\varepsilon(V)$ $\displaystyle=\tau{\mathcal{L}}(V^{\ast},T_{V^{\ast}})-\tau{\mathcal{L}}(V,T_{V})=\int\overline{V_{q}}(x)d\rho_{0}(x)+\int V_{q}(y)d\rho_{1}(y)-\int\langle x,T_{V^{\ast}}(x)\rangle d\rho_{0}(x).$ (29) $$

Now we note that $T_{V^{\ast}}=\nabla\overline{V^{\ast}_{q}}:{\mathcal{X}}\rightarrow\mathbb{R}^{D}$. From the properties of convex functions it follows that $\nabla V^{\ast}_{q}:{\mathcal{X}}\rightarrow\mathbb{R}^{D}$ defines the inverse Optimal Transport mapping. In particular (we assume $V^{\ast}_{q}$ to be convex) $\nabla V_{q}^{\ast}\sharp\rho_{1}=\rho_{0}$. Then, changing the variables in ([29](https://arxiv.org/html/2506.01502v2#A4.E29)) results in:

$$ $\displaystyle\tau\varepsilon(V)$ $\displaystyle=\int\overline{V_{q}}(\nabla V_{q}^{\ast}(y))d\rho_{1}(y)+\int V_{q}(y)d\rho_{1}(y)-\int\langle\nabla V_{q}^{\ast}(y),y\rangle d\rho_{1}(y)$ $\displaystyle=\int_{{\mathcal{X}}}\Big[\overline{V_{q}}(\nabla V_{q}^{\ast}(y))+V_{q}(y)-\langle\nabla V_{q}^{\ast}(y),y\rangle\Big]d\rho_{1}(y).$ (30) $\displaystyle=\int_{{\mathcal{X}}}{\mathcal{D}}_{\overline{V_{q}}}(\nabla V_{q}^{\ast}(y),\nabla V_{q}(y))d\rho_{1}(y),$ (31) $$

where ${\mathcal{D}}_{\overline{V_{q}}}(\,\cdot\,,\,\cdot\,)$ is the Bregman divergence, see (banerjee2005clustering, Def. 1) for the definition and (banerjee2005clustering, Appendix A) for the “Dual Divergences” property used in the transition from ([30](https://arxiv.org/html/2506.01502v2#A4.E30)) to ([31](https://arxiv.org/html/2506.01502v2#A4.E31)).

Since $V_{q}$ is $\frac{1}{\beta}$-smooth, then $\overline{V_{q}}$ is $\beta$-strongly convex. Therefore, by the property of strongly convex functions (we pick $y\in\text{supp}(\rho_{1})$):

$$ $\displaystyle\frac{\beta}{2}\|\nabla V_{q}^{\ast}(y)-\nabla V_{q}(y)\|_{2}^{2}$ $\displaystyle\leq\overline{V_{q}}(\nabla V_{q}^{\ast}(y))-\overline{V_{q}}(\nabla V_{q}(y))-\langle\nabla\overline{V_{q}}(\nabla V_{q}(y)),\nabla\overline{V_{q}^{\ast}}(y)-\nabla V_{q}(y)\rangle$ $\displaystyle={\mathcal{D}}_{\overline{V_{q}}}(\nabla V_{q}^{\ast}(y),\nabla V_{q}(y)).$ (32) $$

Combining ([31](https://arxiv.org/html/2506.01502v2#A4.E31)) and ([32](https://arxiv.org/html/2506.01502v2#A4.E32)) yields:

$$ $\displaystyle\tau\varepsilon(V)\geq\frac{\beta}{2}\int_{{\mathcal{X}}}\|\nabla V_{q}^{\ast}(y)-\nabla V_{q}(y)\|_{2}^{2}d\rho_{1}(y).$ $$

To complete the proof, we are left to note that $\nabla V_{q}^{\ast}(y)-\nabla V_{q}(y)=\tau\big(\nabla V^{\ast}(y)-\nabla V(y)\big)$.
∎

###### Lemma D.1 (Solution to the JKO problem with potential energy is unique) .

###### Proof.

###### Proof.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - In which regimes or datasets is JKO-based modeling justified or superior to flow-matching, neural ODEs, or Schrödinger bridge-based methods?
- How sensitive is iJKOnet to missing, irregularly spaced, or noisy population snapshots? Do learned energy functionals overfit to sampling artifacts in sparse settings?
soundness: 3
strengths: - iJKOnet is different from previous methods that either require potential-only energies or precomputed OT couplings. The formulation of energy functional recovery using inverse optimization within the JKO scheme provides a conceptually clear route to modeling population level dynamics from discrete snapshots, directly leveraging optimal transport geometry.
- Theorem 3.1 gives a non-trivial quality guarantee, explicitly bounding the distance between the gradients of learned and ground-truth potential functions in terms of the inverse optimization gap. This is a concrete advance compared to many contemporary works that lack precise recovery guarantees.
- iJKOnet dispenses with input-convex neural networks for parameterizing transport maps, allowing the use of more scalable and widely applicable neural network architectures such as MLPs and ResNets
summary: This paper introduces iJKOnet, an approach to learning population dynamics by bridging inverse optimization techniques with the Jordan-Kinderlehrer-Otto (JKO) scheme. The method frames the task as a min-max optimization problem for identifying the underlying energy functional that governs observed distributions, enabling end-to-end adversarial training without restrictive neural architectures. Theoretical guarantees are provided for recovery of potential-driven dynamics, and empirical results on both synthetic and real datasets show improved performance over prior JKO-based baselines.
weaknesses: - The main theoretical guarantee (Theorem 3.1) exclusively addresses potential energy functionals with strongly convex smooth potentials, while the practical method is claimed to work for broader energy functionals (e.g., including interaction and internal terms). This leaves a theoretical gap between what the analysis covers and what is demonstrated empirically. Specifically, the inability to provide quality bounds for non-potential and more general functionals limits the rigor of the claims. (This could explain why the performance of iJKOnet is much worse than baseline in Table 2)
- While iJKOnet is, in principle, applicable to full free-energy functionals, all main experiments revert to only learning potential energy, explicitly setting interaction and entropy terms to zero. As a result, the empirical support for learning richer dynamics is missing, and the scalability to very high-dimensional or truly multimodal interaction-dependent data is unproven.
- Some recent and directly related works [1,2,3,4] on similar task with unknown parameters, alternative population-level diffusion modeling, and advanced regularization techniques are not discussed.

ref:

[1] Computational and Statistical Asymptotic Analysis of the JKO Scheme for Iterative Algorithms to update distributions

[2] An Eulerian approach to regularized JKO scheme with low-rank tensor decompositions for Bayesian inversion

[3] WeightFlow: Learning Stochastic Dynamics via Evolving Weight of Neural Network

[4] Modeling Cell Dynamics and Interactions with Unbalanced Mean Field Schrodinger Bridge

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: - Could you explain about the potential stability of the proposed method resulted by the adversarial training procedure, and what further improvements could be made?

- Could you more explicitly summarize the key algorithmic or theoretical innovations that distinguish your approach from its predecessors?
soundness: 3
strengths: 1. This paper is well written.

2. The figures are well-prepared and greatly facilitate the understanding of the content.

3. The theoretical proofs is provided, which give solid theoretical guarantees of the proposed method.

4. The background section is clearly presented and offers helpful context.
summary: In this paper, the author propose  a method to learn the population dynamics, which combine JKO framework with inverse optimization techniques. In addition, this method relays on adversarial training procedure. Theoretical guarantees are also provided in this paper.
weaknesses: 1. The interchange of the min and Σ operators between Equations (10) and (11) lacks justification or analysis. It is unclear under what conditions this exchange is mathematically valid.

2. The method relies on an adversarial training procedure, which is known to introduce instability. It would be important to discuss whether any measures were taken to mitigate this issue.

3. The work builds heavily on existing studies (e.g., Terpin et al., 2024; JKOnet∗). The specific contributions and novelty of the present approach should be more clearly emphasized.

4. Several relevant references appear to be missing from the related work section, such as:

[1] Zhang Z, et al. Modeling Cell Dynamics and Interactions with Unbalanced Mean Field Schr\" odinger Bridge[J]. arXiv preprint arXiv:2505.11197, 2025.

[2] Li R, et al. WeightFlow: Learning Stochastic Dynamics via Evolving Weight of Neural Network[J]. arXiv preprint arXiv:2508.00451, 2025.

[3] Zhang Y, Levin M. Equilibrium flow: From Snapshots to Dynamics[J]. arXiv preprint arXiv:2509.17990, 2025.

---

confidence: 4
contribution: 2
initial_rating: 2
presentation: 3
questions: 1) Can you clarify what are the obstacles to implement an interaction kernel and the entropy term in the algorithm?

2) Can you comment on the difficulty of computing the transport maps?
soundness: 3
strengths: The JKO formulation of gradient flow as a proximal optimzation problem is a fundamental tool in probbaility theory.  It has been found difficult to use it effectively as a generating tool and and the reformulation of the JKO functional in this paper  is interesting and could be promising.
summary: The paper is an attempt to use the (famous) JKO optimization scheme to learn a dynamical system from a time- series $\rho_1, \cdots, \rho_M$  where $\rho_i$ represent the data of the system at time $T_1< T_2 < \cdots < T_M$.  The JKO scheme are expressed as proximal operator of the $W_2^2$ optimal cost with a functional (typically of the form energy plus entropy plus possibly interaction terms described by a kernel).   The approach taken here is to parametrize a transport map (for the $W^2$ term) and to parametrize  the energy terms using neural architectures and to estimate the entropy terms via an entropy estimator.
weaknesses: 1) The main weakness is that the paper does not fully implement the proposed loss functional. Indeed the the interaction terms and the entropy terms are ignored and the same architecture as in JKONet is used.  This reduces the class of PDEs considered to a very small class of linear PDE. The mismatch between the theory which considers general PDEs and the implementation is too great to make the paper convincing.   Ignoring the entropy term,  which I suspect is very difficult to estimate,  is a serious defect of the approach.  

2) The proposed architecture requires the construction of a transport map which is both expensive and only works theoretically when the measures have density.  The authors propose to replace input convex neural networks (since the transport maps are gradient of a convex functions) by a more general architecture but I do not see why this will make the problem more scalable or stable.  Maybe a approach using triangular map as in the Yousef Mazrouk group at MIT may be more scalable?


Decision: Accept (Poster)