Ricardo Baptista

Youssef Marzouk

Olivier Zahm
email: olivier.zahm@inria.fr

On the representation and learning of monotone triangular transport maps

Keywords: Transportation of measure, Knothe-Rosenblatt rearrangement, normalizing flows, monotone functions, infinite-dimensional optimization, adaptive approximation, multivariate polynomials, wavelets, density estimation. MSC2020 subject classifications: 65C20, 49Q22, 62G07

Transportation of measure provides a versatile approach for modeling complex probability distributions, with applications in density estimation, Bayesian inference, generative modeling, and beyond. Monotone triangular transport maps-approximations of the Knothe-Rosenblatt (KR) rearrangement-are a canonical choice for these tasks. Yet the representation and parameterization of such maps have a significant impact on their generality and expressiveness, and on properties of the optimization problem that arises in learning a map from data (e.g., via maximum likelihood estimation). We present a general framework for representing monotone triangular maps via invertible transformations of smooth functions. We establish conditions on the transformation such that the associated infinite-dimensional minimization problem has no spurious local minima, i.e., all local minima are global minima; and we show for target distributions satisfying certain tail conditions that the unique global minimizer corresponds to the KR map. Given a sample from the target, we then propose an adaptive algorithm that estimates a sparse semi-parametric approximation of the underlying KR map. We demonstrate how this framework can be applied to joint and conditional density estimation, likelihood-free inference, and structure learning of directed graphical models, with stable generalization performance across a range of sample sizes.

Introduction

Many sampling, estimation, and inference algorithms seek to characterize a somehow intractable or complex probability distribution µ on R d . Transportation of measure provides a useful and versatile approach to this problem. The underlying idea is to construct a coupling of µ with a tractable "reference" distribution ν on R d -for instance, a standard normal. Formally, one jointly constructs a pair of random variables (X, Z) such that X ∼ µ and Z ∼ ν. A special class of couplings is given by deterministic transformations S : R d → R d such that S(X) = Z in distribution, and a transformation S that satisfies this property is called a transport map [START_REF] Villani | Optimal transport: old and new 338[END_REF]. As a result, if the transport map is invertible, one can generate realizations of X by first simulating Z. If µ and ν have densities π and η with respect to a common base measure, one can also explicitly represent the target density π as a transformation of the reference density η. The construction of such transport maps has found numerous applications: density estimation [START_REF] Tabak | A family of nonparametric density estimation algorithms[END_REF][START_REF] Anderes | A general spline representation for nonparametric and semiparametric density estimates using diffeomorphisms[END_REF][START_REF] Dinh | Density estimation using Real NVP[END_REF][START_REF] Huang | Convex Potential Flows: Universal Probability Distributions with Optimal Transport and Convex Optimization[END_REF]; variational Bayesian inference [START_REF] El Moselhy | Bayesian inference with optimal maps[END_REF][START_REF] Rezende | Variational inference with normalizing flows[END_REF][START_REF] Brennan | Greedy inference with structure-exploiting lazy maps[END_REF]; generative modeling of images, video, and other structured objects [START_REF] Oord | Parallel WaveNet: Fast high-fidelity speech synthesis[END_REF][START_REF] Kingma | Glow: Generative flow with invertible 1x1 convolutions[END_REF]; likelihood-free inference [START_REF] Papamakarios | Fast ε-free inference of simulation models with Bayesian conditional density estimation[END_REF][START_REF] Radev | BayesFlow: Learning complex stochastic models with invertible neural networks[END_REF][START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF]; and beyond.

In general, there exist infinitely many transport maps between two absolutely continuous distributions on R d . In this paper, we focus on a specific, canonical choice: triangular transport maps [START_REF] Bogachev | Triangular transformations of measures[END_REF] of the form

S(x) =      S 1 (x 1 ) S 2 (x 1 , x 2 )
. . .

S d (x 1 , . . . , x d )      , (1)
where each component function S k depends only on the k variables x ≤k := (x 1 , . . . , x k ) and is monotone increasing with respect to the last input x k . In particular, each function S k encodes an increasing rearrangement [START_REF] Villani | Optimal transport: old and new 338[END_REF] between ordered marginal conditionals of µ and ν, i.e., µ(dx k |x <k ) and ν(dz k |z <k ). Later we will discuss properties of such triangular maps-known also as Knothe-Rosenblatt (KR) rearrangements [START_REF] Villani | Optimal transport: old and new 338[END_REF][START_REF] Rosenblatt | Remarks on a multivariate transformation[END_REF][START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]-more precisely, but we comment here on their utility. First of all, triangular structure facilitates computational tractability: S is easy to invert and the determinant of its Jacobian (a lower triangular matrix) is easy to evaluate (see Section 2). Triangular maps have thus been used extensively in the density estimation, inference, and generative modeling applications noted above. For example, triangular maps are the building blocks of many normalizing flows, popularized by the machine learning community [START_REF] Kobyzev | Normalizing flows: An introduction and review of current methods[END_REF][START_REF] Papamakarios | Normalizing flows for probabilistic modeling and inference[END_REF]; specifically, autoregressive normalizing flows define triangular maps [START_REF] Jaini | Sum-of-squares polynomial flow[END_REF] via particular structural choices and parameterizations, and compose these maps to produce more expressive transformations. More fundamentally, because triangular maps expose certain conditionals of µ, they are particularly well suited to conditional density estimation and conditional sampling; we will describe this link in Section 2. Triangular maps also inherit sparsity from the conditional independence properties of µ and ν, as detailed in [START_REF] Spantini | Inference via low-dimensional couplings[END_REF].

It is worth noting, of course, that other canonical choices of transport have different attractive features. For instance, optimal transport maps are invariant under relabeling of the coordinates or more general isometries on R d (unlike triangular maps), and have deep links to partial differential equations [START_REF] Villani | Optimal transport: old and new 338[END_REF]. But optimal transport maps are in general more challenging to represent, evaluate, and estimate from data in the continuous setting; moreover, they do not enjoy such a direct link to conditioning or to graphical models.

Many representations and finite-dimensional parameterizations of monotone triangular maps have been proposed in recent years. These include representations based on polynomials [START_REF] Marzouk | Sampling via Measure Transport: An Introduction In Handbook of Uncertainty Quantification 1-41[END_REF][START_REF] Jaini | Sum-of-squares polynomial flow[END_REF], radial basis functions [START_REF] Tabak | A family of nonparametric density estimation algorithms[END_REF], neural networks of varying capacity [START_REF] Dinh | Density estimation using Real NVP[END_REF][START_REF] Kingma | Glow: Generative flow with invertible 1x1 convolutions[END_REF], and tensor decompositions [START_REF] Cui | Deep composition of Tensor-Trains using squared inverse Rosenblatt transports[END_REF][START_REF] Cui | Conditional Deep Inverse Rosenblatt Transports[END_REF]. A core challenge in this setting is to satisfy the monotonicity constraint ∂ x k S k > 0. For instance, one might enforce the monotonicity constraint at a finite collection of points in the support of π [START_REF] Parno | Transport map accelerated Markov chain Monte Carlo[END_REF], but this approach cannot in general guarantee that S is monotone over the entire support of π. Other efforts have sought to enforce monotonicity by construction-via the parameterization of S itself. For example, [START_REF] Papamakarios | Masked autoregressive flow for density estimation[END_REF] employs map components with affine dependence on the last variable, i.e., S k (x ≤k ) = α(x <k ) + exp(β(x <k ))x k , where α and β are neural networks. While S is then guaranteed to be monotone, it can only represent a restricted class of distributions µ. (If ν is Gaussian, then µ can only be a product of Gaussian marginal conditionals.) Such representations therefore cannot consistently approximate the true KR rearrangement for general µ. Recent work [START_REF] Teshima | Couplingbased invertible neural networks are universal diffeomorphism approximators[END_REF] has shown that a composition of such affine maps, interleaved with rotations and permutations, can approximate a general class of distributions, though approximation rates remain unknown. The required rotations or permutations break the overall triangular structure of the transformation, however. Alternatively, to increase the "expressiveness" of a given triangular function, [START_REF] Wehenkel | Unconstrained monotonic neural networks[END_REF][START_REF] Huang | Neural Autoregressive Flows[END_REF][START_REF] Durkan | Neural spline flows[END_REF] have introduced more general parametric representations of the monotone component functions S k . For distributions with analytic densities on bounded domains, a complete approximation theory for the KR map was recently developed in [START_REF] Zech | Sparse approximation of triangular transports. Part II: the infinite dimensional case[END_REF][START_REF] Zech | Sparse Approximation of Triangular Transports. Part I: The Finite-Dimensional Case[END_REF], using polynomial or ReLU neural network representations of range-constrained monotone triangular functions.

Despite these myriad proposals, relatively little attention has been paid to the structure and tractability of the optimization problem involved in learning triangular maps (e.g., in estimating a map given an i.i.d. sample {X i } from µ). Properties of this optimization problem are intimately tied to the representation and parameterization of the associated map. It is desirable to have a flexible and general representation-one that can consistently recover the KR rearrangement for a broad class of distributions (µ, ν)-that at the same time makes optimization tractable. It is also desirable to have adaptivity: a parameterization whose size or complexity can be adapted to properties of the target distribution and the available sample size, for good empirical statistical performance.

This paper directly addresses these desiderata. We do so by developing and analyzing a functional framework for representing and learning triangular maps. Our main contributions are as follows. We propose a general representation of monotone triangular functions, based on a rectification operator R k that transforms sufficiently smooth non-monotone functions f k : R k → R into monotone component functions S k of a triangular map. This operator takes the form

R k (f k )(x ≤k ) = f k (x <k , 0) + x k 0 g ∂ k f k (x <k , t) dt, (2)
where g : R → R >0 is a positive function. We then analyze the infinite-dimensional optimization problem associated with learning maps from data, recast as optimization over functions {f k } d k=1 . Specifically, we establish conditions on the rectification operator and on the target distribution such that the resulting optimization problem is well-posed, smooth, and has no spurious local minima. Under further conditions on the target distribution (essentially that it has Gaussian tails), we show that the optimization problem has a unique global minimizer corresponding to the KR map.

These theoretical results guarantee, in practice, fast and reliable learning of monotone triangular maps given an appropriate function space V k for each f k . The second main contribution of our paper is algorithmic: given a hierarchical basis (e.g., polynomials or wavelets) for each V k , we propose a greedy adaptive procedure to learn parametric representations of f k . The procedure naturally produces map representations that are sparse and interpretable-in particular, it exploits and implicitly discovers conditional independence. We use these learned maps for density estimation, given an i.i.d. sample {X i } from π. Maintaining a strict triangular structure also exposes marginal conditionals of the target density, thus enabling conditional density estimation. Our numerical experiments show that the algorithm provides robust performance at small-to-moderate sample sizes, and constitutes a semi-parametric approach that naturally links map complexity to the size of the data.

The remainder of the paper is organized as follows. Section 2 recalls properties of triangular transport maps and introduces some estimation problems of interest. Our main theoretical contributions are in Section 3, which introduces a framework for representing monotone triangular maps and analyzes the resulting optimization problem. Section 4 then introduces our greedy adaptive procedure for learning maps, and Section 5 contains numerical experiments. Proofs of certain theoretical results are deferred to the appendix.

Triangular transport for density estimation and simulation

Consider the unsupervised learning problem of approximating a target probability density function π defined on R d , given an i.i.d. sample from π. Our goal is to construct a sufficiently smooth and invertible map S : R d → R d such that the pullback density

S ♯ η(x) = η (S(x)) | det ∇S(x)|, (3)
is a good approximation to π, where η is a simple/tractable probability density function on R d . The choice of η is a degree of freedom of the method, and here we take η(x) ∝ exp(-∥x∥ 2 /2); i.e., η is the density of the standard normal distribution on R d , where

∥ • ∥ is the canonical norm of R d .
To ensure invertibility of S, a common practice is to constrain S to be an increasing lower triangular map of the form [START_REF] Ambrogioni | The kernel mixture network: A nonparametric method for conditional density estimation of continuous random variables[END_REF], where each component

S k : R k → R is such that x k → S k (x <k , x k ) is increasing for all x <k = (x 1 , . . . , x k-1 ) ∈ R k-1
. Such a map is easy to invert1 and has a lower triangular Jacobian ∇S(x) so that

| det ∇S(x)| = d k=1 ∂ k S k (x ≤k
) is readily computable. This structure in fact corresponds to the Knothe-Rosenblatt (KR) rearrangement S KR : the increasing lower triangular map satisfying π(x) = S ♯ KR η(x). For a measure µ on R d that is absolutely continuous with respect to a Gaussian measure ν (and hence has a density π with respect to the Lebesgue measure), the KR rearrangement S KR exists and is the unique map of the form (1) that pulls back η to π (or equivalently pushes forward π to η), up to sets of measure zero [START_REF] Bogachev | Triangular transformations of measures[END_REF].

A useful way to measure the difference between π and its approximation S ♯ η is the Kullback-Leibler (KL) divergence D KL (π||S ♯ η) = log(π/S ♯ η)dπ. As explained below, this choice has direct links to maximum likelihood estimation of S. Furthermore, the following inequality shows that convergence in the KL sense D KL (π||S ♯ η) → 0 implies convergence of S towards S KR in the L 2 π sense. This result is a direct consequence of Corollary 3.10 in [START_REF] Bogachev | Triangular transformations of measures[END_REF]; see Appendix A.1 for a proof. Proposition 1. Let S KR be the KR rearrangement pushing forward a distribution with density π on R d to the standard normal distribution on R d , with density η. For any map S : R d → R d as in (1), we have

∥S KR (x) -S(x)∥ 2 dπ(x) ≤ 2D KL (π||S ♯ η). (4)
Since the standard normal density η is a product of its marginal densities, the KL divergence decomposes as

D KL (π||S ♯ η) = d k=1 J k (S k ) -J k (S KR,k ), (5)
where the functionals J 1 , . . . , J d are given by

J k (s) = 1 2 s(x ≤k ) 2 -log |∂ k s(x ≤k )| π(x)dx. (6)
Minimizing the KL divergence (5) over triangular maps S of the form (1) is therefore equivalent to independently minimizing each objective functional J k to find the associated map component S k [START_REF] Marzouk | Sampling via Measure Transport: An Introduction In Handbook of Uncertainty Quantification 1-41[END_REF]. Solution of these optimization problems is thus embarrassingly parallel. This parallel structure was also exploited for Cholesky factorization via KL minimization in [START_REF] Schäfer | Sparse Cholesky Factorization by Kullback-Leibler Minimization[END_REF]. In addition, minimizing each objective s → J k (s) over functions s : R k → R that are strictly increasing in the last variable is a strictly convex optimization problem; see Lemma 10. Given an i.i.d. sample {X i } n i=1 from π, we can replace the expectation in ( 6) by the sample average, which yields the objective

J k (s) = 1 n n i=1 1 2 s(X i ≤k ) 2 -log ∂ k s(X i ≤k ) . (7)
Minimizing [START_REF] Bogachev | Triangular transformations of measures[END_REF] under the constraint ∂ k s(x ≤k ) > 0 produces an estimator S k of S KR,k . The collection of all such component functions defines an estimator S = ( S 1 , . . . , S d ) of the KR rearrangement, along with an estimate of the density π as π(x) := S ♯ η(x). This S is also the maximum likelihood estimator of S KR , i.e.,

S = arg max

{S∈(1), ∂ k S k >0, k=1,...,d} n i=1 log S ♯ η(X i ),
with the optimization being over the space of monotone increasing and triangular maps of the form [START_REF] Ambrogioni | The kernel mixture network: A nonparametric method for conditional density estimation of continuous random variables[END_REF]. This connection between maximum likelihood estimation and minimization of an empirical forward KL divergence is standard. Furthermore, convexity of the optimization objective is preserved when replacing the expectation in (6) by the sample average in [START_REF] Bogachev | Triangular transformations of measures[END_REF].

A core question when estimating maps and densities by minimizing ( 7) is how to parameterize sufficiently expressive monotone map components S k -i.e., maps capable of representing a wide class of distributions π-while ensuring that the optimization problem can be solved efficiently. As explained in the introduction, this question is intimately tied to the monotonicity constraint ∂ k S k (x ≤k ) > 0. For example, choosing a linear parameterization for S k that admits only affine dependence on x k (to easily enforce monotonicity) allows the map component be identified efficiently through the solution of a least-squares problem (see [START_REF] Spantini | Coupling techniques for nonlinear ensemble filtering[END_REF]Appendix A]), but such maps can only capture distributions that factor into a sequence of Gaussian marginal conditionals. On the other hand, a more complex ansatz for S k will often yield a much more difficult optimization problem. Note that with any nonlinear parameterization of S k , the convexity of ( 6) and [START_REF] Bogachev | Triangular transformations of measures[END_REF] with respect to S k does not in general yield convexity in the parameters. As we shall demonstrate later, many nonlinear parameterizations that enforce monotonicity yield optimization problems that may not even be smooth, and that have many local minima. We will address these issues in Section 3.

Conditional density estimation and sampling

{(Y i , W i )} n
i=1 from the joint density π(y, w). Observe that the KR map immediately has the block structure:

S(y, w) = S Y (y) S W (y, w) , (8)
where y → S Y (y) : R m → R m and w → S W (y * , w) : R p → R p are increasing lower triangular maps, the latter for any y * ∈ R m . Recall that the reference density η is the standard normal on R d , and thus is a product of two standard normals, η 1 (z ′ )η 2 (z), where z ′ ∈ R m and z ∈ R p . Using a KR map of the form (8), we can write the marginal density of Y as π(y) = (S Y KR ) ♯ η 1 (y) and, more interestingly, the conditional density of W as π(w|y) = S W KR (y, •) ♯ η 2 (w). Each of the last p components of the KR rearrangement, S W,k KR (y, w ≤k ), 1 ≤ k ≤ p, can be estimated from a sample {(Y i , W i )} n i=1 by minimizing

J k (s) = 1 n n i=1 1 2 s(Y i , W i ≤k ) 2 -log ∂ m+k s(Y i , W i ≤k ) , (9)
under the constraints ∂ m+k s(y, w ≤k ) > 0. This produces an estimator S W of S W KR , which in turn yields an estimator of the conditional density π(w|y) as π(w|y) := S W (y, •) ♯ η 2 (w). This property has been used to perform conditional density estimation (CDE) in [START_REF] Marzouk | Sampling via Measure Transport: An Introduction In Handbook of Uncertainty Quantification 1-41[END_REF][START_REF] Papamakarios | Fast ε-free inference of simulation models with Bayesian conditional density estimation[END_REF][START_REF] Kovachki | Conditional Sampling With Monotone GANs[END_REF][START_REF] Cui | Conditional Deep Inverse Rosenblatt Transports[END_REF].

One important application of CDE is likelihood-free inference, where W represents a parameter to be inferred and Y are data whose conditional density π(y|w) is computationally intractable or unavailable in closed form. This setting arises in many applications: inference in stochastic models, or inference in the presence of high-dimensional nuisance parameters or latent variables (including parameter inference for state-space models). Given a joint sample {(Y i , W i )} n i=1 , we can estimate the map component S W and use it to simulate from the estimated conditional density π(w|y * ) given any realization of the data y * , simply by sampling z i ∼ η 2 and solving the triangular system S W (y * , w i ) = z i for each w i . Note that we can simulate from or evaluate the estimated conditional densities for multiple realizations of the data y * , including values that are not present in the dataset {Y i } n i=1 . Thus, learning a single map S W parameterized by y is said to amortize the cost of conditional sampling over the data.

Representing and learning continuous monotone functions

In this section, we define a general representation for components of monotone triangular maps and present our main theoretical results. The essential idea, as mentioned in the introduction, is to express each component S k as a nonlinear transformation S k = R k (f ) of a smooth function f , where R k (2) is an operator that enforces the monotonicity constraint by construction. We then identify the map component by solving the re-parameterized optimization problem

min f ∈V k L k (f ) = J k (R k (f )), (10)
where V k is a linear space of functions in which we seek f . With this nonlinear transformation, we lose the convexity of the constrained problem min {s : ∂ k s>0} J k (s), but obtain an unconstrained minimization problem instead. It turns out that, with appropriate conditions on R k , the transformed optimization problem retains many desirable properties. In Section 3.1, we discuss the construction of R k and motivate certain critical choices therein. Then, in Section 3.2, we analyze the regularity of the composed objective functional L k for a certain choice of the function space V k . In Section 3.3 we discuss the existence and uniqueness of solutions to [START_REF] Chang | Chernoff-type bounds for the Gaussian error function[END_REF], and describe conditions under which solving [START_REF] Chang | Chernoff-type bounds for the Gaussian error function[END_REF] exactly recovers the KR rearrangement.

The rectification operator

Recalling (2), for any sufficiently smooth f : R k → R we let R k (f ) : R k → R be the function defined by

R k (f )(x ≤k ) = f (x <k , 0) + x k 0 g ∂ k f (x <k , t) dt,
where g : R → R >0 is a positive function. We call the operator R k : f → R k (f ) a rectifier because it transforms any f into a function which is increasing in the k-th variable, i.e., 1 illustrates the application of the rectifier to a nonlinear univariate function f . Here we choose g(•) = log(1 + exp(•)). S is an increasing transport that pushes forward a onedimensional mixture of Gaussians π(x) = 0.5N (x; -1, 1) + 0.5N (x; 1, 1) to the standard Gaussian reference density η.

∂ k R k (f )(x ≤k ) = g(∂ k f (x ≤k )) > 0. As a simple example, functions of the form f (x ≤k ) = α(x <k ) + β(x <k )x k are transformed into R k (f )(x ≤k ) = α(x <k ) + g(β(x <k ))x k . Figure
The choice of the function g in (2) has a crucial impact on properties of the optimization problem [START_REF] Chang | Chernoff-type bounds for the Gaussian error function[END_REF]. One possible choice, proposed in [START_REF] Jaini | Sum-of-squares polynomial flow[END_REF][START_REF] Brennan | Greedy inference with structure-exploiting lazy maps[END_REF], is the square function g(ξ) = ξ 2 . While this choice permits closed-form computation of the integral in (3.1) when f is polynomial, it yields an optimization problem [START_REF] Chang | Chernoff-type bounds for the Gaussian error function[END_REF] which possesses many spurious local minima and is in fact non-smooth; see Figure 2 for an illustration. This can be explained in part by the fact that this g is not bijective.

Instead, one can choose g to be a bijective function from R to R >0 , such as the soft-plus function,

g(ξ) = log(1 + exp(ξ)), (11)
whose inverse is g -1 (ξ) = log(exp(ξ) -1). Another example of a bijective function, considered in [START_REF] Wehenkel | Unconstrained monotonic neural networks[END_REF], is the shifted exponential linear unit (ELU),

g(ξ) = exp(ξ) ξ < 0 ξ + 1 ξ ≥ 0 , ( 12
)
whose inverse is g

-1 (ξ) = ξ -1 if ξ ≥ 1 and g -1 (ξ) = log(ξ) otherwise.
As a consequence of g being bijective, the inverse of the rectifier R -1 k (s) exists for any sufficiently smooth s : R k → R with ∂ k s(x ≤k ) > 0 and can be written as

R -1 k (s)(x ≤k ) = s(x <k , 0) + x k 0 g -1 ∂ k s(x <k , t) dt. ( 13
)
More importantly, the fact that g is invertible yields an objective function L k that is far better behaved than with g(ξ) = ξ 2 ; see the numerical illustration in Figure 2, with both soft-plus and shifted ELU g. With these choices of g, we observe that L k = J k • R k , though non-convex in general, has no local minima and no saddle points. In the next sections, we will analyze these properties of the optimization problem (10) (e.g., smoothness, existence and uniqueness of solutions) and elucidate their precise dependence on g, V k , and π.

Remark 1. In the finite-dimensional setting, it can easily be shown that composing a smooth and (strictly) convex objective function with a C 1 -diffeomorphic map preserves the (unique) global minima of the objective. Such diffeomorphic maps have been explored for accelerating optimization on finite-dimensional manifolds; see, e.g., [START_REF] Lezcano Casado | Trivializations for gradient-based optimization on manifolds[END_REF]. Establishing that the operator R k in (3.1) is C 1diffeomorphic is, however, non-trivial. We will show in Section 3.3 that, under additional assumptions, our reparameterization yields an infinite-dimensional optimization problem where local minima are global minima. where the rectifier R 1 is defined using the soft-plus g (11) (left), the shifted ELU g (12) (middle), or the square function g(ξ) = ξ 2 (right). Here, π(x) = 1/2N (x; -2, 0.5) + 1/2N (x; 2, 2) is a univariate Gaussian mixture, and we use n = 50 to define J 1 with f 1 represented using a linear combination of Hermite functions up to degree 10. The objective is evaluated along line segments that interpolate between random initial maps (t = 0) and critical points resulting from a gradient-based optimization method (t = 1). Observe that with bijective g (left and middle) the algorithm always arrives at the same optimal value, whereas with the square function g (right) the algorithm gets stuck in local minima and rarely attains the optimal value.

Smoothness of the optimization problem

In this section we give sufficient conditions on the function g and the target density π to guarantee that the objective function L k is smooth over an appropriate function space for f . We introduce this function space V k and show continuity properties of the rectifier that hold for functions f ∈ V k , before stating our main result in Proposition 4. We begin by defining the weighted Sobolev space

V k = f : R k → R such that ∥f ∥ 2 V k := |f (x)| 2 + |∂ k f (x)| 2 η ≤k (x)dx < +∞ , (14)
where η ≤k is the standard normal density on R k . By [START_REF] Kufner | How to define reasonably weighted Sobolev spaces[END_REF]Theorem 1.11], this space is complete and thus is a Hilbert space. The space V k has sufficient regularity for ∂ k f to exist, but also to permit the pointwise evaluation f (x <k , 0), as required in the definition (3.1) of the rectifier. This property is formalized by the following proposition, a trace theorem, which shows that for any f ∈ V k , the function x <k → f (x <k , 0) is a function in L 2 η <k , the η <k -weighted space of square integrable functions.

Proposition 2. There exists a constant C T < ∞ such that for any f ∈ V k f (x <k , 0) 2 η <k (x)dx ≤ C T ∥f ∥ 2 V k . (15)
Proof. See Appendix A.3.

Remark 2. Notice that H 1 η ≤k ⊂ V k ⊂ L 2 η ≤k , where L 2 η ≤k = {f : R k → R : ∥f ∥ 2 L 2 η ≤k := f 2 dη ≤k < ∞} and H 1 η ≤k = {f : R k → R : ∥f ∥ 2 H 1 η ≤k := f 2 + ∥∇f ∥ 2 dη ≤k < ∞} are the standard weighted Sobolev spaces.
Given that the standard normal density factorizes as η ≤k (x) = k i=1 η i (x i ), the space V k admits the following tensor product structure

V k = L 2 η 1 ⊗ • • • ⊗ L 2 η k-1 ⊗ H 1 η k , (16)
and the norm ∥ • ∥ V k is a product norm. 2 This tensor product structure will be used later in Section 4 to construct a numerical scheme for approximating the map components S k .

The following proposition shows that, under mild assumptions on g, the rectifier R k is a Lipschitz continuous operator from V k to itself. The proof relies on the Hardy inequality [START_REF] Muckenhoupt | Hardy's inequality with weights[END_REF] and on Proposition 2. Proposition 3. Let g : R → R >0 be a Lipschitz function, i.e., there exists a constant L < ∞ so that

|g(ξ) -g(ξ ′ )| ≤ L|ξ -ξ ′ |, ( 17
)
holds for any ξ, ξ ′ ∈ R. Then R k (f ) ∈ V k for any f ∈ V k , where R k (f ) is defined in (3.1). Furthermore there exists a constant C < ∞ such that ∥R k (f 1 ) -R k (f 2 )∥ V k ≤ C∥f 1 -f 2 ∥ V k . ( 18
)
holds for any

f 1 , f 2 ∈ V k . Proof. See Appendix A.4.
The next result relies on Proposition 3 to show that the objective functional L k in [START_REF] Chang | Chernoff-type bounds for the Gaussian error function[END_REF], seen as a function from V k to R, is well-defined, continuous, and differentiable.

Theorem 4. Let π be a probability density function on R d satisfying π(x) ≤ C π η(x), (19)
for all x ∈ R d and a constant C π < ∞, and let g : R → R ≥0 be an increasing function where

|g(ξ) -g(ξ ′ )| ≤ L|ξ -ξ ′ |, (20)
| log • g(ξ) -log • g(ξ ′ )| ≤ L|ξ -ξ ′ |, ( 21
) 2 That is ∥v1 ⊗ • • • ⊗ v k ∥V k = ∥v1∥ L 2 η 1 ∥v2∥ L 2 η 2 • • • ∥v k-1 ∥ L 2 η k-1 ∥v k ∥ H 1 η k for any vj ∈ L 2 η k and v k ∈ H 1 η k .
holds for any ξ, ξ ′ ∈ R and a constant L < ∞. Then [START_REF] Bishop | Mixture density networks Technical Report[END_REF] and where R k (f ) is defined as in (3.1). Moreover, the function

L k (f ) := J k (R k (f )) < ∞, for any f ∈ V k , where J k (s) = 1 2 s(x ≤k ) 2 -log |∂ k s(x ≤k )| π(x)dx as in
L k : V k → R is continuous and differentiable at any f ∈ V k and its gradient ∇L k (f ) ∈ V k is given by ⟨∇L k (f ), ε⟩ V k = R k (f )(x) ε(x <k , 0) + x k 0 g ′ ∂ k f (x <k , t) ∂ k ε(x <k , t)dt π(x)dx -(log • g) ′ (∂ k f (x))∂ k ε(x) π(x)dx. ( 22
)
for any ε ∈ V k , where ⟨•, •⟩ V k is the scalar product in V k . Proof. See Appendix A.5.
It is useful to discuss some implications of Theorem 4. A smooth objective function permits us to use deterministic first or second-order optimization algorithms. Furthermore, in Section 4, we exploit the gradients of J k in [START_REF] Jaini | Sum-of-squares polynomial flow[END_REF] to construct an adaptive polynomial (or wavelet) basis for V k . We also note that many functions g satisfy the conditions ( 20) and ( 21) in the proposition. For example, these include the soft-plus [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF] and shifted ELU [START_REF] Cohen | Numerical analysis of wavelet methods[END_REF] functions considered in Figure 2.

Remark 3 (Assumption [START_REF] El Moselhy | Bayesian inference with optimal maps[END_REF] implies Gaussian tails). The condition [START_REF] El Moselhy | Bayesian inference with optimal maps[END_REF] implies that P X∼π (∥X∥ 2 > t) ≤ C π P Z∼η (∥Z∥ 2 > t) holds for any t ≥ 0, and hence that π has sub-Gaussian tails [START_REF] Vershynin | High-dimensional probability: An introduction with applications in data science 47[END_REF]. The reverse, however, does not hold in general. For instance, a compactly supported random variable with density π(x) ∝ 1/ |x|1 {x∈[-1,1]} is sub-Gaussian, but the density is unbounded at x = 0 and thus does not satisfy [START_REF] El Moselhy | Bayesian inference with optimal maps[END_REF]. We also note that (19) can be relaxed with the (less interpretable) assumption that there exist a constant C π < ∞ and diagonal scalings

d k > 0 so that π(x) ≤ C π d k=1 η k (d k x k ) for all x ∈ R d .
In practice, we can always normalize the marginals of π to match the mean and variance of a standard Gaussian reference. Thus, without loss of generality, we will use the assumption above for the remainder of this article.

Remark 4. Under additional assumptions on g, the gradient

f → ∇L k (f ) is locally Lipschitz from V k to V k , where V k = {f ∈ V k , ∂ k f ∈ L ∞ } is the space endowed with the norm ∥f ∥ V k = ∥f ∥ V k +∥∂ k f ∥ L ∞ . More specifically, if g ′ and (log • g) ′ are Lipschitz functions, then there exists a constant M < ∞ such that ∥∇L k (f 1 ) -∇L k (f 2 )∥ V k ≤ M (1 + ∥f 2 ∥ V k )∥f 1 -f 2 ∥ V k , (23)
holds for any f 1 , f 2 ∈ V k . See the derivation of this result in Appendix A.6. Such local Lipschitz regularity is useful to analyze the convergence of backtracking gradient descent procedures, i.e., gradient descent with an inexact line search such as Armijo's rule, to a stationary point; see Theorem 2.1 in [START_REF] Truong | Backtracking Gradient Descent Method and Some Applications in Large Scale Optimisation. Part 2: Algorithms and Experiments[END_REF] and Proposition 2.1.1 in [START_REF] Bertsekas | Nonlinear programming[END_REF]. We leave the extension of such optimization guarantees for solving min f ∈V k L k (f ) to future work.

Existence and uniqueness of solutions

In this section, we show that the optimization problem [START_REF] Chang | Chernoff-type bounds for the Gaussian error function[END_REF] does not admit any spurious local minima, meaning that local minimizers are in fact global minimizers. We also show that problem (10) admits a unique global minimizer which permits us to recover the KR rearrangement. To prove these results, we will need the following proposition, which provides conditions ensuring that the inverse rectifier R -1 k is continuous.

Proposition 5. Let g be a bijective function from R to R >0 such that for any c > 0 there exists a

constant L c < ∞ so that |g -1 (ξ) -g -1 (ξ ′ )| ≤ L c |ξ -ξ ′ |, (24)
holds for any ξ, ξ ′ ≥ c. Then, for any

s ∈ V k such that ess inf ∂ k s > 0, we have R -1 k (s) ∈ V k and ess inf ∂ k R -1 k (s) > -∞. Furthermore for any c > 0, there exists a constant C c < ∞ such that ∥R -1 k (s 1 ) -R -1 k (s 2 )∥ V k ≤ C c ∥s 1 -s 2 ∥ V k , (25)
holds for any s 1 , s 2 ∈ V k such that ess inf ∂ k s i ≥ c.
Proof. See Appendix A.8.

Note that the softplus function [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF] and the shifted exponential linear unit ( 12) satisfy ( 24) with L c = (1 -e -c ) -1 and L c = max{1/c, 1}, respectively.

Local minima are global minima

The following proposition shows that, under certain conditions on the function g, the image set

R k (V k ) = {R k (f ) : f ∈ V k } is convex. Proposition 6. Let g : R → R ≥0 be an increasing function such that ξ → (g -1 (ξ)) 2 is a convex function. Let R k (f ) be defined as in (3.1) and V k as in (14). Then R k (V k ) = {R k (f ) : f ∈ V k } is convex. Proof. See Appendix A.7.
When g is the softplus function [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF] or the exponential linear unit [START_REF] Cohen | Numerical analysis of wavelet methods[END_REF], the mapping ξ → (g -1 (ξ)) 2 is indeed convex. Interestingly, however, the exponential function g(ξ) = exp(ξ)-which has been used to parameterize monotone maps in [START_REF] Marzouk | Sampling via Measure Transport: An Introduction In Handbook of Uncertainty Quantification 1-41[END_REF][START_REF] Spantini | Inference via low-dimensional couplings[END_REF] and for monotone regression in [START_REF] Ramsay | Estimating smooth monotone functions[END_REF][START_REF] Shin | Joint estimation of monotone curves via functional principal component analysis[END_REF]-does not satisfy this property, and there is no guarantee for R k (V k ) to be convex in this case.

An important consequence of Proposition 6 is that the constrained problem min s∈R k (V k ) L k (s) remains convex. Hence, this optimization problem has a unique global minimizer by the strict convexity of L k (see Appendix A.2). The next result shows that the local minima of the unconstrained problem [START_REF] Chang | Chernoff-type bounds for the Gaussian error function[END_REF] are in fact global minima. Theorem 7. Under the assumptions of Propositions 5 and 6, let f * ∈ V k be a local minimizer of f → J k (R k (f )), meaning that there exists ρ > 0 such that

J k (R k (f * )) ≤ J k (R k (f )), ∀f ∈ V k such that ∥f -f * ∥ V k ≤ ρ. ( 26
)
If ess inf ∂ k f * > -∞, then f * is a global minimizer, meaning that J k (R k (f * )) ≤ J k (R k (f )), ∀f ∈ V k . Proof. Let f ∈ V k . For any 0 < t < 1 we let s t = tR k (f ) + (1 -t)R k (f * ). By Proposition 6, R k (V k ) is convex so that s t ∈ R k (V k ). By convexity of J k (c.f. Appendix A.2), we can write J k (s t ) ≤ tJ k (R k (f )) + (1 -t)J k (R k (f * )), or equivalently J k (s t ) -J k (R k (f * )) ≤ t J k (R k (f )) -J k (R k (f * )) . (27)
Next we show that there exists a sufficiently small t > 0 such that the above left-hand side is positive. As a consequence, the right hand side of (27) will be positive, which will conclude the proof.

Let b = ess inf ∂ k f * and c = g(b)/2 > 0. For any t ≤ 1/2, as have ∂ k s t = tg(∂ k f )+(1-t)g(∂ k f * ) ≥ 1/2g(∂ k f * ) so that ess inf ∂ k s t ≥ c. In addition, we have ess inf ∂ k R k (f * ) = ess inf g(∂ k f * ) ≥ c. Thus, by Proposition 5, there exists a constant C c < ∞ such that ∥R -1 k (s t ) -f * ∥ V k = ∥R -1 k (s t ) -R -1 k (R k (f * ))∥ V k ≤ C c ∥s t -R k (f * )∥ V k = tC c ∥R k (f ) -R k (f * )∥ V k . By letting t = ρ/(C c ∥R k (f ) -R k (f * )∥ V k ), we have ∥R -1 k (s t ) -f * ∥ V k ≤ ρ. Therefore, setting f = R -1 k (s t ) in (26) ensures that 0 ≤ J k (s t ) -J k (R k (f * )).
Let us remark that Theorem 7 holds for general target densities π without the bound assumed in Theorem 4. Applying this result, however, depends on the existence of a local minimizer in the function space V k , which must be established for a given target density. The following section studies a relatively broad class of distributions where the k-th component of the KR rearrangement is in R k (V k ), thereby providing a candidate solution to apply Theorem 7.

Uniqueness of the global minimizer and recovery of the KR rearrangement

As discussed in Section 2, the Knothe-Rosenblatt rearrangement S KR is the unique lower triangular and monotone map such that D KL (π||S ♯ KR η) = 0; see [START_REF] Bogachev | Triangular transformations of measures[END_REF]. The decomposition (5) of the KL divergence D KL (π||S ♯ η) thus permits us to write

J k (S KR,k ) ≤ J k (R k (f )), (28)
for any f ∈ V k and for any 1 ≤ k ≤ d. Indeed, by letting S be the map such that S k = R k (f ) and S i = S KR,i for i ̸ = k, (5) yields [START_REF] Kufner | How to define reasonably weighted Sobolev spaces[END_REF]. Thus, if there exists a function

f KR,k ∈ V k such that S KR,k = R k (f KR,k ), then f KR,k is a global minimizer of f → J k (R k (f )) over f ∈ V k .
To show this, we first need the following intermediate result.

Proposition 8. Let π be a probability density function on R d such that

cη(x) ≤ π(x) ≤ Cη(x), (29)
for all x ∈ R d with some constants 0 < c ≤ C < ∞. Then, for all

x <k ∈ R k-1 and k = 1, . . . , d, S KR,k (x <k , x k ) = O(x k ) and ∂ k S KR,k (x <k , x k ) = O(1) as |x k | → ∞. Furthermore, we have S KR,k ∈ V k and ess inf ∂ k S KR,k > 0 for all k = 1, . . . , d.
Proof. See Appendix A.9.

Let us comment on the assumption [START_REF] Lezcano Casado | Trivializations for gradient-based optimization on manifolds[END_REF]. While the results in Section 3.2 assume only an upper bound on the joint density π, here we have both an upper and a lower bound. These bounds together imply that the k-th marginal conditional density, for

1 ≤ k ≤ d, satisfies (c/C)η k (x k ) ≤ π(x k |x <k ) ≤ (C/c)η k (x k ) for all x ≤k ∈ R k .
This means that the marginal conditionals of the target density have Gaussian tails, rather than potentially being lighter than Gaussian (as in assumption [START_REF] El Moselhy | Bayesian inference with optimal maps[END_REF]). This condition guarantees that the asymptotic behavior of each map component S KR,k in its last variable (as |x k | → ∞) is affine.

We now combine our previous results to show that under the assumption (29) on the target density π, we have f KR,k ∈ V k and the existence of a unique solution to

min f ∈V k J k (R k (f )).
Corollary 9. Under the assumptions on π in Proposition 8, let g be a Lipschitz bijection from R to R ≥0 such that ξ → (g -1 (ξ)) 2 is convex and such that, for any c > 0 there exists a constant

L c < ∞ such that |g -1 (ξ) -g -1 (ξ ′ )| ≤ L c |ξ -ξ ′ | for any ξ, ξ ′ ≥ c. Then, for any 1 ≤ k ≤ d, there exists a unique function f KR,k ∈ V k that satisfies ess inf ∂ k f KR,k > -∞ such that R k (f KR,k ) is the k-th component of the KR rearrangement S KR . As a consequence, min f ∈V k J k (R k (f )) admits a unique solution and J k (R k (f KR,k )) = min f ∈V k J k (R k (f )).
Proof. Combining Proposition 8 and Proposition 5, there exists a function

f KR,k ∈ V k that satisfies ess inf ∂ k f KR,k > -∞ such that R k (f KR,k ) = S KR,k .
Then from Theorem 7 and inequality [START_REF] Kufner | How to define reasonably weighted Sobolev spaces[END_REF],

we have that f KR,k is a global minimizer of f → J k (R k (f )) over f ∈ V k .
From the uniqueness of the KR rearrangement S KR,k among monotone triangular maps and the injectivity of

R -1 k , f KR,k = R -1 k (S KR,k ) is unique.

Adaptive parameterization of transport maps

Given an i.i.d. sample {X i } n i=1 ∼ π, we now propose an adaptive algorithm to build f ∈ V k that minimizes the empirical objective

L k := J k • R k (f )
, where J k is as in [START_REF] Bogachev | Triangular transformations of measures[END_REF]. Doing so for each k = 1, . . . , d yields a monotone triangular map, as described in Section 2. For each k, we represent

f ∈ V k with an m-term expansion f (x ≤k ) = α∈Λ c α ψ α (x ≤k ), (30)
where Λ is a set of multi-indices α = (α 1 , . . . , α k ) ∈ N k 0 with #Λ = m, c α ∈ R are coefficients, and ψ α : R k → R are basis functions for V k , constructed as products of univariate functions,

ψ α (x ≤k ) = k j=1 ψ j α j (x j ).
Here, {ψ j α } α∈N 0 is chosen to be a basis of L 2 η j if j < k and of H 1 η k if j = k. Because V k possesses a tensor product structure [START_REF] Dinh | Density estimation using Real NVP[END_REF] and because its norm ∥ • ∥ V k is a product norm, the basis {ψ α } α∈N k 0 is orthonormal if {ψ j α } α∈N 0 is an orthonormal basis for all j ≤ k. Since f depends linearly on the coefficients c α , the smoothness properties of the objective L k transfer to the objective function parameterized by the coefficients c α , α ∈ Λ.

Sections 4.1 and 4.2 describe two different choices of basis: polynomials and wavelets, respectively. Then, in Section 4.3, we present a greedy algorithm for building the multi-index set Λ, applicable to any such hierarchical basis. In addition, we propose a cross validation procedure to determine when to stop enriching Λ in order to avoid over-fitting to the finite sample.

Polynomial basis

Probabilists' Hermite polynomials {φ α } α∈N form an orthogonal basis for L 2 η where η is the univariate standard Gaussian density. φ α is a polynomial of degree α ≥ 0 defined as

φ α (x) = (-1) α exp(x 2 /2) d α dx α exp(-x 2 /2).
Furthermore, we have ⟨φ α , φ β ⟩ L 2 η = α! δ α,β so that {φ α / √ α!} α∈N 0 forms an orthonormal basis for L 2 η . Similarly, ⟨φ α , φ β ⟩ H 1 η = (α + 1)! δ α,β so that {φ α / (α + 1)!} α∈N 0 forms a orthonormal basis for H 1 η ; see [START_REF] Schmuland | Dirichlet forms with polynomial domain[END_REF]Proposition 1.3]. We can thus tensorize these basis functions to obtain an orthonormal basis for the anisotropic space V k [START_REF] Dinh | Density estimation using Real NVP[END_REF].

In practice, we modify the univariate Hermite polynomials to be linear outside of an arbitrary compact domain. Following [START_REF] Parno | Transport map accelerated Markov chain Monte Carlo[END_REF], we let

ψ α j (x j ) = 1 Z α j    φ α j (x a j ) + φ ′ α j (x a j )(x j -x a j ) if x j < x a j φ α j (x j ) if x a j ≤ x ≤ x b j φ α j (x b j ) + φ ′ α j (x b j )(x j -x b j ) if x j > x b j ( 31
)
where Z α j = α j ! for j = 1, . . . , k -1 and Z α j = (α j + 1)! for j = k. In our numerical experiments, we set x a j and x b j to be the 0.01 and 0.99 (empirical) quantiles, respectively, of the marginal sample {X i j } n i=1 . The basis {ψ α j } α j ∈N 0 , albeit not exactly orthonormal, is close to being orthonormal for sufficiently small x a j and large x b j . This is important for the numerical stability of the algorithm presented in the following section.

The basis {ψ α } α generated by the piecewise polynomials (31) has numerical and structural advantages over the Hermite polynomials. First, these functions avoid the fast growth rate of standard polynomials as ∥x∥ → ∞, providing a more stable extrapolation of the map to the tails of π, where there are typically no samples available to otherwise constrain its growth. Second, a map that reverts to linearity far from the origin yields densities in the class considered in Subsection 3.3.2; see condition [START_REF] Lezcano Casado | Trivializations for gradient-based optimization on manifolds[END_REF]. Indeed, each map component

S k = R k (f ) with f ∈ span{ψ α : α ∈ Λ} behaves linearly as |x k | → ∞,
so that the pullback density S ♯ η has Gaussian tails like the reference η.

Wavelet basis

Wavelets are a popular approach for approximating functions that are not uniformly smooth [START_REF] Cohen | Numerical analysis of wavelet methods[END_REF][START_REF] Vidakovic | Statistical modeling by wavelets 503[END_REF]. In particular, wavelet techniques define a multi-resolution approximation to a function that can better capture local features than, for instance, global polynomials. Given a compactly supported function ψ : R → R called the mother wavelet, and indices (l, q) ∈ Z 2 , the function

ψ (l,q) : x → 2 l/2 ψ(2 l x -q),
is the qth wavelet at the level l. Common choices for the mother wavelet include the Haar, the Daubechies, and the Meyer wavelets.

In our numerical experiments, we use the continuous Mexican hat mother wavelet

ψ(x) = 2 √ 3σπ 1/4 (1 -(x/σ) 2 ) exp(-x 2 /(2σ 2 )), (32)
with scale parameter σ = 1 [START_REF] Mallat | A wavelet tour of signal processing[END_REF], which is obtained from the second derivative of the univariate Gaussian density. The Mexican hat function is also commonly referred to as the Ricker wavelet in the geophysics community. We treat this function as having essentially compact support on [-6, 6], up to numerical precision. Tensorizing a univariate wavelet in the kth coordinate direction with modified Hermite polynomials in the first k-1 directions yields an expansion of the form in [START_REF] Lichman | UCI Machine Learning Repository[END_REF], where the indices in α = (α 1 , . . . , α k ) are tuples of the form

α 1 , . . . , α k-1 ∈ N 0 , α k = (l k , q k ) ∈ Z 2 . ( 33
)
In practice, we consider a truncated set of wavelets by first rescaling the mother wavelet ψ to have the same support as the interval between the 0.01 and 0.99 empirical quantiles of the marginal sample {X i j } n i=1 . Then, we constrain l j ≥ 0 and consider l j = 0 to be the coarsest level of the approximation. Next, for any l j , the translation parameter q j is bounded as 0 ≤ q j ≤ 2 l j -1. Thus, we can consider

α 1 , . . . , α k-1 ∈ N 0 , α k = (l k , q k ) ∈ N 2 0 .
rather than [START_REF] Marzouk | Sampling via Measure Transport: An Introduction In Handbook of Uncertainty Quantification 1-41[END_REF]. Lastly, we note that any function f of the form (30) expanded with compactly supported wavelets will decay to 0 as |x k | → ∞. Thus, the map S k = R k (f ) will have asymptotic linear growth as |x k | → ∞, similarly to the modified Hermite polynomial expansions.

Adaptive transport map algorithm

Now we propose a greedy method for constructing the multi-index set Λ = Λ t in [START_REF] Lichman | UCI Machine Learning Repository[END_REF]. For simplicity, we consider only the case of single indices α j ∈ N 0 . The extension to the case α j ∈ N 2 0 (as in the wavelet construction above) is described in Appendix B. At each greedy iteration t, we add one multi-index α * t to Λ t . Starting with Λ 0 = ∅ we let

Λ t+1 = Λ t ∪ {α * t }, (34)
where α * t / ∈ Λ t is a multi-index that yields the best improvement of L k in a sense to be defined below. Borrowing ideas from [START_REF] Migliorati | Adaptive polynomial approximation by means of random discrete least squares[END_REF][START_REF] Migliorati | Adaptive approximation by optimal weighted least-squares methods[END_REF], we constrain the sets {Λ t } t≥0 to be downward closed [START_REF] Chkifa | Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs[END_REF][START_REF] Cohen | Multivariate approximation in downward closed polynomial spaces[END_REF], meaning that they satisfy the property

α ∈ Λ t and α ′ ≤ α ⇒ α ′ ∈ Λ t , (35)
where α ′ ≤ α denotes α ′ i ≤ α i for all 1 ≤ i ≤ k. Intuitively, [START_REF] Migliorati | Adaptive approximation by optimal weighted least-squares methods[END_REF] means that Λ t has a pyramidal shape containing no holes. Downward-closed sets are known to preserve good approximation properties and permit a tractable construction of Λ t ; see [START_REF] Cohen | Multivariate approximation in downward closed polynomial spaces[END_REF]. Indeed, Λ t+1 remains downward-closed if and only if the multi-index α * t is selected from the so-called reduced margin of Λ t , defined by

Λ RM t = {α / ∈ Λ t : α -e i ∈ Λ t for all 1 ≤ i ≤ k with α i ̸ = 0}
, where e i denotes the i-th canonical vector of N k . The reduced margin is a subset of the margin set Λ M t (i.e., multi-indices α ̸ ∈ Λ t such that ∃ i > 0 where αe i ∈ Λ t ); see Figure 3. The size of the reduced margin typically grows more slowly with respect to the dimension k than the margin itself. For instance, if Λ t contains all multi-indices in a hypercube of width p in dimension k, the margin has cardinality (p + 1) k -p k , while the reduced margin has cardinality k.

Denoting by

f t the minimizer of f → L k (f ) over f ∈ span{ψ α , α ∈ Λ t }, we select α * t in the reduced margin Λ RM t with the following heuristic α * t ∈ arg max α∈Λ RM t |∇ α L k (f t )|. (36)
Here, the notation ∇ α L k (f t ) denotes the derivative of c α → L k (f t + c α ψ α ) evaluated at c α = 0. In other words, we select the multi-index α * t by choosing the largest functional derivative of L k at f t with respect to the basis functions not contained in the current expansion. This greedy procedure for learning each map component is presented in detail in Algorithm 1. Update the active set:

Λt+1 = Λt ∪ {α * t } 7:
Update the approximation: ft+1 = arg min f ∈span{ψα:α∈Λ t+1 } L k (f ) 8: end for 9: Output:

S k = R k (fm)
Remark 5. If the objective function L k is the linear least-squares loss, then Algorithm 1 corresponds to orthogonal matching pursuit for sparse linear regression with normalized basis functions. An alternative heuristic for selecting α * t that requires second-order information from the objective is

α * t ∈ arg max α∈Λ RM t |∇ α L(f t )| 2 /|∇ 2 α L(f t )|.
We found this criterion to perform similarly to (36) for the polynomial basis, but it led to faster convergence of the objective when working with the wavelet basis. Remark 6. A potential drawback of Algorithm 1 is that the greedy enrichment procedure cannot "see" behind the reduced margin. For instance, if a relevant multi-index is located far beyond the reduced margin and the gradient ∇ α L k (f t ) vanishes for all α ∈ Λ RM t , then the algorithm will be stuck. [START_REF] Migliorati | Adaptive approximation by optimal weighted least-squares methods[END_REF] suggested a safeguard mechanism to avoid this behaviour: arbitrarily activate the most ancient index from the reduced margin every t sg ∈ N iterations. This modification, however, was not needed in our numerical experiments.

Remark 7. As pointed out in [START_REF] Cohen | Multivariate approximation in downward closed polynomial spaces[END_REF][START_REF] Migliorati | Adaptive polynomial approximation by means of random discrete least squares[END_REF][START_REF] Migliorati | Adaptive approximation by optimal weighted least-squares methods[END_REF], adding multiple multi-indices at each greedy iteration could also yield better performance compared to adding only one multi-index at a time. The so-called "bulk-chasing" procedure identifies a subset λ * t ⊂ Λ RM t of multi-indices that capture a fraction of the L 2 -norm of the gradient ∇ α L k (f t ) along the reduced margin.

To determine the maximal cardinality of Λ t (i.e., when to stop adaptation) we use ν-fold crossvalidation as in [START_REF] Wasserman | All of statistics: a concise course in statistical inference[END_REF][START_REF] Bigoni | Nonlinear dimension reduction for surrogate modeling using gradient information[END_REF]. For each fold, we run Algorithm 1 with ν -1 folds of the training data for m = n ν-1 ν iterations, and evaluate the objective function in [START_REF] Bogachev | Triangular transformations of measures[END_REF] on the held-out test set. We then select the number of terms m * ≤ n in the expansion (30) that minimizes the test error averaged over the ν folds. The full sample is then used to run Algorithm 1 for m * iterations. The complete procedure for learning each map component S k , with cross-validation, is called the Adaptive Transport Map (ATM) algorithm. This procedure is detailed in Algorithm 2. In practice, we also stop the training on each fold early if the log-likelihood of the test set does not continue to increase for more than 20 iterations.

The cross-validation procedure produces an approximation for f with an expansion whose cardinality is adapted to the size of the training sample n. With a larger sample, we observe that the cross-validation procedure reliably adds more parameters, when needed, to reduce the bias in the approximation while controlling the variance of the estimated map parameters. Thus, we consider ATM a semi-parametric approach for approximating monotone triangular transport maps. Adaptation of the map complexity to the sample size n will be demonstrated in the numerical results to follow. Evaluate log-likelihood log( S j,t k ) ♯ η for iterations t = 1, . . . , m with sample χ test j 8: end for 9: Define m * ∈ {1, . . . , m} as the minimizer of negative log-likelihood t → ν j=1 -log( S j,t k ) ♯ η(χ test j ) 10: Estimate map f * using Algorithm 1 for m * iterations with sample χ 11:

Output: S k = R k (f * )

Numerical experiments

In this section, we evaluate the performance of the ATM algorithm for learning monotone triangular transport maps associated with a variety of target distributions. Sections 5.1 and 5.2 evaluate and visualize the approximation power of the proposed method for strongly non-Gaussian targets. Sections 5.3 illustrates the benefits of adaptivity over a range of sample sizes, and Section 5.4 demonstrates how the estimated maps reveal and exploit conditional independence structure in the target distribution. Section 5.5 presents joint and conditional density estimation results for a suite of UCI datasets. Code for reproducing all of these numerical results is available online. 3In all of the numerical experiments, we pre-process the data by standardizing each marginal, i.e., subtracting the empirical mean and dividing each component of X by its empirical standard deviation. Within the rectifier R k , we employ the modified soft-plus function g(ξ) = log(1 + 2 ξ )/ log(2), which satisfies the conditions of Propositions 4, 5, and 6. This function also satisfies g(0) = 1, so that f (x

1:k ) = 0 is transformed into R k (f )(x 1:k ) = x k .
We evaluate the integral in (3.1) numerically using an adaptive quadrature method based on Clenshaw-Curtis points with a relative error of 10 -3 . At each iteration of the ATM algorithm, we optimize over the coefficients c α for α ∈ Λ t using a BFGS (quasi-Newton) method [START_REF] Nocedal | Numerical optimization[END_REF].

One-dimensional examples

We first consider a one-dimensional mixture of Gaussians with density π(x) = 0.5N (x; -2, 0.5) + 0.5N (x; 2, 2). To estimate π, we generate n = 10 4 training samples {X i } n i=1 and use them to estimate the map S KR that pushes forward π to η. In one dimension, the KR rearrangement (which coincides with the optimal transport map [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]) is simply the increasing function S KR = F -1 η • F π , where F η and F π denote the cumulative distribution functions of the reference and target distributions, respectively. We can thus compare our estimated maps to this exact expression. We approximate S as the transformation R(f ) of a smooth non-monotone function f , seeking f in the finite-dimensional space V p 1 ⊂ V 1 spanned by polynomials of degree at most p. Figure 4a plots the estimated pullback densities π = S ♯ η for increasing polynomial degree p, or equivalently an increasing number of terms #Λ in f (30). Figure 4b shows convergence of this approximation in KL divergence, again for increasing p. We also observe convergence of the estimated map S towards the true KR rearrangement in the L 2 π norm, as expected from the upper bound in (4). Figure 5 illustrates the approximation to the map S and the associated non-monotone function f for different choices of basis {ψ α } α . Here we observe that approximating f using the modified Hermite polynomials [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF] very closely tracks the KR rearrangement S KR and the function f KR = R -1 (S KR ) as x → ±∞, as compared to standard Hermite polynomials or Hermite functions [START_REF] Boyd | Asymptotic coefficients of Hermite function series[END_REF]. Enforcing linear asymptotic behavior of the transport map ensures that the approximate distribution has Gaussian tails in regions where there are few samples. Next we consider a Gaussian mixture with density π(x) = 0.5N (x; 0, 1) + 0.5N (x; 0, 0.025).This mixture is a common test case for sampling and density estimation methods, as it evaluates their ability to capture densities with multiple scales [START_REF] Sisson | Sequential Monte Carlo without likelihoods[END_REF]. Given n = 10 4 samples, we compute the approximate map using either modified Hermite polynomials (Section 4.1) or Mexican hat wavelets (Section 4.2). Figure 6(a) plots approximations to the target density using an m = 15 term expansion (30) for f , while Figure 6(b) shows convergence in KL divergence for both the polynomial and wavelet bases. We observe that the polynomial approximation suffers from oscillations, while the wavelets better capture localized features. This results in a much faster convergence of the KL divergence for wavelets than for polynomials.

Two-dimensional datasets

To demonstrate the expressiveness of the modified Hermite polynomial basis [START_REF] Lueckmann | Benchmarking simulation-based inference[END_REF], we use the ATM algorithm to approximate the Knothe-Rosenblatt rearrangement for several two-dimensional distributions with strongly non-Gaussian geometries: the "banana," "funnel," "cosine," "mixture of Gaussians (MoG)," and "ring" distributions π considered in [START_REF] Wenliang | Learning deep kernels for exponential family densities[END_REF][START_REF] Jaini | Sum-of-squares polynomial flow[END_REF]. For each distribution, we generate an i.i.d. sample of size n = 10 4 from π and apply a random rotation to the data using an angle uniformly distributed in [0, π/2]. Figure 7 plots the true densities π and the approximate densities π := S ♯ η found using ATM. We use 5-fold cross-validation to identify the optimal number of elements in the multi-index set for each map component: S 1 and S 2 . Figure 7 indicates the total number of coefficients #Λ m * across both map components, underscoring the parsimony of the approximation. Five different densities π (top row) and their approximations π (bottom row) built using the ATM algorithm with a Hermite polynomial basis for V k , given a sample of size n = 10 4 from π. The total number of coefficients #Λ m * in both map components is indicated below each density.

Truth Banana Funnel Cosine MoG Ring ATM #Λ m * = 44 #Λ m * = 31 #Λ m * = 72 #Λ m * = 10 #Λ m * = 39

Random mixture of Gaussians

In this example, we evaluate the stability and accuracy of the ATM algorithm for learning transport maps in the small sample regime. We consider a three-dimensional distribution π defined as a mixture of Gaussians centered at the eight vertices of the hypercube [-4, 4] 3 , each with covariance matrix I d . The weights of the mixture components are randomly sampled from a uniform distribution U([0, 1]) and then normalized so that R 3 dπ = 1. To estimate π, we generate a training sample {X i } n i=1 ∼ π and use the ATM algorithm to estimate the KR rearrangement that pushes forward π to η using 5-fold cross validation.

Figure 8a plots the KL divergence of π from π averaged over 10 experiments. Here, the training sets used to build π = S ♯ η are of varying size n ∈ [10 1 , 10 4 ] and the reported KL divergence is computed on a common test set of 10 4 samples. Figure 8b plots the total number of coefficients #Λ m * identified by the ATM algorithm. The performance of ATM is compared to a non-adaptive method where Λ = Λ(p) := {α ∈ N k 0 , ∥α∥ 1 ≤ p} is arbitrarily fixed, for p = 1, 3, 5. Note that Λ(p) corresponds to polynomial f with total degree p. For each sample size n, ATM consistently finds a better estimator of π (in the sense of KL divergence) than a non-adaptive method with a fixed degree p, by adaptively identifying the basis {ψ α } α∈Λt to represent the map components. In addition, the ATM estimator achieves smaller KL divergence with fewer total map coefficients, as seen in Figure 8b.

Stochastic volatility

Next we consider data from a Markov process that describes the volatility of the return on a financial asset over time. The model has two hyperparameters µ and ϕ, with a state (Z k ) k that represents the log-volatility at times k = 1, . . . , T . The two hyperparameters are drawn from the distributions µ ∼ N (0, 1) and ϕ = 2 exp(ϕ * )/(1 + exp(ϕ * )) for ϕ * ∼ N (3, 1). The log-volatility obeys the order-one autoregressive process Z k+1 = µ + ϕ(Z k -µ) + ϵ k for all k > 1, where ϵ k ∼ N (0, 1) is independent of all other variables and Z 0 |µ, ϕ ∼ N (µ, 1 1-ϕ 2 ). While the states are Gaussian conditioned on the hyperparameters, the joint distribution of X = (µ, ϕ, Z 1 , . . . , Z T ), is non-Gaussian. In this example, the dimension of X is made arbitrarily large by increasing T .

A key property of triangular transport maps is that they inherit sparsity from the conditional independence structure of π. [START_REF] Spantini | Inference via low-dimensional couplings[END_REF] showed that the Markov structure of π yields a lower bound on the sparsity of the map S (i.e., the functional dependence of each component S k on the input variables x 1:k-1 ). This sparsity was exploited to learn the structure of undirected probabilistic graph models from data in [START_REF] Morrison | Beyond normality: Learning sparse probabilistic graphical models in the non-Gaussian setting[END_REF][START_REF] Baptista | Learning non-Gaussian graphical models via Hessian scores and triangular transport[END_REF]. From the conditional independence properties of the stochastic volatility model described above, we know that the KR rearrangement S KR between the joint distribution of X and a standard Gaussian reference η is sparse. Moreover, the exact sparsity of S KR can be derived from Theorem 5.1 in [START_REF] Spantini | Inference via low-dimensional couplings[END_REF]. We now show that the ATM algorithm is capable of detecting and exploiting this structure-without knowing in advance that it exists.

Figure 9a compares the variable dependence of the true KR rearrangement S KR and the map S learned by the ATM algorithm for a distribution with T = 40 using a sample of size n = 10 3 from π; a non-filled entry (j, k) entry of the plot indicates that k-th map component does not depend on variable x j . The dependence of component S k on (x k , x k-1 ) shows that each state Z k strongly depends on the previous state in time. Most of the map components also show dependence on the hyperparameters (µ, ϕ). The estimated sparsity closely matches the exact sparsity of the KR rearrangement. Furthermore, the sparse variable dependence of the kth map component S k on parent nodes Pa(k) ⊆ {1, . . . , k -1} produces an approximation to the kth marginal conditional density given by π(

x k |x <k ) = π(x k |x Pa(k) ).
By identifying the parent nodes Pa(k) of each variable k, we also learn a sparse Bayesian network or directed acyclic graphical (DAG) model representing the target distribution [START_REF] Koller | Probabilistic graphical models: principles and techniques[END_REF]. As a result, we can see the ATM algorithm as a technique for learning DAGs from samples, given a prescribed variable ordering.

Next, we consider the approximation of S KR for Markov models of increasing state dimension T and hence increasing map dimension d = T + 2; these experiments use a fixed sample size n = 10 3 as before. Figure 9b plots the KL divergence from the ATM approximation to the joint density of X ∈ R d , for increasing dimensions d. We compare ATM with a variant of ATM where the exact sparsity pattern of the the KR rearrangement is provided to the algorithm in advance. This variant, labeled "sparsity-aware ATM" in Figure 9b, differs from the ATM algorithm in that it only activates multi-indices α * t which match the sparsity of S KR (meaning of the form α = (α 1 , α 2 , 0, . . . , 0, α k-1 , α k )). We also compare ATM with non-adaptive maps of degree p = 1 and p = 2; these maps do not exploit the conditional independence structure of π, and hence depend on all input variables. For low dimensions, we see that degree p = 2 maps can better capture the non-Gaussian target distributions than p = 1 linear maps. As the dimension d increases, however, the growing number of coefficients in the p = 2 maps results in higher-variance estimators, outweighing their smaller bias, and hence we see a rapidly increasing KL divergence as the dimension d increases (with crossover around d = 21). In contrast, ATM outperforms the non-adaptive maps for all d and, more importantly, achieves a KL divergence that grows linearly with d; indeed, it performs similarly to the best case "oracle sparsity" of sparsity-aware ATM.

Tabular datasets

Lastly, we evaluate ATM's performance for density estimation on a suite of UCI datasets [START_REF] Lichman | UCI Machine Learning Repository[END_REF] that were also considered in [START_REF] Uria | RNADE: The real-valued neural autoregressive density-estimator[END_REF]. These datasets have dimensionalities between d = 10 and 15 and sample sizes between n = 506 and 5875. We pre-process each dataset to eliminate discrete-valued variables and one variable in every pair that has Pearson correlation coefficient greater than 0.98, following the procedure in [START_REF] Uria | RNADE: The real-valued neural autoregressive density-estimator[END_REF]. We consider 10 splits of the data. For each split, we use one fold (i.e., 10% of the data) as a test set and the remaining 9 folds (i.e., 90% of the data) as the training set to build S. To assess the quality of our estimated map S, we evaluate the negative log-likelihood of the pullback density π = S ♯ η on the test set. The negative log-likelihood is an empirical estimator of -E π [log S ♯ η] = D KL (π|| S ♯ η) -log πdπ and is, up to the unknown constant -log πdπ, the KL divergence from S ♯ η to π. Table 1 presents the mean of the negative log-likelihoods over the 10 splits and a 95% confidence interval for the mean. For all datasets, we observe an improvement using ATM over non-adaptive p = 2 maps and multivariate Gaussian approximations (i.e., p = 1 maps), with a lower number of total coefficients. Lastly, we evaluate ATM's performance for conditional density estimation on a different suite of UCI datasets. We follow a similar procedure as above to pre-process each dataset. Each dataset has a one-dimensional predictor variable X d and covariates X <d of varying dimension. To approximate the conditional density π d (x d |x <d ), we estimate one map component S d as a function of the predictor and the covariates using joint samples {(X i <k , X i k )} n i=1 . Table 2 presents the mean negative conditional log-likelihoods over 10 splits of the training data. The negative conditional log-likelihood is an empirical estimator of

-E π [log S ♯ d η] = E π(x <d ) [D KL (π d (x d |x <d )|| S ♯ d η)] -log π d (x d |x <d )dπ d and is, up to the unknown constant -log π d (x d |x <d )dπ d , the expected KL divergence from S ♯ d η to the marginal conditional π d (x d |x <d ).
In Table 2, we compare ATM to conditional kernel density estimation (CKDE) [START_REF] Silverman | On the estimation of a probability density function by the maximum penalized likelihood method[END_REF], ϵ-neighborhood kernel density estimation (NKDE), and kernel mixture networks (KMN) [START_REF] Ambrogioni | The kernel mixture network: A nonparametric method for conditional density estimation of continuous random variables[END_REF] using the implementation provided by [START_REF] Rothfuss | Conditional Density Estimation with Neural Networks: Best Practices and Benchmarks[END_REF]. We also include results for high-capacity parametric conditional density estimation methods that rely on neural networks: mixture density networks (MDN) [START_REF] Bishop | Mixture density networks Technical Report[END_REF] and conditional normalizing flows (NF) [START_REF] Trippe | Conditional density estimation with Bayesian normalising flows[END_REF]; implementation details of these methods are provided in Appendix C. For all datasets, we observe that ATM has performance comparable to the neural-network based methods and improved performance over the nonparametric approaches (e.g., on the Concrete and Yacht datasets).

In Table 3, we compare the number of coefficients in the ATM and MDN approximations and the computational time required to learn each of these conditional density representations. The runtime of the greedy procedure in Algorithm 1 (sequentially, for all d map components) is generally less than that of MDN, achieving similar approximation performance with about two orders of magnitude fewer coefficients. It is more realistic, however, to consider ATM with the cross-validation procedure of Algorithm 2 and to compare this with the hyperparameter tuning required to achieve the reported performance of MDN. ATM only requires performing cross-validation for a single hyper-parameter, i.e., the total number of coefficients m * (see Algorithm 2). In comparison, tuning several hyperparameters in the neural-network based methods (e.g., the learning rate, and the structure of the networks) by performing cross-validation over a tensor product grid increases runtimes by more than an order of magnitude relative to ATM, as presented below.

Table 2

Mean negative conditional log-likelihood for UCI datasets over 10 sets of training samples.

The method with best performance from both categories is highlighted in bold.

Dataset (d, n) ATM CKDE NKDE MDN KMN NF Boston (12,506) 2.6 ± 0.2 2.6 ± 0.2 3.1 ± 0.2 2.4 ± 0.2 2.7 ± 0.2 2.4 ± 0.1 Concrete (9,1030) 3.1 ± 0.1 3.2 ± 0.1 3.9 ± 0.1 2.9 ± 0.1 3.5 ± 0.1 3.2 ± 0.2 Energy (10,768) 1.5 ± 0.1 1.0 ± 0.1 2.1 ± 0.2 1.2 ± 0.1 1.7 ± 0.1 1.7 ± 0.3 Yacht (7,308) 0.5 ± 0.2 1.1 ± 0.3 3.8 ± 0.2 0.7 ± 0.2 1.8 ± 0.2 1.3 ± 0.5

Conclusions

This paper has presented and analyzed a functional framework for learning monotone triangular transport maps. Our approach represents monotone component functions of a triangular map through the action of an invertible rectification operator R k on smooth, generally non-monotone, functions.

Imposing appropriate structure on this operator and the function space V k on which it acts, along with conditions on the target density π, yields an unconstrained optimization problem for learning the map components that has many desirable and useful properties. First, under certain assumptions on R k and π, we show that the optimization objective is bounded, continuous, and differentiable for all f ∈ V k . This permits the use of deterministic gradient-based optimization methods to find the minimizers. Next, under certain conditions on R k , we show that the optimization problem has no spurious local minima on V k . In practice, this yields important robustness to the initial guess and other parameters of the optimization algorithms. Finally, under the same (cumulative) conditions on R k and some additional assumptions on the target density, we show that the optimization problem has a unique global minimizer on V k that corresponds to the canonical KR rearrangement.

Our functional framework also enables the construction of novel transport map estimators, based on maximum likelihood, given a finite sample from the target. The procedure we develop here is semi-parametric, making use of different hierarchical bases for V k . In particular, we propose a greedy, adaptive algorithm that identifies a sparse set of basis functions, automatically tailored to the target density and to the sample size. Sparsity in the basis selection reflects conditional independence structure in the target. More generally, the maps we build can capture the structure of complex probability distributions with sufficient data, but also are robust in settings with few observations. This is crucial to deploying these algorithms within large-scale applications, such as data assimilation, where triangular transport maps must be learned online [START_REF] Spantini | Coupling techniques for nonlinear ensemble filtering[END_REF]. We demonstrate good computational performance of our algorithm-and the parsimonious representations it produces-in a variety of joint and conditional density estimation problems. We outline some interesting future research directions below.

KR recovery for broader classes of target distributions While our result guaranteeing no spurious local minima (Theorem 7) makes no explicit assumptions on π, our subsequent results on recovery of the KR rearrangement (see Corollary 9) make particular use of Gaussian tail assumptions, relating the target π to the Gaussian reference η. In principle, however, the space V k contains functions R -1 k (S KR,k ) corresponding to a broader class of target densities; V k includes functions that grow up to sub-exponentially fast, as a result of its Gaussian weighting. Hence, it includes the inverse images of KR rearrangements for lighter-tailed π, whose component functions grow faster than linear as x k → ±∞. It may be possible to obtain unique recovery results analogous to Corollary 9 for such densities, although some modifications to our analysis would be needed.

Approximation theory and statistical guarantees

An open research topic, to the best of our knowledge, is to analyze the finite-dimensional approximation of triangular transport maps on unbounded domains. (See [START_REF] Zech | Sparse approximation of triangular transports. Part II: the infinite dimensional case[END_REF][START_REF] Zech | Sparse Approximation of Triangular Transports. Part I: The Finite-Dimensional Case[END_REF] for an approximation theory of triangular transports between analytic densities π, η on bounded domains.) These results would show how approximate maps converge to the KR rearrangement with different parameterizations (e.g., polynomial spaces or neural networks of increasing size), and we expect that our rectification framework could be a useful ingredient of such analyses. A parallel line of work is to develop non-asymptotic statistical convergence results for triangular transports (e.g., in the context of density estimation), to understand how the quality of density estimates and of map estimates depends on the sample size n. Combining these lines of work could provide lower bounds on the number of samples required to learn maps of a given complexity, and also some analytical guidance for how to navigate the bias-variance tradeoff of a map estimator given finite samples.

Other sources of low-dimensional structure The ATM algorithm in this paper identifies maps with sparse variable dependence by exploiting conditional independence structure in the target distribution [START_REF] Spantini | Inference via low-dimensional couplings[END_REF]. In future work it will be interesting to investigate other notions of low-dimensional structure and how they could facilitate the learning of maps from small samples. One promising notion is when the target density π departs from the reference η only along a low-dimensional subspace. In this case, the triangular map S pushing forward π to η can be written as a low-dimensional perturbation of the identity map, after a variable rotation. See [START_REF] Brennan | Greedy inference with structure-exploiting lazy maps[END_REF] for a recent contribution in this direction applied to variational inference.

Map ordering One disadvantage of triangular maps is that the approximation depends on the choice of variable ordering. In particular, each variable ordering yields a different factorization of the target density and a different Knothe-Rosenblatt rearrangement S KR . Thus, it is of interest to develop variable ordering algorithms that minimize the finite-sample error of the estimated transport map S, the estimated pullback density π = S ♯ η, or other goal-oriented metrics. For target distributions that induce sparse variable dependence in the map, one approach is to find the permutation that maximizes the sparsity of S. This is equivalent to finding the the sparsest Bayesian network or directed acyclic graphical (DAG) model for a distribution from samples. While this problem is in general NP-complete, effective algorithms [START_REF] Raskutti | Learning directed acyclic graph models based on sparsest permutations[END_REF] have been proposed. This algorithm reduces to finding a sparse maximum likelihood estimator for the Cholesky factor of the inverse covariance matrix of π in the Gaussian setting. Since a linear transport map is precisely this Cholesky factor for Gaussian π and standard normal η (see [START_REF] Baptista | Learning non-Gaussian graphical models via Hessian scores and triangular transport[END_REF]Section 3]), we expect that such sparse permutation algorithms could be generalized to the nonlinear transport map setting.

Nonparametric methods Instead of finding a particular finite-dimensional basis for the map components S k = R k (f ), nonparametric methods do not limit the functional form of f (e.g., as linear combinations of multivariate polynomials). A broadly useful nonparametric approach involves seeking f in a reproducing kernel Hilbert space (RKHS). Thanks to the representer theorem [START_REF] Schölkopf | A generalized representer theorem[END_REF], the optimal f in the RKHS can still be identified by solving a finite-dimensional optimization problem. Choosing the kernel so that the RKHS is a suitable weighted Sobolev space (see, e.g., [START_REF] Novak | Reproducing kernels of Sobolev spaces on R d and applications to embedding constants and tractability[END_REF]) will yield map representations that fall within the framework of Section 3.2. A related approach uses Gaussian process representations for the map components [START_REF] Katzfuss | Scalable Bayesian transport maps for high-dimensional non-Gaussian spatial fields[END_REF]. It will be interesting to compare the finite-sample performance of such methods to the semi-parametric procedure proposed in this work.

A.3. Proof of Proposition 2

To prove Proposition 2 we need the following lemma.

Lemma 11. Let H 1 ([0, 1]) = f : [0, 1] → R such that ∥f ∥ 2 H 1 ([0,1]) := 1 0 f (t) 2 + f ′ (t) 2 dt . Then |f (0)| ≤ √ 2∥f ∥ H 1 ([0,1]) , (38)
holds for any f ∈ H 1 ([0, 1]).

Proof. Because C ∞ ([0, 1]) is dense in H 1 ([0, 1]
), it suffices to show [START_REF] Nocedal | Numerical optimization[END_REF] for any f ∈ C ∞ ([0, 1]). By the mean value theorem, there exists 0 ≤ z ≤ 1 such that

f (z) = 1 1 -0 1 0 f (t)dt.
Thus we can write

|f (0)| 2 ≤ 2|f (z) -f (0)| 2 + 2|f (z)| 2 = 2 z 0 f ′ (t)dt 2 + 2 1 0 f (t)dt 2 ≤ 2 1 0 f ′ (t) 2 dt + 2 1 0 |f (t)| 2 dt.
This concludes the proof.

We now prove Proposition 2.

Proof. For any f ∈ V k , Lemma 11 permits us to write

|f (x <k , 0)| 2 η <k (x <k )dx <k (38)
≤ 2

1 0 |f (x <k , t)| 2 + |∂ k f (x <k , t)| 2 dt η <k (x <k )dx <k ≤ C T 1 0 |f (x <k , t)| 2 + |∂ k f (x <k , t)| 2 η <k (x <k )η 1 (t)dx <k dt ≤ C T +∞ -∞ |f (x <k , t)| 2 + |∂ k f (x <k , t)| 2 η ≤k (x <k , t)dx <k dt = C T ∥f ∥ 2 V k ,
where C T = 2 sup 0≤t≤1 η 1 (t) -1 .

A.4. Proof of Proposition 3

The proof relies on Proposition 2 and on the following generalized integral Hardy inequality, see [START_REF] Muckenhoupt | Hardy's inequality with weights[END_REF].

Lemma 12. Let η ≤k be the standard Gaussian density on R k . Then there exists a constant C H such that for any v ∈ L 2 η (R k ),

x k 0 v(x <k , t)dt 2 η ≤k (x)dx ≤ C H v(x) 2 η ≤k (x)dx. ( 39
)
Proof of Lemma 12. Let us recall the integral Hardy inequality [START_REF] Muckenhoupt | Hardy's inequality with weights[END_REF].

Theorem 13 (from [START_REF] Muckenhoupt | Hardy's inequality with weights[END_REF]). For weight ρ : R + → R + and u ∈ L 2 ρ (R), there exists a constant

C H < ∞ such that +∞ 0 x 0 u(t)dt 2 ρ(x)dx ≤ C H +∞ 0 u(x) 2 ρ(x)dx (40)
if and only if

sup x>0 +∞ x ρ(t)dt 1/2 x 0 ρ(t) -1 dt 1/2 < +∞. (41)
We apply Theorem 13 with the one-dimensional standard Gaussian density ρ = η for x > 0. In order to check condition [START_REF] Papamakarios | Fast ε-free inference of simulation models with Bayesian conditional density estimation[END_REF], we need to show that

D(x) := +∞ x ρ(t)dt 1/2 x 0 ρ(t) -1 dt 1/2 = +∞ x e -t 2 /2 dt 1/2 x 0 e t 2 /2 dt 1/2
, is bounded. Since x → D(x) is a continuous function with a finite limit as x → 0, it is sufficient to show that D(x) has a finite limit when x → ∞. For x > 1,

+∞ x e -t 2 /2 dt ≤ e -x 2 /2 and D(x) 2 ≤ e -x 2 /2 x 0 e t 2 /2 dt. Furthermore, using integration-by-parts we have

x 0 e t 2 /2 dt = 1 0 e t 2 /2 dt + e x 2 /2 /x -√ e +

x 1 e t 2 /2 /t 2 dt. As x → ∞ the dominating term in the sum is e x 2 /2 /x. Thus, e -x 2 /2 x 0 e t 2 /2 dt behaves asymptotically as O( 1x ), so that D(x) → 0 when x → ∞. Thus, condition (41) is satisfied. Then, by the Hardy inequality in [START_REF] Oord | Parallel WaveNet: Fast high-fidelity speech synthesis[END_REF]

for u ∈ L 2 η (R) we have +∞ 0 x k 0 u(t)dt 2 η(x k )dx k ≤ C H +∞ 0 u(x k ) 2 η(x k )dx k . (42)
For the symmetric density η(x k ) = η(-x k ) we also have

0 -∞ x k 0 u(t)dt 2 η(x k )dx k ≤ C H 0 -∞ u(x k ) 2 η(x k )dx k . (43)
Combining the results in ( 42) and ( 43) we have

+∞ -∞ x k 0 u(t)dt 2 η(x k )dx k ≤ C H +∞ -∞ u(x k ) 2 η(x k )dx k .
Setting u(t) = v(x <k , t) and integrating both sides over x <k ∈ R k-1 with the standard Gaussian weight function η(x <k ) gives the result.

We now prove Proposition 3.

Proof. By Proposition 2, Lemma 12 and by the Lipschitz property of g, we can write

∥R k (f 1 ) -R k (f 2 )∥ 2 L 2 η ≤k ≤ 2 f 1 (x <k , 0) -f 2 (x <k , 0) 2 η ≤k (x)dx + 2 x k 0 g ∂ k f 1 (x <k , t) -g ∂ k f 2 (x <k , t) dt 2 η ≤k (x)dx ≤ 2C T ∥f 1 -f 2 ∥ 2 V k + 2C H ∥g(∂ k f 1 ) -g(∂ k f 2 )∥ 2 L 2 η ≤k ≤ 2C T ∥f 1 -f 2 ∥ 2 V k + 2C H L 2 ∥∂ k f 1 -∂ k f 2 ∥ 2 L 2 η ≤k ≤ 2(C T + C H L 2 )∥f 1 -f 2 ∥ 2 V k , (44)
for any f 1 , f 2 ∈ V k . Furthermore, using the Lipschitz property of g we have

∥∂ k R k (f 1 ) -∂ k R k (f 2 )∥ 2 L 2 η ≤k = g(∂ k f 1 (x <k , t)) -g(∂ k f 2 (x <k , t)) 2 η ≤k (x)dx ≤ L 2 ∂ k f 1 (x <k , t) -∂ k f 2 (x <k , t) 2 η ≤k (x)dx ≤ L 2 ∥f 1 -f 2 ∥ 2 V k . (45)
Combining ( 44) with ( 45), we obtain [START_REF] Durkan | Neural spline flows[END_REF] with

C = 2(C T + C H L 2 ) + L 2 . It remains to show that ∥R k (f )∥ V k < ∞ for any f ∈ V k .
Letting with f 1 = f and f 2 = 0 in (18), the triangle inequality yields

∥R k (f )∥ V k ≤ ∥R k (0)∥ V k + C∥f ∥ V k Because R k (0) is the affine function x → g(0)x k , we have that ∥R k (0)∥ 2 L 2 η ≤k = g(0) 2 x 2 k η(x)dx and ∥∂ k R k (0)∥ 2 L 2 η ≤k = g(0) 2 are finite and so is ∥R k (0)∥ V k . Thus, R k (f ) ∈ V k for all f ∈ V k . A.5. Proof of Proposition 4 Proof. For any f ∈ V k we have |L k (f )| = 1 2 R k (f ) 2 -log(∂ k R k (f )) dπ (19)
≤ C π 2 ∥R k (f )∥ 2 L 2 η ≤k + C π |log(g(∂ k f ))| dη ≤k ≤ C π 2 ∥R k (f )∥ 2 L 2 η ≤k + C π | log(g(0))| + |log(g(∂ k f )) -log(g(0))| dη ≤k (21) ≤ C π 2 ∥R k (f )∥ 2 L 2 η ≤k + C π | log(g(0))| + C π L |∂ k f -0| dη ≤k ≤ C π 2 ∥R k (f )∥ 2 L 2 η ≤k + C π | log(g(0))| + C π L∥f ∥ 2 V k . Because Proposition 3 ensures R k (f ) ∈ V k ⊂ L 2 η ≤k , we have that L k (f ) is finite for any f ∈ V k . Now, for any f 1 , f 2 ∈ V k , we can write |L k (f 1 ) -L k (f 2 )| = 1 2 R k (f 1 ) 2 - 1 2 R k (f 2 ) 2 -log(∂ k R k (f 1 )) + log(∂ k R k (f 2 )) dπ (19) ≤ C π 1 2 R k (f 1 ) 2 -R k (f 2 ) 2 + log(g(∂ k f 1 )) -log(g(∂ k f 2 )) dη (21) ≤ C π 2 ∥R k (f 1 ) + R k (f 2 )∥ L 2 η ≤k ∥R k (f 1 ) -R k (f 2 )∥ L 2 η ≤k + C π L∥∂ k f 1 -∂ k f 2 ∥ L 2 η ≤k (18) ≤ C π ∥R k (f 1 )∥ L 2 η ≤k + ∥R k (f 2 )∥ L 2 η ≤k 2 C∥f 1 -f 2 ∥ V k + C π L∥f 1 -f 2 ∥ V k . This shows that L k : V k → R is continuous. To show that L k is differentiable, we let f, ε ∈ V k so that L k (f + ε) = 1 2 R k (f + ε) 2 -log(∂ k R k (f + ε)) dπ = 1 2 f (x <k , 0) + ε(x <k , 0) + x k 0 g ∂ k (f + ε)(x <k , t) dt 2 -log •g(∂ k (f + ε)(x)) π(x)dx = 1 2 R k (f )(x) 2 + R k (f )(x) ε(x <k , 0) + x k 0 g ′ ∂ k f (x <k , t) ∂ k ε(x <k , t)dt π(x)dx -log •g(∂ k f ) + (log •g) ′ (∂ k f )∂ k εdπ + O(∥ε∥ 2 V k ) = L k (f ) + ℓ(ε) + O(∥ε∥ 2 V k )
where ℓ : V k → R is the linear form defined by

ℓ(ε) = R k (f )(x) ε(x <k , 0) + x k 0 g ′ ∂ k f (x <k , t) ∂ k ε(x <k , t)dt -(log •g) ′ (∂ k f (x))∂ k ε(x) π(x)dx.
If ℓ is continuous, meaning if there exists a constant C ℓ such that |ℓ(ε)| ≤ C ℓ ∥ε∥ V k for any ε ∈ V k , then the Riesz representation theorem states that there exists a vector

∇L k (f ) ∈ V k such that ℓ(ε) = ⟨∇L k (f ), ε⟩ V k .
This proves L k is differentiable everywhere.

To show that ℓ is continuous, we write

|ℓ(ε)| (19) ≤ C π R k (f )(x) ε(x <k , 0) + x k 0 g ′ ∂ k f (x <k , t) ∂ k ε(x <k , t)dt η ≤k (x)dx + C π (log •g) ′ (∂ k f (x))∂ k ε(x) η ≤k (x)dx (21) ≤ C π ∥R k (f )∥ L 2 η ≤k ε(x <k , 0) + x k 0 g ′ ∂ k f (x <k , t) ∂ k ε(x <k , t)dt 2 η ≤k (x)dx + C π L∥∂ k ε∥ L 2 η ≤k (20) ≤ C π ∥R k (f )∥ L 2 η ≤k 2C T ∥ε∥ 2 V k + 2C H L 2 ∂ k ε(x)) 2 η ≤k (x)dx + C π L∥∂ k ε∥ L 2 η ≤k ≤ C π ∥R k (f )∥ L 2 η ≤k 2C T + 2C H L 2 + L ∥ε∥ V k ,
where the second last inequality also uses Proposition 2 and Lemma 12. This concludes the proof.

A.6. Proof of the local Lipschitz regularity (23)

Proposition 14. In addition to the assumptions of Theorem 4, we further assume there exists a constant L < ∞ such that for all ξ, ξ ′ ∈ R we have

|g ′ (ξ) -g ′ (ξ ′ )| ≤ L|ξ -ξ ′ | (46) |(log •g) ′ (ξ) -(log •g) ′ (ξ ′ )| ≤ L|ξ -ξ ′ |. (47)
Then there exists M < ∞ such that

∥∇L k (f 1 ) -∇L k (f 2 )∥ V k ≤ M (1 + ∥R k (f 2 )∥ V k )∥f 1 -f 2 ∥ V k , for any f 1 , f 2 ∈ V k , where V k = {f ∈ V k , ∂ k f ∈ L ∞ } is the space endowed with the norm ∥f ∥ V k = ∥f ∥ V k + ∥∂ k f ∥ L ∞ . Proof. Recall the definition (22) of ∇L k (f ) ⟨∇L k (f ), ε⟩ V k = R k (f )(x) ε(x <k , 0) + x k 0 g ′ ∂ k f (x <k , t) ∂ k ε(x <k , t)dt π(x)dx -(log •g) ′ (∂ k f (x))∂ k ε(x)π(x)dx.
Then for any f 1 , f 2 ∈ V k . we can write

⟨∇L k (f 1 ) -∇L k (f 2 ), ε⟩ V k = A + B + C + D,
where

A = R k (f 1 )(x) -R k (f 2 )(x) ε(x <k , 0)π(x)dx B = R k (f 1 )(x) -R k (f 2 )(x) x k 0 g ′ (∂ k f 1 (x <k , t))∂ k ε(x <k , t)dt π(x)dx C = R k (f 2 )(x) x k 0 g ′ (∂ k f 1 (x <k , t)) -g ′ (∂ k f 2 (x <k , t)) ∂ k ε(x <k , t)dt π(x)dx D = (log •g) ′ (∂ k f 1 (x)) -(log •g) ′ (∂ k f 2 (x)) ∂ k ε(x)π(x)dx.
For the first term A we write

|A| (19) ≤ C π R k (f 1 )(x) -R k (f 2 )(x) |ε(x <k , 0)|η(x)dx ≤ C π ∥R k (f 1 ) -R k (f 2 )∥ V k |ε(x <k , 0)| 2 η <k (x)dx 1/2 (15) ≤ C π C T ∥R k (f 1 ) -R k (f 2 )∥ V k ∥ε∥ V k (18) ≤ C π C T C∥f 1 -f 2 ∥ V k ∥ε∥ V k .
For the second term B we write

|B| (19) ≤ C π ∥R k (f 1 ) -R k (f 2 )∥ V k x k 0 g ′ (∂ k f 1 (x <k , t))∂ k ε(x <k , t)dt 2 η(x)dx 1/2 (39) (18)
≤ C π C H C∥f 1 -f 2 ∥ V k g ′ (∂ k f 1 (x ≤k ))∂ k ε(x ≤k ) 2 η(x)dx 1/2 (20) ≤ C π C H CL∥f 1 -f 2 ∥ V k ∂ k ε(x ≤k ) 2 η(x)dx 1/2 ≤ C π C H CL∥f 1 -f 2 ∥ V k ∥ε∥ V k .
For the third term C we write

|C| (19) ≤ C π ∥R k (f 2 )∥ V k x k 0 g ′ (∂ k f 1 (x <k , t)) -g ′ (∂ k f 2 (x <k , t)) ∂ k ε(x <k , t)dt 2 η(x)dx 1/2 (39) ≤ C π C H ∥R k (f 2 )∥ V k g ′ (∂ k f 1 (x ≤k )) -g ′ (∂ k f 2 (x ≤k )) ∂ k ε(x ≤k ) 2 η(x)dx 1/2 ≤ C π C H ∥R k (f 2 )∥ V k ess sup g ′ • ∂ k f 1 -g ′ • ∂ k f 2 (∂ k ε(x ≤k )) 2 η(x)dx 1/2 (46)
≤ C π C H L∥R k (f 2 )∥ V k ess sup ∂ k f 1 -∂ k f 2 ∥ε∥ V k ≤ C π C H L∥R k (f 2 )∥ V k ∥f 1 -f 2 ∥ V k ∥ε∥ V k .
For the last term D we write

|D| (19) ≤ C π (log •g) ′ (∂ k f 1 (x)) -(log •g) ′ (∂ k f 2 (x)) 2 η(x)dx 1/2 ∥ε∥ V k (47) ≤ C π L ∂ k f 1 (x) -∂ k f 2 (x) 2 η(x)dx 1/2 ∥ε∥ V k ≤ C π L∥f 1 -f 2 ∥ V k ∥ε∥ V k . Thus, because ∥f 1 -f 2 ∥ V k ≤ ∥f 1 -f 2 ∥ V k we obtain |⟨∇L k (f 1 ) -∇L k (f 2 ), ε⟩ V k | ∥ε∥ V k ≤ C π C T C + C H CL + C H L∥R k (f 2 )∥ V k + L ∥f 1 -f 2 ∥ V k ≤ M (1 + ∥R k (f 2 )∥ V k )∥f 1 -f 2 ∥ V k , where M = C π max{ C T C + C H CL + L; C H L}.
This concludes the proof.

A.7. Proof of Proposition 6

Proof.

To show that R k (V k ) = {R(f ) : f ∈ V k } is convex, let f 1 , f 2 ∈ V k and 0 ≤ α ≤ 1.
We need to show that there exists f α ∈ V k such that R(f α ) = S α where

S α := αR k (f 1 ) + (1 -α)R k (f 2 ).
Let f α (x ≤k ) := R -1 (S α )(x ≤k ) = S α (x <k , 0) +

x k 0 g -1 (∂ k S α (x <k , t))dt.

It remains to show that f α ∈ V k , meaning that f α ∈ L 2 η ≤k and ∂ k f α ∈ L 2 η ≤k . By convexity of ξ → g -1 (ξ) 2 , we have

∥∂ k f α ∥ 2 L 2 η ≤k = g -1 α∂ k R k (f 1 ) + (1 -α)∂ k R k (f 2 ) 2 dη ≤k = g -1 αg(∂ k f 1 ) + (1 -α)g(∂ k f 2 ) 2 dη ≤k ≤ αg -1 (g(∂ k f 1 )) 2 + (1 -α)g -1 (g(∂ k f 2 )) 2 dη ≤k = α∥∂ k f 1 ∥ 2 L 2 η ≤k + (1 -α)∥∂ k f 2 ∥ 2 L 2 η ≤k . ( 48
)
Thus ∂ k f α ∈ L 2 η ≤k . Furthermore we have To show that the above quantity is finite, Proposition 2 permits us to write S α (x <k , 0) 2 η <k (x <k )dx = αf 1 (x <k , 0) + (1 -α)f 2 (x <k , 0) , which is finite. We deduce that f α ∈ L 2 η ≤k and therefore that f α ∈ V k .

∥f α ∥ 2
2 η <k (x <k )dx ≤ C T ∥αf 1 + (1 -α)f 2 ∥ 2 V k ,
A.8. Proof of Proposition 5

Proof. Let s 1 , s 2 ∈ V k be strictly increasing functions with respect to x k that satisfy ∂ k s i (x ≤k ) ≥ c for i = 1, 2 and all x ≤k ∈ R k . By the Lipschitz property of g -1 on the domain [c, ∞) with constant L c , we can write

∥∂ k R -1 k (s 1 ) -∂ k R -1 k (s 2 )∥ 2 L 2 η ≤k = (g -1 (∂ k s 1 (x ≤k )) -g -1 (∂ k s 2 (x ≤k ))) 2 η ≤k (x)dx ≤ L 2 c (∂ k s 1 (x ≤k ) -∂ k s 2 (x ≤k )) 2 η ≤k (x)dx ≤ L 2 c ∥s 1 -s 2 ∥ 2 V k . ( 49
)
Applying Proposition 2 to s 1 , s 2 ∈ V k and Lemma 12 to

∂ k R -1 k (s i ) = g -1 (∂ k s i ) ∈ L 2 η ≤k for i = 1, 2
we have

∥R -1 k (s 1 ) -R -1 k (s 2 )∥ 2 L 2
η ≤k ≤ 2 (s 1 (x <k , 0) -s 2 (x <k , 0)) 2 η ≤k (x)dx + 2

x k 0 g -1 (∂ k s 1 ) -g -1 (∂ k s 2 )dt 2 η ≤k (x)dx ≤ 2C T ∥s 1 -s 2 ∥ 2 V k
+ 2C H (g -1 (∂ k s 1 (x ≤k )) -g -1 (∂ k s 2 (x ≤k ))) 2 η ≤k (x)dx

≤ (2C T + 2C H L 2 c )∥s 1 -s 2 ∥ 2 V k , (50)
where the last inequality follows from [START_REF] Rosenblatt | Remarks on a multivariate transformation[END_REF].

It remains to show that ∥R -1 k (s)∥ V k < ∞ for any s ∈ V k such that ess inf ∂ k s > 0. Letting s 1 = s and s 2 = g(0)x k , the triangle inequality combined with (25) yields = 2g(0) 2 is finite, and so

∥R -1 k (s)∥ V k ≤ ∥R -1 k (g(0)x k )∥ V k + C c ∥s -g(0)x k ∥ V k . The function R -1 k (g(0)x k ) is zero. Therefore, ∥R -1 k (s)∥ V k ≤ C c ∥s -g(0)x k ∥ V k ≤ C c (∥s∥ V k + ∥g(0)x k ∥ V k ).
∥R -1 k (s)∥ V k < ∞ for s ∈ V k . Furthermore, if ∂ k s ≥ c > 0, then ∂ k R -1 k (s) = g -1 (∂ k s) ≥ g -1 (c) > -∞ and so ess inf R -1 k (s) > -∞.
A.9. Proof for the KR rearrangement

Proof. Let S KR,k be the kth component of the KR rearrangement, given by composing the inverse CDF of the standard Gaussian marginal F η,k (x k ) with the CDF of the target's kth marginal conditional F π k (x k |x <k ). That is,

S KR,k (x ≤k ) = F -1 η k • F π k (x k |x <k ). ( 51
)
The goal is to show S KR,k ∈ V k , that is, S KR,k ∈ L 2 η ≤k and ∂ k S KR,k ∈ L 2 η ≤k . First we show S KR,k ∈ L 2 η ≤k . From condition (29), we have

F -1 η k (c k F η k (x k )) ≤ S KR,k (x k |x <k ) ≤ F -1 η k (C k F η k (x k )) so that S KR,k (x k |x <k ) 2 ≤ max{F -1 η k (c k F η k (x k )) 2 ; F -1 η k (C k F η k (x k )) 2 }, (52)
for all x <k ∈ R k-1 . To show that S KR,k ∈ L 2 η ≤k , it is sufficient to prove that any function of the form x k → F -1 η k (CF η k (x k )) is in L 2 η ≤k for any C > 0. From Theorems 1 and 2 in [START_REF] Chang | Chernoff-type bounds for the Gaussian error function[END_REF], there exists strictly positive constants α i , β i > 0 for i = 1, 2 such that

α 1 exp(-β 1 x 2 k ) ≤ F η k (x k ) ≤ α 2 exp(-β 2 x 2 k ), (53)
for x k < 0. With a change of variable u = F η k (x k ) we obtain F -1 η k (u) 2 ≤ 1/β 2 log(α 2 /u) for all u < F η k (0) = 1/2. Letting u = CF η k (x k ) yields

F -1 η k (CF η k (x k )) 2 ≤ 1 β 2 log α 2 CF η k (x k ) , (53)
∂ k S KR,k (x ≤k ) = π k (x k |x <k ) η k (S KR,k (x ≤k )) = π k (F -1 π k (F π k (x k |x <k )|x <k )|x <k ) η k (F -1 η k • F π k (x k |x <k )) , (54)
is continuous, where F -1 π k (•|x <k ) denotes the inverse of the map x k → F π k (x k |x <k ) for each x <k ∈ R k-1 . Hence, it is sufficient to show that ∂ k S KR,k goes to a finite limit as |x k | → ∞. For the left-hand limit, we can write

lim x k →-∞ ∂ k S KR,k (x ≤k ) = lim u→0 + π k (F -1 π k (u|x <k )|x <k ) η k (F -1 η k (u)) = lim u→0 + (F -1 η,k ) ′ (u) (F -1 π k ) ′ (u|x <k ) = lim u→0 + F -1 η k (u) F -1 π k (u|x <k ) , (55)
where in the second equality we used the implicit function theorem and the third equality follows from l'Hôpital's rule. To analyze the ratio F -1 η k (u)/F -1 π k (u|x <k ), we combine the lower bound in (29) and the bounds in [START_REF] Schmuland | Dirichlet forms with polynomial domain[END_REF] to get

F -1 η k (u) F -1 π k (u|x <k ) ≤ F -1 η k (u) F -1 η k (u/c k ) ≤ β 2 (log α 1 -log(u)) β 1 (log α 2 -log(u/c k ))
.

Similarly, from the upper bound in [START_REF] Lezcano Casado | Trivializations for gradient-based optimization on manifolds[END_REF] and the bounds in (53), we have

F -1 η k (u) F -1 π k (u|x <k ) ≥ F -1 η k (u) F -1 η k (u/C k ) ≥ β 1 (log α 2 -log(u)) β 2 (log α 1 -log(u/C k ))
.

Thus, ∂ k S KR,k (x ≤k ) = O(1) as x k → -∞, and we have ∂ k S KR,k ∈ L 2 η k . Lastly, taking the limit in [START_REF] Shin | Joint estimation of monotone curves via functional principal component analysis[END_REF] we have lim x k →-∞ ∂ k S KR,k (x ≤k ) ≥ β 1 /β 2 . For a target distribution π with full support, all marginal conditional densities satisfy π k (x k |x <k ) > 0 for each x ≤k ∈ R k . Given that the ∂ k S KR,k does not approach zero as |x k | → ∞, we can find a strictly positive constant c k > 0 such that ∂ k S KR,k (x ≤k ) ≥ c k for all x ≤k ∈ R k . This shows that ess inf ∂ k S KR,k > 0.

  Another important feature of the triangular structure (1) is that each component of the map represents one marginal conditional density of π. More precisely, in the present setting where η is a product density, S KR,k pushes forward the marginal conditional π k (x k |x <k ) to the k-th marginal of the reference η k (z k ) [51, Section 2.3]. Now partition x = (y, w), where y ∈ R m and w ∈ R p , with d = m + p. This property of the KR map lets us estimate the conditional probability density function π(w|y), for any value of y, given a sample

Fig 1 :

 1 Fig 1: The rectifier (2) transforms the non-monotone function f into the monotone function S = R(f ).Here we choose g(•) = log(1 + exp(•)). S is an increasing transport that pushes forward a onedimensional mixture of Gaussians π(x) = 0.5N (x; -1, 1) + 0.5N (x; 1, 1) to the standard Gaussian reference density η.

Fig 2 :

 2 Fig 2: Objective function L 1 = J 1 • R 1 where the rectifier R 1 is defined using the soft-plus g (11) (left), the shifted ELU g (12) (middle), or the square function g(ξ) = ξ 2 (right). Here, π(x) = 1/2N (x; -2, 0.5) + 1/2N (x; 2, 2) is a univariate Gaussian mixture, and we use n = 50 to define J 1 with f 1 represented using a linear combination of Hermite functions up to degree 10. The objective is evaluated along line segments that interpolate between random initial maps (t = 0) and critical points resulting from a gradient-based optimization method (t = 1). Observe that with bijective g (left and middle) the algorithm always arrives at the same optimal value, whereas with the square function g (right) the algorithm gets stuck in local minima and rarely attains the optimal value.

Fig 3 :

 3 Fig 3: A k = 2 dimensional downward-closed active set of multi-indices Λ t with its margin Λ Mt and reduced margin Λ RM t . The margin and reduced margins are plotted before (left) and after (right) adding to α * t = (2, 1), that is denoted with a cross, to Λ t .

Algorithm 2 1 := χ \ χ test j 5 :| iterations with sample χ train j 6 :

 2156 ATM algorithm for learning map component S k Input: Training sample χ = {X i 1:k } n i=1 , number of folds ν 2: Partition data into ν folds of equal size 3: for j = 1, . . . , ν do 4: Partition χ: χ test j is the jth subset of χ and χ train j Estimate ft using Algorithm 1 for m = |χ train j Store iterates S j,1 k := R k (f1), . . . , S j,m k := R k (fm) 7:

Fig 4 :

 4 Fig 4: (a) The pullback density S ♯ η approaches the target Gaussian mixture density π when increasing the maximum polynomial degree p of the space V p 1 ⊂ V 1 . (b) With increasing p, the pullback density converges to π in KL divergence, and the estimated map converges to S KR in L 2 π .

Fig 5 :Fig 6 :

 56 Fig 5: (a) The approximate transport maps S compared to S KR (black), and (b) the corresponding non-monotone functions f compared to f KR := R -1 (S KR ). Both subfigures illustrate different choices of basis {ψ α } α , for the Gaussian mixture target of Figure 4. The modified Hermite polynomial basis provides the closest approximation to S KR and f KR .

Fig 7 :

 7 Fig 7: Five different densities π (top row) and their approximations π (bottom row) built using the ATM algorithm with a Hermite polynomial basis for V k , given a sample of size n = 10 4 from π. The total number of coefficients #Λ m * in both map components is indicated below each density.

Fig 8 :

 8 Fig 8: (a) KL divergence over 10 sets of training samples is lower for ATM than for non-adapted expansions. (b) Trade-off between the approximation quality and number of coefficients using the adaptive algorithm for different n. For any given number of map coefficients, the ATM algorithm finds a representation with similar or lower KL divergence than the lowest achievable error of each non-adaptive approximation (indicated with stars in plot (b)).

Fig 9 :

 9 Fig 9: (a) The sparsity of an ATM map S and of the KR rearrangement S KR for the stochastic volatility model with T = 40 time steps. (b) KL divergence using ATM with increasing dimension d, compared to adaptive map estimators with known sparsity and non-adaptive maps estimators.

L 2 η 2 η

 22 ≤k = S α (x <k , 0) + x k 0 g -1 (∂ k S α (x <k , t))dt ≤k (x)dx ≤ 2 S α (x <k , 0) 2 η <k (x <k )dx + 2 x k 0 g -1 (∂ k S α (x <k , t))dt 2 η ≤k (x)dx.

0 g - 1 ( 2 η 2 L 2 η

 01222 which is finite. Finally, because g -1 (∂ k S α ) = ∂ k f α ∈ L 2η ≤k by (48), Lemma 12 yieldsx k ∂ k S α (x <k , t))dt ≤k (x)dx ≤ C H g -1 (∂ k S α (x ≤k )) 2 η ≤k (x)dx = C H ∥∂ k f α ∥ ≤k

  For a linear function, ∥g(0)x k ∥ 2 V k = ∥g(0)x k ∥ 2

  Λ t with its margin Λ M t and reduced margin Λ RM t . The margin and reduced margins are plotted before (left) and after (right) adding to α * t = (2, 1), that is denoted with a cross, to Λ t .

	Algorithm 1 Estimate map component S k
	1: Input: training sample {X i 1:k } n i=1 , maximum cardinality m for Λt
	2: Initialize Λ0 = ∅, f0 = 0
	3: for t = 0, . . . , m -1 do
	4:	Construct the reduced margin: Λ RM t
	5:	Select the new multi-index: α * t ∈ arg max α∈Λ RM t	|∇α L k (ft)|
	6:

Table 1

 1 Mean negative log-likelihood for UCI datasets over 10 sets of training samples. The estimator with best performance (lowest negative log-likelihood) is highlighted in bold.

	Dataset	(d, n)	Gaussian	ATM	# coefficients p = 2 maps # coefficients
	Boston	(10, 506) 11.3 ± 0.5 3.1 ± 0.6	228 ± 7	6.5 ± 0.4	285
	Red Wine	(11, 1599) 13.2 ± 0.3 9.8 ± 0.4	289 ± 9	10.5 ± 0.2	363
	White Wine (11, 4898) 13.2 ± 0.5 11.0 ± 0.2	342 ± 26	12.0 ± 1.0	363
	Parkinsons (15, 5875) 10.8 ± 0.4 2.8 ± 0.4	783 ± 17	5.1 ± 0.4	815

Table 3

 3 The number of coefficients in the ATM and MDN approximations, as well as the runtimes (in seconds) to identify coefficients via optimization and to perform cross-validation of hyperparameters.

			ATM			MDN
	Dataset	# coefficients runtime (s) runtime w/CV	# coefficients	runtime (s) runtime w/CV
	Boston	31 ± 3	6 ± 1	88 ± 14	(5 ± 2) × 10 3	18 ± 8	1584 ± 13
	Concrete	42 ± 5	14 ± 3	230 ± 18	(2.6 ± 0.3) × 10 3	11 ± 1	1846 ± 10
	Energy	41 ± 7	12 ± 5	193 ± 29	(8 ± 3) × 10 3	32 ± 12	1758 ± 17
	Yacht	28 ± 4	7 ± 2	99 ± 20	(6 ± 2) × 10 3	23 ± 6	1453 ± 15

  for all x k < min{F -1 η k (1/(2C)), 0}. Using the same argument, we obtain a similar bound on F -1 η k (CF η k (x k )) 2 for all x k larger than a certain value. Together with the continuity ofx k → F -1 η k (CF η k (x k )) 2 these bounds ensure that x k → F -1 η k (CF η k (x k )) 2 is in L 2 η ≤k for any C. Then S KR,k ∈ L 2 η ≤k . Furthermore, we have S KR,k (x ≤k ) = O(x k ) as |x k | → ∞. Now we show that ∂ k S KR,k ∈ L 2η ≤k by showing ∂ k S KR,k is a continuous and bounded function. From the absolute continuity of µ and ν, we have that

	≤	1 β 2	log	α 2 Cα 1 exp(-β 1 x 2 k )
	=	1 β 2	log	α 2 Cα 1	+	β 1 β 2	x 2 k

For any z ∈ R d , x = S -1 (z) can be computed recursively by x k = T k (x <k , z k ) for k = 1, . . . , d, where the functionT k (x <k , •) is the inverse of x k → S k (x <k , x k ).

https://github.com/baptistar/ATM

Acknowledgments

RB, YM, and OZ gratefully acknowledge support from the INRIA associate team Unquestionable. RB and YM are also grateful for support from the AFOSR Computational Mathematics program (MURI award FA9550-15-1-0038) and the US Department of Energy AEOLUS center. RB acknowledges support from an NSERC PGSD-D fellowship. OZ also acknowledges support from the ANR JCJC project MODENA (ANR-21-CE46-0006-01).

Proof. Recall that the KR rearrangement S KR is a transport map that satisfies S ♯ KR η = π, where η is the density of the standard Gaussian measure on R d and π is the target density. Corollary 3.10 in [START_REF] Bogachev | Triangular transformations of measures[END_REF] states that for any PDF ϱ on R d of the form ϱ(x)

holds, where T is the KR rearrangement such that T ♯ η = ϱ. Let S be an increasing lower triangular map as in [START_REF] Ambrogioni | The kernel mixture network: A nonparametric method for conditional density estimation of continuous random variables[END_REF] and let ϱ = S ♯ π. Thus we have T = S • S -1 KR and so the left-hand side of (37) becomes

and the right-hand side becomes

which yields (4).

A.2. Convexity of s → J k (s)

Lemma 10. The optimization problem min {s :

For any 0 < t < 1, the function s t = ts 1 + (1 -t)s 2 is also strictly increasing with respect to x k . Finally, because both ξ → 1 2 ξ 2 and ξ → -log(ξ) are strictly convex functions, we have

which shows that J k is strictly convex.

Appendix B: Multi-index refinement for the wavelet basis

In this section we show how to greedily enrich the index set Λ t for a one-dimensional wavelet basis parameterized by the tuple of indices (l, q) representing the level l and translation q of each wavelet ψ (l,q) . To define the allowable indices, we construct a binary tree where each node is indexed by (l, q) and has two children with indices (l + 1, 2q) and (l + 1, 2q + 1). The root of the tree has index (0, 0) and corresponds to the mother wavelet ψ. Analogously to the downward closed property for polynomial indices, we only add nodes to the tree (i.e., indices in Λ t ) if its parents have already been added. Given any set Λ t , we define its reduced margin as

for odd q (l -1, (q -1)/2) ∈ Λ t for even q .

Then, the ATM algorithm with a wavelet basis follows from Algorithm 1 with this construction for the reduced margin at each iteration.

Appendix C: Architecture details of alternative methods

In this section we present the details of the alternative methods to ATM that we consider in Section 5.

For each normalizing flow model, we use the recommended stochastic gradient descent optimizer with a learning rate of 10 -3 . We partition 10% of the samples in each training set to be validation samples and use the remaining samples for training the model. We select the optimal hyper-parameters for each dataset by fitting the density with the training data and choosing the parameters that minimize the negative log-likelihood of the approximate density on the validation samples. We also use the validation samples to set the termination criteria during the optimization.

We follow the implementation of [START_REF] Rothfuss | Conditional Density Estimation with Neural Networks: Best Practices and Benchmarks[END_REF] to define the architectures of these models. The hyperparameters we consider for the neural networks in the MDN and NF models are: 2 hidden layers, 32 hidden units in each layer, {5, 10, 20, 50, 100} centers or flows, weight normalization, and a dropout probability of {0, 0.2} for regularizing the neural networks during training. For CKDE and NKDE we select the bandwidth of the kernel estimators using 5-fold cross-validation.