Marouane Tliba

Aladine Chetouani

Giuseppe Valenzise

Fréderic Dufaux

PCQA-GRAPHPOINT: EFFICIENT DEEP-BASED GRAPH METRIC FOR POINT CLOUD QUALITY ASSESSMENT

Keywords: 3D Point Clouds, Image Quality Assessment, Graph Neural Network, Deep Learning

Following the advent of immersive technologies and the increasing interest in representing interactive geometrical format, 3D Point Clouds (PC) have emerged as a promising solution and effective means to display 3D visual information. In addition to other challenges in immersive applications, objective and subjective quality assessments of compressed 3D content remain open problems and an area of research interest. Yet most of the efforts in the research area ignore the local geometrical structures between points representation. In this paper, we overcome this limitation by introducing a novel and efficient objective metric for Point Clouds Quality Assessment, by learning local intrinsic dependencies using Graph Neural Network (GNN). To evaluate the performance of our method, two well-known datasets have been used. The results demonstrate the effectiveness and reliability of our solution compared to state-of-the-art metrics.

INTRODUCTION

Recently, with the increasing demand for immersive applications and the widespread adoption of 3D acquisition devices, 3D Point Cloud (PC) data becomes more and more popular due to its capability to display the real world in modern immersive communications systems. Point clouds are becoming more widely used in a variety of fields, such as robotics, 3D gaming and telepresence. PCs consist of unordered geometrical coordinates (x, y, z) as well as potentially associated attributes such as color, curvatures, reflectance and normal vectors. Consequently, PCs can accurately represent a 3D model or an entire scene by incorporating thousands or even millions of points, requiring efficient compression schemes. To optimize the coding of point clouds and obtain the best visual quality for a given bitrate, accurate perceptual quality metrics are needed. Hence, point cloud quality assessment (PCQA) has received a great deal of attention from the research community in recent years.

A number of point cloud quality metrics have been developed in the past few years. They can be classified into three main groups: Point-based, Feature-based and Projection-based metrics. Point-based metrics such as Point-to-Point (Po2Po) [START_REF] Tulvan | Evaluation criteria for pcc (point cloud compression)[END_REF], Point-to-Plane (Po2Pl) [START_REF] Tian | Geometric distortion metrics for point cloud compression[END_REF], Plane-to-Plane (Pl2Pl) [START_REF] Alexiou | Point cloud quality assessment metric based on angular similarity[END_REF] and Pointto-Mesh (Po2Mesh) [START_REF] Rochinni | Metro: measuring errors on simplified surfaces[END_REF], predict the quality through point-wise geometric and/or features distance between the reference PC and its distorted version. Po2Po measures the relative distance between point pairs to estimate the final quality, Po2Pl extends Po2Po by projecting the error vector along the local normal, and Pl2Pl quantify the quality through measuring the angular similarity between surfaces associated to the points from the reference and degraded contents. Po2Mesh creates a polygonal mesh from the reference sample and then compute the distance between each distorted point and the corresponding surface. Currently, MPEG is adopting Po2Po MSE and Po2Pl MSE with the associated PSNR as the standard point cloud geometry quality metrics. Feature-based PC quality metrics extract the geometry with the associated attributes from point-wise level in a global or local way. Among those metrics, we can cite PC-MSDM [START_REF] Digne | Pc-msdm: A quality metric for 3d point clouds[END_REF] that extends the 2D SSIM metric [START_REF] Zhou Wang | Image quality assessment: from error visibility to structural similarity[END_REF] to PC by considering local curvature statistics, the Geotex [START_REF] Diniz | Towards a point cloud quality assessment model using local binary patterns[END_REF] metric that exploits the Local Binary Pattern (LBP) [START_REF] Pietikäinen | Two decades of local binary patterns: A survey[END_REF] descriptors, and PCQM [START_REF] Digne | Pcqm: A full-reference quality metric for colored 3d point clouds[END_REF] that combines the geometry and color features. In projection-based PC Quality Metrics, the 3D points or their associated features are projected into 2D regular grids, and later, 2D methods are applied on these views [START_REF] Bourbia | Blind projection-based 3d point cloud quality assessment method using a convolutional neural network[END_REF], or 2D grid of features [START_REF] Chetouani | Deep learning-based quality assessment of 3d point clouds without reference[END_REF] [START_REF] Chetouani | Convolutional neural network for 3d point cloud quality assessment with reference[END_REF]. The above point-based metrics or projection based metrics require a lot of pre-processing time, as well as the reference point clouds. Recently, new efficient deep based metrics such as PointNet-SSNR [START_REF] Tliba | Representation learning optimization for 3d point cloud quality assessment without reference[END_REF] and PointNet-DCCFR [START_REF] Tliba | Point cloud quality assessment using cross-correlation of deep features[END_REF] have been released. They exploit the intrinsic features of point cloud data, by operating on the native point cloud directly without applying any projection or transformation. A main limitation of these methods is that they mainly use point-based convolutions to extract features from points, disregarding the local relationships among points.

To overcome this limitation, we propose in this paper a novel efficient end-to-end deep graph-based metric for no reference point cloud quality assessment that considers the structure of the point cloud. Unlike previous methods, we take into account geometrical local connectivity between points, as well as hidden deep representation dependencies. Our method is designed to capture efficiently geometrical local point cloud structures through a dynamic graph construction at each network layer, and improves the aggregation of distant point dependencies in the la-Fig. 1. General pipeline of our proposed method tent space using cross/self-attention mechanisms. In this paper we demonstrate the effectiveness of our approach in estimating the visual quality of point clouds corrupted by compression artifacts, which is one of the main distortion scenarios in practice. The main contributions of this paper are summarized as follows:

1. We propose a novel efficient end-to-end deep graph based metric for PCQA, that operates on the whole point cloud without any projection or other transformation.

2. Our method is designed to capture efficiently geometrical local point cloud structures through the dynamic graph construction at each layer, and draws better aggregation of distant dependencies on the latent space using cross/selfattention mechanism.

PROPOSED METHOD

The main contribution of this work is to provide an efficient and robust blind metric for estimating the quality of 3D compressed PCs. Specifically, we propose a graph-based metric considering local connectivity information between adjacent point representations. The overall pipeline of our method consists of several steps: pre-processing, features extraction and aggregation using dynamic graph convolutional neural network (DGCNN), attention mechanisms and finally quality estimation. The general pipeline is depicted in Fig. 1.

Pre-processing

The pre-processing stage consists in dividing the PC into vertical slices of points. This partitioning has two goals: on one hand, to enable parallel processing of points, and on the other hand, to meet memory requirements for GPU processing. The number of partitions varies based on the size of the original PC.

In our tests, we use between 8 and 24 partitions. Each new point set (partition) is then divided into subsets in order to form a local graph on each of them. To do so, we select a random number of (distant) centroids, and apply the k-nearest neighbor clustering method to form a patch around each centroid. Each point subset is then fed into our graph model for graph construction and features extraction.

Proposed Graph Model

Our proposed model is mainly inspired from Pointnet++ [START_REF] Qi | Pointnet++: Deep hierarchical feature learning on point sets in a metric space[END_REF] and DGCNN [START_REF] Wang | Dynamic graph cnn for learning on point clouds[END_REF]. In particular, we leverage the information from local structures by building a graph based on neighboring points, and apply convolution using the graph edge information. Also, similar to DGCNN [START_REF] Wang | Dynamic graph cnn for learning on point clouds[END_REF], we update dynamically the graphs after each layer of the network, to better capture the dependencies at different feature levels. However, unlike DGCNN that is designed to extract semantic features for high-level computer vision tasks, our model employs graphs only to extract local features, and uses instead self/cross attention to process semantic information. More specifically, we use two-stream graph neural networks, one for the geometrical information, and one for color information. At a later stage, we aggregate the two streams using a cross-attention mechanism.

In more details, at each layer of the network, we construct a graph using k-nearest neighbors (k = 6) of a point. In the first layer, the graph is built using the geometric distances among points. At subsequent layers, the graph is constructed using distances between embeddings of points at that layer, thus effectively capturing local structures in the feature space. This results in a new adjacency matrix computed at each layer. We use the adjacency matrices of the geometry stream also for the color stream, in order to align the color information propagation with local geometrical structures.

Two-stream graph network

We consider geometry and color information as features, corresponding to two different processing streams. Depending on which stream we consider, input points x i bring different information. For the geometry stream, x i = (x i , y i , z i ) contains the 3-dimensional coordinates of the points, while for the color stream, points represent the RGB attribute information. Notice that it is possible to include additional features in other streams.

We construct two directed graphs by finding the k-nearest neighbors of each point G xyz = (V xyz , E), and G rgb = (V rgb , E), representing point cloud local structures along side the geometric and color features, where V = {1, . . . , k} and E ⊆ V × V are the vertices and edges, respectively. Each graph is also completed with self-loop.

The edge features at each network layer of the two streams are defined as an asymmetric function e ij = h Θ (x i , x i -x j ), where x j is the central vertex, {x j : (i, j) ∈ E} is the neighbor vertex, and h Θ : R 3 × R 3 → R F represents a nonlinear function with learnable parameters Θ. The adopted convolution operation, so-called Edge Conv, as depicted in Fig. 2, is defined by applying a max channel-wise aggregation of the features associated with all the edges resulted from each vertex. Therefore, the output at the i-th vertex of each stream is given as

x ′ i = max j:(i,j)∈E h Θ (x i , x i -x j ). (1)
In order to have faster convergence, after each Edge Conv layer we applied GraphNorm [START_REF] Cai | Graphnorm: A principled approach to accelerating graph neural network training[END_REF] operation, an adapted layer of InstanceNorm [START_REF] Ulyanov | Instance normalization: The missing ingredient for fast stylization[END_REF] for graph normalization with a learnable parameter α that controls how much of the channel-wise average to preserve in the shift operation:

x ′ = x ′ -α • E[x ′ ] Var[x ′ -α • E[x ′ ]] + ϵ • γ + β (2)
where γ and, β are the affine learnable parameters similar to other normalization methods.

After multiple consecutive (Edge Conv, GraphNorm) layers on each stream, a Nodes Aggregation is applied on each graph independently. Here a Max and Mean pooling operation are applied on node channel dimension. Afterwards, the result of the two pooling is concatenated to one vector, so each graph G, transformed into sequence of 2F -dimension vectors, denoted as Hidden.

Two Streams Fusion

The representations produced by the two streams, Hidden rgb and Hidden xyz , are finally fused using a cross-attention layer [START_REF] Chen | Crossvit: Cross-attention multi-scale vision transformer for image classification[END_REF]. More precisely, the fusion is carried out by measuring the similarity between the queries q (a transformed vector of Hidden rgb ) and the keys k (a transformed vector of Hidden xyz ) as well as using it to adjust the values vector v (a transformed vector of Hidden xyz ). Mathematically, the cross-attention can be expressed as follows:

q = Hidden rgb W q , k = Hidden xyz W k , v = Hidden xyz W v (3) A = sof tmax(qk T ) C h (4)
CrossAttention(Hidden xyz , Hidden rgb ) = Av [START_REF] Digne | Pc-msdm: A quality metric for 3d point clouds[END_REF] where W q , W k , W v ∈ R C(C/h) are learnable parameters, C and h are the embedding dimension and number of heads. We refer the reader to [START_REF] Vaswani | Attention is all you need[END_REF][START_REF] Chen | Crossvit: Cross-attention multi-scale vision transformer for image classification[END_REF] for more details about the computation of cross-attention.

Feature Aggregation and Quality Estimation

To better aggregate the fused representation and draw the relationship between distant embedded features on the output of the fusion stage, we employed also a self-attention layer [START_REF] Vaswani | Attention is all you need[END_REF]. Unlike cross attention [START_REF] Chen | Crossvit: Cross-attention multi-scale vision transformer for image classification[END_REF], here the set of k, q, and v are belonging to the same distribution of features [START_REF] Tliba | Satsal: A multi-level self-attention based architecture for visual saliency prediction[END_REF]. Following to that, similar to transformer block [START_REF] Vaswani | Attention is all you need[END_REF], we add skip connections with layer normalization, followed by a non-lineaar function.

Finally, we use a Graphs Aggregation layer, which applies a max poling operation over the vector output by the transformer block (TB), in order to have one vector that represents each PC partition. The later is passed to a very shallow multi-layer perceptron (i.e. regressor) in order to estimate the quality score of each PC partition. The final quality score of the whole PCs is obtained by applying the mean on the sequence of partitions scores.

Training Protocol and Implementation details

We train our model in an end-to-end manner, using the mean square error (MSE) as a loss function. The network is trained to create a mapping function between the input PCs and the mean opinion score (MOS) quality. The loss function is defined as follows:

L = MSE mean i Out i , Y . (6)
where Out i refers to each partition predicted score, and Y refers to the MOS. We implemented our model in PyTorch. The model was trained using Adam optimizer [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] with initial learning rate of 0.0001 and a batch size set to one.

EXPERIMENTS

To validate the effectiveness of our model, we evaluate it on two publicly available benchmarks with subjective scores, which adopt different emergent compression schemes with multiple encoding levels: ICIP20 [START_REF] Perry | Quality evaluation of static point clouds encoded using mpeg codecs[END_REF] and PointXR [START_REF] Alexiou | Pointxr: A toolbox for visualization and subjective evaluation of point clouds in virtual reality[END_REF]. ICIP20 is composed of 6 reference PCs from which 90 degraded versions were derived through three types of compression: V-PCC, G-PCC with triangle soup coding and G-PCC with octree coding [START_REF] Schwarz | Emerging mpeg standards for point cloud compression[END_REF][START_REF] De Queiroz | Compression of 3d point clouds using a region-adaptive hierarchical transform[END_REF][START_REF] Pavez | Dynamic polygon clouds: representa-tion and compression for vr/ar[END_REF]. Each reference PC was compressed using five different PointXR is composed of 5 PCs from which 45 degraded versions were derived through G-PCC with octree coding for geometry compression and, Lifting and RAHT for color compression.

To train and test our model on ICIP20, we split the dataset according to the reference point cloud by applying a 6-fold cross validation protocol where 6 corresponds to the number of reference point cloud samples, i.e., at each iteration 5 reference point cloud samples and their compressed versions are used for training, and one reference point cloud sample and its compressed versions are used for testing. Then, in order to evaluate the generalization capacity of our method to predict the quality on unknown PCs, a cross-dataset evaluation is employed by training our model on ICIP20 and testing it on PointXR datasets. Pearson Correlation Coefficient (PCC) and Spearman Rank-Order Correlation Coefficient (SROCC) are considered as criteria for measuring the quality prediction ability of our method, see e.g. ITU-T Rec. P1401. Notice that the correlations are computed for each fold and the mean over all folds is reported. The results obtained are shown in Tables 1 and2. We compared the performance of our method to state-of-the-art methods by applying strictly the same protocol (i.e. splitting and reporting the mean correlations).

In Table 1, we show the performance reached by our method on ICIP20 dataset for the two streams independently, denoted here as Color stream and Geometry stream, as well as after considering the two streams (i.e. fusion). The results are also compared to state-of-the art methods. As we can see, the two streams are well correlated with the subjective ground truth. The geometry stream predicts the quality much better than the color one, especially according to PLCC. After combining the two streams, the performance improves with a clear gap for both PLCC and SROCC. Also the ablation shows the importance of employing transformer block (TB) before the graph aggregation. Compared to the existing methods, the proposed one outperforms most of them with the best SROCC and the second best PLCC. It is worth noting that our method surpasses most of the fullreference methods. Notice also that all the listed methods in the table are full-reference except PointNet-SSNR [START_REF] Qi | Pointnet++: Deep hierarchical feature learning on point sets in a metric space[END_REF], which is no-reference.

Table 2 shows the performance obtained for the cross-dataset evaluation. As can be seen, high correlations are achieved by our method, outperforming some of the compared ones, showing the generalization capability of our metric to predict the quality of unknown PCs. These results show that the proposed method has consistent performance across different validation sets.

CONCLUSION

In this paper, we proposed a novel learning-based no-reference quality metric for point clouds. Unlike previous learning-based approaches, we employ graph convolutions to extract local features at different network layers. Specifically, we consider a two-stream architecture to process in parallel geometry and color distortion, and we use cross-and self-attention to combine features leveraging long-term dependencies in the point cloud. Our proposed method achieves state-of-the-art correlations on two benchmark datasets. In the future, we believe that this approach should be further validated as larger annotated datasets for point clouds quality become available.

Fig. 2 .

 2 Fig. 2. Left: Computing an edge feature e ij , Right: The Edge-Conv operation [16].

Table 1 .

 1 Results obtained on ICIP20 dataset using 6-fold cross validation

	Model	PLCC ↑	SROCC ↑
	po2point MSE	0.946	0.934
	po2plane MSE	0.959	0.951
	PSNR po2point MSE	0.868	0.855
	PSNR po2point HAU	0.548	0.456
	PSNR po2plane HAU	0.580	0.547
	color Y MSE	0.876	0.892
	color Cb MSE	0.683	0.694
	color Cr MSE	0.594	0.616
	color Y PSNR	0.887	0.892
	color Cb PSNR	0.693	0.694
	color Cr PSNR	0.626	0.617
	pl2plane AVG	0.922	0.910
	pl2plane MSE	0.925	0.912
	GraphSim	0.931	0.893
	PCQM	0.796	0.832
	PointNet-SSNR	0.908	0.955
	PointNet-DCCFR	0.947	0.973
	Ablation (TB Graph Aggregation)	0.835	0.862
	Ablation (Color stream)	0.836	0.925
	Ablation (Geometry stream)	0.925	0.939
	Our	0.946	0.973

Table 2 .

 2 Results obtained by training the model on ICIP20 and testing it on PointXR

	Model	PLCC ↑	SROCC ↑
	po2pointMSE	0.887	0.978
	po2planeMSE	0.855	0.942
	PSNRpo2pointMSE	0.983	0.978
	PSNRpo2planeMSE	0.972	0.950
	PointNet-DCCFR	0.981	0.964
	Our	0.967	0.988


Guillaume Morel
email: guillaume.morel@imt-atlantique.fr

Simon Benaichouche
email: simon.benaichouche@imt-atlantique.fr

Nicolas Courty
email: nicolas.courty@irisa.fr

François Rousseau
email: francois.rousseau@imt-atlantique.fr

Turning Normalizing Flows into Monge Maps with Geodesic Gaussian Preserving Flows

Normalizing Flows (NF) are powerful likelihood-based generative models that are able to trade off between expressivity and tractability to model complex densities. A now well established research avenue leverages optimal transport (OT) and looks for Monge maps, i.e. models with minimal effort between the source and target distributions. This paper introduces a method based on Brenier's polar factorization theorem to transform any trained NF into a more OT-efficient version without changing the final density. We do so by learning a rearrangement of the source (Gaussian) distribution that minimizes the OT cost between the source and the final density. The Gaussian preserving transformation is implemented with the construction of high dimensional divergence free functions and the path leading to the estimated Monge map is further constrained to lie on a geodesic in the space of volume-preserving diffeomorphisms thanks to Euler's equations. The proposed method leads to smooth flows with reduced OT costs for several existing models without affecting the model performance. The code is available here https://github.com/morel-g/GPFlow.

Introduction

Modeling high dimensional data is a central question in data science as they are ubiquitous in applications. Various tasks such as probabilistic inference, density estimation or sampling of new data require accurate probabilistic models that need to be defined efficiently. There exists a large variety of generative models in the literature. Among other approaches, diffusion / score based models (Ho et al., 2020;[START_REF] Sohl-Dickstein | Deep unsupervised learning using nonequilibrium thermodynamics[END_REF][START_REF] Song | Generative modeling by estimating gradients of the data distribution[END_REF], variational autoencoders (VAES) [START_REF] Diederik | Auto-encoding variational bayes[END_REF][START_REF] Jimenez Rezende | Stochastic backpropagation and approximate inference in deep generative models[END_REF] and generative adversarial networks (GAN) [START_REF] Goodfellow | Generative adversarial nets[END_REF] are frequent choices, each with their strengths and weaknesses.

Normalizing flows. A fourth popular class of generative models is Normalizing flows (NF). NF models transform a known probability distribution (Gaussian in most cases) into a complex one allowing for efficient sampling and density estimation. To do so they use a diffeomorphism f : R d → R d which maps a target probability distribution µ to the known source distribution ν = f # µ [START_REF] Dinh | Nice: Non-linear independent components estimation[END_REF][START_REF] Rezende | Variational inference with normalizing flows[END_REF]. In practice the flow satisfies the change of variables formula: log p µ (x) = log p ν (f (x)) + log | det ∇f (x)|.

(1)

There are many possible parameterizations of f , usually relying on automatic differentiation to train their parameters via first order optimization algorithms. For density estimation applications, training is done by maximizing the likelihood of the observed data. The data are generally high dimensional and accessing p µ (x) for a given x requires computing determinant of the Jacobian matrix of f . This operation has a complexity of O(d 3 ) in general and thus a requirement for the flow architecture is to have a tractable determinant of the Jacobian while remaining expressive enough [START_REF] Dinh | Nice: Non-linear independent components estimation[END_REF][START_REF] Durk | Glow: Generative flow with invertible 1x1 convolutions[END_REF][START_REF] Rezende | Variational inference with normalizing flows[END_REF][START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF].

Optimal transport. A diffeomorphism transforming any well-behaved distribution into another always exists in theory [START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF]. However, there can be many ways to transform one probability measure µ into another probability measure ν, and therefore the function f is generally not unique. This has led to many proposed architectures in the literature [START_REF] Durk | Glow: Generative flow with invertible 1x1 convolutions[END_REF][START_REF] Grathwohl | Ffjord: Free-form continuous dynamics for scalable reversible generative models[END_REF][START_REF] Huang | Neural autoregressive flows[END_REF][START_REF] De Cao | Block neural autoregressive flow[END_REF][START_REF] Papamakarios | Masked autoregressive flow for density estimation[END_REF]. The question of choosing the "best" transformation among all existing ones is therefore crucial, independently of how accurately µ models the target distribution. One way to make the architecture unique (under appropriate conditions on the two distributions) is to use optimal transport [START_REF] Hamfeldt | Optimal transport. Youtube videos[END_REF][START_REF] Peyré | Computational optimal transport[END_REF][START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF][START_REF] Villani | Optimal Transport: Old and New[END_REF], that is to choose the one giving the Wasserstein distance between µ and ν, with a squared L 2 ground cost:

W 2 2 (µ, ν) = min f R d |f (x) -x| 2 dµ(x), ν = f # µ, (2)
where f # µ denotes the push forward operator of µ through f . An optimal model in the sense of (2) minimizes the total mass displacement which can be a desirable property even if it is often a difficult task. In particular Brenier's theorem [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] states that the optimal function f is the gradient of a scalar convex function, which is widely used when solving (2).

One key property of OT mappings is that they should better preserve the structure of the distribution compared to non OT transformations. This makes them particularly appealing for machine learning applications, and may also help with generalization performance [START_REF] Karkar | A principle of least action for the training of neural networks[END_REF].

Optimal transport and NF models. Including OT in NF models has recently received much attention with various approaches to obtain a map g which satisfies the property (2). Among all these methods, many use either directly Brenier's theorem [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] or the dynamic OT formulation with the Benamou-Brenier approach [START_REF] Benamou | A computational fluid mechanics solution to the monge-kantorovich mass transfer problem[END_REF]. One important remark is that most of the approaches considered need dedicated architectures in order to satisfy the OT property. For example the transformation f is often written as neural network modeling the gradient of a (possibly convex) scalar function [START_REF] Amos | Input convex neural networks[END_REF]Finlay et al., 2020a;[START_REF] Huang | Convex potential flows: Universal probability distributions with optimal transport and convex optimization[END_REF][START_REF] Onken | Ot-flow: Fast and accurate continuous normalizing flows via optimal transport[END_REF][START_REF] Zhang | Monge-ampere flow for generative modeling[END_REF]. This sometimes requires some particular training process (Finlay et al., 2020a;[START_REF] Huang | Convex potential flows: Universal probability distributions with optimal transport and convex optimization[END_REF][START_REF] Onken | Ot-flow: Fast and accurate continuous normalizing flows via optimal transport[END_REF] and/or the addition of some penalization terms in the loss function [START_REF] Onken | Ot-flow: Fast and accurate continuous normalizing flows via optimal transport[END_REF]Finlay et al., 2020b;[START_REF] Yang | Potential flow generator with l2 optimal transport regularity for generative models[END_REF]. When considering the Benamou-Brenier formulation, the normalizing flow is interpreted as the discretization of a continuous ordinary differential equation (Chen et al., 2018b) and the optimal transport problem is then solved dynamically (Finlay et al., 2020b;[START_REF] Onken | Ot-flow: Fast and accurate continuous normalizing flows via optimal transport[END_REF][START_REF] Zhang | Monge-ampere flow for generative modeling[END_REF].

The methods mentioned above require to constrain the architecture and/or the training procedure which can be problematic in some cases and has already motivated other OT related approaches [START_REF] Uscidda | The monge gap: A regularizer to learn all transport maps[END_REF]. In this paper we present a method which achieve the optimal map of a given NF without constraining the architecture or the training procedure. Polar factorization. An overlooked implication of Brenier's theorem is the so-called polar factorization theorem, that states that the optimal transport map ∇ψ solving (2) can be factorized into the composition of two functions ∇ψ = s • f , the function f being some arbitrary smooth map from µ to ν and s an associated measure preserving function of ν [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]. The idea we exploit is the possibility, from a given flow f (and its corresponding inverse g = f -1 ), to rearrange the distribution ν using s to obtain a new map reducing the OT cost without changing the distribution given by the push-forward µ = g # ν = (g • s -1 ) # ν. The OT-improved map can then be obtained with the composition g • s -1 . Interestingly this property is not as popular as the previous one and to our knowledge is not used when dealing with OT and NF models. Yet normalizing flows can take advantage of this formulation mostly because the distribution ν is known and simple (here and in the following ν is a standard normal) which makes it possible to construct architectures preserving ν. [START_REF] Durk | Glow: Generative flow with invertible 1x1 convolutions[END_REF][START_REF] Grathwohl | Ffjord: Free-form continuous dynamics for scalable reversible generative models[END_REF][START_REF] Huang | Neural autoregressive flows[END_REF][START_REF] De Cao | Block neural autoregressive flow[END_REF][START_REF] Papamakarios | Masked autoregressive flow for density estimation[END_REF] each with their pros and cons which sometimes depend specifically on the test case considered. Our idea is to use Brenier's polar factorization theorem to rearrange the points in the known distribution to obtain the optimal map associated with a given flow. We consider the most common case where the known distribution is a standard normal and call such rearranging maps Gaussian Preserving (GP) flows, see Figure 1. An important point is that by construction our GP map will only change the OT cost of the model. The target density and therefore the training loss given by the model will stay the same. This allows us to take any pre-trained model and compute the associated Monge map, thus improving the model in terms of OT displacement from the source to the target distribution, without changing the modeled density see Figure 2.

Main contributions

Construction of divergence free functions in high dimensions. We show that divergence free functions can be used to model GP flows. We therefore derive an effective construction of divergence free functions in high dimensions and apply it in the latent space of some popular VAE models.

Euler's equations. Since several GP flow models can solve the same OT problem, we also look for a way to find the "best" GP flow. This is somehow similar to the approaches from Finlay et al. (2020b); [START_REF] Onken | Ot-flow: Fast and accurate continuous normalizing flows via optimal transport[END_REF] where the trajectories of a continuous normalizing flow are penalized to be straight lines. This is not strictly needed to find the Monge map but can be interpreted as some geodesic over all the flows which solve the associated OT problem. In this work, we show that the geodesics associated with the OT problem are actually given by solutions to the Euler equations, following a celebrated result by Arnold [START_REF] Arnold | Sur la géométrie différentielle des groupes de lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF].

The penalization of Euler's equations in high dimensions and its practical implementation is therefore also considered, which is to the best of our knowledge an original contribution.

Data structure preservation with optimal transport. Finally we show one potential interest of GP flows by studying the preservation of the data structure experimentally. More specifically we focus on the preservation of disentanglement on the dSprites [START_REF] Matthey | dsprites: Disentanglement testing sprites dataset[END_REF], MNIST [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] and Chairs [START_REF] Mathieu | Seeing 3d chairs: exemplar part-based 2d-3d alignment using a large dataset of cad models[END_REF] datasets in some variational auto-encoder (VAE) latent space. On this particular example we show that OT allows to improve the preservation of the structure of the latent data points which is otherwise destroyed when applying the NF model.

Polar factorization theorem

The main idea is to use the Brenier's polar factorization theorem to construct the Monge map with a rearrangement of the known probability distribution ν.

Notation and measure preserving definition. In the following we will use the shortened notation (X , ν) for a probability space (X , B, ν) where X is a set, B a σ-algebra of its subset and ν a probability distribution.

We recall the definition of measure preserving function use in [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]: a measure preserving mapping from a probability space (X , ν) into itself is a mapping s : X → X such that for every ν measurable subset A of X , s -1 (A) is ν-measurable and ν(s -1 (A)) = ν(A).

Theorem 1 (Brenier's polar factorization [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]). Let (X , ν) be a probability space, X bounded. Then for each g ∈ L p (X , ν, R d ) satisfying the non degeneracy condition

ν(g -1 (E)) = 0 for each Lebesgue negligible subset E of R d ,
there exists a unique convex function ψ : X → R and a unique measure preserving function s : X → X such that g(s(x)) = ∇ψ(x), and s(x) minimizes the cost X |g(s(x)) -x| 2 dν(x).

Note that in his work Brenier mention the maximization of x • g • s(x)dν(x) which is equivalent to the minimization with repect to s of the quadratic cost thanks to the measure preserving property of s. Our goal is to leverage the polar factorization theorem in order to solve the OT problem between ν and µ := g # ν where g is given and ν = N (0, Id), by looking for the rearrangement s via an optimization problem. To do so we need to construct a class of measure preserving maps. Remark 1. Since in practice we consider ν to be a standard normal, the domain X is not bounded and therefore does not strictly satisfy the hypothesis of Theorem 1. We do not investigate this point further and simply quote a remark from Brenier's work [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]: "we believe that the result is still true when X is unbounded, provided that p > 1 and X x q β(x)dx < +∞, where 1/q + 1/p = 1". The function β(x) = e -x 2 /2 is the probability density of ν, and the inequality is therefore satisfied.

Gaussian preserving flows

In order to apply Brenier's polar factorization theorem, it is therefore needed to construct a class of measure preserving maps. Since we consider the case where ν is a standard normal, we call such maps Gaussian preserving (GP). All proofs of the propositions and lemmas are given in Appendix B.

Consider two probability measures α and β with density h α and h β respectively. A map s is measure preserving between α and β if it satisfies the change of variable equality (same as (1) without the log)

h α (x) = h β (s(x))| det(∇s(x))|.
In our case, we want s to be Gaussian preserving therefore h α = h β = e -x 2 /2 and one gets

| det ∇s(x)| = e ( s(x) 2 -x 2 )/2 . ( 3
)
It turns out that Lebesgue preserving functions (i.e. satisfying | det ∇φ| = 1) can be used to construct maps satisfying (3). In the following we will denote erf : R d → R d the distribution function of a one dimensional Gaussian (that is erf(x) = 2 √ π

x 0 e -t 2 dt) applied component wise.

(x) = √ 2 erf -1 •φ • erf ( x √ 2 ), x ∈ R d .
From now on we will focus on the construction of volume and orientation preserving maps (i.e. satisfying det ∇φ = 1) since functions satisfying det ∇φ = -1 can be constructed from them see Appendix B.1.2. Moreover, one has the following result regarding the regularity of GP flows.

Lemma 1. Assume the Monge map and the NF architecture are C 1 diffeomorphisms. Then the corresponding GP flow s is C 1 , the associated function φ is also C 1 and either satisfies det ∇φ(x) = 1 everywhere or det ∇φ(x) = -1 everywhere.

Volume-orientation preserving maps

First we introduce the space SDiff(Ω) we will working with from now on. Let Diff(Ω) be the set of all diffeomorphisms in Ω then SDiff(Ω) := ψ ∈ Diff(Ω) | det(∇ψ)(x) = 1, ∀x ∈ Ω , where Ω = (-1, 1) d . That is we need a transformation which satisfies two properties: 1) the function must be volume and orientation preserving, 2) the solution must stay in the domain (-1, 1) d . Consider the following ODE:

d dt X(t, x) = v(t, X(t, x)), x ∈ Ω, 0 ≤ t ≤ T, X(0, x) = x. (4)
We impose two conditions on the velocity v:

∇ • v = 0, in Ω, ( 5
) v • n = 0, on ∂Ω, ( 6
)
where n is the outward normal at the boundary of Ω. We define φ to be the solution at the final time φ(x) := X(T, x). Property (5) implies that det ∇φ = 1 (this can be checked with the formula d dt det ∇X = div(v) det ∇X), and property (6) ensures that φ has values in Ω. Any function in SDiff(Ω) can be written as a solution to (4) for d ≥ 3 [START_REF] Shnirelman | Attainable diffeomorphisms. Geometric and functional analysis[END_REF], for d = 2 some pathological cases can be constructed [START_REF] Shnirelman | Generalized fluid flows, their approximation and applications[END_REF].

Divergence free vector fields. First we focus on the vector fields satisfying (5) for arbitrary large dimensions. Property (6) can then be incorporated with very little additional work.

Proposition 2. Consider an arbitrary vector field

v : R d → R d . Then ∇ • v = 0 if and only if there exists smooth scalar functions ψ i j : R d → R, with ψ i j = -ψ j i such that v i (x) = d j=1 ∂ xj ψ i j (x), i = 1, ..., d, ( 7
)
where v = (v 1 , ..., v d ).
To impose the boundary conditions (6) one can simply multiply each ψ i j by (x 2 i -1)(x 2 j -1). Lemma 2. Consider the coefficients

ψ i j (x) = h i (x i )h j (x j ) ψ i j (x) ( 8
)
where h i (x i ) = h 1 i (x i -1)h 2 i (x i + 1), h 1 i , h 2 i are functions satisfying h 1 i (0) = h 2 i (0) = 0 and ψ i j (x) : R d → R are bounded functions satisfying ψ i j = -ψ j i . Then the function v defined in Proposition 2 satisfies ∇ • v = 0 and v • n = 0 on ∂Ω.
The function h = (h 1 , ..., h d ) can typically be parametrized as neural networks with no bias but in our experiments we have simply chosen to consider the case h i (x i ) = (x 2 i -1). One drawback of Lemma 2 is that it does not guarantee a universal approximation of divergence free functions near the boundaries. Note that when we are away from the boundaries however Proposition 2 gives this universal approximation result. In our experiments we observed that we can at least significantly reduce the OT cost with the construction (8) and we therefore let the investigations related to the boundary conditions for future work.

The incompressible property (5) and the boundary conditions ( 6) can be exactly implemented in the network in any dimension. Note however that in order to get all the incompressible vector fields (7), we need to construct at least d(d -1)/2 arbitrary scalar functions. See Appendix A for the practical construction of these divergence free functions in high dimensions.

Euler's geodesics

GP flows give a way to compute the Monge map for any trained NF architecture. Many transformations can achieve this goal and the question of finding the best flow among all volume preserving transformations needs to be considered. In the following we will consider a regularization term to smooth the trajectories by minimizing the energy v 2 . While it could be tempting to try to penalize directly v 2 , the energy will unfortunately have an opposite objective from the minimization of the OT cost. Indeed the global minimum for the energy is v 2 = 0 (that is the particles do not move) which is obviously not the velocity field which minimizes the OT cost. In fact the global minimum of a loss composed with the two terms OT + v 2 may not minimizes the OT cost. On the contrary, as we will detail below, Euler's equations minimize the energy on SDiff over all other solutions with the same initial and final states. By finding the correct final state and solving Euler equations it is therefore possible to obtain a solution which both minimizes the OT cost and the energy.

Arnold's theorem

In 1966, Arnold [START_REF] Arnold | Sur la géométrie différentielle des groupes de lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF] showed that the flow described by Euler's equations coincides with the geodesic flow on the manifold of volume preserving diffeomorphisms. This theoretical result therefore gives the reason why regularizing our flows with Euler's equations is a desirable property. Mainly that Euler's equations take the path with the lowest energy to reach the final configuration. Consider the Euler equations:

         ∂ t v + (v • ∇)v = -∇p, t ∈ [0, T ], x ∈ Ω, ∇ • v = 0, t ∈ [0, T ], x ∈ Ω, v • n = 0, t ∈ [0, T ], x ∈ ∂Ω, v(0, •) = v 0 , (9)
where v := v(t, x) is the velocity field, p := p(t, x) the pressure and n := n(x) the outward normal at the boundary of Ω. Here the pressure p ensures that ∂ t v + (v • ∇)v can be written as the gradient of some scalar function which is uniquely defined (up to a constant) thanks to the additional divergence free and boundary conditions on v. Additionally the pressure field can be interpreted as the Lagrange multiplier of the divergence free constraint for the associated variational formulation of Euler's equations. In particular, it may not be needed to compute p when solving numerically (9). We introduce E the energy of a smooth function X(t, •):

E(X) = T 0 Ω 1 2 |∂ t X(t, x)| 2 dxdt, ( 10
)
Now assume φ ∈ SDiff(Ω). Arnold's problem's consists in finding the path X(t, •) t∈[0,T ] in SDiff(Ω) joining the identity to φ which minimizes E:

min X(t,•)∈SDiff(Ω) E(X), X(0, •) = Id, X(T, •) = φ(•). ( 11
)
In other words ( 11) is the geodesic in SDiff(Ω) between Id and φ. Theorem 2 [START_REF] Arnold | Sur la géométrie différentielle des groupes de lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits[END_REF]). Assuming the existence of a solution to Arnold's problem, X is solution to (11) if and only if v(t, x) := ∂ t X(t, x) satisfies Euler's equations (9).

Penalization of Euler's equations in high dimensions

Numerical schemes developed to efficiently solve the Euler equations [START_REF] Canuto | Spectral Methods[END_REF][START_REF] Quarteroni | Numerical models for differential problems[END_REF] (mainly for fluid mechanics problems, i.e. for dimensions up to 3) scale badly when the dimension increases.

In this work, the solution to Euler's equations is interpreted as the geodesic to reach the solution of the OT problem and the dimension can be arbitrary large. Therefore we approach the equation (9) through a penalization procedure which can be carried out in any dimension. As explained in the previous section we notice that the second and third equations in ( 9) are satisfied by construction in the network.

Our remaining goal is to constrain the network to be a smooth solution to

∂ t v + (v • ∇)v = -∇p.
The left hand side can therefore be written as the gradient of a scalar function and we note that a vector satisfies w t,x ∈ R d satisfies w t,x = ∇p(t, x) if and only if its Jacobian is symmetric ∇w t,x = (∇w t,x ) T . In order to solve the first equation in ( 9), we propose to penalize the non-symmetric part of the Jacobian for the total derivative of v. Since a Jacobian-vector product can be efficiently evaluated in high dimensions (unlike the calculation of the full Jacobian which is computationally expensive), we do not calculate directly the Jacobian and use instead the following property of symmetric matrices: M is symmetric if and only if

y T M z -z T M y = 0, ∀y, z ∈ R d .
The idea is to sample random vectors y, z during the training and to penalize this term for the total derivative, that is to minimize:

R(x) := E y,z T 0 y T (∇w t,x )z -z T (∇w t,x )y) 2 dt , ( 12
)
with y, z ∼ N (0, Id) and

w t,x = ∂ t v + (v • ∇)v.
In practice, we do not compute the full time integral in (12) as it would be computationally too expensive but calculate the penalization only at our time steps discretization.

Approximation of the total derivative. To reduce the computational burden, we do not calculate exactly the total derivative w t,x but use an approximation of its Lagrangian formulation instead. More precisely, consider the variable X(t, x) from (4) that is the position of a particle at time t with initial position x. We recall the equality (see Appendix B.3.1) D Dt v(t, X(t, x)) = ∂ t v(t, X(t, x)) + (v(t, X(t, x)) • ∇)v(t, X(t, x)) and therefore choose to approximate the right-hand side by using a first order Taylor expansion of Dv/Dt:

D Dt v(t, X(t, x)) ≈ v(t n+1 , X(t n+1 , x)) -v(t n , X(t n , x)) ∆t , ( 13
)
∆t := t n+1 -t n . In practice, ∆t is set to 2 √ ε where ε is the machine precision. This approximation can be easily computed since it requires only the evaluation of the velocity at two positions of a particle.

To summarize our approach requires to penalize a Jacobian-vector product of the form y T (∇w)zz T (∇w)y where w is given by ( 13). Penalizing Jacobian-vector product has already been done in other contexts and prove to efficiently scale with the dimension [START_REF] Song | Sliced score matching: A scalable approach to density and score estimation[END_REF].

Procedure

In order to solve the optimal transport problem, we can either use the forward NF function f or g := f -1 . Depending on this choice the loss function is then either

E µ(x) x -s • f (x) 2 or E ν(x) x -g • s(x) 2 .
In practice the GP flow is parametrized as a standard residual network (ResNet) with a Runge-Kutta 4 time discretization [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF] (other discretization are possible) and is estimated by minimizing the parameters of the velocity field. When regularizing with Euler's equations, we replace the term ∇w in ( 12) by ( 13) and calculate the Jacobian-vector product with the function torch.autograd from pytorch. A parameter λ > 0 is also added in front of the penalization term that is if we use the function f :

min θ E µ(x) x -s θ • f (x) 2 + λR θ • erf • f (x) √ 2 , ( 14
)
or if the function g is considered instead

min θ E ν(x) x -g • s θ (x) 2 + λR θ • erf • x √ 2 , ( 15
)
where the vector θ denotes the parameters of the velocity field v, f is the NF architecture, s the GP flow and R corresponds to the term penalized with Euler's equations. The subscript θ has been added to highlight the dependence of s and R to the parameters. We emphasize that ( 14)-( 15) are two distinct optimization strategies:

• If we do not want to invert the NF model (for example if it is computationally expensive) we can simply minimize ( 14) over the training data. The probability distribution µ in ( 14) is then the unknown distribution from which the data are taken. This assumes however that the training data are correctly mapped to the standard distribution ν with f .

• If the NF model is cheap to invert we can minimize (15) instead. To do so we optimize over samples from the standard normal distribution ν. In this case, the probability distribution µ is the transformation of ν by the inverse of the NF model µ = g # ν.

If the NF model transforms perfectly the training data over the standard normal ν, these two approaches are equivalent. If this is not the case, we notice that the second approach requires a cheap inverse of the NF model, but has the advantage of not using any training data to train the GP flow and may therefore better generalize.

Results

We apply GP flows on two popular NF models: BNAF, a discrete NF [START_REF] De Cao | Block neural autoregressive flow[END_REF] for two-dimensional test cases and FFJORD, a continuous NF [START_REF] Grathwohl | Ffjord: Free-form continuous dynamics for scalable reversible generative models[END_REF] for higher dimensional cases. Both of these models are solid references among NF and do not incorporate any OT knowledge in their architecture or training procedure. The codes are taken from the official repositories1 . The FFJORD model has an inverse function directly available in the code, which is not the case for the BNAF model. For this reason we consider only the FFJORD model when interpolating in the latent space of the dSprites and MNIST datasets because interpolations require the NF architecture to have an inverse function available. To compare our results we consider the CP-Flow architecture [START_REF] Huang | Convex potential flows: Universal probability distributions with optimal transport and convex optimization[END_REF]. The CP-flow network is constrained to be the composition multiple blocks which are gradient of scalar convex function and therefore converges by construction towards the optimal map when considering only 1 block (provided the optimization reaches a global minimum). This makes CPFlow a good candidate for comparison.

Density estimation on toy 2D data

In this section we perform density estimation on several 2d standard toy distributions [START_REF] Grathwohl | Ffjord: Free-form continuous dynamics for scalable reversible generative models[END_REF][START_REF] Wehenkel | Unconstrained monotonic neural networks[END_REF]. In particular we train our GP flow on a pre-trained BNAF model. We provide two dimensional toy examples for the eight Gaussians, two moons and pinwheel test case on Figure 3. In the case of the pinwheel dataset we use Euler regularization see below for more details. To compute the exact discrete OT projection we use the POT library [START_REF] Flamary | Pot: Python optimal transport[END_REF]. The distribution is colored to compare the transformation with the exact OT map. We observe on each experiment that adding GP flow makes the transformation closer to the Monge map.

Euler regularization. Let's briefly recall the goal when considering Euler's regularization. As discussed in Theorem 1 the measure preserving transformation s which minimizes the OT cost is unique. In this work we have constructed s to be the solution of an ODE with divergence free velocity and therefore even if s (that is the solution of the ODE at the final time) is unique, there are infinitely many trajectories which reach this final configuration. The role of Euler regularization is to obtain smooth trajectories for our ODE, a property which may be helpful during the training process to minimize the OT cost.

To highlight the value of using the Euler's penalization we focus on the pinwheel test case. If the pinwheel experiment from Figure 3 was done with Euler regularization, this is because we notice a bad convergence of the GP flow without it. Figure 4 shows the reason why the training has trouble to converge: the trajectories are very irregular when no penalization is applied. This is an example of the practical utilization of Euler's regularization and as shown in Figure 5 it leads to a better global convergence. We give another illustration of the trajectories for the two moons test case on Figure 6. We observe that the trajectories obtained with Euler regularization are more structured.

Note that Euler regularization is not always needed to reach the OT map and in practice it seems to depend on the case considered. More details on Euler penalization are given in appendix D.

An example of data structure preservation: improving disentanglement preservation with optimal transport

The goal here is to show that OT transformations better preserve the data structure than non OT ones.

Preserving data structure can be very useful for example when trying to interpolate between data points. Indeed, interpolating directly in the data space often leads to poor results mainly because we may interpolate through low density regions. In this case it can be a good idea to map the data density into a Gaussian, interpolate in the Gaussian space and then apply the reverse transformation. However, for the interpolation to be meaningful we expect the NF architecture to preserve as much of the data structure as possible. In order to study high dimensional data structure we therefore seek a data space which satisfy two properties:

-The initial data space should be very structured.

-The interpolations between data points should be easily comparable and in particular we should be able to classify them easily.

Disentangled representations satisfy these two properties and allow to encode the data in a latent space where change in one direction result in the change over one generative factor in the data. Recently the construction of variational auto-encoders (VAE) (Kingma & Welling, 2014) with disentangled latent space has received much attention [START_REF] Higgins | beta-vae: Learning basic visual concepts with a constrained variational framework[END_REF][START_REF] Christopher P Burgess | Understanding disentangling in beta-vae[END_REF]Chen et al., 2018a;[START_REF] Kim | Disentangling by factorising[END_REF]. We therefore propose to experimentally study disentanglement preservation of NF with and without OT. To do this we apply a NF architecture (FFJORD in our case) to the VAE's latent distribution and consider the addition of a GP flow. The goal here is therefore to get as close as possible to the initial interpolations in order to show that even if the initial density is transformed to a Gaussian the data structure is preserved at least to some extent.

For the disentangled interpolation in the NF target (Gaussian) space we consider the same directions which are present in the VAE latent space and are aligned with the axes. This may be a little naive since these directions could be slightly modified depending on the source and target distributions, but it seems a good enough approximation here. Our experiments are run with the β-TCVAE architecture (Chen et al., 2018a) and the latent space dimension is 10.

dSprites dataset. The dSprites dataset [START_REF] Matthey | dsprites: Disentanglement testing sprites dataset[END_REF] is made of 64 × 64 images of 2D shapes procedurally generated from 5 ground truth independent latent factors. These factors are shape, scale, rotation, x and y positions of a sprite. Since the factors are known, we can compute a quantitative evaluation of disentanglement and we choose here to consider the metric from [START_REF] Eastwood | A framework for the quantitative evaluation of disentangled representations[END_REF] on the continuous factors (i.e. all the factors except the shape) for the three criteria: disentanglement, completeness and informativeness. Table 1 shows that FFJORD destroys the latent structure and gives the worst disentanglement, completeness and informativeness scores. Adding GP flows allow to recover the same disentanglement score as the initial latent space and get values closer both for completeness and informativeness. On Table 2 the OT costs are compared, and as expected GP flows allow to reduce the OT cost without changing the loss. Interestingly GP flows with no additional regularization do not converge completely to the Monge map because particles get out of the domain at some point, making it impossible to continue the training process (the training is stopped at ≈ 60 epochs). We conjecture that this may be due to non-smooth trajectories of our GP flows and a regularization is therefore needed. As shown on Table 2 adding Euler regularization fixed this issue and allows to further reduce the OT cost.

Model

To illustrate the preservation of disentanglement, some interpolations are also presented in Figures 9 and10 in Appendix E. The dimensions are sorted with respect to their KL divergence in the initial latent space and therefore only the first dimensions carry information: each of the first 5 lines correspond to a generative factor while the last dimensions leave the image unchanged. This structure is lost when mapping the latent space to a Gaussian distribution with the FFJORD architecture. The addition of a GP flow fixed this issue and the interpolation better match the initial latent one.

MNIST dataset. We also consider the MNIST dataset [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. As opposed to the dSprites test case, GP flows do not seem to have trouble to converge here without Euler penalization and therefore we obtain comparable OT costs with and without Euler regularization see Table 2. Since the generative factors are not known in this case we cannot make a quantitative evaluation of disentanglement as we did for the dSprites dataset. We focus instead on the interpolations presented in Figures 11 and12 in Appendix E. GP flows better preserve the data structure of the initial latent space compare when applying only the FFJORD model. This can be seen in particular on the last two rows of each block which are not changing in the initial latent space. This structure is lost with the FFJORD model and recovered when training a GP flow.

Chairs dataset.

Finally we look at the chairs dataset [START_REF] Mathieu | Seeing 3d chairs: exemplar part-based 2d-3d alignment using a large dataset of cad models[END_REF]. Again, the GP method converges both with and without Euler's penalization and greatly reduced the OT cost see Table 2. Interpolations are presented in Figures 13 and14 in Appendix E. By looking at the last rows we observe once again that GP flows better preserve the data structure of the initial latent space compare when applying only the FFJORD model.

Comparison with CPFlow.

To check that GP flows get close to the Monge map we also compare the OT costs with the CPFlow architecture. Note that the number of parameters as well as the training procedure differ from the FFJORD+GP and CPFlow architectures so we only use CPFlow OT costs as a baseline to ensure that FFJORD+GP is close to the Monge map. The losses are only given here as an indication to show that the probability distribution of both FFJORD and CP flow are close from each other and a comparison between their OT costs is therefore relevant. Finally note that we both test the CP flow architecture with 1 and 3 blocks. The advantage with the 1 block architecture is that CP flow should theoretically converges to the Monge map. However it seems harder to make the network converges in this case and we therefore also show some comparisons with 3 blocks which leads to much lower losses.

GP flows get OT values which are always lower than the CP flow ones for similar losses and we therefore conclude that GP flows are at least as close to the Monge map as CPFlow. The only exception is the chairs dataset where CP flow with one block has a lower OT cost. We note however that in this case the loss for CP flow is much bigger and therefore the comparison of OT costs with GP flow may not be relevant. We also notice that CPFlow may leads to high OT cost when trained on the MNIST dataset. One possible explanation for this is that CPFlow is trained only on the data. Since the OT costs are evaluated with random samples drawn from the standard normal this may show some issue with the generalization during the backward process. On the contrary GP flows are directly trained on random samples from this distribution and may therefore better generalize.

It is therefore possible to better preserve the data structure with GP flows by significantly reducing the OT cost. Note however that since we were not able to make the CP flow architecture converge with only 1 block in our experiments, it could be a good idea to try to assess more precisely how close GP flow is to the Monge Map. To this end there exists many other approaches which could be considered. One could for example use dedicated OT benchmarks [START_REF] Korotin | Do neural optimal transport solvers work? a continuous wasserstein-2 benchmark[END_REF] or try more robust approaches [START_REF] Korotin | Neural optimal transport[END_REF][START_REF] Makkuva | Optimal transport mapping via input convex neural networks[END_REF].

Discussion

This article describes a method to reduce the OT cost of any pre-trained NF model without changing the estimated target density. The proposed method has been tested up to d = 10 and does not require to constrain the architecture of the original model. The procedure relies on building Gaussian preserving flows to rearrange the source distribution to satisfy the OT property. The proposed approach is based on incompressible vector fields which allow to use a nice interpretation of Euler's equations as a geodesic in the group of volume-preserving diffeomorphisms between the identity and the transformation minimizing the OT cost. This original contribution allows to add a regularity condition to the estimated map in addition to simply enforcing the OT property.

Perspectives.

The numerical experiments presented here pave the way to new research perspectives. First compared to other OT approaches in the NF literature, GP flows is to the best of our knowledge the first one which does not constrain the NF architecture to obtain the Monge map. This could be a great advantage when the NF architecture is already constrained for other reasons (for example to satisfy some symmetries or data-related properties). We believe that in this case GP flows may stand out as it could be difficult to further constrain the network to satisfy the OT property with standard approaches. The proposed approach could also be a starting point to investigate other type of (potentially non OT) costs, and more specifically non-quadratic OT costs, such as the L 1 norm which is much less considered in the literature due to the lack of theoretical foundations (the OT map is not the gradient of a convex function anymore). In a GP-based framework, such extension could be easily implemented since we do not rely explicitly on this property. Finally let's mention the case where we consider a map h between two unkown distribution D 1 and D 2 but none of them is Gaussian. To recover the Monge map associated with h a possibility might be to map D 1 to a Gaussian, makes the rearrangement in the Gaussian space with respect to the quadratic cost of h (that is between D 1 and D 2 ) and then map the points back to D 1 . In this case we could recover the OT map between D 1 and D 2 even if none of them is Gaussian.

Limitations. The main limitation of the proposed method is probably related to the number of independent functions require to construct incompressible vector fields in high dimensions. Indeed, as explained in Proposition 2 one needs to construct at least d(d -1)/2 scalar functions to get all the divergence free vector fields in dimension d and in practice our approach requires d -1 vector valued functions in R d . It could therefore be desirable to scale up more efficiently with the dimension. One possible way to overcome this difficulty would be to give up on the exact implementation of divergence free functions. For example one could add a penalization term of the divergence in the loss with an unbiased estimator of the divergence [START_REF] Song | Sliced score matching: A scalable approach to density and score estimation[END_REF]. Also, let us mention that while Euler's regularization adds nice properties to the flows considered, it also increases the computational time required to train the model. Finally the total duration of the training can be impacted by the initial NF chosen since its evaluation is required to compute the OT cost. NF architectures with fast forward (or backward) pass should therefore be preferred.

A Practical construction of incompressible vector fields in high dimensions

The goal here is to have a GPU-friendly construction of the incompressible vector fields given in Proposition 2. In the following we consider stationary divergence free functions but the time variable can be added with no additional work simply by considering functions in R d+1 instead of R d (the gradients are still taken only on the space variables though). All the proofs are given in Appendix B.

Notations.

Regarding the notations we will use the operator diag for two distinct cases: 1) When w is a vector diag(w) denotes the diagonal matrix obtained from the vector w. 2) When W is a matrix diag(W ) denotes the vector obtained from the diagonal of W . The operation • denotes the scalar product between two vectors. Finally we have adopted the convention that when a scalar multiplies a vector it multiplies each of its component.

Practical construction. Let u n : R d → R d-n . We construct a divergence free function with the functions ψ i j defined as

(ψ i j ) n = u n i-n -u n j-n , if i, j ≥ n + 1, 0, otherwise. (16)
To construct this divergence free function we define the matrix (∇u(x)) n ∈ R d×d and the vector

1 n ∈ R d as (∇u(x)) n ij = ∂ i-n u n j-n , if i, j ≥ n + 1 0, otherwise. , ( 17
)
1 n i = 1, if i ≥ n + 1 0, otherwise. . Lemma 3. Let n ∈ N, n ≤ d -2 and consider the function u n : R d → R d-n . Then the vector field v n : R d → R d defined as v n (x) = (∇u) n 1 n -[diag(∇u) n • 1 n ]1 n , ( 18
)
is divergence free.

To construct the functions (18) we need 1) to compute the product between the Jacobian of a vector valued function and a constant vector 2) sum the diagonal elements of the Jacobian matrix. Both of these operations can be done efficiently on GPU. Note that in order to satisfy the boundary conditions v n • n = 0 one can modify the equation ( 18) as in Lemma 2 (recalling we are taking in practice

h i = x 2 i -1) to obtain v n (x) = (x 2 -1) 2M n x + (∇u) n (x)(x 2 -1) -((x 2 -1) • diag(∇u(x)) n )1 n , ( 19
)
where M n ij = u n i -u n j , if i, j ≥ n + 1 and M n ij = 0 otherwise. It is possible to recover all the incompressible functions from Proposition 2 by adding the blocks

v 0 + v 1 + ... + v d-2 . Proposition 3. Let v be a divergence free function in R d . Then there exists d -1 functions v 0 , ..., v d-2 constructed as in Lemma 3 such that v(x) = d-2 n=0 v n (x).
The attentive reader would have noticed that with the vector functions u n , n = 0, ..., d -2 we have a total of (d + 2)(d -1)/2 independent scalar functions while Proposition 2 only requires the construction of d(d -1)/2 scalar functions leaving d -1 additional functions which are not strictly needed to obtain the divergence free vectors. This is due to the vectorized constructions ( 18)-( 19) which allow a fast evaluation of the divergence free functions on GPU. Having d -1 additional scalar functions in return is not a big issue since the general order remains O(d 2 ). Also note that the practical implementation of the equations ( 16)-( 17) requires to find a pythonic way to efficiently pad a group of matrices with different dimensions. We have not yet find such way and therefore have simply chosen in our applications to construct d -1 vector valued functions u n ∈ R d , n = 0, ..., d -2 and fill the appropriate dimensions in ( 16)-( 17) with 0. Again even if not optimal this is not a big issue as it multiplies the number of independent scalar functions by a factor 2, but the general order remains O(d 2 ).

In practice, we have written the functions u n as the output of a big function u : R d → R (d-1)×d allowing to evaluate all the functions u n in a single pass. The vector u is written as the composition of linear functions with some simple non-linearity

u(x) = M n x n + b n , x i = σ(M i-1 x i-1 + b i-1 ), i = 1, ..., n -1, x 0 = x, (20)
where M i are rectangular matrices, b i a vector field and typically we have taken σ = tanh. One big advantage of the formulation ( 20) is that the Jacobian of u (and therefore of all the functions u n ) can be computed analytically ∇u(x) = M n ∇x n ,

∇x i = diag(σ (M i-1 ∇x i-1 + b i-1 ))M i-1 , ( 21
)
for i = 1, ..., n -1. The formulation (21) therefore allows a fast evaluation of the term (18) in particular when summing the diagonal elements of the Jacobian. In our experiments we have noticed that the analytical formulation of the Jacobian (21) was faster than using torch.autograd.

B Technical material

(x) = √ 2 erf -1 •φ • erf ( x √ 2 ), x ∈ R d . ( 22
)
Proof. We recall some basic properties about the distribution function of a one dimensional Gaussian and its inverse. One has for x ∈ R

erf(x) = 2 √ π x 0 e t 2 dt, d dx erf(x) = 2 √ π e -x 2 , d dx erf -1 (x) = √ π 2 e (erf -1 (x)) 2 . ( 23
)
Consider the function φ : (-1, 1) d → (-1, 1) d defined as

φ(x) = erf • s √ 2 • √ 2 erf -1 (x). ( 24
)
The goal here is to show that | det ∇φ| = 1 then equation ( 22) will follows from (24). By definition

| det ∇φ(x)| := | det ∇ erf • s √ 2 • √ 2 erf -1 (x) |.
Applying the equalities ( 23) component wise and denoting x = (x 1 , ..., x d ) one gets

| det ∇φ(x)| = | i 2 √ Π e -( s √ 2 • √ 2 erf -1 (x)) 2 i × det ∇s( √ 2 erf -1 (x)) √ 2 × i √ 2 √ Π 2 e (erf -1 (x)) 2 i |,
since the determinant of the composition is the product of the determinants. That is

| det ∇φ(x)| = | i e -( s √ 2 • √ 2 erf -1 (x)) 2 i × det ∇s( √ 2 erf -1 (x)) × i e (erf -1 (x)) 2 i |,
which can be written

| det ∇φ(x)| = e (-s( √ 2 erf -1 (x)) 2 + √ 2 erf -1 (x) 2 )/2 ×| det ∇s( √ 2 erf -1 (x))|. Finally using | det ∇s(x)| = e ( s(x) 2 -x 2 )/2 one obtains | det ∇φ(x)| = 1.
Following the same kind of arguments one shows the other way around that is if s

(x) = √ 2 erf -1 •φ • erf ( x √ 2 ) with | det ∇φ| = 1 then | det ∇s(x)| = e ( s(x) 2 -x 2 )/2 that is s is Gaussian preserving.

B.1.2 Orientation reversing functions

In the following Lemma we prove that very orientation reversing function satisfying det ∇ψ = -1 can be written as the composition of a volume and orientation preserving function and the function h(x 1 , ..., x d ) = (-x 1 , x 2 , x 3 , ..., x d ). Lemma. Assume ψ : R d → R d is a function satisfying det ∇ψ = -1 and let h(x 1 , ..., x d ) = (-x 1 , x 2 , x 3 , ..., x d ). Then there exists a volume and orientation preserving function φ such that ψ = φ • h.

Proof. We define φ as φ = ψ • h. This function is indeed volume and orientation preserving since it satisfies det ∇φ = det ∇ψ det ∇h = 1. By noticing h • h = Id one gets ψ = φ • h.

B.1.3 Proof of Lemma 1

Lemma. Assume the Monge map m and the NF architecture g are C 1 diffeomorphisms. Then the corresponding GP flow s is C 1 , the associated function φ is also C 1 and either satisfies det ∇φ(x) = 1 everywhere or det ∇φ(x) = -1 everywhere.

Proof. The definition of s := m • g -1 ensures that s is indeed C 1 . Moreover since g is invertible either det ∇g > 0 everywhere or det ∇g < 0 everywhere (if det ∇g(x) = 0 this would mean that g is not invertible at this point) and the same argument applies to g -1 and to the Monge map m. Therefore the equality s = m•g -1 implies that either det ∇s > 0 everywhere or det ∇s < 0 everywhere. The equality φ

= erf • s √ 2 • √ 2 erf -1
shows that φ is also C 1 and that the sign of det ∇φ does not change.

B.2 Divergence free functions

B.2.1 Proof of Proposition 2

Proposition. Consider an arbitrary vector field v : R d → R d . Then ∇ • v = 0 if and only if there exists smooth scalar functions ψ i j : R d → R, with ψ i j = -ψ j i such that

v i (x) = d j=1 ∂ xj ψ i j (x), i = 1, ...d, ( 25
)
where v = (v 1 , ..., v d ).
Proof. 1) First we prove that ∇ • v = 0. The divergence of v can be written

∇ • v = i ∂ xi j ∂ xj ψ i j .
Using ψ i j = -ψ j i one has

∇ • v = i     j j<i ∂ xi ∂ xj ψ i j - j j>i ∂ xi ∂ xj ψ j i     , (26)
For the term i j j<i ∂ xi ∂ xj ψ i j on the left hand side one can sum over the index j first instead of the index i that is

i j j<i ∂ xi ∂ xj ψ i j = j i i>j ∂ xi ∂ xj ψ i j .
Injecting this equality in (26) one gets

∇ • v = j i i>j ∂ xi ∂ xj ψ i j - i j j>i ∂ xi ∂ xj ψ j i = 0.
2) Now we prove that every vector field satisfying ∇ • v = 0 can be written under the form (25). The proof from Stephen Montgomery-Smith is available online2 for completeness we rewrite it here. The proof is made by induction with the following assumption.

Assumption 1. Let k ∈ N, k ≤ d. Given a smooth vector field v such that div k v := k i=1 v i = 0, there exists scalar functions ψ i j : R d → R, 1 ≤ i, j ≤ k with ψ i j = -ψ j i such that v i = j ∂ j ψ i j .
The Assumption 1 is trivial for k = 0. Suppose it is true for k -1 we prove it for k: assume div k v = 0 and let

f 1 (x 1 , ..., x n ) = x1 0 ∂ k v k (ξ, x 2 , ..., x n )dξ. (27) Since ∂ 1 f 1 = ∂ x v k one has ∂ 1 (v 1 + f 1 ) + ∂ 2 v 2 + ... + ∂ k-1 v k-1 = 0.
Thanks to Assumption 1 there exists functions ψ i j with ψ i j = -ψ j i such that

v 1 + f 1 = k-1 j=1 ∂ j ψ 1 j , v i = k-1 j=1 ∂ j ψ i j , for 2 ≤ i ≤ k -1. ( 28
)
Now we define

f 2 (x 1 , ..., x d ) = x1 0 v k (ξ, x 2 , ..., x k-1 , 0, x k+1 ..., x d )dξ - x k 0 f 1 (x 1 , ..., x k-1 , ξ, ..., x d )dξ, (29)
then

∂ k f 2 = -f 1 , (30)
and using ( 27) in ( 29) one gets

∂ 1 f 2 = v k (x 1 , ..., x k-1 , 0, ..., x d ) - x k 0 ∂ k v k (x 1 , ..., x k-1 , ξ, ..., x d )dξ = -v k . ( 31
)
Now we extend the functions

ψ i j , 1 ≤ i, j ≤ k -1 by defining ψ 1 k = -ψ k 1 = f 2 and ψ i k = -ψ k i = 0 for 2 ≤ i ≤ k. Then extending the equations (28) with k one has k j=1 ∂ j ψ 1 j = v 1 + f 1 + ∂ k f 2 = v 1 , k j=1 ∂ j ψ i j = v i ,
for 2 ≤ i ≤ k -1 and where we used the equality (30) in the first equation. Moreover with the definition of the function ψ k j one has

k j=1 ∂ j ψ k j = -∂ 1 f 2 = v k ,
thanks to (31). This proves Assumption 1 for k.

B.2.2 Proof of Lemma 2

Lemma. Let Ω = [-1, 1] d and consider the coefficients Proof. Indeed since ψ i i = 0 there is no ∂ xi term which appear in the sum of (7) for the component v i . The term h i (x i ) can therefore be factored that is

ψ i j (x) = h i (x i )h j (x j ) ψ i j (x) where h i (x i ) = h 1 i (x i -1)h 2 i (x i + 1), h 1 i , h 2 i are functions satisfying h 1 i (0) = h 2 i (0) =
v i = 0 if x i = ±1. Hence v • n = 0 on ∂Ω.

C Experiment details

We run the experiments on two separate GPUs: a NVIDIA Quadro RTX 8000 and a NVIDIA TITAN X.

Our loss is given by the negative log-likelihood (1). The FFJORD model has an inverse function directly available in the code, which is not the case for the BNAF model. Therefore, as explained in the Section 5, the GP flow is trained using f -1 applied on a standard multi-dimensional normal distribution for FFJORD, whereas only the training data are used for BNAF.

Toy datasets. We give the parameters used for the 2D toy experiments in Table 3 3 . We consider a training set of 80K samples and a testing set of 20K samples, 15 time steps with a Runge-Kutta 4 discretization and a GP flow with two intermediate layers of size 15. For the Euler penalization we take the exact same parameters with λ = 5 × 10 -4 which is divided periodically by a factor 2 the number of division is given in table 3.

dSprites, MNIST and chairs datasets. The VAE architecture used in the experiments is taken from the github repository of Yann Dubois [START_REF] Dubois | Disentangling vae[END_REF]. The parameters used for the GP flows are given in Table 4 and the architectures of CPFlow are given in Table 5. For all test cases we consider a GP flow with 15 time steps, three intermediate layers of 50 parameters for the dSprites dataset and four intermediate layers of 100 parameters for the MNIST and chairs dataset. For the Euler penalization we take an initial parameter λ = 5 × 10 -5 which is then divided 10 times periodically by a factor 2 during the training.

D Additional details on 2d examples

We provide additional details on some toy examples presented in Section 6.1. First we consider the eight Gaussian test case, study the transformation of a uniform mesh by the NF model and compare with the CPFlow architecture. As shown in Figure 7, the mesh transformation when adding the GP flow is getting closer to the CP-Flow one, and the OT cost is roughly the same. Another illustration of the transformation of the source distribution is given at the bottom of Figure 7. We clearly see the added value of GP flow, as the points' configuration gets much closer to an isotropic distribution indicating that we get closer to the Monge map. Note that adding the GP flow does not affect the estimated density nor the test loss.

Euler's penalization. We now turn our attention to the two moons dataset and evaluate how close our penalization procedure is from the true solution of Euler equations. To this end we compare our solution with a more standard method where the initial condition is optimized through a numerical scheme (implemented in Pytorch to be differentiable) which directly solves the Euler equations, in the line of "differentiable physics" approaches (de Avila [START_REF] De | End-to-end differentiable physics for learning and control[END_REF]. The baseline numerical scheme is a spectral method, a very efficient and popular numerical method [START_REF] Canuto | Spectral Methods[END_REF]. We do not go into details on how to implement such methods since they are restricted in practice to low dimensions and we refer to [START_REF] Canuto | Spectral Methods[END_REF]; [START_REF] Quarteroni | Numerical models for differential problems[END_REF] for a complete presentation. To have a more precise comparison we use 20 layers and a learning rate of 5×10 -3 . On Figure 8, we compare three different cases: GP alone, GP with Euler penalization and GP with Euler constrain through a spectral method. As indicated by the OT costs the final positions of particles are very similar for each model. We observe that the trajectories, the initial incompressible velocities and the energy Ē = 1

N N i=1
T 0 1 2 v(t, x i ) 2 dt are very similar both for the penalization-based procedure and the spectral method which show that we can correctly solve the Euler equations in 2D with our penalization approach (the additional benefit being the generalization to high dimensions). Without the addition of the Euler constraint however, the initial velocity field and the trajectories look very different resulting in a higher energy, and thus a transformation that does not correspond to a geodesic in SDiff(Ω).

E Interpolation examples

3 The 2D test cases can be found for example here https://github.com/rtqichen/ffjord/blob/master/lib/toy_data.py

Figure 1 :

 1 Figure 1: A GP transformation applied on particles sampled from a two-dimensional Gaussian distribution. The mean and standard deviation stay the same; only the positions of the particles change.

Figure 2 :

 2 Figure 2: Eight Gaussians test case with colored distributions. A GP flow is trained on a pre-trained BNAF model[START_REF] De Cao | Block neural autoregressive flow[END_REF] to reduce the OT cost.

Figure 3 :

 3 Figure 3: Color map for several 2d toy examples. Adding GP flow makes the transformation closer to the OT map.

Figure 4 :Figure 6 :

 46 Figure 4: Trajectories for the pinwheel test case from T = 0 (leftmost image) to T = 1 (rightmost image). Top: Gaussian motion without regularization. Bottom: Gaussian motion with Euler regularization. The latter gives much smoother trajectories. Target distribution GP flow without Euler GP flow + Euler Exact discrete OT map

  0 and ψ i j (x) : R d → R are bounded functions. Then the function v defined in Proposition 2 satisfies ∇ • v = 0 and v • n = 0 on ∂Ω.

Figure 9 :

 9 Figure 9: Examples of interpolation for the dSprites dataset where each block correspond to the interpolation of a different data point along the 10 dimensions axis represented by the rows. The dimensions are sorted with respect to their KL divergence in the VAE latent space, so the higher rows carry more information while the last rows should leave the image unchanged. 26

Figure 10 :

 10 Figure 10: Examples of interpolation for the dSprites dataset where each block correspond to the interpolation of a different data point along the 10 dimensions axis represented by the rows. The dimensions are sorted with respect to their KL divergence in the VAE latent space, so the higher rows carry more information while the last rows should leave the image unchanged. 27

Figure 12 :

 12 Figure 12: Examples of interpolation for the MNIST dataset where each block correspond to the interpolation of a different data point along the 10 dimensions axis represented by the rows. The dimensions are sorted with respect to their KL divergence in the VAE latent space, so the higher rows carry more information while the last rows should leave the image unchanged. 29

Figure 13 :

 13 Figure 13: Examples of interpolation for the chairs dataset where each block correspond to the interpolation of a different data point along the 10 dimensions axis represented by the rows. The dimensions are sorted with respect to their KL divergence in the VAE latent space, so the higher rows carry more information while the last rows should leave the image unchanged. 30

Table 1 :

 1 Quantitative evaluation from[START_REF] Eastwood | A framework for the quantitative evaluation of disentangled representations[END_REF] of disentanglement (higher is better), completeness (higher is better) and informativeness (lower is better) on the dSprites dataset. Adding GP flows make the scores closer to the initial ones. For all models considered here (except in the case of the vanilla FFJORD model) the values represent the mean and standard deviation taken over 3 runs.

	Disentenganlement ↑ Completeness ↑ Informativness ↓
	Init. latent space		0.58	0.81	0.55
	FFJORD		0.39	0.26	0.62
	FFJORD+GP	0.58 ± 0.01	0.66 ± 0.01	0.62 ± 0.
	FFJORD+GP+EULER	0.58 ± 0.01	0.66 ± 0.01	0.58 ± 0.
	Model		dSprites	MNIST	chairs
	OT cost
	FFJORD		10.45	6.81	5.98
	FFJORD+GP		5.69 ± 0.03	3.11 ± 0.01	2.34 ± 0.01
	FFJORD+GP+EULER	5.30 ± 0.01	3.11 ± 0.01	2.36 ± 0.01
	CPFlow (3 blocks)		6.27 ± 0.06	6.01 ± 2.12	2.46 ± 0.14
	CPFlow (1 block)		8.97 ± 3.44	27.92 ± 3.00 2.21 ± 0.50
	Loss
	FFJORD		-16.62	-0.45	5.48
	CPFlow (1 block)		-16.99 ± 0.08	0.95 ± 0.11	6.45 ± 0.01

FFJORD+GP

-16.62 ± 0.

-0.45 ± 0. 5.48 ± 0.

FFJORD+GP+EULER

-16.62 ± 0.

-0.45 ± 0. 5.48 ± 0.

CPFlow (3 blocks)

-19.07 ± 0.06 -0.06 ± 0.05 5.94 ± 0.03

Table 2 :

 2 Losses and mean OT costs. GP flows reduce the OT cost without changing the loss. Adding Euler regularization allows to further reduced the OT cost on the dSprites dataset. For all models considered here (except in the case of the vanilla FFJORD model) the values represent the mean and standard deviation taken over 3 runs.

1 Gaussian preserving flows B.1.1 Proof of Proposition 1

  We recall that here erf : R d → R d is the distribution function of a one dimensional Gaussian applied component wise.Proposition.A map s is a smooth Gaussian preserving function satisfying (3) if and only if there exists φ : (-1, 1) d → (-1, 1) d such that | det ∇φ| = 1 and s

B.

Table 3 :

 3 Parameters used for the training of GP flows on the 2D toy examples.

	Model	nb params (pre-trained model) + GP	epochs nb layers batch size	lr	nb decay euler
	eight Gaussian
	BNAF+GP		(15.4K) + 332		2000		20	1000	10 -2	5
	two moons
	BNAF+GP		(15.4K) + 332		1000		15	1000	2 × 10 -3	5
	pinwheel
	BNAF+GP		(15.4K) + 332		800		15	1000	2 × 10 -3	4
	We recall that the notation d/dt must be understood as deriving the first variable of v(t, X(t, x)) while D/Dt
	represents the derivative in time of the function f (t) := v(t, X(t, x)). By deriving the first and second variable
	with respect to t one has
		D Dt	v(t, X(t, x)) =	d dt	v(t, X(t, x)) + (	d dt	X(t, x) • ∇)v(t, X(t, x)).
	Using d

dt X(t, x) = v(t, X(t, x)) one finally obtains

D Dt v(t, X(t, x)) = d dt v(t, X(t, x)) + (v(t, x • ∇)v(t, X(t, x)).

github.com/rtqichen/ffjord, github.com/nicola-decao/BNAF, github.com/CW-Huang/CP-Flow

https://math.stackexchange.com/questions/578898

Acknowledgments

The research leading to these results has received funding from IMT Atlantique, Cominlabs Labex (DY-NALEARN project) and ANR (AI4CHILD, ANR-19-CHIA-0015-01, project LEMONADE, ANR-21-CE48-0005 and project OTTOPIA ANR-20-CHIA-0030).

B.2.3 Proof of Lemma 3

Lemma. Let n ∈ N, n ≤ d-2 and consider the function u n : R d → R d-n . Then the vector field v n : R d → R d defined as

is divergence free.

Proof. As explain in Appendix A the formulation (32) is equivalent to consider v n under the form (7) with ψ i j defined as in ( 16). From Proposition 2 the function v n is divergence free.

B.2.4 Proof of Proposition 3

Proposition. Let v be a divergence free function in R d . Then there exists d -

Proof. Let v(x) be a divergence free function and denote (ψ i j ) its coefficients from (7). We construct the associate vectors u n of the functions v n procedurally: for k = 1, ..., d -1, we iteratively chose u k-1 1 arbitrarily and define the other components of the vector u k-1 as

note that the sum is well-defined because the components u k are constructed iteratively from k = 1 to k = d -1 and we have adopted the convention -1 n=0 = 0. We claim that with this construction we recover the equality (33). Indeed from (34) one has

Using ( 16) one gets

Again using (16) one has (ψ k j ) n = 0 for n ≥ k and therefore

Since the coefficients ψ k j and (ψ k j ) n are all antisymmetric this equality is also satisfied for j < k + 1. We conclude with the decomposition (7) of the divergence free functions.

B.3 Additional material

B.3.1 Total derivative expansion

In this subsection we recall the expansion of the Lagrangian derivative