# Triangular Flows for Generative Modeling:
Statistical Consistency, Smoothness Classes,
and Fast Rates

Nicholas J. Irons  Department of Statistics, University of Washington  Meyer
Scetbon  Department of Statistics, CREST ENSAE  Soumik Pal  Department of
Mathematics, University of Washington  Zaid Harchaoui  Department of
Statistics, University of Washington

(December 31, 2021)

###### Abstract

Triangular flows, also known as Knöthe-Rosenblatt measure couplings, comprise
an important building block of normalizing flow models for generative modeling
and density estimation, including popular autoregressive flows such as Real
NVP. We present statistical guarantees and sample complexity bounds for
triangular flow statistical models. In particular, we establish the
statistical consistency and the finite sample convergence rates of the
Kullback-Leibler estimator of the Knöthe-Rosenblatt measure coupling using
tools from empirical process theory. Our results highlight the anisotropic
geometry of function classes at play in triangular flows, shed light on
optimal coordinate ordering, and lead to statistical guarantees for Jacobian
flows. We conduct numerical experiments on synthetic data to illustrate the
practical implications of our theoretical findings.

##  1 Introduction

Triangular flows are popular generative models that allow one to define
complex multivariate distributions via push-forwards from simpler multivariate
distributions [25]. Triangular flow models target the _Knöthe-Rosenblatt map_
[33], which originally appeared in two independent papers by Knöthe and
Rosenblatt [24, 30]. The Knöthe-Rosenblatt map is a multivariate function
$S^{*}$ from $\mathbb{R}^{d}$ onto itself pushing a Lebesgue probability
density $f$ onto another one $g$:

$S^{*}\\#f=g$

The Knöthe-Rosenblatt map has the striking property of being _triangular_ in
that its Jacobian is an upper triangular matrix. This map and its properties
have been studied in probability theory, nonparametric statistics, and optimal
transport, under the name of Knöthe-Rosenblatt (KR) coupling or rearrangement.

KR can also be seen as an _optimal transport under anisotropic distortion_
[6], hence shared connections with optimal transport theory [32, 35]. Unlike
the well-known Brenier map, KR has a simple and explicit definition in terms
of conditional densities of distributions at play, sparing the need for
optimization to be defined.

KR can be used to synthesize a sampler of a probability distribution given
data drawn from that probability distribution. Moreover, any probability
distribution can be well approximated via a KR map. This key property has been
an important motivation of a number of flow models [10, 11, 15, 23, 25, 27,
33, 19, 29, 37, 16] in machine learning, statistical science, computational
science, and AI domains. Flow models or normalizing flows have achieved great
success in a number of settings, allowing one to generate images that look
realistic as well as texts that look as if they were written by humans. The
theoretical analysis of normalizing flows is yet a huge undertaking as many
challenges arise at the same time: the recursive construction of a push
forward, the learning objective to estimate the push forward, the neural
network functional parameterization, and the statistical modeling of complex
data.

We focus in this paper on the Knöthe-Rosenblatt coupling, and more generally
triangular flows, as it can be seen as the statistical backbone of normalizing
flows. The authors of [33] showed how to estimate KR from data by minimizing a
Kullback-Leibler objective. The estimated KR can then be used to sample at
will from the probability distribution at hand. The learning objective of [33]
has the benefit of being statistically classical, hence amenable to detailed
analysis compared to adversarial learning objectives which are still subject
to active research. The sample complexity or rate of convergence of this KR
estimator is however, to this day, unknown. On the other hand, KR and its
multiple relatives are frequently motivated from a universal approximation
perspective [19, 25], which can be misleading as it focuses on the
approximation aspect only. Indeed, while universal approximation holds for
many proposed models, slow rates can occur (see Theorem 2.1).

##### Contributions.

We present a theoretical analysis of the Knöthe-Rosenblatt coupling, from its
statistical framing to convergence rates. We put forth a simple example of
slow rates showing the limitations of a viewpoint based on universal
approximation only. This leads us to identify the function classes that the KR
maps belong to and bring to light their anisotropic geometry. We then
establish finite sample rates of convergence using tools from empirical
process theory. Our analysis delineates different regimes of statistical
convergence, depending on the dimension and the sample. Our theoretical
results hold under general conditions. Assuming that the source density is
log-concave, we establish fast rates of Sobolev-type convergence in the smooth
regime. We outline direct implications of our results on Jacobian flows. We
provide numerical illustrations on synthetic data to demonstrate the practical
implications of our theoretical results.

##  2 Triangular Flows

##### Knöthe-Rosenblatt rearrangement.

KR originated from independent works of Rosenblatt and Knöthe and has spawned
fruitful applications in diverse areas. M. Rosenblatt studied the KR map for
statistical purposes, specifically multivariate goodness-of-fit testing [30].
H. Knöthe, on the other hand, elegantly utilized the KR map to extend the
Brunn-Minkowski inequality, which can be used to prove the celebrated
isoperimetric inequality [24]. More recently, triangular flows have been
proposed as simple and expressive building blocks of generative models, which
aim to model a distribution given samples from it. They have been studied and
implemented for the problems of sampling and density estimation, among others
[25, 33, 27, 15]. A triangular flow can be used to approximate the KR map
between a source density and a target density from their respective samples.

Consider two Lebesgue probability densities $f$ and $g$ on $\mathbb{R}^{d}$.
The Knöthe-Rosenblatt map $S^{*}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ between $f$
and $g$ is the unique monotone non-decreasing upper triangular measurable map
pushing $f$ forward to $g$, written $S^{*}\\#f=g$ [33]. The KR map is upper
triangular in the sense that its Jacobian is an upper triangular matrix, since
$S^{*}$ takes the form

$S^{*}(x)=\begin{bmatrix}S^{*}_{1}(x_{1},\ldots,x_{d})\\\
S^{*}_{2}(x_{2},\ldots,x_{d})\\\ \vdots\\\ S^{*}_{d-1}(x_{d-1},x_{d})\\\
S^{*}_{d}(x_{d})\end{bmatrix}.$

Furthermore, $S^{*}$ is monotone non-decreasing in the sense that the
univariate map $x_{k}\mapsto S^{*}_{k}(x_{k},\ldots,x_{d})$ is monotone non-
decreasing for any $(x_{k+1},\ldots,x_{d})\in\mathbb{R}^{d-k}$ for each
$k\in\\{1,\ldots,d\\}=:[d]$. We will sometimes abuse notation and write
$S^{*}_{k}(x)$ for $x=(x_{1},\ldots,x_{d})$, even though $S^{*}_{k}$ depends
only on $(x_{k},\ldots,x_{d})$.

The components of the KR map can be defined recursively via the monotone
transport between the univariate conditional densities of $f$ and $g$. Let
$F_{k}(x_{k}|x_{(k+1):d})$ denote the cdf of the conditional density

$f_{k}(x_{k}|x_{(k+1):d}):=\frac{\int f(x_{1},\ldots,x_{d})dx_{1}\cdots
dx_{k-1}}{\int f(x_{1},\ldots,x_{d})dx_{1}\cdots dx_{k}}.$

Similarly, let $G_{k}(y_{k}|y_{(k+1):d})$ denote the conditional cdf of $g$.
In particular, when $k=d$ we obtain $F_{d}(x_{d})$, the cdf of the $d$-th
marginal density $f_{d}(x_{d})=\int f(x_{1},\ldots,x_{d})dx_{1}\cdots
dx_{d-1}$, and similarly for $g$. Assuming $g(y)>0$ everywhere, the maps
$y_{k}\mapsto G_{k}(y_{k}|y_{(k+1):d})$ are strictly increasing and therefore
invertible. In this case, the $d$th component of $S^{*}$ is defined as the
monotone transport between $f_{d}(x_{d})$ and $g_{d}(y_{d})$ [32],

$S^{*}_{d}(x_{d})=G_{d}^{-1}(F_{d}(x_{d})).$

From here the $k$th component of $S^{*}$ is given by

$S^{*}_{k}(x_{k},\ldots,x_{d})=G_{k}^{-1}\left(F_{k}(x_{k}|x_{(k+1):d})\bigg{|}S^{*}_{(k+1):d}(x_{(k+1):d})\right)$

for $k=1,\ldots,d-1$, where $x_{(k+1):d}=(x_{k+1},\ldots,x_{d})$ and

$S^{*}_{(k+1):d}(x_{(k+1):d})=(S^{*}_{k+1}(x_{(k+1):d}),\ldots,S^{*}_{d}(x_{d})).$

That $S^{*}$ is upper triangular and monotonically non-decreasing is clear
from the construction. Under tame assumptions on $f$ and $g$ discussed below,
$S^{*}$ is invertible. We denote $T^{*}=(S^{*})^{-1}$, which is the KR map
from $g$ to $f$. In this paper, we shall be interested in the asymptotic and
non-asymptotic convergence of statistical estimators towards $S^{*}$ and
$T^{*}$, respectively.

In one dimension ($d=1$), the Knöthe-Rosenblatt rearrangement simply matches
the quantiles of $f$ and $g$ via monotone transport, $S^{*}(x)=(G^{-1}\circ
F)(x)$. In this case, the KR rearrangement is an optimal transport map for any
convex cost that is a function of the displacement $x-S(x)$ (see, e.g.,
Theorem 2.9 in [32]). In higher dimensions, the KR map is not generally an
optimal transport map for the commonly studied cost functions. However, it
does arise as a limit of optimal transport maps under an anisotropic quadratic
cost. We discuss this connection further in Section 4, in which we prove a
theorem concerning densities of anisotropic smoothness that complements the
result of Carlier, Galichon, and Santambrogio [6].

|  |  |
---|---|---|---
Figure 1: Left: Heat map of standard bivariate Gaussian density
$f(x_{1},x_{2})$ with correlation $\rho=0$ (leftmost) and the first component
of the KR map $F_{1}(x_{1}|x_{2})$ from $f$ to $g$, the uniform density on
$[0,1]^{2}$. Right: Bivariate Gaussian density with $\rho=0.7$ and first
component of the KR map $F_{1}(x_{1}|x_{2})$ (rightmost).

##### From the uniform distribution to any distribution.

In his seminal paper on the Knöthe-Rosenblatt rearrangement, Rosenblatt
considered the special case in which $g$ is the uniform density on the
hypercube $[0,1]^{d}$ [30]. This implies that the conditional cdfs are
identity maps, $G_{k}(y_{k}|y_{k+1},\ldots,y_{d})=y_{k}$, and so the KR map
from $f$ to $g$ consists simply of the conditional cdfs

$\displaystyle S^{*}_{d}(x_{d})$ $\displaystyle=P(X_{d}\leq
x_{d})=F_{d}(x_{d}),$ $\displaystyle S^{*}_{d-1}(x_{d-1},x_{d})$
$\displaystyle=P(X_{d-1}\leq x_{d-1}|X_{d}=x_{d})=F_{d-1}(x_{d-1}|x_{d}),$
$\displaystyle\qquad\vdots$ $\displaystyle S^{*}_{1}(x_{1},\ldots,x_{d})$
$\displaystyle=P(X_{1}\leq
x_{1}|X_{2}=x_{2},\ldots,X_{d}=x_{d})=F_{1}(x_{1}|x_{2},\ldots,x_{d}),$

When $d=2$ and $f$ is a bivariate Gaussian $N\left(\begin{bmatrix}\mu_{1}\\\
\mu_{2}\end{bmatrix},\begin{bmatrix}\sigma_{1}^{2}&\rho\sigma_{1}\sigma_{2}\\\
\rho\sigma_{1}\sigma_{2}&\sigma_{2}^{2}\end{bmatrix}\right)$, we have

$\displaystyle F_{2}(x_{2})$
$\displaystyle=\Phi\left(\frac{x_{2}-\mu_{2}}{\sigma_{2}}\right),$
$\displaystyle F_{1}(x_{1}|x_{2})$
$\displaystyle=\Phi\left(\frac{x_{1}-\mu_{1}-\frac{\rho\sigma_{2}}{\sigma_{1}}(x_{2}-\mu_{2})}{\sigma_{1}\sqrt{1-\rho^{2}}}\right),$

where $\Phi$ is the univariate standard normal cdf. Figure 1 exhibits heat
maps of the target density $f(x_{1},x_{2})$ and the first component of the KR
map $F_{1}(x_{1}|x_{2})$ for the choice
$\mu_{1}=\mu_{2}=0,\sigma_{1}=\sigma_{2}=1,$ with $\rho=0$ (left panels) and
$\rho=0.7$ (right panels).

##### Kullback-Leibler objective.

The authors of [33] proposed to estimate the KR map using a minimum distance
approach based on the Kullback-Leibler divergence. Since $S^{*}$ is the unique
monotone upper triangular map satisfying $S\\#f=g$, it is the unique such map
satisfying $\text{KL}(S\\#f|g)=0$, where KL denotes the Kullback-Leibler
divergence,

$\text{KL}(p|q)=\int p(x)\log\frac{p(x)}{q(x)}dx$

for Lebesgue densities $p,q$ on $\mathbb{R}^{d}$. As such, $S^{*}$ uniquely
solves the variational problem $\min_{S\in\mathcal{T}}\text{KL}(S\\#f|g),$
where $\mathcal{T}$ denotes the convex cone of monotone non-decreasing upper
triangular maps $S:\mathbb{R}^{d}\to\mathbb{R}^{d}$. The change of variables
formula for densities states that

$(S^{-1}\\#g)(x)=g(S(x))|\det(JS(x))|,$ (1)

where $JS(x)$ denotes the Jacobian matrix of $S$ evaluated at $x$. Applying
this formula to $\text{KL}(S\\#f|g)=\text{KL}(f|S^{-1}\\#g)$, we rewrite
$\min_{S\in\mathcal{T}}\text{KL}(S\\#f|g)$ as

$\min_{S\in\mathcal{T}}\mathbb{E}\left\\{\log[f(X)/g(S(X))]-\sum_{k=1}^{d}\log
D_{k}S_{k}(X)\right\\},$ (2)

where $X$ is a random variable on $\mathbb{R}^{d}$ with density $f$ and
$D_{k}$ is shorthand for differentiation with respect to the $k$th component
$\frac{\partial}{\partial x_{k}}$. By monotonicity, $D_{k}S_{k}(x)$ is defined
Lebesgue almost everywhere for every $S\in\mathcal{T}$. The relation (2) is
proved in the Supplement.

###  2.1 Statistical estimator of the KR map

We shall study an estimator $S^{n}$ of $S^{*}$ derived from the sample average
approximation to (2), which yields the minimization problem [33]

$\min_{S\in\mathcal{S}}\frac{1}{n}\sum_{i=1}^{n}\left\\{\log[f(X^{i})/g(S(X^{i}))]-\sum_{k=1}^{d}\log
D_{k}S_{k}(X^{i})\right\\},$ (3)

where $X^{1},\ldots,X^{n}$ is an i.i.d random sample from $f$ and
$\mathcal{S}$ is a hypothesis function class. In generative modeling, we have
a finite sample from $f$, perhaps an image dataset, that we use to train a map
that can generate more samples from $f$. In this case, $f$ is the target
density and $g$, the source density, is a degree of freedom in the problem. In
practice, $g$ should be chosen so that it is easy to sample from, e.g., a
multivariate normal density. The target density $f$ could be also unknown in
practice, and if necessary we can omit the terms involving $f$ from the
objective function in (3), since they do not depend on the argument $S$.

With an estimator $S^{n}$ in hand, which approximately solves the sample
average problem (3), we can generate approximate samples from $f$ by pulling
back samples from $g$ under $S^{n}$, or equivalently by pushing forward
samples from $g$ under $T^{n}=(S^{n})^{-1}$. As $S^{n}$ is defined via KL
projection, it can also be viewed as a non-parametric maximum likelihood
estimator (MLE).

In practice, $S^{n}$ can be estimated by parameterizing the space of
triangular maps via some basis expansion, for example, using orthogonal
polynomials [27]. In this case, the problem becomes parametric and $S^{n}$ is
the MLE. Universal approximation is insufficient to reason about nonparametric
estimators, since slow rates can happen, as we shall show next.

Figure 2 illustrates the problem at hand in the case where the target density
$f$ is an unbalanced mixture of three bivariate normal distributions centered
on the vertices of an equilateral triangle with spherical covariance. Panel
(a) of Figure 2 displays level curves of $f$. The source density $g$ is a
standard bivariate normal density. We solved for $S^{n}$ by parametrizing its
components via a Hermite polynomial basis, as described in [27]. Since $g$ is
log-concave, the optimization problem (3) is convex and can be minimized
efficiently using standard convex solvers. By the change of variables formula,
$f(x)=g(S^{*}(x))|\det(JS^{*}(x)|$, where $|\det(JS^{*}(x))|$ denotes the
Jacobian determinant of the KR map $S^{*}$ from $f$ to $g$. Hence we take
$f_{n}(x)=g(S^{n}(x))|\det(JS^{n}(x))|$ as an estimate of the target density.
Panels (b)-(d) of Figure 2 display the level curves of $f_{n}$ as the sample
size $n$ increases from $1500$ to $5000$. The improving accuracy of the
density estimates $f_{n}$ as the sample size grows is consistent with our
convergence results established in Section 3.

|  |  |
---|---|---|---
(a) Truth |  (b) $n=1500$ |  (c) $n=3000$ |  (d) $n=5000$
Figure 2: (a) Level curves of the normal mixture density
$f(x)=g(S^{*}(x))|\det(JS^{*}(x))|$; (b)-(d) Level curves of the density
estimates $f_{n}(x)=g(S^{n}(x))|\det(JS^{n}(x))|$ trained on
$n=1500,3000,5000$ samples, respectively.

##### Slow Rates.

Without combining both a tail condition (e.g., common compact support) and a
smooth regularity condition (e.g., uniformly bounded derivatives) on the
function class $\mathcal{F}$ of the target density $f$, we show that
convergence of any estimator $T^{n}$ of the direct map $T$ from $g$ to $f$ can
occur at an arbitrarily slow rate.

######  Theorem 2.1.

Let $\mathcal{F}$ denote the class of infinitely continuously differentiable
Lebesgue densities supported on the $d$-dimensional hypercube $[0,1]^{d}$ and
uniformly bounded by 2, i.e., $\sup_{f\in\mathcal{F}}\|f\|_{\infty}\leq 2$.
Let $g$ be any Lebesgue density on $\mathbb{R}^{d}$.

For any $n\in\mathbb{N}$, the minimax risk in terms of KL divergence is
bounded below as

$\inf_{T^{n}}\sup_{f\in\mathcal{F}}\mathbb{E}_{f}[\text{KL}(f|f_{n}))]\geq
1/2,$

where $T^{n}:\mathbb{R}^{d}\to[0,1]^{d}$ is any estimate of the KR map from
$g$ to $f$ based on an iid sample of size $n$ from $f$, and $f_{n}=T^{n}\\#g$
is the density estimate of $f$.

Theorem 2.1 underscores the importance of going beyond universal approximation
results to study the sample complexity and statistical performance of KR
estimation. The proof of this “no free lunch” theorem follows an idea of Birgé
[3]; see also [2, 5, 7, 8, 9, 18]. We construct a family of densities in
$\mathcal{F}$ built from rapidly oscillating perturbations of the uniform
distribution on $[0,1]^{d}$. Such densities are intuitively difficult to
estimate. As is evident from the construction, however, a suitable uniform
bound on the derivatives of the functions in $\mathcal{F}$ would preclude the
existence of such pathological examples. As such, in what proceeds we aim to
derive convergence rate bounds under the assumption that the target and source
densities $f$ and $g$, respectively, are compactly supported and sufficiently
regular, in the sense that they lie in a Sobolev space of functions with
continuous and uniformly bounded partial derivatives up to order $s$ for some
$s\geq 1$. For simplicity, our theoretical treatment assumes that $f$ and $g$
are fixed, but our convergence rate bounds as stated also bound the worst-case
KL risk over any $f$ and $g$ lying in the $L^{\infty}$ Sobolev ball
$\mathcal{F}=\left\\{h:\sum_{|\mathbf{\alpha}|\leq
s}\|D^{\mathbf{\alpha}}h\|_{\infty}\leq B\right\\}$ for any fixed $B>0$.

Theorem 2.1 provides a lower bound on the minimax KL risk of KR estimation
over the hypothesis class of target densities $\mathcal{F}$. The following
stronger result, based on work of Devroye [7, 8, 9], demonstrates that
convergence can still occur arbitrarily slowly for the task of estimating a
single target density.

######  Theorem 2.2.

Let $g$ be any Lebesgue density on $\mathbb{R}^{d}$ and
$\\{a_{n}\\}_{n=1}^{\infty}$ any sequence converging to zero with $1/512\geq
a_{1}\geq a_{2}\geq\cdots\geq 0.$ For every sequence of KR map estimates
$T^{n}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ based on a random sample of size $n$,
there exists a target distribution $f$ on $\mathbb{R}^{d}$ such that

$\mathbb{E}_{f}[\text{KL}(f|f_{n})]\geq a_{n},$

where $f_{n}=T^{n}\\#g$ is the density estimate of $f$.

##  3 Statistical Consistency

Setting the stage for the theoretical results, let us introduce the main
assumptions.

######  Assumption 3.1.

The Lebesgue densities $f$ and $g$ have convex compact supports
$\mathcal{X},\mathcal{Y}\subset\mathbb{R}^{d}$, respectively.

######  Assumption 3.2.

The densities $f,g$ are bounded away from 0, i.e.,
$\inf_{x\in\mathcal{X},y\in\mathcal{Y}}\\{f(x),g(y)\\}>0$

######  Assumption 3.3.

Let $s\geq 1$ be a positive integer. The densities $f$ and $g$ are $s$-smooth
on their supports, in the sense that

$D^{\alpha}f(x):=\frac{\partial^{|\alpha|}}{\partial
x_{1}^{\alpha_{1}}\cdots\partial x_{d}^{\alpha_{d}}}f(x)$

is continuous for every multi-index
$\mathbf{\alpha}=(\alpha_{1},\ldots,\alpha_{d})\in\mathbb{Z}^{d}_{+}$
satisfying $|\mathbf{\alpha}|:=\sum_{k=1}^{d}\alpha_{k}\leq s$ and similarly
for $g$.

It is well known that the KR map $S^{*}$ from $f$ to $g$ is as smooth as the
densities $f$ and $g$, but not smoother [32]. As such, under Assumptions
3.1-3.3 we can restrict our attention from $\mathcal{T}$, the set of monotone
non-decreasing upper triangular maps, to the smaller function class of
monotone upper triangular maps that are $s$-smooth, of which the KR map
$S^{*}$ is an element. That is, we can limit our search for an estimator
$S^{n}$ solving (3) to a space of functions with more structure. This
restriction is crucial to establishing a rate of convergence of the estimator
$S^{n}$, as we can quantitatively bound the complexity of spaces of smooth
maps. We discuss these developments in further detail below. Proofs of all
results are included in the Supplement.

###  3.1 Upper bounds on metric entropy

We first derive useful estimates of the metric entropy of function classes
previously introduced. Assumptions 3.1-3.3 allow us focus on smooth subsets of
the class of monotone upper triangular maps $\mathcal{T}$.

######  Definition 3.1.

Let $M>0$. For $\mathfrak{s}\in\mathbb{Z}^{d}_{+}$ and $k\in[d]$, let
$\tilde{\mathfrak{s}}_{k}=(s_{k}+1,s_{k+1},\ldots,s_{d})$. Define
$\mathcal{T}(\mathfrak{s},d,M)\subset\mathcal{T}$ as the convex subset of
strictly increasing upper triangular maps $S:\mathcal{X}\to\mathcal{Y}$
satisfying:

  1. 1.

$\inf_{k\in[d],x\in\mathcal{X}}D_{k}S_{k}(x)\geq 1/M$,

  2. 2.

$\|D^{\alpha}S_{k}\|_{\infty}\leq M$ for all $k\in[d]$ and
$\alpha_{k:d}\preceq\tilde{\mathfrak{s}}_{k}$.

For $s\in\mathbb{N}$, we also define the homogeneous smoothness class
$\mathcal{T}(s,d,M)=\mathcal{T}((s,\ldots,s),d,M)$.

Condition 1 guarantees that the Jacobian term in (3) is bounded, and condition
2 guarantees smoothness of the objective. In Section 4 below we consider
densities with anisotropic smoothness, in which case the number of continuous
derivatives $s_{k}$ varies with the coordinate $x_{k}$. For simplicity and
clarity of exposition, we first focus on the case when $f$ and $g$ are smooth
in a homogeneous sense, as in Assumption 3.3, and we work in the space
$\mathcal{T}(s,d,M)$. As remarked above, the KR map $S^{*}$ from $f$ to $g$
lies in $\mathcal{T}(s,d,M^{*})$ under Assumptions 3.1-3.3 when $M^{*}$ is
sufficiently large. The same is true of the direct map $T^{*}$ from $g$ to
$f$. In fact, all of the results stated here for the sampling map $S^{*}$ also
hold for the direct map $T^{*}$, possibly with minor changes (although the
proofs are generally more involved). For brevity, we mainly discuss $S^{*}$
and direct the interested reader to the Supplement. Henceforth, we consider
estimators $S^{n}$ lying in $\mathcal{T}(s,d,M^{*})$ that minimize the
objective in (3). We leave the issue of model selection, i.e., determining a
sufficiently large $M^{*}$ such that $\mathcal{T}(s,d,M^{*})$ contains the
true KR map $S^{*}$, for future work. In the Supplement, we calculate explicit
quantitative bounds on the complexity of this space as measured by the metric
entropy in the $d$-dimensional $L^{\infty}$ norm
$\|S\|_{\infty,d}:=\max_{k\in[d]}\|S_{k}\|_{\infty}$. The compactness of
$(\overline{\mathcal{T}(s,d,M^{*})},\|\cdot\|_{\infty,d})$ derived as a
corollary of this result is required to establish the convergence of a
sequence of estimators $S^{n}$ to $S^{*}$, and the entropy bound on the
corresponding class of Kullback-Leibler loss functions over
$\mathcal{T}(s,d,M^{*})$ in Proposition 3.3 below allows us to go further by
deriving bounds on the rate of convergence in KL divergence. This result
builds off known entropy estimates for function spaces of Besov type [2, 3,
28].

######  Definition 3.2.

For a map $S\in\mathcal{T}$, define the loss function
$\psi_{S}:\mathcal{X}\to\mathbb{R}$ by

$\displaystyle\psi_{S}(x)$ $\displaystyle=\log f(x)-\log
g(S(x))-\log|\det(JS(x))|$ $\displaystyle=\log f(x)-\log
g(S(x))-\sum_{k=1}^{d}\log D_{k}S_{k}(x).$

We also define the class $\Psi(s,d,M^{*})$ of loss functions over
$\mathcal{T}(s,d,M^{*})$ as

$\Psi(s,d,M^{*}):=\\{\psi_{S}:S\in\mathcal{T}(s,d,M^{*})\\}.$

By (2) we have $\mathbb{E}[\psi_{S}(X)]=\text{KL}(S\\#f|g)$, where $X\sim f$.
Similarly, the sample average of $\psi_{S}$ is the objective in (3). Hence, to
derive finite sample bounds on the expected KL loss, we must study the sample
complexity of the class $\Psi(s,d,M^{*})$. Define
$N(\epsilon,\Psi(s,d,M^{*}),\|\cdot\|_{\infty})$ as the $\epsilon$-covering
number of $\Psi(s,d,M^{*})$ with respect to the uniform norm
$\|\cdot\|_{\infty}$, and the metric entropy

$H(\epsilon,\Psi(s,d,M^{*}),\|\cdot\|_{\infty})=\log
N(\epsilon,\Psi(s,d,M^{*}),\|\cdot\|_{\infty}).$

######  Proposition 3.3.

Under Assumptions 3.1-3.3, the metric entropy of $\Psi(s,d,M^{*})$ in
$L^{\infty}(\mathcal{X})$ is bounded as

$H(\epsilon,\Psi(s,d,M^{*}),\|\cdot\|_{\infty})\lesssim\epsilon^{-d/s}.$

Consequently, $\Psi(s,d,M^{*})$ is totally bounded and therefore precompact in
$L^{\infty}(\mathcal{X})$.

Here, for functions $a(\epsilon),b(\epsilon)$ (or sequences $a_{n},b_{n}$) we
write $a(\epsilon)\lesssim b(\epsilon)$ (resp. $a_{n}\lesssim b_{n}$) if
$a(\epsilon)\leq Cb(\epsilon)$ (resp. $a_{n}\leq Cb_{n}$) for all $\epsilon$
(resp. $n$) for some constant $C>0$. For brevity, in this result and those
that follow, we suppress scalar prefactors that do not depend on the sample
size $n$. As our calculations in the Supplement demonstrate, the constant
prefactors in this and subsequent bounds are polynomial in the
$\|\cdot\|_{\infty}$ radius $M^{*}$ and exponential in the dimension $d$. This
dependence is similar to others in the literature concerning sample complexity
of transport map estimators [20].

###  3.2 Statistical consistency

For the sake of concision, we introduce empirical process notation [13, 14,
34, 36]. Let $\mathcal{H}$ be a collection of functions from
$\mathcal{X}\subseteq\mathbb{R}^{d}\to\mathbb{R}$ measurable and square
integrable with respect to $P$, a Borel probability measure on
$\mathbb{R}^{d}$. Let $P_{n}$ denote the empirical distribution of an iid
random sample $X_{1},\ldots,X_{n}$ drawn from $P$. For a function
$h\in\mathcal{H}$ we write $Ph:=\mathbb{E}[h(X)]$,
$P_{n}h:=\frac{1}{n}\sum_{i=1}^{n}h(X_{i})$, and
$\|P_{n}-P\|_{\mathcal{H}}:=\sup_{h\in\mathcal{H}}|(P_{n}-P)h|$. As discussed
in [34], measurability of the supremum process is a delicate subject. For
simplicity, we gloss over this technicality in our theoretical treatment, but
note that all results are valid as stated with outer expectation
$\mathbb{E}^{*}$ replacing $\mathbb{E}$ where relevant.

Let $P$ denote the probability measure with density $f$. With these new
definitions, the sample average minimization objective in (3) can be expressed
$P_{n}\psi_{S}$, while the population counterpart in (2) reads as
$P\psi_{S}=\text{KL}(S\\#f|g)$. Suppose the estimator $S^{n}$ of the sampling
map $S^{*}$ is a random element in $\mathcal{T}(s,d,M^{*})$ obtained as a
near-minimizer of $P_{n}\psi_{S}$. Let

$R_{n}=P_{n}\psi_{S^{n}}-\inf_{S\in\mathcal{T}(s,d,M)}P_{n}\psi_{S}\geq 0$

denote the approximation error of our optimization algorithm. Our goal is to
bound the loss $P\psi_{S^{n}}$. Fix $\epsilon>0$ and let $\tilde{S}$ be any
deterministic element of $\mathcal{T}(s,d,M^{*})$ that nearly minimizes
$P\psi_{S}$, i.e., suppose

$P\psi_{\tilde{S}}\leq\inf_{S\in\mathcal{T}(s,d,M^{*})}P\psi_{S}+\epsilon.$

It follows that

$\displaystyle P\psi_{S^{n}}-\inf_{S\in\mathcal{T}(s,d,M^{*})}P\psi_{S}$
$\displaystyle\leq P\psi_{S^{n}}-P\psi_{\tilde{S}}+\epsilon$
$\displaystyle=[P\psi_{S^{n}}-P\psi_{\tilde{S}}]+[P_{n}\psi_{\tilde{S}}-P_{n}\psi_{S^{n}}]+R_{n}+\epsilon$
$\displaystyle\leq 2\|P_{n}-P\|_{\Psi(s,d,M^{*})}+R_{n}+\epsilon.$

As $\epsilon>0$ was arbitrary, we conclude that

$P\psi_{S^{n}}-\inf_{S\in\mathcal{T}(s,d,M^{*})}P\psi_{S}\leq
2\|P_{n}-P\|_{\Psi(s,d,M^{*})}+R_{n}.$ (4)

Controlling the deviations of the empirical process
$\|P_{n}-P\|_{\Psi(s,d,M^{*})}$ as in Lemma 3.4 allows us to bound the loss of
the estimator $S^{n}$ and establish consistency and a rate of convergence in
KL divergence.

######  Lemma 3.4.

Under Assumptions 3.1-3.3, we have

$\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M^{*})}\lesssim\begin{cases}n^{-1/2},&d<2s,\\\
n^{-1/2}\log n,&d=2s,\\\ n^{-s/d},&d>2s.\end{cases}$

######  Remark 3.1.

Consider the case where $2s>d$, for example when both $f,~{}g$ are the
densities of the standard normal distribution, then we have by the central
limit theorem that for any $S\in\mathcal{T}(s,d,M^{*})$,
$\sqrt{n}\left[P_{n}\psi_{S}-P\psi_{S}\right]$ converge in law towards a
centered Gaussian distribution. Therefore we have that for any
$S\in\mathcal{T}(s,d,M^{*})$

$\displaystyle\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M^{*})}\geq\sup\limits_{S\in\mathcal{T}(s,d,M^{*})}\mathbb{E}\left[|P_{n}\psi_{S}-P\psi_{S}|\right]\gtrsim
n^{-1/2}$

which shows that our results are tight at least in the smooth regime.

The proof of Lemma 3.4 relies on metric entropy integral bounds established by
Dudley in [13] and van der Vaart and Wellner in [34]. Although we have phrased
the sample complexity bounds in Lemma 3.4 in terms of the expectation of the
empirical process $\|P_{n}-P\|_{\Psi(s,d,M)}$, high probability bounds can be
obtained similarly [36]. Hence, the following KL consistency theorem is
obtained as a direct result of Lemma 3.4 and the risk decomposition (4).

######  Theorem 3.5.

Suppose Assumptions 3.1-3.3 hold. Let $S^{n}$ be a near-optimizer of the
functional $S\mapsto P_{n}\psi_{S}$ on $\mathcal{T}(s,d,M^{*})$ with remainder
$R_{n}$ given by

$R_{n}=P_{n}\psi_{S^{n}}-\inf_{S\in\mathcal{T}(s,d,M^{*})}P_{n}\psi_{S}=o_{P}(1).$

Then $P\psi_{S^{n}}\stackrel{{\scriptstyle p}}{{\to}}P\psi_{S^{*}}=0$, i.e.,
$S^{n}$ is a consistent estimator of $S^{*}$ with respect to KL divergence.

Moreover, if $R_{n}$ is bounded in expectation as

$\mathbb{E}[R_{n}]\lesssim\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M^{*})},$

then the expected KL divergence of $S^{n}$ is bounded as

$\mathbb{E}[P\psi_{S^{n}}]\lesssim\begin{cases}n^{-1/2},&d<2s,\\\ n^{-1/2}\log
n,&d=2s,\\\ n^{-s/d},&d>2s.\end{cases}$

######  Remark 3.2.

Note that, as we work on a compact set, we also obtain the rates of
convergence of $S^{n}$ with respect to the Wasserstein metric:

$\mathbb{E}[W(S^{n}\\#f,g)]\lesssim\begin{cases}n^{-1/4},&d<2s,\\\
n^{-1/4}\log n,&d=2s,\\\ n^{-s/2d},&d>2s.\end{cases}$

where for any probability measures $\mu$ and $\nu$ on $\mathbb{R}^{d}$ with
finite first moments, we denote:

$\displaystyle W(\mu,\nu):=\sup\left\\{\left|\int hd\mu-\int
hd\nu\right|~{}:|h(x)-h(y)|\leq\|x-y\|\right\\}.$

###### Proof.

Let $\mu$ and $\nu$ be two probability measures on
$\mathcal{Z}:=\mathcal{X}\cup\mathcal{Y}\subset\mathbb{R}^{d}$. Then as both
$\mathcal{X}$ and $\mathcal{Y}$ are compact, we have the following well-known
inequalities [35]:

$\displaystyle
W(\mu,\nu)\leq\text{diam}(\mathcal{Z})\text{TV}(\mu,\nu)\leq\frac{\text{diam}(\mathcal{Z})}{\sqrt{2}}\sqrt{\text{KL}(\mu|\nu)}\;,$

then applying Jensen’s inequality to our results yields control in the
Wasserstein sense. ∎

##### Uniform convergence.

Although Theorem 3.5 only establishes a weak form of consistency in terms of
the KL divergence, we leverage this result to prove strong consistency, in the
sense of uniform convergence of $S^{n}$ to $S^{*}$ in probability, in Theorem
3.6. The proof requires understanding the regularity of the KL divergence with
respect to the topology induced by the $\|\cdot\|_{\infty,d}$ norm. In the
Supplement, we establish that KL is lower semicontinuous in
$\|\cdot\|_{\infty,d}$. Our result relies on the weak lower semicontinuity of
KL proved by Donsker and Varadhan using their well known dual representation
in Lemma 2.1 of [12].

######  Theorem 3.6.

Suppose Assumptions 3.1-3.3 hold. Let $S^{n}$ be any near-optimizer of the
functional $S\mapsto P_{n}\psi_{S}$ on $\mathcal{T}(s,d,M^{*})$, i.e., suppose

$P_{n}\psi_{S^{n}}=\inf_{S\in\mathcal{T}(s,d,M^{*})}P_{n}\psi_{S}+o_{P}(1).$

Then $\|S^{n}-S^{*}\|_{\infty,d}\stackrel{{\scriptstyle p}}{{\to}}0$, i.e,
$S^{n}$ is a consistent estimator of $S^{*}$ with respect to the uniform norm
$\|\cdot\|_{\infty,d}$ .

###### Proof.

Recall that $\overline{\mathcal{T}(s,d,M^{*})}$ is compact with respect to
$\|\cdot\|_{\infty,d}$. Together, lower semicontinuity of KL in
$\|\cdot\|_{\infty,d}$ and compactness guarantee that the KR map $S^{*}$,
which is the unique minimizer of $S\mapsto P\psi_{S}$ over
$\overline{\mathcal{T}(s,d,M^{*})}$, is well-separated in
$\mathcal{T}(s,d,M^{*})$. In other words, for any $\epsilon>0$,

$P\psi_{S^{*}}<\inf_{S\in\mathcal{T}(s,d,M^{*}):\|S-S^{*}\|_{\infty,d}\geq\epsilon}P\psi_{S}.$

Indeed, suppose to the contrary that we can find a deterministic sequence
$\tilde{S}^{n}\in\mathcal{T}(s,d,M^{*})$ satisfying
$\|\tilde{S}^{n}-S^{*}\|_{\infty,d}\geq\epsilon$ such that
$P\psi_{\tilde{S}^{n}}\to P\psi_{S^{*}}$. Since $\mathcal{T}(s,d,M^{*})$ is
precompact with respect to $\|\cdot\|_{\infty,d}$, we can extract a
subsequence $\tilde{S}^{n_{k}}$ converging to some
$\tilde{S}^{*}\in\overline{\mathcal{T}(s,d,M^{*})}$, which necessarily
satisfies $\|\tilde{S}^{*}-S^{*}\|_{\infty,d}\geq\epsilon$. By lower
semicontinuity of $S\mapsto P\psi_{S}$ with respect to $\|\cdot\|_{\infty,d}$,
it follows that

$P\psi_{\tilde{S}^{*}}\leq\liminf_{k\to\infty}P\psi_{\tilde{S}^{n_{k}}}=P\psi_{S^{*}}.$

We have a contradiction, since the KR map $S^{*}$ is the unique minimizer of
$S\mapsto P\psi_{S}$. This proves the claim that $S^{*}$ is a well-separated
minimizer.

Now fix $\epsilon>0$ and define

$\delta=\inf_{S\in\mathcal{T}(s,d,M^{*}):\|S-S^{*}\|_{\infty,d}\geq\epsilon}P\psi_{S}-P\psi_{S^{*}}>0.$

It follows that

$\\{\|S^{n}-S^{*}\|_{\infty,d}\geq\epsilon\\}\subseteq\\{P\psi_{S^{n}}-P\psi_{S^{*}}\geq\delta\\}.$

We have shown $P\psi_{S^{n}}\stackrel{{\scriptstyle p}}{{\to}}P\psi_{S^{*}}$
in Theorem 3.5. As a consequence,

$P(\|S^{n}-S^{*}\|_{\infty,d}\geq\epsilon)\leq
P(P\psi_{S^{n}}-P\psi_{S^{*}}\geq\delta)\to 0.$

As $\epsilon>0$ was arbitrary, we have
$\|S^{n}-S^{*}\|_{\infty,d}\stackrel{{\scriptstyle p}}{{\to}}0$. ∎

##### Inverse consistency.

We have proved consistency of the estimator $S^{n}$ of the sampling map
$S^{*}$, pushing forward the target $f$ to the source $g$. We can also get the
consistency and an identical rate of convergence of $T^{n}=(S^{n})^{-1}$
estimating $T^{*}=(S^{*})^{-1}$, although the proof of the analog to Theorem
3.6 establishing uniform consistency of $T^{n}$ is much more involved. We
defer to the Supplement for details.

##  4 Log-Concavity, Dimension Ordering, Jacobian Flows

##### Sobolev-type rates under log-concavity.

Suppose the source density $g$ is log-concave. Then
$\min_{S\in\mathcal{T}(s,d,M)}P_{n}\psi_{S}$ is a convex problem; moreover if
$g$ is strongly log-concave, strong convexity follows. The user can choose a
convenient $g$, such as a multivariate Gaussian with support truncated to a
compact convex set. In this case, we can establish a bound on the rate of
convergence of $S^{n}$ to $S^{*}$ in the $L^{2}$ Sobolev-type norm

$\|S\|^{2}_{H_{f}^{1,2}(\mathcal{X})}:=\sum_{k=1}^{d}\left\\{\|S_{k}\|^{2}_{L^{2}_{f}(\mathcal{X})}+\|D_{k}S_{k}\|^{2}_{L^{2}_{f}(\mathcal{X})}\right\\}.$

Here $\|\cdot\|_{L^{2}_{f}(\mathcal{X})}$ denotes the usual $L^{2}$ norm
integrating against the target density $f$.

######  Theorem 4.1.

Suppose Assumptions 3.1-3.3 hold. Assume further that $g$ is $m$-strongly log-
concave. Let $S^{n}$ be a near-optimizer of the functional $S\mapsto
P_{n}\psi_{S}$ on $\mathcal{T}(s,d,M^{*})$ with remainder $R_{n}$ satisfying

$\mathbb{E}[R_{n}]=\mathbb{E}\left\\{P_{n}\psi_{S^{n}}-\inf_{S\in\mathcal{T}(s,d,M^{*})}P_{n}\psi_{S}\right\\}\lesssim\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M^{*})}.$

Then $S^{n}$ converges to the true sampling map $S^{*}$ with respect to the
norm $H^{1,2}_{f}(\mathcal{X})$ norm with rate

$\mathbb{E}\|S^{n}-S^{*}\|_{H^{1,2}_{f}(\mathcal{X})}^{2}\lesssim\begin{cases}n^{-1/2},&d<2s,\\\
n^{-1/2}\log n,&d=2s,\\\ n^{-s/d},&d>2s,\end{cases}$

With more work, we can establish Sobolev convergence rates of the same order
in $n$ for $T^{n}=(S^{n})^{-1}$ to $T^{*}=(S^{*})^{-1}$, yet now in the
appropriate norm $\|\cdot\|_{H^{1,2}_{g}(\mathcal{Y})}$; see details in the
Supplement.

######  Remark 4.1.

Note that here we obtain the rates with respect to the
$H^{1,2}_{f}(\mathcal{X})$ norm from which we deduce immediately similar rates
for the $H^{1,p}_{f}(\mathcal{X})$ norm with $p\geq 2$. In addition, for fixed
$p$, under Assumptions 3.1-3.3, the $H^{1,p}_{f}(\mathcal{X})$ norm is
equivalent to $H^{1,p}(\mathcal{X})$ norm with the Lebesgue measure as the
reference measure. Finally, here we show the rates of the strong consistency
of our estimator with respect $H^{1,\cdot}_{f}(\mathcal{X})$ norm as we do not
require any additional assumptions on the higher order derivatives of the maps
living in $\mathcal{T}(s,d,M^{*})$. Additional assumptions on these
derivatives may lead to the convergence rates of higher order derivatives of
our estimate with respect to the $H^{k,p}_{f}(\mathcal{X}),k>1$ norm which is
beyond the scope of the present paper.

##### Dimension ordering.

Suppose now that the smoothness of the target density $f$ is anisotropic. That
is, assume $f(x_{1},\ldots,x_{d})$ is $s_{k}$-smooth in $x_{k}$ for each
$k\in[d]$. As there are $d!$ possible ways to order the coordinates, the
question arises: how should we arrange $(x_{1},\ldots,x_{d})$ such that the
estimator $S^{n}$ converges to the true KR map $S^{*}$ at the fastest possible
rate? Papamakarios et al. (2017) provide a discussion of this issue in the
context of autoregressive flows in Section 2.1 of [29]; they construct a
simple 2D example in which the model fails to learn the target density if the
wrong order of the variables is chosen.

This relates to choices made in neural architectures for normalizing flows on
images and texts. Our results suggest here that one would rather start with
the coordinates (i.e. data parts) that are the least smooth and make their way
through the triangular construction to the most smooth ones.

We formalize the anisotropic smoothness of the KR map as follows.

######  Assumption 4.1.

Let $\mathfrak{s}=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}_{+}$ be a multi-index
with $s_{k}\geq 1$ for all $k\in[d]$. The density $f$ is $\mathfrak{s}$-smooth
on its support, in the sense that

$D^{\alpha}f(x):=\frac{\partial^{|\alpha|}}{\partial
x_{1}^{\alpha_{1}}\cdots\partial x_{d}^{\alpha_{d}}}f(x)$

exists and is continuous for every multi-index
$\alpha=(\alpha_{1},\ldots,\alpha_{d})\in\mathbb{Z}^{d}_{+}$ satisfying
$\alpha\preceq\mathfrak{s}$, i.e., $\alpha_{k}\leq s_{k}$ for every $k\in[d]$.
Furthermore, we assume that $g(y)$ is
$(\|\mathfrak{s}\|_{\infty},\ldots,\|\mathfrak{s}\|_{\infty})$-smooth with
respect to $y=(y_{1},\ldots,y_{d})$ on $\mathcal{Y}$.

As the source density $g$ is a degree of freedom in the problem, we are free
to impose this assumption on $g$. Note that
$(\|\mathfrak{s}\|_{\infty},\ldots,\|\mathfrak{s}\|_{\infty})$-smoothness of
$g$ as defined in Assumption 4.1 is equivalent to
$\|\mathfrak{s}\|_{\infty}$-smoothness of $g$ as defined in Assumption 3.3.
The results that follow are slight variations on those in Sections 3.1 and 3.2
adapted to the anisotropic smoothness of the densities posited in Assumption
4.1.

Under Assumptions 3.1, 3.2, and 4.1, there exists some $M^{*}>0$ such that the
KR map $S^{*}$ from $f$ to $g$ lies in $\mathcal{T}(\mathfrak{s},d,M^{*})$;
see the Supplement for a proof. We also define the class of loss functions

$\Psi(\mathfrak{s},d,M^{*})=\\{\psi_{S}:S\in\mathcal{T}(\mathfrak{s},d,M^{*})\\},$

which appear in the objective $P_{n}\psi_{S}$. Hence, we can proceed as above
to bound the metric entropy and obtain uniform convergence and Sobolev-type
rates for estimators in the function class
$\mathcal{T}(\mathfrak{s},d,M^{*})$. Appealing to metric entropy bounds for
anisotropic smoothness classes [3, Proposition 2.2], we have the following
analog of Lemma 3.4.

######  Lemma 4.2.

For $k\in[d]$, let $d_{k}=d-k+1$ and
$\sigma_{k}=d_{k}\left(\sum_{j=k}^{d}s_{j}^{-1}\right)^{-1}$. Under
Assumptions 3.1, 3.2, and 4.1, we have

$\mathbb{E}\|P_{n}-P\|_{\Psi(\mathfrak{s},d,M^{*})}\lesssim\sum_{k=1}^{d}c_{n,k}.$

where we define

$c_{n,k}=\begin{cases}n^{-1/2},&d_{k}<2\sigma_{k},\\\ n^{-1/2}\log
n,&d_{k}=2\sigma_{k},\\\ n^{-\sigma_{k}/d_{k}},&d_{k}>2\sigma_{k}.\end{cases}$

Lemma 4.2 is proved via a chain rule decomposition of relative entropy which
relies upon the triangularity of the hypothesis maps
$\mathcal{T}(\mathfrak{s},d,M^{*})$. From here we can repeat the analysis in
Section 3.2 to obtain consistency and bounds on the rate of convergence of the
estimators $S^{n}$ and $T^{n}=(S^{n})^{-1}$ of the sampling map $S^{*}$ and
the direct map $T^{*}=(S^{*})^{-1}$, respectively, in the anisotropic
smoothness setting of Assumption 4.1. All the results in these Sections are
true with $\sum_{k}c_{n,k}$ replacing the rate under isotropic smoothness.

In order to minimize this bound to obtain an optimal rate of convergence, we
should order the coordinates $(x_{1},\ldots,x_{d})$ such that $\sigma_{k}$ is
as large as possible for each $k\in[d]$. Inspecting the definition of
$\sigma_{k}$ in Lemma 4.2, we see that this occurs when $s_{1}\leq\cdots\leq
s_{d}.$

######  Theorem 4.3.

The bound on the rate of convergence $\sum_{k}c_{n,k}$ is minimized when
$s_{1}\leq\cdots\leq s_{d}$, i.e., when the smoothness of the target density
$f$ in the direction $x_{k}$ increases with $1\leq k\leq d$.

Our result on the optimal ordering of coordinates complements the following
theorem of Carlier, Galichon, and Santambrogio adapted to our setup.

######  Theorem 4.4 (Theorem 2.1; Carlier, Galichon, and Santambrogio (2008)
[6]).

Let $f$ and $g$ be compactly supported Lebesgue densities on $\mathbb{R}^{d}$.
Let $\epsilon>0$ and let $\gamma^{\epsilon}$ be an optimal transport plan
between $f$ and $g$ for the cost

$c_{\epsilon}(x,y)=\sum_{k=1}^{d}\lambda_{k}(\epsilon)(x_{k}-y_{k})^{2},$

for some weights $\lambda_{k}(\epsilon)>0$. Suppose that for all
$k\in\\{1,\ldots,d-1\\}$, $\lambda_{k}(\epsilon)/\lambda_{k+1}(\epsilon)\to 0$
as $\epsilon\to 0$. Let $S^{*}$ be the Knöthe-Rosenblatt map between $f$ and
$g$ and $\gamma^{*}=(\text{id}\times S^{*})\\#f$ the associated transport
plan. Then $\gamma^{\epsilon}\rightsquigarrow\gamma^{*}$ as $\epsilon\to 0$.
Moreover, should the plans $\gamma^{\epsilon}$ be induced by transport maps
$S^{\epsilon}$, then these maps would converge to $S^{*}$ in $L^{2}(f)$ as
$\epsilon\to 0$.

With this theorem in mind, the KR map $S^{*}$ can be viewed as a limit of
optimal transport maps $S^{\epsilon}$ for which transport in the $d$th
direction is more costly than in the $(d-1)$th, and so on. The anisotropic
cost function $c_{\epsilon}(x,y)$ inherently promotes increasing regularity of
$S^{\epsilon}$ in $x_{k}$ for larger $k\in[d]$. Theorem 4.3 establishes the
same heuristic for learning triangular flows based on Knöthe-Rosenblatt
rearrangement to build generative models.

In particular, our result suggests that we should order the coordinates such
that $f$ is smoothest in the $x_{d}$ direction. This makes sense intuitively
because the estimates of the component maps $S^{*}_{k},k\in[d]$ all depend on
the $d$th coordinate of the data. As such, we should leverage smoothness in
$x_{d}$ to stabilize all of our estimates as much as possible. If not, we risk
propagating error throughout the estimation. In comparison, only the estimate
of the first component $S^{*}_{1}$ depends on $x_{1}$. Since our estimator
depends comparatively little on $x_{1}$, we should order the coordinates so
that $f$ is least smooth in the $x_{1}$ direction.

##### Jacobian Flows.

Suppose now we solve

$\inf_{S}\text{KL}(S\\#f|g)=\inf_{S}P\psi_{S}$

where the candidate map $S$ is a composition of smooth monotone increasing
upper triangular maps $U^{j}$ and orthogonal linear transformations
$\Sigma^{j}$ for $j\in[m]$, i.e.,

$S(x)=U^{m}\circ\Sigma^{m}\circ\cdots\circ U^{1}\circ\Sigma^{1}(x).$ (5)

We call $S$ a Jacobian flow of order $m$. This model captures many popular
normalizing autoregressive flows [25], such as Real NVP [11], in which the
$\Sigma^{j}$ are “masking” permutation matrices.

For simplicity, in this section we assume that $f$ and $g$ are supported on
the unit ball $B_{d}(0,1)\subset\mathbb{R}^{d}$ centered at the origin. Since
$\Sigma^{j}$ are orthogonal, we have $\Sigma^{j}(B_{d}(0,1))=B_{d}(0,1)$.
Hence, to accommodate the setup of the preceding sections, we can guarantee
that $S$ maps from $\mathcal{X}=B_{d}(0,1)$ to $\mathcal{Y}=B_{d}(0,1)$ by
requiring $\sup_{x\in B_{d}(0,1)}\|U^{j}(x)\|_{2}\leq 1$ for $j\in[m]$.

######  Definition 4.5.

Define the class $\mathcal{J}_{m}(\Sigma,s,M)$ of $s$-smooth Jacobian flows of
order $m$ to consist of those maps $S$ of the form (5) such that

  1. 1.

$\Sigma=(\Sigma^{1},\ldots,\Sigma^{m})$ are fixed orthogonal matrices,

  2. 2.

$\sup_{x\in B_{d}(0,1)}\|U^{j}(x)\|_{2}\leq 1$ for $j\in[m]$,

  3. 3.

$U^{j}\in\mathcal{T}(s,d,M)$ for $j\in[m]$.

We also define
$\Psi_{m}(\Sigma,s,M)=\\{\psi_{S}:S\in\mathcal{J}_{m}(\Sigma,s,M)\\}$.

By expanding our search to $\mathcal{J}_{m}(\Sigma,s,M)$, we are not targeting
the KR map $S^{*}$ and we are no longer guaranteed uniqueness of a minimizer
of the KL. Nevertheless, we can study the performance of estimators
$S^{n}\in\mathcal{J}_{m}(\Sigma,s,M)$ as compared to minimizers of KL in
$\overline{\mathcal{J}_{m}(\Sigma,s,M)}$, which are guaranteed to exist by
compactness of $\overline{\mathcal{J}_{m}(\Sigma,s,M)}$ and lower
semicontinuity of KL in $\|\cdot\|_{\infty,d}$.

Since the $\Sigma^{j}$ are orthogonal, we have
$|\det(J\Sigma^{j}(x))|=|\det(\Sigma^{j})|=1$, and therefore

$\displaystyle\psi_{S}(x)=$
$\displaystyle\log[f(x)/g(S(x))]-\sum_{j=1}^{m}\sum_{k=1}^{d}\log
D_{k}U_{k}^{j}(x^{j}),$

where we define

$x^{j}=\Sigma^{j}\circ U^{j-1}\circ\Sigma^{j-1}\circ\cdots\circ
U^{1}\circ\Sigma^{1}(x),\quad j\in[m].$

Hence, with a loss decomposition mirroring that in Definition 3.2, we can
apply the methods of the preceding sections to establish quantitative limits
on the loss incurred by the estimates $S^{n}$ in finite samples. See the
Supplement for details.

######  Theorem 4.6.

Suppose $f,g$ are $s$-smooth and supported on $B_{d}(0,1)$. Let $S^{n}$ be a
near-optimizer of the functional $S\mapsto P_{n}\psi_{S}$ on
$\mathcal{J}_{m}(\Sigma,s,M)$ with

$R_{n}=P_{n}\psi_{S^{n}}-\inf_{S\in\mathcal{J}_{m}(\Sigma,s,M)}P_{n}\psi_{S}=o_{P}(1).$

Further, let $S^{0}$ be any minimizer of $S\mapsto P\psi_{S}$ on
$\overline{\mathcal{J}_{m}(\Sigma,s,M)}$. It follows that

$P\psi_{S^{n}}\stackrel{{\scriptstyle p}}{{\to}}P\psi_{S^{0}}.$

Moreover, if $R_{n}$ is bounded in expectation as

$\mathbb{E}[R_{n}]\lesssim\mathbb{E}\|P_{n}-P\|_{\Psi_{m}(\Sigma,s,M)},$

then the expected KL divergence of $S^{n}$ is bounded as

$\mathbb{E}[P\psi_{S^{n}}]-P\psi_{S^{0}}\lesssim\begin{cases}n^{-1/2},&d<2s,\\\
n^{-1/2}\log n,&d=2s,\\\ n^{-s/d},&d>2s.\end{cases}$

##  5 Numerical Illustrations

|  |  |
---|---|---|---
|  |  |
|  |  |
|  |  |
Figure 3: Experimental comparison of convergence rates on 2D examples (top:
density; bottom: rates) for the 8 datasets considered in [37, 17].

We conducted numerical experiments to illustrate our theoretical results. To
estimate the KR map, we used Unconstrained Monotonic Neural Networks Masked
Autoregressive Flows (UMNN-MAF), a particular triangular flow introduced in
[37], with code to implement the model provided therein. UMNN-MAF learns an
invertible monotone triangular map targeting the KR rearrangement via KL
minimization

$S(x;\theta)=\begin{bmatrix}S_{1}(x_{1},\ldots,x_{d};\theta)\\\
S_{1}(x_{2},\ldots,x_{d};\theta)\\\ \vdots\\\ S_{d-1}(x_{d-1},x_{d};\theta)\\\
S_{d}(x_{d};\theta)\end{bmatrix},$

where each component $S_{k}(x_{k:d};\theta)$ is parametrized as the output of
a neural network architecture that can learn arbitrary monotonic functions.
Specifically, we have

$S_{k}(x_{k:d},\theta)=\int_{0}^{x_{k}}f_{k}(t,h_{k}(x_{(k+1):d};\phi_{k});\psi_{k})dt+\beta_{k}(h_{k}(x_{(k+1):d};\phi_{k});\psi_{k}),$
(6)

where $h_{k}(\cdot;\phi_{k}):\mathbb{R}^{d-k-1}\to\mathbb{R}^{q}$ is a
$q$-dimensional neural embedding of $x_{(k+1):d}$ and
$\beta_{k}(\cdot;\psi_{k}):\mathbb{R}^{q}\to\mathbb{R}$ is parametrized by a
neural network. Each $f_{k}$ is a strictly positive function parametrized by a
neural network, which guarantees that $\int_{0}^{x_{k}}f_{k}$ is increasing in
$x_{k}$. Here the total parameter $\theta$ is defined as
$\theta=\\{(\psi_{k},\phi_{k})\\}_{k=1}^{d}$. Further details of the model are
provided in [37]. We note that UMNN-MAF is captured by the model setup in our
theoretical treatment of KR map estimation.

The model was trained via log-likelihood maximization using minibatch gradient
descent with the Adam optimizer [22] with minibatch size 64, learning rate
$10^{-4}$, and weight decay $10^{-5}$. The integrand network architecture
defined in equation (6) consisted of 4 hidden layers of width 100. Following
Wehenkel and Louppe [37], the architecture of the embedding networks is the
best performing MADE network [16] used in NAF [19]. We used 20 integration
steps to numerically approximate the integral in equation (6). The source
density $g$ is a bivariate standard normal distribution. The population
negative log-likelihood loss, which differs from the KL objective by a
constant factor (namely, the negative entropy of the target density $f$), was
approximated by the empirical negative log-likelihood on a large independently
generated test set of size $N=10^{5}$.

Figure 3 exhibits our results for UMNN-MAF trained on 8 two-dimensional
datasets considered in [37] and [17]. Heatmaps of the target densities are
displayed in the top rows, while the bottom rows show the log-likelihood
convergence rates as the sample size increases from $n=10^{3}$ to $n=10^{4}$
on a log-log scale. For each training sample size, we repeated the experiment
100 times. We report the mean of the loss over the 100 replicates with 95%
error bars.

These experiments highlight the impact of the ordering of coordinates on the
convergence rate, as predicted by Theorem 4.3. The blue curves correspond to
first estimating the KR map along the horizontal $x_{1}$ axis, then the
vertical axis conditional on the horizontal, $x_{2}|x_{1}$. The orange curves
show the reverse order, namely estimating the KR map along the vertical
$x_{2}$ axis first. The 5 densities in the top rows (2 spirals, pinwheel,
moons, and banana) and the bottom right (swiss roll) are asymmetric in
$(x_{1},x_{2})$ (i.e., $f(x_{1},x_{2})$ is not exchangeable in
$(x_{1},x_{2})$). These densities exhibit different convergence rates
depending on the choice of order. The remaining 3 densities in the bottom rows
(8 gaussians, 2 circles, checkerboard) are symmetric in $(x_{1},x_{2})$ and do
not exhibit this behavior, which is to be expected. We also note that a linear
trend (on the log-log scale) is visible in the plots for smaller $n$, which
aligns with the convergence rates established in Theorem 3.5. For larger $n$,
however, approximation error dominates and the loss plateaus.

As an illustrative example, we focus in on the top right panel of Figure 3,
which plots the banana density $f(x_{1},x_{2})$ corresponding to the random
variables

$\displaystyle X_{2}$ $\displaystyle\sim N(0,1)$ $\displaystyle X_{1}|X_{2}$
$\displaystyle\sim N(X_{2}^{2}/2,1/2).$

It follows that $f$ is given by

$\displaystyle f(x_{1},x_{2})$ $\displaystyle=f(x_{1}|x_{2})f(x_{2})$
$\displaystyle\propto\exp\left\\{-\left(x_{1}-x_{2}^{2}/2\right)^{2}\right\\}\cdot\exp\left\\{-x_{2}^{2}/2\right\\}.$

Intuitively, estimating the normal conditional
$f(x_{1}|x_{2})=N(x_{2}^{2}/2,1/2)$ and the standard normal marginal
$f(x_{2})$ should be easier than estimating $f(x_{2}|x_{1})$ and $f(x_{1})$.
Indeed, as $x_{1}$ increases, we see that $f(x_{2}|x_{1})$ transitions from a
unimodal to a bimodal distribution. As such, we expect that estimating the KR
map $S^{(21)}$ from $f(x_{2},x_{1})$ to the source density $g$ should be more
difficult than estimating the KR map $S^{(12)}$ from $f(x_{1},x_{2})$ to $g$.
This is because the first component of $S^{(21)}$ targets the conditional
distribution $f(x_{2}|x_{1})$ and the second component targets $f(x_{1})$,
while the first component of $S^{(12)}$ targets $f(x_{1}|x_{2})$ and the
second component targets $f(x_{2})$. Indeed, this is what we see in the top
right panels of Figure 3, which shows the results of fitting UMNN-MAF to
estimate $S^{(12)}$ (orange) and $S^{(21)}$ (blue). As expected, we see that
estimates of $S^{(21)}$ converge more slowly than those of $S^{(12)}$. These
results are consistent with the findings of Papamakarios et al. (2017) [29],
who observed this behavior in estimating the banana density with MADE [16].

In the first 3 panels of Figure 4, we repeat the experimental setup on the 3
normal mixture densities considered by Kong and Chaudhuri in [26]. The
conclusions drawn from Figure 3 are echoed here. We observe no dependence in
the convergence rates on the choice of variable ordering, since the target
densities $f(x_{1},x_{2})$ are exchangeable in $(x_{1},x_{2})$. Furthermore,
we observe a linear convergence rate as predicted by Theorem 3.5.

|  |  |
---|---|---|---
|  |  |
Figure 4: Experimental comparison of convergence rates on 2D examples (top:
density; bottom: rates) for the 3 datasets considered in [26] and the sine
density with varying frequency $k_{2}\in\\{3,5,7\\}$.

Inspired by the pathological densities constructed in the proof of the “no
free lunch” Theorem 2.1, which are rapidly oscillating perturbations of the
uniform density on the hypercube, we now consider the sine density on the
hypercube, defined as

$f(x_{1},\ldots,x_{d})=1+\prod_{j=1}^{d}\sin(2\pi k_{j}x_{j})$

where $k_{j}\in\mathbb{Z}$ for $j\in[d]$, and
$x=(x_{1},\ldots,x_{d})\in[0,1]^{d}$. The smoothness of $f$, as measured by
any $L^{p}$ norm of its derivative(s), decreases as the frequency $|k_{j}|$
increases. As such, $f$ parametrizes a natural family of functions to test our
theoretical results concerning the statistical performance of KR map
estimation as a function of the sample size $n$, the smoothness of the
underlying target density, and the order of coordinates. The rightmost panels
of Figure 4 plot the sine density with $k_{1}=1,k_{2}=3$ (top row) and
convergence rates for the choices $k_{1}=1$ and $k_{2}\in\\{3,5,7\\}$ (bottom
row). The dashed lines correspond to estimating the marginal $x_{1}$ first,
followed by $x_{2}|x_{1}$; the solid lines indicate the reverse order. We
again see an effect of coordinate ordering on convergence rates. It is also
apparent that convergence slows down as $k_{2}$ increases and $f$ becomes less
smooth.

##  6 Discussion

##### Related work.

Previous work on KR couplings has focused on existence and approximation
questions in relation to universal approximation [4, 1, 19]. This universal
approximation property of KR maps has been a frequent motivation from a
theoretical viewpoint of normalizing flows; see Sec. 3 in [25]. KR maps have
been used for various learning and inference problems [25, 33, 27, 15].

Kong and Chaudhuri [26] study the expressiveness of basic types of normalizing
flows, such as planar flows, Sylvester flows, and Householder flows. In the
one-dimensional setting, they show that such flows are universal
approximators. However, when the distributions lives in a $d$-dimensional
space with $d\geq 2$, the authors provide a partially negative answer to the
universal approximation power of these flows. For example, they exhibit cases
where Sylvester flows cannot recover the target distributions. Their results
can be seen as complementary to ours as we give examples of arbitrary slow
statistical rates and we develop the consistency theory of KR-type flows [33].

Jaini et al. [21] investigate the properties of the increasing triangular map
required to push a tractable source density with known tails onto a desired
target density. Then they consider the general $d$-dimensional case and show
similarly that imposing smoothness condition on the increasing triangular map
will result in a target density with the same tail properties as the source.
Such results suggest that without any assumption on the target distribution,
the transport map might be too irregular to be estimated. These results echo
our assumptions on the target to obtain fast rates and complement ours by
focusing on the tail behavior while we focus on the consistency and the rates.

##### Conclusion.

We have established the uniform consistency and convergence rates of
statistical estimators of the Knöthe-Rosenblatt rearrangement, highlighting
the anisotropic geometry of function classes at play in triangular flows. Our
results shed light on coordinate ordering and lead to statistical guarantees
for Jacobian flows.

##### Acknowledgements.

The authors would like to thank the authors of [26] for sharing their code to
reproduce their results. N. Irons is supported by a Shanahan Endowment
Fellowship and a Eunice Kennedy Shriver National Institute of Child Health and
Human Development training grant, T32 HD101442-01, to the Center for Studies
in Demography & Ecology at the University of Washington, and NSF CCF-2019844.
M. Scetbon is supported by "Chaire d’excellence de l’IDEX Paris Saclay". S.
Pal is supported by NSF DMS-2052239 and a PIMS CRG (PIHOT). Z. Harchaoui is
supported by NSF CCF-2019844, NSF DMS-2134012, NSF DMS-2023166, CIFAR-LMB, and
faculty research awards. Part of this work was done while Z. Harchaoui was
visiting the Simons Institute for the Theory of Computing. The authors would
like to thank the Kantorovich Initiative of the Pacific Institute for the
Mathematical Sciences (PIMS) for supporting this collaboration. This work was
first presented at the Joint Statistical Meetings in August 2021.

## References

  * [1] D. E. Alexandrova.  Convergence of triangular transformations of measures.  Theory of Probability & Its Applications, 50(1):113–118, 2006\.
  * [2] L. Birgé.  Approximation dans les espaces métriques et théorie de l’estimation.  Z. Wahrscheinlichkeitstheor. Verw. Geb., 65:181–237, 1983.
  * [3] L. Birgé.  On estimating a density using Hellinger distance and some other strange facts.  Probab. Th. Rel. Fields, 71:271–291, 1986.
  * [4] V. I. Bogachev, A. V. Kolesnikov, and K. V. Medvedev.  Triangular transformations of measures.  Sbornik: Mathematics, 196(3):309, 2005.
  * [5] O. Bousquet, S. Boucheron, and G. Lugosi.  Introduction to statistical learning theory.  In O. Bousquet, U. von Luxburg, and G. Rätsch, editors, Advanced Lectures on Machine Learning, pages 169–207. Springer, Berlin, Heidelberg, 2004.
  * [6] G. Carlier, A. Galichon, and F. Santambrogio.  From Knöthe’s transport to Brenier’s map and a continuation method for optimal transport.  SIAM Journal on Mathematical Analysis, 41(6):2554–2576, 2010.
  * [7] L. Devroye.  On arbitrarily slow rates of global convergence in density estimation.  Z. Wahrscheinlichkeitstheorie verw. Gebiete, 62:475–483, 1983\.
  * [8] L. Devroye.  Another proof of a slow convergence result of Birgé.  Statistics and Probability Letters, 23:63–67, 1995.
  * [9] L. Devroye, L. Györfi, and G. Lugosi.  A Probabilistic Theory of Pattern Recognition.  Springer, 1996.
  * [10] L. Dinh, D. Krueger, and Y. Bengio.  NICE: Non-linear Independent Components Estimation.  In ICLR Workshop, 2015.
  * [11] L. Dinh, J. Sohl-Dickstein, and S. Bengio.  Density estimation using Real NVP.  In ICLR, 2017.
  * [12] M. Donsker and S. Varadhan.  Asymptotic evaluation of certain markov process expectations for large time, I.  Communications on Pure and Applied Mathematics, 28(1):1–47, 1975\.
  * [13] R. Dudley.  The sizes of compact subsets of Hilbert space and continuity of Gaussian processes.  Journal of Functional Analysis, 1:290–330, 1967.
  * [14] R. Dudley.  The speed of mean Glivenko-Cantelli convergence.  Ann. Math. Statist., 40:40–50, 1968.
  * [15] T. A. El Moselhy and Y. M. Marzouk.  Bayesian inference with optimal maps.  Journal of Computational Physics, 231(23):7815–7850, 2012.
  * [16] M. Germain, K. Gregor, I. Murray, and H. Larochelle.  MADE: Masked Autoencoder for Distribution Estimation.  In ICML, 2015.
  * [17] W. Grathwohl, R. Chen, J. Bettencourt, I. Sutskever, and D. Duvenaud.  FFJORD: Free-form continuous dynamics for scalable reversible generative models.  In ICLR, 2019.
  * [18] L. Györfi, M. Kohler, A. Krzyżak, and H. Walk.  A distribution-free theory of nonparametric regression, volume 1.  Springer, 2002.
  * [19] C.-W. Huang, D. Krueger, A. Lacoste, and A. Courville.  Neural autoregressive flows.  In ICML, 2018.
  * [20] J.-C. Hütter and P. Rigollet.  Minimax estimation of smooth optimal transport maps.  The Annals of Statistics, 49(2):1166–1194, 2021.
  * [21] P. Jaini, I. Kobyzev, Y. Yu, and M. Brubaker.  Tails of Lipschitz triangular flows.  In International Conference on Machine Learning, pages 4673–4681. PMLR, 2020.
  * [22] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In ICLR, 2015.
  * [23] D. Kingma and P. Dhariwal.  Glow: Generative flow with invertible 1x1 convolutions.  In Advances in Neural Information Processing Systems, pages 10215––10224, 2018.
  * [24] H. Knöthe.  Contributions to the theory of convex bodies.  The Michigan Mathematical Journal, 4(1):39–52, 1957.
  * [25] I. Kobyzev, S. Prince, and M. Brubaker.  Normalizing flows: An introduction and review of current methods.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2020\.
  * [26] Z. Kong and K. Chaudhuri.  The expressive power of a class of normalizing flow models.  In International Conference on Artificial Intelligence and Statistics, pages 3599–3609. PMLR, 2020.
  * [27] Y. Marzouk, T. Moselhy, M. Parno, and A. Spantini.  An introduction to sampling via measure transport.  arXiv: 1602.05023, 2016.
  * [28] R. Nickl and B. Pötscher.  Bracketing metric entropy rates and empirical central limit theorems for function classes of Besov- and Sobolev-Type.  J Theor Probab, 20:177–199, 2007.
  * [29] G. Papamakarios, T. Pavlakou, and I. Murray.  Masked autoregressive flow for density estimation.  In NeurIPS, 2017.
  * [30] M. Rosenblatt.  Remarks on a multivariate transformation.  The Annals of Mathematical Statistics, 23(3):470–472, 1952.
  * [31] W. Rudin.  Principles of Mathematical Analysis.  McGraw-Hill, 3rd edition, 1976.
  * [32] F. Santambrogio.  Optimal Transport for Applied Mathematicians.  Birkhäuser, 2015.
  * [33] A. Spantini, D. Bigoni, and Y. Marzouk.  Inference via low-dimensional couplings.  The Journal of Machine Learning Research, 19(1):2639–2709, 2018\.
  * [34] A. van der Vaart and J. Wellner.  Weak Convergence and Empirical Processes.  Springer, New York, 1996.
  * [35] C. Villani.  Topics in Optimal Transportation.  American Mathematical Society, 2003.
  * [36] M. Wainwright.  High-Dimensional Statistics: A Non-Asymptotic Viewpoint.  Cambridge University Press, 2019.
  * [37] A. Wehenkel and G. Louppe.  Unconstrained monotonic neural networks.  In Advances in Neural Information Processing Systems, pages 1543–1553, 2019.

This Supplement collects the detailed proofs of the theoretical results stated
in the main text and additional materials. Sec. A.1 details the derivations of
the Kullback-Leibler objective. Sec. A.2 details the proof of the slow rates
in Sec. 2 of the main text. Sec. A.3 provides estimates of metric entropy from
Sec. 3 of the main text. Sec. A.4 provides the proofs of the statistical
consistency from Sec. 3 of the main text. Sec. A.5 provides the proofs of the
Sobolev rates under log-concavity of $g$ from Sec. 4 of the main text. Sec.
A.6 expands on the dimension ordering from Sec. 4 of the main text. Sec. A.7
expands on the extension to Jacobian flows from Sec. 4 of the main text. Sec.
A.8 discusses conditions under which the KL objective (3) defining the KR
estimator $S^{n}=(S^{n}_{1},\ldots,S^{n}_{d})$ is separable in the components
$S^{n}_{k}$.

##  Appendix A Detailed proofs

###  A.1 Kullback-Leibler objective

###### Derivation of (2).

By the change of variables formula (1), the density $S\\#f$ is given by

$(S\\#f)(y)=f(S^{-1}(y))|\det(J(S^{-1})(y))|.$

Consequently, $\text{KL}(S\\#f|g)$ rewrites as

$\displaystyle\text{KL}(S\\#f|g)$
$\displaystyle=\int_{Y}(S\\#f)(y)\log\left(\frac{(S\\#f)(y)}{g(y)}\right)dy$
$\displaystyle=\int_{\mathcal{Y}}f(S^{-1}(y))|\det(J(S^{-1})(y))|\log\left(\frac{f(S^{-1}(y))|\det(J(S^{-1})(y))|}{g(y)}\right)dy$
$\displaystyle=\int_{\mathcal{Y}}f(S^{-1}(y))|\det(JS(S^{-1}(y)))|^{-1}\log\left(\frac{f(S^{-1}(y))}{g(y)|\det(JS(S^{-1}(y)))|}\right)dy$
(inverse function theorem)
$\displaystyle=\int_{\mathcal{X}}f(x)|\det(JS(x))|^{-1}\log\left(\frac{f(x)}{g(S(x))|\det(JS(x))|}\right)\cdot|\det(JS(x))|\
dx$ ($x:=S^{-1}(y)$)
$\displaystyle=\int_{\mathcal{X}}f(x)\log\left(\frac{f(x)}{g(S(x))|\det(JS(x))|}\right)dx$
$\displaystyle=\mathbb{E}_{X\sim f}\left\\{\log f(X)-\log
g(S(X))-\log|\det(JS(x))|\right\\}$ $\displaystyle=\mathbb{E}_{X\sim
f}\left\\{\log f(X)-\log g(S(X))-\sum_{k=1}^{d}\log D_{k}S_{k}(X)\right\\}.$

The last line follows because $S$ is assumed to be upper triangular and
monotone non-decreasing, and therefore

$|\det(JS(x))|=\prod_{k=1}^{d}D_{k}S_{k}(x).$

Note that in the above calculation, we have also established that

$\text{KL}(S\\#f|g)=\text{KL}(f|S^{-1}\\#g),$

for any diffeomorphism $S:\mathbb{R}^{d}\to\mathbb{R}^{d}$, since

$\displaystyle\text{KL}(S\\#f|g)$
$\displaystyle=\int_{\mathcal{X}}f(x)\log\left(\frac{f(x)}{g(S(x))|\det(JS(x))|}\right)dx$
$\displaystyle=\int_{\mathcal{X}}f(x)\log\left(\frac{f(x)}{(S^{-1}\\#g)(x)}\right)dx$
(change of variables) $\displaystyle=\text{KL}(f|S^{-1}\\#g).$

This completes the proof. ∎

###  A.2 Slow rates

####  A.2.1 Proof of Theorem 2.1

Our proof follows the argument in Section V of [3].

###### Proof of Theorem 2.1.

For fixed $\epsilon\in(0,1)$, let $\tilde{h}(x;\epsilon)$ be a $C^{\infty}$
bump function on $\mathbb{R}$ satisfying

  1. 1.

$0\leq\tilde{h}(x;\epsilon)\leq 1\qquad\forall x\in\mathbb{R}$,

  2. 2.

$\tilde{h}(x;\epsilon)=1$ on the interval $[\epsilon/4,1/2-\epsilon/4]$,

  3. 3.

$\tilde{h}(x;\epsilon)=0$ outside of the interval $[0,1/2]$.

Now for $r\in\mathbb{N}$, define the function
$h_{\epsilon,r}:[0,1]^{d}\to[-1,1]$ by

$h_{\epsilon,r}(x)=\tilde{h}(x_{1}r;\epsilon)-\tilde{h}(x_{1}r-1/2;\epsilon).$

It is clear that $h_{\epsilon,r}$ is smooth,
$\sup_{x\in[0,1]^{d}}|h_{\epsilon,r}(x)|\leq 1$, and $\int
h_{\epsilon,r}(x)dx=0$. Therefore, $1+h_{\epsilon,r}\in\mathcal{F}$ is a
smooth Lebesgue density on $[0,1]^{d}$ uniformly bounded by 2. Also note that
$|h_{\epsilon,r}(x)|=1$ whenever $[\epsilon/4\leq x_{1}r\leq 1/2-\epsilon/4]$
or $[1/2+\epsilon/4\leq x_{1}r\leq 1-\epsilon/4]$. Furthermore, the support of
$h_{\epsilon,r}$ is contained in the set $[0,1/r]\times[0,1]^{d-1}$. It
follows that

$\displaystyle\text{TV}(1+h_{\epsilon,r},1-h_{\epsilon,r})$
$\displaystyle=\frac{1}{2}\int|(1+h_{\epsilon,r}(x))-(1-h_{\epsilon,r}(x))|dx$
$\displaystyle=\int|h_{\epsilon,r}(x)|dx$
$\displaystyle\geq\int_{[\epsilon/4\leq x_{1}r\leq
1/2-\epsilon/4]\cup[1/2+\epsilon/4\leq x_{1}r\leq 1-\epsilon/4]}1\ dx_{1}$
$\displaystyle=\frac{1-\epsilon}{r},$

and

$\displaystyle\text{TV}(1+h_{\epsilon,r},1-h_{\epsilon,r})$
$\displaystyle=\int|h_{\epsilon,r}(x)|dx$ $\displaystyle\leq\int_{[0,1/r]}1\
dx_{1}$ $\displaystyle=1/r.$

As we will see, these bounds on the total variation imply that the
perturbations $1+h_{\epsilon,r},1-h_{\epsilon,r}$ are sufficiently similar to
make identification a challenging task, but sufficiently different to incur
significant loss when mistaken for each other.

Now define the translates
$h_{i}(x;\epsilon,r)=h_{\epsilon,r}\left(x_{1}-\frac{i-1}{r},x_{2},\ldots,x_{d}\right)$,
which are disjointedly supported with support contained in
$H_{i}=[(i-1)/r,i/r]\times[0,1]^{d-1}$ for $i=1,\ldots,r$. Hereafter, we
suppress dependence of $h_{i}$ on $\epsilon,r$ for notational convenience.
Consider the family of densities

$\mathcal{F}(\epsilon,r)=\left\\{1+\sum_{i=1}^{r}\delta_{i}h_{i}:\delta_{i}=\pm
1\right\\}\subset\mathcal{F}$

with cardinality $2^{r}$. For $\delta\in\\{\pm 1\\}^{r}$ we write

$f_{\delta}=1+\sum_{i=1}^{r}\delta_{i}h_{i}.$

The worst-case KL risk on $\mathcal{F}$ of any density estimate
$f_{n}=T^{n}\\#g$ derived from a KR map estimate $T^{n}$ can be bounded below
as

$\displaystyle\sup_{f\in\mathcal{F}}\mathbb{E}_{f}[\text{KL}(f|f_{n})]$
$\displaystyle\geq\sup_{f\in\mathcal{F}(\epsilon,r)}\mathbb{E}_{f}[\text{KL}(f|f_{n})]$
$\displaystyle\geq\sup_{f\in\mathcal{F}(\epsilon,r)}\mathbb{E}_{f}[2\text{TV}(f,f_{n})^{2}]$
(Pinsker’s inequality)
$\displaystyle\geq\sup_{f\in\mathcal{F}(\epsilon,r)}2\mathbb{E}_{f}[\text{TV}(f,f_{n})]^{2}.$
(Jensen’s inequality)

We aim to lower bound the total variation risk on $\mathcal{F}(\epsilon,r)$.
We have

$\displaystyle\sup_{f\in\mathcal{F}(\epsilon,r)}\mathbb{E}_{f}[\text{TV}(f,f_{n})]$
$\displaystyle\geq
2^{-r}\sum_{f\in\mathcal{F}(\epsilon,r)}\mathbb{E}_{f}[\text{TV}(f,f_{n})]$
$\displaystyle=2^{-r}\sum_{\delta\in\\{\pm
1\\}^{r}}\mathbb{E}_{f_{\delta}}\left[\text{TV}\left(f_{\delta},f_{n}\right)\right],$

i.e, the worst-case risk is larger than the Bayes risk associated to the
uniform prior on $\mathcal{F}(\epsilon,r)$. Now note that

$\displaystyle\text{TV}\left(f_{\delta},f_{n}\right)$
$\displaystyle=\frac{1}{2}\int|f_{n}(x)-f_{\delta}(x)|dx$
$\displaystyle=\frac{1}{2}\int\left|f_{n}(x)-\left(1+\sum_{i=1}^{r}\delta_{i}h_{i}(x)\right)\right|dx$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{r}\int_{H_{i}}|f_{n}(x)-(1+\delta_{i}h_{i}(x))|dx.$
($\\{H_{i}\\}$ are disjoint)

Define

$\ell_{i}(f_{n})=\frac{1}{2}\int_{H_{i}}|f_{n}(x)-(1+h_{i}(x))|dx,\qquad\ell_{i}^{\prime}(f_{n})=\frac{1}{2}\int_{H_{i}}|f_{n}(x)-(1-h_{i}(x))|dx$

and note that, by the triangle inequality,

$\ell_{i}(f_{n})+\ell_{i}^{\prime}(f_{n})\geq\frac{1}{2}\int_{H_{i}}|(1+h_{i}(x))-(1-h_{i}(x))|dx=\int|h_{i}(x)|dx\geq\frac{1-\epsilon}{r}.$

Writing $F_{\delta}^{n}$ to denote the cdf of an iid sample of size $n$ from
$f_{\delta}$, we then have

$\displaystyle
2^{-r}\sum_{\delta}\mathbb{E}_{f_{\delta}}\left[\text{TV}\left(f_{\delta},f_{n}\right)\right]$
$\displaystyle=2^{-r}\sum_{i=1}^{r}\left\\{\sum_{\delta_{i}=1}\mathbb{E}_{f_{\delta}}[\ell_{i}(f_{n})]+\sum_{\delta_{i}=-1}\mathbb{E}_{f_{\delta}}[\ell_{i}^{\prime}(f_{n})]\right\\}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{r}\left\\{\int\ell_{i}(f_{n})d\left[2^{1-r}\sum_{\delta_{i}=1}F_{\delta}^{n}\right]+\int\ell_{i}^{\prime}(f_{n})d\left[2^{1-r}\sum_{\delta_{i}=-1}F_{\delta}^{n}\right]\right\\}$
$\displaystyle\geq\frac{1}{2}\sum_{i=1}^{r}\left\\{\int[\ell_{i}(f_{n})+\ell_{i}^{\prime}(f_{n})]d\left(2^{1-r}\sum_{\delta_{i}=1}F_{\delta}^{n}\wedge
2^{1-r}\sum_{\delta_{i}=-1}F_{\delta}^{n}\right)\right\\}$
$\displaystyle\geq\frac{1-\epsilon}{2r}\sum_{i=1}^{r}\int
d\left(2^{1-r}\sum_{\delta_{i}=1}F_{\delta}^{n}\wedge
2^{1-r}\sum_{\delta_{i}=-1}F_{\delta}^{n}\right)$
$\displaystyle:=\frac{1-\epsilon}{2r}\sum_{i=1}^{r}\pi\left(2^{1-r}\sum_{\delta_{i}=1}F_{\delta}^{n},2^{1-r}\sum_{\delta_{i}=-1}F_{\delta}^{n}\right),$

where $x\wedge y=\min(x,y)$. In the last line we defined the testing affinity
$\pi$ between two distribution functions $F_{p},F_{q}$ with Lebesgue densities
$p,q$,

$\pi(F_{p},F_{q})=\int d(F_{p}\wedge F_{q})=\int(p\wedge q)dx,$

which satisfies the well known identity

$\pi(F_{p},F_{q})=1-\text{TV}(p,q).$

Since $\min$ is concave, Jensen’s inequality implies that

$\displaystyle\pi\left(2^{1-r}\sum_{\delta_{i}=1}F_{\delta}^{n},2^{1-r}\sum_{\delta_{i}=-1}F_{\delta}^{n}\right)$
$\displaystyle\geq
2^{1-r}\sum_{(\delta,\delta^{\prime})\in\Delta_{i}}\pi(F_{\delta}^{n},F_{\delta^{\prime}}^{n}),$

where
$\Delta_{i}=\\{(\delta,\delta^{\prime}):\delta_{i}=1,\delta^{\prime}_{i}=-1,\delta_{j}=\delta^{\prime}_{j}\quad\forall
j\neq i\\}$. For any $(\delta,\delta^{\prime})\in\Delta_{i}$, we have

$\displaystyle\pi(F_{\delta}^{n},F_{\delta^{\prime}}^{n})$
$\displaystyle=1-\text{TV}(f_{\delta},f_{\delta}^{\prime})$
$\displaystyle=1-\frac{1}{2}\int|f_{\delta}(x)-f_{\delta^{\prime}}(x)|dx$
$\displaystyle=1-\int|h_{i}(x)|dx$ $\displaystyle\geq 1-\frac{1}{r}.$

Hence, we conclude that

$\pi\left(2^{1-r}\sum_{\delta_{i}=1}F_{\delta}^{n},2^{1-r}\sum_{\delta_{i}=-1}F_{\delta}^{n}\right)\geq
2^{1-r}\sum_{\Delta_{i}}\left(1-\frac{1}{r}\right)=1-\frac{1}{r}.$

Thus, putting this all together, we have shown that

$\displaystyle\sup_{f\in\mathcal{F}(\epsilon,r)}\mathbb{E}_{f}[\text{TV}(f,f_{n})]$
$\displaystyle\geq\frac{1-\epsilon}{2r}\sum_{i=1}^{r}\left(1-\frac{1}{r}\right)$
$\displaystyle=\frac{1-\epsilon}{2}\left(1-\frac{1}{r}\right).$

Sending $\epsilon\to 0$ and $r\to\infty$, it follows that

$\sup_{f\in\mathcal{F}(\epsilon,r)}\mathbb{E}_{f}[\text{TV}(f,f_{n})]\geq
1/2,$

and hence

$\sup_{f\in\mathcal{F}}\mathbb{E}_{f}[\text{KL}(f|f_{n})]\geq\sup_{f\in\mathcal{F}(\epsilon,r)}2\mathbb{E}_{f}[\text{TV}(f,f_{n})]^{2}\geq
1/2.$

This completes the proof. ∎

####  A.2.2 Proof of Theorem 2.2

###### Proof.

Let $h_{n}$ denote the first marginal of the density estimate
$f_{n}=T^{n}\\#g$. Note that $h_{n}$ is a density on $\mathbb{R}$. Defining
$\pi_{1}:\mathbb{R}^{d}\to\mathbb{R}$ to be the projection along the first
factor, we have

$h_{n}=\pi_{1}\\#f_{n}=(\pi_{1}\circ T^{n})\\#g.$

By Problem 7.5 in Devroye et al. (1996) [9], for any positive sequence
$1/16\geq b_{1}\geq b_{2}\cdots$ converging to zero and any density estimate
$h_{n}$ there exists a density $h$ on $\mathbb{R}$ such that

$\mathbb{E}\left\\{\int|h(x)-h_{n}(x)|dx\right\\}\geq b_{n}.$

Letting TV denote the total variation distance, this inequality can be
rewritten as

$2\mathbb{E}[\text{TV}(h,h_{n})]\geq b_{n}.$

Setting $a_{n}=b_{n}^{2}/2$, which satisfies $a_{1}\leq 1/512$, we find that

$\displaystyle\mathbb{E}[\text{KL}(h|h_{n}))]$
$\displaystyle\geq\mathbb{E}[2\cdot\text{TV}(h,h_{n})^{2}]$ (Pinsker’s
inequality) $\displaystyle\geq 2\mathbb{E}[\text{TV}(h,h_{n})]^{2}$ (Jensen’s
inequality) $\displaystyle\geq b_{n}^{2}/2$ $\displaystyle=a_{n}.$

Finally, let $f=h^{\otimes d}$ be the density on $\mathbb{R}^{d}$ defined as a
$d$-fold product of $h$. By the chain rule of relative entropy, it follows
that

$\displaystyle\mathbb{E}[\text{KL}(f|f_{n})]\geq\mathbb{E}[\text{KL}(h|h_{n})]\geq
a_{n}.$

This completes the proof. ∎

###  A.3 Upper bounds on metric entropy

We begin by defining relevant Sobolev function spaces, for which metric
entropy bounds are known.

######  Definition A.1.

For $\mathcal{X}\subseteq\mathbb{R}^{d}$, define the function space

$D_{s}(\mathcal{X})=\\{\phi:D^{\alpha}\phi\text{ are uniformly continuous for
all }|\alpha|\leq s\\}.$

and its subset

$C_{s}(\mathcal{X})=\left\\{\phi:\mathcal{X}\to\mathbb{R}:\sum_{0\leq|\alpha|\leq
s}\|D^{\alpha}\phi\|_{\infty}<\infty\right\\}\cap D_{s}(\mathcal{X}).$

endowed with the Sobolev norm

$\|\phi\|_{H^{s,\infty}(\mathcal{X})}=\sum_{|\alpha|\leq
s}\|D^{\alpha}\phi\|_{\infty}.$

######  Proposition A.2 (Corollary 3, Nickl and Pötscher (2007) [28]).

Assume $\mathcal{X}\subset\mathbb{R}^{d}$ is compact and let $\mathcal{F}$ be
a bounded subset of $C_{s}(\mathcal{X})$ with respect to
$\|\cdot\|_{H^{s,\infty}(\mathcal{X})}$ for some $s>0$. The metric entropy of
$\mathcal{F}$ in the $L^{\infty}$ norm is bounded as

$H(\epsilon,\mathcal{F},\|\cdot\|_{\infty})\lesssim\epsilon^{-d/s}.$

With this result in hand, we can proceed to the proof of Proposition 3.3.

###### Proof of Proposition 3.3.

This is a direct consequence of Proposition A.2. Indeed, under Assumptions
3.1-3.3 and Definition 3.1, for every $S\in\mathcal{T}(s,d,M)$, every term in
$\exp(\psi_{S}(x))$ is bounded away from 0 and $s$-smooth with uniformly
bounded derivatives. Since $\log$ is smooth away from 0, it follows that that
every $\psi_{S}\in\Psi(s,d,M)$ is $s$-smooth with uniformly bounded
derivatives. Consequently, $\Psi(s,d,M)$ is a bounded subset of
$C_{s}(\mathcal{X})$ for every $M>0$. ∎

We also provide a metric entropy bound for the space
$\mathcal{T}(\mathfrak{s},d,M)$, which in particular establishes compactness
of $\overline{\mathcal{T}(\mathfrak{s},d,M)}$ with respect to
$\|\cdot\|_{\infty,d}$ (although this can also be proved with the Arzelà-
Ascoli theorem).

######  Proposition A.3.

Let $\mathfrak{s}=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}_{+}$, $d_{k}=d-k+1$,
and
$\tilde{\sigma}_{k}=d_{k}\left((s_{k}+1)^{-1}+\sum_{j=k+1}^{d}s_{j}^{-1}\right)^{-1}$
for $k\in[d]$. Under Assumptions 3.1, 3.2, and 4.1 the space
$\mathcal{T}(\mathfrak{s},d,M)$ is totally bounded (and therefore precompact)
with respect to the uniform norm $\|\cdot\|_{\infty,d}$ with metric entropy
satisfying

$\displaystyle
H(\epsilon,\mathcal{T}(\mathfrak{s},d,M),\|\cdot\|_{\infty,d})\leq$
$\displaystyle\sum_{k=1}^{d}c_{k}(\epsilon/2M)^{-d_{k}/\tilde{\sigma}_{k}}$

for some positive constants $c_{k},k\in[d]$ independent of $\epsilon$ and $M$.

This result relies on known metric entropy bounds for anisotropic smoothness
classes.

######  Proposition A.4 (Prop. 2.2, Birgé (1986) [3]).

Let $\mathfrak{s}=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}_{+}$ and
$\sigma=d\left(\sum_{j=1}^{d}s_{j}^{-1}\right)^{-1}$. Assume that $\Phi$ is a
family of functions $\mathbb{R}^{d}\to\mathbb{R}$ with common compact convex
support of dimension $d$ and satisfying

$\sup_{\phi\in\Phi,\alpha\preceq\mathfrak{s}}\|D^{\alpha}\phi\|_{\infty}<\infty.$

The metric entropy of $\Phi$ in the $L^{\infty}$ norm is bounded as

$H(\epsilon,\Phi,\|\cdot\|_{\infty})\lesssim\epsilon^{-d/\sigma}.$

We now proceed to the proof of Proposition A.3.

###### Proof of Proposition A.3.

For every $k\in[d]$, define the set of functions
$\mathcal{X}_{k:d}\to\mathbb{R}$ given by

$\mathcal{T}_{k}=\\{S_{k}:S\in\mathcal{T}(\mathfrak{s},d,M)\\}.$

By Definition 3.1, for each $k\in[d]$ we have that $\mathcal{T}_{k}/M$
satisfies the assumptions of Proposition A.4, and hence

$n_{k}:=H(\epsilon,\mathcal{T}_{k},L^{\infty})\leq
c_{k}(\epsilon/M)^{-d_{k}/\tilde{\sigma}_{k}},$

for some $c_{k}>0$ independent of $\epsilon$ and $M$.

Now note that
$\mathcal{T}(\mathfrak{s},d,M)\subseteq\prod_{k=1}^{d}\mathcal{T}_{k}$. For
each $k\in[d]$, let $\\{g_{k,1},\ldots,g_{k,n_{k}}\\}$ be a minimal
$\epsilon$-cover of $\mathcal{T}_{k}$ with respect to the $L^{\infty}$ norm,
and define the subset

$\mathcal{E}=\\{f_{i_{1},\ldots,i_{d}}=(g_{1,i_{i}},\ldots,g_{d,i_{d}}):1\leq
i_{k}\leq n_{k}\\}\subset\prod_{k=1}^{d}\mathcal{T}_{k}$

which has cardinality $\prod_{k=1}^{d}n_{k}$. Now fix an arbitrary
$f\in\prod_{k=1}^{d}\mathcal{T}_{k}$. For each $k\in[d]$, we can find some
$g_{k,i_{k}}$ such that $\|f_{k}-g_{k,i_{k}}\|_{\infty}\leq\epsilon$. It
follows that

$\|f-f_{i_{1},\ldots,i_{d}}\|_{\infty,d}=\max_{k\in[d]}\|f_{k}-g_{k,i_{k}}\|_{\infty}\leq\epsilon.$

Hence, $\mathcal{E}$ is an $\epsilon$-cover of
$\prod_{k=1}^{d}\mathcal{T}_{k}$ and so

$H\left(\epsilon,\prod_{k=1}^{d}\mathcal{T}_{k},\|\cdot\|_{\infty,d}\right)\leq\log\left(\prod_{k=1}^{d}n_{k}\right)=\sum_{k=1}^{d}n_{k}.$

To conclude the proof, we claim that
$\mathcal{T}(\mathfrak{s},d,M)\subseteq\prod_{k=1}^{d}\mathcal{T}_{k}$ implies
that

$H\left(\epsilon,\mathcal{T}(\mathfrak{s},d,M),\|\cdot\|_{\infty,d}\right)\leq
H\left(\epsilon/2,\prod_{k=1}^{d}\mathcal{T}_{k},\|\cdot\|_{\infty,d}\right).$

Indeed, suppose
$\\{f_{1},\ldots,f_{m}\\}\subseteq\prod_{k=1}^{d}\mathcal{T}_{k}$ is a finite
$(\epsilon/2)$-cover of $\prod_{k=1}^{d}\mathcal{T}_{k}$. For each
$j=1,\ldots,m$, if $\mathcal{T}(\mathfrak{s},d,M)\cap
B_{\infty}(f_{j},\epsilon/2)$ is non-empty, we define $g_{j}$ to be an element
of this intersection. Here $B_{\infty}(f_{j},\epsilon/2)$ denotes the ball of
radius $\epsilon/2$ in $\prod_{k=1}^{d}\mathcal{T}_{k}$ centered at $f_{j}$
with respect to the norm $\|\cdot\|_{\infty,d}$. Now let
$g\in\mathcal{T}(\mathfrak{s},d,M)$ be arbitrary. Since
$g\in\prod_{k=1}^{d}\mathcal{T}_{k}$, there is some $j\in\\{1,\ldots,m\\}$
such that $\|g-f_{j}\|_{\infty,d}\leq\epsilon/2$. This implies that $g_{j}$ is
defined and hence

$\|g-g_{j}\|_{\infty,d}\leq\|g-f_{j}\|_{\infty,d}+\|f_{j}-g_{j}\|_{\infty,d}\leq\epsilon/2+\epsilon/2=\epsilon.$

It follows that $\\{g_{j}\\}$ is a finite $\epsilon$-cover of
$\mathcal{T}(\mathfrak{s},d,M)$ with respect to $\|\cdot\|_{\infty,d}$, which
establishes the claim and completes the proof.

∎

###  A.4 Statistical consistency

####  A.4.1 Proof of Lemma 3.4

###### Proof of Lemma 3.4.

Assume $d<2s$. By Theorem 2.14.2 in van der Vaart and Wellner (1996) [34],
since functions in $\Psi(s,d,M)$ are by definition uniformly bounded there
exist positive constants $C_{1},C_{2}$ such that

$\displaystyle\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M)}$
$\displaystyle\leq\frac{C_{1}C_{2}}{\sqrt{n}}\int_{0}^{1}\sqrt{1+H(C_{2}\epsilon,\Psi(s,d,M),\|\cdot\|_{\infty})}d\epsilon$
$\displaystyle\lesssim\frac{1}{\sqrt{n}}\int_{0}^{1}\sqrt{1+\epsilon^{-d/s}}d\epsilon$
(Proposition 3.3) $\displaystyle\lesssim n^{-1/2}.$ ($d<2s$)

The last line follows since $d<2s$ implies that the integral on the right side
is finite.

When $d\geq 2s$, the metric entropy integral above is no longer finite. In
this case, we appeal to Dudley’s metric entropy integral bound [13] (see also
Theorem 5.22 in [36]), which states that there exists positive constants
$C_{3},D>0$ for which

$\displaystyle\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M)}$
$\displaystyle\leq\min_{\delta\in[0,D]}\left\\{\delta+\frac{C_{3}}{\sqrt{n}}\int_{\delta}^{D}\sqrt{H(\epsilon,\Psi(s,d,M),\|\cdot\|_{\infty})}d\epsilon\right\\}$
$\displaystyle\leq\min_{\delta\in[0,D]}\left\\{\delta+\frac{C_{3}\sqrt{c}}{\sqrt{n}}\int_{\delta}^{D}\epsilon^{-d/2s}d\epsilon\right\\}.$
(Proposition 3.3)

First assume $d=2s$. Evaluating the integral in Dudley’s bound, we obtain

$\displaystyle\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M)}$
$\displaystyle\leq\min_{\delta\in[0,D]}\left\\{\delta+\frac{C_{3}\sqrt{c}}{\sqrt{n}}[\log
D-\log\delta]\right\\}.$

To minimize the expression on the right side in $\delta$, we differentiate
with respect to $\delta$ and find where the derivative vanishes. The bound is
optimized by choosing $\delta$ proportional to $n^{-1/2}$, which implies that
$\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M)}\lesssim n^{-1/2}\log n$.

Now assume $d>2s$. Evaluating the integral in Dudley’s bound, we obtain

$\displaystyle\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M)}$
$\displaystyle\leq\min_{\delta\in[0,D]}\left\\{\delta+\frac{C_{3}\sqrt{c}}{\sqrt{n}}\int_{\delta}^{\infty}\epsilon^{-d/2s}d\epsilon\right\\}$
$\displaystyle=\min_{\delta\in[0,D]}\left\\{\delta+\frac{C_{4}}{\sqrt{n}}\delta^{1-d/2s}\right\\}.$

We optimize this bound by choosing $\delta$ proportional to $n^{-s/d}$, which
implies that $\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M)}\lesssim n^{-s/d}$. ∎

####  A.4.2 Proof of Theorem 3.5

Before continuing on to the proof of Theorem 3.5, we first prove that indeed
$S^{*}\in\mathcal{T}(s,d,M^{*})$ for $M^{*}>0$ sufficiently large, and
similarly for $T^{*}$.

######  Lemma A.5.

Let $\mathfrak{s}=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}_{+}$. Under
Assumptions 3.1, 3.2, and 4.1, there exists some $M^{*}>0$ such that the KR
map $S^{*}$ from $f$ to $g$ lies in $\mathcal{T}(\mathfrak{s},d,M)$ for all
$M\geq M^{*}$.

By symmetry, if we switch the roles of $f$ and $g$ in Assumption 4.1, we see
that the same is true for the KR map $T^{*}$ from $g$ to $f$. In particular,
if $\mathfrak{s}=(s,\ldots,s)$ for some $s\in\mathbb{N}$, then we are in the
case of Assumption 3.3. It then follows from Lemma A.5 that
$S^{*},T^{*}\in\mathcal{T}(s,d,M^{*})$ for $M^{*}>0$ sufficiently large.

The proof of Lemma A.5 requires an auxiliary result establishing the
regularity of marginal and conditional pdfs and quantile functions associated
to a smooth density.

######  Lemma A.6.

Let $h$ be a density on $\mathbb{R}^{d}$ with compact support $\mathcal{Z}$ .
Assume further that $h$ is strictly positive on $\mathcal{Z}$ and $s$-smooth
in the sense of Assumption 4.1, where
$s=(s_{1},\ldots,s_{d})\in\mathbb{Z}^{d}_{+}$. The following hold.

  1. 1.

For any index set $J\subseteq[d]$, the marginal density
$h(x_{J})=\int_{\mathbb{R}^{d-|J|}}h(x_{J},x_{-J})dx_{-J}$ is $s_{J}$-smooth.
Similarly, the conditional density
$h(x_{J}|x_{-J})=\frac{h(x_{J},x_{-J})}{h(x_{-J})}$ is $s$-smooth as a
function of $x=(x_{J},x_{-J})$.

  2. 2.

For any $J\subseteq[d],k\in J^{c}$, the univariate conditional cdf
$H(x_{k}|x_{J})=\int_{-\infty}^{x_{k}}h(y_{k}|x_{J})dy_{k}$ is
$(s_{k}+1,s_{J})$-smooth as a function of $(x_{k},x_{J})$.

  3. 3.

Assume further that $\mathcal{Z}$ is convex. Then the conditional quantile
function $H^{-1}(p_{k}|x_{J})$ is $(s_{k}+1,s_{J})$-smooth as a function of
$(p_{k},x_{J})$ for any $J\subseteq[d],k\in J^{c},p_{k}\in[0,1]$.

###### Proof.

For (1), the regularity assumptions on $h$ imply that we can differentiate
under the integral. For any multi-index $\alpha\preceq s$ satisfying
$\alpha_{-J}=0$ we have

$\displaystyle D^{\alpha}h(x_{J})$ $\displaystyle=\int
D^{\alpha}h(x_{J},x_{-J})dx_{-J}.$

Since $D^{\alpha}h(x_{J},x_{-J})$ is continuous and compactly supported by
hypothesis, the dominated convergence theorem implies that for any sequence
$x_{J,n}\to x_{J}$ we have

$D^{\alpha}h(x_{J,n})=\int D^{\alpha}h(x_{J,n},x_{-J})dx_{-J}\to\int
D^{\alpha}h(x_{J},x_{-J})dx_{-J}=D^{\alpha}h(x_{J}).$

This proves that $h(x_{J})$ is $s_{J}$-smooth. Since $h>0$ on its support and
$x\mapsto 1/x$ is a smooth function for $x>0$, this implies that $1/h(x_{-J})$
is $s_{-J}$-smooth. As products of differentiable functions are
differentiable, we conclude that $h(x_{J}|x_{-J})=h(x_{J},x_{-J})/h(x_{-J})$
is $s$-smooth.

For (2), the result from (1) implies that the univariate conditional density
$h(x_{k}|x_{J})$ is $(s_{k},s_{J})$-smooth. For any multi-index $\alpha\preceq
s$ satisfying $\alpha_{k}=0$, we can repeat the differentiation-under-the-
integral argument above, combined with the dominated convergence theorem, to
see that $D^{\alpha}H(x_{k}|x_{J})$ is continuous. If $\alpha_{k}\neq 0$, we
apply the fundamental theorem of calculus to take care of one derivative in
$x_{k}$ and then recall that $h(x_{k}|x_{J})$ is $(s_{k},s_{J})$-smooth to
complete the proof, appealing to Clairaut’s theorem to exchange the order of
differentiation. Note that this shows that $H(x_{k}|x_{J})$ is
$(s_{k}+1,s_{J})$-smooth as a function of $(x_{k},x_{J})$.

For (3), the assumption that $\mathcal{Z}$ is convex, combined with the fact
that $h>0$ on its support, implies that $H(x_{k}|x_{J})$ is strictly
increasing in $x_{k}$. As such, $H^{-1}(p_{k}|x_{J})$ is defined, strictly
increasing, and continuous in $p_{k}$. Since $H(x_{k}|x_{J})$ is
differentiable in $x_{k}$, we have

$\frac{\partial}{\partial
p_{k}}H^{-1}(p_{k}|x_{J})=\frac{1}{\frac{\partial}{\partial
x_{k}}H(H^{-1}(p_{k}|x_{J})|x_{J})}=\frac{1}{h(H^{-1}(p_{k}|x_{J})|x_{J})},$

which is continuous in $p_{k}$ as a composition of continuous functions.
Continuing in this way, we see that $H^{-1}(p_{k}|x_{J})$ is
$(s_{k}+1)$-smooth in $p_{k}$.

For the other variables, note that the following relation holds generally:

$\displaystyle H(H^{-1}(p_{k}|x_{J})|x_{J})$ $\displaystyle=p_{k}.$

Defining the function

$u(p_{k},x_{J},x_{k})=H(x_{k}|x_{J})-p_{k},$

we see that $u$ is $(\infty,s_{J},s_{k}+1)$-smooth in its arguments
$(p_{k},x_{J},x_{k})$ and satisfies

$u(p_{k},x_{J},H^{-1}(p_{k}|x_{J}))=0.$

Furthermore $\frac{\partial}{\partial x_{k}}u>0$ by assumption, since $h>0$ on
its support. Hence we can appeal to the implicit function theorem (stated
precisely below) to conclude that $H^{-1}(p_{k}|x_{J})$ is
$(s_{k}+1,s_{J})$-smooth in $(p_{k},x_{J})$. For example, for any $j\in J$ we
can evaluate the partial derivative in $x_{j}$ by implicitly differentiating
the above relation. Letting $x_{k}(p_{k},x_{J})=H^{-1}(p_{k}|x_{J})$ we obtain

$\displaystyle\frac{\partial}{\partial x_{j}}H(x_{k}(p_{k},x_{J})|x_{J})$
$\displaystyle=\frac{\partial H}{\partial
x_{k}}(x_{k}(p_{k},x_{J})|x_{J})\cdot\frac{\partial x_{k}}{\partial
x_{j}}(p_{k},x_{J})+\frac{\partial H}{\partial x_{j}}(x_{k}|x_{J})$
$\displaystyle=\frac{\partial p_{k}}{\partial x_{j}}$ $\displaystyle=0.$

Rearranging terms and noting that $\frac{\partial x_{k}}{\partial
x_{j}}(p_{k},x_{J})=\frac{\partial}{\partial x_{j}}H^{-1}(p_{k}|x_{J})$, we
find that

$\frac{\partial}{\partial x_{j}}H^{-1}(p_{k}|x_{J})=-\frac{\frac{\partial
H}{\partial x_{j}}(H^{-1}(p_{k}|x_{J})|x_{J})}{h(H^{-1}(p_{k}|x_{J})|x_{J})}.$

∎

######  Theorem A.7 (Implicit function theorem [31]).

Let $u:\mathbb{R}^{n+m}\to\mathbb{R}^{m}$ be continuously differentiable and
suppose $u(a,b)=0$ for some $(a,b)\in\mathbb{R}^{n}\times\mathbb{R}^{m}$. For
$(x,y)\in\mathbb{R}^{n}\times\mathbb{R}^{m}$ define the partial Jacobians

$J_{n}u(x,y)=\left[\frac{\partial}{\partial
x_{j}}u_{i}(a,b)\right]_{(i,j)\in[m]\times[n]}\qquad
J_{m}u(x,y)=\left[\frac{\partial}{\partial
y_{j}}u_{i}(a,b)\right]_{(i,j)\in[m]\times[m]},$

and suppose that $J_{m}u(a,b)$ is invertible. Then there exists an open
neighborhood $V\subseteq\mathbb{R}^{n}$ of $a$ such that there exists a unique
continuously differentiable function $v:V\to\mathbb{R}^{m}$ satisfying
$v(a)=b$, $u(x,v(x))=0$ for all $x\in V$, and

$Jv(x)=-J_{m}u(x,v(x))^{-1}J_{n}u(x,v(x)).$

We now proceed to the proof of Lemma A.5.

###### Proof of Lemma A.5.

By definition of the KR map, for each $k\in[d]$ we have

$S^{*}_{k}(x_{k}|x_{(k+1):d})=G_{k}^{-1}(F_{k}(x_{k}|x_{(k+1):d})|S^{*}_{(k+1):d}(x)).$

The requisite smoothness in Definition 3.1 of $S^{*}_{k:d}$ then follows from
the chain rule of differentiation, appealing to Assumption 4.1, and an
application of Lemma A.6. Note also that $x_{k}\mapsto
G_{k}^{-1}(F_{k}(x_{k}|x_{(k+1):d})|S^{*}_{(k+1):d}(x))$ is strictly
increasing as a composition of strictly increasing functions, since $f,g$ are
bounded away from 0 on their supports by assumption. Defining
$\mathcal{T}(\mathfrak{s},d)\subset\mathcal{T}$ as the subset of strictly
increasing triangular maps that are $\mathfrak{s}$-smooth (although not
necessarily with uniformly bounded derivatives), we see that
$S^{*}\in\mathcal{T}(\mathfrak{s},d)$. Since
$\mathcal{T}(\mathfrak{s},d)=\cup_{M>0}\mathcal{T}(\mathfrak{s},d,M)$, there
exists some $M^{*}>0$ for which $S^{*}\in\mathcal{T}(\mathfrak{s},d,M^{*})$.
Since
$\mathcal{T}(\mathfrak{s},d,M_{1})\subseteq\mathcal{T}(\mathfrak{s},d,M_{2})$
for all $M_{1}\leq M_{2}$ by definition, we conclude that
$S^{*}\in\mathcal{T}(\mathfrak{s},d,M)$ for all $M\geq M^{*}$. ∎

With Lemma A.5 in hand, we can now prove Theorem 3.5.

###### Proof of Theorem 3.5.

The theorem is a direct result of inequality (4) and Lemma 3.4. Indeed, since
$S^{*}\in\mathcal{T}(s,d,M^{*})$ by Lemma A.5, we have
$\inf_{S\in\mathcal{T}(s,d,M^{*})}P\psi_{S}=P\psi_{S^{*}}=0$, and therefore

$\displaystyle\mathbb{E}[P\psi_{S^{n}}]$
$\displaystyle=\mathbb{E}\left\\{P\psi_{S^{n}}-\inf_{S\in\mathcal{T}(s,d,M^{*})}P\psi_{S}\right\\}$
$\displaystyle\leq\mathbb{E}\left\\{2\|P_{n}-P\|_{\Psi(s,d,M^{*})}+R_{n}\right\\}$
(inequality (4))
$\displaystyle\lesssim\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M^{*})}.$

Appealing to Lemma 3.4 establishes the stated bound on
$\mathbb{E}[P\psi_{S^{n}}]$. To conclude that
$P\psi_{S^{n}}\stackrel{{\scriptstyle p}}{{\to}}0$, we apply Markov’s
inequality:

$P(P\psi_{S^{n}}\geq\epsilon)\leq\epsilon^{-1}\mathbb{E}[P\psi_{S^{n}}]\to 0.$

As $\epsilon>0$ was arbitrary, this completes the proof. ∎

####  A.4.3 Proof of KL lower semicontinuity

Although Theorem 3.5 only establishes a weak form of consistency in terms of
the KL divergence, we leverage this result to prove strong consistency, in the
sense of uniform convergence of $S^{n}$ to $S^{*}$ in probability, in Theorem
3.6. The proof requires understanding the regularity of the KL divergence with
respect to the topology induced by the $\|\cdot\|_{\infty,d}$ norm. Lemma A.8
establishes that KL is lower semicontinuous with respect to this topology. It
relies on the weak lower semicontinuity of KL proved by Donsker and Varadhan
using their dual representation in Lemma 2.1 of [12].

######  Lemma A.8.

Under Assumptions 3.1-3.3, the functional $S\mapsto P\psi_{S}$ on the domain
$\mathcal{T}$ is lower semicontinuous with respect to the uniform norm
$\|\cdot\|_{\infty,d}$.

###### Proof.

Assume $\|S^{n}-S\|_{\infty,d}\to 0$ for fixed maps $S^{n},S\in\mathcal{T}$.
We first claim that $\nu_{n}:=S^{n}\\#\mu\rightsquigarrow S\\#\mu=:\nu$ for
any probability measure $\mu$. We will show that $\int h\ d\nu_{n}\to\int h\
d\nu$ for all bounded continuous functions $h\in C_{b}(\mathbb{R}^{d})$. Note
that $\int h\ d\nu_{n}=\int h\circ S^{n}\ d\mu$. By assumption, $h\circ S^{n}$
is bounded and measurable with $\|h\circ S^{n}\|_{\infty}\leq\|h\|_{\infty}$
for all $n$. Furthermore, $h\circ S^{n}\to h\circ S$ pointwise, since $h$ is
continuous. Then by the dominated convergence theorem, $\int|h\circ
S^{n}-h\circ S|\ d\mu\to 0$. In particular, this implies that $\int h\
d\nu_{n}=\int h\circ S^{n}\ d\mu\to\int h\circ S\ d\mu=\int h\ d\nu$. This
proves the claim.

Combining this fact with the weak lower semicontinuity of KL divergence [12]
proves the lemma. Indeed, since $\|S^{n}-S\|_{\infty,d}\to 0$ implies that
$S^{n}\\#\mu\rightsquigarrow S\\#\mu$ for any $\mu$, it follows that

$\displaystyle P\psi_{S}$ $\displaystyle=\text{KL}(S\\#f|g)$
$\displaystyle\leq\liminf_{n\to\infty}\text{KL}(S^{n}\\#f|g)$ (KL is weakly
lower semicontinuous) $\displaystyle=\liminf_{n\to\infty}P\psi_{S^{n}}.$

Thus, $S\mapsto P\psi_{S}$ on $\mathcal{T}$ is lower semicontinuous with
respect to $\|\cdot\|_{\infty,d}$. ∎

As a corollary, we also obtain an existence result for minimizers of $S\mapsto
P\psi_{S}$ on $\overline{\mathcal{T}(s,d,M)}$.

######  Corollary A.9.

Under Assumptions 3.1-3.3, for any $M>0$ the minimum
$\inf_{S\in\overline{\mathcal{T}(s,d,M)}}P\psi_{S}$ is attained.

###### Proof.

This follows from the direct method of calculus of variations, since
$\overline{\mathcal{T}(s,d,M)}$ is compact with respect to
$\|\cdot\|_{\infty,d}$ (Proposition A.3) and $S\mapsto P\psi_{S}$ is bounded
below and lower semicontinuous with respect to $\|\cdot\|_{\infty,d}$. ∎

####  A.4.4 Uniform consistency of the inverse map

We have proved consistency of the estimator $S^{n}$ of the sampling map
$S^{*}$, which pushes forward the target density $f$ to the source density
$g$. However, we require knowledge of the direct map $T^{*}=(S^{*})^{-1}$ to
generate new samples from $f$ by pushing forward samples from $g$ under
$T^{*}$. In this section we prove consistency and a rate of convergence of the
estimator $T^{n}=(S^{n})^{-1}$ of the direct map $T^{*}$.

First note that $\text{KL}(S\\#f|g)=\text{KL}(f|S^{-1}\\#g)$ for any
diffeomorphism $S:\mathbb{R}^{d}\to\mathbb{R}^{d}$, as we proved in Section
A.1. As such, the consistency and rate of convergence of
$\text{KL}(S^{n}\\#f|g)$ obtained in Theorem 3.5 yield the same results for
the estimator $T^{n}$ in terms of $\text{KL}(f|T^{n}\\#g)$ under identical
assumptions. We can also establish uniform consistency of $T^{n}$ as we did
for $S^{n}$ in Theorem 3.6, although the proof of this fact requires a bit
more work.

######  Theorem A.10.

Suppose Assumptions 3.1-3.3 hold. Let $S^{n}$ be any near-optimizer of the
functional $S\mapsto P_{n}\psi_{S}$ on $\mathcal{T}(s,d,M^{*})$, i.e., suppose

$P_{n}\psi_{S^{n}}=\inf_{S\in\mathcal{T}(s,d,M^{*})}P_{n}\psi_{S}+o_{P}(1).$

Let $T^{n}=(S^{n})^{-1}$. Then
$\|T^{n}-T^{*}\|_{\infty,d}\stackrel{{\scriptstyle p}}{{\to}}0$, i.e, $T^{n}$
is a uniformly consistent estimator of $T^{*}$.

The proof of Theorem A.10 relies heavily on the bounds on $S^{n}$ and its
derivatives posited in Definition 3.1, which we utilize in conjunction with
the inverse function theorem to uniformly bound the Jacobian
$JT^{n}(y)=(JS^{n}(T^{n}(y)))^{-1}$ over $y\in\mathcal{Y}$ and
$n\in\mathbb{N}$, thereby establishing uniform equicontinuity of the family of
estimators $\\{T^{n}\\}_{n=1}^{\infty}$. We combine this uniform
equicontinuity with the uniform consistency of $S^{n}$ from Theorem 3.6 to
complete the proof.

First we establish a lemma that allows us to bound the derivatives of the
inverse map estimates $T^{n}$.

######  Lemma A.11.

Suppose $A\in\mathbb{R}^{d\times d}$ is an invertible upper triangular matrix
satisfying

$\max_{i,j\in[d],i<j}|A_{ij}|\leq
L\qquad\text{and}\qquad\min_{j\in[d]}|A_{jj}|\geq 1/M$

for some positive $L,M>0$. Then $A^{-1}$ is upper triangular and the diagonal
entries are bounded as

$\max_{j\in[d]}|A^{-1}_{jj}|\leq M.$

Furthermore, the superdiagonal terms $i,j\in[d]$ with $i<j$ are bounded as

$|A_{ij}^{-1}|\leq M^{2}L(ML+1)^{j-i-1}.$

###### Proof.

Let $D=\text{diag}(A)$ denote the matrix with diagonal entries $D_{jj}=A_{jj}$
for $j\in[d]$ and zeros elsewhere. Note that $D$ is invertible since $A$ is,
and $D^{-1}$ is diagonal with entries $D^{-1}_{jj}=1/A_{jj}$ for $j\in[d]$.
Let $U=A-D$ denote the strictly upper triangular part of $A$. Now note that

$A^{-1}=(D+U)^{-1}=[D(I+D^{-1}U)]^{-1}=(I+D^{-1}U)^{-1}D^{-1}.$

To calculate $(I+D^{-1}U)^{-1}$ we make use of the matrix identity

$(I-X)\left[\sum_{k=0}^{d-1}X^{k}\right]=I-X^{d}.$

We plug in $X=-D^{-1}U$ and note that $D^{-1}U$ is strictly upper triangular,
which implies that it is nilpotent of degree at most $d$, i.e.,
$(D^{-1}U)^{d}=0$. Hence, we obtain

$(I+D^{-1}U)\left[\sum_{k=0}^{d-1}(-D^{-1}U)^{k}\right]=I-(-D^{-1}U)^{d}=I,$

which implies that $(I+D^{-1}U)^{-1}=\sum_{k=0}^{d-1}(-D^{-1}U)^{k}$. It
follows that

$A^{-1}=\sum_{k=0}^{d-1}(-D^{-1}U)^{k}D^{-1}.$

This shows that $A^{-1}$ is upper triangular as a sum of products of upper
triangular matrices. Note also that the terms in the sum with $k>0$ are
strictly upper triangular. Hence, we see that $\text{diag}(A^{-1})=D^{-1}$ and
therefore the diagonal entries of $A^{-1}$ satisfy
$|A^{-1}_{jj}|=\frac{1}{|A_{jj}|}\leq M$ by hypothesis.

Now we bound the superdiagonal entries. By repeated application of the
triangle inequality and appealing to the bounds on $A_{ij}$, we can bound
$A^{-1}_{ij}$ for $i<j$ as

$\displaystyle|A^{-1}_{ij}|$
$\displaystyle=\left|\sum_{k=0}^{d-1}[(D^{-1}U)^{k}D^{-1}]_{ij}\right|$
$\displaystyle=\left|\sum_{k=1}^{d-1}[(D^{-1}U)^{k}D^{-1}]_{ij}\right|$
($i<j$) $\displaystyle\leq\sum_{k=1}^{d-1}|[(D^{-1}U)^{k}D^{-1}]_{ij}|$
$\displaystyle\leq\sum_{k=1}^{d-1}|[(MU)^{k}M]_{ij}|$
($|D^{-1}_{jj}|=|A^{-1}_{jj}|\leq M$)
$\displaystyle=\sum_{k=1}^{d-1}M^{k+1}|(U^{k})_{ij}|.$

Now let $V$ denote the $d\times d$ matrix with ones strictly above the
diagonal and zeros elsewhere, i.e.,

$V_{ij}=\begin{cases}1,&i<j,\\\ 0,&i\geq j.\end{cases}$

We now claim that for $k\geq 1$,

$(V^{k})_{ij}=\begin{cases}\binom{j-i-1}{k-1},&k\leq j-i\\\
0,&\text{else}.\end{cases}$

Our proof of the claim proceeds by induction. The statement is clearly true
for the base case $k=1$ by definition of $V$. Suppose the claim holds up to
$k$. It follows that

$\displaystyle(V^{k+1})_{ij}$ $\displaystyle=(V^{k}\cdot V)_{ij}$
$\displaystyle=\sum_{\ell=1}^{d}(V^{k})_{i\ell}V_{\ell j}$
$\displaystyle=\sum_{\ell=1}^{d}\binom{\ell-i-1}{k-1}1_{[k\leq\ell-i]}\cdot
1_{[\ell<j]}$ (inductive step) $\displaystyle=1_{[k+1\leq
j-i]}\sum_{\ell=i+k}^{j-1}\binom{\ell-i-1}{k-1}$ $\displaystyle=1_{[k+1\leq
j-i]}\sum_{p=k-1}^{j-i-2}\binom{p}{k-1}$ ($p:=\ell-i-1$)
$\displaystyle=1_{[k+1\leq
j-i]}\left\\{1+\sum_{p=k}^{j-i-2}\binom{p}{k-1}\right\\}$
$\displaystyle=1_{[k+1\leq
j-i]}\left\\{1+\sum_{p=k}^{j-i-2}\left[\binom{p+1}{k}-\binom{p}{k}\right]\right\\}$
(Pascal’s formula) $\displaystyle=1_{[k+1\leq
j-i]}\left\\{1+\binom{j-i-1}{k}-\binom{k}{k}\right\\}$ (telescoping sum)
$\displaystyle=1_{k+1\leq j-i}\binom{j-i-1}{k}.$

This proves the claim for $k+1$. The result then follows by induction.

Note now that $|U_{ij}|\leq LV_{ij}$ for all $i,j\in[d]$ by hypothesis. Again
applying the triangle inequality and the assumption $\max_{i<j}|A_{ij}|\leq
L$, it follows that

$\displaystyle|A^{-1}_{ij}|$
$\displaystyle\leq\sum_{k=1}^{d-1}M^{k+1}|(U^{k})_{ij}|$
$\displaystyle\leq\sum_{k=1}^{d-1}M^{k+1}|((LV)^{k})_{ij}|$
$\displaystyle=\sum_{k=1}^{d-1}M^{k+1}L^{k}\binom{j-i-1}{k-1}1_{[k\leq j-i]}$
$\displaystyle=\sum_{k=1}^{j-i}M^{k+1}L^{k}\binom{j-i-1}{k-1}$
$\displaystyle=M^{2}L\sum_{\ell=0}^{j-i-1}(ML)^{\ell}1^{j-i-1-\ell}\binom{j-i-1}{\ell}$
($\ell:=k-1$) $\displaystyle=M^{2}L(ML+1)^{j-i-1}.$ (binomial formula)

This completes the proof. ∎

We now use this lemma to establish strong consistency of the direct map
estimator.

###### Proof of Theorem A.10.

First note that convexity of $\mathcal{X}$ along with positivity $f>0$ on
$\mathcal{X}$ implies that the inverse KR map $T^{*}=(S^{*})^{-1}$ is well-
defined and continuous, as noted in Lemma A.5. Furthermore, by definition of
$\mathcal{T}(s,d,M^{*})$, the inverse maps $T^{n}=(S^{n})^{-1}$ exist and are
continuous also.

We first prove that the sequence $\\{T^{n}\\}_{n=1}^{\infty}$ is uniformly
equicontinuous with respect to the $\ell_{\infty}$ norm $\|\cdot\|_{\infty}$
on $\mathbb{R}^{d}$. For a function $S:\mathcal{X}\to\mathcal{Y}$ let $JS(x)$
denote the Jacobian matrix at $x\in\mathcal{X}$. For a matrix
$A:\mathbb{R}^{d}\to\mathbb{R}^{d}$ let $\|A\|_{p}$ denote the operator norm
induced by the $\ell_{p}$ norm on $\mathbb{R}^{d}$, i.e.,

$\|A\|_{p}=\sup\left\\{\|Ax\|_{p}:x\in\mathbb{R}^{d},\|x\|_{p}=1\right\\}.$

When $p=\infty$, this norm is simply the maximum absolute row sum of the
matrix:

$\|A\|_{\infty}=\max_{1\leq i\leq d}\sum_{j=1}^{d}|A_{ij}|.$

We aim to bound $\|JT^{n}(y)\|_{\infty}$ uniformly in
$y\in\mathcal{Y},n\in\mathbb{N}$. Since
$\\{S^{n}\\}\subset\mathcal{T}(s,d,M^{*})$, the Jacobian matrix $JS^{n}(x)$ is
upper triangular for every $x\in\mathcal{X}$ and

$\max_{i<j}\sup_{x\in\mathcal{X}}|[JS^{n}(x)]_{ij}|=\max_{i<j}\|D_{j}S^{n}_{i}\|_{\infty}\leq
M^{*}.$

Similarly, we also have that

$\min_{j\in[d]}\inf_{x\in\mathcal{X}}|[JS^{n}(x)]_{jj}|=\min_{j\in[d]}\inf_{x\in\mathcal{X}}|D_{j}S^{n}_{j}(x)|\geq
1/M^{*}.$

It follows that for each $x\in\mathcal{X}$, the Jacobian $JS^{n}(x)$ satisfies
the hypotheses of Lemma A.11 with $L=M^{*}$. Applying the inverse function
theorem, we conclude from the lemma that the entries of
$JT^{n}(y)=(JS^{n}(T^{n}(y)))^{-1}$ are bounded as

$\displaystyle\max_{i<j}\sup_{y\in\mathcal{Y}}|[JT^{n}(y)]_{ij}|$
$\displaystyle\leq(M^{*})^{3}((M^{*})^{2}+1)^{j-i-1},$
$\displaystyle\max_{j\in[d]}\sup_{y\in\mathcal{Y}}|[JT^{n}(y)]_{jj}|$
$\displaystyle\leq M^{*}.$

It follows that the $\ell_{\infty}$ operator norm is bounded as

$\displaystyle\sup_{y\in\mathcal{Y}}\|JT^{n}(y)\|_{\infty}\leq$
$\displaystyle\ M^{*}+\sum_{j=2}^{d}(M^{*})^{3}((M^{*})^{2}+1)^{j-2}$
$\displaystyle=$
$\displaystyle\begin{cases}M^{*}+(M^{*})^{3}\cdot\frac{1-((M^{*})^{2}+1)^{d-2}}{1-((M^{*})^{2}+1)},&d\geq
2,\\\ M^{*},&d=1.\end{cases}$ (partial geometric series) $\displaystyle=$
$\displaystyle\begin{cases}M^{*}((M^{*})^{2}+1)^{d-2},&d\geq 2,\\\
M^{*},&d=1.\end{cases}$ $\displaystyle=$ $\displaystyle\
M^{*}\cdot\max\\{1,((M^{*})^{2}+1)^{d-2}\\}$ $\displaystyle:=$ $\displaystyle\
C(d,M^{*}).$

Here we have used the convention that $\sum_{j=2}^{1}a_{j}=0$ for any sequence
$\\{a_{j}\\}$. Now we apply the mean value inequality for vector-valued
functions to deduce that the $\\{T^{n}\\}$ are uniformly equicontinuous.
Indeed, for any $y_{1},y_{2}\in\mathcal{Y}$, we have

$\displaystyle\|T^{n}(y_{1})-T^{n}(y_{2})\|_{\infty}$
$\displaystyle\leq\sup_{y\in\mathcal{Y}}\|JT^{n}(y)\|_{\infty}\|y_{1}-y_{2}\|_{\infty}$
$\displaystyle\leq C(d,M^{*})\|y_{1}-y_{2}\|_{\infty}.$

Now let $y\in\mathcal{Y}$ and note that $y=S^{*}(x)$ for some
$x\in\mathcal{X}$. We then have

$\displaystyle\|T^{n}(y)-T^{*}(y)\|_{\infty}$
$\displaystyle=\|T^{n}(y)-T^{*}(S^{*}(x))\|_{\infty}$
$\displaystyle=\|T^{n}(y)-x\|_{\infty}$
$\displaystyle=\|T^{n}(S^{*}(x))-T^{n}(S^{n}(x))\|_{\infty}$
$\displaystyle\leq C(d,M^{*})\|S^{*}(x)-S^{n}(x)\|_{\infty}$
$\displaystyle\leq C(d,M^{*})\|S^{*}-S^{n}\|_{\infty,d}.$

Since $y\in\mathcal{Y}$ was arbitrary, we can take the supremum over $y$ on
the left side to obtain the desired result:

$\|T^{n}-T^{*}\|_{\infty,d}\leq
C(d,M^{*})\|S^{*}-S^{n}\|_{\infty,d}\stackrel{{\scriptstyle p}}{{\to}}0.$

∎

###  A.5 Sobolev rates under log-concavity

####  A.5.1 Proof of Theorem 4.1

Suppose the source density $g$ is log-concave, which implies that $S\mapsto
P\psi_{S}$ and $S\mapsto P_{n}\psi_{S}$ are strictly convex functionals. Since
$\mathcal{T}(s,d,M)$ is convex,

$\min_{S\in\mathcal{T}(s,d,M)}P\psi_{S},\qquad\min_{S\in\mathcal{T}(s,d,M)}P_{n}\psi_{S}$

are convex optimization problems. If in addition $g$ is strongly log-concave,
we obtain strong convexity of the objective, as we show in Lemma A.12.

######  Lemma A.12.

Suppose Assumptions 3.1-3.3 hold. Assume further that the source density $g$
is $m$-strongly log-concave for some $m>0$:

$[\nabla\log g(y_{1})-\nabla\log g(y_{2})]^{T}(y_{1}-y_{2})\leq
m\|y_{1}-y_{2}\|_{2}^{2}\qquad\forall y_{1},y_{2}\in\mathcal{Y}.$

Then the map $S\mapsto P\psi_{S}$ on $\mathcal{T}(s,d,M)$ is
$\min\\{m,M^{-2}\\}$-strongly convex with respect to the $L^{2}$ Sobolev-type
norm

$\displaystyle\|S\|_{H^{1,2}_{f}(\mathcal{X})}^{2}:=$ $\displaystyle\
\|S\|^{2}_{L^{2}_{f}(\mathcal{X})}+\sum_{k=1}^{d}\|D_{k}S_{k}\|_{L^{2}_{f}(\mathcal{X})}^{2}.$

###### Proof.

We first calculate the Gâteaux derivative of $S\mapsto P\psi_{S}$ in the
direction $A\in\mathcal{T}(s,d,M)$.

$\displaystyle\nabla P\psi_{S}(A)$ $\displaystyle=\lim_{t\to
0}\frac{P\psi_{S+tA}-P\psi_{S}}{t}$ $\displaystyle=\lim_{t\to
0}-t^{-1}\mathbb{E}\bigg{\\{}\left[\log g((S+tA)(X))-\log g(S(X))\right]$
$\displaystyle\hskip 85.35826pt+\sum_{k=1}^{d}\left[\log
D_{k}(S_{k}+tA_{k})(X)-\log D_{k}S_{k}(X)\right]\bigg{\\}}$
$\displaystyle=-\mathbb{E}\left\\{\nabla\log
g(S(X))^{T}A(X)+\sum_{k}\frac{D_{k}A_{k}(X)}{D_{k}S_{k}(X)}\right\\}.$

We can differentiate under the integral by the dominated convergence theorem,
since the integrand is smooth and compactly supported by Assumptions 3.1-3.3.

Now note that $\nabla P\psi_{S}(A)$ is a bounded linear operator in $A$.
Furthermore, since the KR map $S^{*}$ is the global minimizer of $S\mapsto
P\psi_{S}$, we have $\nabla P\psi_{S^{*}}(A)=0$ for all
$A\in\mathcal{T}(s,d,M)$ satisfying $S^{*}+tA\in\mathcal{T}(s,d,M)$ for all
$t$ sufficiently small. We now check the strong convexity condition. Assume
$A,B\in\mathcal{T}(s,d,M)$ for some $M>0$. We have

$\displaystyle(\nabla P\psi_{A}-\nabla$ $\displaystyle P\psi_{B})(A-B)=\nabla
P\psi_{A}(A-B)-\nabla P\psi_{B}(A-B)$ $\displaystyle=\mathbb{E}\\{[\nabla\log
g(B(X))-\nabla\log g(A(X))]^{T}(A(X)-B(X))$
$\displaystyle\qquad+\sum_{k}\left[\frac{D_{k}(A_{k}-B_{k})(X)}{D_{k}B_{k}(X)}-\frac{D_{k}(A_{k}-B_{k})(X)}{D_{k}A_{k}(X)}\right]\\}$
$\displaystyle\geq\mathbb{E}\left\\{m\|A(X)-B(X)\|_{2}^{2}+\sum_{k}\left[\frac{D_{k}(A_{k}-B_{k})(X)}{D_{k}B_{k}(X)}-\frac{D_{k}(A_{k}-B_{k})(X)}{D_{k}A_{k}(X)}\right]\right\\}$
(strong log-concavity of $g$)
$\displaystyle=\mathbb{E}\left\\{m\|A(X)-B(X)\|_{2}^{2}+\sum_{k}\frac{[D_{k}(A_{k}-B_{k})(X)]^{2}}{D_{k}A_{k}(X)D_{k}B_{k}(X)}\right\\}$
$\displaystyle\geq\mathbb{E}\left\\{m\|A(X)-B(X)\|_{2}^{2}+\frac{1}{M^{2}}\sum_{k}[D_{k}(A_{k}-B_{k})(X)]^{2}\right\\}$
($A,B\in\mathcal{T}(s,d,M)$)
$\displaystyle=m\|A-B\|_{L^{2}(f)}^{2}+\frac{1}{M^{2}}\sum_{k}\|D_{k}(A_{k}-B_{k})\|_{L^{2}(f)}^{2}$
$\displaystyle\geq\min\\{m,M^{-2}\\}\|A-B\|_{H^{1}(f)}^{2}.$

Hence, $S\mapsto P\psi_{S}$ satisfies the first-order strong convexity
condition. ∎

We now proceed to prove Theorem 4.1.

###### Proof of Theorem 4.1.

By Lemma A.12, strong convexity of $S\mapsto P\psi_{S}$ with respect to
$\|\cdot\|_{H^{1,2}_{f}(\mathcal{X})}$ implies that

$\displaystyle P\psi_{S^{n}}$ $\displaystyle=P\psi_{S^{n}}-P\psi_{S^{*}}$
$\displaystyle\geq\nabla
P\psi_{S^{*}}(S^{n}-S^{*})+\frac{\min\\{m,(M^{*})^{-2}\\}}{2}\|S^{n}-S^{*}\|_{H^{1,2}_{f}(\mathcal{X})}^{2}$
$\displaystyle=\frac{\min\\{m,(M^{*})^{-2}\\}}{2}\|S^{n}-S^{*}\|_{H^{1,2}_{f}(\mathcal{X})}^{2},$

since $\nabla P\psi_{S^{*}}(S^{n}-S^{*})=0$, as $S^{*}$ minimizes $P\psi_{S}$.
We complete the proof by appealing to the bound on $\mathbb{E}[P\psi_{S^{n}}]$
established in Theorem 3.5. ∎

####  A.5.2 Sobolev rates for the inverse map

Now we prove a rate of convergence of the inverse map estimator in the $L^{2}$
Sobolev norm $\|\cdot\|_{H^{1,2}_{g}(\mathcal{Y})}$ assuming strong log-
concavity, as in Theorem 4.1.

######  Theorem A.13.

Suppose Assumptions 3.1-3.3 hold. Assume further that $g$ is $m$-strongly log-
concave. Let $S^{n}$ be a near-optimizer of the functional $S\mapsto
P_{n}\psi_{S}$ on $\mathcal{T}(s,d,M^{*})$ with remainder $R_{n}$ satisfying

$\mathbb{E}[R_{n}]=\mathbb{E}\left\\{P_{n}\psi_{S^{n}}-\inf_{S\in\mathcal{T}(s,d,M^{*})}P_{n}\psi_{S}\right\\}\lesssim\mathbb{E}\|P_{n}-P\|_{\Psi(s,d,M^{*})}.$

Then $T^{n}=(S^{n})^{-1}$ converges to $T^{*}=(S^{*})^{-1}$ with respect to
the norm $\|\cdot\|_{H^{1,2}_{g}(\mathcal{Y})}$ with rate

$\frac{\min\\{m,(M^{*})^{-2}\\}}{2}\mathbb{E}\|T^{n}-T^{*}\|_{H^{1}(g)}^{2}\lesssim\begin{cases}n^{-1/2},&d<2s,\\\
n^{-1/2}\log n,&d=2s,\\\ n^{-s/d},&d>2s.\end{cases}$

###### Proof.

We aim to bound $\|T^{n}-T^{*}\|_{H^{1,2}_{g}(\mathcal{Y})}$ by a multiple of
$\|S^{n}-S^{*}\|_{H^{1,2}_{f}(\mathcal{X})}$, which will establish the same
rate of convergence for the inverse map (up to constant factors) as derived
for $S^{n}$ in Theorem 4.1.

In the proof of Theorem A.10 we showed that we can bound the $\ell_{\infty}$
matrix norm of the Jacobian $JT^{n}(y)$ uniformly as

$\sup_{y\in\mathcal{Y}}\|JT^{n}(y)\|_{\infty}\leq
C(d,M^{*}):=M^{*}\cdot\max\\{1,((M^{*})^{2}+1)^{d-2}\\}.$

Since $JT^{n}(y)$ is upper triangular, we can use the exact same argument to
arrive at the same bound on the $\ell_{1}$ matrix norm, which equals the
maximum absolute column sum of the matrix,

$\|A\|_{1}=\max_{1\leq j\leq d}\sum_{i=1}^{d}|A_{ij}|.$

Hence we conclude that $\sup_{y\in\mathcal{Y}}\|JT^{n}(y)\|_{1}\leq
C(d,M^{*})$. Now we apply Hölder’s inequality for matrix norms to conclude
that

$\sup_{y\in\mathcal{Y}}\|JT^{n}(y)\|_{2}\leq\sup_{y\in\mathcal{Y}}\sqrt{\|JT^{n}(y)\|_{1}\|JT^{n}(y)\|_{\infty}}\leq
C(d,M^{*}).$

Consequently, we can bound the first term in
$\|\cdot\|_{H^{1,2}_{g}(\mathcal{Y})}$ as

$\displaystyle\|T^{n}-T^{*}\|_{L^{2}_{g}(\mathcal{Y})}^{2}$
$\displaystyle=\int\|T^{n}(y)-T^{*}(y)\|_{2}^{2}\ g(y)dy$
$\displaystyle=\int\|T^{n}(S^{*}(x))-T^{*}(S^{*}(x))\|_{2}^{2}\ |\det
JS^{*}(x)|g(S^{*}(x))dx$ (Change of variables)
$\displaystyle=\int\|T^{n}(S^{*}(x))-T^{*}(S^{*}(x))\|_{2}^{2}\ f(x)dx$
($g=T_{\\#}f$) $\displaystyle=\int\|T^{n}(S^{*}(x))-T^{n}(S^{n}(x))\|_{2}^{2}\
f(x)dx$
$\displaystyle\leq\int\sup_{y\in\mathcal{Y}}\|JT^{n}(y)\|_{2}^{2}\|S^{*}(x)-S^{n}(x)\|_{2}^{2}\
f(x)dx$ $\displaystyle\leq C(d,M^{*})^{2}\int\|S^{*}(x)-S^{n}(x)\|_{2}^{2}\
f(x)dx$
$\displaystyle=C(d,M^{*})^{2}\|S^{*}-S^{n}\|_{L^{2}_{f}(\mathcal{X})}^{2}.$

To bound the deviations of the first derivatives, note that

$\displaystyle|D_{k}T^{n}_{k}(y)-D_{k}T^{*}_{k}(y)|$
$\displaystyle=\left|\frac{1}{D_{k}S_{k}^{n}(T^{n}(y))}-\frac{1}{D_{k}S^{*}_{k}(T^{*}(y))}\right|$
$\displaystyle=\left|\frac{D_{k}S_{k}^{n}(T^{n}(y))-D_{k}S^{*}_{k}(T^{*}(y))}{[D_{k}S_{k}^{n}(T^{n}(y))][D_{k}S^{*}_{k}(T^{*}(y))]}\right|$
$\displaystyle\leq(M^{*})^{2}\left|D_{k}S_{k}^{n}(T^{n}(y))-D_{k}S^{*}_{k}(T^{*}(y))\right|$
($S,S^{n}\in\mathcal{T}(s,d,M^{*})$)
$\displaystyle\leq(M^{*})^{2}\left|D_{k}S_{k}^{n}(T^{n}(y))-D_{k}S_{k}^{n}(T^{*}(y))\right|$
$\displaystyle\quad+(M^{*})^{2}\left|D_{k}S_{k}^{n}(T^{*}(y))-D_{k}S^{*}_{k}(T^{*}(y))\right|$
$\displaystyle\leq(M^{*})^{2}\sup_{x\in\mathcal{X}}\|\nabla(D_{k}S_{k}^{n})(x)\|_{2}\|T^{n}(y)-T^{*}(y)\|_{2}$
$\displaystyle\quad+(M^{*})^{2}\left|D_{k}S_{k}^{n}(T^{*}(y))-D_{k}S^{*}_{k}(T^{*}(y))\right|.$
(mean value inequality)

Now note that $S_{k}^{n}$ depends only on $(x_{k},\ldots,x_{d})$, and
therefore, since $S^{n}\in\mathcal{T}(s,d,M^{*})$ implies that
$\max_{j\in[d]}\|D_{j}D_{k}S_{k}^{n}\|_{\infty}\leq M^{*}$, we have

$\displaystyle\sup_{x\in\mathcal{X}}\|\nabla(D_{k}S_{k}^{n})(x)\|_{2}$
$\displaystyle=\sup_{x\in\mathcal{X}}\sqrt{\sum_{j=1}^{d}|D_{j}D_{k}S_{k}^{n}(x)|^{2}}$
$\displaystyle=\sup_{x\in\mathcal{X}}\sqrt{\sum_{j=k}^{d}|D_{j}D_{k}S_{k}^{n}(x)|^{2}}$
$\displaystyle\leq\sqrt{\sum_{j=k}^{d}(M^{*})^{2}}$
$\displaystyle=M^{*}\sqrt{d-k+1}.$

Hence, we conclude that

$\displaystyle|D_{k}T^{n}_{k}(y)-D_{k}T^{*}_{k}(y)|^{2}$
$\displaystyle\leq\\{(M^{*})^{3}\sqrt{d-k+1}\|T^{n}(y)-T^{*}(y)\|_{2}$
$\displaystyle\qquad+(M^{*})^{2}\left|D_{k}S_{k}^{n}(T^{*}(y))-D_{k}S^{*}_{k}(T^{*}(y))\right|\\}^{2}$
$\displaystyle\leq 2(M^{*})^{6}(d-k+1)\|T^{n}(y)-T^{*}(y)\|_{2}^{2}$
$\displaystyle\qquad+2(M^{*})^{4}\left|D_{k}S_{k}^{n}(T^{*}(y))-D_{k}S^{*}_{k}(T^{*}(y))\right|^{2}.$
($(a+b)^{2}\leq 2(a^{2}+b^{2})$)

Summing over $k$ and integrating against the density $g$, we obtain

$\displaystyle\sum_{k=1}^{d}\|D_{k}T^{n}_{k}-$ $\displaystyle
D_{k}T^{*}_{k}\|_{L^{2}_{g}(\mathcal{Y})}^{2}=\sum_{k=1}^{d}\int|D_{k}T^{n}_{k}(y)-D_{k}T^{*}_{k}(y)|^{2}g(y)dy$
$\displaystyle\leq\sum_{k=1}^{d}\int
2(M^{*})^{6}(d-k+1)\|T^{n}(y)-T^{*}(y)\|_{2}^{2}\ g(y)dy$
$\displaystyle\quad+\sum_{k=1}^{d}\int
2(M^{*})^{4}\left|D_{k}S_{k}^{n}(T^{*}(y))-D_{k}S^{*}_{k}(T^{*}(y))\right|^{2}g(y)dy$
$\displaystyle=(M^{*})^{6}d(d+1)\|T^{n}-T^{*}\|^{2}_{L^{2}_{g}(\mathcal{Y})}$
$\displaystyle\quad+2(M^{*})^{4}\sum_{k=1}^{d}\int\left|D_{k}S_{k}^{n}(x)-D_{k}S^{*}_{k}(x)\right|^{2}g(S^{*}(x))|\det
JS^{*}(x)|dx$
$\displaystyle=(M^{*})^{6}d(d+1)\|T^{n}-T^{*}\|^{2}_{L^{2}_{g}(\mathcal{Y})}$
$\displaystyle\quad+2(M^{*})^{4}\sum_{k=1}^{d}\int\left|D_{k}S_{k}^{n}(x)-D_{k}S^{*}_{k}(x)\right|^{2}f(x)dx$
$\displaystyle=(M^{*})^{6}d(d+1)\|T^{n}-T^{*}\|^{2}_{L^{2}_{g}(\mathcal{Y})}+2(M^{*})^{4}\sum_{k=1}^{d}\|D_{k}(S_{k}^{n}-S^{*}_{k})\|^{2}_{L^{2}_{f}(\mathcal{X})}$
$\displaystyle\leq(M^{*})^{6}d(d+1)C(d,M^{*})^{2}\|S^{n}-S^{*}\|^{2}_{L^{2}_{f}(\mathcal{X})}+2(M^{*})^{4}\sum_{k=1}^{d}\|D_{k}(S_{k}^{n}-S^{*}_{k})\|^{2}_{L^{2}_{f}(\mathcal{X})}$

Putting all of these calculations together, we have shown that

$\displaystyle\|T^{n}-T^{*}\|^{2}_{H^{1,2}_{g}(\mathcal{Y})}$
$\displaystyle=\|T^{n}-T^{*}\|^{2}_{L^{2}_{g}(\mathcal{Y})}+\sum_{k=1}^{d}\|D_{k}(T^{n}_{k}-T^{*}_{k})\|^{2}_{L^{2}_{g}(\mathcal{Y})}$
$\displaystyle\leq
C(d,M^{*})^{2}\|S^{n}-S^{*}\|^{2}_{L^{2}_{f}(\mathcal{X})}+(M^{*})^{6}d(d+1)C(d,M^{*})^{2}\|S^{n}-S^{*}\|^{2}_{L^{2}_{f}(\mathcal{X})}$
$\displaystyle\quad+2(M^{*})^{4}\sum_{k=1}^{d}\|D_{k}(S_{k}^{n}-S^{*}_{k})\|^{2}_{L^{2}_{f}(\mathcal{X})}$
$\displaystyle\leq\tilde{C}(d,M^{*})^{2}\|S^{n}-S^{*}\|^{2}_{H^{1,2}_{f}(\mathcal{X})},$

where we define

$\tilde{C}(d,M)^{2}=\max\\{C(d,M)^{2}[1+M^{6}d(d+1)],2M^{4}\\}.$

Finally, we appeal to the bound on
$\|S^{n}-S^{*}\|^{2}_{H^{1,2}_{f}(\mathcal{X})}$ derived in Theorem 4.1 to
conclude the proof:

$\displaystyle\frac{\min\\{m,(M^{*})^{-2}\\}}{2}\mathbb{E}\|T^{n}-T^{*}\|^{2}_{H^{1,2}_{g}(\mathcal{Y})}$
$\displaystyle\leq\frac{\min\\{m,(M^{*})^{-2}\\}}{2}\tilde{C}(d,M^{*})^{2}\mathbb{E}\|S^{n}-S^{*}\|_{H^{1,2}_{f}(\mathcal{X})}^{2}$
$\displaystyle\lesssim\begin{cases}n^{-1/2},&d<2s,\\\ n^{-1/2}\log n,&d=2s,\\\
n^{-s/d},&d>2s.\end{cases}$

∎

###  A.6 Dimension ordering

####  A.6.1 Proof of Lemma 4.2

Now we establish a rate of convergence in the anisotropic smoothness setting.
Define

$\psi_{S}^{k}(x)=\log f_{k}(x_{k}|x_{(k+1):d})-\log
g_{k}(S_{k}(x)|S_{(k+1):d}(x))+\log D_{k}S_{k}(x),$

which is $(s_{k},\ldots,s_{d})$-smooth in $(x_{k},\ldots,x_{d})$ whenever each
$S_{k}$ is $(s_{k}+1,s_{k+1},\ldots,s_{d})$-smooth in
$(x_{k},x_{k+1},\ldots,x_{d})$ by Assumption 4.1 and Lemma A.6. Since

$f(x)=\prod_{k=1}^{d}f_{k}(x_{k}|x_{(k+1):d})$

by the chain rule of densities, and similarly for $g$, we have

$(P_{n}-P)\psi_{S}=\sum_{k=1}^{d}(P_{n}-P)\psi_{S}^{k}.$

Note that $\psi_{S}^{k}$ is a function of $S_{k:d}$ and the $d_{k}$ variables
$(x_{k},\ldots,x_{d})$ only. Defining

$\Psi_{k}(\mathfrak{s},d,M)=\\{\psi_{S}^{k}:S\in\mathcal{T}(\mathfrak{s},d,M)\\},$

we have the following analog of Proposition 3.3.

######  Proposition A.14.

Under Assumptions 3.1, 3.2, and 4.1, the metric entropy of
$\Psi_{k}(\mathfrak{s},d,M)$ in the $L^{\infty}$ norm is bounded as

$H(\epsilon,\Psi_{k}(\mathfrak{s},d,M),\|\cdot\|_{\infty})\lesssim\epsilon^{-d_{k}/{\sigma}_{k}}.$

Consequently, $\Psi_{k}(\mathfrak{s},d,M)$ is totally bounded and therefore
precompact in $L^{\infty}(\mathcal{X}_{k:d})$.

Hence, we obtain in Lemma 4.2 a bound on the supremum process analogous to
Lemma 3.4 now applied to the family $\Psi_{k}(\mathfrak{s},d,M)$ of
$(s_{k},\ldots,s_{d})$-smooth maps $\psi_{S}^{k}$.

###### Proof of Lemma 4.2.

The metric entropy integral bounds utilized in the proof of Lemma 3.4,
combined with the entropy bound on $\Psi_{k}(\mathfrak{s},d,M)$ derived in
Proposition A.14, yield the following rate

$\mathbb{E}\|P_{n}-P\|_{\Psi_{k}(\mathfrak{s},d,M)}=\mathbb{E}\left\\{\sup_{S\in\mathcal{T}(\mathfrak{s},d,M)}|(P_{n}-P)\psi_{S}^{k}|\right\\}\lesssim
c_{n,k}.$

Applying the triangle inequality then yields

$\displaystyle\mathbb{E}\|P_{n}-P\|_{\Psi(\mathfrak{s},d,M)}$
$\displaystyle=\mathbb{E}\left\\{\sup_{S\in\mathcal{T}(\mathfrak{s},d,M)}|(P_{n}-P)\psi_{S}|\right\\}$
$\displaystyle=\mathbb{E}\left\\{\sup_{S\in\mathcal{T}(\mathfrak{s},d,M)}\left|\sum_{k=1}^{d}(P_{n}-P)\psi^{k}_{S}\right|\right\\}$
$\displaystyle\leq\sum_{k=1}^{d}\mathbb{E}\left\\{\sup_{S\in\mathcal{T}(\mathfrak{s},d,M)}\left|(P_{n}-P)\psi^{k}_{S}\right|\right\\}$
$\displaystyle=\sum_{k=1}^{d}\mathbb{E}\|P_{n}-P\|_{\Psi_{k}(\mathfrak{s},d,M)}$
$\displaystyle\lesssim\sum_{k=1}^{d}c_{n,k}.$

∎

###  A.7 Jacobian flows

###### Proof of Theorem 4.6.

The proof is practically identical to that of Theorem 3.5, since the functions
in $\Psi_{m}(\Sigma,s,M)$ are $s$-smooth with uniformly bounded derivatives
(analogous to $\Psi(s,d,M)$) by definition of $\mathcal{J}_{m}(\Sigma,s,M)$,
the chain rule of differentiation, and the relation

$\displaystyle\psi_{S}(x)=$
$\displaystyle\log[f(x)/g(S(x))]-\sum_{j=1}^{m}\sum_{k=1}^{d}\log
D_{k}U_{k}^{j}(x^{j}),$

where we define

$x^{j}=\Sigma^{j}\circ U^{j-1}\circ\Sigma^{j-1}\circ\cdots\circ
U^{1}\circ\Sigma^{1}(x),\quad j\in[m].$

Hence, the entropy estimates for $\Psi(s,d,M)$ in Proposition 3.3 hold also
for $\Psi_{m}(\Sigma,s,M)$. Thus, we obtain similar bounds on
$\mathbb{E}\|P_{n}-P\|_{\Psi_{m}(\Sigma,s,M)}$ as in Lemma 3.4. Combining this
with the risk decomposition (4) and the argument in Theorem 3.5 completes the
proof. ∎

###  A.8 On separability

Suppose the source $g$ is a product density, i.e.,
$g(y)=\prod_{k=1}^{d}g_{k}(y_{k})$ for some smooth densities
$g_{k}:\mathbb{R}\to\mathbb{R}$. As the source density $g$ is a degree of
freedom in our problem, we are free to choose $g$ to factor as such. For
example, $g$ could be the standard normal density in $d$-dimensions or the
uniform density on a box $\mathcal{Y}\subset\mathbb{R}^{d}$. We will show that
the task of estimating the KR map $S^{*}$ is amenable to distributed
computation in this case.

######  Assumption A.1.

The source density $g$ factors as a product:
$g(y)=\prod_{k=1}^{d}g_{k}(y_{k})$.

Recall the minimization objective defining our estimator:

$-\frac{1}{n}\sum_{i=1}^{n}\left[\log g(S(X^{i}))+\sum_{k=1}^{d}\log
D_{k}S_{k}(X^{i})\right],$

where $X^{i}=(X^{i}_{1},\ldots,X^{i}_{d}),i=1,\ldots,n$ is an iid random
sample from $f$. Here we omit the entropy term involving $f$ that does not
depend on $S$ without loss of generality. For a general density $g$, we can
simplify the above expression by appealing to the chain rule for densities

$g(y)=\prod_{k=1}^{d}g_{k}(y_{k}|y_{k+1},\ldots,y_{d}),$

where

$g_{k}(y_{k}|y_{k+1},\ldots,y_{d})=\frac{\int
g(y_{1},\ldots,y_{d})dy_{1}\cdots dy_{k-1}}{\int
g(y_{1},\ldots,y_{d})dy_{1}\cdots dy_{k}}.$

The objective then becomes

$-\frac{1}{n}\sum_{i=1}^{n}\sum_{k=1}^{d}\left[\log
g_{k}(S_{k}(X^{i})|S_{k+1}(X^{i}),\ldots,S_{d}(X^{i}))+\log
D_{k}S_{k}(X^{i})\right].$

When $g$ is a product density we have
$g_{k}(y_{k}|y_{k+1},\ldots,y_{d})=g_{k}(y_{k})$ and we obtain

$\sum_{k=1}^{d}\left\\{-\frac{1}{n}\sum_{i=1}^{n}\left[\log
g_{k}(S_{k}(X^{i}))+\log D_{k}S_{k}(X^{i})\right]\right\\},$

which is a separable objective over the component maps $S_{k}$. In this case,
we can find our estimator $S^{n}=(S^{n}_{1},\ldots,S^{n}_{d})$ by solving for
the components $S^{n}_{k}$ in parallel.