{'fraction_non_alphanumeric': 0.0449696799887926, 'fraction_numerical': 0.01949286529109212, 'mean_word_length': 4.137569401603948, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Building Normalizing Flows with Stochastic Interpolants

A simple generative model based on a continuous-time normalizing flow between any pair of base and target probability densities is proposed. The velocity field of this flow is inferred from the probability current of a time-dependent density that interpolates between the base and the target in finite time. Unlike conventional normalizing flow inference methods based the maximum likelihood principle, which require costly backpropagation through ODE solvers, our interpolant approach leads to a simple quadratic loss for the velocity itself which is expressed in terms of expectations that are readily amenable to empirical estimation. The flow can be used to generate samples from either the base or target, and to estimate the likelihood at any time along the interpolant. In addition, the flow can be optimized to minimize the path length of the interpolant density, thereby paving the way for building optimal transport maps. The approach is also contextualized in its relation to diffusions. In particular, in situations where the base is a Gaussian density, we show that the velocity of our normalizing flow can also be used to construct a diffusion model to sample the target as well as estimating its score. This allows one to map methods based on stochastic differential equations to those using ordinary differential equations, simplifying the mechanics of the model, but capturing equivalent dynamics. Benchmarking on density estimation tasks illustrates that the learned flow can match and surpass maximum likelihood continuous flows at a fraction of the conventional ODE training costs.


INTRODUCTION
Contemporary generative models have primarily been designed around the construction of a map between two probability distributions that transform samples from the first into samples from the second. While progress has been from various angles with tools such as implicit maps (Goodfellow et al., 2014;Brock et al., 2019), and autoregressive maps (Menick & Kalchbrenner, 2019;Razavi et al., 2019;Lee et al., 2022), we focus on the case where the map has a clear associated probability flow. Advances in this domain, namely from flow and diffusion models, have been governed by the introduction of algorithms or inductive biases that make learning this map, and the Jacobian of the associated change of variables, more tractable. The challenge is to choose what structure to impose on the transport to best reach a complex target distribution from a simple one used as base, while maintaining computational efficiency.
In the continuous time perspective, this problem can be framed as the design of a time-dependent map, X t (x) with t ∈ [0, 1], which functions as the push-forward of the base distribution at time t = 0 onto some time-dependent distribution that reaches the target at time t = 1. Assuming that these distributions have densities supported on Ω ⊆ R d , say ρ 0 for the base and ρ 1 for the target, this amounts to constructing X t : Ω → Ω such that if x ∼ ρ 0 then X t (x) ∼ ρ t for some density ρ t such that ρ t=0 = ρ 0 and ρ t=1 = ρ 1 .
(1) One convenient way to represent this time-continuous map is to define it as the flow associated with the ordinary differential equation (ODE)  Figure 1: The density ρ t (x) produced by the stochastic interpolant based on (5) between a standard Gaussian density and a Gaussian mixture density with three modes. Also shows in white are the flow lines of the map X t (x) our method produces.
where the dot denotes derivative with respect to t and v t (x) is the velocity field governing the transport. This is equivalent to saying that the probability density function ρ t (x) defined as the pushforward of the base ρ 0 (x) by the map X t satisfies the continuity equation (see e.g. (Villani, 2009;Santambrogio, 2015) and Appendix A) ∂ t ρ t + ∇ · (v t ρ t ) = 0 with ρ t=0 = ρ 0 and ρ t=1 = ρ 1 , and the inference problem becomes to estimate a velocity field such that 3 holds.
In this paper we propose a simple solution to this problem based on introducing a time-differentiable interpolant I t : Ω × Ω → Ω such that I t=0 (x 0 , x 1 ) = x 0 and I t=1 (x 0 , x 1 ) = x 1 (4) A useful instance of such an interpolant that we will employ is I t (x 0 , x 1 ) = cos( 1 2 πt)x 0 + sin( 1 2 πt)x 1 , (5) though we stress the framework we propose applies to any I t (x 0 , x 1 ) satisfying (4) under mild additional assumptions on ρ 0 , ρ 1 , and I t specified below. Given this interpolant, we then construct the stochastic process x t by sampling independently x 0 from ρ 0 and x 1 from ρ 1 , and passing them through I t : x t = I t (x 0 , x 1 ), x 0 ∼ ρ 0 , x 1 ∼ ρ 1 independent.
(6) We refer to the process x t as a stochastic interpolant. Under this paradigm, we make the following key observations as our main contributions in this work: • The probability density ρ t (x) of x t , henceforth referred to as the interpolant density, satisfies (3) with a velocity v t (x) which is the unique minimizer of a simple quadratic objective. This result is the content of Proposition 1 below, and it can be leveraged to estimate v t (x) in a parametric class (e.g. using deep neural networks) to construct a generative model through the solution of the probability flow equation (2). • By specifying an interpolant density, the method therefore separates the tasks of minimizing the objective from discovering a path between the base and target densities. This is in contrast with conventional maximum likelihood (MLE) training of flows where one is forced to couple the choice of path in the space of measures to maximizing the objective. • We show that the Wasserstein-2 (W 2 ) distance between the target density ρ 1 and the densityρ 1 obtained by transporting ρ 0 using an approximate velocityv t in (2) is controlled by our objective function. We also show that the value of the objective onv t during training can be used to check convergence of this learned velocity field towards the exact v t . • We show that our approach can be generalized to shorten the path length of the interpolant density and optimize the transport by additionally maximizing our objective over the interpolant I t (x 0 , x 1 ) and/or adjustable parameters in the base density ρ 0 . • By choosing ρ 0 to be a Gaussian density and using (5) as interpolant, we show that the score of the interpolant density, ∇ log ρ t , can be explicitly related to the velocity field v t . This allows us to draw connection between our approach and score-based diffusion models, providing theoretical groundwork for future exploration of this duality. • We demonstrate the feasibility of the method on toy and high dimensional tabular datasets, and show that the method matches or supersedes conventional ODE flows at lower cost, as it avoids the need to backpropagate through ODE solves.

RELATED WORKS
Early works on exploiting transport maps for generative modeling go back at least to Chen & Gopinath (2000), which focuses on normalizing a dataset to infer its likelihood. This idea was brought closer to contemporary use cases through the work of Tabak & Vanden-Eijnden (2010) and Tabak & Turner (2013), which devised to expressly map between densities using simple transport functions inferred through maximum likelihood estimation (MLE). These transformations were learned in sequence via a greedy procedure. We detail below how this paradigm has evolved in the case where the map is represented by a neural network and optimized accordingly.
Discrete and continuous time flows. Instead of working with the continuous map X t (x) as described in Section 1, the first success of normalizing flows with neural network parametrizations follow the work of Tabak & Turner (2013) with a finite set of steps along the map. By imposing structure on the transformation so that it remains an efficiently invertible diffeomorphism, the models of Rezende & Mohamed (2015); Dinh et al. (2017); Huang et al. (2018); Durkan et al. (2019) can be optimized through maximum likelihood estimation. However, this comes at the cost of limiting the expressive power of the representation, as the Jacobian of the map must be kept simple to calculate the likelihood. Extending this to the continuous case allowed the Jacobian to be unstructured yet still estimable through trace estimation techniques (Chen et al., 2018;Grathwohl et al., 2019;Hutchinson, 1989). However learning this map through MLE requires costly backpropagation through numerical integration. The research of Finlay et al. (2020); Onken et al. (2021) shows that the number of solver calls can be reduced with the addition of some path optimality regularization, though this does not alleviate the main structural challenge of the optimization. Our work uses a continuous map X t as well but allows for direct estimation of the underlying velocity.
Score-based flows. Adjacent research has made use of diffusion processes, commonly the Ornstein-Uhlenbeck (OU) process, to connect the target ρ 1 to the base ρ 0 . In this case the transport is governed by a stochastic differential equation (SDE) that is evolved for infinite time, and the challenge of learning a generative model can be framed as fitting the reverse time evolution of the SDE from Gaussian noise back to ρ 1 (Sohl-Dickstein et al., 2015;Ho et al., 2020;Song et al., 2021b). Doing so indirectly learns the velocity field by means of learning the score function ∇ log ρ t (x), instead of through MLE. While this approach has shown great promise to model high dimensional distributions (Rombach et al., 2022;Hoogeboom et al., 2022), particularly in the case of text-to-image generation (Ramesh et al., 2022;Saharia et al., 2022), there is an absence of theoretical motivation for the SDEflow framework and the complexity it induces. Namely, the SDE must evolve for infinite time to connect the distributions, the parameterization of the time steps remains heuristic Feller (1949); Xiao et al. (2022), and the criticality of noise, as well as the score, is not absolutely apparent Bansal et al. (2022); Lu et al. (2022). This motivates further research into effective, ODE-driven, approaches to learning the map. Indeed, under certain conditions, the primary SDE in score-based diffusion can already be associated to its ODE probability flow through a lower bound (Song et al., 2021a). Our approach can be viewed as an alternative to score-based diffusion models in which the ODE velocity is learned through the stochastic interpolant x t rather than the solution of an OU process.
Interpolants and optimal transport. Data interpolants were recently used in (Choi et al., 2022) for the purposes of density ratio estimation. Here we make use of them to directly discern a velocity field behind a generative model obeying the continuity equation. Co-incident work by  derives an analogous optimization to us, with a focus on straight interpolants, also contrasting it with score-based methods. Interestingly, they describe an iterative way of rectifying the interpolant path, which can be shown to arrive at an optimal transport map when the procedure is repeated ad infinitum (Liu, 2022). Here we propose an alternative and more direct solution to the problem of optimal transport that involves optimizing our objective over the stochastic interpolant. This is discussed in more detail in Section 2.2.

NOTATIONS AND ASSUMPTIONS
We assume that the base and the target distribution are both absolutely continuous with respect to the Lebesgue measure on R d , with densities ρ 0 and ρ 1 , respectively. We do not require these densities to be positive everywhere on R d , but we assume that ρ 0 (x) and ρ 1 (x) are continuously differentiable in x.
Regarding the interpolant I t : R d × R d → R d , we assume that it is surjective for all t ∈ [0, 1] and satisfies (4). We also assume that I t (x 0 , x 1 ) is continuously differentiable in (t, x 0 , x 1 ), and that it is such that E |∂ t I t (x 0 , x 1 )| 2 < ∞ (7) A few additional technical assumptions on ρ 0 , ρ 1 and I t are listed in Appendix B.
Given any function f t (x 0 , x 1 ) we denote its expectation over t, x 0 , and x 1 drawn independently from the uniform density on [0, 1], ρ 0 , and ρ 1 , respectively. We use ∇ to denote the gradient operator and the dot to denote the derivative with respect to time t.

STOCHASTIC INTERPOLANTS AND ASSOCIATED FLOWS
Our first main theoretical result can be phrased as follows: Proposition 1. The stochastic interpolant x t defined in (6) with I t (x 0 , x 1 ) satisfying (4) has a probability density ρ t (x) that satisfies the continuity equation (3) with a velocity v t (x) which is the unique minimizer overv t (x) of the objective In addition the minimum value of this objective is given by Both these statements also hold if G(v) is minimized over velocities that are gradient fields, in which case the minimizer is of the form v t (x) = ∇φ t (x) for some potential φ t : R d → R uniquely defined up to a constant.
Proposition 1 is proven in Appendix B under Assumption B.1. As this proof shows, the statement of the proposition remains true if the expectation over t is performed using any probability density ω(t) > 0, which may prove useful in practice. A similar result is stated in in  using the linear interpolant x t = x 0 (1 − t) + x 1 t. We now describe some primary facts resulting from this proposition, itemized for clarity: • The objective G(v) is given in terms of an expectation that is amenable to empirical estimation given samples t, x 0 , and x 1 drawn from ρ 0 , ρ 1 and U ([0, 1]). Below, we will exploit this property to propose a numerical scheme to perform the minimization of G(v). • While the minimizer of the objective G(v) is not available analytically in general, a notable exception is when ρ 0 and ρ 1 are Gaussian mixture densities and we use the trigonometric interpolant (5) or generalization thereof, as discussed in Appendix C.
• Minimizing G(v) over gradient fieldsv t (x) = ∇φ t (x) guarantees that the minimizer does not contain any componentṽ t (x) such that ∇ · (ṽ t (x)ρ t (x)) = 0. Such a component of the velocity has no effect on the evolution of ρ t (x) but affects the map X t , the solution of (2). We stress however that: even if we minimize G(v) over velocities that are not gradient fields, the minimizer v t (x) produces a map X t via (2) that satisfies the pushforward condition in (1). • The minimal value in (10) achieved by the objective implies that a necessary (albeit not sufficient) (11) In our numerical experiments we will monitor this quantity. This minimal value also suggests to maximize G(v) = minv G(v) with respect to additional control parameters (e.g. the interpolant) to shorten the W 2 length of the path {ρ t (x) : t ∈ [0, 1]}. In Section 2.2, we show that this procedure achieves optimal transport under minimal assumptions. • The last bound in (10), which is proven in in Lemma B.2, implies that the path length is always finite, even if it is not the shortest possible.
Let us now provide some intuitive derivation of the statements of Proposition 1: Continuity equation. By definition of the stochastic interpolant x t we can express its density ρ t (x) using the Dirac delta distribution as Since I t=0 (x 0 , x 1 ) = x 0 and I t=1 (x 0 , x 1 ) = x 1 by definition, we have ρ t=0 = ρ 0 and ρ t=1 = ρ 1 , which means that ρ t (x) satisfies the boundary conditions at t = 0, 1 in (3). Differentiating (12) in time using the chain rule gives where we defined the probability current Therefore if we introduce the velocity we see that we can write (13) as the continuity equation in (3).
Variational formulation. Using the expressions in (12) and (14) for ρ t (x) and j t (x) shows that we can write the objective (9) as Since ρ t (x) and j t (x) have the same support, the minimizer of this quadratic objective is unique for all (t, x) where ρ t (x) > 0 and given by (15).
Minimum value of the objective. Let v t (x) be given by (15) and consider the alternative objective where we expanded the square and used the identity v t (x)ρ t (x) = j t (x) to get the second equality. The objective G(v) given in (16) can be written in term of H(v) as The equality in (10) follows by evaluating (18) Gradient field velocity. Another way to write (13) as a continuity equation is to define the potential φ t : R d → R as the solution to with ρ t (x) and j t (x) given by (12) and (14), respectively. This is a Poisson equation for φ t (x) which has a unique (up to a constant) solution by the Fredholm alternative since ρ t (x) and j t (x) have the same support and R d ∇ · j t (x)dx = 0. In terms of φ t (x), (13) can therefore be written as The velocity ∇φ t (x) is also the unique minimizer of the objective (9) over gradient fields since if we setv t (x) = ∇φ t (x) and optimize G(∇φ) overφ, the Euler-Lagrange equation for the minimizer is precisely the Poisson equation (19). The lower bound on the objective evaluated at ∇φ t (x) still holds since in the argument above involving (17) and (18) can be made by replacing the (15) and φ t (x) defined as the solution to (19), we have consistent with the fact that the gradient field ∇φ t (x) is the velocity that minimizes given by (12).
Optimality gap. The argument above shows that we can use as alternative objective to G(v) since their first variations coincide. However, it should be stressed that this quadratic objective remains strictly positive atv = v in general so it offers no baseline measure of convergence. To see why, complete the square in G(v) to write (18) as where we used the shorthand notation However we stress that this inequality is not saturated We will come back to this inequality in Section 2.2 when positioning our approach in the context of optimal transport.

WASSERSTEIN BOUNDS
The following result shows that the objective in (17) controls the Wasserstein distance between the target density ρ 1 and the the densityρ 1 obtained as the pushforwrad of the base density ρ 0 by the mapX t=1 associated with the velocityv t : Proposition 2. Let ρ t (x) be the exact interpolant density defined in (B.5) and, given a velocity field v t (x), let us defineρ t (x) as the solution of the initial value problem The proof of Proposition 2 is given in Appendix D: it leverages the following bound on the square of W-2 distance where X t is the flow map solution of (2) with the exact v t (x) defined in (15) andX t is the flow map obtained by solving (2) with v t (x) replaced byv t (x). Note that the proposition suggests to regularize G(v) using e.g.
with some small λ > 0. In the numerical results presented below no such regularization was included.

OPTIMIZING THE TRANSPORT
Using the velocity v t (x) in (15) that minimizes the objective (9) in the ODE (2) gives an exact transport map T = X t=1 from ρ 0 to ρ 1 . However, this map is not optimal in general, in the sense that it does not minimize over allT such thatT ρ 0 = ρ 1 . It is easy to understand why: In their seminal paper Benamou & Brenier (2000) showed that finding the optimal map requires solving the minimization problem As also shown in (Benamou & Brenier, 2000), the velocity minimizing (30) In contrast, in our construction the interpolant density ρ t (x) is fixed by the choice of interpolant (9) is not the minimum in (30).
Minimizing (9) over gradient fields reduces the value of the objective in (30), but it does not yield an optimal map either-indeed the gradient velocity ∇φ t (x) with the potential φ t (x) solution to (19) only minimizes the objective in (30) over allv t (x) withρ t (x) = ρ t (x) fixed, as explained after (21). It is natural to ask whether our stochastic interpolant construction can be amended or generalized to derive optimal maps. This question is discussed next, from two complementary perspectives: via optimization of the interpolant I t (x 0 , x 1 ), and/or via optimization of the base density ρ 0 , assuming that we have some leeway to choose this density.
Optimizing the interpolant. Since (10) indicates that the minimum of G(v) is lower bounded by minus the value of the objective in (30), one way to optimize the transport is to maximize this minimum over the interpolant. Under some assumption on the Benamou-Brenier density ρ * t (x) solution of (31), this procedure works, as we show now. Let us begin with a definition: Definition 1 (Interpolable density). We say that one-parameter family of probability densities ρ t (x) with t ∈ [0, 1] is interpolable (in short: ρ t (x) is interpolable) if there exists a one-parameter family of invertible maps T t : R d → R d with t ∈ [0, 1], continuously-differentiable in time and space, such that ρ t is the pushforward by T t of the Gaussian density with mean zero and covariance identity, i.e.
Interpolable densities form a wide class, as discussed e.g. in Mikulincer & Shenfeld (2022). These densities also are the ones that can be learned by score-based diffusion modeling discussed in Section 2.3. They are useful for our purpose because of the following result showing that any interpolable density can be represented as an interpolant density: Proposition 3. Let ρ t (x) be an interpolable density in the sense of Definition 1. Then satisfies (4) and is such that the stochastic interpolant defined in (6) satisfies x t ∼ ρ t .
The proof of this proposition is given in E. We stress that the interpolant in (32) is in general not the only one giving the interpolable density ρ t (x), and the actual value of the map T t plays no role in the results below. Proposition 3 implies that: Proposition 4. Assume that (i) the optimal density function ρ * t (x) minimizing (30) is interpolable and (ii) (31) has classical solution. Consider the max-min problem where G(v) is the objective in (9) and the maximum is taken over interpolants satisfying (4). Then a maximizer of (33) exists, and any maximizer I * t (x 0 , x 1 ) is such that the probability density function of The proof of Proposition 4 is given in Appendix E. It relies on Proposition 3 as well as Lemma E.1 that reformulates (30) in a way that shows that this problem is equivalent to the max-min problem in (33) for interpolable densities.
Since our primary aim here is to construct a map T = X t=1 that pushes forward ρ 0 onto ρ 1 , but not necessarily to identify the optimal one, we can perform the maximization overÎ t (x 0 , x 1 ) in a restricted class (though of course the corresponding map is no longer optimal in that case). We investigate this option in the numerical examples below, using interpolants of the type (B.2) and maximizing over the functions a t and b t , subject to a 0 = b 1 = 1, a 1 = b 0 = 0. In Section 2.4 we also discussion generalizations of the interpolant that can render the optimization of the transport easier to perform. We leave the full investigation of the consequences of Proposition 4 for future work.
Optimizing the base density. Wile our construction allows one to connect any pair of densities ρ 0 and ρ 1 , the typical situation of interest is when ρ 1 is a complex target density and we wish to construct a generative model for this density by transporting samples from a simple base density ρ 0 . In this case it is natural to adjust parameters in ρ 0 to optimize the transport via maximization of G(v) = minv G(v) over these parameters-indeed changing ρ 0 also affects the stochastic interpolant x t defined in (6), and hence both the interpolant density ρ t (x) and the velocity v t (x). For example we can take ρ 0 to be a Gaussian with mean m ∈ R d and covariance C ∈ R d × R d , and maximize G(v) = minv G(v) over m and C. This construction is tested in the numerical examples treated below. We also discuss how to generalize it in Section 2.4.
Optimizing ρ 0 only makes practical sense if we we do so in a restricted class, like that of Gaussian densities that we just discussed. Still, we may wonder whether optimizing ρ 0 over all densities would automatically give ρ 0 = ρ 1 and v t = 0. If the interpolant is fixed, the answer is no, in general. Indeed, even if we set ρ 0 = ρ 1 , the interpolant density will still evolve, i.e. ρ t = ρ 0 except at t = 0, 1, in general. This indicates that optimizing the interpolant in concert with ρ 0 is necessary if we want to optimize the transport.
Recifying the map. In Liu, 2022), where an approach similar to ours using the linear interpolant x t = x 0 (1 − t) + x 1 t is introduced, an alternative procedure is proposed to optimize the transport. Specifically, it is suggested to rectify the map T = X t=1 learned by repeating the procedure using the new interpolant x t = x 0 (1 − t) + T (x 0 )t with x 0 ∼ ρ 0 . As shown in (Liu, 2022) iterating on this rectification step yields successive maps that are getting closer to optimality. The main drawback of this approach is that it requires each of these maps (including the first one) to be learned exactly, i.e. we must have T ρ 0 = ρ 1 to use the interpolant x t above. If the maps are not exact, which is unavoidable in practice, the procedure introduces a bias whose amplitude will grow with the iterations, leading to instability (as noted in ; Liu (2022)).

LINK WITH SCORE-BASED GENERATIVE MODELS
The following results shows that if ρ 0 is a Gaussian density, the velocity v t (x) can be related to the score of the density ρ t (x): Proposition 5. Assume that the base density ρ 0 (x) is a standard Gaussian density N (0, Id) and suppose that the interpolant I t (x 0 , x 1 ) is given by (5). Then the score ∇ log ρ t (x) is related to velocity v t (x) as The proof of this proposition is given in Appendix F. The first formula for t ∈ [0, 1) is based on a direct calculation using Gaussian integration by parts; the second formula at t = 1 is obtained by taking the limit of the first using v t=1 (x) = 0 from (B.20) and l'Hôpital's rule. It shows that we can in principle resample ρ t at any t ∈ [0, 1] using the stochastic differential equation in artificial time τ whose drift is the scoreŝ t (x) obtained by evaluating (34) on the estimatedv t (x): Similarly, the scoreŝ t (x) estimated this way could in principle be used in score-based diffusion models, as explained in Appendix F. We stress however that while our velocity is well-behaved for all times in [0, 1], as shown in (B.20), the drift and diffusion coefficient in the associated SDE are singular at t = 0, 1.

GENERALIZED INTERPOLANTS
Our construction can be easily generalized in various ways, e.g. by making the interpolant depend on additional latent variables to be averaged upon. This enlarges the class of interpolant density we can construct, which may prove useful to get simpler (or more optimal) velocity fields in the continuity equation (3). Let us consider one specific generalization of this type: Factorized interpolants. Suppose that we decompose both ρ 0 and ρ 1 as where K ∈ N, ρ k 0 and ρ k 1 are normalized PDF for each k = 1, . . . , K, and p k > 0 with K k=1 p k = 1. We can then introduce K interpolants I k t (x 0 , x 1 ) with k = 1, . . . , K, each satisfying (4), and define the stochastic interpolant This corresponds to splitting the samples from ρ 0 and ρ 1 into K (soft) clusters, and only interpolating between samples in cluster k in ρ 0 and samples in cluster k in ρ 1 . This clustering can either be done beforehand, based on some prior information we may have about ρ 0 and ρ 1 , or be learned (more on this point below).
It is easy to see that the PDF ρ t of x t is formally given by and that this density satisfies the continuity equation (B.11) for the current Therefore this equation can be written as (3) with the velocity v t (x) which is the unique minimizer of a generalization of the objective (9). We sate this as: Proposition 6. The stochastic interpolant x t defined in (37) with I k t (x 0 , x 1 ) satisfying (4) for each k = 1, . . . , K has a probability density ρ t (x) that satisfies the continuity equation (3) with a velocity v t (x) which is the unique minimizer overv t (x) of the objective where we used the shorthand notations I k t = I k t (x 0 , x 1 ) and ∂ t I k t = ∂ t I k t (x 0 , x 1 ). In addition the minimum value of this objective is given by and both these statements remain true if G(v) is minimized over velocities that are gradient fields, in which case the minimizer is of the form v t (x) = ∇φ t (x) for some potential φ t : R d → R uniquely defined up to a constant.
We omit the proof of this proposition since it is similar to the one of Proposition 1. The advantage of this construction is that it gives us the option to make the transport more optimal by maximizing G K (v) over I k t and/or the partitioning used to define ρ k 0 , ρ k 1 , and p k . For example, if we know that the target density ρ 1 has K clusters with relative mass p k , we can define ρ 0 as a Gaussian mixture with K modes, set the weight of mode k to p k , and maximize G K (v) = minv G K (v) over the mean m k ∈ R d and and the covariance C k ∈ R d × R d of each mode k = 1, . . . , K in the mixture density ρ 0 .

PRACTICAL IMPLEMENTATION AND NUMERICAL EXPERIMENTS
The objective detailed in Section 2 is amenable to efficient empirical estimation, which we utilize to experimentally validate the method. Let n samples from the target density ρ 1 , and {t k } K k=1 K samples from the uniform density on [0, 1]. Then an empirical estimate of the objective function in (9) is given by This calculation is fortunately parallelizable. Moreover, it is appealing to consider a neural network parameterization of the velocity field. In this case, the parameters of the modelv can be optimized through stochastic gradient descent (SGD) or its variants, namely the Adam optimizer (Kingma & Ba, 2015). We take this approach in the following experiments to demonstrate the validity and efficacy of the method using very simple network parameterizations, the purpose of which is to state that sensitive fine tuning is not necessary to demonstrate efficacy.
Following the literature in the advent of new recent density estimation methods, we benchmark the method on visualizable yet complicated densities that display multimodality, as well as higher dimensional tabular data initially provided in (Papamakarios et al., 2017) and tested in other works such as (Grathwohl et al., 2019). We take the 2D test case to demonstrate the ability to flow between empirical densities with no known analytic form while also maintaining a benchmark of convergence. In both instances, we compare both the computational cost as well as algorithmic efficiency to the conventional means of training ODE flows and show orders of magnitude speed up.
In all cases, numerical integration for sampling is done with the Dormand-Prince, explicit Runge-Kutta of order (4)5 (Dormand & Prince, 1980). In Sections 3.1-3.4 the choice of interpolant for experimentation was selected to be that of (5), as it is the one used to draw connections to the technique of score based diffusions in Proposition 5. In Section 3.5 we use the interpolant (B.2) and optimize a t and b t to investigate the possibility and impact of shortening the path length.

2D DENSITY ESTIMATION
An intuitive first test to benchmark the validity of the method is sampling a target density whose analytic form is known or whose density can be visualized for comparison. To this end, we follow the practice of choosing a few complicated 2-dimensional toy datasets, namely those from (Grathwohl et al., 2019), which were selected to differentiate the flexibility of continuous flows from discrete time flows, which cannot fully separate the modes. We consider anisotropic curved densities, a mixture of 8 separated Gaussians, and a checkerboard density.
The velocity field of the interpolant flow is parameterized by a simple feed forward neural network with ReLU Nair & Hinton (2010) activation functions. The network for each model has 3 layers, each of width equal to 256 hidden units. Adam optimizer is performed on G(v) for 10k epochs.
The results are presented in Figure 2. We plot a kernel density estimate over 80k samples from the flow as well as from the true distribution for comparison. The bounds from (11)

Flow True
Figure 2: Comparison of density estimation on highly multimodal, structured 2-D distributions.
As described in Section 2, the velocity field associated to the flow can be inferred from arbitrary densities ρ 0 , ρ 1 -this deviates from the score-based diffusion perspective, in which one distribution must be taken to be Gaussian for the training paradigm to be tractable.
In Figure 3, we illustrate this capacity by learning the velocity field connecting the anisotropic swirls distribution to that of the checkerboard. The interpolant formulation allows us to draw samples from ρ t at any time t ∈ [0, 1], which we exploit to check that the velocity field is empirically correct at all times on the interval, rather than just at the end points. This aspect of interpolants is also noted in (Choi et al., 2022), but for the purpose of density ratio estimation. The above observation highlights an intrinsic difference of the proposed method compared to MLE training of flows, where the map that is the minimizer of G(v) is not empirically known. We stress that query access to ρ 0 or ρ 1 is not needed to use our interpolation procedure since it only uses samples from these densities.

TABULAR DATA FOR HIGHER DIMENSIONAL TESTING
A set of tabular datasets introduced by (Papamakarios et al., 2017) has served as a consistent test bed for demonstrating flow-based sampling and its associated density estimation capabilities. We continue that practice here to provide a benchmark of the method on models which provide an exact likelihood, separating and comparing to exemplary discrete and continuous flows: MADE (Germain et al., 2015), Real NVP (Dinh et al., 2017), Convex Potential Flows (CPF) (Huang et al., 2021), Neural Spline Flows (NSP) Durkan et al. (2019), Free-form continuous flows (FFJORD) (Grathwohl et al., 2019), and OT-Flow (Finlay et al., 2020). Our primary point of comparison is to other continuous time models, so we sequester them in benchmarking.
We train the interpolant flow model on each target dataset listed in Table 1, choosing the reference distribution of the interpolant ρ 0 to be a Gaussian density with mean zero and variance I d , where d is the data dimension. The architectures and hyperparameters are given in Appendix G. We highlight  Figure 3: Demonstration of learning a flow map between distributions when neither have analytically known densities. Samples from the intermediate density are still available by the interpolant function x t = I t (x 0 , x 1 ), emitting an empirical diagnostic. some of the main characteristics of the models here. All models are feed-forward ReLU networks, except for the network trained on the BSDS300 dataset, which used ELU activations (Clevert et al., 2016). In each case, sampling of the time t was reweighted according to a beta distribution, with parameters α, β provided in the same appendix. These were selected to sample times closer to the t = 1 endpoint of the transport to ensure that the velocity at times close to the target is well-inferred.
Results from the tabular experiments are displayed in Table 1, in which the negative log-likelihood averaged over a test set of held out data is computed. We note that the interpolant flow achieves better or equivalent held out likelihoods on all ODE based models, except BSDS300, in which the FFJORD outperforms the interpolant by ∼ 0.6%. We note upwards of 30% improvements compared to baselines. We observed during model selection that increasing batch size generally improved performance, but this was not consistent across datasets, and the improvements were marginal. For example, a large target batch size proved useful on the POWER dataset, but the same did not hold for the BSDS300 dataset, for which there was no correlation between batch size and performance above a minimal threshold. Note that these likelihoods are achieved without direct optimization it.

COMPUTATIONAL EFFICIENCY AND CONVERGENCE
In the following, we show that the results achieved in the previous section are driven by a model that can train significantly more efficiently than the maximum likelihood approach to ODEs. We also provide an illustration of the convergence requirements on the objective defined in (11), which helps track that the learned model is meeting necessary conditions of the minimizer of G(v).   Figure 4 shows that, overall, the cost per iteration on a batch size of 50 is constant for the interpolant flow, while it grows for MLE based approaches as the ODE gets more complicated to solve. Moreover, the two examples show a 40x and 400x speedup for the interpolant flow as compared to FFJORD, suggesting continued gains in higher dimensions. The right figure shows that, under the same training conditions (optimizer, batch size, and learning rate), the interpolant flow can converge faster in number of training steps, in addition to being cheaper per step. The extent of this benefit is dependent on the set of hyperparameters necessary to achieve good performance, as well as the dataset, so a general statement about convergence speeds is difficult to make. For the sake of this comparison, we averaged over 5 trials for each model and dataset, the variance of which is shown shaded around the curves.
As described in Section 1, the minimum of G(v) is bounded by the square of the path taken by the map X t (x). The shifted value of the objectiveG(v) in (11) can be tracked to ensure that the model velocityv t (x) meets the requirement of the objective. It must be the case thatG(v) = 0 ifv t (x) is taken to be the minimizer of G(v), so we can look for this signature during the training of the interpolant flow. Figure 5 displays this phenomenon for an interpolant flow trained on the POWER dataset. Here, the shifted loss converges to the targetG(v) = 0 and remains there throughout training. This suggests that the dynamics of the stochastic optimization of G(v) are dual to the squared path length of the map.

PATH OPTIMALITY
As discussed in Section 2.2, the minimizer of the objective in equation (10) can be maximized with respect to the interpolant I t as a route toward optimal transport. We motivate this by choosing a parametric class for the interpolant I t and demonstrating that the max-min optimization in equation (33) can give rise to velocity fields which are easier to learn, resulting in better likelihood estimation. The 2D checkerboard example is an appealing test case because the transport is nontrivial. In this case, we train the same flow as in Section 3.1, with and without optimizing the interpolant. We choose a simple parametric class for the interpolant given by a Fourier series expansion This is but one set of possible parameterizations. Another, for example, could be a rational quadratic spline, so that the endpoints forâ t ,b t can be properly constrained as they are above in the Fourier expansion. In Figure 6, we plot the log likelihood, the learned interpolantsâ t ,b t compared to their initializations, as well as the path length as it evolves over training epochs. With the learned interpolants, the path length is reduced, and the resultant velocity field under the same number of training epochs endows a model with better likelihood estimation, as shown in the left plot of the figure. This suggests that minimizing the transport cost can create models which are, practically, objective function. This objective has an easily accessible diagnostic which can verify whether a proposed minimizer of the loss is a valid minimizer.
While our method is computationally more efficient, the new objective introduces an additional sum over time slices, which increases memory requirements. However, as shown in the experiments, these batch sizes need not be prohibitively large when stochastically optimizing G(v). This additional time sum also appears in score-based diffusion approaches and demonstrably poses no issue.
One salient feature of the proposed method is that choosing an interpolant I t (x 0 , x 1 ) decouples the optimization problem from that of also choosing a transport path, which is implicit in the MLE approach. This separation offers the opportunity of having explicit control on both. In future work, we will investigate more thoroughly how to optimize the transport by learning the interpolant I t in addition to minimizing G(v), as discussed in Secs. 2.2 and 2.4. It would also be interesting to approach this question from the perspective of path regularization as suggested by (28).
The intrinsic connection to score-based diffusion presented in Proposition 5 may be fruitful ground for understanding the benefits and tradeoffs of SDE vs ODE approaches to generative modeling. Exploring this relation is already underway (Lu et al., 2022;Boffi & Vanden-Eijnden, 2022), and can hopefully provide theoretical insight into designing more effective models.
where X s,t is the probability flow solution to In addition, given any test function φ : Ω → R, we have In words, Lemma A.1 states that an evaluation of the PDF ρ t at a given point x may be obtained by evolving the probability flow equation (2) backwards to some earlier time t to find the point x that evolves to x at time t, assuming that ρ t (x ) is available. In particular, for t = 0, we obtain and Proof. The assumed C 1 and globally Lipschitz conditions on v t guarantee global existence (on t ≥ 0) and uniqueness of the solution to (2). Differentiating ρ t (X t ,t (x)) with respect to t and using (2) and (A.1) we deduce Integrating this equation in t from t = t to t = t gives Evaluating this expression at x = X t,t (x) and using the group properties (i) X t ,t (X t,t (x)) = x and (ii) X t ,s (X t,t (x)) = X t,s (x) gives (A.2). Equation (A.4) can be derived by using (A.2) to express ρ t (x) in the integral at the left hand-side, changing integration variable x → X t ,t (x) and noting that the factor exp − t t ∇ · v s (X t,s (x))ds is precisely the Jacobian of this change of variable. The result is the integral at the right hand-side of (A.4).

B PROOF OF PROPOSITION 1
We will work under the following assumption: Assumption B.1. The densities ρ 0 (x) and ρ 1 (x) are continuously differentiable in x; I t (x 0 , x 1 ) is continuously differentiable in (t, x 0 , x 1 ) and satisfies (4) and (7); and for all t ∈ [0, 1] we have This structural assumption guarantees that the stochastic interpolant x t has a probability density and a probability current. As shown in Appendix C it is satisfied e.g. if ρ 0 and ρ 1 are Gaussian mixture densities and where a t and b t are C 1 function of t ∈ [0, 1] satisfyinġ The interpolant (4) is in this class for the choice Our proof of Proposition 1 will rely on the following result that quantifies the probability density and the probability current of the stochastic interpolant x t defined in (6): Lemma B.1. If Assumption B.1 holds, then the stochastic interpolant x t defined in (6) has a probability density function ρ t (x) given by and it satisfies the continuity equation withe the probability current j t (x) given by In addition the action of ρ t (x) and j t (x) against any test function φ : R d → R can be expressed as Note that (B.8) and (B.9) can be formally rewritten as (12) and (14) using the Dirac delta distribution.
Proof. By definition of x t in (6), the characteristic function of this random variable is Under Assumption B.1, the Fourier inversion theorem implies that x t has a density ρ t (x) given by (B.5). Taking the time derivative of this density gives with j t (x) given by (B.7).
Lemma B.1 shows that ρ t (x) satisfies the continuity equation (3) with the velocity field v t (x) defined in (15). It also show that the objective function G(v) in (9) is well-defined, which implies the first part Proposition 1.
For the second part, we will need Lemma B.2. If Assumption B.1 holds, then Then, using the pointwise identity v t (x)ρ t (x) = j t (x) as well as (B.8) and (B.9), we can write (B.14) where we use the shorthand I t = I t (x 0 , x 1 ) and ∂ t I t = ∂ t I t (x 0 , x 1 ). Therefore and this quantity is finite by assumption we deduce that lim K→∞ E φ K t (I t )|v t (I t )| 2 exists and is bounded by E |∂ t I t | 2 .
Lemma B.2 implies that the objective H(v) in (17) is well-defined, and the second part of statement of Proposition 1 follows from the argument given after (17). For the third part of the proposition we can then proceed as explained in main text, starting from the Poisson equation (19). Remark B.3. Let us show on a simple example that the inequality E |v t (I t (x 0 , x 1 ))| 2 ≤ E |∂ t I t (x 0 , x 1 )| 2 is not saturated in general. Assume that ρ 0 (x) is a Gaussian density of mean zero and variance one, and ρ 1 (x) a Gaussian density of mean m ∈ R and variance one. In this case, if we use the trigonometric interpolant (5), (C.4) indicates that ρ t (x) is a Gaussian density with mean sin( 1 2 πt)m and variance one, and (C.6) simplifies to v t (x) = 1 4 π 2 m 2 cos 2 ( 1 2 πt) , so that R |v t (x)| 2 ρ t (x)dx = 1 4 π 2 m 2 cos 2 ( 1 2 πt) (B.16) At the same time
Notice that we can evaluate v t=0 (x) and v t=1 (x) more explicitly. For example, with the linear interpolant (B.2) we have For the trigonometric interpolant that uses (B.4) these reduce to

C THE CASE OF GAUSSIAN MIXTURE DENSITIES
Here we consider the case where ρ 0 and ρ 1 are both Gaussian mixture densities. We denote the Gaussian probability density with mean vector m ∈ R d and positive-definite symmetric covariance matrix C = C T ∈ R d × R d . We assume that positivedefinite, and similarly for p 1 i , m 1 i , and C 1 i . We assume that the interpolant is of the form (B.2) and we denote Note that if all the covariance matrices are the same, C 0 i = C 1 j = C, with the trigonometric interpolant in (5) we have C ij t = C, which justifies this choice of interpolant. We have: Proposition C.1. The interpolant density ρ t (x) obtained by connecting the probability densities in (C.2) using the linear interpolant (B.2) is given by and it satisfies the continuity equation ∂ t ρ t (x) + ∇ · j t (x) = 0 with the current To use the bound in (27), let us consider the evolution of Therefore, by Gronwall's inequality and since Q 0 = 0 we deduce Since W 2 2 (ρ 1 ,ρ 1 ) ≤ Q 1 by (27) so that the boundary condition in (4) are satisfied. Then observe that, by definition of T t , This implies that, if x 0 ∼ ρ 0 , x 1 ∼ ρ 1 , and they are independent, T −1 0 (x 0 ) cos( 1 2 πt) + T −1 1 (x 1 ) sin( 1 2 πt) (E.3) is a Gaussian random variable with mean zero and covariance Id cos 2 ( 1 2 πt) + Id sin 2 ( 1 2 πt) = Id, (E.4) i.e. it is a sample from N (0, Id). By definition of T t , this implies that x t = I t (x 0 , x 1 ) = T t T −1 0 (x 0 ) cos( 1 2 πt) + T −1 1 (x 1 ) sin( 1 2 πt) (E.5) is a sample from ρ t if x 0 ∼ ρ 0 , x 1 ∼ ρ 1 , and they are independent.
The proof of Proposition 4 relies on the following simple reformulation of (30): Lemma E.1. The Benamou-Brenier minimization problem in (30) is equivalent to the min-max problem subject to: ∂ tρt + ∇ · t = 0,ρ t=0 = ρ 0 ,ρ t=1 = ρ 1 (E.6) In particular, under the conditions on ρ 0 and ρ 1 such that (30) has a minimizer, the optimizer (ρ * Solving this expression in ∇ log ρ t (x) and specializing it to the trigonometric interpolant with a t , b t given in (B.4) gives the first equation in (34). The second one can be obtained by taking the limit of this first equation using v t=1 (x) = 0 from (B.20) and l'Hôpital's rule.
Note that (F.3) shows that, when the interpolant is of the type (B.2) and ρ 0 (x 0 ) = (2π) −d/2 e − 1 2 |x0| 2 , the continuity equation (3) can also be written as the diffusion equation Since we assume thatȧ t ≤ 0 andḃ t ≥ 0 (see (B.3), the diffusion coefficient in this equation is negative This means that (F.5) is well-posed backward in time, i.e. it corresponds to backward diffusion from ρ t=1 = ρ 1 to ρ t=0 = ρ 0 = (2π) −d/2 e − 1 2 |x0| 2 . Therefore, reversing this backward diffusion, similar to what is done in score-based diffusion models, gives an SDE that transforms samples from ρ 0 to ρ 1 . Interestingly, these forward and backward diffusion processes arise on the finite time interval t ∈ [0, 1]; notice however that both the drift and the diffusion coefficient are singular at t = 1. This is unlike the velocity v t (x) which is finite at t = 0, 1 and is given by (B.19  The architectural information and hyperparameters of the models for the resulting likelihoods in Table 1 is presented in Table 2. ReLU (Nair & Hinton, 2010) activations were used throughout, barring the BSDS300 dataset, where ELU (Clevert et al., 2016) was used. Table formatting based on (Durkan et al., 2019).
In addition, reweighting of the uniform sampling of time values in the empirical calculation of (42) was done using a Beta distribution under the heuristic that the flow should be well trained near the target. This is in line with the statements under Proposition 1 that any weight w(t) maintains the same minimizer.
G.1 PATH OPTIMALITY AND PARAMETERIZING ρ 0 For the path optimality experiments, M = 7 Fourier coefficients were used to parameterize the interpolant.
In addition to parameterizing the interpolant, one can also parameterize the base density ρ 0 in some simple parametric class. We show that including this in the min-max optimization given in equation (33) can further reduce the transport cost and improve the final log likelihood. The results are given in Figure G.1. We train an interpolant just as described above, but also allow the base density ρ 0 to be parameterized as a Gaussian with meanμ and covarianceΣ. The inclusion of the learnable base density results in a significantly reduced path length, thereby bringing the model closer to optimal transport.