# Convergence of Continuous Normalizing Flows for Learning Probability
Distributions

Yuan<EMAIL_ADDRESS>[ Jian
<EMAIL_ADDRESS>[ Yuling
<EMAIL_ADDRESS>[ Shurong
<EMAIL_ADDRESS>[ Department of Applied Mathematics, The
Hong Kong Polytechnic University, Hong Kong SAR, Chinapresep=, ]e1  Department
of Applied Mathematics, The Hong Kong Polytechnic University, Hong Kong SAR,
Chinapresep=, ]e2  School of Mathematics and Statistics and Hubei Key
Laboratory of Computational Science, Wuhan University, Wuhan, Chinapresep=,
]e3  School of Mathematics and Statistics, Northeast Normal University,
Changchun, Chinapresep=, ]e4

###### Abstract

Continuous normalizing flows (CNFs) are a generative method for learning
probability distributions, which is based on ordinary differential equations.
This method has shown remarkable empirical success across various
applications, including large-scale image synthesis, protein structure
prediction, and molecule generation. In this work, we study the theoretical
properties of CNFs with linear interpolation in learning probability
distributions from a finite random sample, using a flow matching objective
function. We establish non-asymptotic error bounds for the distribution
estimator based on CNFs, in terms of the Wasserstein-2 distance. The key
assumption in our analysis is that the target distribution satisfies one of
the following three conditions: it either has a bounded support, is strongly
log-concave, or is a finite or infinite mixture of Gaussian distributions. We
present a convergence analysis framework that encompasses the error due to
velocity estimation, the discretization error, and the early stopping error. A
key step in our analysis involves establishing the regularity properties of
the velocity field and its estimator for CNFs constructed with linear
interpolation. This necessitates the development of uniform error bounds with
Lipschitz regularity control of deep ReLU networks that approximate the
Lipschitz function class, which could be of independent interest. Our
nonparametric convergence analysis offers theoretical guarantees for using
CNFs to learn probability distributions from a finite random sample.

62G05; ,

68T07,

Generative learning,

Lipschitz regularity,

non-asymptotic error bound,

neural network approximation,

velocity fields,

###### keywords:

[class=MSC]

###### keywords:

, , and

##  1 Introduction

Let $\\{\mathsf{X}_{i}\\}_{i=1}^{n}$ be independent and identically
distributed (i.i.d.) random variables drawn from an underlying probability
distribution $\nu$ with support in ${\mathbb{R}}^{d}$. The task of generative
learning is to learn $\nu$ from the data $\\{\mathsf{X}_{i}\\}_{i=1}^{n}$ by
generating new samples (Salakhutdinov, 2015). Several generative learning
methods have been developed during the recent decade, including generative
adversarial networks (GANs) (Goodfellow et al., 2014), variational
autoencoders (Kingma and Welling, 2014), diffusion models (Sohl-Dickstein et
al., 2015; Ho, Jain and Abbeel, 2020; Song et al., 2021), and normalizing
flows (Tabak and Turner, 2013; Rezende and Mohamed, 2015; Chen et al., 2018).
Deep neural networks (LeCun, Bengio and Hinton, 2015), as a powerful modeling
tool, have played an important role in the development of these methods. Among
the generative learning methods, continuous normalizing flows (CNFs) use
ordinary differential equations (ODEs) to determine a stochastic process for
transporting a Gaussian distribution to the target distribution, achieving the
goal of generative learning.

Simulation-free CNFs have achieved impressive empirical performance across
various applications. These applications include large-scale image synthesis
(Ma et al., 2024), protein structure prediction (Jing, Berger and Jaakkola,
2023), and 3D molecule generation (Song et al., 2023). Rectified flow (Liu,
Gong and Liu, 2023), a CNF model that linearly interpolates Gaussian noise and
data, has been recently implemented in the large image model Stable Diffusion
3 (Esser et al., 2024). Simulation-free CNFs that use flow matching to learn
probability distributions have been the focus of much recent attention
(Albergo and Vanden-Eijnden, 2023; Lipman et al., 2023; Liu, Gong and Liu,
2023; Neklyudov et al., 2023).

An early model of CNFs was proposed by Chen et al. (2018). This model is based
on neural ODEs and employs a simulation-based maximum likelihood method to
estimate velocity fields. However, simulation-based CNFs are computationally
demanding in large-scale applications. To address the computational challenges
of simulation-based CNFs, significant efforts have been made to develop
simulation-free CNFs, where velocity fields can be represented in terms of
conditional expectations. Noteworthy examples include the probability flows of
diffusion models (Song et al., 2021) and denoising diffusion implicit models
(Song, Meng and Ermon, 2021), which are trained using a denoising score
matching objective function. In contrast, the flow matching method solves a
least squares problem to estimate the conditional expectation that represents
the velocity field (Albergo and Vanden-Eijnden, 2023; Lipman et al., 2023;
Liu, Gong and Liu, 2023).

When a random sample from the target distribution is available, the process of
learning simulation-free CNFs involves statistically and numerically solving a
class of ODE-based initial value problems (IVPs). Specifically, consider the
IVP on the unit time interval

$\frac{\mathrm{d}X_{t}}{\mathrm{d}t}(x)=v(t,X_{t}(x)),\quad
X_{0}(x)=x\sim\mu,\quad(t,x)\in[0,1]\times{\mathbb{R}}^{d},$ (1.1)

where $v$ represents the velocity field, which can be estimated based on data.
The solution to the IVP (1.1) is a family of flow maps $(X_{t})_{t\in[0,1]}$
indexed by the time variable, which generates a smoothing path between the
source and target distributions. Gao, Huang and Jiao (2023) have studied the
mathematical properties of these ODE-based IVPs under the framework of
Gaussian interpolation flow. This defines a Gaussian smoothing path as a
transport map that pushes forward a Gaussian distribution onto the target
distribution in terms of measure transport. Simulation-free CNFs adopt a two-
step “estimation-then-simulation” approach to learning the desired transport
map based on a random sample. In the estimation stage, a deep learning model
is trained to estimate the velocity field without simulating the ODE that
defines the CNF. During the simulation stage, numerical solvers simulate the
numerical solution of the ODE associated with the estimated velocity field,
and the generated data is collected at the end time point.

In this work, we study the theoretical properties of simulation-free CNFs. We
develop a general framework for error analyses of CNFs with flow matching for
learning probability distributions based on a random sample. Central to
simulation-free CNFs, deep ReLU networks are employed for function
approximation and nonparametric estimation of the velocity field. We establish
the approximation properties of deep ReLU networks with Lipschitz regularity
control, which is essential for analyzing the impact of the estimated velocity
field on the distribution of the data generated through the flow. In
particular, it is crucial to control the Lipschitz regularity of the estimated
velocity field to ensure that the associated IVP is well-posed.

###  1.1 Preview of main results

The following informal descriptions provide a preview of our main results. Our
analysis is based on the key assumption that the target distribution satisfies
one of the three conditions: it either has a bounded support, is strongly log-
concave, or is a finite or infinite mixture of Gaussian distributions. The
assumption of bounded support is practical, as all measurements can only be
taken within a finite range. However, it is theoretically interesting to
consider distributions with unbounded support. In such cases, we impose
additional structures on the distributions, including mixtures of Gaussians
and strong log-concavity. Mixtures of Gaussian distributions are multimodal
and have a universal approximation capacity for smooth probability density
functions. The assumption of strong log-concavity is crucial for the
convergence analyses of Langevin Monte Carlo and Metropolis-Hastings
algorithms when sampling from unnormalized distributions (Dalalyan, 2017;
Durmus and Moulines, 2017; Dwivedi et al., 2019).

Our first main result concerns the regularity of the velocity fields of the
CNFs constructed with linear interpolation. For a detailed definition of such
CNFs, see Lemma 3.1 and equation (3.4).

######  Theorem 1.1 (Informal).

Assume that the target distribution either has a bounded support, is strongly
log-concave, or is a mixture of Gaussians. Let $0<\underline{t}\ll 1$. The
velocity fields of the CNFs with linear interpolation have the following
regularity properties:

  * (i)

The velocity field $v^{*}$ is Lipschitz continuous in the space variable $x$
for $(t,x)\in[0,1]\times{\mathbb{R}}^{d}$, where the Lipschitz constant is
uniformly bounded;

  * (ii)

The velocity field $v^{*}$ is Lipschitz continuous in the time variable $t$
for $(t,x)\in[0,1-\underline{t}]\times{\mathbb{R}}^{d}$, where the Lipschitz
constant grows at the order of ${\mathcal{O}}(\underline{t}^{-2})$ as
$\underline{t}\downarrow 0$;

  * (iii)

The velocity field $v^{*}$ spatially has a linear growth on ${\mathbb{R}}^{d}$
for each $t\in[0,1]$.

######  Remark 1.

The regularity properties of the velocity fields stated in Theorem 1.1 are
derived from the assumptions made about the underlying target distribution.
These properties are essential for studying the distributions generated by the
corresponding CNFs.

######  Theorem 1.2 (Informal).

Suppose that the target distribution has a bounded support, or is strongly
log-concave, or is a mixture of Gaussians. Let $n$ be the sample size and
$0<\underline{t}\ll 1$ satisfying $\underline{t}\asymp n^{-1/(d+5)}$. By
properly setting the deep ReLU network structure and the forward Euler
discretization step sizes, the distribution estimation error of the CNFs
learned with linear interpolation and flow matching is evaluated by

$\displaystyle\mathbb{E}\mathcal{W}_{2}(\hat{\nu}_{1-\underline{t}},\nu)=\widetilde{\mathcal{O}}(n^{-\frac{1}{d+5}}),$
(1.2)

where the expectation is taken with respect to all random samples,
$\hat{\nu}_{1-\underline{t}}$ is the law of generated data, $\nu$ is the law
of target data, $\mathcal{W}_{2}(\cdot,\cdot)$ is the Wasserstein-2 distance,
and a polylogarithmic prefactor in $n$ is omitted.

######  Remark 2.

As can be seen from Theorem 1.1 or Theorem 5.1 below, the velocity fields
associated with the CNFs based on linear interpolation may be singular in the
time variable at $t=1$ due to the exploding Lipschitz constant bound. This
singularity affects the convergence rate in (1.2). Without the time
singularity of the velocity field, the distribution estimation error would be
bounded by $\widetilde{{\mathcal{O}}}(n^{-1/(d+3)})$. The time singularity
leads to a necessary trade-off regarding $\underline{t}$ between the error due
to velocity estimation and the early stopping error. The trade-off reduces the
nonparametric convergence rate of the distribution estimator to
$\widetilde{{\mathcal{O}}}(n^{-1/(d+5)})$. See also Remark 10 for additional
comments and explanation.

###  1.2 Our contributions

We present a comprehensive error analysis of simulation-free CNFs with linear
interpolation, trained using flow matching. To the best of our knowledge, this
is the first analysis of its kind in the context of simulation-free CNFs. Our
results are based solely on assumptions about the target distribution, and all
regularity conditions are rigorously derived from these assumptions. Our
analysis accurately reflects the practical computational implementation of
flow matching for learning simulation-free CNFs. Although our focus is on CNFs
based on linear interpolation due to their widespread use and for the sake of
simplicity, our analytical framework can be applied to CNFs based on other
types of interpolation as well.

We summarize our main contributions into four points:

  * (1)

We establish non-asymptotic error bounds for distribution estimators based on
simulation-free CNFs with linear interpolation and flow matching. These bounds
apply to distributions that satisfy one of the following conditions: (a) they
have a bounded support, (b) they are strongly log-concave, or (c) they are
mixtures of Gaussians. We present a convergence analysis framework that
encompasses the error due to velocity estimation, the discretization error,
and the early stopping error. We show that the nonparametric convergence rate
of the distribution estimator is $\widetilde{{\mathcal{O}}}(n^{-1/(d+5)})$ up
to a polylogarithmic prefactor in the sample size $n$.

  * (2)

We derive regularity properties for the velocity field of the CNF with linear
interpolation. We demonstrate the Lipschitz regularity properties of the
velocity field in both the space and time variables and establish bounds for
the Lipschitz constants. We also show that the velocity field grows at most
linearly with respect to the space variable.

  * (3)

We establish error bounds for deep ReLU network approximation within the
Lipschitz function class, demonstrating that the constructed approximation
function maintains Lipschitz regularity. We also derive time-space
approximation bounds for approximating the velocity field in both the time and
space variables. These time-space approximation bounds are novel in three
respects. Firstly, the approximation bounds are derived in terms of the
$L^{\infty}$ norm. Secondly, we demonstrate that the constructed time-space
approximation function is Lipschitz in both the time and space variables, with
Lipschitz constants of the same order as those of the target function. Lastly,
the time-space approximation function can exhibit different Lipschitz
regularity in the time and space variables. These neural network approximation
results, which maintain Lipschitz regularity, could be of independent
interest.

  * (4)

We establish the statistical consistency of the flow matching estimator for
the velocity field. By rigorously bounding the stochastic and approximation
errors, we show that the convergence rate of the flow matching estimator
coincides with the minimax optimal rate of nonparametric estimation of
regression functions belonging to the Sobolev space $W^{1,\infty}([0,1]^{d})$.

The remainder of this paper is organized as follows. In Section 2, we present
the preliminary materials required for subsequent sections. In Section 3, we
describe simulation-free CNFs and outline the steps for using these CNFs for
generative learning. In Section 4, we derive our main result concerning the
error bounds for the distribution estimator based on CNFs with linear
interpolation. In Section 5, we first present some useful regularity
properties of the velocity field and establish error bounds for the estimated
velocity fields through flow matching. Section 6 contains discussions on
related works in the existing literature. We conclude with remarks and
discussions on potential future problems in Section 7.

##  2 Preliminaries

In this section, we summarize several useful definitions, assumptions, the
background of CNFs and deep ReLU networks.

### Notation

Let $\mathbb{N}_{0}$ and $\mathbb{N}$ denote the set of non-negative integers
and the set of positive integers, respectively, that is,
$\mathbb{N}=\\{1,2,3,\cdots\\}$ and $\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$. We
use $\|x\|_{p}$ to denote the $\ell_{p}$-norm of a vector
$x\in{\mathbb{R}}^{d}$ for $p\in[1,\infty]$. Let $\mathrm{I}_{d}$ denote the
$d\times d$ identity matrix. For
$\Omega_{1}\subseteq{\mathbb{R}}^{k},\Omega_{2}\subseteq{\mathbb{R}}^{d},n\in\mathbb{N}$,
we denote by $C^{n}(\Omega_{1};\Omega_{2})$ the space of continuous functions
$f:\Omega_{1}\to\Omega_{2}$ that are $n$ times differentiable and whose
partial derivatives of order $n$ are continuous. For any $f(x)\in
C^{2}({\mathbb{R}}^{d};{\mathbb{R}})$, let $\nabla_{x}f$ and $\dot{f}$ denote
its gradient, and we use $\nabla^{2}_{x}f$ to denote its Hessian. For
simplicity, let $\gamma_{d,\sigma^{2}}:=N(0,\sigma^{2}\mathrm{I}_{d})$ and
$\gamma_{d}:=N(0,\mathrm{I}_{d})$. Let $f:{\mathbb{R}}^{k}\to{\mathbb{R}}^{d}$
be a measurable mapping and $\mu$ be a probability distribution on
${\mathbb{R}}^{k}$. The push-forward distribution $f_{\\#}\mu$ of a measurable
set $A$ is defined as $f_{\\#}\mu:=\mu(f^{-1}(A))$. We use $X\lesssim Y$ and
$Y\gtrsim X$ to denote $X\leq CY$ for some constant $C>0$. The notation
$X\asymp Y$ indicates that $X\lesssim Y\lesssim X$. For a set
$\Omega\subset{\mathbb{R}}^{d}$, let
$\Omega^{c}:=\\{x\in{\mathbb{R}}^{d}:x\notin\Omega\\}$, and we use
$\mathds{1}_{\Omega}:{\mathbb{R}}^{d}\to\\{0,1\\}$ to denote the indicator
function of $\Omega$.

###  2.1 Definitions

We list a few useful definitions in this subsection.

######  Definition 2.1 (Cattiaux and Guillin, 2014).

A probability distribution $\mu(\mathrm{d}x)=\exp(-U)\mathrm{d}x$ is
$\kappa$-semi-log-concave for some $\kappa\in{\mathbb{R}}$ if its support
$\Omega\subseteq{\mathbb{R}}^{d}$ is convex and its potential function $U\in
C^{2}(\Omega;{\mathbb{R}})$ satisfies

$\nabla^{2}_{x}U(x)\succeq\kappa\mathrm{I}_{d},\quad\forall x\in\Omega.$

Note that here $\kappa$ is allowed to be negative.

######  Definition 2.2 (Eldan and Lee, 2018).

A probability distribution $\mu(\mathrm{d}x)=\exp(-U)\mathrm{d}x$ is
$\beta$-semi-log-convex for some $\beta>0$ if its support
$\Omega\subseteq{\mathbb{R}}^{d}$ is convex and its potential function $U\in
C^{2}(\Omega;{\mathbb{R}})$ satisfies

$\nabla^{2}_{x}U(x)\preceq\beta\mathrm{I}_{d},\quad\forall x\in\Omega.$

The rectified linear unit (ReLU) activation function is defined as
$\varrho(x):=\max(0,x)$ for $x\in{\mathbb{R}},$ which also operates
coordinate-wise on elements of $x\in{\mathbb{R}}^{d}$. For $k\geq 2$, the
ReLUk activation function is given by $\varrho_{k}(x):=(\max(0,x))^{k}$.

######  Definition 2.3 (Deep ReLU networks).

Deep ReLU networks stand for a class of feed-forward artiﬁcial neural networks
defined with ReLU activation functions. The function
$f_{\theta}(x):{\mathbb{R}}^{\mathtt{k}}\to{\mathbb{R}}^{\mathtt{d}}$
implemented by a deep ReLU network with parameter $\theta$ is expressed as
composition of a sequence of functions

$\displaystyle f_{\theta}(x):=l_{\mathtt{D}}\circ\varrho\circ
l_{\mathtt{D}-1}\circ\varrho\circ\cdots\circ l_{1}\circ\varrho\circ l_{0}(x)$

for any $x\in{\mathbb{R}}^{\mathtt{k}}$, where $\varrho(x)$ is the ReLU
activation function and the depth $\mathtt{D}$ is the number of hidden layers.
For $i=0,1,\cdots,\mathtt{D}$, the $i$-th layer is represented by
$l_{i}(x):=A_{i}x+b_{i}$, where $A_{i}\in{\mathbb{R}}^{d_{i+1}\times d_{i}}$
is the weight matrix, $b_{i}\in{\mathbb{R}}^{d_{i+1}}$ is the bias vector, and
$d_{i}$ is the width of the $i$-th layer. The network $f_{\theta}$ contains
$\mathtt{D}+1$ layers in all. We use a $(\mathtt{D}+1)$-dimension vector
$(d_{0},d_{1},\cdots,d_{\mathtt{D}})^{\top}$ to describe the width of each
layer. In particular, $d_{0}=\mathtt{k}$ is the dimension of the domain and
$d_{\mathtt{D}}=\mathtt{d}$ is the dimension of the codomain. The width
$\mathtt{W}$ is defined as the maximum width of hidden layers, that is,
$\mathtt{W}=\max\left\\{d_{1},d_{2},\cdots,d_{\mathtt{D}}\right\\}$. The size
$\mathtt{S}$ denotes the total number of nonzero parameters in the network
$f_{\theta}$. The bound $\mathtt{B}$ denotes the $L^{\infty}$ bound of
$f_{\theta}$, that is,
$\sup_{x\in{\mathbb{R}}^{\mathtt{k}}}\|f_{\theta}(x)\|_{\infty}\leq\mathtt{B}$.
We denote the function class
$\\{f_{\theta}:{\mathbb{R}}^{\mathtt{k}}\to{\mathbb{R}}^{\mathtt{d}}\\}$
implemented by deep ReLU networks with size $\mathtt{S}$, width $\mathtt{W}$,
depth $\mathtt{D}$, and bound $\mathtt{B}$ as
$\mathcal{NN}(\mathtt{S},\mathtt{W},\mathtt{D},\mathtt{B},\mathtt{k},\mathtt{d})$.

######  Definition 2.4 (Wasserstein-$2$ distance).

The Wasserstein-$2$ distance between two probability distributions on
${\mathbb{R}}^{d}$ is the $L^{2}$ optimal transportation cost defined by

$\displaystyle\mathcal{W}_{2}(\mu,\nu):=\inf_{\pi\in\Pi(\mu,\nu)}\left(\mathbb{E}_{(\mathsf{X},\mathsf{Y})\sim\pi}\|\mathsf{X}-\mathsf{Y}\|_{2}^{2}\right)^{1/2},$

where $\Pi(\mu,\nu)$ denotes the set of all joint probability distributions
$\pi$ whose marginals are respectively $\mu$ and $\nu$. A distribution
$\pi\in\Pi(\mu,\nu)$ is called a coupling of $\mu$ and $\nu$.

###  2.2 Assumptions

We state the assumptions on the target distribution $\nu$ on
${\mathbb{R}}^{d}$ to support the main results. Below we use
$\mathrm{supp}(\nu)$ to denote the support of a probability distribution
$\nu$. We also use $\mathrm{diam}(\Omega)$ to represent the diameter of a set
$\Omega\subset{\mathbb{R}}^{d}$.

######  Assumption 1.

The probability distribution $\nu$ is absolutely continuous with respect to
the Lebesgue measure and has a zero mean.

######  Assumption 2.

Let $D:=(1/\sqrt{2})\mathrm{diam}(\mathrm{supp}(\nu))$. The probability
distribution $\nu$ satisfies any one of the following conditions:

  * (i)

$\nu$ is $\kappa$-semi-log-concave for some $\kappa>0$ with $D\in(0,\infty]$
and $\beta$-semi-log-convex for some $\beta>0$;

  * (ii)

$\nu$ is $\kappa$-semi-log-concave for some $\kappa\leq 0$ with
$D\in(0,\infty)$ and $\beta$-semi-log-convex for some $\beta>0$;

  * (iii)

$\nu=\gamma_{d,\sigma^{2}}*\rho$ where $\rho$ is a probability distribution
supported on a Euclidean ball of radius $R$ on $\mathbb{R}^{d}$.

######  Remark 3 (Distribution classes).

Assumption 2 covers three classes of distributions that are of great interest
in the literature on generative learning and sampling. Let us briefly remark
on the properties of the distributions:

  * (1)

Distributions with a bounded support are covered in Assumption 2-(ii). The
boundedness provides a great source of regularity, which frees us from the
requirement that the Hessian matrix of the potential function is positively
lower bounded. A negative lower bound of the Hessian matrix allows the
distribution has a multimodal landscape.

  * (2)

Strongly log-concave distributions are considered in Assumption 2-(i). For
distributions in this class, the Hessian matrix of the potential function is
both positively lower bounded and positively upper bounded.

  * (3)

Mixtures of Gaussians are considered in Assumption 2-(iii). A mixture of
Gaussians is notably a multimodal probability distribution, which is neither
strongly log-concave nor bounded.

######  Remark 4 (Score Lipschitzness).

For any $\beta>0,\kappa\in{\mathbb{R}},$ and $\kappa\leq\beta$, the
$\beta$-semi-log-convexity and $\kappa$-semi-log-concavity measure the
Lipschitzness of the smooth score function
$S(x):=\nabla_{x}\log\frac{\mathrm{d}\nu}{\mathrm{d}x}(x)$ in the sense that
$\kappa\mathrm{I}_{d}\preceq\nabla_{x}S(x)\preceq\beta\mathrm{I}_{d}$. The
Lipschitzness of the score function for a probability distribution is a common
assumption in the literature studying convergence properties of Langevin Monte
Carlo algorithms and score-based diffusion models (cf. Dalalyan (2017); Durmus
and Moulines (2017); Chen et al. (2023a); Chen, Lee and Lu (2023)).

######  Remark 5 (Sub-Gaussianity).

The probability distribution $\nu$ considered in Assumption 2 satisfies the
log-Sobolev inequality with a finite constant $C_{\mathrm{LSI}}$ depending on
$\kappa>0$ or $(\kappa,D)$ or $(\sigma,R)$ (Mikulincer and Shenfeld, 2021; Dai
et al., 2023). Let $\mathsf{V}\sim\nu$ and
$\mathsf{V}=[\mathsf{V}_{1},\mathsf{V}_{2},\cdots,\mathsf{V}_{d}]^{\top}$.
Then a standard Herbst’s argument shows that $\mathsf{V}_{i}$ is sub-Gaussian,
and its sub-Gaussian norm
$\|\mathsf{V}_{i}\|_{\psi_{2}}\asymp\sqrt{C_{\mathrm{LSI}}}$ for $1\leq i\leq
d$ owing to Ledoux (2001, Theorem 5.3). In addition, a sub-Gaussian random
variable has a finite fourth moment.

##  3 Simulation-free continuous normalizing flows

The basic idea of simulation-free CNFs is to construct an ODE-based IVP with a
tractable velocity field. The flow map of the IVP pushes forward a simple
source distribution onto the underlying target distribution. It is essential
to be able to efficiently estimate the velocity field with a random sample
from the target distribution. Since CNFs use ODEs to model a target
distribution, the corresponding velocity fields depend on the target
distribution and may have a complex, unknown structure. Therefore, it is
natural to employ nonparametric methods with deep neural networks to estimate
velocity fields.

In Table 1, we summarize the four steps of using simulation-free CNFs for
generative learning, and we discuss each step in detail below.

Table 1: Four steps to conduct generative learning via simulation-free CNFs.
Task: Generating samples from a distribution approximating the target
distribution $\nu$.
---
Step 1. Construct a simulation-free CNF defined by the IVP such that
$\nu={X_{T}}_{\\#}\mu$:
$\displaystyle\frac{\mathrm{d}X_{t}}{\mathrm{d}t}(x)=v(t,X_{t}(x)),\quad
X_{0}(x)=x\sim\mu,\quad(t,x)\in[0,T]\times{\mathbb{R}}^{d},$ where
$T>0,v:[0,T]\times{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$ is the velocity field,
and $\mu$ is a simple source distribution. Step 2. Estimate the velocity field
$v(t,x)$ with a deep neural network $\hat{v}(t,x)$. Step 3. Use a proper
numerical solver to solve the IVP associated with $\hat{v}(t,x)$, and return
the numerical solution $\hat{Y}_{T}(x)$ at time $t=T$:
$\displaystyle\frac{\mathrm{d}Y_{t}}{\mathrm{d}t}(x)=\hat{v}(t,Y_{t}(x)),\quad
Y_{0}(x)=x\sim\mu,\quad(t,x)\in[0,T]\times{\mathbb{R}}^{d}.$ Step 4. Generate
samples from ${\hat{Y}_{T}{}}_{\\#}\mu$, which approximates
$\nu={X_{T}}_{\\#}\mu$.

###  3.1 Construction of simulation-free CNFs

Based on the concept of stochastic interpolation, Liu, Gong and Liu (2023) and
Albergo and Vanden-Eijnden (2023) proposed a new class of simulation-free
CNFs. Gao, Huang and Jiao (2023) analyzed probability flows of diffusion
models and denoising diffusion implicit models in the framework of stochastic
interpolation. Let
$a_{t}:[0,1]\to{\mathbb{R}}_{+},b_{t}:[0,1]\to{\mathbb{R}}_{+}$ satisfy the
following conditions:

$\displaystyle\begin{aligned} &\dot{a}_{t}\leq 0,\quad\dot{b}_{t}\geq 0,\quad
a_{0}>0,\quad b_{0}\geq 0,\quad a_{1}=0,\quad b_{1}=1,\\\ &a_{t}>0\ \
\text{for any $t\in(0,1)$},\quad b_{t}>0\ \ \text{for any $t\in(0,1)$},\\\
&a_{t},b_{t}\in C^{2}([0,1)),\quad a_{t}^{2}\in C^{1}([0,1]),\quad b_{t}\in
C^{1}([0,1]).\end{aligned}$

Then a general class of simulation-free CNFs is constructed in Lemma 3.1 and
satisfies the requirements of Step 1 in Table 1. For simplicity, we use
$\mathrm{Law}(\mathsf{X})$ to denote the distribution of a random variable
$\mathsf{X}$.

######  Lemma 3.1 (Theorem 38 in Gao, Huang and Jiao (2023)).

Suppose that a probability distribution $\nu$ satisfies Assumptions 1 and 2.
Let $\mu=\mathrm{Law}(a_{0}\mathsf{Z}+b_{0}\mathsf{X}_{1})$ with
$\mathsf{Z}\sim\gamma_{d},\mathsf{X}_{1}\sim\nu$,
$\mathsf{X}_{t}:=a_{t}\mathsf{Z}+b_{t}\mathsf{X}_{1}$ for any $t\in(0,1).$ Let
the velocity field $v(t,x)$ be defined by

$\displaystyle v(t,x)$
$\displaystyle:=\mathbb{E}[\dot{a}_{t}\mathsf{Z}+\dot{b}_{t}\mathsf{X}_{1}|\mathsf{X}_{t}=x],\quad(t,x)\in(0,1)\times{\mathbb{R}}^{d},$
(3.1) $\displaystyle v(0,x)$ $\displaystyle:=\lim_{t\downarrow 0}v(t,x),\quad
v(1,x):=\lim_{t\uparrow 1}v(t,x),~{}~{}x\in{\mathbb{R}}^{d}.$ (3.2)

Then there exists a unique solution $(X_{t})_{t\in[0,1]}$ to the IVP

$\frac{\mathrm{d}X_{t}}{\mathrm{d}t}(x)=v(t,X_{t}(x)),\quad
X_{0}(x)=x\sim\mu,\quad(t,x)\in[0,1]\times{\mathbb{R}}^{d}.$ (3.3)

Moreover, the push-forward distribution satisfies
${X_{t}}_{\\#}\mu=\mathrm{Law}(a_{t}\mathsf{Z}+b_{t}\mathsf{X}_{1})$ with
$\mathsf{Z}\sim\gamma_{d},\mathsf{X}_{1}\sim\nu$.

Lemma 3.1 shows that the velocity field $v(t,x)$ of the simulation-free CNF
(3.3) takes the form of conditional expectations. As a result, the least
squares method, also known as the flow matching method for simulation-free
CNFs (Lipman et al., 2023), is an effective approach to estimating the
velocity field $v(t,x)$.

Among various choices of the coefficients $a_{t}$ and $b_{t}$, the linear
interpolation scenario, where

$\displaystyle a_{t}=1-t,\ b_{t}=t,$ (3.4)

has been shown to have excellent properties for generative learning tasks
(Liu, Gong and Liu, 2023; Albergo, Boffi and Vanden-Eijnden, 2023). A CNF
model with linear interpolation has been used in the implementation of a large
image generative model (Esser et al., 2024). The coefficients of the linear
interpolation are the same as those of the displacement interpolation in
optimal transport (McCann, 1997; Villani, 2009). In this paper, we focus on
the regularity and convergence properties of the CNF with linear interpolation
(3.4).

######  Corollary 3.2.

Suppose that a probability distribution $\nu$ satisfies Assumptions 1 and 2.
Let $\mathsf{X}_{0}=\mathsf{Z}\sim\gamma_{d}$, $\mathsf{X}_{1}\sim\nu$.
Consider the linear interpolation
$\mathsf{X}_{t}:=(1-t)\mathsf{Z}+t\mathsf{X}_{1}$ for any $t\in(0,1)$ and the
velocity field $v^{*}(t,x)$ defined by

$\displaystyle v^{*}(t,x)$
$\displaystyle:=\mathbb{E}[\mathsf{X}_{1}-\mathsf{Z}|\mathsf{X}_{t}=x],\quad(t,x)\in[0,1]\times{\mathbb{R}}^{d}.$
(3.5)

Then there exists a unique solution $(X_{t})_{t\in[0,1]}$ to the IVP

$\frac{\mathrm{d}X_{t}}{\mathrm{d}t}(x)=v^{*}(t,X_{t}(x)),\quad
X_{0}(x)=x\sim\gamma_{d},\quad(t,x)\in[0,1]\times{\mathbb{R}}^{d},$ (3.6)

and the push-forward distribution satisfies
${X_{t}}_{\\#}\gamma_{d}=\mathrm{Law}(\mathsf{X}_{t})$.

Corollary 3.2 implies that the push-forward distributions
$({X_{t}}_{\\#}\gamma_{d})_{t\in[0,1]}$ coincide with the marginal
distributions of the Gaussian channel
$\mathsf{X}_{t}=(1-t)\mathsf{Z}+t\mathsf{X}_{1}$. The connections between the
velocity fields and Gaussian channels have inspired moment expressions for the
derivatives of the velocity field. We show these expressions in Lemmas 3.3 and
3.4 for the purpose of examining the regularity properties of velocity fields.

######  Remark 6.

Since $x=\mathbb{E}[(1-t)\mathsf{Z}+t\mathsf{X}_{1}|\mathsf{X}_{t}=x]$ for
$(t,x)\in[0,1]\times{\mathbb{R}}^{d}$, an alternative expression of the
velocity field is given by

$v^{*}(t,x)=-\frac{1}{1-t}x+\frac{1}{1-t}\mathbb{E}[\mathsf{X}_{1}|\mathsf{X}_{t}=x],\quad(t,x)\in[0,1)\times{\mathbb{R}}^{d}.$
(3.7)

Expression (3.7) shows that the velocity field $v^{*}(t,x)$ only depends on
the conditional expectation $\mathbb{E}[\mathsf{X}_{1}|\mathsf{X}_{t}=x].$

######  Lemma 3.3 (Lemma 26 in Gao, Huang and Jiao (2023)).

The Jacobian matrix $\nabla_{x}v^{*}(t,x)$ has a covariance expression as
follows

$\nabla_{x}v^{*}(t,x)=\frac{t}{(1-t)^{3}}\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)-\frac{1}{1-t}\mathrm{I}_{d},\quad(t,x)\in[0,1)\times{\mathbb{R}}^{d}.$
(3.8)

######  Remark 7.

The covariance expression (3.8) has been used to derive regularity properties
of the velocity field $v^{*}(t,x)$ in the space variable $x$. For example, see
Proposition 29 in Gao, Huang and Jiao (2023).

######  Lemma 3.4 (Proposition 63 in Gao, Huang and Jiao (2023)).

The time derivative of the velocity field $v^{*}(t,x)$ has a moment expression
for any $t\in[0,1)$ as follows

$\displaystyle\partial_{t}v^{*}(t,x)$
$\displaystyle=-\frac{1}{(1-t)^{2}}x+\frac{1}{(1-t)^{2}}M_{1}+\frac{t+1}{(1-t)^{4}}M^{c}_{2}x-\frac{t}{(1-t)^{4}}(M_{3}-M_{2}M_{1}),$
(3.9)

where
$M_{1}:=\mathbb{E}[\mathsf{X}_{1}|\mathsf{X}_{t}=x],M_{2}:=\mathbb{E}[\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x],M^{c}_{2}:=\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)$,
and
$M_{3}:=\mathbb{E}[\mathsf{X}_{1}\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x]$
with omitted dependence on $(t,x)$.

######  Remark 8.

To quantify the regularity of the velocity field $v^{*}(t,x)$ in the time
variable $t$, one can try to bound the moments in Eq. (3.9) defined in Lemma
3.4. Following this idea, we conduct regularity analyses on the velocity field
in Appendix A and summarize the results in Theorem 5.1.

###  3.2 Flow matching

This subsection concerns Step 2 in Table 1, which estimates the velocity field
of a CNF with linear interpolation. As shown in (3.5), the velocity field
$v^{*}(t,x)=\mathbb{E}[\mathsf{X}_{1}-\mathsf{Z}|\mathsf{X}_{t}=x]$ for each
$t\in[0,1].$ For notational simplicity, let
$\mathsf{Y}:=\mathsf{X}_{1}-\mathsf{Z}$, and note that
$\mathsf{X}_{t}=(1-t)\mathsf{Z}+t\mathsf{X}_{1}$. Given $\tau\in(0,1]$, we
consider the time interval $[0,\tau]$. For $t\in[0,\tau]$, we denote
$\mathsf{X}_{t}\sim p_{t}$. When the time is a random variable distributed as
a uniform distribution on $[0,\tau]$, that is,
$\mathsf{t}\sim\mathrm{U}(0,\tau)$, we denote
$\mathsf{X}_{\mathsf{t}}|\mathsf{t}=t\sim p_{t}$. Then the flow matching
method (Lipman et al., 2023; Liu, Gong and Liu, 2023) solves a nonlinear least
squares problem for estimating
$v^{*}(t,x)=\mathbb{E}[\mathsf{Y}|\mathsf{X}_{t}=x]$ on the domain
$[0,\tau]\times{\mathbb{R}}^{d}$:

$\displaystyle
v^{*}\in\operatorname*{arg\,min}_{v}\Big{\\{}\mathcal{L}(v):=\mathbb{E}_{\mathsf{t}\sim\mathrm{U}(0,\tau)}\mathbb{E}_{\mathsf{X}_{1}\sim\nu,\mathsf{Z}\sim\gamma_{d}}\|v(\mathsf{t},\mathsf{X}_{\mathsf{t}})-\mathsf{Y}\|_{2}^{2}\Big{\\}}.$
(3.10)

In practice, $\tau$ is often taken as $1$. Here, we leave $\tau$ as a quantity
adaptive to the time regularity of the velocity field $v^{*}$. We will analyze
the time regularity of $v^{*}$ in Subsection 5.1.

Let $\\{\mathsf{Z}_{i}\\}_{i=1}^{n}$, $\\{\mathsf{X}_{1,i}\\}_{i=1}^{n}$, and
$\\{\mathsf{t}_{i}\\}_{i=1}^{n}$ be i.i.d. random samples from $\gamma_{d}$,
$\nu$, and $\mathrm{U}(0,\tau)$, respectively. For $i=1,2,\cdots,n$, we denote
$\mathsf{X}_{\mathsf{t}_{i}}:=(1-\mathsf{t}_{i})\mathsf{Z}_{i}+\mathsf{t}_{i}\mathsf{X}_{1,i}$
and $\mathsf{Y}_{i}:=\mathsf{X}_{1,i}-\mathsf{Z}_{i}$. The population risk
$\mathcal{L}(v)$ defined in (3.10) leads to the empirical risk

$\displaystyle\mathcal{L}_{n}(v):=\frac{1}{n}\sum_{i=1}^{n}\|v(\mathsf{t}_{i},\mathsf{X}_{\mathsf{t}_{i}})-\mathsf{Y}_{i}\|_{2}^{2}.$
(3.11)

We approximate the velocity field $v^{*}$ using deep neural networks. We
consider a deep ReLU network class with the input dimension $\mathtt{k}=d+1$
and the output dimension $\mathtt{d}=d$. Let
$\mathcal{F}_{n}:=\mathcal{NN}(\mathtt{S},\mathtt{W},\mathtt{D},\mathtt{B},\mathtt{L}_{x},\mathtt{L}_{t},d+1,d)$
denote a function class
$\\{f_{\theta}(t,x):{\mathbb{R}}^{d+1}\to{\mathbb{R}}^{d}\\}$ implemented by
deep ReLU networks with size $\mathtt{S}$, width $\mathtt{W}$, depth
$\mathtt{D}$, bound $\mathtt{B}$, and Lipschitz constants at most
$\mathtt{L}_{x}$ in $x$ and $\mathtt{L}_{t}$ in $t$ over
$(t,x)\in[0,\tau]\times{\mathbb{R}}^{d}$. The network parameters can depend on
the sample size $n$, and we show the fact by making $\mathcal{F}_{n}$ depend
on $n$. For any $f\in\mathcal{F}_{n}$, the Lipschitz continuity of $f$ implies
that $\|f(t,x)-f(s,x)\|_{\infty}\leq\mathtt{L}_{t}|t-s|$ and
$\|f(t,x)-f(t,y)\|_{\infty}\leq\mathtt{L}_{x}\|x-y\|_{\infty}$ for any
$s,t\in[0,\tau]$ and $x,y\in{\mathbb{R}}^{d}$. It is easy to see that
$\mathcal{F}_{n}\subseteq\mathcal{NN}(\mathtt{S},\mathtt{W},\mathtt{D},\mathtt{B},d+1,d)$,
that is a deep ReLU network class without Lipschitz regularity control. To
estimate the velocity field within the hypothesis class $\mathcal{F}_{n}$, we
consider the empirical risk minimization problem

$\displaystyle\hat{v}_{n}\in\operatorname*{arg\,min}_{v\in\mathcal{F}_{n}}\mathcal{L}_{n}(v),$
(3.12)

where $\hat{v}_{n}$ is a deep ReLU network estimator for the velocity field
$v^{*}$. We call $\hat{v}_{n}$ the flow matching estimator because it is a
minimizer of the empirical flow matching objective.

###  3.3 Forward Euler discretization

In this subsection, we proceed to Step 3 in Table 1, where a numerical solver
is used to solve the IVP:

$\displaystyle\frac{\mathrm{d}\tilde{X}_{t}}{\mathrm{d}t}(x)=\hat{v}_{n}(t,\tilde{X}_{t}(x)),~{}~{}\tilde{X}_{0}(x)\sim\gamma_{d},~{}~{}(t,x)\in[0,\tau]\times{\mathbb{R}}^{d}.$
(3.13)

The forward Euler method is a first-order numerical procedure for solving ODE-
based IVPs, which is commonly used in algorithms for sampling and generative
learning. First, we set a time grid of $[0,\tau]$ as
$t_{0}=0<t_{1}<\cdots<t_{K}=\tau\leq 1$. The forward Euler method for solving
the IVP (3.13) yields the numerical iterations:

$\displaystyle\frac{\mathrm{d}\hat{X}_{t}}{\mathrm{d}t}(x)=\hat{v}_{n}(t_{k-1},\hat{X}_{t_{k-1}}(x)),~{}~{}\hat{X}_{0}(x)\sim\gamma_{d},~{}~{}t\in[t_{k-1},t_{k}],\
k=1,2,\cdots,K.$ (3.14)

Finally, Step 4 in Table 1 is accomplished by drawing samples from
$\gamma_{d}$ and running the numerical iterations given in (3.14).

##  4 Main result: Error bounds for distribution estimation

In this section, we derive our main result, Theorem 4.4, on the error bounds
for the distribution estimator based on CNFs with linear interpolation.

###  4.1 Error decomposition

We begin with three IVPs (3.6), (3.13), and (3.14) in Section 3. The IVP (3.6)
defines the true process without approximation of the velocity field $v^{*}$
and discretization in time. The IVP (3.13) defines the neural process
resulting from replacing the velocity field $v^{*}$ with the flow matching
estimator $\hat{v}_{n}$. The IVP (3.14) is the forward Euler discretization
counterpart of the IVP (3.13).

Let $0<\underline{t}\ll 1$. As shown in Theorem 1.1 or Theorem 5.1 below, the
Lipschitz constant bound of the velocity field in the time variable is of the
order ${\mathcal{O}}(\underline{t}^{-2})$ on the time interval
$[0,1-\underline{t}]$. This order shows that the bound explodes at the time
$t=1$. To maintain the Lipschitz regularity in the time variable for flow
matching, we need to consider an early stopping time by letting the end time
$\tau=1-\underline{t}$.

Solving the IVPs (3.6), (3.13), and (3.14), we obtain the flow maps
$(X_{t})_{t\in[0,1]},(\tilde{X}_{t})_{t\in[0,1-\underline{t}]}$, and
$(\hat{X}_{t})_{t\in[0,1-\underline{t}]}$. To simplify the notations, we
denote the push-forward distributions
${X_{t}}_{\\#}\gamma_{d},{\tilde{X}_{t}{}}_{\\#}\gamma_{d},{\hat{X}_{t}{}}_{\\#}\gamma_{d}$
by $p_{t},\tilde{p}_{t},\hat{p}_{t}$, respectively. We summarize the three
processes (3.6), (3.13), and (3.14), their corresponding flow maps and push-
forward distributions defined by the three IVPs, their corresponding velocity
fields and density functions, and sources of errors in Table 2.

Table 2: A list of three IVPs and related notations defining the generative
learning process. In the first column, we present three processes defined by
the three IVPs referred in the second column. The corresponding notations are
given in the following columns. Particularly, we show the error source for
each process in the last column.

Process |  IVP |  Flow map |  Velocity field |  Push-forward distribution |  Density |  Error source
---|---|---|---|---|---|---
True process |  (3.6) |  $X_{t}(x)$ |  $v^{*}(t,X_{t}(x))$ |  ${X_{t}}_{\\#}\gamma_{d}$ |  $p_{t}$ |  Early stopping
Neural process |  (3.13) |  $\tilde{X}_{t}(x)$ |  $\hat{v}_{n}(t,\tilde{X}_{t}(x))$ |  ${\tilde{X}_{t}{}}_{\\#}\gamma_{d}$ |  $\tilde{p}_{t}$ |  Velocity estimation
Discrete process |  (3.14) |  $\hat{X}_{t}(x)$ |  $\hat{v}_{n}(t_{k-1},\hat{X}_{t_{k-1}}(x))$ |  ${\hat{X}_{t}{}}_{\\#}\gamma_{d}$ |  $\hat{p}_{t}$ |  Discretization

There are three sources of errors introduced in the generative learning
process (3.6), (3.13), and (3.14). The discretization error comes from the
forward Euler discretization. The error due to velocity estimation results
from flow matching with deep ReLU networks. The early stopping error is due to
the time singularity of the velocity field at time $t=1$. We use the
Wasserstein-2 distance ${\mathcal{W}}_{2}$ to measure the difference between
the estimated generative distribution $\hat{p}_{1-\underline{t}}$ and the
target distribution $p_{1}$. We derive an upper bound for
$\mathcal{W}_{2}(\hat{p}_{1-\underline{t}},p_{1}),$ which takes into account
all the three sources of error.

It is important to consider the trade-off between the different sources of
errors. The early stopping error is reduced when the parameter $\underline{t}$
gets smaller. However, a smaller value of $\underline{t}$ increases the time
singularity of the velocity field on the time interval $[0,1-\underline{t}]$,
thus leads to a larger error due to velocity estimation.

Keeping the error trade-off in mind, we consider a basic decomposition of the
total error in terms of the Wasserstein-2 distance as follows:

$\displaystyle\mathcal{W}_{2}(\hat{p}_{1-\underline{t}},p_{1})\leq\underbrace{\mathcal{W}_{2}(\hat{p}_{1-\underline{t}},\tilde{p}_{1-\underline{t}})}_{\text{discretization}}+\underbrace{\mathcal{W}_{2}(\tilde{p}_{1-\underline{t}},p_{1-\underline{t}})}_{\text{velocity
estimation}}+\underbrace{\mathcal{W}_{2}(p_{1-\underline{t}},p_{1})}_{\text{early
stopping}}.$ (4.1)

In (4.1), the first term
$\mathcal{W}_{2}(\hat{p}_{1-\underline{t}},\tilde{p}_{1-\underline{t}})$
measures the discretization error, the second term
$\mathcal{W}_{2}(\tilde{p}_{1-\underline{t}},p_{1-\underline{t}})$ measures
the error due to velocity estimation, and the third term
$\mathcal{W}_{2}(p_{1-\underline{t}},p_{1})$ measures the early stopping
error. We evaluate each error term in Lemmas 4.1, 4.2, and 4.3 below,
respectively.

######  Lemma 4.1 (Discretization error).

Suppose that Assumptions 1 and 2 hold. Let $\Upsilon\equiv t_{k}-t_{k-1}$ for
$k=1,2,\cdots,K$. Then the discretization error is evaluated by

$\displaystyle\mathcal{W}_{2}(\hat{p}_{1-\underline{t}},\tilde{p}_{1-\underline{t}})=\mathcal{O}\left(\sqrt{d}e^{\mathtt{L}_{x}}(\mathtt{L}_{x}\mathtt{B}+\mathtt{L}_{t})\Upsilon\right).$

Lemma 4.1 shows that the error due to the forward Euler discretization is well
controlled when the discretization step size $\Upsilon$ is sufficiently small.
We use the perturbation analysis of ODE flows to derive the error bound in
Lemma 4.1, and the proof can be found in Appendix D.

######  Remark 9.

Lemma 4.1 considers a uniform step size for the forward Euler discretization.
For more general choices of the step size, we need the general condition

$\displaystyle\Upsilon^{2}=\sum\nolimits_{k=1}^{K}(t_{k}-t_{k-1})^{3}$ (4.2)

to ensure that the error bound in Lemma 4.1 holds. This can be shown following
the proof of Lemma 4.1.

In the sequel, we frequently take expectations over
$(\mathsf{t},\mathsf{X}_{\mathsf{t}})$ whose joint distribution is specified
by $\mathsf{t}\sim\mathrm{U}(0,1-\underline{t})$ and
$\mathsf{X}_{\mathsf{t}}|\mathsf{t}=t\sim p_{t}$. For ease of notation, we may
omit the notation of the joint distribution when taking expectations over
$(\mathsf{t},\mathsf{X}_{\mathsf{t}}).$

######  Lemma 4.2 (Error due to velocity estimation).

Suppose that Assumptions 1 and 2 hold, and let $\hat{v}_{n}\in\mathcal{F}_{n}$
satisfy
$\|\hat{v}_{n}(t,x)-\hat{v}_{n}(t,y)\|_{\infty}\leq\mathtt{L}_{x}\|x-y\|_{\infty}$
for any $t\in[0,1-\underline{t}]$ and $x,y\in{\mathbb{R}}^{d}$. Then the error
due to velocity estimation is bounded by

$\displaystyle\mathcal{W}_{2}^{2}(\tilde{p}_{1-\underline{t}},p_{1-\underline{t}})\leq\exp(2\mathtt{L}_{x}+1)\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|\hat{v}_{n}(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})\|_{2}^{2}.$
(4.3)

Lemma 4.2 states that the error due to velocity estimation is controlled by
the excess risk of the flow matching estimator $\hat{v}_{n}$ when the
Lipschitz constant $\mathtt{L}_{x}$ is bounded. We will analyze the excess
risk of the flow matching estimator in Section 5.

By combining the excess risk bound for $\hat{v}_{n}$ and the $\mathcal{W}_{2}$
distance bound (4.3), we can deduce the error bound attributable to the
estimated velocity field. Generally, bounding the error due to velocity
estimation involves establishing a perturbation error bound for the ODE flow
associated with the velocity field $\hat{v}_{n}$ or $v^{*}$, as well as an
estimation error bound for the velocity field $v^{*}$. We will use the
Grönwall’s inequality to establish the perturbation error bound (4.3) based on
the $\mathcal{W}_{2}$ distance. Similar perturbation error bounds have been
obtained in Albergo and Vanden-Eijnden (2023, Proposition 3), Benton,
Deligiannidis and Doucet (2024, Theorem 1), and Gao, Huang and Jiao (2023,
Proposition 54). The proof of Lemma 4.2 is given in Appendix D.

######  Lemma 4.3 (Early stopping error).

Suppose that Assumptions 1 and 2 hold. The early stopping error is evaluated
by

$\displaystyle\mathcal{W}_{2}(p_{1-\underline{t}},p_{1})\lesssim\underline{t},$

where we omit a polynomial prefactor in $d$ and
$\mathbb{E}[\|\mathsf{X}_{1}\|_{2}^{2}]$.

The $\mathcal{W}_{2}$ distance bound in Lemma 4.3 formalizes the intuition
that the early stopping error scales with the early stopping parameter
$\underline{t}$. The proof of Lemma 4.3 uses a coupling argument, and is given
in Appendix D.

###  4.2 Error bounds for the estimated distribution

We now apply the error bounds in the preceding subsection to derive error
bounds for the distribution estimator
$\hat{\nu}_{1-\underline{t}}(\mathrm{d}x)=\hat{p}_{1-\underline{t}}(x)\mathrm{d}x$.

By Lemma 4.1, it is clear that the discretization error can be controlled by
choosing the step size $\Upsilon$ properly. Lemma 4.2 shows that the error due
to velocity estimation is upper bounded by the excess risk of the flow
matching estimator $\hat{v}_{n}$. Furthermore, we will provide a detailed
nonparametric analysis of the flow matching estimator $\hat{v}_{n}$ in Section
5.

Before presenting our main result, we first describe the trade-off between the
different sources of errors. By Theorem 5.12, the excess risk of the flow
matching estimator $\hat{v}_{n}$ satisfies

$\displaystyle\mathbb{E}_{\mathbb{D}_{n}}\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|\hat{v}_{n}(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})\|_{2}^{2}\lesssim(n\underline{t}^{2})^{-2/(d+3)}\mathrm{polylog}(n)\log(1/\underline{t}),$

where $\mathrm{polylog}(n)$ stands for a polylogarithmic prefactor in $n$.
Consequently, the error due to velocity estimation satisfies

$\displaystyle\mathbb{E}_{\mathbb{D}_{n}}\mathcal{W}_{2}(\tilde{p}_{1-\underline{t}},p_{1-\underline{t}})\lesssim(n\underline{t}^{2})^{-1/(d+3)}\mathrm{polylog}(n)\log(1/\underline{t}),$

where we use Lemma 4.2 and the bound $\mathtt{L}_{x}\asymp\log(\log n)$.
According to Lemma 4.3, the early stopping error is upper bounded by
$\mathcal{W}_{2}(p_{1-\underline{t}},p_{1})\lesssim\underline{t}$. By
substituting the error bounds into the error decomposition (4.1), it follows
that

$\displaystyle\mathbb{E}_{\mathbb{D}_{n}}\mathcal{W}_{2}(\hat{p}_{1-\underline{t}},p_{1})$
$\displaystyle\leq\mathcal{W}_{2}(\hat{p}_{1-\underline{t}},\tilde{p}_{1-\underline{t}})+\mathbb{E}_{\mathbb{D}_{n}}\mathcal{W}_{2}(\tilde{p}_{1-\underline{t}},p_{1-\underline{t}})+\mathcal{W}_{2}(p_{1-\underline{t}},p_{1})$
$\displaystyle\lesssim\Big{\\{}\underbrace{e^{\mathtt{L}_{x}}(\mathtt{L}_{x}\mathtt{B}+\mathtt{L}_{t})\Upsilon}_{\text{controlled
by
}\Upsilon}+\underbrace{(n\underline{t}^{2})^{-1/(d+3)}+\underline{t}}_{\text{trade-
off on }\underline{t}}\Big{\\}}\mathrm{polylog}(n)\log(1/\underline{t})$
$\displaystyle\lesssim\Big{\\{}\underline{t}^{-2}\Upsilon+(n\underline{t}^{2})^{-1/(d+3)}+\underline{t}\Big{\\}}\mathrm{polylog}(n)\log(1/\underline{t}),$
(4.4)

where we use the following bounds in deriving (4.4)

$\displaystyle\mathtt{L}_{x}\asymp\log(\log n),\quad
e^{\mathtt{L}_{x}}\asymp\log n,\quad\mathtt{B}\asymp\log(\log
n),\quad\mathtt{L}_{t}\asymp\log(\log n)\underline{t}^{-2}.$

Our main result stated below is obtained by balancing the error terms on the
right-hand side of (4.4).

######  Theorem 4.4 (Distribution estimation error).

Suppose that Assumptions 1 and 2 hold. Let us set $A\asymp\log(\log n)$,
$NL\asymp n^{d/(2d+10)}$, $\underline{t}\asymp n^{-1/(d+5)}$, and
$\Upsilon\lesssim n^{-3/(d+5)}$. We consider the deep ReLU network class
${\mathcal{F}}_{n}=\mathcal{NN}(\mathtt{S},\mathtt{W},\mathtt{D},\mathtt{B},\mathtt{L}_{x},\mathtt{L}_{t},d+1,d)$
whose parameters satisfy the following bounds

$\displaystyle\mathtt{S}\asymp\underline{t}^{-2}(NL)^{2/d}(N\log N)^{2}L\log
L,\quad\mathtt{W}\asymp\underline{t}^{-2}(NL)^{2/d}N\log N,$
$\displaystyle\mathtt{D}\asymp L\log L,\quad\mathtt{B}\asymp
A,\quad\mathtt{L}_{x}\asymp A,\quad\mathtt{L}_{t}\asymp A\underline{t}^{-2}.$

For any random sample
$\mathbb{D}_{n}:=\\{(\mathsf{Z}_{i},\mathsf{X}_{1,i},\mathsf{t}_{i})\\}_{i=1}^{n}$
satisfying $n\geq\mathrm{Pdim}(\mathcal{F}_{n})$, the distribution estimation
error of the CNF learned with linear interpolation and flow matching is upper
bounded as

$\displaystyle\mathbb{E}_{\mathbb{D}_{n}}\mathcal{W}_{2}(\hat{p}_{1-\underline{t}},p_{1})=\widetilde{\mathcal{O}}(n^{-\frac{1}{d+5}}),$

where we omit a polylogarithmic prefactor in $n$.

The proof of Theorem 4.4 is given in Appendix D.

######  Remark 10.

As shown in (4.4), we need to consider a trade-off on the early stopping
parameter $\underline{t}$ as well as an appropriate order of the step size
$\Upsilon$. By setting $\underline{t}\asymp n^{-1/(d+5)}$ and
$\Upsilon\lesssim n^{-3/(d+5)}$, we attain a concrete convergence rate
$\widetilde{\mathcal{O}}(n^{-1/(d+5)})$ of the distribution estimator
$\hat{p}_{1-\underline{t}}.$

######  Remark 11.

We consider the uniform step size $\Upsilon$ in deriving the distribution
estimation error bound in Theorem 4.4. The uniform step size is common in
implementing numerical solvers for ODE flows. Additionally, the condition
(4.2) in Remark 9 provides a guideline on general settings of the
discretization step size.

##  5 Error analysis of flow matching

In this section, we first present some useful regularity properties of the
velocity field $v^{*}$, which are essential to the convergence analysis of the
flow matching estimator $\hat{v}_{n}$ defined in (3.12). The error from the
flow matching estimation constitutes a main source of the total error for the
distribution estimation given in Theorem 4.4.

###  5.1 Regularity of velocity fields

The regularity properties of the velocity field are needed in studying the
nonparametric estimation error of flow matching. We summarize the regularity
results of the velocity field in Theorem 5.1.

######  Theorem 5.1.

Suppose that Assumptions 1 and 2 are satisfied, and let $0<\underline{t}\ll
1$. Then the velocity field
$v^{*}(t,x):[0,1]\times{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}^{d}$ has the
following regularity properties:

  1. (1)

For any $s,t\in[0,1-\underline{t}]$ and $x\in{\mathbb{R}}^{d}$,
$\|v^{*}(t,x)-v^{*}(s,x)\|_{\infty}\leq L_{t}|t-s|$ with
$L_{t}\lesssim\underline{t}^{-2}$;

  2. (2)

For any $x,y\in{\mathbb{R}}^{d}$ and $t\in[0,1]$,
$\|v^{*}(t,x)-v^{*}(t,y)\|_{\infty}\leq L_{x}\|x-y\|_{\infty}$ with
$L_{x}\lesssim 1$;

  3. (3)

$\sup_{(t,x)\in[0,1]\times\Omega_{A}}\|v^{*}(t,x)\|_{\infty}\leq B$ with
$B\lesssim A$,

where we omit constants in $d,\kappa,\beta,\sigma,D,R$ and denote
$\Omega_{A}:=[-A,A]^{d}$.

Theorem 5.1 states that the Lipschitz regularity of the velocity field $v^{*}$
holds in the time variable $t$ and the space variable $x$. Moreover, the
Lipschiz constant in $x$ is uniformly bounded for any
$(t,x)\in[0,1]\times{\mathbb{R}}^{d}$, and the Lipschitz constant in $t$ is
bounded for any $(t,x)\in[0,1-\underline{t}]\times{\mathbb{R}}^{d}$ but
depends on $\underline{t}$. Due to the uniform Lipschitzness in $x$, the
velocity field $v^{*}$ further satisfies the linear growth property (3).

######  Remark 12.

We note that the velocity field may be singular at time $t=1$, since the
Lipschitz constant bound of $L_{t}$ explodes at time $t=1$. We quantify the
time singularity through the upper bound $L_{t}\lesssim\underline{t}^{-2}$
where $0<\underline{t}\ll 1$. Taking the time singularity into account, we set
the end time $\tau=1-\underline{t}$ in Subsection 4.1.

######  Remark 13.

The global Lipschitz continuity of the velocity field in $x$ ensures that the
associated IVP has a unique solution, according to the Cauchy-Lipschitz
theorem (Hartman, 2002).

Proof idea of Theorem 5.1. The proof of Theorem 5.1 can be found in Appendix
A.5. As shown in Lemmas 3.3 and 3.4, the derivatives of the velocity field
$v^{*}$ can be expressed in terms of the moments of
$\mathsf{X}_{1}|\mathsf{X}_{t}$. The key idea of the proof is to bound these
moments. Under Assumptions 1 and 2, Gao, Huang and Jiao (2023) have shown that
the covariance matrix $\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t})$ is both
lower and upper bounded uniformly in $t\in[0,1]$ (cf. Lemma A.1 in Appendix
A). As a result, we can prove that the Lipschitz property (1) holds. The
linear growth property (3) follows from the Lipschitz property (1). To prove
the Lipschitz property (2), we derive bounds for the moments
$M_{1}=\mathbb{E}[\mathsf{X}_{1}|\mathsf{X}_{t}],M_{2}=\mathbb{E}[\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}]$,
and
$M_{3}=\mathbb{E}[\mathsf{X}_{1}\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}]$
by transforming the bounds of the covariance matrix
$M^{c}_{2}=\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t})$. In particular, the
bound transformations can be derived based on the Hatsell-Nolte identity
(Dytso, Poor and Shamai Shitz, 2023, Proposition 1), the Brascamp-Lieb
inequality (Brascamp and Lieb, 1976), and a basic inequality on
$M_{3}-M_{2}M_{1}$. Moreover, we validate the sharpness of moment bounds using
a Gaussian example.

###  5.2 Error decomposition of flow matching

The starting point of our analysis is the decomposition of the excess risk of
$\hat{v}_{n}$ below.

######  Lemma 5.2.

Let $\mathsf{t}\sim\mathrm{U}(0,1-\underline{t})$. For any random sample
$\mathbb{D}_{n}:=\\{(\mathsf{Z}_{i},\mathsf{X}_{1,i},\mathsf{t}_{i})\\}_{i=1}^{n}$,
the excess risk of the flow matching estimator $\hat{v}_{n}$ satisfies

$\displaystyle\mathbb{E}_{\mathbb{D}_{n}}\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|\hat{v}_{n}(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})\|_{2}^{2}$
$\displaystyle=\mathbb{E}_{\mathbb{D}_{n}}[\mathcal{L}(\hat{v}_{n})-\mathcal{L}(v^{*})]$
$\displaystyle\leq\mathcal{E}_{\mathrm{stoc}}+2\mathcal{E}_{\mathrm{appr}},$
(5.1)

where the stochastic error
$\mathcal{E}_{\mathrm{stoc}}:=\mathbb{E}_{\mathbb{D}_{n}}[\mathcal{L}(v^{*})-2\mathcal{L}_{n}(\hat{v}_{n})+\mathcal{L}(\hat{v}_{n})]$
and the approximation error
$\mathcal{E}_{\mathrm{appr}}:=\inf_{v\in\mathcal{F}_{n}}\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|v(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})\|_{2}^{2}$.

The proof of Lemma 5.2 is given in Appendix C.1. The decomposition (5.1) of
the excess risk can be considered a bias-variance decomposition. The
stochastic error $\mathcal{E}_{\mathrm{stoc}}$ bounds the variance term of the
flow matching estimator, and the approximation error
$\mathcal{E}_{\mathrm{appr}}$ represents the bias term of the flow matching
estimator. We derive bounds for $\mathcal{E}_{\mathrm{stoc}}$ and
$\mathcal{E}_{\mathrm{appr}}$. Then the best bound for the excess risk under
the decomposition (5.1) is obtained by balancing these two error bounds.

###  5.3 Approximation error

We derive the error bounds for approximating Lipschitz functions using deep
ReLU networks with Lipschitz regularity. The results are presented in Theorem
5.4, which is crucial for bounding the approximation error of the flow
matching estimator $\hat{v}_{n}$. To address the challenge posed by the
unbounded support of the velocity field $v^{*}$ in the space variable $x$, we
use the standard technique of truncated approximation. This allows us to
divide the approximation error into the truncated approximation error and the
truncation error.

######  Lemma 5.3.

For $\bar{v}\in\mathcal{F}_{n}$ and any $A>0$, the approximation error
satisfies a basic inequality as follows:

$\displaystyle\mathcal{E}_{\mathrm{appr}}=\inf_{v\in\mathcal{F}_{n}}\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|v(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})\|_{2}^{2}\lesssim\mathcal{E}_{\mathrm{appr}}^{\mathrm{trunc}}+\mathcal{E}_{\mathrm{trunc}},$
(5.2)

where the truncated approximation error

$\displaystyle\mathcal{E}_{\mathrm{appr}}^{\mathrm{trunc}}:=\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|[\bar{v}(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})]\mathds{1}_{\Omega_{A}}(\mathsf{X}_{\mathsf{t}})\|_{2}^{2},$

and the truncation error

$\displaystyle\mathcal{E}_{\mathrm{trunc}}:=\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|[\bar{v}(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})]\mathds{1}_{\Omega_{A}^{c}}(\mathsf{X}_{\mathsf{t}})\|_{2}^{2}.$

Lemma 5.3 follows from the triangle inequality and the inequality $2ab\leq
a^{2}+b^{2}$ for any $a,b\in{\mathbb{R}}$. We bound the truncated
approximation error $\mathcal{E}_{\mathrm{appr}}^{\mathrm{trunc}}$ by
considering deep ReLU network approximation of the velocity field on the
$(d+1)$-dimensional hypercube. The truncation error
$\mathcal{E}_{\mathrm{trunc}}$ measures how fast the approximation error
decays according to the tail property of the probability distribution $p_{t}$
with $t\in[0,1-\underline{t}]$.

Approximation with Lipschitz regularity control. We study the capacity of an
approximation function $\bar{v}(t,x)$ implemented by a deep ReLU network with
Lipschitz regularity for approximating the velocity field $v^{*}(t,x).$ For
balancing the approximation error with the stochastic and discretization
errors to obtain an overall error bound for the distribution estimation error,
we construct the approximation function $\bar{v}(t,x)$ so that it satisfies
the following three requirements:

  * (a)

Good approximation power under the sup norm over the hypercube
$\Omega_{\underline{t},A}:=[0,1-\underline{t}]\times[-A,A]^{d}$,

  * (b)

Lipschitz continuity with respect to both the time variable $t$ and the space
variable $x$,

  * (c)

Independent regularity in the time variable $t$ and the space variable $x$.

Let us briefly comment on each of these requirements. Requirement (a) is
needed for bounding the approximation error of the flow matching estimator
$\hat{v}_{n}.$ The time-space Lipschitz regularity of the approximation
function $\bar{v}(t,x)$ required in (b) is essential to bounding the
discretization error and the error due to velocity estimation. Requirement (c)
stems from the time singularity of the velocity field at $t=1$ and the
different roles of the time regularity and the space regularity in the error
analysis.

######  Theorem 5.4.

For any $N,L\in{\mathbb{N}}$, there exists a function $\bar{v}(t,x)$
implemented by a deep ReLU network with width
$\mathcal{O}(\underline{t}^{-2}(NL)^{2/d}N\log N)$, depth $\mathcal{O}(L\log
L)$, and size $\mathcal{O}(\underline{t}^{-2}(NL)^{2/d}(N\log N)^{2}L\log L)$
such that the following properties hold simultaneously:

  1. (i)

Boundedness and Lipschitz regularity: for any $s,t\in[0,1-\underline{t}]$ and
any $x,y\in{\mathbb{R}}^{d}$,

$\displaystyle\sup_{(t,x)\in[0,1-\underline{t}]\times{\mathbb{R}}^{d}}\|\bar{v}(t,x)\|_{\infty}$
$\displaystyle\lesssim A,$
$\displaystyle\sup_{x\in{\mathbb{R}}^{d}}\|\bar{v}(t,x)-\bar{v}(s,x)\|_{\infty}$
$\displaystyle\lesssim A\underline{t}^{-2}|t-s|,$
$\displaystyle\sup_{t\in[0,1-\underline{t}]}\|\bar{v}(t,x)-\bar{v}(t,y)\|_{\infty}$
$\displaystyle\lesssim A\|x-y\|_{\infty}.$

  2. (ii)

Approximation error bound:

$\displaystyle\sup_{(t,x)\in\Omega_{\underline{t},A}}\|\bar{v}(t,x)-v^{*}(t,x)\|_{\infty}$
$\displaystyle\lesssim A(NL)^{-2/d}.$

Note that we omit some prefactors in $d,\kappa,\beta,\sigma,D,R$ and denote
$\Omega_{\underline{t},A}:=[0,1-\underline{t}]\times[-A,A]^{d}$.

The proof of Theorem 5.4 is given in Appendix B.3. Let $l\in\mathbb{N}$ and
$\Omega\subset{\mathbb{R}}^{l}$ denote a subset of $\mathbb{R}^{d}$. We denote
by $L^{p}(\Omega)$ the standard Lebesgue space on $\Omega$ with the Lebesgue
norm $\|\cdot\|_{L^{p}(\Omega)}$ for $p\in[1,\infty]$. Let $k\in\mathbb{N}$.
We show the definitions of the Sobolev space $W^{k,\infty}(\Omega)$, the
Sobolev norm $\|\cdot\|_{W^{k,\infty}(\Omega)}$, and the Sobolev semi-norm
$|\cdot|_{W^{k,\infty}(\Omega)}$ in Appendix E.

Proof idea of Theorem 5.4. To fulfill Requirement (c), we use different
approaches to approximation in the time variable and the space variable. Our
approximation approaches can ensure the constructed approximation function to
be global Lipschitz for fulfilling Requirement (b). We take four steps to
construct the time-space approximation function $\bar{v}$ with Lipschitz
regularity control.

The first step is to derive an $L^{\infty}([0,1]^{d})$ error bound of using
deep ReLU networks for approximating a Lipschitz function in the space
variable. We show that the constructed deep neural approximation function is
globally Lipschitz in the space variable. The approximation results are
presented in Appendix B.1, and we summarize them here. Lemma B.1 and Corollary
B.2 show that for any $N,L\in\mathbb{N}$ and any $f\in
W^{1,\infty}((0,1)^{d})$, there exists a function $\phi$ implemented by a deep
ReLU network with width $\mathcal{O}(2^{d}dN\log N)$ and depth
$\mathcal{O}(d^{2}L\log L)$ such that
$|\phi|_{W^{1,\infty}((0,1)^{d})}\lesssim|f|_{W^{1,\infty}((0,1)^{d})}$ and
that $\|\phi-f\|_{L^{\infty}([0,1]^{d})}\lesssim(NL)^{-2/d},$ omitting the
prefactors depending only on $d$.

The second step is to derive an $L^{\infty}([0,1])$ approximation error bound
of deep ReLU networks for approximating a Lipschitz function in the time
variable. We establish an $L^{\infty}([0,1])$ approximation bound with
Lipschitz regularity control in Appendix B.2. Lemma B.8 states the main
results for approximation in time in such a way: for any $M\in\mathbb{N}$ and
any $f\in W^{1,\infty}((0,1))$, there exists a function $\xi$ implemented by a
deep ReLU network with width $\mathcal{O}(M)$, depth $\mathcal{O}(1)$, and
size ${\mathcal{O}}(M)$ such that
$|\xi|_{W^{1,\infty}((0,1))}\lesssim|f|_{W^{1,\infty}((0,1))}$ and that
$\|\xi-f\|_{L^{\infty}([0,1])}\lesssim|f|_{W^{1,\infty}((0,1))}/M$.

In the third step, we combine the constructed approximation in Lemmas B.1 and
B.8 to establish an $L^{\infty}(\Omega_{\underline{t},A})$ approximation bound
for the time-space approximation of the velocity field $v^{*}$. This
guarantees that Requirement (a) is fulfilled.

The last step is to show that the constructed time-space approximation
satisfies the remaining Requirements (b) and (c). We summarize these
discussions in Theorem 5.4 and present detailed construction and derivations
in the proof of Theorem 5.4.

######  Remark 14 (Optimality).

Under the assumption of continuous parameter selection, DeVore, Howard and
Micchelli (1989, Theorem 4.2) and Yarotsky (2017, Theorem 3) provided a lower
bound $\Omega(\epsilon^{-d/k})$ on the number of parameters for parametric
approximations in the Sobolev space $W^{k,\infty}([0,1]^{d})$, using the
approach of continuous nonlinear widths, when the $L^{\infty}$ approximation
error is no more than $\epsilon$. Our approximation rate for the time variable
in the Sobolev space $W^{1,\infty}([0,1])$ matches this lower bound in the
sense that a deep ReLU network with size ${\mathcal{O}}(S)$ can yield an
$L^{\infty}$ approximation error no more than ${\mathcal{O}}(1/S)$. Suppose
that the deep ReLU network has width ${\mathcal{O}}(N)$, depth
${\mathcal{O}}(L)$, and size ${\mathcal{O}}(S)$ with $S\asymp N^{2}L$. The
approximation rate ${\mathcal{O}}(S^{-k/d})$ in $W^{k,\infty}([0,1]^{d})$ can
be improved to the nearly optimal rate
${\mathcal{O}}((NL)^{-2k/d}\mathrm{polylog}(NL))$ with the bit-extraction
technique (Bartlett, Maiorov and Meir, 1998; Bartlett et al., 2019; Lu et al.,
2021). Our approximation rate for the space variable in the Sobolev space
$W^{1,\infty}([0,1]^{d})$ is nearly optimal in the sense that a deep ReLU
network with width ${\mathcal{O}}(N)$ and depth ${\mathcal{O}}(L)$ can yield
an $L^{\infty}$ approximation error no more than
${\mathcal{O}}((NL)^{-2/d}(\log(NL))^{2/d})$.

The $L^{\infty}(\Omega_{\underline{t},A})$ approximation error bound of
Theorem 5.4 implies the following $L^{2}$ bound of the truncated approximation
error for analyzing the flow matching estimator $\hat{v}_{n}$.

######  Corollary 5.5.

The truncated approximation error satisfies

$\displaystyle\mathcal{E}_{\mathrm{appr}}^{\mathrm{trunc}}=\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|[\bar{v}(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})]\mathds{1}_{\Omega_{A}}(\mathsf{X}_{\mathsf{t}})\|_{2}^{2}\lesssim
A^{2}(NL)^{-4/d},$

where we omit a constant in $d,\kappa,\beta,\sigma,D,R$.

As elaborated in the proof of Theorem 5.4 given in Appendix B, the deep ReLU
network implementing $\bar{v}$ consists of
${\mathcal{O}}(\underline{t}^{-2}(NL)^{2/d})$ parallel subnetworks which have
width ${\mathcal{O}}(N\log N)$ and depth ${\mathcal{O}}(L\log L)$. We take
advantage of the parallel structure and the construction of each subnetwork to
estimate the complexity of the deep ReLU network class ${\mathcal{F}}_{n}$
implementing $\bar{v}$. We derive the complexity of the deep ReLU network
class ${\mathcal{F}}_{n}$ in Lemma 5.6.

######  Lemma 5.6.

Suppose that Assumptions 1 and 2 hold. The complexity of the deep ReLU network
class ${\mathcal{F}}_{n}$ implementing $\bar{v}$ is quantified by

$\displaystyle\mathtt{S}\asymp\underline{t}^{-2}(NL)^{2/d}(N\log N)^{2}L\log
L,\quad\mathtt{W}\asymp\underline{t}^{-2}(NL)^{2/d}N\log N,$
$\displaystyle\mathtt{D}\asymp L\log L,\quad\mathtt{B}\asymp
A,\quad\mathtt{L}_{x}\asymp A,\quad\mathtt{L}_{t}\asymp A\underline{t}^{-2},$

where we omit some prefactors in $d,\kappa,\beta,\sigma,D,R$.

Lemma 5.6 follows from the bounds for the number of parameters in the deep
ReLU network implementing $\bar{v}$ in Theorem 5.4.

Through Theorem 5.4 and Corollary 5.5, we have established bounds for the
truncated approximation error $\mathcal{E}_{\mathrm{appr}}^{\mathrm{trunc}}$.
In what follows, we focus on the truncation error
$\mathcal{E}_{\mathrm{trunc}}$. We show the sub-Gaussian property of
$\mathsf{X}_{t}\sim p_{t}$ in Lemma 5.7 under Assumptions 1 and 2. In Lemma
5.8, we prove that the truncation error $\mathcal{E}_{\mathrm{trunc}}$ decays
very fast in the parameter $A$, as a result of the sub-Gaussian property of
$p_{t}$.

######  Lemma 5.7 (Tail probability).

Let $\mathsf{X}_{t}=(1-t)\mathsf{Z}+t\mathsf{X}_{1}$ with
$\mathsf{Z}\sim\gamma_{d}$, $\mathsf{X}_{1}\sim\nu$, and $t\in[0,1]$. Suppose
that Assumptions 1 and 2 are satisfied. For any $A>0$, it holds that

$\displaystyle\sup_{t\in[0,1]}\mathbb{P}(\mathsf{X}_{t}\in\Omega_{A}^{c})\leq
2d\exp\left(-\frac{C_{2}A^{2}}{C_{\mathrm{LSI}}}\right),$ (5.3)

where $C_{2}$ is a universal constant, and $C_{\mathrm{LSI}}>0$ depends on
$\kappa,\beta,\sigma,D$, and $R$.

######  Lemma 5.8 (Truncation error).

Suppose that Assumptions 1 and 2 are satisfied. For any $A>0$, the truncation
error satisfies

$\displaystyle\mathcal{E}_{\mathrm{trunc}}=\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|[\bar{v}(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})]\mathds{1}_{\Omega_{A}^{c}}(\mathsf{X}_{\mathsf{t}})\|_{2}^{2}\lesssim
A^{2}\exp(-C_{3}A^{2}/C_{\mathrm{LSI}}),$

where $C_{3}$ is a universal constant, and we omit a constant in
$d,\kappa,\beta,\sigma,D,R$, and the fourth moment of the target
$\mathsf{X_{1}}$.

The proofs of Lemmas 5.7 and 5.8 are given in Appendix C.2. We are now ready
to provide an upper bound for the approximation error
$\mathcal{E}_{\mathrm{appr}}.$

######  Corollary 5.9.

Suppose that Assumptions 1 and 2 hold. For any $N,L\in{\mathbb{N}}$ and $A>0$,
the approximation error is evaluated by

$\displaystyle\mathcal{E}_{\mathrm{appr}}\lesssim
A^{2}(NL)^{-4/d}+A^{2}\exp(-C_{3}A^{2}/C_{\mathrm{LSI}}).$

Corollary 5.9 holds by combining (5.2) in Lemma 5.3, Lemma 5.8, and Corollary
5.5.

###  5.4 Stochastic error

We now establish upper bounds for the stochastic error of the estimated
velocity field based on a class of deep ReLU networks.

######  Lemma 5.10.

Consider the flow matching model and the hypothesis class
$\mathcal{F}_{n}\subseteq\mathcal{NN}(\mathtt{S},\mathtt{W},\mathtt{D},\mathtt{B},d+1,d)$.
For any $n\in\mathbb{N}$ satisfying $n\geq\mathrm{Pdim}(\mathcal{F}_{n})$, the
stochastic error satisfies

$\displaystyle\mathcal{E}_{\mathrm{stoc}}=\mathbb{E}_{\mathbb{D}_{n}}[\mathcal{L}(v^{*})-2\mathcal{L}_{n}(\hat{v}_{n})+\mathcal{L}(\hat{v}_{n})]\lesssim\frac{1}{n}(\log
n)^{4}dA^{4}\mathtt{S}\mathtt{D}\log(\mathtt{S})\log(An^{2}).$

The proof of Lemma 5.10 is given in Appendix C.3.

######  Corollary 5.11.

Suppose that Assumptions 1 and 2 hold. The stochastic error satisfies

$\displaystyle\mathcal{E}_{\mathrm{stoc}}\lesssim\frac{1}{n}\underline{t}^{-2}(NL)^{2+2/d}(\log
N\log L)^{2}A^{4}\log(A)\log(\underline{t}^{-2}(NL)^{2/d}(N\log N)^{2}L\log
L),$

where we omit a polylogarithmic prefactor in $n$ and a prefactor in
$d,\kappa,\beta,\sigma,D,R$.

Corollary 5.11 follows from Lemmas 5.6 and 5.10.

###  5.5 Overall error bound for the estimated velocity field

By Lemma 5.2, the overall error for the estimated velocity field is bounded by
the sum of the approximation error $\mathcal{E}_{\mathrm{appr}}$ and the
stochastic error $\mathcal{E}_{\mathrm{stoc}}$. The approximation error is
analyzed in Subsection 5.3 and an upper bound is given in Corollary 5.9. It
decreases at a fast rate as the depth and width of the deep ReLU networks
grow. The stochastic error $\mathcal{E}_{\mathrm{stoc}}$ is analyzed in
Subsection 5.4 and its upper bound is provided in Corollary 5.11. The
stochastic error increases when the size and depth of the deep ReLU networks
grow, as a result of the increasing complexity of the hypothesis class
${\mathcal{F}}_{n}$. By balancing the bounds for $\mathcal{E}_{\mathrm{appr}}$
and $\mathcal{E}_{\mathrm{stoc}}$, we obtain the best error bound for the flow
matching estimator $\hat{v}_{n}$ under the error decomposition in Lemma 5.2.

######  Theorem 5.12 (Flow matching error).

Suppose that Assumptions 1 and 2 are satisfied. Let
$NL\asymp(n\underline{t}^{2})^{d/(2d+6)}$ and $A\asymp\log(\log n)$. Then the
excess risk of flow matching satisfies

$\displaystyle\mathbb{E}_{\mathbb{D}_{n}}\mathbb{E}_{(\mathsf{t},\mathsf{X}_{\mathsf{t}})}\|\hat{v}_{n}(\mathsf{t},\mathsf{X}_{\mathsf{t}})-v^{*}(\mathsf{t},\mathsf{X}_{\mathsf{t}})\|_{2}^{2}\lesssim(n\underline{t}^{2})^{-2/(d+3)},$

where we omit a polylogarithmic prefactor in $n$, a prefactor in
$\log(1/\underline{t})$, and a prefactor in $d,\kappa,\beta,\sigma,D,$ and
$R.$

The proof of Theorem 5.12 is given in Appendix C.4.

######  Remark 15.

In Theorem 5.12, the polynomial prefactor in $1/\underline{t}$ is due to the
singularity of $v^{*}$ in the time variable $t$. Without the singularity at
time $t=1$, the convergence rate of the flow matching error in Theorem 5.12
becomes $n^{-2/(d+3)}\mathrm{polylog}(n)$, which is nearly minimax optimal for
nonparametric least squares regression in the Sobolev space
$W^{1,\infty}([0,1]^{d+1})$ according to Stone (1982).

##  6 Related work

Process-based generative models aim to construct a stochastic process that
transports an easy-to-sample source probability distribution to the target
distribution. This goal is achieved by estimating a nonlinear transport map
implemented through deep neural networks based on a random sample from the
target distribution. CNFs and diffusion models are two prominent approaches
that have been developed for deep generative learning. Many researchers have
considered the theoretical properties of various generative learning methods.
In this section, we discuss the connections and differences between our work
and the existing studies. We focus on the studies concerning CNFs and
diffusion models that are most relevant to our work. We also discuss the
differences between the neural network approximation theory developed in this
work and those in the existing literature, focusing on the regularity
properties of the neural network functions. In particular, we highlight the
fact that our approximation results concern velocity field functions that have
different regularities in the space and time variables, while the existing
results are only applicable to functions with the same regularity in all the
variables.

###  6.1 Continuous normalizing flows

CNFs are an ODE-based generative learning approach which estimates a
stochastic process for sampling from the target distribution. Marzouk et al.
(2023) conducted a nonparametric statistical convergence analysis for
simulation-based CNF distribution estimators trained through likelihood
maximization. However, this analysis does not extend to simulation-free CNFs.
Probability flow ODEs (Song et al., 2021), denoising diffusion implicit models
(DDIMs) (Song, Meng and Ermon, 2021), and flow matching methods (Liu, Gong and
Liu, 2023; Albergo and Vanden-Eijnden, 2023; Lipman et al., 2023) all fall
under the category of simulation-free CNFs. In these models, either the score
function or the velocity field is estimated. The overall error analysis needs
to address both the estimation error of the velocity field (or score function)
and the discretization error.

In existing literature, it is typical to assume strong regularity conditions
directly on the velocity field (or score function) and its estimator.
Furthermore, current studies often only consider certain sources of errors,
neglecting either the discretization error or the estimation error of the
velocity field (or score function). In contrast, our results are derived based
on assumptions about the target distribution. Additionally, our analysis
encompasses the error due to velocity estimation, the discretization error of
the forward Euler solver, and the early stopping error. These errors are
included in the overall error bound. We provide a summary of the comparison
between our work and relevant existing studies in Table 3, with more detailed
commentary provided below.

Table 3: Comparison of convergence analyses of simulation-free CNFs. We use
$W_{2}$, KL, TV to represent the Wasserstein-2 distance, the Kullback-Leibler
divergence, and the total variation distance. We say a numerical sampler is
“mixed” if it is a combination of deterministic and stochastic samplers. For
assumptions on velocity fields or velocity field estimators, we mark “Yes” if
the assumptions are required and “No” if not. Since the unknown nonlinear part
of the velocity field is a score function, assumptions and estimation error
bounds on score functions or assumptions on score estimators can be regarded
as those on velocity fields or velocity field estimators.

|  Distribution metric |  Numerical sampler |  Estimation error bound of velocity fields |  Perturbation error bound |  Discretization error bound |  Assumptions on velocity fields |  Assumptions on estimated velocity fields
---|---|---|---|---|---|---|---
Albergo and Vanden-Eijnden (2023) |  $W_{2}$ |  Deterministic |  ✗ |  ✓ |  ✗ |  No |  Yes
Chen, Daras and Dimakis (2023) |  KL |  Deterministic |  ✗ |  ✗ |  ✓ |  Yes |  Yes
Albergo, Boffi and Vanden-Eijnden (2023) |  KL |  Deterministic |  ✗ |  ✓ |  ✗ |  No |  No
Chen et al. (2023b) |  TV |  Mixed |  ✗ |  ✓ |  ✓ |  Yes |  Yes
Benton, Deligiannidis and Doucet (2024) |  $W_{2}$ |  Deterministic |  ✗ |  ✓ |  ✗ |  Yes |  Yes
Li et al. (2024) |  TV |  Deterministic |  ✗ |  ✓ |  ✓ |  No |  Yes
Gao and Zhu (2024) |  $W_{2}$ |  Deterministic |  ✗ |  ✓ |  ✓ |  Yes |  Yes
This work |  $W_{2}$ |  Deterministic |  ✓ |  ✓ |  ✓ |  No |  No

Albergo and Vanden-Eijnden (2023) derived a perturbation error bound similar
to that in Lemma 4.3 for the CNF distribution estimator, under a Lipschitz
assumption for the estimated velocity field. Chen, Daras and Dimakis (2023)
conducted convergence analyses of DDIM-type samplers with the Kullback-Leibler
(KL) divergence, assuming second-order smoothness in the space variable and
Hölder-type regularity in the time variable for the score function, while
ignoring the score estimation error. Albergo, Boffi and Vanden-Eijnden (2023)
also derived a new perturbation error bound on the CNF distribution estimator
using the KL divergence.

Chen et al. (2023b) provided polynomial-time convergence guarantees for
distribution estimation using the probability flow ODE trained with denoising
score matching and simulated with additional randomness. To derive these
convergence rates, Chen et al. (2023b) assumed that the score function and the
score estimator both have Lipschitz regularity in the space variable, and that
the score estimation error is sufficiently small in the $L^{2}$ distance.

Benton, Deligiannidis and Doucet (2024) studied the distribution estimation
error of the flow matching method, but their results rely on the small $L^{2}$
estimation error assumption, the existence and uniqueness of smooth flows
assumption, and the spatial Lipschitzness of estimated velocity field
assumption. Li et al. (2024) derived convergence rates of probability flow
ODEs in the total variation distance, and their results depend on a small
$L^{2}$ score estimation error assumption and a small $L^{2}$ Jacobian
estimation error assumption.

Cui et al. (2023) studied the problem of learning a high-dimensional mixture
of two Gaussians with the flow matching method, in which the velocity field is
parametrized by a two-layer auto-encoder. Furthermore, Cui et al. (2023)
conducted convergence analyses of the Gaussian mixture distribution estimator
in the asymptotic limit $d\to\infty$.

Cheng et al. (2023) presented a theoretical analysis of the distribution
estimator defined by Jordan-Kinderleherer-Otto (JKO) flow models, which
implements the JKO scheme in a normalizing flow network. Gao and Zhu (2024)
assumed a small $L^{2}$ score estimation error, Lipschitz-type time regularity
of the score function, and a smooth log-concave data distribution, and then
studied the distribution estimation error for a general class of probability
flow distribution estimators in the Wasserstein-2 distance.

Finally, Chang et al. (2024) considered a conditional generative learning
model, in which the predictor $\mathsf{X}$ and the response $\mathsf{Y}$ are
both random variables with bounded support. They provided an error analysis
for learning the conditional distribution of $\mathsf{Y}|\mathsf{X}$ via the
Föllmer flow.

In this study, we derive non-asymptotic error bounds for the estimated
velocity fields and discretization error bounds for the forward Euler sampler.
These error bounds are incorporated into the end-to-end convergence analysis
of the CNF distribution estimator with flow matching. Furthermore, we only
stipulate general assumptions on the target distribution, rather than making
assumptions on the velocity field (or score function) and its estimator. We
believe that these theoretical contributions set our work apart from previous
studies.

###  6.2 Diffusion models

Diffusion models (Sohl-Dickstein et al., 2015; Song and Ermon, 2019; Ho, Jain
and Abbeel, 2020; Song et al., 2021; Song, Meng and Ermon, 2021) have emerged
as a powerful SDE-based framework for deep generative learning. The diffusion
model estimators share a deep connection with the CNF distribution estimators
due to the correspondence between an Itô SDE and its probability flow ODE.
There has been a growing interest in the statistical analyses of diffusion
model estimators, as evidenced by the works of Lee, Lu and Tan (2022); Bortoli
(2022); Chen et al. (2023a); Lee, Lu and Tan (2023); Chen, Lee and Lu (2023);
Oko, Akiyama and Suzuki (2023); Li et al. (2024), among others.

Unlike the deterministic sampler of CNF distribution estimators, diffusion
model estimators employ a stochastic sampler (such as the Euler-Maruyama
method) to simulate the time-reversed Itô SDEs. This stochasticity plays a
crucial role in the discretization error analysis of diffusion model
estimators and leads to the development of useful techniques such as
Girsanov’s theorem (Chen et al., 2023a), a chain rule-based variant (Chen, Lee
and Lu, 2023) of the interpolation technique (Vempala and Wibisono, 2019), and
the stochastic interpolation formula (Bortoli, 2022). However, it remains
uncertain whether these techniques can be generalized for analyzing the CNF
distribution estimators.

Compared to the CNF distribution estimators, the diffusion model estimators
have been extensively investigated from a statistical perspective. For
instance, the estimation error bounds of the score function have been
established by Oko, Akiyama and Suzuki (2023); Chen et al. (2023c); Huang et
al. (2023); Cole and Lu (2024). There is also a vast body of literature on
analyzing the discretization error of diffusion model estimators, including
works by Wibisono and Yang (2022); Benton et al. (2024); Pedrotti, Maas and
Mondelli (2023); Gao, Nguyen and Zhu (2023); Bruno et al. (2023); Shah, Chen
and Klivans (2023) and others. However, the absence of stochasticity presents
significant challenges when attempting to analyze the ODE-based CNF
distribution estimators using techniques developed for diffusion model
estimators.

###  6.3 Neural network approximation with Lipschitz regularity control

The approximation theory of deep ReLU networks has developed rapidly since the
seminal work of Yarotsky (2017). Previous studies have shown that deep ReLU
networks can efficiently approximate functions in a smooth function class,
such as the Hölder class, the Sobolev class, and the Besov class, under the
$L^{\infty}$ norm (Yarotsky, 2017; Petersen and Voigtlaender, 2018; Suzuki,
2019; Yarotsky, 2018; Gühring, Kutyniok and Petersen, 2020; DeVore, Hanin and
Petrova, 2021; Daubechies et al., 2022; Lu et al., 2021; Jiao et al., 2023;
Siegel, 2023). Recent works have also considered nonparametric or
semiparametric estimation using deep ReLU networks, including least squares
regression (Bauer and Kohler, 2019; Schmidt-Hieber, 2020; Nakada and Imaizumi,
2020; Kohler and Langer, 2021; Suzuki and Nitanda, 2021; Chen et al., 2022;
Jiao et al., 2023), quantile regression (Shen et al., 2022; Padilla, Tansey
and Chen, 2022) semiparametric inference (Farrell, Liang and Misra, 2021),
factor augmented sparse throughput models (Fan and Gu, 2023), among others. In
the convergence analysis of these models, it is sufficient to know the error
bounds of using deep neural networks for approximating smooth functions.

Analyzing deep generative distribution estimators becomes more challenging as
it requires not only approximation error bounds but also additional regularity
properties of the constructed neural network approximation functions. For
instance, the error analysis of Wasserstein GANs necessitates an upper bound
of the Lipschitz constant of the discriminator network (Chen et al., 2020;
Huang et al., 2022). Chen et al. (2020) demonstrated that the wide and shallow
ReLU network constructed by Yarotsky (2017), for which the depth grows
logarithmically but the width grows polynomially, can approximate 1-Lipschitz
functions with a uniformly bounded Lipschitz constant. Huang et al. (2022)
provided a Lipschitz constant bound for the deep ReLU network approximation
function proposed by Lu et al. (2021). However, this bound increases with the
width and depth of the network. Furthermore, Jiao, Wang and Yang (2023)
succeeded in controlling the Lipschitz constant of deep ReLU networks by
enforcing a norm constraint on the neural network weights, and applied the
approximation bound to analyze the distribution estimation error of GANs. In
addition to the error analyses of GANs, the convergence analysis of
simulation-based CNFs by Marzouk et al. (2023), also requires a Lipschitz
regularity control of the constructed approximation function to ensure the
CNFs are well-posed.

In the current context, the Lipschitz regularity of the neural network
approximation functions is crucial for analyzing the behavior of the estimated
velocity field. Indeed, a key step in our error analysis involves constructing
deep ReLU networks to approximate the Lipschitz velocity field $v^{*}(t,x)$
for $(t,x)\in[0,1-\underline{t}]\times\mathbb{R}^{d}.$ To achieve this target,
we need to derive an $L^{\infty}$ bound of the approximation error and
demonstrate that the Lipschitz constant of the constructed deep ReLU network
is uniformly bounded, regardless of the varying width and depth of the neural
network. Establishing the Lipschitz regularity of the neural network
approximation functions, in addition to the approximation error bounds, is a
more challenging task that requires different techniques. Specifically, our
uniform bounds of the Lipschitz constants are sharper than those obtained by
Huang et al. (2022) for varying width and depth of the deep ReLU network.
Compared to the approximation bound of Chen et al. (2020), our approximation
bound is valid for any network width and depth specified by the parameters $N$
and $L.$ Marzouk et al. (2023) considered the Lipschitz regularity of deep
neural networks activated by the smooth function ReLUk with $k\geq 2$, which
is based on spline approximation and technically differs from our work.

##  7 Conclusion and discussion

We have established non-asymptotic error bounds for the CNF distribution
estimator trained via flow matching, using the Wasserstein-$2$ distance.
Assuming that the target distribution belongs to several rich classes of
probability distributions, we have established Lipschitz regularity properties
of the velocity field for simulation-free CNFs defined with linear
interpolation. To meet the regularity requirements of flow matching
estimators, we have developed $L^{\infty}$ approximation bounds of deep ReLU
networks for Lipschitz functions, along with Lipschitz regularity control of
the constructed deep ReLU networks. By integrating the regularity results, the
deep approximation bounds, and perturbation analyses of ODE flows, we have
shown that the convergence rate of the CNF distribution estimator is
$\widetilde{{\mathcal{O}}}(n^{-1/(d+5)}),$ up to a polylogarithmic prefactor
of $n$. Our error analysis framework can be extended to study more general
CNFs based on interpolation, beyond the CNFs constructed with linear
interpolation.

There are several questions deserving further investigation. Firstly, it would
be interesting to consider target distributions with general smoothness
properties and investigate the resulting smoothness properties of the
corresponding velocity fields. Secondly, the inevitability of the time
singularity of the velocity field remains unclear and warrants further
analysis, as we have not provided a lower bound on the Lipschitz constant in
the time variable. This is a challenging problem that requires more effort and
careful analyses. Lastly, it would be interesting to derive general non-
asymptotic error bounds and convergence rates for CNF distribution estimators
under general smoothness conditions. For this purpose, we need to combine
general smoothness properties of velocity fields with the deep neural network
approximation theory.

##  Appendix A Regularity of the velocity field

In this appendix, we study the regularity properties of the velocity field and
present necessary lemmas, theorems, propositions, and their proofs.

We first list some useful notations. The space ${\mathbb{R}}^{d}$ is endowed
with the Euclidean metric and we denote by $\|\cdot\|_{2}$ and
$\langle\cdot,\cdot\rangle$ the corresponding norm and inner product. Let
$\mathbb{S}^{d-1}:=\\{x\in{\mathbb{R}}^{d}:\|x\|_{2}=1\\}$,
$\mathbb{B}^{d}(x_{0},r,\|\cdot\|_{p}):=\\{x\in{\mathbb{R}}^{d}:\|x-x_{0}\|_{p}<r\\}$,
and
$\bar{\mathbb{B}}^{d}(x_{0},r,\|\cdot\|_{p}):=\\{x\in{\mathbb{R}}^{d}:\|x-x_{0}\|_{p}\leq
r\\}$. The spectral norm of a matrix $A\in{\mathbb{R}}^{k\times d}$ is denoted
by $\|A\|_{2,2}:=\sup_{x\in\mathbb{S}^{d-1}}\|Ax\|_{2}$ and $A^{\top}$ is the
transpose of $A$. For $X,Y\in{\mathbb{R}}$, we denote $X\vee
Y:=\max\\{X,Y\\}$. For two random variables $\mathsf{X}$ and $\mathsf{Y}$, let
$\mathsf{X}\overset{d}{=}\mathsf{Y}$ mean that $\mathsf{X}$ and $\mathsf{Y}$
have the same distribution.

We introduce several auxiliary conditions to assist studying the regularity
properties of the velocity field. These conditions are covered in the three
cases of Assumption 2.

######  Condition 1 (Semi-log-concavity).

Let $\nu(\mathrm{d}x)=\exp(-U(x))\mathrm{d}x$. The potential function $U(x)$
is of class $C^{2}$ and satisfies $\nabla^{2}U(x)\succeq\kappa\mathrm{I}_{d}$
for some $\kappa\in{\mathbb{R}}$.

######  Condition 2 (Bounded support).

The target distribution $\nu$ has bounded support, that is, $D<\infty$.

######  Condition 3 (Gaussian smoothing).

The target distribution $\nu=\gamma_{d,\sigma^{2}}*\rho$ where $\rho$ is a
probability distribution supported on a Euclidean ball of radius $R$ on
$\mathbb{R}^{d}$.

######  Lemma A.1 (Proposition 29 in Gao, Huang and Jiao (2023)).

Let $\nu(\mathrm{d}y)=p(y)\mathrm{d}y$ be a probability distribution on
${\mathbb{R}}^{d}$ with $D:=(1/\sqrt{2})\mathrm{diam}(\mathrm{supp}(\nu))$.

  * (1)

For any $t\in(0,1)$,

$-\frac{1}{1-t}\mathrm{I}_{d}\preceq\nabla_{x}v^{*}(t,x)\preceq\left\\{\frac{t}{(1-t)^{3}}D^{2}-\frac{1}{1-t}\right\\}\mathrm{I}_{d},\quad\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)\preceq
D^{2}\mathrm{I}_{d}.$

  * (2)

Suppose that $p$ is $\beta$-semi-log-convex with $\beta>0$. Then for any
$t\in(0,1)$,

$\nabla_{x}v^{*}(t,x)\succeq\frac{(\beta+1)t-\beta}{\beta(1-t)^{2}+t^{2}}\mathrm{I}_{d},\quad\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)\succeq\frac{(1-t)^{2}}{\beta(1-t)^{2}+t^{2}}\mathrm{I}_{d}.$

  * (3)

Suppose that $p$ is $\kappa$-semi-log-concave with $\kappa\in{\mathbb{R}}$.
Then for any $t\in(t_{0},1)$,

$\nabla_{x}v^{*}(t,x)\preceq\frac{(\kappa+1)t-\kappa}{\kappa(1-t)^{2}+t^{2}}\mathrm{I}_{d},\quad\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)\preceq\frac{(1-t)^{2}}{\kappa(1-t)^{2}+t^{2}}\mathrm{I}_{d},$

where $t_{0}$ is the root of the equation $\kappa+t^{2}/(1-t)^{2}=0$ over
$t\in(0,1)$ if $\kappa<0$ and $t_{0}=0$ if $\kappa\geq 0$.

  * (4)

Fix a probability distribution $\rho$ on ${\mathbb{R}}^{d}$ supported on a
Euclidean ball of radius $R$, and let $\nu:=\gamma_{d,\sigma^{2}}*\rho$ with
$\sigma>0$. Then for any $t\in(0,1)$,

$\displaystyle\frac{(\sigma^{2}+1)t-1}{(1-t)^{2}+\sigma^{2}t^{2}}\mathrm{I}_{d}\preceq\nabla_{x}v^{*}(t,x)\preceq\left\\{\frac{t(1-t)}{((1-t)^{2}+\sigma^{2}t^{2})^{2}}R^{2}+\frac{(\sigma^{2}+1)t-1}{(1-t)^{2}+\sigma^{2}t^{2}}\right\\}\mathrm{I}_{d},$
$\displaystyle~{}~{}\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)\preceq\left\\{\left(\frac{(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\right)^{2}R^{2}+\frac{\sigma^{2}(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\right\\}\mathrm{I}_{d},$
$\displaystyle~{}~{}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)\overset{d}{=}\frac{(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\mathsf{Q}+\sqrt{\frac{\sigma^{2}(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}}\mathsf{Z}+\frac{\sigma^{2}t^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}x$

where $\mathsf{Q}\sim\tilde{\rho}$ is supported on the same ball as $\rho$,
$\mathsf{Z}\sim\gamma_{d}$, and $\mathsf{Q},\mathsf{Z}$ are independent.

In Lemma A.2 below, we show that the velocity field and its spatial derivative
is (locally) bounded under mild regularity conditions. The boundedness of the
spatial derivative directly follows from Lemma A.1. Since a Lipschitz property
results in a linear growth property, we obtain the the velocity field is
locally bounded. For ease of presentation, let us define two parameter sets by

$\displaystyle\mathcal{S}_{1}:=\begin{cases}\\{\kappa,\beta\\}&\text{if
Assumption \ref{assump:target}-(i) holds},\\\ \\{D,\kappa,\beta\\}&\text{if
Assumption \ref{assump:target}-(ii) holds},\\\ \\{R,\sigma\\}&\text{if
Assumption \ref{assump:target}-(iii) holds},\end{cases}$
$\displaystyle\mathcal{S}_{2}:=\begin{cases}\\{d,\kappa,\beta\\}&\text{if
Assumption \ref{assump:target}-(i) holds},\\\ \\{d,D,\kappa,\beta\\}&\text{if
Assumption \ref{assump:target}-(ii) holds},\\\ \\{d,R,\sigma\\}&\text{if
Assumption \ref{assump:target}-(iii) holds}.\end{cases}$

We say a prefactor scales polynomially with $\mathcal{S}_{1}$ if it scales
polynomially with parameters in $\mathcal{S}_{1}$.

######  Lemma A.2.

Suppose that Assumptions 1 and 2 hold. Then it holds that

$\sup_{(t,x)\in[0,1]\times\Omega_{A}}\|v^{*}(t,x)\|_{2}\lesssim
A,\quad\sup_{(t,x)\in[0,1]\times\mathbb{R}^{d}}\|\nabla_{x}v^{*}(t,x)\|_{2,2}\lesssim
1,$

where we omit prefactors scaling polynomially with $\mathcal{S}_{2}$.

###### Proof.

Under Assumptions 1 and 2, Lemma A.1 shows that

$\displaystyle
C_{1}(\mathcal{S}_{1})\mathrm{I}_{d}\preceq\nabla_{x}v^{*}(t,x)\preceq
C_{2}(\mathcal{S}_{1})\mathrm{I}_{d},$

where $C_{1}(\mathcal{S}_{1})$ and $C_{2}(\mathcal{S}_{1})$ are constants
scaling polynomially with $\mathcal{S}_{1}$. It further yields that

$\sup_{(t,x)\in[0,1]\times\mathbb{R}^{d}}\|\nabla_{x}v^{*}(t,x)\|_{2,2}\lesssim
1,$ (A.1)

when we omit a prefactor scaling polynomially with $\mathcal{S}_{1}$. Notice
that for any $t\in(0,1)$, it holds that

$\displaystyle v^{*}(t,0)$
$\displaystyle=\frac{1}{1-t}\mathbb{E}[\mathsf{X}_{1}|\mathsf{X}_{t}=0]=\frac{1}{1-t}\int_{\mathbb{R}^{d}}yq(y|t,0)\mathrm{d}y$
$\displaystyle\lesssim\frac{1}{1-t}\int_{{\mathbb{R}}^{d}}yp(y)(1-t)^{-d}\exp\left(-\frac{t^{2}\|y\|^{2}_{2}}{2(1-t)^{2}}\right)\mathrm{d}y,$

which implies $\|v^{*}(t,0)\|_{2}<\infty$ due to fast growth of the
exponential function. Besides,
$v^{*}(0,0)=\mathbb{E}[\mathsf{X}_{1}],v^{*}(1,0)=0$. Then by the boundedness
of $\|v^{*}(t,0)\|_{2}$ over $[0,1]$ and (A.1), we bound $v^{*}(t,x)$ as
follows

$\displaystyle\|v^{*}(t,x)\|_{2}$
$\displaystyle\leq\|v^{*}(t,0)\|_{2}+\|v^{*}(t,x)-v^{*}(t,0)\|_{2}$
$\displaystyle\leq\|v^{*}(t,0)\|_{2}+\left\\{\sup_{(t,y)\in[0,1]\times{\mathbb{R}}^{d}}\|\nabla_{y}v^{*}(t,y)\|_{2,2}\right\\}\|x\|_{2}$
$\displaystyle\lesssim\|x\|_{2}\vee 1,$

where we omit a prefactor scaling polynomially with $\mathcal{S}_{1}$. It
further yields that

$\sup_{(t,x)\in[0,1]\times\Omega_{A}}\|v^{*}(t,x)\|_{2}\lesssim A$

by omitting a prefactor scaling polynomially with $\mathcal{S}_{2}$. This
completes the proof. ∎

###  A.1 Control with semi-log-concavity

We derive moment bounds under Condition 1. The moment bounds are useful to
estimate the time regularity of the velocity field.

######  Lemma A.3 (Moment bounds).

Suppose that Condition 1 holds. Let $\eta\in(0,1)$ be a constant. Let $t_{1}$
be the root of the equation $\kappa(1-t)^{2}+t^{2}=\eta$ over $t\in(0,1)$ if
$\kappa\leq 0$ or $t_{1}=0$ if $\kappa>0$. Then for any
$t\in[t_{1},1-\underline{t}]$, it holds that

$\displaystyle\sup_{x\in\Omega_{A}}\|M_{1}\|_{2}\lesssim
A,\quad\sup_{x\in{\mathbb{R}}^{d}}\|M^{c}_{2}\|_{2,2}\lesssim(1-t)^{2},\quad\sup_{x\in\Omega_{A}}\|M_{3}-M_{2}M_{1}\|_{2}\lesssim
A(1-t)^{2},$

where we omit polynomial prefactors in $d,\kappa,\eta$.

###### Proof.

First, we bound $M^{c}_{2}$. According to Lemma A.1, the following covariance
bound holds for any $t\in[t_{1},1)$

$0\mathrm{I}_{d}\preceq\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)\preceq\frac{(1-t)^{2}}{\kappa(1-t)^{2}+t^{2}}\mathrm{I}_{d}\quad\text{with
}\kappa(1-t)^{2}+t^{2}\geq\begin{cases}C_{\kappa},\quad&\text{if
$\kappa>0$},\\\ \eta,\quad&\text{if $\kappa\leq 0$},\end{cases}$

where $C_{\kappa}:=\kappa/(\kappa+1)$. Then it implies that for any
$t\in[t_{1},1-\underline{t}]$,
$\sup_{x\in{\mathbb{R}}^{d}}\|M^{c}_{2}\|_{2,2}\lesssim(1-t)^{2}$ with
omitting a polynomial prefactor in $\kappa,\eta$.

Then, we bound $M_{1}$ and $M_{2}$. By the Hatsell-Nolte identity (Dytso, Poor
and Shamai Shitz, 2023, Proposition 1), we obtain
$\nabla_{x}M_{1}(t,x)=(t/(1-t)^{2})M^{c}_{2}$ which implies that

$\sup_{(t,x)\in[t_{1},1-\underline{t}]\times{\mathbb{R}}^{d}}\|\nabla_{x}M_{1}(t,x)\|_{2,2}=\sup_{(t,x)\in[t_{1},1-\underline{t}]\times{\mathbb{R}}^{d}}\frac{t}{(1-t)^{2}}\|M^{c}_{2}\|_{2,2}\lesssim
1$ (A.2)

with a polynomial prefactor in $\kappa,\eta$ hidden. Notice that for any
$t\in(0,1)$, it holds that

$\displaystyle
M_{1}(t,0)=\int_{\mathbb{R}^{d}}yq(y|t,0)\mathrm{d}y\lesssim\int_{{\mathbb{R}}^{d}}yp(y)(1-t)^{-d}\exp\left(-\frac{t^{2}\|y\|^{2}_{2}}{2(1-t)^{2}}\right)\mathrm{d}y,$

which implies $\|M_{1}(t,0)\|_{2}<\infty$ due to fast growth of the
exponential function. Besides, $M_{1}(0,0)=\mathbb{E}_{p}[\mathsf{X}_{1}]$ and
$M_{1}(1,0)=0$. By the boundedness of $\|M_{1}(t,0)\|_{2}$ for any $t\in[0,1]$
and (A.2), we further bound $M_{1}(t,x)$ for any
$(t,x)\in[t_{1},1-\underline{t}]\times{\mathbb{R}}^{d}$ as follows

$\displaystyle\|M_{1}(t,x)\|_{2}$
$\displaystyle\leq\|M_{1}(t,0)\|_{2}+\|M_{1}(t,x)-M_{1}(t,0)\|_{2}$
$\displaystyle\leq\|M_{1}(t,0)\|_{2}+\left\\{\sup_{(t,y)\in[t_{1},1-\underline{t}]\times{\mathbb{R}}^{d}}\|\nabla_{y}M_{1}(t,y)\|_{2,2}\right\\}\|x\|_{2}$
$\displaystyle\lesssim\|x\|_{2}\vee 1,$

where a polynomial prefactor in $\kappa,\eta$ is hidden. It further yields
that

$\sup_{(t,x)\in[t_{1},1-\underline{t}]\times\Omega_{A}}\|M_{1}\|_{2}\lesssim
A$

when omitting a polynomial prefactor in $d,\kappa,\eta$. Moreover, notice that
$M_{2}=\operatorname{Tr}(M^{c}_{2})+\|M_{1}\|_{2}^{2}$, which further yields
that

$\sup_{(t,x)\in[t_{1},1-\underline{t}]\times\Omega_{A}}|M_{2}|\lesssim A^{2}$

with an omitted polynomial prefactor in $d,\kappa,\eta$.

Lastly, we bound $M_{3}-M_{2}M_{1}$. For any $i\in\\{1,2,\cdots,d\\}$, let
$\mathsf{X}_{1,i}$ denote the $i$-th element of $\mathsf{X}_{1}$. Then it
holds that

$\displaystyle~{}~{}~{}~{}~{}~{}\|M_{3}-M_{2}M_{1}\|_{2}^{2}$
$\displaystyle=\sum\nolimits_{i=1}^{d}\left(\mathbb{E}[\mathsf{X}_{1,i}\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x]-\mathbb{E}[\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x]\mathbb{E}[\mathsf{X}_{1,i}|\mathsf{X}_{t}=x]\right)^{2}$
$\displaystyle=\sum\nolimits_{i=1}^{d}\left(\mathrm{Cov}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1},\mathsf{X}_{1,i}|\mathsf{X}_{t}=x)\right)^{2}$
$\displaystyle\leq\sum\nolimits_{i=1}^{d}\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)\mathrm{Var}(\mathsf{X}_{1,i}|\mathsf{X}_{t}=x)$
(By the Cauchy-Schwarz inequality)
$\displaystyle=\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)\sum\nolimits_{i=1}^{d}\mathrm{Var}(\mathsf{X}_{1,i}|\mathsf{X}_{t}=x)$
$\displaystyle=\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)\operatorname{Tr}(M^{c}_{2})$
$\displaystyle\leq d\
\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)\|M^{c}_{2}\|_{2,2}.$

Let $\mathsf{X}_{1}\sim p(y)$ be $\kappa$-semi-log-concave for some
$\kappa\in{\mathbb{R}}$. Then for any $t\in[0,1)$,
$\mathsf{X}_{1}|\mathsf{X}_{t}\sim q(y|t,x)$ is
$(\kappa+t^{2}/(1-t)^{2})$-semi-log-concave because

$-\nabla^{2}_{y}\log q(y|t,x)=-\nabla^{2}_{y}\log p(y)-\nabla^{2}_{y}\log
q(t,x|y)\succeq\left(\kappa+\frac{t^{2}}{(1-t)^{2}}\right)\mathrm{I}_{d}.$

When $t\in\left\\{t:\kappa+t^{2}/(1-t)^{2}>0,t\in(0,1)\right\\}$, by the
Brascamp-Lieb inequality (Brascamp and Lieb, 1976), it yields that

$\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)\leq
4M_{2}\left(\kappa+\frac{t^{2}}{(1-t)^{2}}\right)^{-1}=4M_{2}\frac{(1-t)^{2}}{\kappa(1-t)^{2}+t^{2}}.$

Analogous to the control of $\|M^{c}_{2}\|_{2,2}$, we further obtain that for
any $(t,x)\in[t_{1},1-\underline{t}]\times\Omega_{A}$,

$\displaystyle\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)\lesssim\begin{cases}A^{2}(1-t)^{2}/C_{\kappa},\quad&\text{if
$\kappa>0$},\\\ A^{2}(1-t)^{2}/\eta,\quad&\text{if $\kappa\leq
0$}.\end{cases}$

Hence, we deduce that for any $t\in[t_{1},1-\underline{t}]$,

$\sup_{x\in\Omega_{A}}\|M_{3}-M_{2}M_{1}\|_{2}\lesssim A(1-t)^{2},$

where we omit a polynomial prefactor in $d,\kappa,\eta$. This completes the
proof. ∎

######  Lemma A.4.

Suppose that Condition 1 holds. Then it holds that

$\sup_{(t,x)\in[t_{1},1-\underline{t}]\times\Omega_{A}}\|\partial_{t}v^{*}(t,x)\|_{2}\lesssim
A/\underline{t}^{2},$

where we omit a polynomial prefactor in $d,\kappa,\eta$.

###### Proof.

By Lemma 3.4, it holds that

$\displaystyle\|\partial_{t}v^{*}(t,x)\|_{2}\leq$
$\displaystyle\frac{1}{(1-t)^{2}}\|x\|_{2}+\frac{1}{(1-t)^{2}}\|M_{1}\|_{2}+\frac{1}{(1-t)^{4}}\|M^{c}_{2}\|_{2,2}\cdot\|x\|_{2}$
$\displaystyle+\frac{1}{(1-t)^{4}}\|M_{3}-M_{2}M_{1}\|_{2}.$

Applying Lemma A.3, we obtain

$\displaystyle\sup_{(t,x)\in[t_{1},1-\underline{t}]\times\Omega_{A}}\|\partial_{t}v^{*}(t,x)\|_{2}\lesssim\frac{A}{\underline{t}^{2}},$

where we omit a polynomial prefactor in $d,\kappa,\eta$. ∎

###  A.2 Control with bounded support

We derive moment bounds under Condition 2. The moment bounds are useful to
estimate the time regularity of the velocity field.

######  Lemma A.5 (Moment bounds).

Suppose that Condition 2 holds. Then it holds that

$\displaystyle\sup_{(t,x)\in[0,1]\times\mathbb{R}^{d}}\|M_{1}\|_{2}\lesssim
D,\sup_{(t,x)\in[0,1]\times\mathbb{R}^{d}}\|M^{c}_{2}\|_{2,2}\lesssim
D^{2},\sup_{(t,x)\in[0,1]\times\mathbb{R}^{d}}\|M_{3}-M_{2}M_{1}\|_{2}\lesssim
D^{3}.$

###### Proof.

The desired bounds hold naturally due to the boundedness of $\mathsf{X}_{1}$.
∎

To ease the notation, we denote
$\Omega_{\underline{t},A}:=[0,1-\underline{t}]\times[-A,A]^{d}$.

######  Lemma A.6.

Suppose that Condition 2 holds. Then it holds that

$\displaystyle\sup_{(t,x)\in\Omega_{\underline{t},A}}\|v^{*}(t,x)\|_{2}\lesssim
A/\underline{t},\sup_{(t,x)\in[0,1-\underline{t}]\times\mathbb{R}^{d}}\|\nabla_{x}v^{*}(t,x)\|_{2,2}\lesssim
1/\underline{t}^{3},$
$\displaystyle\sup_{(t,x)\in\Omega_{\underline{t},A}}\|\partial_{t}v^{*}(t,x)\|_{2}\lesssim
A/\underline{t}^{4},$

where we omit polynomial prefactors in $d,D$.

###### Proof.

By Eq. (3.7), for $t\in(0,1)$, it implies that

$\displaystyle\|v^{*}(t,x)\|_{2}\leq\frac{1}{1-t}\|x\|_{2}+\frac{1}{1-t}\|M_{1}\|_{2}\leq\frac{\|x\|_{2}\vee\|M_{1}\|_{2}}{1-t},$

which yields that

$\sup_{(t,x)\in\Omega_{\underline{t},A}}\|v^{*}(t,x)\|_{2}\lesssim\frac{A}{\underline{t}}$

with a polynomial prefactor in $D$ omitted. By Lemma 3.3, for $t\in(0,1)$, it
holds that

$\displaystyle\|\nabla_{x}v^{*}(t,x)\|_{2,2}\leq\frac{t}{(1-t)^{3}}\|M^{c}_{2}\|_{2,2}+\frac{1}{1-t}\lesssim\frac{1}{(1-t)^{3}},$

which implies that

$\sup_{(t,x)\in[0,1-\underline{t}]\times\mathbb{R}^{d}}\|\nabla_{x}v^{*}(t,x)\|_{2,2}\lesssim\frac{1}{\underline{t}^{3}}$

with a polynomial prefactor in $D$ omitted. By Lemma 3.4, the following bounds
hold

$\displaystyle\|\partial_{t}v^{*}(t,x)\|_{2}\leq$
$\displaystyle\frac{1}{(1-t)^{2}}\|x\|_{2}+\frac{1}{(1-t)^{2}}\|M_{1}\|_{2}+\frac{1}{(1-t)^{4}}\|M^{c}_{2}\|_{2,2}\cdot\|x\|_{2}$
$\displaystyle+\frac{1}{(1-t)^{4}}\|M_{3}-M_{2}M_{1}\|_{2}$
$\displaystyle\lesssim$
$\displaystyle\frac{1}{(1-t)^{4}}(\|x\|_{2}+\|M_{1}\|_{2}+\|M^{c}_{2}\|_{2,2}\cdot\|x\|_{2}+\|M_{3}-M_{2}M_{1}\|_{2}),$

which implies that

$\sup_{(t,x)\in\Omega_{\underline{t},A}}\|\partial_{t}v^{*}(t,x)\|_{2}\lesssim\frac{A}{\underline{t}^{4}}$

with a polynomial prefactor in $d,D$ omitted. This completes the proof. ∎

###  A.3 Control with Gaussian smoothing

We derive moment bounds under Condition 3. The moment bounds are useful to
estimate the time regularity of the velocity field.

######  Lemma A.7 (Moment bounds).

Suppose that Condition 3 holds. Then for any $t\in[0,1-\underline{t}]$, it
holds that

$\displaystyle\sup_{x\in\Omega_{A}}\|M_{1}\|_{2}\lesssim
A,\quad\sup_{x\in{\mathbb{R}}^{d}}\|M^{c}_{2}\|_{2,2}\lesssim(1-t)^{2},\quad\sup_{x\in\Omega_{A}}\|M_{3}-M_{2}M_{1}\|_{2}\lesssim
A(1-t)^{2},$

where we omit polynomial prefactors in $d,R,\sigma$.

###### Proof.

The proof idea is partially similar to that of Lemma A.3.

First, we bound $M^{c}_{2}$. According to Lemma A.1, the following covariance
bound holds for any $t\in[0,1)$,

$0\mathrm{I}_{d}\preceq\mathrm{Cov}(\mathsf{X}_{1}|\mathsf{X}_{t}=x)\preceq(1-t)^{2}\left\\{\frac{R^{2}(1-t)^{2}}{((1-t)^{2}+\sigma^{2}t^{2})^{2}}+\frac{\sigma^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\right\\}\mathrm{I}_{d}.$

Notice that

$\displaystyle\frac{R^{2}(1-t)^{2}}{((1-t)^{2}+\sigma^{2}t^{2})^{2}}+\frac{\sigma^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\leq\left(1+\frac{1}{\sigma^{2}}\right)^{2}R^{2}+\sigma^{2}+1.$

It implies that for any $t\in[0,1-\underline{t}]$,

$\sup_{x\in{\mathbb{R}}^{d}}\|M^{c}_{2}\|_{2,2}\lesssim(1-t)^{2}$ (A.3)

with omitting a polynomial prefactor in $R,\sigma$.

Then, we bound $M_{1}$ and $M_{2}$. Again, by the Hatsell-Nolte identity
(Dytso, Poor and Shamai Shitz, 2023, Proposition 1), we obtain that

$\sup_{(t,x)\in[0,1-\underline{t}]\times{\mathbb{R}}^{d}}\|\nabla_{x}M_{1}(t,x)\|_{2,2}=\sup_{(t,x)\in[0,1-\underline{t}]\times{\mathbb{R}}^{d}}\frac{t}{(1-t)^{2}}\|M^{c}_{2}\|_{2,2}\lesssim
1$ (A.4)

with a polynomial prefactor in $R,\sigma$ hidden. Identical to how we proceed
in the proof of Lemma A.3, we have

$\sup_{(t,x)\in\Omega_{\underline{t},A}}\|M_{1}\|_{2}\lesssim A$

when omitting a polynomial prefactor in $d,R,\sigma$.

Finally, we bound $M_{3}-M_{2}M_{1}$. Recall that we have deduced the
following inequality in the proof of Lemma A.3

$\displaystyle\|M_{3}-M_{2}M_{1}\|_{2}^{2}\leq d\
\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)\|M^{c}_{2}\|_{2,2}.$
(A.5)

We next focus on bounding
$\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)$. By Lemma
A.1-(4), it is shown that

$\displaystyle(\mathsf{X}_{1}|\mathsf{X}_{t}=x)\overset{d}{=}\mathsf{P}_{x}:=\frac{(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\mathsf{Q}+\sqrt{\frac{\sigma^{2}(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}}\mathsf{Z}+\frac{\sigma^{2}t^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}x$

where $\mathsf{Q}\sim\tilde{\rho}$ is supported on the same ball as $\rho$,
$\mathsf{Z}\sim\gamma_{d}$, and $\mathsf{Q},\mathsf{Z}$ are independent. In
the expression above, we note that the denominator $(1-t)^{2}+\sigma^{2}t^{2}$
is lower bounded by $\sigma^{2}/(\sigma^{2}+1)$ over $t\in[0,1]$. Let
$\mathsf{R}_{x,y}:=\mathsf{P}_{x}^{\top}\mathsf{P}_{x}|\mathsf{Q}=y$. Then by
the law of total variance, it yields that

$\displaystyle\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)=\mathrm{Var}(\mathsf{P}_{x}^{\top}\mathsf{P}_{x})=\mathbb{E}[\mathrm{Var}(\mathsf{R}_{x,y})]+\mathrm{Var}(\mathbb{E}[\mathsf{R}_{x,y}]).$
(A.6)

We claim that $\mathsf{R}_{x,y}/\eta$ is distributed as a noncentral chi-
squared distribution with degrees of freedom $d$ and the noncentrality
parameter $\xi_{x,y}$ where

$\displaystyle\eta=\frac{\sigma^{2}(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}},\quad\xi_{x,y}=\frac{1}{\eta}\left\|\frac{(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}y+\frac{\sigma^{2}t^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}x\right\|_{2}^{2}.$

By properties of the noncentral chi-squared distribution, it holds that

$\displaystyle\mathbb{E}(\mathsf{R}_{x,y})=\eta(d+\xi_{x,y}),\quad\mathrm{Var}(\mathsf{R}_{x,y})=2\eta^{2}(d+2\xi_{x,y}).$

Then we bound the first term in the variance decomposition (A.6) as follows

$\displaystyle\mathbb{E}[\mathrm{Var}(\mathsf{R}_{x,y})]=2\eta\mathbb{E}\left[\eta
d+2\eta\xi_{x,y})\right]\lesssim(1-t)^{2}(\|x\|_{2}^{2}\vee 1)$

where we omit a polynomial prefactor in $d,R,\sigma$. To bound the second
term, we do the following calculations

$\displaystyle\mathrm{Var}(\mathbb{E}[\mathsf{R}_{x,y}])$
$\displaystyle=\mathrm{Var}(\eta d+\eta\xi_{x,y})=\mathrm{Var}\left(\eta
d+\left\|\frac{(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\mathsf{Q}+\frac{\sigma^{2}t^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}x\right\|_{2}^{2}\right)$
$\displaystyle=\mathrm{Var}\left(\left\|\frac{(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\mathsf{Q}\right\|_{2}^{2}+2\left\langle\frac{(1-t)^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}\mathsf{Q},\frac{\sigma^{2}t^{2}}{(1-t)^{2}+\sigma^{2}t^{2}}x\right\rangle\right)$
$\displaystyle\lesssim(1-t)^{4}(\|x\|_{2}^{2}\vee 1),$

where we omit a polynomial prefactor in $d,R,\sigma$. Combining the control of
the two terms, we obtain that

$\displaystyle\mathrm{Var}(\mathsf{X}_{1}^{\top}\mathsf{X}_{1}|\mathsf{X}_{t}=x)\lesssim(1-t)^{2}(\|x\|_{2}^{2}\vee
1)$ (A.7)

by omitting a polynomial prefactor in $d,R,\sigma$. Therefore, using (A.3),
(A.5), and (A.7), the bound of $M_{3}-M_{2}M_{1}$ is deduced for any
$t\in[0,1-\underline{t}]$ by

$\displaystyle\sup_{x\in\Omega_{A}}\|M_{3}-M_{2}M_{1}\|_{2}\lesssim
A(1-t)^{2}$

with a polynomial prefactor in $d,R,\sigma$ hidden. This completes the proof.
∎

######  Lemma A.8.

Suppose that Condition 3 holds. Then it holds that

$\sup_{(t,x)\in[0,1-\underline{t}]\times\Omega_{A}}\|\partial_{t}v^{*}(t,x)\|_{2}\lesssim
A/\underline{t}^{2},$

where we omit a polynomial prefactor in $d,\kappa,\eta$.

###### Proof.

Based on Lemma A.7, the proof is almost identical to that of Lemma A.4. ∎

###  A.4 Sharpness of moment bounds

The moment bounds in Lemmas A.3 and A.7 are sharp in $(t,x)$ because of a
Gaussian example.

######  Proposition A.9.

Let $\mathsf{X}_{1}\sim\gamma_{d}$. The conditional distribution of
$\mathsf{X}_{1}|\mathsf{X}_{t}$ has the following explicit expression

$\displaystyle\mathsf{X}_{1}|\mathsf{X}_{t}=x\sim
N\left(\frac{t}{(1-t)^{2}+t^{2}}x,\frac{(1-t)^{2}}{(1-t)^{2}+t^{2}}\mathrm{I}_{d}\right).$

Moreover, for any $t\in(0,1]$, the moment bounds are given by

$\displaystyle\sup_{x\in\Omega_{A}}\|M_{1}\|_{2}\asymp A,\
\sup_{x\in{\mathbb{R}}^{d}}\|M^{c}_{2}\|_{2,2}\asymp(1-t)^{2},\
\sup_{x\in\Omega_{A}}\|M_{3}-M_{2}M_{1}\|_{2}\asymp A(1-t)^{2},$

where we omit polynomial prefactors in $d$.

###### Proof.

By Bayes’ rule, for $\mathsf{X}_{1}\sim\gamma_{d}$, it implies that

$\displaystyle\mathsf{X}_{1}|\mathsf{X}_{t}=x\sim
N\left(\frac{t}{(1-t)^{2}+t^{2}}x,\frac{(1-t)^{2}}{(1-t)^{2}+t^{2}}\mathrm{I}_{d}\right).$

By properties of the Gaussian distribution, the desired moment bounds hold. ∎

###  A.5 Proof of Theorem 5.1

###### Proof of Theorem 5.1.

By Lemma A.2, it holds that for any $x,y\in{\mathbb{R}}^{d}$ and $t\in[0,1]$,
$\|v^{*}(t,x)-v^{*}(t,y)\|_{\infty}\lesssim\|x-y\|_{\infty}$, and that
$\sup_{(t,x)\in[0,1]\times\Omega_{A}}\|v^{*}(t,x)\|_{\infty}\lesssim A$, where
we omit constants in $d,\kappa,\beta,\sigma,D,R$.

Then we show that the Lipschitz continuity of $v^{*}(t,x)$ in $t$. Concretely,
for any $s,t\in[0,1-\underline{t}]$ and $x\in{\mathbb{R}}^{d}$,
$\|v^{*}(t,x)-v^{*}(s,x)\|_{\infty}\leq L_{t}|t-s|$ with
$L_{t}\lesssim\underline{t}^{-2}$ by omitting a constant in
$d,\kappa,\beta,\sigma,D,R$. We analyze the cases in Assumption 2 one by one
as follows:

  * •

Suppose that Assumptions 1 and 2-(i) holds. Condition 1 holds as well. We use
controls with semi-log-concavity in Lemma A.4 and derive the desired Lipschitz
continuity in $t$.

  * •

Suppose that Assumptions 1 and 2-(ii) holds. Conditions 1 and 2 hold as well.
We use controls with bounded support in Lemma A.6 for $t\in[0,t_{1}]$ and use
controls with semi-log-concavity in Lemma A.4 for
$t\in[t_{1},1-\underline{t}]$. Then we derive the desired Lipschitz continuity
in $t$.

  * •

Suppose that Assumptions 1 and 2-(iii) holds. Condition 3 holds as well. We
use controls with with Gaussian smoothing in Lemma A.8, and the Lipschitz
continuity in $t$ follows.

We complete the proof. ∎

##  Appendix B Approximation error of the velocity field

In this section, we analyze the approximation error of the velocity field by a
constructive approach.

Before proceeding, we present a few useful notations for the Sobolev function
class. A $d$-dimensional multi-index is a $d$-tuple
$\alpha=(\alpha_{1},\alpha_{2},\cdots,\alpha_{d})^{\top}\in\mathbb{N}_{0}^{d}$.
We define $\|\alpha\|_{1}=\sum_{i=1}^{d}\alpha_{i}$ and
$\partial^{\alpha}:=\partial_{1}^{\alpha_{1}}\partial_{2}^{\alpha_{2}}\cdots\partial_{d}^{\alpha_{d}}$
to represent the partial derivative of a $d$-dimensional function. We also use
$D$ to denote the weak derivative of a single variable function and
$D^{\alpha}$ to denote the partial derivative
$D^{\alpha_{1}}_{1}D^{\alpha_{2}}_{2}\cdots D^{\alpha_{d}}_{d}$ of a
$d$-dimensional function with $\alpha_{i}$ as the order of derivative $D_{i}$
in the $i$-th variable.

###  B.1 Approximation in space with Lipschitz regularity

In this subsection, we study the approximation capacity of deep ReLU networks
joint with an estimate of the Lipschitz regularity. The strong expressive
power of deep ReLU networks has been studied with the localized or averaged
Taylor polynomials. We follow the localized approximation approach, and
establish the global Lipschitz continuity and non-asymptotic approximation
estimate of deep ReLU networks.

######  Lemma B.1.

Given any $f\in W^{1,\infty}((0,1)^{d})$ with
$\|f\|_{W^{1,\infty}((0,1)^{d})}\leq 1$, for any $N,L\in\mathbb{N}$, there
exists a function $\phi$ implemented by a deep ReLU network with width
$\mathcal{O}(2^{d}dN\log N)$ and depth $\mathcal{O}(d^{2}L\log L)$ such that
$\|\phi\|_{W^{1,\infty}((0,1)^{d})}\lesssim 1$ and

$\displaystyle\|\phi-f\|_{L^{\infty}([0,1]^{d})}\lesssim(NL)^{-2/d},$

where we omit some prefactors depending only on $d$.

######  Corollary B.2.

Given any $f\in W^{1,\infty}((0,1)^{d})$ with
$\|f\|_{W^{1,\infty}((0,1)^{d})}<\infty$, for any $N,L\in\mathbb{N}$, there
exists a function $\phi$ implemented by a deep ReLU network with width
$\mathcal{O}(2^{d}dN\log N)$ and depth $\mathcal{O}(d^{2}L\log L)$ such that
$\|\phi\|_{W^{1,\infty}((0,1)^{d})}\lesssim\|f\|_{W^{1,\infty}((0,1)^{d})}$
and

$\displaystyle\|\phi-f\|_{L^{\infty}([0,1]^{d})}\lesssim\|f\|_{W^{1,\infty}((0,1)^{d})}(NL)^{-2/d},$

where we omit some prefactors depending only on $d$.

######  Remark 16.

The approximation rate is nearly optimal for the unit ball of functions in
$W^{1,\infty}((0,1)^{d})$ according to Shen, Yang and Zhang (2020, 2022) and
Lu et al. (2021).

Proof sketch of Lemma B.1. The proof idea is similar to that of Yang, Yang and
Xiang (2023, Theorem 3), and we divide the proof into three steps.

Step 1. Discretization. We use a partition of unity to discretize the set
$(0,1)^{d}$. As in Definitions B.3 and B.4, we construct a partition of unity
$\\{g_{m}\\}_{m\in\\{1,2\\}^{d}}$ on $(0,1)^{d}$ with
$\mathrm{supp}(g_{m})\cap(0,1)^{d}\subset\Omega_{m}$ for any
$m\in\\{1,2\\}^{d}$. Then we approximate the partition of unity
$\\{g_{m}\\}_{m\in\\{1,2\\}^{d}}$ by a collection of deep ReLU networks
$\\{\phi_{m}\\}_{m\in\\{1,2\\}^{d}}$ as in Lemma B.5.

Step 2. Approximation on $\Omega_{m}$. Given any $m\in\\{1,2\\}^{d}$, for each
subset $\Omega_{m}\subset[0,1]^{d}$, we find a piecewise constant function
$f_{K,m}$ satisfying

$\displaystyle\|f_{K,m}-f\|_{W^{1,\infty}(\Omega_{m})}\lesssim
1,\quad\|f_{K,m}-f\|_{L^{\infty}(\Omega_{m})}\lesssim 1/K,$

where we omit constants in $d$. Piecewise constant functions can be
approximated by deep ReLU networks. Then, following Lu et al. (2021) and Yang,
Yang and Xiang (2023), we construct a deep ReLU network $\psi_{m}$ with width
$\mathcal{O}(2^{d}dN\log N)$ and depth $\mathcal{O}(d^{2}L\log L)$ such that

$\displaystyle\|\psi_{m}-f\|_{W^{1,\infty}(\Omega_{m})}\lesssim
1,\quad\|\psi_{m}-f\|_{L^{\infty}(\Omega_{m})}\lesssim(NL)^{-2/d},$

where we omit constants in $d$.

Step 3. Approximation on $[0,1]^{d}$. Combining the approximations on each
subset $\Omega_{m}$ properly, we construct an approximation of the target
function $f$ on the domain $[0,1]^{d}$. That is, for any $N,L\in\mathbb{N}$,
there exists a function $\phi$ implemented by a deep ReLU network with width
$\mathcal{O}(N\log N)$ and depth $\mathcal{O}(L\log L)$ such that

$\displaystyle\|\phi-f\|_{L^{\infty}([0,1]^{d})}\lesssim(NL)^{-2/d}\quad\text{with}\quad\|\phi\|_{W^{1,\infty}((0,1)^{d})}\lesssim
1,$

where we omit constants in $d$.

######  Definition B.3.

Given $K,d\in\mathbb{N}$, and for any
$m=[m_{1},m_{2},\cdots,m_{d}]^{\top}\in\\{1,2\\}^{d}$, we define
$\Omega_{m}:=\prod_{i=1}^{d}\Omega_{m_{j}}$ where
$\Omega_{1}:=\bigcup_{i=1}^{K-1}\left[\frac{i}{K},\frac{i}{K}+\frac{3}{4K}\right]$
and
$\Omega_{2}:=\bigcup_{i=0}^{K}\left[\frac{i}{K}-\frac{1}{2K},\frac{i}{K}+\frac{1}{4K}\right]\cap[0,1]$.

Figure 1: Functions $g_{1}$ and $g_{2}$ for defining a partition of unity.

######  Definition B.4.

Given $K,d\in\mathbb{N}$, for any integer $i\in\mathbb{Z}$, we define

$\displaystyle
g_{1}(x):=\begin{cases}1,&x\in\left[\frac{i}{K}+\frac{1}{4K},\frac{i}{K}+\frac{1}{2K}\right],\\\
0,&x\in\left[\frac{i}{K}+\frac{3}{4K},\frac{i}{K}+\frac{1}{K}\right],\\\
4K\left(x-\frac{i}{K}\right),&x\in\left[\frac{i}{K},\frac{i}{K}+\frac{1}{4K}\right],\\\
-4K\left(x-\frac{i}{K}-\frac{3}{4K}\right),&x\in\left[\frac{i}{K}+\frac{1}{2K},\frac{i}{K}+\frac{3}{4K}\right],\end{cases}\quad
g_{2}(x):=g_{1}\left(x+\frac{1}{2K}\right).$

For any $m=[m_{1},m_{2},\cdots,m_{d}]^{\top}\in\\{1,2\\}^{d}$, we further
define $g_{m}(x):=\prod_{j=1}^{d}g_{m_{j}}(x_{j})$ where
$x=[x_{1},x_{2},\cdots,x_{d}]^{\top}$.

######  Lemma B.5 (Proposition 1 in Yang, Yang and Xiang (2023)).

Given any $N,L\in\mathbb{N}$ and any $m\in\\{1,2\\}^{d}$, for $K=\lfloor
N^{1/d}\rfloor^{2}\lfloor L^{2/d}\rfloor$, there exists a function $\phi_{m}$
implemented by a deep ReLU network with width $\mathcal{O}(dN)$ and depth
$\mathcal{O}(d^{2}L)$ such that

$\displaystyle\|\phi_{m}-g_{m}\|_{W^{1,\infty}((0,1)^{d})}\leq
50d^{5/2}(N+1)^{-4dL}.$

######  Lemma B.6.

Let $K\in\mathbb{N}$. For any $f\in W^{1,\infty}((0,1)^{d})$ with
$\|f\|_{W^{1,\infty}((0,1)^{d})}\leq 1$ and $m\in\\{1,2\\}^{d}$, there exists
a piecewise constant function $f_{K,m}$ on $\Omega_{m}$ satisfying

$\displaystyle\|f_{K,m}-f\|_{W^{1,\infty}(\Omega_{m})}\lesssim
1,\quad\|f_{K,m}-f\|_{L^{\infty}(\Omega_{m})}\lesssim 1/K$

with prefactors in $d$ omitted.

###### Proof.

We leverage approximation properties of averaged Taylor polynomials (Brenner
and Scott, 2008, Definition 4.1.3) and the Bramble-Hilbert Lemma (Brenner and
Scott, 2008, Lemma 4.3.8) to deduce local estimates and then combine them
through a partition of unity to obtain a global estimate. The key observation
is that the $L^{\infty}$ approximation bound can be established while
uniformly controlling the Lipschitz constant of the piecewise constant
function with a mild regularity assumption on the target function such as
$f\in W^{1,\infty}((0,1)^{d})$.

Without loss of generality, let us assume $m=m_{*}:=[1,1,\cdots,1]^{\top}$.
Following the proofs of Gühring, Kutyniok and Petersen (2020, Lemma C.4) and
Yang, Yang and Xiang (2023, Theorem 6), we first define an extension operator
$E:W^{1,\infty}((0,1)^{d})\to W^{1,\infty}({\mathbb{R}}^{d})$ to handle the
boundary. Accordingly, let $\tilde{f}:=Ef$ and $C_{E}$ be the norm of the
extension operator. Then for any $\Omega\subset{\mathbb{R}}^{d}$, it holds
that

$\displaystyle|\tilde{f}|_{W^{1,\infty}(\Omega)}\leq\|\tilde{f}\|_{W^{1,\infty}(\Omega)}\leq
C_{E}\|f\|_{W^{1,\infty}((0,1)^{d})}\leq C_{E}.$

Next, we define an average Taylor polynomial of order $1$ over
$B_{i,K}:=\mathbb{B}^{d}(\frac{8i+3}{8K},\frac{1}{4K},\|\cdot\|_{2})$ by

$\displaystyle
p_{f,i}(x):=\int_{B_{i,K}}T_{y}^{1}\tilde{f}(x)\phi_{K}(y)\mathrm{d}y$

where $\phi_{K}$ is a cut-off function supported on $\bar{B}_{i,K}$ as given
in Example 1. By Definition E.4, $T_{y}^{1}\tilde{f}(x)=\tilde{f}(y)$. Then,
it implies that $p_{f,i}(x)$ is a constant function as

$\displaystyle p_{f,i}(x)=\int_{B_{i,K}}\tilde{f}(y)\phi_{K}(y)\mathrm{d}y.$

Step 1. Get local estimates. For any
$i=[i_{1},i_{2},\cdots,i_{d}]^{\top}\in\\{0,1,\cdots,K\\}^{d}$, we would like
to employ the Bramble-Hilbert lemma E.8 on the subset

$\displaystyle\Omega_{m_{*},i}=\bar{\mathbb{B}}^{d}\left(\frac{8i+3}{8K},\frac{3}{8K},\|\cdot\|_{\infty}\right)=\prod_{j=1}^{d}\left[\frac{i_{j}}{K},\frac{3+4i_{j}}{4K}\right].$

It is easy to check the conditions of the Bramble-Hilbert lemma are fulfilled
as

$\displaystyle\frac{1}{4K}\geq\frac{1}{2}\times\frac{3}{8K}=\frac{1}{2}r_{\max}(\Omega_{m_{*},i}),\quad\gamma(\Omega_{m_{*},i})=\frac{d_{\Omega_{m_{*},i}}}{r_{\max}(\Omega_{m_{*},i})}=2\sqrt{d}.$

Hence, by the Bramble-Hilbert Lemma E.8, it yields that

$\displaystyle\|\tilde{f}-p_{f,i}\|_{L^{\infty}(\Omega_{m_{*},i})}\leq
C_{1}(d)|\tilde{f}|_{W^{1,\infty}(\Omega_{m_{*},i})}/K,$
$\displaystyle|\tilde{f}-p_{f,i}|_{W^{1,\infty}(\Omega_{m_{*},i})}\leq
C_{1}(d)|\tilde{f}|_{W^{1,\infty}(\Omega_{m_{*},i})}.$

Combining $|\tilde{f}|_{W^{1,\infty}(\Omega_{m_{*},i})}\leq C_{E}$ and the
inequalities above, it implies that

$\displaystyle\|\tilde{f}-p_{f,i}\|_{L^{\infty}(\Omega_{m_{*},i})}\leq
C_{1}(d)C_{E}/K,$ (B.1)
$\displaystyle\|\tilde{f}-p_{f,i}\|_{W^{1,\infty}(\Omega_{m_{*},i})}\leq
C_{1}(d)C_{E}.$ (B.2)

Step 2. Define a partition of unity. We construct a partition of unity in
order to combine the local estimates. Let $K\in\mathbb{N}$. For any $0\leq
i\leq K$, we define $h_{i}:{\mathbb{R}}\to{\mathbb{R}}$ by

$\displaystyle h_{i}(x)=h\left(4K\left(x-\frac{8i+3}{8K}\right)\right)\text{
where }h(x)=\begin{cases}1,&|x|<3/2,\\\ 0,&|x|>2,\\\ 4-2|x|,&3/2\leq|x|\leq
2.\end{cases}$

One can verify that $\\{h_{i}\\}_{i=1}^{K}$ is a partition of unity of $[0,1]$
and $h_{i}(x)=1$ for any $x\in\left[\frac{i}{K},\frac{3+4i}{4K}\right]$.
Considering the multidimensional case, for any
$x=[x_{1},x_{2},\cdots,x_{d}]^{\top}\in{\mathbb{R}}^{d}$ and any
$i=[i_{1},i_{2},\cdots,i_{d}]^{\top}\in\\{0,1,\cdots,K\\}^{d}$, let us define

$\displaystyle h_{i}(x):=\prod\nolimits_{j=1}^{d}h_{i_{j}}(x_{j}).$

Then a partition of unity of $[0,1]^{d}$ is defined by
$\\{h_{i}:i\in\\{0,1,\cdots,K\\}^{d}\\}$. Moreover, $h_{i}(x)=1$ for any
$x\in\Omega_{m_{*},i}=\prod_{j=1}^{d}\left[\frac{i_{j}}{K},\frac{3+4i_{j}}{4K}\right]$
and $i=[i_{1},i_{2},\cdots,i_{d}]^{\top}\in\\{0,1,\cdots,K\\}^{d}$. By the
definition of $h_{i}(x)$ on $\Omega_{m_{*},i}$ and (B.1), (B.2) , it yields
that

$\displaystyle\|h_{i}(\tilde{f}-p_{f,i})\|_{L^{\infty}(\Omega_{m_{*},i})}\leq\|\tilde{f}-p_{f,i}\|_{L^{\infty}(\Omega_{m_{*},i})}\leq
C_{1}(d)C_{E}/K,$
$\displaystyle\|h_{i}(\tilde{f}-p_{f,i})\|_{W^{1,\infty}(\Omega_{m_{*},i})}\leq\|\tilde{f}-p_{f,i}\|_{W^{1,\infty}(\Omega_{m_{*},i})}\leq
C_{1}(d)C_{E}.$

Step 3. Get global estimates. To deduce the global estimates, we start with
defining $f_{K,m_{*}}$ over $\Omega_{m_{*}}$ by

$\displaystyle
f_{K,m_{*}}:=\sum\nolimits_{i\in\\{0,1,\cdots,K\\}^{d}}h_{i}p_{f,i}.$

The error bounds follow that

$\displaystyle\|f_{K,m_{*}}-f\|_{L^{\infty}(\Omega_{m_{*}})}\leq\max_{i\in\\{0,1,\cdots,K\\}^{d}}\|h_{i}(\tilde{f}-p_{f,i})\|_{L^{\infty}(\Omega_{m_{*},i})}\leq
C_{1}(d)C_{E}/K,$
$\displaystyle\|f_{K,m_{*}}-f\|_{W^{1,\infty}(\Omega_{m_{*}})}\leq\max_{i\in\\{0,1,\cdots,K\\}^{d}}\|h_{i}(\tilde{f}-p_{f,i})\|_{W^{1,\infty}(\Omega_{m_{*},i})}\leq
C_{1}(d)C_{E}.$

This completes the proof. ∎

######  Lemma B.7.

Given any $f\in W^{1,\infty}((0,1)^{d})$ with
$\|f\|_{W^{1,\infty}((0,1)^{d})}\leq 1$, for any $N,L\in\mathbb{N}$ and any
$m\in\\{1,2\\}^{d}$, there exists a deep ReLU network $\psi_{m}$ with width
$\mathcal{O}(N\log N)$ and depth $\mathcal{O}(L\log L)$ such that

$\displaystyle\|\psi_{m}-f\|_{W^{1,\infty}(\Omega_{m})}\lesssim
1,\quad\|\psi_{m}-f\|_{L^{\infty}(\Omega_{m})}\lesssim(NL)^{-2/d},$

where we omit constants in $d$.

###### Proof.

The idea of proof is similar to those of Hon and Yang (2022, Theorem 3.1) and
Yang, Yang and Xiang (2023, Theorem 7). For completeness, we provide a
concrete proof in the following. Without loss of generality, we consider
$m=m_{*}:=[1,1,\cdots,1]^{\top}$. Given $K=\lfloor N^{1/d}\rfloor^{2}\lfloor
L^{2/d}\rfloor$, by Lemma B.6, we have

$\displaystyle\|f_{K,m_{*}}-f\|_{W^{1,\infty}(\Omega_{m_{*}})}\lesssim 1,$
$\displaystyle\|f_{K,m_{*}}-f\|_{L^{\infty}(\Omega_{m_{*}})}\lesssim
1/K\lesssim(NL)^{-2/d},$

where $f_{K,m_{*}}$ is a constant function for
$x\in\prod_{j=1}^{d}\left[\frac{i_{j}}{K},\frac{3+4i_{j}}{4K}\right]$ and
$i=[i_{1},i_{2},\cdots,i_{d}]^{\top}\in\\{0,1,\cdots,K-1\\}^{d}$. The insight
is to approximate $f_{K,m_{*}}$ with deep ReLU networks. Let
$\delta=1/(4K)\leq 1/(3K)$ in Lemma F.9. Then by Lemma F.9, there exists a
deep ReLU network $\phi_{1}(x)$ with width $4N+5$ and depth $4L+4$ such that

$\displaystyle\phi_{1}(x)=k,\quad
x\in\left[\frac{k}{K},\frac{k+1}{K}-\frac{1}{4K}\right],\quad
k=0,1,\cdots,K-1.$

We further define

$\displaystyle\phi_{2}(x)=\left[\frac{\phi_{1}(x_{1})}{K},\frac{\phi_{1}(x_{2})}{K},\cdots,\frac{\phi_{1}(x_{d})}{K}\right]^{\top}.$

For each $p=0,1,\cdots,K^{d}-1$, there exists a bijection

$\displaystyle\eta(p)=[\eta_{1},\eta_{2},\cdots,\eta_{d}]^{\top}\in\\{0,1,\cdots,K-1\\}^{d}$

satisfying $\sum_{j=1}^{d}\eta_{j}K^{j-1}=p$. We also define

$\displaystyle\xi_{p}=\frac{f_{K,m_{*}}(\eta(p)/K)+C_{2}(d)}{2C_{2}(d)}\in[0,1],$

where $|f_{K,m_{*}}|<C_{2}(d):=1+C_{1}(d)C_{E}$. Then, due to Lemma F.10,
there exists a deep ReLU network $\tilde{\phi}$ with width
$16(N+1)\log_{2}(8N)$ and depth $(5L+2)\log_{2}(4L)$ such that
$|\tilde{\phi}(p)-\xi_{p}|\leq(NL)^{-2}$ for $p=0,1,\cdots,K^{d}-1$. Let us
define

$\displaystyle\phi(x):=2C_{2}(d)\tilde{\phi}\left(\sum\nolimits_{j=1}^{d}\eta_{j}K^{j}\right)-C_{2}(d).$

Then it is clear that

$\displaystyle|\phi(\eta(p)/N)-f_{K,m_{*}}(\eta(p)/N)|=2C_{2}(d)|\tilde{\phi}(x)-\xi_{p}|\leq
2C_{2}(d)(NL)^{-2}.$

Furthermore, let $\psi_{m_{*}}(x):=\phi\circ\phi_{2}(x)$ for any
$x\in\Omega_{m_{*}}$. Since $\psi_{m_{*}}-f_{K,m_{*}}$ is a step function
whose first-order weak derivative is $0$ over $\Omega_{m_{*}}$, then it
implies that

$\displaystyle\|\psi_{m_{*}}-f_{K,m_{*}}\|_{W^{1,\infty}(\Omega_{m_{*}})}=\|\psi_{m_{*}}-f_{K,m_{*}}\|_{L^{\infty}(\Omega_{m_{*}})}\leq
2C_{2}(d)(NL)^{-2}.$

By the triangle inequalities for $\|\cdot\|_{L^{\infty}(\Omega_{m_{*}})}$ and
$\|\cdot\|_{W^{1,\infty}(\Omega_{m_{*}})}$, it is easy to derive that

$\displaystyle\|\psi_{m_{*}}-f\|_{W^{1,\infty}(\Omega_{m_{*}})}\leq\|\psi_{m_{*}}-f_{K,m_{*}}\|_{W^{1,\infty}(\Omega_{m_{*}})}+\|f_{K,m_{*}}-f\|_{W^{1,\infty}(\Omega_{m_{*}})}\lesssim
1,$
$\displaystyle\|\psi_{m_{*}}-f\|_{L^{\infty}(\Omega_{m_{*}})}\leq\|\psi_{m_{*}}-f_{K,m_{*}}\|_{L^{\infty}(\Omega_{m_{*}})}+\|f_{K,m_{*}}-f\|_{L^{\infty}(\Omega_{m_{*}})}\lesssim(NL)^{-2/d}.$

Lastly, we calculate the width and depth of the deep ReLU network to implement
$\psi_{m_{*}}=\phi\circ\phi_{2}$. Because that $\phi$ has width
$\mathcal{O}(N\log N)$ and depth $\mathcal{O}(L\log L)$ and $\phi_{2}$ has
width $\mathcal{O}(N)$ and depth $\mathcal{O}(L)$, the deep ReLU network of
$\psi_{m_{*}}$ is constructed with width $\mathcal{O}(N\log N)$ and depth
$\mathcal{O}(L\log L)$. This completes the proof. ∎

###### Proof of Lemma B.1.

We proceed in a similar way as the proof of Yang, Yang and Xiang (2023,
Theorem 3). By Lemma B.5, there exists a sequence of deep ReLU networks
$\\{\phi_{m}\\}_{m\in\\{1,2\\}^{d}}$ such that for any $m\in\\{1,2\\}^{d}$,

$\displaystyle\|\phi_{m}-g_{m}\|_{W^{1,\infty}((0,1)^{d})}\leq
50d^{5/2}(N+1)^{-4dL}.$

Each $\phi_{m}$ is implemented by a deep ReLU network with width
$\mathcal{O}(dN)$ and depth $\mathcal{O}(d^{2}L)$. By Lemma B.7, there exists
a collection of deep ReLU networks $\\{\psi_{m}\\}_{m\in\\{1,2\\}^{d}}$ such
that for any $m\in\\{1,2\\}^{d}$,

$\displaystyle\|\psi_{m}-f\|_{W^{1,\infty}(\Omega_{m})}\lesssim
1,\quad\|\psi_{m}-f\|_{L^{\infty}(\Omega_{m})}\lesssim(NL)^{-2/d},$

where we omit constants in $d$. Each $\psi_{m}$ is implemented by a deep ReLU
network with width $\mathcal{O}(N\log N)$ and depth $\mathcal{O}(L\log L)$.
Before proceeding, it is useful to estimate
$\|\phi_{m}\|_{L^{\infty}(\Omega_{m})}$,
$\|\phi_{m}\|_{W^{1,\infty}(\Omega_{m})}$,$\|\psi_{m}\|_{L^{\infty}(\Omega_{m})}$,