Importance Corrected Neural JKO Sampling

## Abstract

Abstract In order to sample from an unnormalized probability density function, we propose to combine continuous normalizing flows (CNFs) with rejection-resampling steps based on importance weights. We relate the iterative training of CNFs with regularized velocity fields to a JKO scheme and prove convergence of the involved velocity fields to the velocity field of the Wasserstein gradient flow (WGF).
The alternation of local flow steps and non-local rejection-resampling steps allows to overcome local minima or slow convergence of the WGF for multimodal distributions.
Since the proposal of the rejection step is generated by the model itself,
they do not suffer from common drawbacks of classical rejection schemes.
The arising model can be trained iteratively, reduces the reverse Kulback-Leibler (KL) loss function in each step, allows to generate iid samples and moreover allows for evaluations of the generated underlying density.
Numerical examples show that our method yields accurate results on various test distributions including high-dimensional multimodal targets and outperforms the state of the art in almost all cases significantly.

## 1 Introduction

We consider the problem of sampling from an unnormalized probability density function. That is, we are given an integrable function $g\colon\mathbb{R}^{d}\to\mathbb{R}_{>0}$ and we aim to generate samples from the probability distribution $\nu$ given by the density $q(x)=g(x)/Z_{g}$, where the normalizing constant $Z_{g}=\int_{\mathbb{R}^{d}}g(x)\mathrm{d}x$ is unknown.
Many classical sampling methods are based on Markov chain Monte Carlo (MCMC) methods like the overdamped Langevin sampling, see, e.g., Welling & Teh (). The generated probability path of the underlying stochastic differential equation follows the Wasserstein-$2$ gradient flow of the reverse KL divergence $\mathcal{F}(\mu)=\mathrm{KL}(\mu,\nu)$.
Over the last years, generative models like normalizing flows (Rezende & Mohamed, ) or diffusion models (Ho et al., ; Song et al., ) became more popular for sampling, see, e.g., Phillips et al. (); Vargas et al. (). Also these methods are based on the reverse KL divergence as a loss function.
While generative models have successfully been applied in data-driven setups, their application to the problem of sampling from arbitrary unnormalized densities is not straightforward. This difficulty arises from the significantly harder nature of the problem, even in moderate dimensions, particularly when dealing with target distributions that exhibit phenomena such as concentration effects, multi-modalities, heavy tails, or other issues related to the curse of dimensionality.

In particular, the reverse KL is non-convex in the Wasserstein space as soon as the target density $\nu$ is not log-concave which is for example the case when $\nu$ consists of multiple modes.
In this case generative models often collapse to one or a small number of modes. We observe that for continuous normalizing flows (CNFs, Chen et al., ; Grathwohl et al., ) this can be prevented by regularizing the $L^{2}$-norm of the velocity field as proposed under the name OT-flow by Onken et al. (). In particular, this regularization converts the objective functional into a convex one.
However, the minimizer of the regularized loss function is no longer given by the target measure $\nu$ but by the Wasserstein proximal mapping of the objective function applied onto the latent distribution. Considering that the Jordan-Kinderlehrer-Otto (JKO) scheme (Jordan et al., ) iteratively applies the Wasserstein proximal mapping and converges to the Wasserstein gradient flow, several papers proposed to approximate the steps of the scheme by generative models, see Altekrüger et al. (); Alvarez-Melis et al. (); Fan et al. (); Lambert et al. (); Mokrov et al. (); Vidal et al. (); Xu et al. (). We will refer to this class of method by the name neural JKO. Even though this approximates the same gradient flow of the Langevin dynamics these approaches have usually the advantage of faster inference (once they are trained) and additional possibly allow for density evaluations.

However, already the time-continuous Wasserstein gradient flow suffers from the non-convexity of the reverse KL loss function by getting stuck in local minima or by very slow convergence times.
In particular, it is well-known that Langevin-based sampling methods often do not distribute the mass correctly onto multimodal target distributions.
As a remedy, Neal () proposed importance sampling, i.e., to reweight the sample based on quotients of the target distribution and its current approximation which leads to an unbiased estimator of the Monte Carlo integral. However, this estimator may lead to highly imbalanced weights between these samples and the resulting estimator might have a large variance. Moreover, the density of the current approximation has to be known up to a possible multiplicative constant, which is not the case for many MCMC methods like Langevin sampling or Hamiltonian Monte Carlo. In Del Moral et al. () the authors propose a scheme for iteratively applying the importance sampling procedure. However, as a downside this strategy might lead to samples which might not be iid.

#### Contributions

Figure: Figure 1: Iterative application of neural JKO steps and rejection steps for a shifted mixture target distribution. The red samples are rejected in the next following rejection step and the green samples are the resampled points. The latter approach enables for the correction of wrong mode weights introduced by the underlying WGF. See Figure 2 for more steps.
Refer to caption: extracted/5762551/introfig/intro0.png

In this paper, we propose a sampling method which combines neural JKO steps based on CNFs with importance based rejection steps.
While the CNFs adjust the position of the generated samples locally, the rejection steps readjusts the inferred distribution non-locally based on the quotient of generated and target distribution. Then, in each rejection step we resample the rejected points based on the current constructed generative model.
We illustrate this procedure in Figure 1.

Our methods generates independent samples and allows to evaluate the density of the generated distribution.
In our numerical examples, we apply our method to common test distributions up to the dimension $d=1600$ which are partially highly multimodal. We show that our importance corrected neural JKO sampling (neural JKO IC) achieves significantly better results than the comparisons(^1^11The code is available at [https://github.com/johertrich/neural_JKO_ic](https://github.com/johertrich/neural_JKO_ic)).

From a theoretical side, we prove that the velocity fields from a sequence of neural JKO steps strongly converges to the velocity field of the corresponding Wasserstein gradient flow and that the reverse KL loss function decreases throughout the importance-based rejection steps.

#### Related Work

Common methods for sampling of unnormalized densities are often based on Markov Chain Monte Carlo (MCMC) methods, see e.g. (Gilks et al., ). In particular first order based variants, such as the Hamiltonian Monte Carlo (HMC)
(Betancourt, ; Hoffman & Gelman, )
and the Metropolis Adjusted Langevin Algorithm (MALA Girolami & Calderhead, ; Rossky et al., ; Roberts & Tweedie, ) are heavily used in practice.
The viewpoint of these samplers
as sample space description of gradient flows defined in a metricized probability space then allows for extensions such as interacting particle systems (Chen et al., ; Eigel et al., ; Garbuno-Inigo et al., ; Wang & Li, ).
However, since these algorithms are based on local transformations of the samples, they are unable to distribute the mass correctly among different modes, which can partially be corrected by importance sampling (Neal, ) and sequential Monte Carlo samplers (Del Moral et al., ) as described above.

As an alternative to MCMC methods one can use generative models for sampling including VAEs (Kingma & Welling, ), normalizing flows (Rezende & Mohamed, ) or diffusion models based on variational approaches (Blessing et al., ; Phillips et al., ; Vargas et al., ) or by computing the score by solving a partial differential equation
(Richter & Berner, ; Sommer et al., ). These methods usually provide much faster sampling time and are often used in combination with some conditioning parameter for inverse problems, where a (generative) prior is combined with a known likelihood term (Ardizzone et al., ; Altekrüger & Hertrich, ; Andrle et al., ; Denker et al., ).
These generative models can be adapted to follow a Wasserstein gradient flow, by mimicking a JKO scheme with generative models (Altekrüger et al., ; Alvarez-Melis et al., ; Fan et al., ; Lambert et al., ; Mokrov et al., ; Vidal et al., ; Xu et al., ) or directly follow the velocity field of a kernel-based functional (Arbel et al., ; Hagemann et al., ; Hertrich et al., ; Liu & Wang, ; Liutkus et al., ).
Combinations of generative models with stochastic sampling steps were considered in the literature under the names stochastic normalizing flows (Hagemann et al., ; ; Wu et al., ) and annealed flow transport Monte Carlo (Arbel et al., ; Matthews et al., ).

#### Outline

The paper is organized as follows.
In Section 2, we recall the fundamental concepts which will be required. Afterwards, we consider neural JKO schemes more detailed in Section 3.
We introduce our importance-based rejection steps in Section 4. Finally, we evaluate our model numerically and compare it to existing methods in Section 5. Conclusions are drawn in Section 6. Additionally, proofs, further numerical and technical details are presented in Appendix A- E.

## 2 Preliminaries

In this section, we provide a rough overview of the required concepts for this paper. To this end, we first revisit the basic definitions of Wasserstein gradient flows, e.g., based on Ambrosio et al. (). Afterwards we recall continuous normalizing flows with OT-regularizations.

### 2.1 Wasserstein Spaces and Absolutely Continuous Curves

#### Wasserstein Distance

Let $\mathcal{P}(\mathbb{R}^{d})$ be the space of probability measures on $\mathbb{R}^{d}$ and denote by $\mathcal{P}_{2}(\mathbb{R}^{d})\coloneqq\{\mu\in\mathcal{P}(\mathbb{R}^{d}):
\int_{\mathbb{R}^{d}}\|x\|^{2}\mathrm{d}x<\infty\}$ the subspace of probability measures with finite second moment. Let $\mathcal{P}_{2}^{\mathrm{ac}}(\mathbb{R}^{d})$ be the subspace of absolutely continuous measures from $\mathcal{P}_{2}(\mathbb{R}^{d})$.
Moreover, we denote for $\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ by $\Gamma(\mu,\nu)\coloneqq\{\bm{\pi}\in\mathcal{P}_{2}(\mathbb{R}^{d}\times
\mathbb{R}^{d}):{\pi_{1}}_{\#}\bm{\pi}=\mu,{\pi_{2}}_{\#}\bm{\pi}=\nu\}$ the set of all transport plans with marginals $\mu$ and $\nu$, where $\pi_{i}\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ defined by $\pi_{i}(x_{1},x_{2})=x_{i}$ is the projection onto the $i$-th component for $i=1,2$.
Then, we equip $\mathcal{P}_{2}(\mathbb{R}^{d})$ with the Wasserstein-2 metric defined by

$$ $W_{2}^{2}(\mu,\nu)=\inf_{\bm{\pi}\in\Gamma(\mu,\nu)}\int_{\mathbb{R}^{d}\times \mathbb{R}^{d}}\|x-y\|^{2}\mathrm{d}\bm{\pi}(x,y).$ $$

The minimum is attained, but not always unique. We denote the subset of minimizers by $\Gamma^{\mathrm{opt}}(\mu,\nu)$.
If $\mu\in\mathcal{P}_{2}^{\mathrm{ac}}(\mathbb{R}^{d})$, then the optimal transport plan is unique
and is given by a so-called transport map.

###### Theorem 1 ( Brenier ( 1987 ) ) .

Let $\mu\in\mathcal{P}_{2}^{\mathrm{ac}}(\mathbb{R}^{d})$ and $\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})$.
Then there is a unique transport plan $\bm{\pi}\in\Gamma^{\rm{opt}}(\mu,\nu)$
which is induced by a unique
measurable optimal transport map $T\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$, i.e.,
$\bm{\pi}=({\rm{Id}},T)_{\#}\mu$
and

$$ $W_{2}^{2}(\mu,\nu)=\min_{T\colon\mathbb{R}^{d}\to\mathbb{R}^{d}}\int_{\mathbb{ R}^{d}}\|T(x)-x\|_{2}^{2}\,\mathrm{d}\mu(x)\quad\text{subject \; to}\quad T_{ \#}\mu=\nu.$ (1) $$

###### Theorem 1 ( Brenier ( 1987 ) ) .

#### Absolutely Continuous Curves

A curve $\gamma\colon I\to\mathcal{P}_{2}(\mathbb{R}^{d})$
on the interval $I\subseteq\mathbb{R}$ is called *absolutely continuous*
if there exists a Borel velocity field $v\colon\mathbb{R}^{d}\times I\to\mathbb{R}^{d}$ with
$\int_{I}\|v(\cdot,t)\|_{L_{2}(\gamma(t),\mathbb{R}^{d})}\mathrm{d}t<\infty$
such that the continuity equation

$$ $\partial_{t}\gamma(t)+\nabla\cdot(v(\cdot,t)\gamma(t))=0$ (2) $$

is fulfilled on $I\times\mathbb{R}^{d}$ in a weak sense.
Then, any velocity field $v$ solving the continuity equation (2) for fixed $\gamma$ characterizes $\gamma$ as
$\gamma(t)=z(\cdot,t)_{\#}\gamma(t_{0}),$
where $z$ is the solution of the ODE $\dot{z}(x,t)=v(z(x,t),t)$ with $z(x,t_{0})=x$ and $t_{0}\in I$.
It can be shown that for an absolutely continuous curve there exists a unique solution of minimal norm which is equivalently characterized by the so-called regular tangent space $\mathrm{T}_{\gamma(t)}\mathcal{P}_{2}(\mathbb{R}^{d})$, see Appendix A for details.
An absolute continuous curve is a geodesic if there exists some $c>0$ such that $W_{2}(\gamma(s),\gamma(t))=c|s-t|$.

The following theorem formulates a dynamic version of the Wasserstein distance based minimal energy curves in the Wasserstein space.

###### Theorem 2 ( Benamou & Brenier ( 2000 ) ) .

Assume that $\mu,\nu\in\mathcal{P}_{2}^{\mathrm{ac}}(\mathbb{R}^{d})$. Then, it holds

$$ $W_{2}^{2}(\mu,\nu)=\inf_{\begin{subarray}{c}v\colon\mathbb{R}^{d}\times[0,1] \to\mathbb{R}^{d},\\ \dot{z}(x,t)=v(z(x,t),t),\\ z(x,0)=x,\,z(\cdot,1)_{\#}\mu=\nu\end{subarray}}\int_{0}^{1}\int_{\mathbb{R}^{ d}}\|v(z(x,t),t)\|^{2}\mathrm{d}\mu(x)\mathrm{d}t.$ $$

Moreover, there exists a unique minimizing velocity field $v$ and the curve defined by $\gamma(t)=z(\cdot,t)_{\#}\mu$ with $t\in[0,1]$ and $\dot{z}(x,t)=v(z(x,t),t)$, $z(x,0)=x$ is a geodesic which fulfills the continuity equation
$\partial_{t}\gamma(t)+\nabla\cdot(v(\cdot,t)\gamma(t))=0.$

Let $\tau>0$. Then, by substitution of $t$ by $t/\tau$ and rescaling $v$ in the time variable, this is equal to

$$ $\displaystyle W_{2}^{2}(\mu,\nu)$ $\displaystyle=\inf_{\begin{subarray}{c}v\colon\mathbb{R}^{d}\times[0,\tau]\to \mathbb{R}^{d},\\ \dot{z}(x,t)=v(z(x,t),t),\\ z(x,0)=x,\,z(\cdot,\tau)_{\#}\mu=\nu\end{subarray}}\tau\int_{0}^{\tau}\int_{ \mathbb{R}^{d}}\|v(z(x,t),t)\|^{2}\mathrm{d}\mu(x)\mathrm{d}t.$ (3) $$

###### Theorem 2 ( Benamou & Brenier ( 2000 ) ) .

#### Wasserstein Gradient Flows

An absolutely continuous curve $\gamma\colon(0,\infty)\to\mathcal{P}_{2}(\mathbb{R}^{d})$
with velocity field $v_{t}\in\mathrm{T}_{\gamma(t)}\mathcal{P}_{2}(\mathbb{R}^{d})$ is a *Wasserstein gradient flow with respect to* $\mathcal{F}\colon\mathcal{P}_{2}(\mathbb{R}^{d})\to(-\infty,\infty]$
if

$$ $v_{t}\in-\partial\mathcal{F}(\gamma(t)),\quad\text{for a.e. }t>0,$ (4) $$

where $\partial\mathcal{F}(\mu)$ denotes the reduced Fréchet subdiffential at $\mu$, see Appendix A for a definition.

To compute Wasserstein gradient flows numerically, we can use the generalized minimizing movements or Jordan-Kinderlehrer-Otto (JKO) scheme (Jordan et al., ). To this end, we consider the Wasserstein proximal mapping defined as

$$ $\operatorname{prox}_{\tau\mathcal{F}}(\hat{\mu})=\operatorname*{arg\,min}_{\mu \in\mathcal{P}_{2}(\mathbb{R}^{d})}\{\frac{1}{2}W_{2}^{2}(\mu,\hat{\mu})+\tau \mathcal{F}(\mu)\}.$ $$

Then, define as $\mu_{\tau}^{k}$ for $k\in\mathbb{N}$ the steps of the minimizing movements scheme, i.e.,

$$ $\mu_{\tau}^{0}=\mu^{0},\qquad\mu_{\tau}^{k+1}=\operatorname{prox}_{\tau \mathcal{F}}(\mu_{\tau}^{k}).$ (5) $$

We denote the piecewise constant interpolations
$\tilde{\gamma}_{\tau}\colon[0,\infty)\to\mathcal{P}_{2}(\mathbb{R}^{d})$
of the minimizing movement scheme by

$$ $\tilde{\gamma}_{\tau}(k\tau+t\tau)=\mu_{\tau}^{k},\quad t\in[0,1).$ (6) $$

Then, the following convergence result holds true.

###### Theorem 3 .

(Ambrosio et al., , Thm 11.2.1)
Let $\mathcal{F}\colon\mathcal{P}_{2}(\mathbb{R}^{d})\to(-\infty,+\infty]$ be
proper,
lsc,
coercive,
and $\lambda$-convex along generalized geodesics,
and let $\mu^{0}\in\overline{\operatorname{dom}\mathcal{F}}$.
Then the curves $\gamma_{\tau}$ defined via the minimizing movement scheme (6)
converge for $\tau\to 0$ locally uniformly to
a locally Lipschitz curve $\gamma\colon(0,+\infty)\to\mathcal{P}_{2}(\mathbb{R}^{d})$
which is the unique Wasserstein gradient flow of $\mathcal{F}$ with $\gamma(0+)=\mu^{0}$.

###### Theorem 3 .

### 2.2 Continuous Normalizing Flows and OT-Flows

The concept of normalizing flows first appeared in Rezende & Mohamed ().
It follows the basic idea to approximate a probability distribution $\nu$ by considering a simple latent distribution $\mu_{0}$ (usually a standard Gaussian) and to construct a diffeomorphism $\mathcal{T}_{\theta}\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ depending on some parameters $\theta$ such that $\nu\approx{\mathcal{T}_{\theta}}_{\#}\mu_{0}$.
In practice, the diffeomorphism can be approximated by coupling-based neural networks (Dinh et al., ; Kingma & Dhariwal, ), residual architectures (Behrmann et al., ; Chen et al., ; Hertrich, ) or autoregressive flows (De Cao et al., ; Durkan et al., ; Huang et al., ; Papamakarios et al., ).
In this paper, we mainly focus on *continuous normalizing flows* proposed by Chen et al. (); Grathwohl et al. (), see also Ruthotto & Haber () for an overview. Here the diffeomorphism $\mathcal{T}_{\theta}$ is parameterized as neural ODE.
To this end let $v_{\theta}\colon\mathbb{R}^{d}\times\mathbb{R}\to\mathbb{R}^{d}$ be a neural network with parameters $\theta$ and let
$z_{\theta}\colon\mathbb{R}^{d}\times[0,\tau]\to\mathbb{R}^{d}$ for fixed $\tau>0$ be the solution of
$\dot{z}_{\theta}=v_{\theta}$, with initial condition $z_{\theta}(x,0)=x$.
Then, we define $\mathcal{T}_{\theta}$ via the solution $z_{\theta}$ as $\mathcal{T}_{\theta}(x)=z_{\theta}(x,\tau)$.

The density $p_{\theta}$ of ${\mathcal{T}_{\theta}}_{\#}\mu_{0}$ can be described by the change-of-variables formula

$$ $p_{\theta}(x)=\frac{p_{0}(x)}{|\mathrm{det}(\nabla\mathcal{T}_{\theta}(x))|},$ $$

where $p_{0}$ is the density of the latent distribution $\mu_{0}$. In the case of continuous normalizing flows, it can be shown that the denominator can be computed as

$$ $\log(|\mathrm{det}(\nabla\mathcal{T}_{\theta}(x))|)=\ell_{\theta}(x,\tau), \quad\partial_{t}\ell_{\theta}(x,t)=\mathrm{trace}(\nabla v_{\theta}(z_{\theta }(x,t),t)),\quad\ell_{\theta}(\cdot,0)=0.$ $$

In order to train a normalizing flow, one usually uses the Kullback-Leibler divergence. If $\nu$ is given by a density $q(x)=Z_{g}g(x)$, where $Z_{g}$ is an unknown normalizing constant, then this amounts to the reverse KL loss function

$$ $\displaystyle\mathcal{L}(\theta)$ $\displaystyle=\mathrm{KL}({\mathcal{T}_{\theta}}_{\#}\mu_{0},\nu)=\mathbb{E}_{ x\sim\mu_{0}}[-\log(q(\mathcal{T}_{\theta}(x)))+\log(p_{\theta}(x))]$ (7) $\displaystyle\propto\mathbb{E}_{x\sim\mu_{0}}[-\log(q(\mathcal{T}_{\theta}(x)) )-\ell_{\theta}(x,\tau)],$ (8) $$

where $\propto$ indicates equality up to a constant.

In order to stabilize the training, Onken et al. () propose to regularize the velocity field $v_{\theta}$ by its expected squared norm. More precisely, they propose to add the regularizer

$$ $\mathcal{R}(\theta)=\tau\int_{0}^{\tau}\|v_{\theta}(z_{\theta}(x,t),t)\|^{2} \mathrm{d}t$ $$

to the loss function.
Following Theorem 2, the authors of Onken et al. () note that for $\beta>0$ the functional $\mathcal{L}(\theta)+\alpha\mathcal{R}(\theta)$ has the same minimizer as the functional $\mathcal{L}(\theta)+\beta W_{2}^{2}(\mu_{0},{\mathcal{T}_{\theta}}_{\#}\mu_{0})$, which relates to the JKO scheme as pointed out by Vidal et al. ().

## 3 Neural JKO Scheme

In the following, we learn the steps (5) of the JKO scheme by neural ODEs. While similar schemes were already suggested in several papers (Altekrüger et al., ; Alvarez-Melis et al., ; Fan et al., ; Lambert et al., ; Mokrov et al., ; Vidal et al., ; Xu et al., ), we are particularly interested in the convergence properties of the corresponding velocity fields.
In Subsection 3.1, we introduce the general scheme and derive its properties. Afterwards, in Subsection 3.2, we describe the corresponding neural network approximation.

Throughout this section, we consider the following assumptions on the objective functional $\mathcal{F}$ and our initialization $\mu_{0}$.

###### Assumption 4 .

Let $\mathcal{F}\colon\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}\cup\{\infty\}$ be proper, lower semi-continuous with respect to narrow convergence, coercive, $\lambda$-convex along generalized geodesics and bounded from below. Moreover, assume that $\operatorname{dom}(|\partial\mathcal{F}|)\subseteq\mathcal{P}_{2}^{\mathrm{ac}
}(\mathbb{R}^{d})$ and that $\mathcal{F}$ has finite metric derivative $|\partial\mathcal{F}|(\mu_{0})<\infty$ at the initialization $\mu_{0}\in\mathcal{P}_{2}^{\mathrm{ac}}(\mathbb{R}^{d})$.

This assumption is fulfilled for many important divergences and loss functions $\mathcal{F}$. We list some examples in Appendix B.1.
We will later pay particular attention to the reverse Kullback-Leibler divergence
$\mathcal{F}(\mu)=\mathrm{KL}(\mu,\nu)=\int p(x)\log\left(\frac{p(x)}{q(x)}
\right)\mathrm{d}x$,
where $\nu$ is a fixed target measure and $p$ and $q$ are the densities of $\mu$ and $\nu$ respectively. This functional fulfills Assumption 4 if $q$ is $\lambda$-convex.

###### Assumption 4 .

### 3.1 Piecewise Geodesic Interpolation

In order to represent the JKO scheme by neural ODEs, we first reformulate it based on Benamou-Brenier (Theorem 2).
To this end, we insert the dynamic formulation of the Wasserstein distance in the Wasserstein proximal mapping defining the steps in (5). For any $\mu\in\mathcal{P}(\mathbb{R}^{d})$ this leads to

$$ $\frac{1}{2\tau}W_{2}^{2}(\mu,\mu_{\tau}^{k})+\mathcal{F}(\mu)=\min_{\begin{ subarray}{c}v\colon\mathbb{R}^{d}\times[0,\tau]\to\mathbb{R}^{d}\\ \dot{z}(x,t)=v(z(x,t),t),\\ z(x,0)=x,\,z(\cdot,\tau)_{\#}\mu_{\tau}^{k}=\mu\end{subarray}}\frac{1}{2}\int_ {0}^{\tau}\int_{\mathbb{R}^{d}}\|v(z(x,t),t)\|^{2}\mathrm{d}\mu_{\tau}^{k}(x) \mathrm{d}t+\mathcal{F}(\mu).$ (9) $$

Hence from a minimizer perspective, we obtain that
$\mu_{\tau}^{k+1}=z_{\tau}^{k}(\cdot,T)_{\#}\mu_{\tau}^{k},$
where

$$ $(v_{\tau,k},z_{\tau,k})\in\operatorname*{arg\,min}_{\begin{subarray}{c}v\colon \mathbb{R}^{d}\times[0,\tau]\to\mathbb{R}^{d}\\ \dot{z}(x,t)=v(z(x,t),t),\,z(x,0)=x\end{subarray}}\left\{\frac{1}{2}\int_{0}^{ \tau}\int_{\mathbb{R}^{d}}\|v(z(x,t),t)\|^{2}\mathrm{d}\mu_{\tau}^{k}(x) \mathrm{d}t+\mathcal{F}(z(\cdot,\tau)_{\#}\mu_{\tau}^{k})\right\}.$ (10) $$

Finally, we concatenate the velocity fields of all steps and obtain the ODE

$$ $v_{\tau}|_{(k\tau,(k+1)\tau]}=v_{\tau,k},\qquad\dot{z}_{\tau}(x,t)=v_{\tau}(z_ {\tau}(x,t),t),\qquad z_{\tau}(x,0)=x.$ $$

As a straightforward observation, we obtain that the curve defined by the velocity field $v_{\tau}$ is the geodesic interpolation between the points from JKO scheme (5).
We state a proof in Appendix B.2.

###### Corollary 5 .

Under Assumption 4 the following holds true.

- (i)
It holds that
$W_{2}^{2}(\mu_{\tau}^{k},\mu_{\tau}^{k+1})=\tau\int_{0}^{\tau}\int_{\mathbb{R}
^{d}}\|v_{\tau,k}(z_{\tau,k}(x,t),t)\|^{2}\mathrm{d}\mu_{\tau}^{k}(x)\mathrm{d}t$,
i.e., $v_{\tau,k}$ is the optimal velocity field from the theorem of Benamou-Brenier.
- (ii)
The curve $\gamma_{\tau}(t)\coloneqq z(\cdot,t)_{\#}\mu_{0}$ fulfills
$\gamma_{\tau}(k\tau+t\tau)=((1-t)I+t\mathcal{T}_{\tau}^{k})_{\#}\mu_{\tau}^{k}$ for $t\in[0,1],$
where $\mathcal{T}_{\tau}^{k}$ is the optimal transport map between $\mu_{\tau}^{k}$ and $\mu_{\tau}^{k+1}$.
- (iii)
$v_{\tau}$ and $\gamma_{\tau}$ solve the continuity equation
$\partial_{t}\gamma_{\tau}(t)+\nabla\cdot(v_{\tau}(\cdot,t)\gamma(t))=0.$

Analogously to Theorem 3 one can show that also the curves $\gamma_{\tau}$ are converging locally uniformly to the unique Wasserstein gradient flow (see, e.g., the proof of Ambrosio et al., , Thm 11.1.6).
In the next subsection, we will approximate the velocity fields $v_{\tau,k}$ by neural networks.
In order to retain the stability of the resulting scheme, the next theorem states that also the velocity fields $v_{\tau}$ converge strongly towards the velocity field from the Wasserstein gradient flow. The proof is given in Appendix B.3.

###### Theorem 6 .

Suppose that Assumption 4 is fulfilled and
let $(\tau_{l})_{l}\subseteq(0,\infty)$ with $\tau_{l}\to 0$. Then, $(v_{\tau_{l}})_{l}$ converges strongly to the velocity field $\hat{v}\in L_{2}(\gamma,\mathbb{R}^{d}\times[0,T])$ of the Wasserstein gradient flow $\gamma\colon(0,\infty)\to\mathcal{P}_{2}(\mathbb{R}^{d})$ of $\mathcal{F}$ starting in $\mu^{0}$.

In some cases the limit velocity field $v$ can be stated explicitly, even though its direct computation is intractable. For details, we refer to Appendix B.1.

###### Corollary 5 .

###### Theorem 6 .

### 3.2 Neural JKO Sampling

In the following, we learn the velocity fields $v_{\tau,k}$ as neural ODEs in order to sample from a target measure $\nu$ given by the density $q(x)=\frac{1}{Z_{g}}g(x)$ with unknown normalization constant $Z_{g}=\int_{\mathbb{R}^{d}}g(x)\mathrm{d}x$.
To this end, we consider the Wasserstein gradient flow with respect to the reverse KL loss function $\mathcal{F}(\mu)=\mathrm{KL}(\mu,\nu)$ which has the unique minimizer $\mu=\nu$.

Then, due to (10) the loss function from the JKO steps for the training of the velocity field $v_{\theta}$ reads as

$$ $\mathcal{L}(\theta)=\frac{1}{2}\int_{0}^{\tau}\int_{\mathbb{R}^{d}}\|v_{\theta }(z_{\theta}(x,t),t)\|^{2}\mathrm{d}\mu_{\tau}^{k}(x)\mathrm{d}t+\mathcal{F}(z _{\theta}(\cdot,\tau)_{\#}\mu_{\tau}^{k}),$ $$

where $z_{\theta}$ is the solution of $\dot{z}_{\theta}(x,t)=v_{\theta}(z_{\theta}(x,t),t)$ with $z_{\theta}(x,0)=x$.
Now, following the derivations of continuous normalizing flows, cf. Section 2.2, the second term of $\mathcal{L}$ can be rewritten (up to an additive constant) as

$$ $\mathcal{F}(z_{\theta}(\cdot,\tau)_{\#}\mu_{\tau}^{k})\propto\mathbb{E}_{x\sim \mu_{\tau}^{k}}[-\log(g(z_{\theta}(x,\tau)))-\ell_{\theta}(x,\tau)],$ $$

where $\ell_{\theta}$ is the solution of $\dot{\ell}_{\theta}(x,t)=\mathrm{trace}(\nabla v_{\theta}(z_{\theta}(x,t),t))$ with $\ell_{\theta}(\cdot,0)=0$.
Moreover, we can rewrite the first term of $\mathcal{L}$ based on

$$ $\int_{0}^{\tau}\int_{\mathbb{R}^{d}}\|v_{\theta}(z_{\theta}(x,t),t)\|^{2} \mathrm{d}\mu_{\tau}^{k}(x)\mathrm{d}t=\mathbb{E}_{x\sim\mu_{\tau}^{k}}[\int_{ 0}^{\tau}\|v_{\theta}(z_{\theta}(x,t),t)\|^{2}\mathrm{d}t]=\mathbb{E}_{x\sim \mu_{\tau}^{k}}[\omega_{\theta}(x,\tau)],$ $$

where $\omega_{\theta}$ is the solution of $\dot{\omega}_{\theta}(x,t)=\|v_{\theta}(z_{\theta}(x,t),t)\|^{2}$.
Hence, we can represent $\mathcal{L}$ up to an additive constant as

$$ $\mathcal{L}(\theta)=\mathbb{E}_{x\sim\mu_{\tau}^{k}}[-\log(g(z_{\theta}(x,\tau )))-\ell_{\theta}(x,\tau)+\omega_{\theta}(x,\tau)],$ (11) $$

where $(z_{\theta},\ell_{\theta},\omega_{\theta})$ solves the ODE system

$$ $\left(\begin{array}[]{c}\dot{z}_{\theta}(x,t)\\ \dot{\ell}_{\theta}(x,t)\\ \dot{\omega}_{\theta}(x,t)\end{array}\right)=\left(\begin{array}[]{c}v_{\theta }(x,t)\\ \mathrm{trace}(\nabla v_{\theta}(z_{\theta}(x,t),t))\\ \|v_{\theta}(z_{\theta}(x,t),t)\|^{2}\end{array}\right),\qquad\left(\begin{ array}[]{c}z_{\theta}(x,0)\\ \ell_{\theta}(x,0)\\ \omega_{\theta}(x,0)\end{array}\right)=\left(\begin{array}[]{c}x\\ 0\\ 0\end{array}\right).$ (12) $$

In particular, the loss function $\mathcal{L}$ can be evaluated and differentiated based on samples from $\mu_{\tau}^{k}$.
Once the parameters $\theta$ are optimized, we can evaluate the JKO steps in the same way as standard continuous normalizing flows. We summarize training and evaluation of the JKO steps in Algorithm 1 and 2 in Appendix D.1. In practice, the density values are computed and stored in log-space for numerical stability.

## 4 Importance-Based Rejection Steps

While a large number of existing sampling methods rely on Wasserstein gradient flows with respect to some divergences, it is well known that these loss functions are non-convex. This leads to very slow divergence speeds or convergence to suboptimal local minima.
In particular, if the target distribution is multi-modal the modes often do not have the right mass assigned.

#### Annealed Importance Sampling

As a remedy, Neal () proposed to use annealed importance sampling.
This means that we assign an importance weight $w_{i}$ to each generated sample $x_{i}$, which is computed by
$w_{i}=\frac{q(x)}{p(x)}$, where $p$ is some proposal distribution and $q$ is the density of the target distribution $\nu$.
Then, it is straightforward to show that independent samples $x_{i}$ together with the weights $w_{i}$ follow the target distribution $\nu$ in the sense that for any $\nu$-integrable function $f\colon\mathbb{R}^{d}\to\mathbb{R}$ it holds that
$\frac{\sum_{i=1}^{N}w_{i}f(x_{i})}{\sum_{i=1}^{N}w_{i}}$
is an unbiased estimator of $\int_{\mathbb{R}^{d}}f(x)\mathrm{d}\nu(x)$.
Note that annealed importance sampling is very sensitive with respect to the proposal $p$ which needs to be designed carefully and problem adapted.

#### Rejection Steps

Inspired by annealed importance sampling, we propose to use importance-based rejection steps.
More precisely, let $\mu$ be a proposal distribution where we can sample from with density $p(x)=f(x)/Z_{f}$ and denote by $\nu$ the target distribution with density $q(x)=g(x)/Z_{g}$. In the following, we assume that we have access to the unnormalized densities $f$ and $g$, but not to the normalization constants $Z_{f}$ and $Z_{g}$
Then, for a random variable $X\sim\mu$, we now generate a new random variable $\tilde{X}$ by the following procedure: First, we compute the importance based acceptance probability $\alpha(X)=\min\left\{1,\frac{q(X)}{\tilde{c}p(X)}\right\}=\min\left\{1,\frac{g
(X)}{cf(X)}\right\}$, where $c>0$ is a positive hyperparameter and $\tilde{c}=cZ_{f}/Z_{g}$. Then, we set $\tilde{X}=X$ with probability $\alpha(X)$ and choose $\tilde{X}=X^{\prime}$ otherwise, where $X^{\prime}\sim\mu$ and $X$ are independent.

###### Remark 7 .

This is a one-step approximation of the classical rejection sampling scheme (Von Neumann, ), see also (Andrieu et al., ) for an overview.
More precisely, we arrive at the classical rejection sampling scheme by choosing $\tilde{c}>\sup_{x}q(x)/p(x)$ and redo the procedure when $X$ is rejected instead of choosing $\tilde{X}=X^{\prime}$.

Similarly to importance sampling, the rejection sampling algorithm is highly sensitive towards the proposal distribution $p$. In particular, if $p$ is chosen as a standard normal distribution, it suffers from the curse of dimensionality. We will tackle this problem later in the section by choosing $p$ already close to the target density $q$.

The following proposition describes the density of the distribution $\tilde{\mu}$ of $\tilde{X}$. Moreover, it ensures that the KL divergence to the target distribution decreases. We include the proof in Appendix C.1

###### Proposition 8 .

Let $\tilde{\mu}$ be the distribution of $\tilde{X}$. Then, the following holds true.

- (i)
$\tilde{\mu}$ admits the density $\tilde{p}$ given by
$\tilde{p}(x)=p(x)(\alpha(x)+1-\mathbb{E}[\alpha(X)])$.
In particular, we have $\tilde{p}(x)=\tilde{f}(x)/Z_{\tilde{f}}$ with $\tilde{f}(x)=f(x)(\alpha(x)+1-\mathbb{E}[\alpha(X)])$.
- (ii)
It holds that $\mathrm{KL}(\tilde{\mu},\nu)\leq\mathrm{KL}(\mu,\nu)$.

Note that the value $\mathbb{E}[\alpha(X)]$ can easily be estimated based on samples during the training phase. Indeed, given $N$ iid copies $X_{1},...,X_{N}$ of $X$, we obtain that $\mathbb{E}[\alpha(X)]\approx\frac{1}{N}\sum_{i=1}^{N}\alpha(X_{N})$ is an unbiased estimator fulfilling the error estimate from the following corollary. The proof is a direct consequence of Hoeffding’s inequality and given in Appendix C.2.

###### Corollary 9 .

Let $X_{1},...,X_{N}$ be iid copies of $X$. Then, it holds

$$ $\mathbb{E}\left[\left|\mathbb{E}[\alpha(X)]-\frac{1}{N}\sum_{i=1}^{N}\alpha(X_ {N})\right|\right]\leq\frac{\sqrt{2\pi}}{\sqrt{N}}\in O\left(\frac{1}{\sqrt{N} }\right).$ $$

###### Remark 10 (Choice of $c$ ) .

We choose the hyperparameter $c$ such that a constant ration $r>0$ of the samples will be resampled, i.e., that $\mathbb{E}[\alpha(X)]\approx 1-r$. To this end, we assume that we are given samples $x_{1},...,x_{N}$ from $X$ and approximate

$$ $\mathbb{E}[\alpha(X)]\approx\frac{1}{N}\sum_{i=1}^{N}\alpha(x_{i})=\frac{1}{N} \sum_{i=1}^{N}\min\left\{1,\frac{g(x)}{cp(x)}\right\}.$ $$

Note that the right side of this formula depends monotonically on $c$ such that we can find $c>0$ such that $\mathbb{E}[\alpha(X)]=1-r$ by a bisection search. In our numerical experiments, we set $r=0.2$.

For the parameter selection and application of a rejection step we refer to Algorithms 3 and 4 in Appendix D.1.

###### Remark 7 .

###### Proposition 8 .

###### Corollary 9 .

###### Remark 10 (Choice of c 𝑐 c italic_c ) .

#### Neural JKO Sampling with Importance Correction

Finally, we combine the neural JKO scheme from the previous section with our rejection steps to obtain a sampling algorithm.
More precisely, we start with a simple latent random variable $X_{0}$ following the probability distribution $\mu_{0}$ with known density and where we can sample from.
In our numerical experiments this will be a standard Gaussian distribution.
Now, we iteratively generate random variables $X_{k}$ following the distribution $\mu_{k}$, $k=1,...,K$ by applying either neural JKO steps as described in Algorithm 1 and 2 or importance-based rejection steps as described in Algorithm 3 and 4.
We call the resulting Markov chain $(X_{0},...,X_{K})$ an importance corrected neural JKO model.
During the sampling process we can maintain the density values $p^{k}(x)$ of the density $p^{k}$ of $\mu_{k}$ for the generated samples by the Algorithms 2 and 4.
Vice versa, we can also use Proposition 8 to evaluate the density $p^{k}$ at some arbitrary point $x\in\mathbb{R}^{d}$. We outline this density evaluation process in Appendix D.2.

###### Remark 11 (Runtime Limitations) .

The sampling time of our importance corrected neural JKO sampling depends exponentially on the number of rejection steps since in each rejection step we resample a constant fraction of the samples. However, due to the moderate base of $1+r$ we will see in the numerical part that we are able to perform a significant number of rejection steps in a tractable time.

###### Remark 11 (Runtime Limitations) .

## 5 Numerical Results

We compare our method with classical Monte Carlo samplers like a Metropolis adjusted Langevin sampling (MALA) and Hamiltonian Monte Carlo (HMC), see e.g., Betancourt (); Roberts & Tweedie ().
Additionally, we compare with two recent deep-learning based sampling algorithms, namely denoising diffusion samplers (DDS, Vargas et al., ) and continual repeated annealed flow transport Monte Carlo (CRAFT, Matthews et al., ). We evaluate all methods on a set of common test distributions which is described in detail in Appendix E.1. Moreover, we report the error measures for our importance corrected neural JKO sampler (neural JKO IC). Additionally, we emphasize the importance of the rejection steps by reporting values for the same a neural JKO scheme without rejection steps (neural JKO).

For evaluating the quality of our results, we use two different metrics. First, we evaluate the energy distance (Székely, ). It is a kernel metrics which can be evaluated purely based on two sets of samples from the model and the ground truth. Moreover it encodes the geometry of the space such that a slight perturbation of the samples only leads to a slight change in the energy distance.
Second, we estimate the log-normalizing constant which is equivalent to approximating the reverse KL loss of the model. A higher estimate of the log-normalizing constant corresponds to a smaller KL divergence between generated and target distribution and therefore to a higher similarity of the two measures. Since this requires the density of the model, this approach is not applicable for MALA and HMC.
The results are given in Table 1 and 2. We can see, that importance corrected neural JKO sampling significantly outperforms the comparison for all test distributions. In particular, we observe that for shifted 8 modes, shifted 8 peaky and GMM-$d$ the energy distance between neural JKO IC and ground truth samples is in the same order of magnitude as the energy distance between to different sets of ground truth samples. This implies that the distribution generated by neural JKO IC is indistinguishable from the target distribution in the energy distance. For these examples, the $\log(Z)$ esitmate is sometimes slightly larger than the ground truth, which can be explained by numerical effects, see Remark 16 for a detailed discussion.
We point out to a precise description of the metrics, the implementation details and additional experiments and figures in Appendix E.

**Table 1: Energy distance. We run each method $5$ times and state the average value and corresponding standard deviations. The rightmost column shows the reference sampling error, i.e., the lower bound magnitude of the average energy distance between two different sets of samples drawn from the ground truth. A smaller energy distance indicates a better result.**
|  |  | Sampler |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Distribution |  | MALA | HMC | DDS | CRAFT | Neural JKO | Neural JKO IC (ours) |  | Sampling Error |
| Mustache |  | $$4.6\text{\times}{10}^{-2}$\pm$1.6\text{\times}{10}^{-3}$$ | $$1.7\text{\times}{10}^{-2}$\pm$4.3\text{\times}{10}^{-4}$$ | $$6.9\text{\times}{10}^{-2}$\pm$1.8\text{\times}{10}^{-3}$$ | $$9.2\text{\times}{10}^{-2}$\pm$9.9\text{\times}{10}^{-3}$$ | $$1.8\text{\times}{10}^{-2}$\pm$2.0\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$2.9\text{\times}{10}^{-3}$\pm$4.4\text{\times}{10}^<br>{-4}$$ |  | $8.6\text{\times}{10}^{-5}$ |
| shifted 8 Modes |  | $$5.3\text{\times}{10}^{-3}$\pm$4.9\text{\times}{10}^{-4}$$ | $$4.1\text{\times}{10}^{-5}$\pm$3.3\text{\times}{10}^{-5}$$ | $$1.2\text{\times}{10}^{-2}$\pm$4.1\text{\times}{10}^{-3}$$ | $$5.2\text{\times}{10}^{-2}$\pm$1.1\text{\times}{10}^{-2}$$ | $$1.3\text{\times}{10}^{-1}$\pm$3.8\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$1.2\text{\times}{10}^{-5}$\pm$5.1\text{\times}{10}^<br>{-6}$$ |  | $2.6\text{\times}{10}^{-5}$ |
| shifted 8 Peaky |  | $$1.3\text{\times}{10}^{-1}$\pm$3.2\text{\times}{10}^{-3}$$ | $$1.2\text{\times}{10}^{-1}$\pm$2.4\text{\times}{10}^{-3}$$ | $$1.1\text{\times}{10}^{-2}$\pm$3.4\text{\times}{10}^{-3}$$ | $$5.2\text{\times}{10}^{-2}$\pm$2.2\text{\times}{10}^{-2}$$ | $$1.3\text{\times}{10}^{-1}$\pm$2.2\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$3.4\text{\times}{10}^{-5}$\pm$8.2\text{\times}{10}^<br>{-6}$$ |  | $2.4\text{\times}{10}^{-5}$ |
| Funnel |  | $$1.2\text{\times}{10}^{-1}$\pm$3.1\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$3.1\text{\times}{10}^{-3}$\pm$3.2\text{\times}{10}^<br>{-4}$$ | $$2.6\text{\times}{10}^{-1}$\pm$2.6\text{\times}{10}^{-2}$$ | $$7.4\text{\times}{10}^{-2}$\pm$2.8\text{\times}{10}^{-3}$$ | $$4.6\text{\times}{10}^{-2}$\pm$1.6\text{\times}{10}^{-3}$$ | $$1.4\text{\times}{10}^{-2}$\pm$8.2\text{\times}{10}^{-4}$$ |  | $3.4\text{\times}{10}^{-4}$ |
| GMM-10 |  | $$1.2\text{\times}{10}^{-2}$\pm$5.5\text{\times}{10}^{-3}$$ | $$1.2\text{\times}{10}^{-2}$\pm$5.2\text{\times}{10}^{-3}$$ | $$3.7\text{\times}{10}^{-3}$\pm$1.6\text{\times}{10}^{-3}$$ | $$1.8\text{\times}{10}^{-1}$\pm$6.6\text{\times}{10}^{-2}$$ | $$1.1\text{\times}{10}^{-2}$\pm$5.6\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$5.3\text{\times}{10}^{-5}$\pm$1.7\text{\times}{10}^<br>{-5}$$ |  | $4.6\text{\times}{10}^{-5}$ |
| GMM-20 |  | $$9.1\text{\times}{10}^{-3}$\pm$2.8\text{\times}{10}^{-3}$$ | $$9.1\text{\times}{10}^{-3}$\pm$2.7\text{\times}{10}^{-3}$$ | $$5.0\text{\times}{10}^{-3}$\pm$1.5\text{\times}{10}^{-3}$$ | $$5.4\text{\times}{10}^{-1}$\pm$1.4\text{\times}{10}^{-1}$$ | $$1.0\text{\times}{10}^{-2}$\pm$2.8\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$1.1\text{\times}{10}^{-4}$\pm$3.4\text{\times}{10}^<br>{-5}$$ |  | $6.4\text{\times}{10}^{-5}$ |
| GMM-50 |  | $$2.4\text{\times}{10}^{-2}$\pm$7.5\text{\times}{10}^{-3}$$ | $$2.4\text{\times}{10}^{-2}$\pm$7.5\text{\times}{10}^{-3}$$ | $$2.3\text{\times}{10}^{-2}$\pm$1.1\text{\times}{10}^{-2}$$ | $$1.8\text{\times}{10}^{0}$\pm$1.7\text{\times}{10}^{-1}$$ | $$2.7\text{\times}{10}^{-2}$\pm$7.8\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$1.0\text{\times}{10}^{-4}$\pm$4.6\text{\times}{10}^<br>{-5}$$ |  | $1.1\text{\times}{10}^{-4}$ |
| GMM-100 |  | $$3.6\text{\times}{10}^{-2}$\pm$1.6\text{\times}{10}^{-2}$$ | $$3.7\text{\times}{10}^{-2}$\pm$1.7\text{\times}{10}^{-2}$$ | $$3.9\text{\times}{10}^{-2}$\pm$2.1\text{\times}{10}^{-2}$$ | $$2.8\text{\times}{10}^{1}$\pm$1.0\text{\times}{10}^{-1}$$ | $$4.7\text{\times}{10}^{-2}$\pm$2.2\text{\times}{10}^{-2}$$ | $\pagecolor{green!15!white}$6.0\text{\times}{10}^{-4}$\pm$3.4\text{\times}{10}^<br>{-4}$$ |  | $1.5\text{\times}{10}^{-4}$ |
| GMM-200 |  | $$6.4\text{\times}{10}^{-2}$\pm$2.1\text{\times}{10}^{-2}$$ | $$6.6\text{\times}{10}^{-2}$\pm$1.9\text{\times}{10}^{-2}$$ | $$9.8\text{\times}{10}^{-2}$\pm$3.1\text{\times}{10}^{-2}$$ | $$3.9\text{\times}{10}^{0}$\pm$1.6\text{\times}{10}^{-1}$$ | $$8.9\text{\times}{10}^{-2}$\pm$2.7\text{\times}{10}^{-2}$$ | $\pagecolor{green!15!white}$3.3\text{\times}{10}^{-3}$\pm$1.9\text{\times}{10}^<br>{-3}$$ |  | $2.0\text{\times}{10}^{-4}$ |

**Table 2: Estimated $\log(Z)$. We run each method $5$ times and state the average value and corresponding standard deviations. The DDS values for LGCP are taken from Vargas et al. (). Higher values of $\log(Z)$ estimates correspond to better results.**
|  |  | Sampler |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Distribution |  | DDS | CRAFT | Neural JKO | Neural JKO IC (ours) |  | Ground Truth |
| Mustache |  | $$-1.5\text{\times}{10}^{-1}$\pm$2.7\text{\times}{10}^{-2}$$ | $$-6.5\text{\times}{10}^{-2}$\pm$5.5\text{\times}{10}^{-2}$$ | $$-3.0\text{\times}{10}^{-2}$\pm$2.6\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$-7.3\text{\times}{10}^{-3}$\pm$8.2\text{\times}{10}<br>^{-4}$$ |  | 0 0 |
| shifted 8 Modes |  | $-$5.7\text{\times}{10}^{-2}$\pm$2.0\text{\times}{10}^{-2}$$ | $-$1.2\text{\times}{10}^{-2}$\pm$1.4\text{\times}{10}^{-3}$$ | $-$3.4\text{\times}{10}^{-1}$\pm$3.1\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}+$5.1\text{\times}{10}^{-6}$\pm$2.4\text{\times}{10}<br>^{-3}$$ |  | 0 0 |
| shifted 8 Peaky |  | $-$1.2\text{\times}{10}^{-1}$\pm$2.2\text{\times}{10}^{-2}$$ | $-$1.8\text{\times}{10}^{-3}$\pm$2.6\text{\times}{10}^{-3}$$ | $-$3.5\text{\times}{10}^{-1}$\pm$3.1\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}-$2.1\text{\times}{10}^{-3}$\pm$3.2\text{\times}{10}<br>^{-3}$$ |  | 0 0 |
| Funnel |  | $-$1.8\text{\times}{10}^{-1}$\pm$6.8\text{\times}{10}^{-2}$$ | $-$1.2\text{\times}{10}^{-1}$\pm$7.9\text{\times}{10}^{-3}$$ | $-$1.4\text{\times}{10}^{-1}$\pm$1.6\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}-$7.1\text{\times}{10}^{-3}$\pm$1.9\text{\times}{10}<br>^{-3}$$ |  | 0 0 |
| GMM-10 |  | $$-2.3\text{\times}{10}^{-1}$\pm$1.0\text{\times}{10}^{-1}$$ | $$-8.5\text{\times}{10}^{-1}$\pm$1.7\text{\times}{10}^{-1}$$ | $$-4.3\text{\times}{10}^{-1}$\pm$5.1\text{\times}{10}^{-2}$$ | $\pagecolor{green!15!white}+$3.5\text{\times}{10}^{-3}$\pm$2.0\text{\times}{10}<br>^{-3}$$ |  | 0 0 |
| GMM-20 |  | $$-5.1\text{\times}{10}^{-1}$\pm$6.0\text{\times}{10}^{-2}$$ | $-$1.5\text{\times}{10}^{0}$\pm$1.7\text{\times}{10}^{-1}$$ | $-$6.3\text{\times}{10}^{-1}$\pm$2.7\text{\times}{10}^{-2}$$ | $\pagecolor{green!15!white}+$6.4\text{\times}{10}^{-3}$\pm$3.8\text{\times}{10}<br>^{-3}$$ |  | 0 0 |
| GMM-50 |  | $-$1.3\text{\times}{10}^{0}$\pm$3.3\text{\times}{10}^{-1}$$ | $-$2.3\text{\times}{10}^{0}$\pm$1.5\text{\times}{10}^{-3}$$ | $-$9.3\text{\times}{10}^{-1}$\pm$4.6\text{\times}{10}^{-2}$$ | $\pagecolor{green!15!white}+$1.1\text{\times}{10}^{-2}$\pm$3.9\text{\times}{10}<br>^{-3}$$ |  | 0 0 |
| GMM-100 |  | $-$3.0\text{\times}{10}^{0}$\pm$7.3\text{\times}{10}^{-1}$$ | $-$2.3\text{\times}{10}^{0}$\pm$8.6\text{\times}{10}^{-2}$$ | $-$1.8\text{\times}{10}^{0}$\pm$9.6\text{\times}{10}^{-2}$$ | $\pagecolor{green!15!white}-$3.9\text{\times}{10}^{-2}$\pm$7.8\text{\times}{10}<br>^{-3}$$ |  | 0 0 |
| GMM-200 |  | $-$9.4\text{\times}{10}^{0}$\pm$7.2\text{\times}{10}^{-1}$$ | $-$6.3\text{\times}{10}^{0}$\pm$1.5\text{\times}{10}^{-1}$$ | $-$5.2\text{\times}{10}^{0}$\pm$2.5\text{\times}{10}^{-1}$$ | $\pagecolor{green!15!white}-$5.6\text{\times}{10}^{-2}$\pm$1.3\text{\times}{10}<br>^{-2}$$ |  | 0 0 |
| LGCP |  | $503.0\pm$7.7\text{\times}{10}^{-1}$$ | $507.6\pm$3.2\text{\times}{10}^{-1}$$ | $499.9\pm$1.7\text{\times}{10}^{-1}$$ | $\pagecolor{green!15!white}508.2\pm$1.0\text{\times}{10}^{-1}$$ |  | not available |

## 6 Conclusions

#### Methodology

We proposed a novel and expressive generative method that enables the efficient and accurate sampling from a prescribed unnormalized target density which is empirically confirmed in numerical examples.
To this end, we combine local sampling steps, relying on piecewise geodesic interpolations of the JKO scheme realized by CNFs, and non-local rejection and resampling steps based on importance weights.
Since the proposal of the rejection step is generated by the model itself,
they do not suffer from the curse of dimensionality as opposed to classical variants of rejection sampling.
The proposed approach provides the advantage that we can draw independent samples while correcting imbalanced mode weights, iteratively refine the current approximation and evaluate the density of the generated distribution.
The latter property is a consequence of the density value propagation through CNFs and Proposition 8 and enables possible further post-processing steps that require density evaluations of the approximated sample process.

#### Outlook

Our method allows for the pointwise access to the approximated target density and the log normalization constant. These quantities can be used for the error monitoring during the training and hence provide guidelines for the adaptive design of the emulator in terms of CNF -or rejection/resampling steps and provide a straightforward stopping criterion.

#### Limitations

In the situation, when the emulator is realised through a stack of underlying rejection/resampling steps, the sample generation process time is negatively affected, see Remark 11. In order to resolve the drawback we plan to utilize diffusion models for the sample generation. This is part of ongoing and future work by the authors.

## Acknowledgements

JH acknowledges funding by the EPSRC programme grant The Mathematics of Deep Learning with reference EP/V026259/1.
RG acknowledges support by the DFG MATH+ project AA5-5 (was EF1-25) -
Wasserstein Gradient Flows for Generalised Transport in Bayesian Inversion.

## Appendix A Background on Wasserstein Spaces

We give some more backgrounds on Wasserstein gradient flows extending Section 2. In what follows we refer to $L_{2}(\mu):=L_{2}(\mu,\mathbb{R}^{d})$ as the set of square integrable measurable functions on $\mathbb{R}^{d}$ with respect to a given measure $\mu\in\mathcal{P}(\mathbb{R}^{d})$.

For an absolutely continuous curve $\gamma$, there exists a unique minimal norm solution $v_{t}$ of the continuity equation (2) in the sense that any solution $\tilde{v}$ of (2) fulfills $\|\tilde{v}(\cdot,t)\|_{L_{2}(\gamma(t))}\geq\|v(\cdot,t)\|_{L_{2}(\gamma(t))}$ for almost every $t$. This is the unique solution of (2) such that $v(\cdot,t)$ is contained in the regular tangent space

$$ $\displaystyle{\mathrm{T}}_{\mu}\mathcal{P}_{2}(\mathbb{R}^{d})$ $\displaystyle\coloneqq\overline{\left\{\lambda(T-\text{Id}):(\text{Id},T)_{\#} \mu\in\Gamma^{\text{opt}}(\mu,T_{\#}\mu),\;\lambda>0\right\}}^{L_{2}(\mu)}.$ (13) $$

An absolutely continuous curve $\gamma\colon(0,\infty)\to\mathcal{P}_{2}(\mathbb{R}^{d})$
with velocity field $v_{t}\in\mathrm{T}_{\gamma(t)}\mathcal{P}_{2}(\mathbb{R}^{d})$ is a *Wasserstein gradient flow with respect to* $\mathcal{F}\colon\mathcal{P}_{2}(\mathbb{R}^{d})\to(-\infty,\infty]$
if

$$ $v_{t}\in-\partial\mathcal{F}(\gamma(t)),\quad\text{for a.e. }t>0,$ (14) $$

where $\partial\mathcal{F}(\mu)$ denotes the reduced Fréchet subdiffential at $\mu$ defined as

$$ $\partial\mathcal{F}(\mu)\coloneqq\left\{\xi\in L_{2}(\mu):\mathcal{F}(\nu)- \mathcal{F}(\mu)\geq\inf_{\pi\in\Gamma^{\mathrm{opt}}(\mu,\nu)}\int\limits_{ \mathbb{R}^{d}\times\mathbb{R}^{d}}\langle\xi(x),y-x\rangle\,\mathrm{d}\pi(x,y )+o(W_{2}(\mu,\nu))\;\forall\nu\in\mathcal{P}_{2}(\mathbb{R}^{d})\right\}.$ (15) $$

The norm $\|v_{t}\|_{L_{2}(\gamma(t))}$ of the velocity field of a Wasserstein gradient flow coincides for almost every $t$ with the metric derivative

$$ $|\partial\mathcal{F}|(\mu)=\inf_{\nu\to\mu}\frac{\mathcal{F}(\mu)-\mathcal{F}( \nu)}{W_{2}(\mu,\nu)}.$ $$

For the convergence result from Theorem 3, we need two more definitions.
First, we need some convexity assumption.
For $\lambda\in\mathbb{R}$, $\mathcal{F}\colon\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}\cup\{+\infty\}$ is called
*$\lambda$-convex along geodesics* if, for every
$\mu,\nu\in\operatorname{dom}\mathcal{F}\coloneqq\{\mu\in\mathcal{P}_{2}(
\mathbb{R}^{d}):\mathcal{F}(\mu)<\infty\}$,
there exists at least one geodesics $\gamma\colon[0,1]\to\mathcal{P}_{2}(\mathbb{R}^{d})$
between $\mu$ and $\nu$ such that

$$ $\mathcal{F}(\gamma(t))\leq(1-t)\,\mathcal{F}(\mu)+t\,\mathcal{F}(\nu)-\tfrac{ \lambda}{2}\,t(1-t)\,W_{2}^{2}(\mu,\nu),\qquad t\in[0,1].$ $$

To ensure uniqueness and convergence of the JKO scheme, a slightly stronger condition, namely being
*$\lambda$-convex along generalized geodesics*
will be in general needed.
Based on the set of three-plans with base
$\sigma\in\mathcal{P}_{2}(\mathbb{R}^{d})$ given by

$$ $\Gamma_{\sigma}(\mu,\nu)\coloneqq\bigl{\{}\bm{\alpha}\in\mathcal{P}_{2}( \mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d}):(\pi_{1})_{\#}\bm{ \alpha}=\sigma,(\pi_{2})_{\#}\bm{\alpha}=\mu,(\pi_{3})_{\#}\bm{\alpha}=\nu \bigr{\}},$ $$

the so-called
*generalized geodesics* $\gamma\colon[0,\epsilon]\to\mathcal{P}_{2}(\mathbb{R}^{d})$
joining $\mu$ and $\nu$ (with base $\sigma$) is defined as

$$ $\gamma(t)\coloneqq\bigl{(}(1-\tfrac{t}{\epsilon})\pi_{2}+\tfrac{t}{\epsilon} \pi_{3}\bigr{)}_{\#}\bm{\alpha},\qquad t\in[0,\epsilon],$ (16) $$

where $\bm{\alpha}\in\Gamma_{\sigma}(\mu,\nu)$
with
$(\pi_{1,2})_{\#}\bm{\alpha}\in\Gamma^{\mathrm{opt}}(\sigma,\mu)$
and
$(\pi_{1,3})_{\#}\bm{\alpha}\in\Gamma^{\mathrm{opt}}(\sigma,\nu)$, see Definition 9.2.2 in Ambrosio et al. ().
The plan $\bm{\alpha}$ may be interpreted as transport from $\mu$ to $\nu$ via $\sigma$.
Then
a function $\mathcal{F}\colon\mathcal{P}_{2}(\mathbb{R}^{d})\to(-\infty,\infty]$ is called
*$\lambda$-convex along generalized geodesics* (see Ambrosio et al., , Definition 9.2.4),
if for every $\sigma,\mu,\nu\in\operatorname{dom}\mathcal{F}$,
there exists at least one generalized geodesics $\gamma\colon[0,1]\to\mathcal{P}_{2}(\mathbb{R}^{d})$
related to some $\bm{\alpha}$ in (16) such that

$$ $\mathcal{F}(\gamma(t))\leq(1-t)\,\mathcal{F}(\mu)+t\,\mathcal{F}(\nu)-\tfrac{ \lambda}{2}\,t(1-t)\,W_{\bm{\alpha}}^{2}(\mu,\nu),\qquad t\in[0,1],$ (17) $$

where

$$ $W_{\bm{\alpha}}^{2}(\mu,\nu)\coloneqq\int_{\mathbb{R}^{d}\times\mathbb{R}^{d} \times\mathbb{R}^{d}}\|y-z\|_{2}^{2}\,\mathrm{d}\bm{\alpha}(x,y,z).$ $$

Every function being $\lambda$-convex along generalized geodesics
is also $\lambda$-convex along geodesics
since generalized geodesics with base $\sigma=\mu$ are actual geodesics.
Second, a $\lambda$-convex functional $\mathcal{F}\colon\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}\cup\{+\infty\}$ is called coercive, if there exists some $r>0$ such that

$$ $\inf\{\mathcal{F}(\mu):\mu\in\mathcal{P}_{2}(\mathbb{R}^{d}),\int_{\mathbb{R}^ {d}}\|x\|^{2}\mathrm{d}\mu(x)\leq r\}>-\infty,$ $$

see (Ambrosio et al., , eq. (11.2.1b)).
In particular, any functional which is bounded from below is coercive.

If $\mathcal{F}$ is proper, lower semicontinuous, coercive and $\lambda$-convex along generalized geodesics, one can show that the $\operatorname{prox}_{\tau\mathcal{F}}(\mu)$ is non-empty and unique for $\tau$ small enough (see Ambrosio et al., , page 295).

## Appendix B Proofs and Examples from Section 3

### B.1 Examples fulfilling Assumption 4

Assumption 4 is fulfilled for many important divergences and loss functions $\mathcal{F}$. We list some examples below. While it is straightforward to check that they are proper, lower semicontinuous, coercive and bounded from below, the convexity is non-trivial. However, the conditions under which these functionals are $\lambda$-convex along generalized geodesics are well investigated in (Ambrosio et al., , Section 9.3).
In the following we denote the Lebesque measure as $\lambda$. For a measure $\mu\in\mathcal{P}^{\mathrm{ac}}(\mathbb{R}^{d})$, we denote by $\mathrm{d}\mu/\mathrm{d}\lambda$ the Lebesque density of $\mu$ if it exists.

- -
Let $\nu\in\mathcal{P}_{2}^{\mathrm{ac}}(\mathbb{R}^{d})$ with Lebesque density $q$ and define the *forward Kullback-Leibler (KL) loss function*
$\mathcal{F}(\mu)=\mathrm{KL}(\nu,\mu)\coloneqq\begin{cases}\int_{\mathbb{R}^{d
}}q(x)\log\left(\frac{q(x)}{p(x)}\right)\mathrm{d}x,&\text{if }\exists\mathrm{
d}\mu/\mathrm{d}\lambda=p\text{ and }\mathrm{d}\nu/\mathrm{d}\lambda=q,\\
+\infty&\text{otherwise.}\end{cases}$
By (Ambrosio et al., , Proposition 9.3.9), we obtain that $\mathcal{F}$ fulfills Assumption 4.
- -
We can also derive a functional, be reversing the arguments in the KL divergence. Then, we arrive at the *reverse KL loss function* given by
$\mathcal{F}(\mu)=\mathrm{KL}(\mu,\nu)\coloneqq\begin{cases}\int_{\mathbb{R}^{d
}}p(x)\log\left(\frac{p(x)}{q(x)}\right)\mathrm{d}x,&\text{if }\exists\mathrm{
d}\mu/\mathrm{d}\lambda=p\text{ and }\mathrm{d}\nu/\mathrm{d}\lambda=q,\\
+\infty&\text{otherwise.}\end{cases}$
Given that $q$ is $\lambda$-convex, we obtain that $\mathcal{F}$ fulfills Assumption 4, see (Ambrosio et al., , Proposition 9.3.2).
- -
Finally, we can define $\mathcal{F}$ based on the Jensen-Shannon divergence. This results into the function
$\mathcal{F}(\mu)=\mathrm{JS}(\mu,\nu)\coloneqq\frac{1}{2}\left[\mathrm{KL}
\left(\mu,\tfrac{1}{2}(\mu+\nu)\right)+\mathrm{KL}\left(\nu,\tfrac{1}{2}(\mu+
\nu)\right)\right].$
Assume $\mu$ and $\nu$ admit Lebesque densities $p$ and $q$ respectively. Then, combining the two previous statements, this fulfills Assumption 4 whenever $p$ and $q$ are $\lambda$-convex.

All of these functionals are integrals of a smooth Lagrangian functional, i.e., there exists some smooth $F\colon\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$ such that

$$ $\mathcal{F}(\mu)=\begin{cases}\int_{\mathbb{R}^{d}}F(x,p(x),\nabla p(x)) \mathrm{d}x,&\text{if }\exists\mathrm{d}\mu/\mathrm{d}\lambda=p,\\ \infty&\text{otherwise.}\end{cases}$ (18) $$

In this case, the limit velocity field from Theorem 3 (which appears as a limit in Theorem 6) can be expressed analytically as the gradient of the so-called variational derivative of $\mathcal{F}$, which is given by

$$ $\frac{\delta\mathcal{F}}{\delta\gamma(t)}(x)=-\nabla\cdot\partial_{3}F(x,p(x), \nabla p(x))+\partial_{2}F(x,p(x),\nabla p(x))$ $$

where $\partial_{i}F$ is the derivative of $F$ with respect to the $i$-th argument and $\gamma$ is the Wasserstein gradient flow, see (Ambrosio et al., , Example 11.1.2).
For the above divergence functionals, computing these terms lead to a (weighted) difference of the Stein scores of the input measure $\gamma(t)$ and the target measure $\nu$ which is a nice link to score-based methods. More precisely, denoting the density of $\gamma(t)$ by $p_{t}$, we obtain the following limiting velocity fields.

- -
For the forward KL loss function we have that $F(x,y,z)=q(x)\log\left(\frac{q(x)}{y}\right)$.
Thus, we have that $\frac{\delta\mathcal{F}}{\delta\mu}(x)=-\frac{q(x)}{p(x)}$.
Hence, the velocity field $v(\cdot,t)=\nabla\frac{\delta\mathcal{F}}{\delta\gamma(t)}$ is given by
$v(x,t)=\frac{q(x)}{p_{t}(x)}\frac{\nabla p_{t}(x)}{p_{t}(x)}-\frac{q(x)}{p_{t}
(x)}\frac{\nabla q(x)}{q(x)}=\frac{q(x)}{p_{t}(x)}\left(\nabla\log(p_{t}(x))-
\nabla\log(q(x))\right).$
- -
For the reverse KL loss function the Lagrangian is given by $F(x,y,z)=y\log\left(\frac{y}{q(x)}\right)$.
Thus, we have that
$\frac{\delta\mathcal{F}}{\delta\mu}(x)=\log\left(\frac{p(x)}{q(x)}\right)+1=
\log(p(x))-\log(q(x))+1$.
Hence, the velocity field $v(\cdot,t)=\nabla\frac{\delta\mathcal{F}}{\delta\gamma(t)}$ is given by
$v(x,t)=\nabla(\log(p_{t})(x))-\nabla(\log(q)(x)).$
- -
For the Jensen-Shannon divergence the Lagrangian is given by $F(x,y,z)=\frac{1}{2}\left(y\log\left(\frac{2y}{y+q(x)}\right)+q(x)\left(\frac{
2q(x)}{y+q(x)}\right)\right)$.
Thus, we have that
$\frac{\delta\mathcal{F}}{\delta\mu}(x)=\frac{1}{2}\log\left(\frac{p(x)}{p(x)+q
(x)}\right).$
Hence, the velocity field $v(\cdot,t)=\nabla\frac{\delta\mathcal{F}}{\delta\gamma(t)}$ is given by
$\displaystyle v(x,t)$
$\displaystyle=\frac{1}{2}\frac{p_{t}(x)+q(x)}{p_{t}(x)}\frac{\nabla p_{t}(x)(p
_{t}(x)+q(x))-p_{t}(x)(\nabla p_{t}(x)+\nabla q(x))}{(p_{t}(x)+q(x))^{2}}$
(19)
$\displaystyle=\frac{1}{2}\frac{\nabla p_{t}(x)q(x)-p_{t}(x)\nabla q(x)}{(p_{t}
(x)+q(x))p_{t}(x)}$
(20)
$\displaystyle=\frac{q(x)}{2(p_{t}(x)+q(x))}\left[\frac{\nabla p_{t}(x)}{p_{t}(
x)}-\frac{\nabla q(x)}{q(x)}\right]$
(21)
$\displaystyle=\frac{q(x)}{2(p_{t}(x)+q(x))}\left[\nabla(\log(p_{t})(x))-\nabla
(\log(q)(x))\right].$
(22)

Note that computing the score $\nabla\log(p_{t})$ of the current approximation is usually intractable, such that these limits cannot be inserted into a neural ODE directly.

### B.2 Proof of Corollary 5

Using similar arguments as in Altekrüger et al. (); Mokrov et al. (); Onken et al. (); Xu et al. (), let $v\colon\mathbb{R}^{d}\times[0,\tau]\to\mathbb{R}^{d}$ such that the solution of

$$ $\dot{z}(x,t)=v(z(x,t),t),\quad z(x,0)=x,$ $$

fulfills $z(\cdot,\tau)_{\#}\mu_{\tau}^{k}=\mu_{\tau}^{k+1}$.
Since $(v_{\tau,k},z_{\tau,k})$ is a minimizer of (10), we obtain that

$$ $\displaystyle\phantom{=}\,\,\,\frac{1}{2}\int_{0}^{\tau}\int_{\mathbb{R}^{d}} \|v_{\tau,k}(z_{\tau,k}(x,t),t)\|^{2}\mathrm{d}\mu_{\tau}^{k}(x)\mathrm{d}t+ \mathcal{F}(z_{\tau,k}(\cdot,\tau)_{\#}\mu_{\tau}^{k})$ (23) $\displaystyle\leq\frac{1}{2}\int_{0}^{\tau}\int_{\mathbb{R}^{d}}\|v(z(x,t),t) \|^{2}\mathrm{d}\mu_{\tau}^{k}(x)\mathrm{d}t+\mathcal{F}(z(\cdot,\tau)_{\#}\mu _{\tau}^{k}).$ (24) $$

Observing that $\mathcal{F}(z_{\tau,k}(\cdot,\tau)_{\#}\mu_{\tau}^{k})=\mathcal{F}(z(\cdot,
\tau)_{\#}\mu_{\tau}^{k})=\mathcal{F}(\mu_{\tau}^{k+1})$, we obtain that

$$ $\tau\int_{0}^{\tau}\int_{\mathbb{R}^{d}}\|v_{\tau,k}(z_{\tau,k}(x,t),t)\|^{2} \mathrm{d}\mu_{\tau}^{k}(x)\leq\tau\int_{0}^{\tau}\int_{\mathbb{R}^{d}}\|v(z(x ,t),t)\|^{2}\mathrm{d}\mu_{\tau}^{k}(x)\mathrm{d}t.$ $$

Since $v$ was chose arbitrary, we obtain that $v_{\tau,k}$ is the optimal velocity field from the theorem of Benamou-Brenier, which now directly implies part (ii) and (iii). $\Box$

### B.3 Proof of Theorem 6

In order to prove convergence of the velocity fields $v_{\tau}$, we first introduce some notations. To this end, let $\mathcal{T}_{\tau}^{k}$ be the optimal transport maps between $\mu_{\tau}^{k}$ and $\mu_{\tau}^{k+1}$ and define by $v_{\tau}^{k}=(\mathcal{T}_{\tau}^{k}-I)/\tau$ the corresponding discrete velocity fields.
Then, the velocity fields $v_{\tau}$ can be expressed as

$$ $v_{\tau}(x,k\tau+t\tau)=v_{\tau}^{k}(((1-t)I+t\mathcal{T}_{\tau}^{k})^{-1}(x)).$ (25) $$

Further, we denote the piece-wise constant concatenation of the discrete velocity fields by

$$ $\tilde{v}_{\tau}(x,k\tau+t\tau)=v_{\tau}^{k},\quad t\in(0,1).$ $$

Note that for any $\tau$ it holds that $v_{\tau}\in L_{2}(\gamma_{\tau},\mathbb{R}^{d}\times[0,T])$ and $\tilde{v}_{\tau}\in L_{2}(\tilde{\gamma}_{\tau},\mathbb{R}^{d}\times[0,T])$.
To derive limits of these velocity fields, we recall the notion of convergence from (Ambrosio et al., , Definition 5.4.3) allowing that the iterates are not defined on the same space. In this paper, we stick to square integrable measurable functions defined on finite dimensional domains, which slightly simplifies the definition.

###### Definition 12 .

Let $\Omega\subseteq\mathbb{R}^{d}$ be a measurable domain, assume that $\mu_{k}\in\mathcal{P}_{2}(\Omega)$ converges weakly to $\mu\in\mathcal{P}_{2}(\Omega)$ and let $f_{k}\in L_{2}(\mu_{k},\Omega)$ and $f\in L_{2}(\mu,\Omega)$. Then, we say that $f_{k}$ converges weakly to $f$, if

$$ $\int_{\Omega}\langle\phi(x),f_{k}(x)\rangle\mathrm{d}\mu_{k}(x)\to\int_{\Omega }\langle\phi(x),f(x)\rangle\mathrm{d}\mu(x),\quad\text{as }k\to\infty$ $$

for all test functions $\phi\in C_{c}^{\infty}(\Omega)$.
We say that $f_{k}$ converges strongly to $f$ if

$$ $\limsup\limits_{k\to\infty}\|f_{k}\|_{L_{2}(\mu_{k},\Omega)}\leq\|f\|_{L_{2}( \mu,\Omega)}.$ (26) $$

Note that (Ambrosio et al., , Theorem 5.4.4 (iii)) implies that formula (26) is fulfilled with equality for any strongly convergent sequence $f_{k}$ to $f$.
Moreover it is known from the literature that subsequences of the piece-wise constant velocity admit weak limits.

###### Theorem 13 ( Ambrosio et al., 2005 , Theorem 11.1.6 ) .

Suppose that Assumption 4 is fulfilled for $\mathcal{F}$. Then, for any $\mu^{0}\in\operatorname{dom}(\mathcal{F})$ and any sequence $(\tau_{l})_{l}\subset(0,\infty)$, there exists a subsequence (again denoted by $(\tau_{l})_{l}$) such that

- -
The piece-wise constant curve $\tilde{\gamma}_{\tau_{l}}(t)$ narrowly converges to some limit curve $\hat{\gamma}(t)$ for all $t\in[0,\infty)$.
- -
The velocity field $\tilde{v}_{\tau_{l}}\in L_{2}(\tilde{\gamma}_{\tau_{l}}(t),\mathbb{R}^{d}
\times[0,T])$ weakly converges to some limit $\hat{v}\in L_{2}(\hat{\gamma}(t),\mathbb{R}^{d}\times[0,T])$ according to Definition 12 for any $T>0$.
- -
The limit $\hat{v}$ fulfills the continuity equation with respect to $\hat{\gamma}$, i.e.,
$\partial_{t}\hat{\gamma}(t)+\nabla\cdot(\hat{v}_{\tau_{l}}(\cdot,t)\hat{\gamma
}(t))=0.$

In order to show the desired result, there remain the following questions, which we answer in the rest of this section:

- -
Does Theorem 13 also hold for the velocity fields $v_{\tau_{l}}$ defined in (25) belonging to the geodesic interpolations?
- -
Can we show strong convergence for the whole sequence $v_{\tau_{l}}$?
- -
Does the limit $\hat{v}$ have the norm-minimizing property that $\|\hat{v}(\cdot,t)\|_{L_{2}(\hat{\gamma}(t))}=|\partial\mathcal{F}|(\hat{
\gamma}(t))$?

To address the first of these questions, we show that weak limits of $\tilde{v}_{\tau}$ and $v_{\tau}$ coincide. The proof is a straightforward computation. A similar statement in a slightly different setting was proven in (Santambrogio, , Section 8.3).

###### Lemma 14 .

Suppose that Assumption 4 is fulfilled and let $(\tau_{l})_{l}\subseteq(0,\infty)$ be a sequence with $\tau_{l}\to 0$ as $l\to\infty$ such that $\tilde{v}_{\tau_{l}}\in L_{2}(\tilde{\gamma}_{\tau_{l}},\mathbb{R}^{d}\times[0
,T])$ converges weakly to some $\hat{v}\in L_{2}(\hat{\gamma},\mathbb{R}^{d}\times[0,T])$. Then, also $v_{\tau_{l}}\in L_{2}(\gamma_{\tau_{l}},\mathbb{R}^{d}\times[0,T])$ converges weakly to $\hat{v}$.

###### Proof.

Let $\phi\in C_{c}^{\infty}(\mathbb{R}^{d}\times[0,T])$. In particular, $\phi$ is Lipschitz continuous with some Lipschitz constant $L<\infty$.
Then, it holds that

$$ $\displaystyle\phantom{\leq}\hskip 3.80005pt\left|\int_{\mathbb{R}^{d}\times[0, T]}\langle\phi,v_{\tau_{l}}\rangle\mathrm{d}\gamma_{\tau_{l}}-\int_{\mathbb{R} ^{d}\times[0,T]}\langle\phi,v\rangle\mathrm{d}\hat{\gamma}\right|$ $\displaystyle\leq\left|\int_{\mathbb{R}^{d}\times[0,T]}\langle\phi,v_{\tau_{l} }\rangle\mathrm{d}\gamma_{\tau_{l}}-\int_{\mathbb{R}^{d}\times[0,T]}\langle \phi,\tilde{v}_{\tau_{l}}\rangle\mathrm{d}\tilde{\gamma}_{\tau_{l}}\right|+ \left|\int_{\mathbb{R}^{d}\times[0,T]}\langle\phi,\tilde{v}_{\tau_{l}}\rangle \mathrm{d}\tilde{\gamma}_{\tau_{l}}-\int_{\mathbb{R}^{d}\times[0,T]}\langle \phi,v\rangle\mathrm{d}\hat{\gamma}\right|.$ $$

Since $\tilde{v}_{\tau_{l}}$ converges weakly to $\hat{v}$, the second term converges to zero as $l\to\infty$. Thus in order to show that also $v_{\tau_{l}}$ converges weakly to $\hat{v}$, it remains to show that also the first term converges to zero.
By using the notations $T_{l}=\min\{k\tau_{l}:k\tau_{l}\geq T,k\in\mathbb{Z}_{\geq 0}\}$ and $K_{l}=\lceil\tfrac{T}{\tau_{l}}\rceil=\tfrac{T_{l}}{\tau_{l}}$, we can estimate

$$ $\displaystyle\phantom{=}\hskip 3.80005pt\left|\int_{\mathbb{R}^{d}\times[0,T]} \langle\phi,v_{\tau_{l}}\rangle\mathrm{d}\gamma_{\tau_{l}}-\int_{\mathbb{R}^{d }\times[0,T]}\langle\phi,\tilde{v}_{\tau_{l}}\rangle\mathrm{d}\tilde{\gamma}_{ \tau_{l}}\right|$ $\displaystyle=\left|\int_{0}^{T}\left(\int_{\mathbb{R}^{d}}\langle\phi(x,t),v_ {\tau_{l}}(x,t)\rangle\mathrm{d}\gamma_{\tau_{l}}(t)(x)-\int_{\mathbb{R}^{d}} \langle\phi(x,t),\tilde{v}_{\tau_{l}}(x,t)\rangle\mathrm{d}\tilde{\gamma}_{ \tau_{l}}(t)(x)\right)\mathrm{d}t\right|$ $\displaystyle\leq\int_{0}^{T_{l}}\left|\int_{\mathbb{R}^{d}}\langle\phi(x,t),v _{\tau_{l}}(x,t)\rangle\mathrm{d}\gamma_{\tau_{l}}(t)(x)-\int_{\mathbb{R}^{d}} \langle\phi(x,t),\tilde{v}_{\tau_{l}}(x,t)\rangle\mathrm{d}\tilde{\gamma}_{ \tau_{l}}(t)(x)\right|\mathrm{d}t$ $\displaystyle\leq\sum_{k=0}^{K_{l}-1}\int_{k\tau_{l}}^{(k+1)\tau_{l}}\left| \int_{\mathbb{R}^{d}}\langle\phi(x,t),v_{\tau_{l}}(x,t)\rangle\mathrm{d}\gamma _{\tau_{l}}(t)(x)-\int_{\mathbb{R}^{d}}\langle\phi(x,t),\tilde{v}_{\tau_{l}}(x ,t)\rangle\mathrm{d}\tilde{\gamma}_{\tau_{l}}(t)(x)\right|\mathrm{d}t$ $\displaystyle=\tau_{l}\sum_{k=0}^{K_{l}-1}\int_{0}^{1}\left|\int_{\mathbb{R}^{ d}}\langle\phi(x,k\tau_{l}+t\tau_{l}),v_{\tau_{l}}(x,k\tau_{l}+t\tau_{l}) \rangle\mathrm{d}\gamma_{\tau_{l}}(k\tau_{l}+t\tau_{l})(x)\right.$ $\displaystyle\qquad\qquad\quad\left.-\int_{\mathbb{R}^{d}}\langle\phi(x,k\tau_ {l}+t\tau_{l}),\tilde{v}_{\tau_{l}}(x,k\tau_{l}+t\tau_{l})\rangle\mathrm{d} \tilde{\gamma}_{\tau_{l}}(k\tau_{l}+t\tau_{l})(x)\right|\mathrm{d}t.$ $$

By denoting with $\mathcal{T}_{\tau_{l}}^{k}$ the optimal transport map from $\mu_{\tau_{l}}^{k}$ to $\mu_{\tau_{l}}^{k+1}$, we have for $t\in(0,1)$ that

$$ $\gamma_{\tau_{l}}(k\tau_{l}+t\tau_{l})=((1-t)I+t\mathcal{T}_{\tau_{l}}^{k})_{ \#}\mu_{\tau_{l}}^{k},\quad\tilde{\gamma}_{\tau_{l}}(k\tau_{l}+t\tau_{l})=\mu_ {\tau_{l}}^{k},$ $$

and the velocity fields satisfy

$$ $v_{\tau_{l}}(x,k\tau_{l}+t\tau_{l})=v_{\tau_{l}}^{k}(((1-t)I+t\mathcal{T}_{ \tau_{l}}^{k})^{-1}(x)),\quad\tilde{v}_{\tau_{l}}(x,k\tau_{l}+t\tau_{l})=v_{ \tau_{l}}^{k}(x).$ $$

Then, the above term becomes

$$ $\displaystyle\quad\tau_{l}\sum_{k=0}^{K_{l}-1}\int_{0}^{1}\left|\int_{\mathbb{ R}^{d}}\langle\phi(x,k\tau_{l}+t\tau_{l}),v_{\tau_{l}}^{k}(((1-t)I+t\mathcal{T }_{\tau_{l}}^{k})^{-1}(x))\rangle\mathrm{d}((1-t)I+t\mathcal{T}_{\tau_{l}}^{k} )_{\#}\mu_{\tau_{l}}^{k}(x)\right.$ $\displaystyle\qquad\qquad\quad\left.-\int_{\mathbb{R}^{d}}\langle\phi(x,k\tau_ {l}+t\tau_{l}),v_{\tau_{l}}^{k}(x)\rangle\mathrm{d}\mu_{\tau_{l}}^{k}(x)\right |\mathrm{d}t$ $\displaystyle=\tau_{l}\sum_{k=0}^{K_{l}-1}\int_{0}^{1}\left|\int_{\mathbb{R}^{ d}}\langle\phi((1-t)x+t\mathcal{T}_{\tau_{l}}^{k}(x)),k\tau_{l}+t\tau_{l}),v_{ \tau_{l}}^{k}(x)\rangle\mathrm{d}\mu_{\tau_{l}}^{k}(x)\right.$ $\displaystyle\qquad\qquad\quad\left.-\int_{\mathbb{R}^{d}}\langle\phi(x,k\tau_ {l}+t\tau_{l}),v_{\tau_{l}}^{k}(x)\rangle\mathrm{d}\mu_{\tau_{l}}^{k}(x)\right |\mathrm{d}t$ $\displaystyle\leq\tau_{l}\sum_{k=0}^{K_{l}-1}\int_{0}^{1}\int_{\mathbb{R}^{d}} \left|\langle\phi((1-t)x+t\mathcal{T}_{\tau_{l}}^{k}(x)),k\tau_{l}+t\tau_{l})- \phi(x,k\tau_{l}+t\tau_{l}),v_{\tau_{l}}^{k}(x)\rangle\right|\mathrm{d}\mu_{ \tau_{l}}^{k}(x)\mathrm{d}t.$ $$

Using Hölders’ inequality and we obtain that this is smaller or equal than

$$ $\displaystyle\tau_{l}\sum_{k=0}^{K_{l}-1}\int_{0}^{1}\left(\int_{\mathbb{R}^{d }}\|\phi((1-t)x+t\mathcal{T}_{\tau_{l}}^{k}(x)),k\tau_{l}+t\tau_{l})-\phi(x,k \tau_{l}+t\tau_{l})\|^{2}\mathrm{d}\mu_{\tau_{l}}^{k}(x)\int_{\mathbb{R}^{d}} \|v_{\tau_{l}}^{k}(x)\|^{2}\mathrm{d}\mu_{\tau_{l}}^{k}(x)\right)^{1/2}\mathrm {d}t$ $$

By the Lipschitz continuity of $\phi$ and inserting the definition of $v_{\tau_{l}}^{k}=\frac{\mathcal{T}_{\tau_{l}}^{k}-I}{\tau_{l}}$ this is smaller or equal than

$$ $\displaystyle\phantom{=}\hskip 3.80005pt\tau_{l}\sum_{k=0}^{K_{l}-1}\int_{0}^{ 1}\left(t^{2}\frac{L^{2}}{\tau_{l}^{2}}\int_{\mathbb{R}^{d}}\|\mathcal{T}_{ \tau_{l}}^{k}(x)-x\|^{2}\mathrm{d}\mu_{\tau_{l}}^{k}(x)\int_{\mathbb{R}^{d}}\| \mathcal{T}_{\tau_{l}}^{k}(x)-x\|^{2}\mathrm{d}\mu_{\tau_{l}}^{k}(x)\right)^{1 /2}\mathrm{d}t$ (27) $\displaystyle=L\sum_{k=0}^{K_{l}-1}\left(\int_{0}^{1}t\mathrm{d}t\right)\left( \int_{\mathbb{R}^{d}}\|\mathcal{T}_{\tau_{l}}^{k}(x)-x\|^{2}\mathrm{d}\mu_{ \tau_{l}}^{k}(x)\right)$ (28) $\displaystyle=\frac{L}{2}\sum_{k=0}^{K_{l}-1}W_{2}^{2}(\mu_{\tau_{l}}^{k},\mu_ {\tau_{l}}^{k+1})$ (29) $\displaystyle\leq\frac{L}{2}\sum_{k=0}^{\infty}W_{2}^{2}(\mu_{\tau_{l}}^{k}, \mu_{\tau_{l}}^{k+1})$ (30) $$

Finally, we have by the definition of the minimizing movements scheme that

$$ $\frac{1}{2\tau_{l}}W_{2}^{2}(\mu_{\tau_{l}}^{k},\mu_{\tau_{l}}^{k+1})\leq \mathcal{F}(\mu_{\tau_{l}}^{k})-\mathcal{F}(\mu_{\tau_{l}}^{k+1}).$ $$

Summing up for $k=0,1,...$ we finally arrive at the bound

$$ $\left|\int_{\mathbb{R}^{d}\times[0,T]}\langle\phi,v_{\tau_{l}}\rangle\mathrm{d }\gamma_{\tau_{l}}-\int_{\mathbb{R}^{d}\times[0,T]}\langle\phi,v\rangle\mathrm {d}\hat{\gamma}\right|\leq\tau_{l}L\left(\mathcal{F}(\mu^{0})-\inf_{\mu\in \mathcal{P}_{2}(\mathbb{R}^{d})}\mathcal{F}(\mu)\right).$ $$

Since $\mathcal{F}$ is bounded from below the upper bound converges to zero as $l\to\infty$.
This concludes the proof.
∎

Finally, we employ the previous results to show Theorem 6 from the main part of the paper. That is, we show that for any $(\tau_{l})_{l}\subseteq(0,\infty)$ with $\tau_{l}\to 0$ the whole $v_{\tau_{l}}$ converges strongly to $v$.

###### Theorem 15 (Theorem 6 ) .

Suppose that Assumption 4 is fulfilled and
let $(\tau_{l})_{l}\subseteq(0,\infty)$ with $\tau_{l}\to 0$. Then, $(v_{\tau_{l}})_{l}$ converges strongly to the velocity field $\hat{v}\in L_{2}(\gamma,\mathbb{R}^{d}\times[0,T])$ of the Wasserstein gradient flow $\gamma\colon(0,\infty)\to\mathcal{P}_{2}(\mathbb{R}^{d})$ of $\mathcal{F}$ starting in $\mu^{0}$.

###### Proof.

We show that any subsequence of $\tau_{l}$ admits a subsequence converging strongly to $\hat{v}$. Using the sub-subsequence criterion this yields the claim.

By Theorem 13 and Lemma 14 we know that any subsequence of $\tau_{l}$ admits a weakly convergent subsequence. In an abuse of notations, we denote it again by $\tau_{l}$ and its limit by $\tilde{v}$. Then, we prove that the convergence is indeed strong and that $\tilde{v}=\hat{v}$.

Step 1: Bounding $\limsup\limits_{l\to\infty}\int_{0}^{T}\|v_{\tau_{l}}(\cdot,t)\|^{2}_{L_{2}(
\gamma_{\tau_{l}}(t),\mathbb{R}^{d})}\mathrm{d}t$ from above.
Since $\lambda$-convexity with $\lambda\geq 0$ implies $\lambda$-convexity with $\lambda=-1$, we can assume without loss of generality that $\lambda<0$.
Then, by (Ambrosio et al., , Lemma 9.2.7, Theorem 4.0.9), we know that for any $\tau>0$ it holds

$$ $W_{2}(\tilde{\gamma}_{\tau}(t),\gamma(t))\leq\tau C(\tau,t),\qquad C(\tau,t)= \frac{(1+2|\lambda|t_{\tau})|\partial\mathcal{F}|(\mu_{0})}{\sqrt{2}(1+\lambda \tau)}\exp\left(-\frac{\log(1+\lambda\tau)}{\tau}t\right),$ $$

where $t_{\tau}=\min\{k\tau:k\tau\geq t,k\in\mathbb{N}_{0}\}$.
For simplicity, we use the notations $\tau_{\max}=\max\{\tau_{l}:l\in\mathbb{N}_{0}\}$, $K_{l}=\lceil\frac{T}{\tau_{l}}\rceil$ and $T_{\max}=T+\tau_{\max}$.
Since $\lambda<0$, we have that $C(\tau,t)\leq C(\tau,T_{\max})\eqqcolon C(\tau)$ for all $t\in[0,T_{\max}]$. Moreover, we have that $C(\tau)\to\frac{(1+2|\lambda|T_{\mathrm{max}}|\partial\mathcal{F}|(\mu_{0})}{
\sqrt{2}}\exp(-\lambda T_{\max})$ as $\tau\to 0$, such that the sequence $(C(\tau_{l}))_{l}$ is bounded. In particular, there exists a $C>0$ such that

$$ $W_{2}(\tilde{\gamma}_{\tau_{l}}(t),\gamma(t))\leq\tau_{l}C,\quad\text{for all} \quad t\in[0,T_{\max}].$ $$

Inserting $t=k\tau_{l}$ for $k=0,...,K_{l}$ gives

$$ $W_{2}(\mu_{\tau_{l}}^{k},\gamma(k\tau_{l}))\leq\tau_{l}C,\quad\text{for all} \quad k=0,...,K_{l}.$ (31) $$

Now, we can conclude by the theorem of Benamou-Brenier that

$$ $\displaystyle\int_{0}^{T}\|v_{\tau_{l}}(\cdot,t)\|^{2}_{L_{2}(\gamma_{\tau_{l} }(t),\mathbb{R}^{d})}\mathrm{d}t$ $\displaystyle\leq\int_{0}^{K_{l}\tau_{l}}\|v_{\tau_{l}}(\cdot,t)\|^{2}_{L_{2}( \gamma_{\tau_{l}}(t),\mathbb{R}^{d})}\mathrm{d}t$ (32) $\displaystyle=\sum_{k=0}^{K_{l}-1}\int_{k\tau_{l}}^{(k+1)\tau_{l}}\|v_{\tau_{l }}(\cdot,t)\|^{2}_{L_{2}(\gamma_{\tau_{l}}(t),\mathbb{R}^{d})}\mathrm{d}t$ (33) $\displaystyle=\sum_{k=1}^{K_{l}-1}W_{2}^{2}(\mu_{\tau_{l}}^{k},\mu_{\tau_{l}}^ {k+1}).$ (34) $$

Now applying the triangular inequality for any $k=1,\ldots,K_{l}-1$

$$ $W_{2}^{2}(\mu_{\tau_{l}}^{k},\mu_{\tau_{l}}^{k+1})\leq\left(W_{2}(\mu_{\tau_{l }}^{k},\gamma(\tau_{l}k))+W_{2}(\gamma(\tau_{l}k),\gamma(\tau_{l}(k+1))+W_{2}( \gamma(\tau_{l}(k+1)),\mu_{\tau_{l}}^{k+1})\right)^{2}$ $$

and the estimate from (31) yields

$$ $\displaystyle\int_{0}^{T}\|v_{\tau_{l}}(\cdot,t)\|^{2}_{L_{2}(\gamma_{\tau_{l} }(t),\mathbb{R}^{d})}\mathrm{d}t$ $\displaystyle\leq\sum_{k=0}^{K_{l}-1}\left(2\tau C+W_{2}(\gamma(\tau_{l}k), \gamma(\tau_{l}(k+1))\right)^{2}$ $$

By Jensens’ inequality the right hand side can be bounded from above by

$$ $4K_{l}\tau^{2}C^{2}+2\tau K_{l}C\left(\sum_{k=0}^{K_{l}-1}\frac{1}{K_{l}}W_{2} ^{2}(\gamma(\tau_{l}k),\gamma(\tau_{l}(k+1)))\right)^{1/2}+\sum_{k=0}^{K_{l}-1 }W_{2}^{2}(\gamma(\tau_{l}k),\gamma(\tau_{l}(k+1)))$ $$

Finally, again Benamou-Brenier gives that $W_{2}^{2}(\gamma(\tau_{l}k),\gamma(\tau_{l}(k+1)))\leq\int_{k\tau_{l}}^{(k+1)
\tau_{l}}\|\hat{v}(\cdot,t)\|^{2}_{L_{2}(\gamma(t))}\mathrm{d}t$ such that the above formula is smaller or equal than

$$ $4K_{l}\tau_{l}^{2}C^{2}+2\tau_{l}\sqrt{K_{l}}C\left(\int_{0}^{K_{l}\tau_{l}}\| \hat{v}(\cdot,t)\|^{2}_{L_{2}(\gamma(t))}\mathrm{d}t\right)^{1/2}+\int_{0}^{K_ {l}\tau_{l}}\|\hat{v}(\cdot,t)\|^{2}_{L_{2}(\gamma(t))}\mathrm{d}t.$ $$

Since $K_{l}\tau_{l}\to T$ and $\tau_{l}\to 0$ as $l\to\infty$ this converges to $\int_{0}^{T}\|\hat{v}(\cdot,t)\|^{2}_{L_{2}(\gamma(t))}\mathrm{d}t$ such that we can conclude

$$ $\limsup\limits_{l\to\infty}\int_{0}^{T}\|v_{\tau_{l}}(\cdot,t)\|^{2}_{L_{2}( \gamma_{\tau_{l}}(t),\mathbb{R}^{d})}\mathrm{d}t\leq\int_{0}^{T}\|\hat{v}( \cdot,t)\|^{2}_{L_{2}(\gamma(t))}\mathrm{d}t.$ $$

Step 2: Strong convergence.
By (Ambrosio et al., , Theorem 8.3.1, Proposition 8.4.5) we know that for any $v$ fulfilling the continuity equation

$$ $\partial_{t}\gamma(t)+\nabla\cdot(v(\cdot,t)\gamma(t))=0,$ $$

it holds that $\|v(\cdot,t)\|_{L_{2}(\gamma(t),\mathbb{R}^{d}}\geq\|\hat{v}(\cdot,t)\|_{L_{2}
(\gamma(t),\mathbb{R}^{d}}$.
Since $\tilde{v}$ fulfills the continuity equation by Theorem 13, this implies that

$$ $\int_{0}^{T}\|\tilde{v}(\cdot,t)\|^{2}_{L_{2}(\gamma(t),\mathbb{R}^{d}}\mathrm {d}t\geq\int_{0}^{T}\|\hat{v}(\cdot,t)\|^{2}_{L_{2}(\gamma(t),\mathbb{R}^{d}} \mathrm{d}t=\limsup\limits_{l\to\infty}\int_{0}^{T}\|v_{\tau_{l}}(\cdot,t)\|^{ 2}_{L_{2}(\gamma_{\tau_{l}}(t),\mathbb{R}^{d})}\mathrm{d}t.$ $$

In particular, $v_{\tau_{l}}\in L_{2}(\gamma_{\tau_{l}},\mathbb{R}^{d}\times[0,T])$ converges strongly to $\tilde{v}$ such that by (Ambrosio et al., , Theorem 5.4.4 (iii)) it holds equality in the above equation, i.e.,

$$ $\int_{0}^{T}\|\tilde{v}(\cdot,t)\|^{2}_{L_{2}(\gamma(t),\mathbb{R}^{d}}\mathrm {d}t=\int_{0}^{T}\|\hat{v}(\cdot,t)\|^{2}_{L_{2}(\gamma(t),\mathbb{R}^{d}} \mathrm{d}t$ $$

Using again (Ambrosio et al., , Theorem 8.3.1, Proposition 8.4.5) this implies that $\tilde{v}=\hat{v}$.
∎

###### Definition 12 .

###### Theorem 13 ( Ambrosio et al., 2005 , Theorem 11.1.6 ) .

###### Lemma 14 .

###### Proof.

###### Theorem 15 (Theorem 6 ) .

###### Proof.

## Appendix C Proofs from Section 4

### C.1 Proof of Proposition 8

- (i)
By the law of total probability, we have
$\tilde{p}(x)=P(\tilde{X}=x)=P(\tilde{X}=x\text{ and }X\text{ was accepted})+P(
\tilde{X}=x\text{ and }X\text{ was rejected}).$
Since it holds $\tilde{X}=X$ if $X$ and $\tilde{X}=X^{\prime}$ if $X$ is rejected this can be reformulated as
$\displaystyle\quad P(X=x\text{ and }X\text{ was accepted})+P(X^{\prime}=x\text
{ and }X\text{ was rejected})$
(35)
$\displaystyle=P(X=x)P(X\text{ was accepted}|X=x)+P(X^{\prime}=x)P(X\text{ was
rejected}),$
(36)
where we used the definition of conditional probabilities and the fact that $X$ and $X^{\prime}$ are independent.
Because $X,X^{\prime}\sim\mu$, we now have $P(X=x)=P(X^{\prime}=x)=p(x)$ and, by definition, that $P(X\text{ was accepted}|X=x)=\alpha(x)$.
Finally, it holds that
$\displaystyle P(X\text{ was rejected})$
$\displaystyle=\int_{\mathbb{R}^{d}}P(X\text{ was rejected}|X=x)p(x)\mathrm{d}x$
(37)
$\displaystyle=\int_{\mathbb{R}^{d}}(1-\alpha(x))p(x)\mathrm{d}x=1-\mathbb{E}[
\alpha(X)].$
(38)
Inserting these terms in (36) yields the claim.
- (ii)
Note, since $X\sim p$, we have
$\int_{\mathbb{R}^{d}}\frac{p(x)\alpha(x)}{\mathbb{E}[\alpha(X)]}\mathrm{d}x=
\int_{\mathbb{R}^{d}}\frac{p(x)\alpha(x)}{\int_{\mathbb{R}^{d}}p(y)\alpha(y)
\mathrm{d}y}\mathrm{d}x=1,$
in particular $\frac{p(x)\alpha(x)}{\mathbb{E}[\alpha(X)]}$ defines a density. Now let $\eta\in\mathcal{P}_{2}(\mathbb{R}^{d})$ be the corresponding probability measure.
Then, it holds by part (i) that $\tilde{\mu}=\mathbb{E}[\alpha(X)]\eta+(1-\mathbb{E}[\alpha(X)])\mu$.
We will show that $\mathrm{KL}(\eta,\nu)\leq\mathrm{KL}(\mu,\nu)$.
Due to the convexity of the KL divergence in the linear space of measures this implies the claim.
We denote $Z=\mathbb{E}[\alpha(X)]=\int_{\mathbb{R}^{d}}\min(p(x),\tfrac{q(x)}{\tilde{c}}
)\mathrm{d}x$. Then it holds
$\displaystyle\mathrm{KL}(\eta,\nu)$
$\displaystyle=\int_{\mathbb{R}^{d}}\frac{p(x)\alpha(x)}{Z}\log\left(\frac{p(x)
\alpha(x)}{Zq(x)}\right)\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\frac{\min(p(x),\frac{q(x)}{\tilde{c}})}{Z}
\log\left(\frac{\tilde{c}p(x)\alpha(x)}{q(x)}\right)\mathrm{d}x-\log(Z\tilde{c
})\int_{\mathbb{R}^{d}}\frac{\min(p(x),\frac{q(x)}{\tilde{c}})}{Z}\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\frac{\min(p(x),\frac{q(x)}{\tilde{c}})}{Z}
\log\left(\min\left(\frac{\tilde{c}p(x)}{q(x)},1\right)\right)\mathrm{d}x-\log
(Z)-\log(\tilde{c})$
$\displaystyle=\int_{\mathbb{R}^{d}}\min\left(\frac{\min(p(x),\frac{q(x)}{
\tilde{c}})}{Z}\log\left(\frac{\tilde{c}p(x)}{q(x)}\right),0\right)\mathrm{d}x
-\log(Z)-\log(\tilde{c}).$
Since $\log\left(\frac{\tilde{c}p(x)}{q(x)}\right)\leq 0$ if and only if $p(x)\leq q(x)/\tilde{c}$ this can be reformulated as
$\displaystyle\quad\int_{\mathbb{R}^{d}}\min\left(\frac{p(x)}{Z}\log\left(\frac
{\tilde{c}p(x)}{q(x)}\right),0\right)\mathrm{d}x-\log(Z)-\log(\tilde{c})$
(39)
$\displaystyle\leq\int_{\mathbb{R}^{d}}\min\left(p(x)\log\left(\frac{\tilde{c}p
(x)}{q(x)}\right),0\right)\mathrm{d}x-\log(Z)-\log(\tilde{c}),$
(40)
where the inequality comes from the fact that $Z=\mathbb{E}[\alpha(X)]\in[0,1]$.
Moreover, it holds by Jensen’s inequality that
$\displaystyle-\log(Z)=-\log(\mathbb{E}[\alpha(X)])$
$\displaystyle=-\log\left(\int_{\mathbb{R}^{d}}p(x)\min\left(\frac{q(x)}{\tilde
{c}p(x)},1\right)\mathrm{d}x\right)$
(41)
$\displaystyle\leq-\int_{\mathbb{R}^{d}}p(x)\log\left(\min\left(\frac{q(x)}{
\tilde{c}p(x)},1\right)\right)\mathrm{d}x$
(42)
$\displaystyle=-\int_{\mathbb{R}^{d}}p(x)\min\left(\log\left(\frac{q(x)}{\tilde
{c}p(x)}\right),0\right)\mathrm{d}x$
(43)
$\displaystyle=\phantom{-}\int_{\mathbb{R}^{d}}p(x)\max\left(\log\left(\frac{
\tilde{c}p(x)}{q(x)}\right),0\right)\mathrm{d}x$
(44)
$\displaystyle=\phantom{-}\int_{\mathbb{R}^{d}}\max\left(p(x)\log\left(\frac{
\tilde{c}p(x)}{q(x)}\right),0\right)\mathrm{d}x.$
(45)
Thus, we obtain that (40) can be bounded from above by
$\int_{\mathbb{R}^{d}}\min\left(p(x)\log\left(\frac{\tilde{c}p(x)}{q(x)}\right)
,0\right)\mathrm{d}x+\int_{\mathbb{R}^{d}}\max\left(p(x)\log\left(\frac{\tilde
{c}p(x)}{q(x)}\right),0\right)\mathrm{d}x-\log(\tilde{c})$
which equals
$\int_{\mathbb{R}^{d}}p(x)\log\left(\frac{\tilde{c}p(x)}{q(x)}\right)\mathrm{d}
x-\log(\tilde{c})=\int_{\mathbb{R}^{d}}p(x)\log\left(\frac{p(x)}{q(x)}\right)
\mathrm{d}x=\mathrm{KL}(\mu,\nu).$
In summary, we have $\mathrm{KL}(\eta,\nu)\leq\mathrm{KL}(\mu,\nu)$, which implies the assertion.

. $\Box$

### C.2 Proof of Corollary 9

Let $Y=\Big{|}\mathbb{E}[\alpha(X)]-\frac{1}{N}\sum_{i=1}^{N}\alpha(X_{N})\Big{|}$ denote the random variable representing the error.
Since $\alpha(X_{N})\in[0,1]$, Hoeffding’s inequality Hoeffding () yields that
$P(Y>t)\leq 2\exp(-\frac{Nt^{2}}{2})$.
Consequently, we have that

$$ $\mathbb{E}[Y]=\int_{0}^{\infty}P(Y>t)\mathrm{d}t\leq 2\int_{0}^{\infty}\exp \left(-\frac{Nt^{2}}{2}\right)\mathrm{d}t=\frac{\sqrt{2\pi}}{\sqrt{N}}.$ $$

. $\Box$

## Appendix D Algorithms

### D.1 Training and Evaluation Algorithms

We summarize the training and evaluation procedures for the neural JKO steps in Algorithm 1 and 2 and for the importance-based rejection steps in Algorithm 3 and 4.

Figure: Algorithm 1 Training of neural JKO steps

### D.2 Density Evaluation of Importance Corrected Neural JKO Models

Let $(X_{0},...,X_{K})$ be an importance corrected neural JKO model. We aim to evaluate the density $p^{K}$ of $X_{K}$ at some given point $x\in\mathbb{R}^{d}$. Using Proposition 8, this can be done recursively by Algorithm 5. Note that the algorithm always terminates since $K$ is reduced by one in each call.

Figure: Algorithm 5 Density Evaluation of Importance Corrected Neural JKO Models

## Appendix E Additional Numerical Results and Implementation Details

### E.1 Test Distributions

We evaluate our method on the following test distributions.

- -
Mustache:
The two-dimensional log-density is given as $\log\mathcal{N}(0,\Sigma)\circ T$ with $\Sigma=[1,\sigma;\sigma,1]$ and $T(x_{1},x_{2})=(x_{1},(x_{2}-(x_{1}^{2}-1)^{2}))$. Note that $\mathrm{det}(\nabla T(x))=1$ for all $x$. In particular, we obtain directly by the transformation formula that the normalization constant is one.
Depending on $\sigma\in[0,1)$ close to $1$, this probability distribution has very long and narrow tails making it hard for classical MCMC methods to sample them. In our experiments we use $\sigma=0.9$.
- -
Shifted 8 Modes: A two-dimensional Gaussian mixture model with 8 equal weighted modes and covariance matrix $$1\text{\times}{10}^{-2}$I$. The modes are placed in a circle with radius 1 and center $(-1,0)$. Due to the shifted center classical MCMC methods have difficulties to distribute the mass correctly onto the modes.
- -
Shifted 8 Peaky: This is the same distribution as the shifted 8 Modes with the difference that we reduce the width of the modes to the covariance matrix $$5\text{\times}{10}^{-3}$I$. Since the modes are disconnected, it becomes harder to sample from them.
- -
Funnel: We consider the (normalized) probability density function given by $q(x)=\mathcal{N}(x_{1}|0,\sigma_{f}^{2})\mathcal{N}(x_{2:10}|0,\exp(x_{1})I)$, where $\sigma_{f}^{2}=9$. This example was introduced by Neal (). Similarly to the mustache example, this distribution has a narrow funnel for small values $x_{1}$ which can be hard to sample.
- -
GMM-$d$: A $d$-dimensional Gaussian mixture model with 10 equal weighted modes with covariance matrix $$1\text{\times}{10}^{-2}$I$ and means drawn randomly from a uniform distribution on $[-1,1]^{d}$. This leads to a peaky high-dimensional and multimodal probability distribution which is hard to sample from.
- -
LGCP: This is a high dimensional standard example taken from Arbel et al. (); Matthews et al. (); Vargas et al. (). It describes a Log-Gaussian Cox process on a $40\times 40$ grid as arising from spatial statistics Møller et al. (). This leads to a $1600$-dimensional probability distribution with the unnormalized density function
$q(x)\propto\mathcal{N}(x|\mu,K)\prod_{i\in\{1,...,40\}^{2}}\exp(x_{i}y_{i}-a
\exp(x_{i}))$, where $\mu$ and $K$ are a fixed mean and covariance kernel, $y_{i}$ is some data and $a$ is a hyperparameter. For details on this example and the specific parameter choice we refer to Matthews et al. ().

### E.2 Error Measures

We use the following error measures.

- -
The energy distance was proposed by Székely () and is defined by
$D(\mu,\nu)=-\frac{1}{2}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\|x-y\|
\mathrm{d}(\mu-\nu)(x)\mathrm{d}(\mu-\nu)(y).$
It is the maximum mean discrepancy with the negative distance kernel $K(x,y)=-\|x-y\|$ (see Sejdinovic et al., ) and can be estimated from below and above by the Wasserstein-1 distance (see Hertrich et al., ).
It is a metric on the space of probability measures. Consequently, a smaller energy distance indicates a higher similarity of the input distributions.
By discretizing the integrals it can be easily evaluated based on $N\in\mathbb{N}$ samples $\bm{x}=(x_{i})_{i}\sim\mu^{\otimes N}$, $\bm{y}=(y_{i})_{i}\sim\nu^{\otimes N}$ as
$D(\bm{x},\bm{y})=\sum_{i,j=1}^{N}\|x_{i}-y_{j}\|-\frac{1}{2}\sum_{i,j=1}^{N}\|
x_{i}-x_{j}\|-\frac{1}{2}\sum_{i,j=1}^{N}\|y_{i}-y_{j}\|.$
We use $N=50000$ samples in Table 1.
- -
We estimate the log normalizing constant (short $\log(Z)$ estimation) which is used as a benchmark standard in various references (e.g., in Arbel et al., ; Matthews et al., ; Phillips et al., ; Vargas et al., ). More precisely, for the generated distribution $\mu$ with normalized density $p$ and target measure $\nu$ with density $q(x)=g(x)/Z_{g}$ we evaluate the term
$\mathbb{E}_{x\sim\mu}\left[\log\left(\frac{g(x)}{p(x)}\right)\right]=\log(Z_{g
})-\mathbb{E}_{x\sim\mu}\left[\log\left(\frac{p(x)}{q(x)}\right)\right]=\log(Z
_{g})-\mathrm{KL}(\mu,\nu).$
Due to the properties of the KL divergence, a higher $\log(Z)$ estimate implies a lower KL divergence between $\mu$ and $\nu$ and therefore a higher similarity of generated and target distribution. In our experiments we compute the $\log(Z)$ estimate based on $N=50000$ samples. The results are given in Table 2.
- -
To quantify how well the mass is distributed on different modes for the mixture model examples (shifted 8 Modes, shifted 8 Peaky, GMM-$d$), we compute the mode weights. That is, we generate $N=50000$ samples and assign each generated samples to the closest mode of the GMM. Afterwards, we compute for each mode the fraction of samples which is assigned to each mode. To evaluate this distribution quantitatively, we compute the mean square error (MSE) between the mode weights of the generated samples and the ground truth weights from the GMM. We call this error metric the mode MSE, give the results are in Table 3.

###### Remark 16 (Bias in $\log(Z)$ Computation) .

In the cases, where the importance corrected neural JKO sampling fits the target distribution almost perfectly, we sometimes report in Table 2 $\log(Z)$ estimates which are slightly larger than the ground truth. This can be explained by the fact that the density evaluation of the continuous normalizing flows uses the Hutchinson trace estimator for evaluating the divergence and a numerical method for solving the ODE. Therefore, we have a small error in the density propagation of the neural JKO steps. This error is amplified by the rejection steps since samples with underestimated density are more likely to be rejected than samples with overestimated density.

We would like to emphasize that this effect only appears for examples, where the energy distance between generated and ground truth samples is in the same order of magnitude like the average distance between two different sets of ground truth samples (see Table 1). This means that in the terms of the energy distance the generated and ground truth distribution are indistinguishable. At the same time the bias in the $\log(Z)$ estimate appears at the third or fourth relevant digit meaning that it is likely to be negligible.

**Table 3: We report the mode MSEs for the different methods for all examples which can be represented as mixture model. A smaller mode MSE indicates a better result.**
|  |  | Sampler |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Distribution |  | MALA | HMC | DDS | CRAFT | Neural JKO | Neural JKO IC (ours) |
| shifted 8 Modes |  | $$3.2\text{\times}{10}^{-3}$\pm$1.3\text{\times}{10}^{-4}$$ | $$1.9\text{\times}{10}^{-5}$\pm$1.1\text{\times}{10}^{-5}$$ | $$8.8\text{\times}{10}^{-3}$\pm$2.2\text{\times}{10}^{-3}$$ | $$3.2\text{\times}{10}^{-2}$\pm$6.1\text{\times}{10}^{-3}$$ | $$8.3\text{\times}{10}^{-2}$\pm$1.0\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$1.3\text{\times}{10}^{-5}$\pm$4.4\text{\times}{10}^<br>{-6}$$ |
| shifted 8 Peaky |  | $$8.3\text{\times}{10}^{-2}$\pm$3.7\text{\times}{10}^{-4}$$ | $$7.8\text{\times}{10}^{-2}$\pm$9.9\text{\times}{10}^{-4}$$ | $$7.9\text{\times}{10}^{-3}$\pm$2.2\text{\times}{10}^{-3}$$ | $$3.3\text{\times}{10}^{-2}$\pm$1.3\text{\times}{10}^{-2}$$ | $$8.4\text{\times}{10}^{-2}$\pm$6.7\text{\times}{10}^{-4}$$ | $\pagecolor{green!15!white}$1.5\text{\times}{10}^{-5}$\pm$3.2\text{\times}{10}^<br>{-6}$$ |
| GMM-10 |  | $$1.2\text{\times}{10}^{-2}$\pm$5.5\text{\times}{10}^{-3}$$ | $$1.0\text{\times}{10}^{-2}$\pm$5.8\text{\times}{10}^{-3}$$ | $$3.6\text{\times}{10}^{-3}$\pm$1.9\text{\times}{10}^{-3}$$ | $$1.4\text{\times}{10}^{-1}$\pm$4.7\text{\times}{10}^{-2}$$ | $$1.1\text{\times}{10}^{-2}$\pm$5.6\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$2.1\text{\times}{10}^{-5}$\pm$3.0\text{\times}{10}^<br>{-6}$$ |
| GMM-20 |  | $$6.6\text{\times}{10}^{-3}$\pm$2.2\text{\times}{10}^{-3}$$ | $$6.6\text{\times}{10}^{-3}$\pm$2.4\text{\times}{10}^{-3}$$ | $$3.6\text{\times}{10}^{-3}$\pm$1.5\text{\times}{10}^{-3}$$ | $$3.8\text{\times}{10}^{-1}$\pm$3.9\text{\times}{10}^{-2}$$ | $$7.3\text{\times}{10}^{-3}$\pm$2.8\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$2.4\text{\times}{10}^{-5}$\pm$9.8\text{\times}{10}^<br>{-6}$$ |
| GMM-50 |  | $$1.0\text{\times}{10}^{-2}$\pm$3.8\text{\times}{10}^{-3}$$ | $$9.8\text{\times}{10}^{-3}$\pm$3.6\text{\times}{10}^{-3}$$ | $$9.6\text{\times}{10}^{-3}$\pm$4.9\text{\times}{10}^{-3}$$ | $$9.0\text{\times}{10}^{-1}$\pm$6.4\text{\times}{10}^{-5}$$ | $$1.2\text{\times}{10}^{-2}$\pm$4.4\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$2.6\text{\times}{10}^{-5}$\pm$6.0\text{\times}{10}^<br>{-6}$$ |
| GMM-100 |  | $$1.1\text{\times}{10}^{-2}$\pm$5.4\text{\times}{10}^{-3}$$ | $$1.1\text{\times}{10}^{-2}$\pm$5.4\text{\times}{10}^{-3}$$ | $$1.1\text{\times}{10}^{-2}$\pm$5.8\text{\times}{10}^{-3}$$ | $$9.0\text{\times}{10}^{-1}$\pm$4.8\text{\times}{10}^{-8}$$ | $$1.4\text{\times}{10}^{-2}$\pm$7.0\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$1.5\text{\times}{10}^{-4}$\pm$8.5\text{\times}{10}^<br>{-5}$$ |
| GMM-200 |  | $$1.4\text{\times}{10}^{-2}$\pm$5.1\text{\times}{10}^{-3}$$ | $$1.4\text{\times}{10}^{-2}$\pm$4.6\text{\times}{10}^{-3}$$ | $$2.2\text{\times}{10}^{-2}$\pm$8.6\text{\times}{10}^{-3}$$ | $$9.0\text{\times}{10}^{-1}$\pm$5.8\text{\times}{10}^{-8}$$ | $$1.9\text{\times}{10}^{-2}$\pm$6.3\text{\times}{10}^{-3}$$ | $\pagecolor{green!15!white}$6.8\text{\times}{10}^{-4}$\pm$4.0\text{\times}{10}^<br>{-4}$$ |

###### Remark 16 (Bias in log ⁡ ( Z ) 𝑍 \log(Z) roman_log ( italic_Z ) Computation) .

### E.3 Additional Figures and Evaluations

Additionally to the results from the main part of the paper, we provide the following evaluations.

#### Visualization of the Rejection Steps

In Figure 1 and with involved steps in Figure 2, we visualize the different steps of our importance corrected neural JKO model on the shifted 8 Peaky example. Due to the shift of the modes, the modes on the right attract initially more mass than the ones on the left. In the rejection layers, we can see that samples are mainly rejected in oversampled regions such that the mode weights are equalized over time. This can also be seen in Figure 3, where we plot the weights of the different modes over time. We observe, that these weights are quite imbalanced for the latent distribution but are equalized over the rejection steps until they reach the ground truth value.

#### (Marginal) Plots of Generated Samples

Despite the quantitative comparison of the methods in the Tables 1, 2 and 3, we also plot the first to coordinates of generated samples for the different test distributions for all methods for a visual comparison.

In Figure 4, we plot samples of the shifted 8 Modes example. We can see that all methods roughly cover the ground truth distribution even though CRAFT and DDS have slight and the uncorrected neural JKO scheme has a severe imbalance in the assigned mass for the different modes.

For the shifted 8 Peaky example in Figure 5, we see that this imbalance increases drastically for CRAFT, HMC, MALA and the uncorrected neural JKO. Also DDS has has a slight imbalance, while the importance corrected neural JKO scheme fits the ground truth almost perfectly.

The Funnel distribution in Figure 6 has two difficult parts, namely the narrow but high-density funnel on the left and the wide moderate density fan on the right. We can see that DDS does cover none of them very well. Also MALA, CRAFT and the uncorrected neural JKO do not cover the end of the funnel correctly and have also difficulties to cover the fan. HMC covers the fan well, but not the funnel. Only our importance corrected neural JKO scheme covers both parts in a reasonable way.

For the Mustache distribution in Figure 7, the critical parts are the two heavy but narrow tails. We observe that CRAFT and DDS are not able to cover them at all, while HMC and MALA and uncorrected neural JKO only cover them only partially. The importance corrected neural JKO covers them well.

Finally, for the GMM-$d$ example we consider the dimensions $d=10$ and $d=200$ in the Figures 8 and 9. We can see that CRAFT mode collapses, i.e., for $d=10$ it already finds some of the modes and for $d=200$ it only finds one mode. DDS, HMC, MALA and the uncorrected neural JKO find all modes but do not distribute the mass correctly onto all modes. While this already appears for $d=10$ it is more severe for $d=200$. The importance corrected neural JKO sampler finds all modes and distributes the mass accurately.

#### Development of Error Measures over the Steps

We plot how the quantities of interest decrease over the application of the steps of our model. The results are given in the Figure 10.
It can be observed that the different steps may optimize different metrics. While the rejection steps improve the $\log(Z)$ estimate more significantly, the JKO steps focus more on the minimization of the energy distance.
Overall, we see that in all figures the errors decrease over time.

### E.4 Implementation Details

To build our importance corrected neural JKO model, we first apply $n_{1}\in\mathbb{N}$ JKO steps followed by $n_{2}\in\mathbb{N}$ blocks consisting out of a JKO step and three importance-based rejection steps. The velocity fields of the normalizing flows are parameterized by a dense three-layer feed-forward neural network. For the JKO steps, we choose an initial step size $\tau_{0}>0$ and then increase the step size exponentially by $\tau_{k+1}=4\tau_{k}$. The choices of $n_{1}$, $n_{2}$, $\tau_{0}$ and the number of hidden neurons from the networks is given in Table 4 together with the execution times for training and sampling. For evaluating the density propagation through the CNFs, we use the Hutchinson trace estimator with $5$ Rademacher distributed random vectors whenever $d>5$ and the exact trace evaluation otherwise. For implementing the CNFs, we use the code from Ffjord (Grathwohl et al., ) and the torchdiffeq library by Chen (). We provide the code online at

[https://github.com/johertrich/neural_JKO_ic](https://github.com/johertrich/neural_JKO_ic).

For the MALA and HMC we run an independent chain for each generated samples and perform $50000$ steps of the algorithm. For HMC we use $5$ momentum steps and set the step size to $0.1$ for 8Modes and Funnel and step size $0.01$ for mustache. To stabilize the first iterations of MALA and HMC we run the first iterations with smaller step sizes ($0.01$ times the final step size in the first $1000$ iterations and $0.1$ times the final step size for the second $1000$ iterations).
For MALA we use step size $0.001$ for all examples. For DDS and CRAFT we use the official implementations. In particular for the test distributions such as Funnel and LGCP we use the hyperparameters suggested by the original authors. For the other examples we optimized them via grid search.

**Table 4: Parameters, training and sampling times for the different examples. For the sampling time we draw $N=50000$ samples once the method is trained. The execution times are measured on a single NVIDIA RTX $4090$ GPU with $24$ GB memory.**
|  |  | Distribution |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Parameter |  | Mustache | shifted 8 Modes | shifted 8 Peaky | Funnel | GMM-10 | GMM-20 | GMM-50 | GMM-100 | GMM-200 | LGCP |
| Dimension |  | 2 | 2 | 2 | 10 | 10 | 20 | 50 | 100 | 200 | 1600 |
| Number $n_{1}$ of flow steps |  | 6 | 2 | 2 | 6 | 4 | 4 | 4 | 4 | 5 | 3 |
| Number $n_{2}$ of rejection blocks |  | 6 | 4 | 4 | 6 | 6 | 6 | 7 | 8 | 8 | 6 |
| Initial step size $\tau_{0}$ |  | 0.05 | 0.01 | 0.01 | 5 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.001 | 5 |
| Hidden neurons |  | 54 | 54 | 54 | 256 | 70 | 90 | 150 | 250 | 512 | 1024 |
| Training time (min) |  | 38 | 20 | 21 | 107 | 33 | 34 | 44 | 53 | 80 | 163 |
| Sampling time (sec) |  | 15 | 2 | 3 | 79 | 22 | 23 | 72 | 129 | 473 | 535 |

Figure: Figure 2: Visualization of the steps of the importance corrected neural JKO model for the shifted 8 Peaky example. We start at the top left with the latent distribution and apply in each image one step from the model. The red color indicates samples which are rejected in the next step and the green color marks the resampled points.
Refer to caption: extracted/5762551/img/plotdata_8peaky_2.png

Figure: Figure 4: Sample generation with various methods for shifted $8$ Modes example with ground truth samples and generated samples for each associated method. While most methods recover the distribution well, we can see a imbalance in the modes for the uncorrected neural JKO, CRAFT and DDS.
Refer to caption: extracted/5762551/img/plot_methodSamples_8mixtures_2.png

Figure: Figure 6: Marginalized sample generation with various methods for the $d=10$ funnel distribution with ground truth samples and generated samples for each associated method. We observe that only the importance corrected neural JKO covers the thin part of the funnel well.
Refer to caption: extracted/5762551/img/plot_methodSamples_funnel_10.png

Figure: Figure 7: Sample generation with various methods for the $d=2$ mustache distribution with ground truth samples and generated samples for each associated method. We can see that MALA, CRAFT and DDS have difficulties to model the long tails of the distribution properly.
Refer to caption: extracted/5762551/img/plot_methodSamples_schnauzbart_2.png

Figure: Figure 8: Marginalized sample generation with various methods for the GMM-$10$ distribution with ground truth samples and generated samples for each associated method. We observe that CRAFT mode collapses and that only the importance corrected neural JKO model distributes the mass correctly onto the modes.
Refer to caption: extracted/5762551/img/plot_methodSamples_mixtures_10.png

Figure: Figure 10: We plot the energy distance (left) and $\log(Z)$ estimate (right) over the steps of our importance corrected neural JKO method for different examples. We observe that the error measures decrease in the beginning and then saturate at some value.
Refer to caption: x2.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: 1. Please give details about the importance sampling scheme: the introduction in Sec. 4 could be shortened so that the **Neural JKO Sampling with Importance Correction** is expanded and the scene proposed clarified. 

2. The paper blends CNF to perform steepest descent in continuous-time via the solution of ODE, with the JKO scheme which is a discrete-time proximal scheme. Doing so therefore requires discretizing the ODE, which of course is a step that always needs to be done, but the details of this discretization step should be discussed more carefully. In particular, how is the parameter $\tau$ in the JKO scheme chosen: is the same at every step?  Is the discretization step used to solve the CNF (and related ODE) taken to be $\tau$, or are multiple integration steps performed for every $\tau$? Can the time step by adapted along the way?

3. The method seems to require restarting the integration of the CNF from the beginning each time a rejection step is performed. Is this indeed the case? This seems quite costly in practice. 

4. Several papers have recently appeared that propose similar schemes, with an importance sampling step also added, see in particular:

https://arxiv.org/abs/2410.02711

https://arxiv.org/abs/2102.07501

https://arxiv.org/abs/2111.15141

https://arxiv.org/abs/2307.01050

These papers should be discussed here, by giving a careful comparison between the methods and their numerical results.

5. More details about the numerics should be given.
soundness: 3
strengths: The problem studied is an important one, and the algorithm proposed may be useful to address it.
summary: The paper considers the problem of sampling a probability distribution known up to a normalization constant, a classical question of interest in statistical mechanics and bayesian inference. It is proposed to do so by combining JKO steps using continuous normalizing flows (CNF) based on neural ODE with importance based rejection steps to correct for the bias introduced by errors in the CNF. More specifically, the JKO scheme is used to build a sequence of piecewise geodesic interpolations between densities performing steepest descent in Wasserstein-2 metric using as objective function the reversed KL divergence of the evolving measure from the target as energy. Since this KL divergence can be expressed as an expectation over the evolving measure of the unnormalized target distribution, it can be written explicitly. The optimization can thereby be performed in steps, which is practice can be done by propagating samples from one measure along the sequence to the next while solving some ODE along the way. A key aspect of the method is claimed to be the addition of an importance sampling step to accept or reject these samples along te way, in order to remove statistical errors that could accumulate as the walkers used in the sampling slowly drift from the sequence of measure generated.
weaknesses: The idea of building JKO schemes using is not new, see in particular the paper by Xu, Cheng, and Tie (cited here),. It would therefore be helpful to better emphasizes what is new here (which seems to be the importance sampling step, and demonstrate via example that it solves the slow convergence and instability issues that otherwise could arise. In this respect, the author could summarize (or relegate to the appendix) the background material on Wasserstein spaces, CNF and OT,  the JKO flow, and AIS as well as Metropolis-Hastings: while this material is explained well here, this exposition takes a lot of space and as a result the details on the actual scheme being proposed are rather sparse (in particular in the way the importance sampling bit of the scheme is implemented), and so are the numerical experiments reported.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: To the reviewer, the rejection sampling step seems more like sequential monte-carlo than it does annealed importance sampling. These methods share a lot of the same conceptual reasoning, and the reviewer thinks in fact that the resampling procedure in SMC could be interpreted in this context too. Could the authors comment on that? Do they think it would be useful to make this connection?
soundness: 3
strengths: - The authors exposit their method very nicely, and in addition make useful connections to various facets of the transport problem in the literature that are necessary for understanding the context of their method -- e.g. the motivation for the JKO scheme and its relation to regularized-OT solutions. 

- The experiments are comprehensive as with those in the literature.

- Having a correctable JKO scheme is important. Given the reviewer's current knowledge, this was absent from previously proposed solutions, and for sampling algorithms, being able to show that your method is unbiased (or has controlled bias) is essential.
summary: This paper combines ideas from JKO sampling and annealed importance sampling to established an unbiased algorithm for implementing the Neural JKO protocol. They establish that using a learned velocity field to perform the transport monotonically decreases the KL divergence between the model density at time $t$ to the target density (as we should expect from a transport with an implied interpolating density). As the author understands it, the main contribution of this paper is the establishment of importance weights that allow to correct errors arising from the approximate velocity field in the neural JKO scheme.
weaknesses: - As the reviewer understands it, the proposed method requires backpropagation through the solution of the ODE for training. Could the authors verify that this is true? This has been a perennial issue in using KL-like optimization for training samplers (it does not scale well with the dimensionality of the problem), and it's not clear to the reviewer that this primary obstacle is avoided in the OT regularized version.  

- Moreover, this scales even worse when the ODE solution requires the computation of the trace of the divergence of a vector field. A measure of computational cost compared to the relevant literature / experiments seems important. In high dimensional sampling problems, I imagine the memory footprint of training probably has to be quite large. 

- These issues are left out of the limitation discussion, which seems important, in addition to what they already remark on that the number of discretization steps (and therefore opportunity for rejection sampling) may be large if the flow is not learned well. Adding these reflections would be useful (and does not retract from the contribution of the authors!), and would allow the reviewer to bump their score a bit!

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 4
questions: - How does the neural JKO approach differ from the rectified flow / flow-matching objective?
- Is it possible to avoid the adjoint formulation, which is widely known to be inefficient?
- Some of the initial motivation focuses on mode collapse. Do any of your examples clearly indicate that this method avoids mode collapse where standard CNFs fail? A simple example in 2d with well-separated modes is typically enough to observe mode collapse when optimizing the reverse KL.
soundness: 3
strengths: - The paper is technically strong throughout and provides a clear assessment of the algorithm on a suite of relatively standard examples for the MCMC literature.
	- The clear connections to optimal transport provide a nice framework to enable connections to other algorithms.
summary: This paper describes a sampling algorithm using proposals from a continuous normalizing flow with an additional reweighting / accept-reject step, which the authors claim enables the use of a data-free variational inference objective that reduces the reverse KL divergence. The methodology is analyzed and studied numerically on simple examples, where it outperforms a number of standard methods and more recent ML-based strategies like CRAFT and DDS.
weaknesses: - It is somewhat difficult to identify the contribution / novelty of this work. The basic idea of using reweighting and or metropolization has been used heavily with normalizing flows. E.g., https://arxiv.org/abs/2105.12603 or https://arxiv.org/abs/1812.01729. Much of the paper is standard background material. 
- The connection between the objectives written in section 3.2 and the extremely widely used flow-matching / rectified flow objectives is not explored at all. While both of the aforementioned algorithms use (or can be directly formulated with) a Benamou-Brenier OT objective, the authors choose to instead use the more tedious and less efficient adjoint state framework for CNFs. 
- The importance sampling algorithm is not very clearly explained. Because it seems to be the primary novel content in the paper, it would be beneficial to focus more heavily on this. 
- Not at all clear from the paper how this scheme avoids mode collapse.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: Besides the questions in the "Weaknesses" section, I have a few other minor ones that I hope the author(s) can clarify.

1. In Theorem 3, the curve is denoted as $\gamma_\tau$, but I only see $\tilde{\gamma}_\tau$ in equation (3). Are they the same?

2. In Section 4, the "Annealed Importance Sampling" paragraph, I would like to point out that the self-normalized estimator $\sum_{i=1}^N w_i f(x_i)/(\sum_{i=1}^N w_i)$ is in general not unbiased, since the denominator is also a random variable. The numerator itself $\sum_{i=1}^N w_i f(x_i)$ is actually unbiased. The self-normalized version is used when $w_i$ is computed from unnormalized densities.

3. In the introduction section, I see "this regularization converts the objective functional into a convex one". Is there any theoretical guarantee for this?
soundness: 3
strengths: Sampling is a fundamental task in machine learning, and this article tackles this problem by utilizing some recent advances in deep generative modeling. The use of JKO scheme to iteratively optimizing the model parameters is also interesting.
summary: In this article, the author(s) consider the problem of sampling from unnormalized density functions, and the basic idea is to implicitly construct a diffeomorphism to transport a simple distribution to the target distribution, based on continuous normalizing flows (CNFs). To further improve the quality of sampling, a rejection-resampling step is also used.
weaknesses: 1. Overall I like the idea of this article, but one of my majors concerns on this article is the motivation of using a continuous normalizing flow (CNF) rather than the classical discrete-time normalizing flow. If my understanding is correct, CNF uses neural networks to model the velocity field $v_\theta(\cdot)$, and then the transport map $\mathcal{T}\_\theta(\cdot)$ is implicitly defined by $v_\theta(\cdot)$. However, when minimizing the (reverse) KL loss function, it would be easier to work on $\mathcal{T}\_\theta(\cdot)$, which is what classical normalizing flows do. Also, the sampling of classical normalizing flows should be faster than CNFs, since they do not need to solve an ODE.

2. Moreover, for classical normalizing flows, the evaluation of the (log-)density function is *exact*, which is important for rejection sampling. CNF needs to approximate the Jacobian term in evaluating the density.

3. There are some notable existing works on sampling using diffeomorphisms/normalizing flows, for example [1][2][3]. When comparing deep-learning-based sampling algorithms for the experiments, I would suggest at least adding one of them, since they are closer to the method used in this article.

4. In the last paragraph of the introduction section, I think it is helpful to distinguish *modeling* methods (learning from data) with *sampling* methods (learning from a density function). The latter is closer to the topic of this article.

5. For the evaluation metric in the experiments, isn't Wasserstein distance a more natural choice? It has the benefit of avoiding selecting the kernel as a hyperparameter.

[1] Marzouk, Y., Moselhy, T., Parno, M., & Spantini, A. (2016). An introduction to sampling via measure transport. arXiv preprint arXiv:1602.05023.

[2] Hoffman, M., Sountsov, P., Dillon, J. V., Langmore, I., Tran, D., & Vasudevan, S. (2019). Neutra-lizing bad geometry in hamiltonian monte carlo using neural transport. arXiv preprint arXiv:1903.03704.

[3] Qiu, Y., & Wang, X. (2024). Efficient multimodal sampling via tempered distribution flow. Journal of the American Statistical Association.

Meta-review: This paper proposes a novel method for sampling from unnormalized distributions by combining continuous normalizing flows (CNFs) with rejection sampling based on importance weights to correct biases. The CNF framework enables fast sampling and density evaluation using neural networks, while the biases introduced by CNFs are appropriately corrected through rejection sampling.

All reviewers acknowledged the significant contribution of the JKO scheme in this study. However, concerns were raised regarding the lack of discussion on importance sampling, computational costs, and parameter tuning. As the paper is borderline, I conducted a detailed review and reached a similar conclusion. While the mathematical contributions are substantial, as an algorithmic proposal, the paper requires major revisions, as also highlighted by Reviewer wGJ8.

Specifically, although the authors added a pseudo code (Algorithm 1) in response to reviewer comments, the algorithm calls several additional sub-algorithms (Algorithms 2–4 in the Appendix), which significantly impacts readability. Moreover, the critical aspect of how the JKO scheme is combined with importance sampling is not explicitly detailed in the pseudo code. For instance, the description, "In practice, we first utilize CNF layers only, then use several blocks consisting of a single CNF layer composed with 3 rejection steps", remains vague. While Appendix E.2 provides details about the combination, this should be clearly presented in the main text as an explicit algorithm.

Additionally, there is no discussion on the stability of the proposed method with respect to how the JKO scheme and importance sampling are combined. The pseudo code also lacks details about key hyperparameters, such as step sizes, which are critical for the algorithm’s implementation. Multiple reviewers raised concerns about hyperparameters, and while these were partially addressed during the discussion phase, such discussions should also be included in the main text due to their importance.

Another significant aspect of the study is the ability to calculate densities, which is an important advantage of the proposed method. However, the specific method for this calculation is only described in the pseudo code in the Appendix. This too should be presented in the main text for clarity and completeness.

Addressing these issues requires substantial revisions to ensure the algorithmic aspects are adequately described. While the mathematical contributions of this study are significant, there is a notable discrepancy in quality between the theoretical and algorithmic components. Given the inadequacies in the latter, I recommend rejection at this time, with the hope that the paper will be resubmitted after major revisions.

Decision: Reject