Lagrangian Flow Networks for Conservation Laws
Authors: Fabricio Arend Torres
Marcello Massimo Negri
Marco Inversi
Jonathan Aellen
Volker Roth
Venue: ICLR 2024 spotlight (2024)
Keywords: Physics-informed Neural Network
Fluid Dynamics
Conservation Law
Partial Differential Equation
Conditional Normalizing Flows
Bird-Migration
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: An approach for modeling fluid densities and velocities continuously in space and time while satisfying the continuity equation by construction.

Abstract: We introduce Lagrangian Flow Networks (LFlows) for modeling fluid densities and velocities continuously in space and time.
By construction, the proposed LFlows satisfy the continuity equation,
a PDE describing mass conservation in its differential form. 
Our model is based on the insight that solutions to the continuity equation can be expressed as
time-dependent density transformations via differentiable and invertible maps.
This follows from classical theory of the existence and uniqueness of Lagrangian flows for smooth vector fields.
Hence, we model fluid densities by transforming a base density with parameterized diffeomorphisms conditioned on time.
The key benefit compared to methods relying on numerical ODE solvers or PINNs is that the analytic expression of the velocity is always consistent with changes in density.
Furthermore, we require neither expensive numerical solvers, nor additional penalties to enforce the PDE.
LFlows show higher predictive accuracy in density modeling tasks compared to competing models in 2D and 3D,
while being computationally efficient.
As a real-world application, we model bird migration based on sparse weather radar measurements.

LAGRANGIAN FLOW NETWORKS FOR CONSERVATION LAWS
Fabricio Arend Torres, Marcello M. Negri, Marco Jonathan Aellen & Volker Roth Department of Mathematics and Computer Science
University of Basel
marcellomassiin megiim...nggi,,... marco. inversi, jonathan.aellen, volker roth woth
ceabricio.arenttooo
ABSTRACT
We introduce Lagrangian Flow Networks (LFlows) for modeling fluid densities and velocities continuously in space and time. By construction, the proposed LFlows satisfy the continuity equation, PDE describing mass conservation its differential form. Our model is based on the insight that solutions to the con- tinuity equation can expressed as time-dependent density transformations via differentiable and invertible maps. This follows from classical theory of the ex- istence and uniqueness of Lagrangian flows for smooth vector fields. Hence, we model fluid densities by transforming base density with parameterized diffeo- morphisms conditioned on time. The key benefit compared to methods relying numerical ODE solvers PINNs is that the analytic expression of the veloc- ity is always consistent with changes in density. Furthermore, we require neither expensive numerical solvers, nor additional penalties to enforce the PDE. LFlows show higher predictive accuracy in density modeling tasks compared to competing models 2D and 3D, while being computationally efficient. As a real-world ap- plication, we model bird migration based on sparse weather radar measurements.
INTRODUCTION
The development of physics-informed Machine Learning (PI-ML) (Karniadakis et al. (02T) opens new opportunities to combine the power of modern ML methods with physical constraints that serve as meaningful regularizers. These constraints might, for example, be available in the form of partial differential equations (PDEs). Within PI-ML consider hydrodynamic flow problems governed by the physical law of mass conservation. This law is described its local and differential form by PDE commonly known as the continuity equation (CE)
{ O.p V V. (vp) 0
x Ω (1, E (lo, x E
p(to, x) = (2))
For any time E [to, I') the function p(t,-) can be thought of as the density of parcels advected by the velocity field u, with initial density Pto Here, [to, T) x Ω is the space-time domain, which is
subset R The partial derivative w.r.t. time is denoted by θ and b V Di 0b.
dz.
is the spatial divergence of a d dimensional vector field b [to, 7 x Ω + Rd
Unlike classical initial value problems, we consider challenging settings where exact boundary and initial conditions are unknown. That is, the continuous density and velocity fields have to be inferred from sparse and noisy data. The important physical constraint is that the solution must comply with the CE in Eq. I To this end we propose neural network based model that fulfills the CE by construction and provides physically consistent velocity and density fields.
We specifically consider two distinct application settings. In both settings we are restricted to sparse and noisy measurements. setting (i) we measure the fluid density p and velocity v, without knowing any additional equations aside from the CE. This occurs for example within the area of radar ornithology (Chilson et al. 2017), where density and velocity of birds can be inferred from radar data. Such radar-based measurements to date the only practical high-throughput data

---

source for birds. The goal is a spatiotemporal estimate of bird migration densities (Nussbaumer et al. 2019). Since more migration-specifie dynamics are unknown, mass conservation is usually considered as a physical constraint (Nussbaumer et al., 2021 Lippert al.] 2022a b). In setting (ii) we exclusively have sparse density observations, but we know additional equations constraining the velocity. This might, for example, occur in dynamical optimal transport problems. Here, two densities are to be interpolated while adhering to the CE. Therefore, minimizing the total transport cost constrains the velocity field. More generally, setting (ii) includes wide range of compressible fluids dynamic problems. In both settings we are mainly interested in accurately modeling the density, with the velocity measurements or additional equations serving as an informative prior. We provide code at ptttps://gtt thub. com bmdaa Sonda-uniiia www.-uuiiaaaaagggaa www.baaaaaaaaaaggiaaattt
woNetworks
Main contributions. The main contributions of this paper are as follows:
We outline a fundamental link between densities modeled by conditional Normalizing Flows and spatiotemporal density fields that satisfy the CE.
We leverage this link to introduce models for ill-posed hydrodynamic flow problems that always satisfy the CE by construction, coined Lagrangian Flow Networks (LFlows). We do so without requiring an explicit representation of the initial density.
We provide way to calculate the velocity without inverting the conditional Normalizing Flow, enabling the use of flexible bijective layers with otherwise costly inverses.
We assess LFlows in multiple application settings and show better predictive performance than existing methods while staying computationally feasible and physically consistent.
RELATED WORK
Neural Networks for Conservation Laws. Physics-informed neural networks (PINNs) (Raissi et al. 2019) enforce PDEs in neural networks by introducing an additional PDE-loss that penalizes pointwise deviations from the PDE. The PDE-Ioss is enforced collocation points that are sampled in the signal domain. The accuracy of PINNs is thus limited by the amount (and distribution) of collocation points, as well the dimension of the signal domain. For conservation laws, recent improvements in PINNs use more sophisticated sampling approaches (Arend Torres et al. (2022), or introduce domain decompositions (Jagtap et al., 2020). Although this alleviates scaling prohlems, the fundamental limitation due to the number of collocation points (or subdomains) still remains.
In contrast, Richter--Poelle et al. (2022) propose parameterization of neural networks that enforces mass conservation by design, which we refer to as divergence-free Neural Networks (DFNNs). Solutions to the CE in Eq. I are represented as divergence-free (d + 1) dimensional vector field b pe) with augmented (d+11 1) dimensional input space = (t,
Op V2 (pv)
d+1 bb
= V. (P))
- bb
+
Σ 1=1 Os;
(2)
Ōt
The generalization of divergence-free vector fields to higher dimensions is achieved through the concept of differential forms. The resulting parameterization, however, heavily relies on expensive higher-order automatic differentiation, posing limitations in terms of scalability.
Concurrent to our work, Li al. (2023) (TIPF) proposes the use of Lagrangian flow maps to model continuous probability flows that fulfill the CE. Unlike LFlows, TIPF considers well-posed PDE settings with known initial conditions. Specifically, they focus on Fokker-Planck equations and Wasserstein gradient flows by additionally introducing an unbiased self-consistency loss. contrast, we consider ill-posed data-assimilation for physical problems with sparse and noisy data where initial conditions are unknown. Additionally, I.Flows stand out as they don't 't require the inversion of bijective layers for computing the velocity, allowing for more expressive bijections.
Deep operator learning (Lu et al. 2019 Liet et al. 2020) was proposed as a general approach to learn dynamics from dense observations. In these settings, the PDE is not provided, but has to be learned from large amounts of dense data often provided through simulations. As such, it is not applicable to setting with spatially sparse data obtained at irregular time steps. Further notable mentions Lagrangian and Hamiltonian neural networks (Cranmer et al, 2020; Greyddnus et al., 2019), which can learn conservation laws from observed trajectories of individual particles.
2

---

Data assimilation with the Adjoint. The adjoint method (Cacuci 1981a Pontryagin 1987
allows to differentiate through numerical solvers by integrating adjoint equations backward in time. By minimizing an objective function, it is then possible to infer the initial conditions and parameters of a dynamical system from data. Within adjoint methods, the well-established semi-Laggaanian data assimilation (SLDA) approach is the closest one conceptually
model and setting (Robert 1982 Staniforth & Côté. [191]] Diamantakis & Magnusson] (0166). SLDA is widely used for for integrating transport equations into atmospheric models (Diamantakis) 2013 Hersbach al. 2020).
It is based on the Lagrangian form of the CE, which can be written as an ODE:
<00, z) In p(0, 2)
] = In (n)
d x x(t, x) df In p(t, x(t,
]
vlt, (t, 2)) att, z)) }
(3)
where E Rd v: |4o, T) IR - R4, E TT) Given the initial position and (log-) density In Po(z) of the parcel, the ODE describes their temporal evolution according to the CE. The density at the departure points (i.e. the initial time) is represented by an interpolated mesh. The data-loss for the density is computed by mapping observations backwards in time with Eq. 3 and then comparing it to the interpolated initial density.
An efficient autograd implementation of the adjoint method was presented by Chen et al. (2018). The implementation enables black-box differentiation for numerically solved ODEs, and further allows to specify the dynamics of ODEs with a neural network. The introduced Continuous Normalizing Flows can be regarded the velocity with a neural network, while fixing the density at the departure points. A limiting factor of neural adjoint based methods is their computational cost, since input derivatives of the velocity v are repeatedly evaluated for every step of the solver. Furthermore, dynamics given by neural networks can become stiff during training. To avoid these issues Bilos et al. (2021) propose time- dependent bijections instead of neural ODEs for modeling time series data.
special case of the SLDA that propagates probability densities and models
(Conditional) Normalizing Flows Normalizing Flows (NFs) are general approach warping a simple probability distribution into more complex target distribution via invertible and differ- entiable transformations, i.e. diffeomorphisms. Let R Rdd be random variable with a known density function R ~ (R( and let Y = T(R), where T is a diffeomorphism with trainable pa- rameters. With change of variables the probability density of Y can be expressed in terms of the base density pR. the map T, and its Jacobian:
(y) = PR(T"(())) det 'T))
(4)
NFs usually rely on transformations for which the Jacobian determinant can be efficiently and easily calculated. A parameterization for conditional distributions py (y|c) can be obtained by additionally
conditioning the parameters of T on another variable e through hypernetwork (Ha et al. 2017 This is commonly called a conditional Normalizing Flow (Atanov et al. 2011:; Kobyzev et al.
2020
For a review of NFs, we refer to Kobyzev et al. (2020) and Papamakarios et al. (2021).
LAGRANGIAN FLOW NETWORKS
We first present some key results of the classical theory of Lagrangian flows for smooth vector fields. These provide us a framework for evolving densities and velocities that always fulfill the CE. We then propose parameterizations that result in simple expressions for the velocity and density. The resulting LFlow models densities and velocities by building upon conditional Normalizing Flows.
3.1 FLOW MAPS AND THE CONTINUITY EQUATION
The Lagrangian view describes fluids from the perspective of moving fluid parcels, i.e. infinitesimal volumes with constant mass. From this point of view the CE states that density changes of the fluid are described by volume changes of parcels. That is, spatial contraction increases the density of a parcel, and expansion decreases it. In order to compute the density of any parcel, we then only need to know its initial density, and how much its volume was distorted.
More formally, let Ito denote the initial position of a parcel at time to- In addition, let X Q + Ω for a fixed t [to. T] be a diffeomorphism that maps 2t to the parcel position time t:
X, (x,),,)) Ez,
(5)
3

---

Phase O (0,
9, (x)
X.
0.3
a(t.a)
6((.z)
(sr)
X.(x))
(+)
X
Figure 1: Illustration of the transformations and involved fields for modeling the temporal evolution of 2D density with LFlows. The red lines indicate trajectories of fluid parcels.
with Xto being the identity map. That is, X; provides the continuous trajectory of the parcel Eta- We further assume basic regularity of X, and X,' such as smoothness and globally bounded derivatives. Since X, provides the trajectory of parcel, the velocity of a given parcel at position æ; and time follows naturally. First, the parcel is mapped back to its initial position with X, The velocity is then the change in position along its trajectory with respect time l:
0X. v(l, (() at
(x,"(x)).
Following classical Cauchy Lipschitz theory, it is known that such a map X, is the unique flow map of a starting time to. Specifically, for any E Q the curve + X, (x) is the unique solution to the Cauchy Problem
0 (x) = v (1. (xx)) L [lo, T), i
A more complete statement is given in A.1.I Theorem and we refer to Hartman (2002) for an extensive description of the theory of ordinary differential equations.
With the velocity given by Eq. 6 we further need define density describe a fluid. Let Ptp Q + R, be the (known) initial fluid density at time ty. We can then define the time-evolved density as transformation of Ptg using the change of variables formula:
p(t,x) = Pso (x-'( (xx) |ett (x))-
(8)
Given the velocity in Eq. 6 and smooth Pto Eq. 8 is solution to the continuity equation (see A.1.I Theorem 2). Proofs Tor this statement vary significantly their complexity and depend on the regularity assumptions for the velocity. That is, they range from classical theory to current mathematical research. For the sake of completeness, we provide precise statements, assumptions and proofs in the Appendix Section A.1 The appended proofs based the well-established classical theory for the existence and uniqueness of Lagrangian flows for smooth vector fields and refer to Ambrosio & Crippa (0088) for basic and advanced results.
3.2 LAGRANGIAN FLOW NETWORKS
We can now exploit the derived connection between the CE and time-evolving diffeomorphisms to model densities and velocities that satisfy the CE by construction. Instead of directly parameterizing both X, and Ptg (as in et al. (2023)), we model the density at each time Pt. including Pt. as transformation of simple fixed density Poase- This requires only one cime-conditioned bijection 4,- We call the resulting model Lagrangian Flow Networks (LFlows), which we illustrate in Figure [
Let 4, Ω + nd be a learnable diffcomorphism with t [T. We propose to parameterize X; as the composition
x, $ "(Dr (x)).
(9)
4

---

In practice, implement L as an invertible neural network with its parameters conditioned on time. We further define the initial density as a transformation of a simple base density:
Pto (2))
Prase
Peo (ar)) |det J 1 (a))
(10)
The base density Prase R H R. is an unnormalized probability density:
Ponse (z) = c NN I), where e R2 is the total mass of the system and freely learnable parameter.
(11)
We now substitute X, in Eq. 8 with the parameterized X, from Eq. 2 We further substitute the density Pto with the parameterized Pto (Eq. I0). The modeled density then simplifies to
(t. Pro
{x (3x) |det JX, (ax) Phase (f.(2))
(12)
We refer to the Appendix A for the intermediate steps. Note that the resulting expression co- incides with the change of variable formula for probability densities in Eq. 4 This allows to elegantly model the evolving density through a conditional normalizing flow with unnormalized base density Phase and bijective layers conditioned on time +,.
The parameterization of X, with X, in Eq. also results in a simple expression for the velocity:
é(t, x)) ax.
x (x)) (J)4
0 (a)
(13)
We provide the explicit steps for arriving at Eq. [3] in the Appendix AA2222 Note that in order to evaluate ý and ē we now require only the forward map Φ;, but not inverse. This proves useful for layers with expensive inverses, or if the inverse is unknown. An illustration that unifies the
Lagrangian view and the provided parameterization of LFlows is given in Figure [
Limitations. LFlows model fluid densities and velocities by transforming a base density with bijective layers. As such, similar limitations as for Normalizing Flows apply. If the target density has disconnected modes, the base density must have the same number of disconnected modes due to topological constraints (Papamakarios et al.) 202T). If not, the space in between disconnected modes will be covered by small but non non-zero density. Furthermore, LFlows power of the bijective layers. Even though state-of- bijective layers are highly flexible, each layer might still be limited in terms of the number of modes that can be modeled (Liao & He 2021
limited by the expressive
4 IMPLEMENTATION
To implement LFlow as outlined in Section 3 we require conditional bijective layers. That is, the diffeomorphism o required for Eq. 12 and Eq. 13 has to be conditioned on time To allow for a flexible parameterization, we first embed into a higher-dimensional space with an embedding network and then condition the bijections on this embedding, i.e. Φ, :- ((:; fee (t)) with fe : (4. To & The k-dimensional embedding is shared between the individual layers. A high-level visualization of the network architecture provided in the Appendix Figure A.6 We implement
ft( as MLPs with residual skip connections and swish activations (Elfwing et al. 2018).
We implement flexible & (æ;; (t)) with Sipschitz-constrannce invertible densenets (Perugachi-De Diaz et al., 2021), which have free-form Jacobians. For a conditional variant of these layers we pass the embedding (t) as an additional input. That is, each invertible densenet layer is a function g(x, t) = x + h(#, (t)) with h Rd+ + Rd and Lip(h) < 1. The activations of h are sinusoidal (xx) = sin(wrr)u with w E R,, where the division by w ensures Lip(s) = 1. Invertible densenets provide numerical inversion via fixed-point iterations. In addition to densenet blocks, we employ (unconditional) intermediate activation normalizations (Kingma & Dhariwal, (0118) and (conditional) SVD layers. The orthogonal components of the SVD layers parameterized by conditional Householder reflections.
In all our settings the total mass, i.e. the normalization constant in Eq. [T] is not known. Therefore, we treat c as a freely learnable hyperparameter, which we initialize based on validation set. We further encourage solutions with small total mass during training with the penalty Luaes = w. c, where w, E Ry0 is a hyperparameter. In practice, this penalty discourages learning significant densities in areas where there measurements.
Normalization Constant.
5

---

Graundtruth
SLDA
Dereity R2
0.75
0.50
LFOOw
SLDA
DFIW
PANN
LFOW
SLD4
DFNN
PINN
120
Q1
SMAPE Consistencs
0.4
0.2
1=0000
0.0
LFOow
SLDA
PINN
LFOW
SLDA
PINN
Figure 2:
Evaluation of the predicted density on the zy-plane with 0 for different methods.
Arrows indicate velocity, except for the DFNNs, where the normalized flux shown. The lower left plot shows the spatial splitting into train (green), validation (yellow), and test (purple) subsets.
Baseline Models. For implementation details all competing methods, we refer to the Appendix Section A.3 and provided code.
5 EXPERIMENTS
We showcase LFlows for two distinct settings. In setting (i), density and velocity are observed, but no equations are available (Snetion In setting (ii) only the density is observed but further equations for the velocity known (Section 5.2). For details on the data, architecture, and code for each experiment, we refer to the Appendix A.4 to An6 the supplementary material. We will compare LFlows with methods that enforce the CE through different means, namely divergence-free neural networks (DFNNs), physics-informed informedd neural networks (PINNs), and -agrangggn data assimilation (SLDA). For details on the implementation we refer the Appendix A.3
Numerical Evaluation of Physical Consistency. Physical consistency in terms of the CE implies that the predicted density at any time coincides with the initial density transformed forward in time by the learned velocity field. An inconsistent model would imply that the predicted velocity field fundamentally disagrees with the predicted density movements. This would make any downstream interpretation of the two learned fields futile. We quantitatively evaluate this potential inconsistency
in the following experiments. We compare the predicted density (model (t, x) with the numerical
solution of the IVP defined by Pmodel((o, x). (t, and Eq. 3 We evaluate the symmetric (b..aa) (o
mean absolute percentage error sMAPE -. -
Prede (1,, [2. (c, that ranges from 0 to 1.
5.1 SIMULATED FLUID FLOW
For synthetic example of setting (i) we simulate densities in 2D and 3D over time by trans- forming a mixture of four unnormalized Gaussians. We parameterize time-dependent bijections in t E [0, 1.2], S (-4, 4)d, which provide analytical forms for the densities and velocities. During training, only sub-regions the domain Q are observed. The dynamics in 3D are similar to the 2D setting, with the zy-velocity being the same for all values and the velocity being 0. The only added difficulty is a higher-dimensional domain. We limit all models to the computing resources of NVIDIA Titan Pascal, and optimize based the explained variance (R22) (R2) of the density on the validation set. For the PINN this resulted in 26 collocation points. We do not include consistency results for DFNNs due to numerical instabilities. As DFNNs only provide access to the flux F. the velocity F/O becomes numerically unstable low-density regions, which are abundant in this experiment.
R
MSE
with R² indicating perfect reconstruction.
(Hoh)

---

Results for 10 random seeds provided in Figure In addition, snapshots of the 3D prediction for 2 = are shown. All methods aside from the PINN have a low consistency sMAPE on the order of le-4 or lower. This is expected, as LFlows enforce the CE by construction. Furthermore, SLDA computes the density similarly to our numerical reference, although with lower order ODE solver. Low error tolerances low-order solvers for SLLA would of course still result in inconsistencies. Looking at the predictive performance, LFlows show the highest average R² for the density in both 2D and 3D. While PINNs perform competitively in 2D, they severely degrade in 3D, the number of collocation points used (limited by GPU memory) is insufficient to enforce the PDE. This is also reflected their increase of the consistency sMAPE in 3D. Finally, DFNNs and SLDA are roughly comparable in terms of predictive accuracy, with DFNNS being in our experience most prone to overfitting. We note that small density displacements can already lead large differences in R².
5.2 DYNAMICAL OPTIMAL TRANSPORT
As example of setting (ii), in which no velocity is observed but additional equations dictating the dynamics are known, we consider dynamical optimal transport problems. Our experimental setting closely follows Richter-Powell al. (2022). Specifically, we consider the Kenamou--rentie formulation of the optimal transport problem between two densities pt, and Pt,, In this case the optimal transport map is the solution map X, of flow that is defined by the vector field and minimizes the following objective:
/ |u(t, x))) x)) dr dt
min 0.p
(14)
subject the constraints p(lo, .²) = Pto (a)) and p(11, x) = pt, (ar). Furthermore, p and e are subject to the continuity equation Oep = - (pr).
Both DFNNS and LFlows can solve the minimization problem Eq. without needing separate estimation of Instead, one fits the densities at to and f and additionally minimize Eq. [4] How- ever, to obtain the transport map from the learned velocity field, DFNNs need to numerically solve the Cauchy problem in Eq. 7 An ODE solver might however struggle due to the numerical instabil- ity of the DFNNs velocity in low-density regions. In contrast, LFlows elegantly provide an analytical form for the continuous transport map through the learned bijections, i.e. (x(x) = Φ, (ot (x:))
LFlow
6180
0.969
8.015
LFlow
Pinwheel Caussians
DFNN
0.20
2300
8240
0
Minimax
Gaussiens Circles
Discrete
Discrete
0.4
0.17
0.50
Figure 3: Left: Approximations of the 2D optimal transport map with LFlows, DFNN and a discrete reference. Right: Corresponding estimated Wasserstein distances for different methods for single run. The red vertical lines denote the minimum, median, and maximum estimates of 5 runs with discrete OT solver. Datasets: Circles < Pinwheel, Pinwheel 8Gaussians, 8Gaussians < Circles.
We train the models by optimizing
min XEg, [],(0, x) (=(=) + ΛEg [l)(1,x) p1(x))] ( ((tt,))) x) dx dt (15)
A.6
7

---

where data is drawn from pi, which is a mixture of p; and uniform density (for i = 0. 1); À is hyperparameter. At test time we empirically estimate the W distance by mapping 5000 samples of Po from 0 to t 1. Different to Pochter-Powell et al. (2022) repeat this estimate 50 times. We compare the 3 estimates of (i) LFlows, (ii) DFNNS and (iii) a minimax formulation of the optimal transport map learned via input convex neural network (Makkuva al. 2020 For
DFNNs the samples are transported through the estimated velocity with an ODE solver. We further compare to the W estimates of a discrete OT-solver (Bonneel et al. 2011) from the por library (Flamary et al. 2021 based on 50000 samples. The À for LFlows is chosen by matching the W2 distance between the moons and swissroll datasets with the discrete estimate. For DFNNs a À value is provided by Powwll et al. (2022). In this experiment we restrict ourselves to methods that can exactly enforce the CE and thus exclude PINNs.
We considered three different pairs of toy 2D distributions. Figure shows the approximated optimal transport maps learned with LFlows and the W3 estimates of the different methods. We verified that the DFNN and the LFlow fit the target densities well, with a test MSEs below 8e-5 for all settings. We excluded SLDA, as it was unstable and did not consistently result in low errors for the two target densities. The LFlow estimates of W? are closest to the range of discrete estimates (Figure 3). In contrast, the minimax model underestimates the distance which is consistent with the results of RRihhter-Powell et al. (2022). DFNNS significantly overestimated W and we were unable fully reproduce results obtained by |Rihher--Poell et al. (2022). We assume that this is due to the ODE solver struggling with the unstable velocity calculation in low-density regions.
5.3 MODELING BIRD MIGRATION
As real-world application for setting (i) we model bird migration within Europe based on weather radar measurements. The data provided by Nusssbumeer et al. (202T) contains estimated bird densi- taken from 37 weather radar stations in France, Germany, and the Netherlands at up to 5-minute intervals at 200m altitude bins, reaching up to 5km. The velocity data does not include z-axis component. We model the bird migration of 3 subsequent nights of April 2018. We assume that the mass is mostly conserved within the three nights during migration. We test the predictions on radars located in the center of the covered region, which were excluded during training (see Figure 4.. Hyperparameters of all models are selected by minimizing the density MSE on three nights of March 2018 As a baseline, we compare the re- sults with 10-layer multilayer perceptron (MLP) with skip connections, 256 hidden units per layer, ReLU activations, and Batch normalization. We additionally evaluate the PINN, DFNN, and SLDA. The PINN has the same general architecture as the MLP but additionally minimizes the penalty A- 1 + V (po) ll2 evaluated on 100 000 collocation points, where À E R., is a hyperparameter.
ties (birds/ (hirds/km²) and velocities (m/s). Measurements
54"
04 02/0.00
Binge
12
32'
60%
4479
4276
Figure 4: Snapshots of predicted bird density at three consecutive nights within central Europe. The 2D projection was obtained by integrating over altitudes covered by the radars. Orange lines indicate 2D (ry) projections of 3D trajectories from to to using randomly sampled departure points.
Figure 4 shows snapshots of the vertically integrated density and flux predicted by LFlows. Predic- tions of the other models are provided in the Appendix [A.6.2] Aside from the velocity and density, LFlows readily provide the trajectories (shown in orange). In practice, experts could compare these with the migration paths taken by individual birds, which are for example obtained from bird-ringing studies. We further explore the role of the total mass penalty w. in Appendix AA6.3
We compare test density errors for each model in Figure 5 (c). Methods enforcing the CE result in lower error than the baseline MLP. The consistency sMAPE shows the PINN and MLP lead to
8

---

inconsistent density and velocity. LFlows, SLDA and DFNNS again have sMAPE that is the order le-4 or lower. Figure S (a) and (b) show the pointwise sMAPE for the PINN and LFlow.
While DFNNs and SLDA nearly competitive with LFlows in terms of the density MSE, they suffer from high computational costs. SLDA requires to evaluate neural network many times to numerically solve the ODE. In addition, SLDA gets slower during training because the dynamics given by the neural network become more stiff. This is shown by the increasing runtime for each epoch in Figure 5 (e) and is known limitation of models of the neural ODE family (Biloš al. 2021). DFNNs on the other hand result in huge memory requirements due to the required second order derivatives. Fig shows the peak memory use in terms of VRAM for varying minibatch sizes of the models. In practice, high memory requirements result in small minibatch sizes, which ultimately lead to slower training and inference pipeline (Shallue al. 2018). The high peak memory and runtime of PINNs is due to the large amount of collocation points.
sMAPE Consistency
0.23
0.5
0.75
LFlow
SLDA
Secords
DFNN
(e)
40
PINN
20
MLP
0.4
0.6
0.8
1.0
100
200
300
MSE of loglp Density
Training Epoch
S.DM
LFlow
Poak Memory (GB)
30
SLDA
20
S
DFNN
PINN
20
MLP
0.0
0.2
0.4
0.6
256
1024
2048
3072
LFlow
SMAPE of Consistency
Minibatch Size
Results of the bird migration experiment. (a) Consistency sMAPE of PINNs and (b) LFlows at fixed time (2018-044 04-08 00:00). (c) Test MSE of the (loglp) density and (d) consistency sMAPE evaluated at multiple timesteps. (e) Time per training epoch during training and (Π) peak memory usage during training GB VRAM for varying minibatch sizes.
Figure 5:
6 CONCLUSION
We introduced LFlows for modeling densities and velocities that adhere the continuity equation by construction. We did so by establishing link between time-conditioned diffeomorphisms and Lagrangian solution maps for the continuity equation. The resulting parameterization allows to elegantly model time evolving densities with a single time-conditioned Normalizing Flow. Further- more, can calculate the velocity without inverting the conditional bijections, allowing the use of expressive bijective layers. We showed that LFlows can be applied to settings where we have sparse data on both density and velocity, and to settings where we have no data on velocity, but instead enforce additional equations.
In terms of density prediction L.llws outperform competing models on both synthetic and real experiments. Different to methods like PINNs, which weakly enforce the continuity equation, LFlows always provide physically consistent predictions. In addition, LFlows avoid scaling lim- itations of DFNNs (peak memory usage) and neural adjoint based methods (training time). For downstream tasks, LFlows directly provide Lagrangian maps without the need for additional numer- ical solvers. In dynamical optimal transport settings, LFlows directly provide the analytic expression of the transport map. When modeling bird migration, the Lagrangian maps provide access to trajec- tories, which could be compared migration paths obtained from different data modalities.

---

7 REPRODUCIBILITY STATEMENT
We provide general information on our implementation of LFlows, DFNNs, SLDA and PINNs in Section 4 and in the Appendix A.3] For each experiment we provide further details in the Appendix A.4 to A.6, where we state the layers, settings, and computational ressources (GPU) used. We also describe how we generated the introduced synthetic dataset. For the real-world bird-migration dataset we provided references for accessing the data, as well as information on the preprocessing in the Appendix [A.6] The supplementary material contains code for all experiments and models. The readme... mdd file lists the random seeds used for experiments and plots. The selected hyperparam- eters for each experiment are provided. The code further includes the generation of the synthetic datasets, as well as automated scripts for downloading and preprocessing the bird-migration data. A cleaned anaconda environment file for reproducing the python environment is provided.
We furthermore provide access to the used conditional bijective layers in separate python package called Flow Conductor This package includes conditional i-DenseNets with sinusoidal activations, well the conditional SVD layers.
REFERENCES
Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. Optuna: A next- next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pp. 2623-22222222 2631, 2019.
Luigi Ambrosio and Gianluca Crippa. Existence, uniqueness, stability and differentiability proper- ties of the flow associated to weakly differentiable vector fields. In Transport equations and multi- D hyperbolic conservation laws, volume 5 of Lect. Notes Unione Mat. Ital., pp. 3-57. Springer, 1.
Berlin, 2008. doi: 10. 101100977---554--77777
Fabricio Arend Torres, Marcello Massimo Negri, Monika Nagy-Huber, Maxim Samarin, and Volker Roth. Mesh-free eulerian physics-informed neural networks. arXiv preprint arXiv:2206.01543 2022.
Andrei Atanov, Alexandra Volokhova, Arsenii Ashukha, Ivan Sosnovik, and Dmitry Vetrov, Semi- conditional normalizing flows for semi-supervised learning. arXiv preprint urXiv:1905.000 2019.
Marin Biloš, Johanna Sommer, Syama Sundar Rangapuram, Tim Januschowski, and Stephan Günnemann. Neural flows: Efficient alternative to neural odes. Advances in neural informa- fion processing systems, 34:21325-21333 2021.
Nicolas Bonneel, Michiel Van De Panne, Sylvain Paris, and Wolfgang Heidrich. Displacement interpolation using lagrangian mass transport. In Proceedings of the 201 SIGGRAPH Asia con- ference, pp. 1-12, 2011.
Haim Brezis. Functional analysis, Sobolev spuces and partial differential equations. Universitext. Springer, New York, 2011. ISBN 978-0-387-70913
Dan G Cacuci. Sensitivity theory for nonlinear systems. nonlinear functional analysis approach. Journal of Mathematical Physics, 22(12):2794-288) 2802, 1981a
Dan G Cacuci. Sensitivity theory for nonlinear systems. ii extensions to additional classes of re- sponses. Journal of Mathematical Physics, 22(12):2033-2812, 1981b
Ricky T. Q. Chen.
torchdiffeq. 2018.
URL hatp:::/gittabb.cc//rrrt
Corchdiffedt
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Ricky TQ Chen, Brandon Amos, and Maximilian Nickel. Neural spatio-temporal point processes. In International Conference on Learning Representationns,2222 2020.
10

---

Phillip B Chilson, Phillip M Stepanian, and Jeffrey F Kelly. Radar aeroecology. Aeroecology, pp- 277-309, 2017.
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. arXiv preprint arXiv:2003.0465 2020.
Michail Diamantakis. The semi-lagrangian technique in atmospheric modelling: current status and futtrre challenges. In ECMWF Seminar in numerical methods for atmosphere and ocean mod- elling, pp. 183-200, 2013.
Michail Diamantakis and Linus Magnusson. Sensitivity of the eemwf model to semi-lagrangian departure point iterations. Monthly Weather Review, 144(9):3233-3250, 2016.
Adriaan M Dokter, Felix Liechti, Herbert Stark, Laurent Delobbe, Pierre Tabary, and Iwan Holle- man. Bird migration flight altitudes studied by network of operational weather radars. Journal of the Royal Society Interface, 8(54):30-43, 2011.
Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios. nflows: normalizing flows in PyTorch, November 2020. URL
22281/zenodo 4296287
Stefan Elfwing, Eiji Uchibe, and Kenji Doya. Siggoid-weighted linear units for neural network function approximation in reinforcement learning. Neural networks, 107:3-11, 2018.
Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanis- las Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, et al. Pot: Python optimal transport. The Journal of Machine Learning Research, 22(1):3571-3578 2021.
Will Grathwohl, Ricky T. Q. Chen, Jesse Bettencourt, and David Duvenaud. Ffjord: Free-form continuous dynamics for scalable reversible generative models. In International Conference on Learning Representations, 2019.
Samuel Greydanus. Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. Advances in neural information processing systems, 32, 2019.
David Ha, Andrew M. Dai, and Quoc V. Le. Hypemetworks. In International Conference on Learning Representations, 2017.
Philip Hartman. Ordinary differential equations, volume 38 of Classics Applied Mathemat- ics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002. ISBN 971119898888899222... Corrected reprint of the second (1982) edition
0-89871-510-5. doi: [Birkhäuser, Boston, MA, MR0658490 (83e:34002),, With foreword by Peter Bates.
Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, et al. The era5 global reanalysis. Quarterly Journal of the Royal Meteorological Society,
2020.
Ameya D Jagtap, Ehsan Kharazmi, and George Em Karniadakis. Conservative physics-informed neural networks on discrete domains for conservation laws: Applications forward and inverse problems. Computer Methods in Applied Mechanics and Engineering, 365: 13028, 2020.
George Em Karniadakis, Ioannis Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning, Nature Reviews Physics, 3(6))422-440, 2021.
Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible Ixl con- volutions. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Gamett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. URL https: di htttps://prrceeedagss..eoec www.:eewwee.....................
files/ pdf.
Ivan Kobyzev, Simon JD Prince, and Mareus A Brubaker. Normalizing flows: An introduction and review of current methods. IEEE transactions on pattern analysis and machine intelligence, 43 (11):3964-397,,2) 2020.

---

Lingxiao Li, Samuel Hurault, and Justin Solomon. Self-consistent velocity matching of probability flows. arXiv preprint arXiv:2301.13737, 2023.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08892 2020.
Huadong Liao and Jiawei He. Jacobian determinant of normalizing flows, 2021.
Fiona Lippert, Bart Kranstauber, Patrick D Forré, and E Emiel van Loon. Learning to predict spatiotemporal movement dynamics from weather radar networks. Methods in Ecology and Evo- lution, (3(1))2281 1-2826, 2022a.
Fiona Lippert, Bart Kranstauber, E. Emiel Loon, and Patrick Forré. Physics-informed inference of acrial animal movements from weather radar data. In NeurIPS 2022 Al for Science: Progress and Promises, 2022b.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for iden- tifying differential equations based on the universal approximation theorem of operators. arXiv preprint urXiv:1910.031990 2019.
Ashok Makkuva, Amirhossein Taghvaei, Sewoong Oh, and Jason Lee. Optimal transport mapping via input convex neural networks. In International Conference Machine Learning, pp. 6672- 6681. PMLR, 2020.
Raphaël Nussbaumer, Lionel Benoit, Grégoire Mariethoz, Felix Liechti, Silke Bauer, and Baptiste Schmid. A geostatistical approach estimate high resolution nocturnal bird migration densities from weather radar network. Remote Sensing, 11(19):223,, 2019.
Raphaël Nussbaumer, Silke Bauer, Lionel Benoit, Grégoire Mariethoz, Felix Liechti, and Baptiste Schmid. Quantifying year-round nocturnal bird migration with a fluid dynamics model. Journal of the Royal Society Interface, 18(179):2011014, 2021.
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lak- shminarayanan. Normalizing flows for probabilistic modeling and inference. The Journal of Machine Learning Research, 22(1):2617-2680, 2021.
Yura Perugachi-DDe -Daa,, Jakub Tomezak, and Sandjai Bhulai. Invertible densenets with concatenated lipswish. Advances in Neural Information Processing Systems, 1111762612227, 2021.
Lev Semenovich Pontryagin. Mathematical theory of optimal processes. CRC press, 1987.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A
deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Jack Richter-Powelle Yaron Lipman, and Ricky TQ Chen. Neural conservation laws: divergence- free perspective. In Advances in Neural Information Processing Systems, 2022.
Andre Robert. A semi-lagrangian and semi-im cemi-implicit numerical integration scheme for the primitive meteorological equations. Journal of the Meteorological Society of Japan. Ser. II, 60(1):319-355, 1982.
Christopher Shallue, Jaehoon Lee, Joseph Antognini, Jascha Sohl-Dickstein, Roy Frostig, and George E Dahl. Measuring the effects of data parallelism on neural network training. arXiv preprint FrXXiv:8111.03000 2018
Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit neural representations with periodic activation functions. Advances in Neural Information Processing Systems, 33:7462-7473. 2020.
Andrew Staniforth and Jean Côté. Semi-lagrangian integration schemes for atmospheric models-4 review. Monthly weather review, (19(9):206-2225) 1991.
12

---

A APPENDIX
A.1 THEORETICAL BACKGROUND
In this section we provide proofs and theoretical background for the method in Section While the underlying theory is well established, were unable find single source that concisely contained all required statements written in an accessible and directly citeable manner. Section 1.] contains the statements rely in Section 3 The proof's for these theorems provided in Sections A.1.2
to A.14
A.1.1 TIME DEPENDENT BIJETIONS AND THE CONTINUITY EQUATION
Theorem I provides velocity field given time-dependent bijections. Theorem 2 links the push- forward of the density to the solution of the continuity equation defined by an initial density and the velocity given by Theorem I
Theorem 1. Let 0 S < T and let Ω C IR be convex open set. Let X lo, T] x Ω 9 be a family of maps such that X₂ : Ω + a bijection for any t E [to. and (x,(x) ) x for any C ((to, T] x 12;21 with globally bounded derivatives.
Ω. Assume that X and X
(x,'(x)) is C In particular, satisfies the assumptions of the Cauchy-Lipsche Theorem and is the unique flow map of u starting at time to. Specifically, for any Ω the curve t + X2 X((x) is the unique solution to the Cauchy Problem
Then, the velocity field v(t, x)
aX.
{
0 (a) v (t, Xx())) E [t,
(16)
X to (x) x
Proof. See Appendix Section A.1.3 Theorem 2. Let Ω, T, lo, X be as in Theorem 7 Given an initial density Pho E LI))), we define
p(t, - Pto (X, (x)) JX, (x,'(()))
(17)
Then p(t, x) is distributional solution to the continuity equation in Eq. according to Definition 1 i.e. the following condition is satisfied for any test function & E C (coo x Q):
["' 10:0 - Vo)p dt = =-f ( x)) dr.
(18)
Moreover, if PL. E (R), then p E C (1to, x)) and i point-wise solution to the continuity equation Eq. 7 If we assume in addition that Pio (x) > 0 for any x E Ω. then the same holds for p(t, x) for any (t, x) E [to, T) Ω and p satisfies the log-density formula of the continuity equation
d log(p(t,X((x)))
- v(t,
(19)
dt
Proof. See Appendix Section A.1.2 and A.1.4
A.1.2 THE FLOW ASSOCIATED TO LIPSCHITZ VECTOR FIELD
We recall the setting the classical Cauchy-Lipschitz Theorem. For simplicity, Ω C Rd convex open set. Given 0 < to < T, let [to. T] x Ω Rd be a bounded vector field. We say that v is Lipschitz continuous in space uniformly in time if there exists a constant >> such that
[t x)- - vt, 4)1 < L x yl
Vt E [to, T Va, y E Ω.
(20)
Throughout this section, we consider maps X : [4o, T] x Ω + Ω with the following properties:
for any æ E Ω the map t + (x))) is 0 ([to: 71) with uniform bounds, namely there exists constant M > such that
[0.X((m)) S M VI E [lo, Vir Ω
(21)
13

---

for any t [to, T the map + (z:) is Lipschitz with uniform bounds, namely there exists a constant L > 0 such that
(,xx(()) (,)) - L| WE [oo T] Va, y Q
(22)
for any 1 E [l,, T] the map æ + (x) a bilipschitz transformation of Ω uniformly in time, namely X, is bijection of 52 and there exists constant C' > 0 such that
C x y < (x) X. S C y Væ, y E Ω VI E [40, T7.
(23)
Theorem in the case of IR and for the forward flow. We refer to Hartman (2002) for an extensive description the theory of ordinary differential equations, as well any book in basic differential calculus.
We state the Cauchy-
Theorem 3. Given T > 0 and to E [o, T), let v: [to, T] x IR + Rd be a bounded vector field that satisfies Eq. 20 for some constant L > 0. For any E I2 there exists a unique trajectory t ++ X((x) solving the Cauchy problem
{ 0 (a)) = v (X)(
E (to,
(24)
(1)(2) = x
in the integral sense. The map Xx [to: T) Rdd Rd is the flow of v starting time to and satisfies Eq. [7 Eq. [2] Eq. 23
Remark 1. Under the assumptions of Theorem point out that the maps X,, X, are Lipschitz continuous for any time. Thus, given time slice t E [fo, T", infer that X X, are differen- tiable almost everywhere in Idd, Hence, the Jacobian matrices J (x)) JX, (xx) are well defined for almost every x E Rdd and we have that
JX,(X, E (x)) = (x))- for almost every a E Bd.
Moreover, there exists a constant M > 0 such that
[JX,(x)) + JX, (a)) S M for almost every a E Rd
Here, |-| a given matrix norm (recall that all norms are equivalent in finite dimensional vector spaces). We point that the uniqueness part of Theorem 3 and the regularity properties of the flow, as well as the existence of the Jacobian matrix, extend to Lipschitz vector field defined on a general convex domain Qt as soon as we assume that the trajectories do not touch the boundary of Ω. In this case, we get that the flow map is a bilipschitz transformation of Q for any time slice.
A.1.3 VELOCITIES FROM I-PARAMETER GROUPS OF DIFFEOMORPHISM
Theorem M is a particular case of the following more general result. For simplicity, we consider convex domains.
Theorem 4. Let 0 < to T. let Ω C Rd be a convex open set and let X Noo
x S Ω
be family of maps on Ω such that X (x) = for any x E Q. Assume that X satisfies Eq. ax. St (X, (æ)) satisfies the assumptions of Theorem 3 and X is the unique flow map of starting at time to. Specifically, for any E Ω the curve t + X (x) ) is the unique solution to the Cauchy problem Eq. 16
27
Eq. 22 and Eg. 23 Then, the velocity field v(t,
the Cauchy-LLpsschn
The reader might note that Theorem iis the inverse of Theorem Indeed, given map satisfying all the properties of a flow map, it is natural to ask whether we can find velocity field v within the Cauchy-Lipschizz framework whose flow starting time tg X.
Proof of Theorem 7 Given (l, E [to, 7) x Ω we define
(t, m) lim XX+,(X, (x)) Xx(X, (z)) xx,
(25)
h-+0
h
Since the map - X. (x. (a)) ) is Ct ((to,, TJ) by assumption for any E Ω the velocity field
well defined. Moreover, by Eq. 25 the curve + X((x) solves the Cauchy problem Eq. 24 with given by Eq. 25 We study the regularity of to ensure that X is the unique flow associated to
v starting at time to. To begin, we remark that v is uniformly bounded by Eq. 21 Moreover, since
x, X, Lipschitz maps for any time slice t Ito, T by Eq. 22 and Eq. 23 we infer that
satisfies Eq. 20 Thus, a satisfies the assumptions of Theorem In particular, X Is the unique flow starting at time to of the vector field v.
14

---

A.1.4 THE CONTINUITY EQUATION
Let C Rd be an open set, let T and En E T). Given vector field v iu. T) x S IR', we shall consider the continuity equation Eq. on (to. x Q with initial condition Pto To deal with irregular vector fields and densities, we consider solutions to Eq. in the sense of distributions. We refer to Ambrosio & Crippa (2008) for some basic and advanced results on the theory of continuity equation.
Definition 1. Let Q C Rdd be an open set and 0 S to <T. Let o E 1. T) XQ; Rd) be a vector field and let Pto E (12). We say that p (Lto, L)(L)) is distributional solution Eq. if the following condition is satisfied for any test function @ C r[to,o x Q):
['' [ (0.00 Vo)p dt = - Pto (io(z) x) dr.
(26)
Remark 2. We point out that Definition [is well posed without differentiability assumptions on v. P. However, if v, p are C functions time and space and Pto is a continuous function, after integrating by parts the left hand side of Eq. 26 for any test function & E C ([fo, T) x 92) we obtain that
T S'
v. (up))odo d + J. p(tu x)o(lo do = - Pto (to()))((os x)) dx.
Since 0.p + V. (vp) is continuous function, by the so-called Fundamental Lemina of Calculus of Variations (see Brezis (2017). .) for instance) we infer that Dep+ + V (up) for any (t (to. T) Q and Pto (z) p(to: x) for any I E Q thus recovering the pointwise formulation of Eq. 7
Solving the Continuity Equation. The proof of the first part of Theorem 2 is corollary of the following general statement.
Ito. T] x Ω - Id be globally bounded velocity field that satisfies Eq. 20 Assume that the flow of u starting at time lo, denoted by X, is well defined in [ta. T x Q and that X, Ω - Q is a bilipschitz transformation of Q. Letting p(x; t) be defined by Eq. 17 then p is a distributional solution to the continuity equation Eq. according to Definition
Theorem 5.
Let Ω C. Rd be an open set and let 0 < to < T. Let U
We remark that, under the assumptions of Theorem 2 the flow map X is given and we build velocity field v that has as unique flow map. Hence, by construction, 12 invariant under X, for any t E (o, Tj. Thus, the assumptions of Theorem 5 are satisfied.
Proof of Theorem 5 By Theorem 3. the flow map X starting at time to associated to v well defined. Hence, defining p by Eq. 17 for any / E [L,, T] we have that p(t,-) is defined almost everywhere in Rdd We check that is a distributional solution to Eq. according to Definition To begin, we check that p ([to, T); L (IR5). Indeed, after the change of variables X, (x) y. using the Area formula with dy = |det Jx (dx),d, have that
[lott. x)) dx =
= 1 (X, (x)) J Jx, (x)) dr (1o)
Therefore, the total mass is preserved along the time evolution. Then, fix test function C (Ito, T) x 12) and, performing again the change of variables X, (x), we have that
/' T f 1000 + v. Volp dt
t0
= I.¹ / (0,((t,)) v(t, v(t, (x)) JX (x)) dr dt
1"' / [1h0(t. X.(v)) + v(l. (x)) Voli, X.(y)) (v) dy dt.
15

---

Recalling that X₂ is the flow map generated by the velocity field v starting at time to. the latter is equal to
('e (0rstt, Xi(v)) (,X,(y) Vo(t, ((y))) dy dt
= -f 1." Hlt, olt, (y)) dt dy,
dt
after using Fubini's Theorem and the chain rule for the derivatives. Thus, by the Fundamental Theorem of Calculus, we have that
I (ro(v)) S(" dt 011 (((,)(())) di dy
Pro (2) (T, XT - <00, X (v))) dy
o(to y)pro (y) dy,
12
since (TT, = 0 and X is the identity map.
Finally, we are able to conclude the proof of Theorem
Conclusion of the proof of Theorem 2 We discussed the fact that p is pointwise solution to the Indeed, by the explicit formula Eq. [7 is clear that is
continuity equation Eq. in Remark 2 O ([to. 7) x Q2) and that p nonnegative whenever Pto nonnegative. Thus, we check that Eq. 19 is satisfied. Indeed, by the chain rule we have that
Vp(t, Oeplt, x,()) Vp(t, X,(t, d x.()))
dt log(p(t, Xx(x))) = ptt, xx))))
dt
1iplt, Xi(x)) Vp(t, X, (x))) r(t, X,(x))
p(l, (x()))
Dep(1, X,(x)) v (o(1, {x(x)))(1 xx(x)))
p(t, x,()))
- (V oit, (x))) Xx
(z))),
since p satisfies the continuity equation at any point (t) E (to, T) Ω
A.2
CALCUUAATING FOR THE DENSITY AND VELOCITY
In this subsection we explicitly provide the steps to get to Eq. 12 and Eq. 13
A.2.1 CALCULATING THE DENSITY
We now report explicitly the steps to get to Eq. 12
Let f and g be diffeomorphisms on Rd, and A and B positive definite matrices. We denote with (g))) the Jacobian of f evaluated at the point g(z). In the following part we make use of identities that follow from the chain rule, the inverse function theorem, and properties of the deter- minant:
J(f o (z) = ((f)((())) (z))())
(27)
det (())(x) |det ((()(()
(28)
(eet(A B) |det(A)) /det (det())) (B)-
(29)
x(t,x)) Pro (X,"' (x)) det det J (x))
(30)
Prese ((De, X,, (zx) (4x, det ((x,))
(31)
16

---

Poaso (0r, o xx)))(z) = Prase (Br, P7v [bi](z)) = Poaso P (l((x))
(32)
det J(fto((X, (ax))
det ((Pto) ( ( -
(33)
((P,,)) (0,, $1o (1)(2)))
(34)
det (0,,,,))((()))
(35)
det J(X, (xx)
= |det |det.[( J[ to o (2))
(36)
det J (( (Bx(z)) det (x))
(37)
Combing the three terms we get:
p(t, Pro (X; (x)) det J (X, (x)]
(38)
Prace p (z) del J (0x))) |det ((0%))(0,))))
|det. J (x))
(39)
Prass (x)) det [()))(=))
(40)
A.2.2 CALCULATING THE VELOCITY WITHOUT INVERTING THE FLOW.
We now report explicitly the steps to get to Eq. 13 Firstly, we show that the velocity can be expressed in terms of the flow bijection &, and its inverse $7
ox,
d(t,x)
X (x,((x))
(41)
0 1
o Dto)
(10 (x)(x))
(42)
01
a(
(Di Dio (x)
(43)
Øt
0
(z((z))
(44)
In second step, the velocity can be written without the need for explicit inversion of the bijective layer tt, allowing for efficient computation in practice.
Let O and F," be the maps from 2; to z and vice versa respectively.
L = (z) = x2
(45)
and z does not depend on time, i.e. de 2 O. We can now explicitly compute the total derivative and find formulation for the velocity that requires inverting a Jacobian instead of inverting the map +,.
Clearly, P. (0,1
d ( p (0, (0,) dt
= dt 0
(46)
(oe)
(.(m))
0b,
Od; (Φ. (2))
04,-"
(-1(2)) (2)) + gz:
=
(47)
0t
Øt
0Φ-1
(0)()))) - [ )
(II)
(48)
0t
Øt
17

---

A.3 IMPLEMENTATION
This section provides high-level descriptions of the implementations of the models used in all exper- iments. For details on experiment-specific implementations, we refer to the corresponding sections
t.4 A.6
LFlows. Our code for LFlows is based on the nffows library for bijective neural networks (Durkan al. 2020). We extended nflows by providing range of additional (conditional) transformations, such as invertible densenets, addition to the modifications required for the LFlows. The code provided as supplementary material.
Embedding Network fe
( fa(t))
(1)
ou
m1(=** sa(t))
0cs(2 fo(t))
4 (x)
Figure A.6: General Architecture used for the conditional bijective layers based on embedding networks. Dotted lines indicate the inverse direction).
SLDA. We represent the density Ptg the departure time with bilinear interpolation of learnable equi-distant mesh in 2D or 3D. We parameterize the velocity with a neural network with smooth activation functions. In 2D settings, we evaluate the divergence Eq. exactly via autograd, whereas we stochasticly estimate it in 3D during training (Grathwohl et al., 2019). To solve multiple leverage the rescaling trick discussed in Appendix of Chen et al. (2020). In all settings we use a Dormand- Prince solver of order 5 with absolute and relative tolerance of 1e-5. The adjoint is computed with the torchdiffeq library (Chen 2018).
ODEs with differing
times, in parallel
DFNNs We rely on the vector-based parameterization of DFNNs given in Section 7.1 of Richter- Powell al. (2022), which ensures non -egnaivte densities.
We implement standard PINNs without additional resampling schemes or loss terms, and use either ReLu (Section 5.1) sinusoidal activations (Sitzmann al., (020) (Section (.3). The collocation points are resampled each training iteration with quasi-random Sobol sequences.
PINNs.
A.4 EXPERIMENT: SIMULATED FLUID FLOW
Data Generation. We generate the data by defining a Lagrangian flow map for an initial unnor- malized Gaussian mixture density. The initial density the simulated problem is based on a mixture of 4 independent Gaussians arranged around the origin with varying radii and a standard deviation of 0. We defined the density be restricted to [-4, 4] with (pv)lsa O:
po(x)= Ein W 0.1Z) if E (-4,4) else
(49)
The changes density on the zy axes are simulated by directly parameterizing the Lagrangian solution map X((xo) (-4,4) (-4,4) for t [0, 1.2]:
(0.5t + 1)
Aroe (t) Ascate (4 atanh(0.25an) + shhit(t))
Xe(xo) = 4- taanh
(50)
4
18

---

with
ccos(2
(2T()
Aro = sin()
(51)
cos (2T()
(scale(t) = + 0.14
(52)
1 0.1e
[ 0.6 0.6
shift(t) = sit(
(53)
In 3D, the z component of X, is always an identity map, limiting the dynamics to the ay axis.
The density and velocity are then given by:
p(t, x) = Po X, (x))
|det det JX,
(54)
(t, 0x, (x, (2:)
(55)
Öt
where all involved derivatives are computed using automatic differentiation. Observations are avail- able at 21 equidistant timesteps in the range '0, 1]. The test data set covers the time range [0, 1.2]. To simulate noisy mcasurements, additional Gaussian noise added to the observed velocities and log densities during training.
A.4.1 TRAINING AND ARCHITECTURE DETAILS
Hyperparameter Optimization. We optimize each model based on the explained variance (R²²) of the density on validation data. Firstly, we manually performed a general architecture selection. Subsequently, we tuned parameters such as the number of layers, units, learning rate, loss- and regularization weights with the black-box optimization framework Optuna Akiba et al. 2019). using the default Tree-structured Parzen Estimator as sampler. We trained the I.Flows and PINNs on a minibatch size of 16384 and the SLDA on a minibatch size 4096. As the 2nd order derivatives of the DFNNs require lot of GPU memory, DFNNs were limited to a minibatch size of 2048. For the optimized hyperparameters of each model we refer to the provided code in the supplementary material.
LFlows. For the initial layer we first rescale the domain and then use a atanh bijection for restrict- ing the domain to [-4, 4]d For the remaining layers we used blocks consisting of invertible Dense Nets (i-DenseNet) (Perugachi-Diaz et al., 2021) followed SVD layer conditioned on the embedding. For the i-DenseNet we rely on sinusoidal activations g(x) sin(15 z)/15. Each i-DenseNet has a depth of 3 and before each block we use an Activation Normalization layer. We enforce the Lipschitz constant of the i-DenseNet be 0.97.
DFNNs. The Divergence-Free Neural Networks do not directly provide access to the velocity v, but only the flux F (pv). Hence, calculating the velocity e F/p low-density regions leads to numerical issues. To avoid this, we train DFNNs directly on the flux instead of the velocity. Furthermore, we require non-negative densities, so we use the parameterization with subharmonic functions discussed in Section 7.1 of Riihhtrr-Poell et al. (022)). As the predicted densities can still zero, we train the DFNNS on the MSE the densities (instead of log densities).
To facilitate training of the PINN, we sinusoidal activation functions presented by Sitzzaan et al. 2020). We use a frequency multiplier of wp = 12 in the first layer. Collocation points are sampled within the full domain, uniformly distributed in [0, 1.2) x Ω. Instead of purely random sampler, we rely on quasi-random lowddiscrpy samples obtained via Sobol Sequences.
PINNs.
In each minibatch, 216 collocation points are sampled. The minimized PDE loss is L(I, [[0,p(1, x) + V (p(t, x (o, 22 averaged over the collocation points.
SLDA Relative and absolute tolerances of the solver during hyperparameter optimization were 10-3 and for the final run with the tuned hyperparameters 10-5 Lower tolerances during hyper- parameter search were not feasible, as they significantly increased the runtime for the problem at
https://optuna.oogg
19

---

hand. For stable dynamics, the hypernetworks provided by the code of Grathwohl et al. (2019) were necessary. These hyper networks are conditioned on time and provide a network taking the space coordinates imput, i.. ve(s) (2) with è being the hyper network. Using instead fully connected layers (that still fulfill the smoothness requirements) led in our experience to difficult dynamics for the adaptive ODE solvers.
Boundary Conditions. For PINNs, SLDA, and DENNS the boundary condition p(zr)n(=)) is enforced via an additional penalty on points sampled at the boundary, For LFlows, the boundary was enforced via a bijection $2 852.
Numerical Evaluation of Physical Consistency We compare the predicted density p(t, x) with the numerical solution of the initial value problem uniquely defined by p(to, æ), 0(t, x) and Eq. 3 We obtain the numerical solution with an 8th order Dormand- Prince ODE solver with an absolute and relative tolerance of le-5. We set to O. For the locations where we evaluate the sMAPE, we
consider 10 equidistant timesteps in [0.1, 1.2]. For each timestep, we first randomly sample locations (4, 4)d, and then we randomly subsample 2500 location with groundtruth densities larger than threshold of 0.1 in 2D or 0.01 in 3D. With this procedure we avoid regions with near-zero density.
Total Mass Regularization. We penalize the total of the system for all methods except PINNs. For the LFlows and SLDA, we penalize the learnable normalization constant. For the DFNNs, no equivalent to the normalization constant is available. We instead introduce the penalty at points sampled on the domain (0. 1.2] x Ω. For PINNs additional regularization is not necessary. This is due the side effect of the PDE loss being numerically small for small densities, leading to an automatic built-in penalty for large total mass.
Computational Resources. Each individual experiment for the synthetic data was individ- ual NVIDIA TITAN X GPUs (I2GB VRAM), using 20 CPU cores and 20GB RAM. To speed up hyperparameter tuning, up experiments were in parallel using a SLURM-based compute cluster.
A.5 EXPERIMENT: DYNAMICAL OPTIMAL TRANSPORT
In Figure A.5 and Figure A.S interpolated densities are shown for the (Caussians++Cirdd dataset and the Sewweele-&Gauseeeee dataset.
LFlow
DFNN
Discrete
Figure A.7: Approximations of the 2D optimal transport map for the 8Gaussians+CCrrlls dataset with LFlows, DFNNs, and a discrete reference.
A.5.1 EVALUATING THE INTEGRAL.
For evaluating the integral required for the objective in Eq. 14 we reuse the code and method provided by RRihherr- Powell et al. (2022). The integral is estimated via importance sampling, with the sampling distribution q(1, m) (()) Samples in time drawn uniformly with q())
20

---

5
LFloo
DFNN
Discrete
Figure A.8: Approximations of the 2D optimal transport map for the Pinwheel=+&Ganssians dataset with LFlows, DFNNs, and a discrete reference.
utto f1). Samples in space are drawn from a uniform mixture
1
(0)) = pP (x) + pPh (ax) +
3 (u(((r)
(56)
with Un being the uniform distribution on the domain.
A.5.2 DISCRETE ESTIMATE OF W2
For estimating the W2 distance with discrete reference method we rely on 50000 samples. In our setting, more samples lead to convergence issues of the numerical solver.
As the discrete W2 estimate depends the number of samples n, explore increasing number of samples and repeated runs. For each fixed n we repeat the estimate 10 times with different seeds and visualize the resulting box plots in A.5.2 We note that the differences with increasing sample size are comparatively small relative to the differences between the continuous methods in Section
Circles
Pinwheol Gaussians
Gaussians
0.12
0.27
0.170
0.11
0.26
0.165
0.10
0.160
0.25
we
00009
0.155
0.24
0.07
0.23
0.145
0.06
0.22
0.140
30000
50000
30000 #samples
10000
50000
#samples
#samples
Figure A.9: Change of the discrete W estimate based increasing number of samples. Each boxplot is based on 10 repeated runs with different seeds.
A.5.3 IMPLEMENTATION
LFlows. For LFlows we use 10 blocks of i-DenseNets, each preceded by Activation Normal- ization layer. Each i-DenseNet has a depth of with sinusoidal activations g(x) = sin(15 * x1/15. We enforce the Lipschitz constant of the i-DenseNet to be 0.97. For the embedding of the time condition we use a residual neural network with hidden layer, swish activations, and 128 hidden features. The dimension of the time embedding, i.e. the output of the embedding network. is of dimension 10. The model was trained with the ADAM optimizer for 5000 iterations with a learning rate of 2e-3 and 2048 data points per iteration.
21

---

DFNNS. For DFNNs we directly used the experiment code provided by Richter- Powell et al. 2022 That is, the DFNNS are based on the parameterization with subharmonic functions dis- cussed in Section 7.1 of Rochie et al. (2022). The 128 mixtures of the subharmonic function are parameterized with a layer neural network with 96 hidden features and swish activations. A fixed À weight of 50 is set for the OTT-peaalty for all datasets. The DFNN is trained with the ADAM optimizer and a learning rate of le-3 over 10 000 iterations with 256 data points per iteration.
Computational Resources. The experiment was run on individual NVIDIA TITAN X GPUs (12GB VRAM), using 20 CPU cores and 20GB RAM.
A.6 EXPERIMENT: BIRD MIGRATION
About the data. The data provided by Nussbaumer et al. (2021) is originally based weather radar measurements made available by the European Operational Program for Exchange of Weather Radar Information (EUMETNET/OPERA, The vertical profiles, i.e. density and velocity estimates at different altitude levels, were provided by the European Network for the Radar surveillance of Animal Movement (ENRAM), based on vol2bird an algorithm for preprocessing raw radar scans. The raw data consists of volume scans from Doppler radars, measuring reflectivity and radial veloc- ity of the surroundings. By filtering out biological and environmental scatters, it possible to retain scans that mostly contain bird movements. Based on the reflectivity and radial velocity, the average bird density and velocity within a 15km radius is estimated for multiple altitude bins. For details about this process, we refer to Dokter et al. (201T). The final density and velocity measurements we use are openly available
Preprocessing. The positions of the radars are given in the WGS84 coordinate reference system. We project it to the Cartesian reference system ETRS89-extended (EPSG:3035), effectively pro- jecting longitude/latitude to x and y coordinates given in meters. As an additional preprocessing step, we excluded velocities that were measured together with near-zero density. To generate the data set used in our experiment, we directly concatenated multiple nights and remove the daytime during which measurements are available. In the supplementary material we provide code for downloading and preprocessing the data.
Numerical Evaluation of Physical Consistency We compare the predicted density p(t, a)) with the numerical solution of the initial value problem uniquely defined by p(to, x), é(t, x) and Eq. 3 The time of the earliest available datapoint within the three nights serves as to- We evaluated the error at spatial equidistant grid of sidelength 50 in the zy dimension. The coordinate is randomly sampled for each of the grid entries. We evaluate the consistency loss at these zy² coordinates 10 time steps between the start and end of the three selected nights, and average the sMAPE over all zyzz coordinates. For all models we use a 8th order Dormand-Prince ODE solver with an absolute and relative tolerance of 1e-5, which is one order higher than the ODE solver used for SLDA. The MLP has no z-component of the velocity, as no measurements of it are in the data. We thus set this component of the MLP zero for computing the numerical ODE solution. The other models indirectly learn z-component by enforcing the CE in 3D.
A.6.1 TRAINING AND ARCHITECTURE DETAILS
LFlows. As first layer of the LFlow we use alanh bijection that constrains Ω to rectangular volume that is multiple times larger than the spatial extent of the radar positions. The following layers consist of 10 blocks of invertible Dense Nets (i-DenseNet) (Perugachi-Diaz where we instead use sinusoidal activations g(x) = sin(10 * z Each block has depth of 5 and an Activation Normalization layer. We enforce the Lipschitz constant of
al., (021I)
before each block the N-DeneeNet to be 0.97. We use a 2 layer residual network with a width of 128 for embedding the time before passing it to the i-DenseNet.
We initialized the log of the normalization constant In(c) with 18.2 and set the weight of the total mass penalty le-3. We trained for 50 epochs with minibatch size of 16384 using the ADAM
httpss//anndd.rrg/reer////33333
22

---

optimizer with learning rate of 1e-2, weight decay of 2e-3 and cosine annealing learning rate schedule.
SLDA. For the SLDA represent the initial density with grid of size 20 20 20. Values in between mesh points are evaluated with bilinear interpolation. We parameterize the velocity network with hypemetwork as provided by the code of Grathwohl simulated fluid flow experiment. The network consists of 5 layers with 512 hidden features and swish activations. We trained for 300 epochs with a minibatch size of 16384 using the ADAM optimizer with a learning rate of 1e-3, weight decay of 5e-3 and a cosine annealing learning rate schedule.
al. (2019), similar to the
the parameterization for non-negative densities based subhar- monic functions (Section 7.1 Richter-Powell et al. (2022). The model consists of 4 layers with 256 hidden features and Swish activations, which parameterize 64 mixture components for the sub- harmonic function. We used a minibatch size of 2048 due to memory constraints. The model was trained for 100 epochs with cosine annealing learning rate schedule.
DFNN. For the DFNN
MLP. We trained 10 layer residual neural network with ReLu activations, intermediate batch normalization and 256 hidden features. The model was trained for 100 epochs with minibatch size of 16384 using the ADAM optimizer with a learning rate of 1e-3, a weight decay of 1e-3 and cosine annealing learning rate schedule.
the same architecture as for the MLP, but with additional PDE loss. We use 100 000 collocation points sampled from a quasi-random Sobol sequence and weight the PDE loss with 7c-4. The PINN was trained for 300 epochs with minibatch size of 16384 using the ADAM optimizer with a learning rate of le-2, weight decay of le-2, and a cosine annealing learning rate schedule
PINN. For the PINN
Computational Ressources. The experiment was run A100 GPU (40GB VRAM), using 20 CPUs and 30GB RAM. Repeated experiments were parallelized on a SLURM based compute cluster.
23

---

A.6.2 PREDICTIONS OF ALL MODELS
2018 06 00:00:00
2018 04 07 00:00.00
2018-04 00:00:00
LFlow
2018 06 00:00.00
2018-85-07 00:00000
SLDA
2018 00:00.00
2018-04-07 00:02000
2018 04 OB 00:00.00
DFNN
PINN
2010-04-06 00:00:00
201D-0N-77 8O108:30
2015-08-28 DC:00-00
MLP
Figure A.10: Predicted density and normalized velocity integrated over z-axis for all considered models.
A.6.3 TOTAL MASS PENALTY
To motivate the total mass penalty we showcase its effect on bird migration predictions using LFlows with varying penalty weights. The different predictions shown in Figure A.11
24

---

Without any penalty, the network freely explain the training data with large amounts of total mass. This leads to significant and never observed densities outside the observed region. However, with active total mass penalty effectively being prior the network is instead encouraged to decrease the total mass. The total mass penalty weight allows to explore this aspect of the ill- posedness of the problem.
Without Penalty 2013 05 02.00
Penalty weight 1e-2
Penalty weight 6e-2
to
Figure A.11: Predictions of the LFlow trained without (left) and with (right) total mass penalty. While predictions at observed radar stations barely change, the total mass outside of the observed region is significantly reduced.
The variation in predictions can be further summarized in a single visualization. Similar to before, train different models for varying total mass penalty weights. We calculate the relative standard deviation the predicted flux (i.e. the product density and velocity) each spatial location for fixed time. Areas with higher relative standard deviations correspond areas with largely varying explanations. Figure A212 shows the resulting map for a single time frame for the bird migration setting. The areas which are never observed and are completely dependent highest variation. Areas closer to the train radar stations (light green points) have a lower variation.
prior show the
0.5
Figure 121 Relative standard deviation of predicted flux with varying mass penalties.
25


Reviews:
[{"review_id": "ICFJDgfNOh", "reviewer": "Reviewer_1yfJ", "summary": "The paper introduces a neural network-based model designed to adhere to the continuity equation, even when dealing with scenarios where precise boundary and initial conditions are unknown. The core technical foundation relies on the classical theory of Lagrangian flows, providing a toolbox that ensures the evolution of density and velocities complies with the continuity equation.", "questions": "Could the authors comment on the challenges they anticipate when scaling this approach to larger problems, such as modeling large turbulent flows in aerodynamics? How might the complexity of LFlows impact scalability in such scenarios?\n\nHow does the computational complexity of LFlows compare to existing methods for solving hydrodynamic flow problems, especially in high-dimensional scenarios? Are there specific trade-offs or advantages in terms of computational efficiency?\n\nDoes the methodology proposed in this paper have the potential to be generalized to tackle other partial differential equations (PDEs) beyond the continuity equation? If so, what challenges or modifications might be necessary when applying it to different PDEs?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper significantly contributes to the field of neural PDE methods by focusing on enforcing constraints. It marries neural networks with conditioning normalization flows, resulting in LFlows that satisfy the continuity equation by design.\n\nThe method is validated across various 2D and 3D problems, including a real-world bird migration dataset, which could potentially serve as a standard benchmark dataset for similar problems. Unlike previous methods that often require known PDEs, initial, and boundary conditions, this approach tackles data assimilation-style problems where many system parameters are unknown. This kind of problem could be particularly useful for biomedical imaging applications.\n\nThe paper rigorously compares its approach to Richter-Powell's divergence-free neural networks, showcasing the computational efficiency of its method without the need for higher-order autodiff. Additionally, it compares favorably against other standard baselines, such as PINNs, especially in spatially sparse settings.", "weaknesses": "The proposed method inherits limitations associated with normalizing flows and bijective layers, such as difficulties in handling discontinuities.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8gx1x6FYaf", "reviewer": "Reviewer_fhDc", "summary": "This paper expresses the continuity equation into the time-dependent density transformation to model the fluid densities and velocities. By introducing differentiable and invertible maps, the proposed LFlows can ensure the continuity equation inherently. The authors provide complete proof to support their insights. Typical experiments are also included. Experimentally, LFlows performs best in both performance and efficiency.", "questions": "All the questions are listed above.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.\tThe paper starts from an interesting insight. The proposed method is reasonable and well-supported.\n\n2.\tThe authors conduct extensive experiments, covering both simulated and real-world dataset. \n\n3.\tThe propose LFlows performs well in both performance and efficiency.", "weaknesses": "1.\tAbout the presentation.\n\nFor me, the writings in Section 3.1 and 4 are confusing and hard-to-read.\n\nIn Section 3.1, the authors use several theorems in an informal way and sometimes not give the exact formalization of theorem. I prefer to place the definitions of theorems in the main text and leave the proof in appendix, such as Theorem 1,2 and the existence and uniqueness of Lagrangian flows for smooth vector fields.\n\nIn Section 4, I don’t think the implementation of LFlows are well described. The model architecture in Appendix A.3 should be placed in main text. Besides, the implementation for baselines should be deferred into the Appendix.\n\n2.\tHow to adopt LFlows to process the spatially sparse data? More details are expected.\n\n3.\tAbout the noise observations in the real-world dataset. Since LFlows enable the strict continuity equation, does it come across problems when the input data is noisy even wrong? Some discussions are expected.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "M1iXKsk5fI", "reviewer": "Reviewer_RQ15", "summary": "This paper proposes LFlow, a neural parameterization of both the density and the velocity field that adhere to the continuity equation. This is achieved with a link between time-conditioned diffeomorphisms and Lagrangian solution maps for the continuity equation. Two application settings were considered, where in one setting the user has sparse data on both density and velocity, and in another setting, the user has no data but knows the equations. Low-dimensional experiments are done to demonstrate the effectiveness of the proposed parameterization.", "questions": "1. At the bottom of the first page, it reads \"in setting (ii) we measure only the density of the fluid but we know additional equations defining the velocity.\" If density is already known for all time, then isn't it already enough to reconstruct the velocity without needing to know any additional equations, by simply solving the continuity equation?\n2. In the evaluation in Section 5.2, the gold standard of $W_2^2$ is computed using discrete estimates (and then averaged). However, the discrete estimate of W2 is known to be biased. How can one be sure that the red bars on the right of Figure 3 are accurate? Would the bar shift if more samples were used?\n3. In (15), how is the double integral computed in practice? It looks like it requires uniform samples from the domain $\\Omega$, which could be very ineffective if the flow only concentrates on a small portion of the domain.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "* The writing of the paper is clear.\n* Both settings are interesting and the proposed method has been shown to outperform the alternatives (PINN, DFNN, SLDA).", "weaknesses": "* To my knowledge, the idea of using the Lagrangian view of the continuity equation to build a neural parameterization that gives both velocity and density access is not new. \n  * In Section 3.1 and 3.2 of [1], the TIPF parameterization appears to be identical to LFlow proposed in the present paper, with the only difference being that in (13) of the present paper, $\\Phi_t$ appears instead of the inverse, which is the case of TIPF.\n  * This is a minor difference for TIPF, since in [1] the conditional normalizing flow used is RealNVP [2], which has ready access to the inverse flow map.\n  * Could the authors elaborate on the advantages of avoiding the inverse here?\n  * Could the authors also comment on the choice of the Lipschitz-constrained invertible densenets as the backbone architecture? What is the advantage/disadvantage compared to e.g. RealNVP [2]?\n* The experiments section is weak in my opinion, and weaker than that in [Richter-Powell 2022]. If I understand correctly, all experiments are in 2D. How would the parameterization fare in higher dimensions? Higher dimensional experiments are done in [Richter-Powell 2022]. In [1] it is found that TIPF does not work well in dimension >= 10. I suspect the same behavior occurs here, due to the reason that invertible pushforward networks might not be expressive (or might be hard to train).\n\n[1] Li, Lingxiao, Samuel Hurault, and Justin Solomon. \"Self-consistent velocity matching of probability flows.\" arXiv preprint arXiv:2301.13737 (2023).\n\n[2] Dinh, Laurent, Jascha Sohl-Dickstein, and Samy Bengio. \"Density estimation using real nvp.\" arXiv preprint arXiv:1605.08803 (2016).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)