A Likelihood Based Approach to Distribution Regression Using Conditional Deep Generative Models

## Abstract

Abstract In this work, we explore the theoretical properties of conditional deep generative models under the statistical framework of distribution regression where the response variable lies in a high-dimensional ambient space but concentrates around a potentially lower-dimensional manifold. More specifically, we study the large-sample properties of a likelihood-based approach for estimating these models. Our results lead to the convergence rate of a sieve maximum likelihood estimator (MLE) for estimating the conditional distribution (and its devolved counterpart) of the response given predictors in the Hellinger (Wasserstein) metric. Our rates depend solely on the intrinsic dimension and smoothness of the true conditional distribution. These findings provide an explanation of why conditional deep generative models can circumvent the curse of dimensionality from the perspective of statistical foundations and demonstrate that they can learn a broader class of nearly singular conditional distributions. Our analysis also emphasizes the importance of introducing a small noise perturbation to the data when they are supported sufficiently close to a manifold. Finally, in our numerical studies, we demonstrate the effective implementation of the proposed approach using both synthetic and real-world datasets, which also provide complementary validation to our theoretical findings.

## 1 Introduction

Conditional distribution estimation provides a principled framework for characterizing the dependence relationship between a response variable $Y$ and predictors $X$, with the primary goal of estimating the distribution of
$Y$ conditional on $X$ through learning the (conditional) data-generating process. Conditional distribution estimation allows one to regress the entire distribution of
$Y$ on $X$, which provides much richer information than the traditional mean regression and plays a central role in various important areas ranging from causal inference (Pearl,, ; Spirtes,, ), graphical models (Jordan,, ; Koller and Friedman,, ), representation learning (Bengio et al.,, ), dimension reduction (Carreira-Perpinán,, ; Van Der Maaten et al.,, ), to model selection (Claeskens and Hjort,, ; Ando,, ). Their applications span across diverse domains such as forecasting (Gneiting and Katzfuss,, ), biology (Krishnaswamy et al.,, ), energy (Jeon and Taylor,, ), astronomy (Zhao et al.,, ), and industrial engineering (Simar and Wilson,, ), among others.

There is a rich literature in statistics and machine learning on conditional distribution estimation including both frequentist and Bayesian methods (Hall and Yao,, ; Norets and Pati,, ). Traditional methods, however, suffer from the curse of dimensionality and often struggle to adapt to the intricacies of modern data types such as the ones with lower-dimensional manifold structures.

Recent methodologies that leverage deep generative models have demonstrated significant advancements in complex data generation. Instead of explicitly modeling the data distribution, these approaches implicitly estimate it through learning the corresponding data sampling scheme. Commonly, these implicit distribution estimation approaches can be broadly categorized into three types. The first one is likelihood-based with notable examples including Kingma and Welling, (), Rezende et al., (), and Burda et al., (). The second approach, based on adversarial learning, matches the empirical distribution of the data with a distribution estimator using an adversarial loss. Representative examples include Goodfellow et al., (), Arjovsky et al., (), and Mroueh et al., (), among others. The third approach, which is more recent, reduces the problem of distribution estimation to score estimation through certain time-discrete or continuous dynamical systems.
The idea of score matching was first proposed in Hyvärinen and Dayan, () and Vincent, (). More recently, score-based diffusion models
have achieved state-of-the-art performance in many applications (Sohl-Dickstein et al.,, ; Nichol and Dhariwal,, ; Song et al.,, ).

On the theoretical front, recent works such as Liu et al., (), Chae et al., (), Altekrüger et al., (), and Tang and Yang, () demonstrate that distribution estimation based on deep generative models can adapt to the intrinsic geometry of the data, with convergence rates dependent on the intrinsic dimension of the data, thus potentially circumventing the curse of dimensionality.
Such advancement has naturally motivated us to employ and investigate *conditional deep generative model* for conditional distribution estimation.
Specifically, we explore and study the theoretical properties of a new likelihood-based approach to conditional sampling using deep generative models for data potentially residing on a low-dimensional manifold corrupted by full-dimensional noise.
More concretely, we consider the following *conditional distributional regression* problem:

$$ $Y|X=V|X+\varepsilon,$ (1) $$

where $X$ serves as a predictor in ${\mathbb{R}}^{\mathfrak{p}}$, $V|X$ represents the (uncorrupted) underlying response supported on a manifold of dimension ${\mathfrak{d}}\leq D$, $Y|X$ represents the observed response, and $\varepsilon\sim{\mathsf{N}}(0,\,\sigma_{*}^{2}I_{D})$ denotes the noise residing in the ambient space ${\mathbb{R}}^{D}$. Our deep generative model focuses on the conditional distribution $V|X$ by using a (conditional) generator of the form $G_{*}(Z,X)$, where $G_{*}$ is a function of a random seed $Z$ and the covariate information $X$. This approach is termed ‘conditional deep generative’ because the conditional generator is modeled using deep neural networks (DNNs). Observe that, when ${\mathfrak{d}}<D$, the distribution of $G_{*}(Z,X)$ is supported on a lower-dimensional manifold, making it singular with respect to the Lebesgue measure in the $D$-dimensional ambient space. We study the statistical convergence rate of sieve MLEs in the conditional deep general model setup and investigate its dependence on the intrinsic dimension, structure properties of the model as well as the noise level of the data.

### 1.1 List of contributions

We briefly summarise the main contributions made in this paper.

- •
To the best of our knowledge, our study is the first attempt to explore the likelihood-based approach for distributional regression using a conditional deep generative model, considering full-dimensional noise and the potential presence of singular underlying support.
We provide a solid statistical foundation for the approach by proving the near-optimal convergence rates for this proposed estimator.
- •
We derive the convergence rates for the conditional density estimator of the corrupted data $Y$ with respect to the Hellinger distance and specialize the obtained rate for two popular deep neural network classes: the sparse and fully connected network classes.
Furthermore, we characterize the Wasserstein convergence rates for the induced intrinsic conditional distribution estimator
on the manifold (i.e., a deconvolution problem). Both rates turn out to depend only on the intrinsic dimension and smoothness of the true conditional distribution.
- •
Our analysis in Corollary 2 suggests the need to inject a small amount of noise into the data when they are sufficiently close to the manifold. Intuitively, this observation validates the underlying structural challenges in related manifold estimation problems with noisy data, as outlined by Genovese et al., ().
- •
We show that the class of learnable (conditional) distributions of our method is broad.
It encompasses not only the smooth distributions class, but also extends to the general (nearly) singular distributions with manifold structures, with minimal assumptions.

### 1.2 Other relevant literature

The problem of non-parametric conditional density estimation has been extensively explored in statistical literature.
Hall and Yao, (), Bott and Kohler, (), and Bilodeau et al., () directly tackle this problem with smoothing and local polynomial-based methods. Fan and Yim, () and Efromovich, () explore suitably transformed regression problems to address this challenge. Other notable approaches include the nearest neighbor method (Izbicki et al.,, ; Bhattacharya and Gangopadhyay,, ), basis function expansion (Sugiyama et al.,, ; Izbicki and Lee,, ), and tree-based boosting (Pospisil and Lee,, ; Gao and Hastie,, ), among others.

In the context of conditional generation, we highlight recent work by Zhou et al., () and Liu et al., (). In Zhou et al., (), GANs were employed to investigate conditional density estimation. While this work offers a consistent estimator, it lacks statistical rates or convergence analysis, and its focus is on a low-dimensional setup. In Liu et al., (), conditional density estimation supported on a manifold using Wasserstein-GANs was examined. However, their setup does not account for smoothness across either covariates or responses, nor do they address how deep generative models specifically tackle the challenges of high-dimensionality.
Moreover, their assumption that the data lies exactly on the manifold can be restrictive.
Our study shares some commonalities with the work of Chae et al., (), as both investigate sieve maximum likelihood estimators (MLEs). However, the fundamental problems addressed and the methodologies employed differ significantly, and our work involves technical challenges that span multiple scales. While Chae et al., () concentrates exclusively on unconditional distribution estimation, our theoretical analysis necessitates much more nuanced techniques due to the conditional nature of our setup. This shift is noteworthy because it demands a more refined analysis of entropy bounds, considering two potential sources of smoothness - across the regressor and the response variables. Furthermore, our setting accommodates the possibility of an infinite number of $x$ values, which gives rise to a dynamic manifold structure, further compounding the intricacy of the problem at hand.

## 2 Conditional deep generative models for distribution regression

We consider the following probabilistic conditional generative model, where for a given predictor value $x$, the response $Y$ is generated by

$$ $Y=\,G_{*}(Z,x)+\varepsilon,\quad x\in{\mathcal{X}}\subset\mathbb{R}^{\mathfrak {p}}.$ (2) $$

Here, $G_{*}(\cdot,x):{\mathcal{Z}}\to{\mathcal{M}}_{x}$ is the unknown generator function, $Z$ a latent variable with a known distribution $P_{Z}$ and support $\mathcal{Z}\subset{\mathbb{R}}^{\mathfrak{d}}$ independent of the predictor $X$. The existence of the generator $G_{*}$ directly follows from Noise Outsourcing Lemma 3. This lemma enables the transfer of randomness into the covariate and an orthogonal (independent) component through a generating function for any regression response. We denote ${\mathcal{M}}:\,=\cup_{x\in{\mathcal{X}}}{\mathcal{M}}_{x}\subset\mathbb{R}^{D}$ as the support of the image of $G_{*}({\mathcal{Z}},{\mathcal{X}})$ such as a (union of) $d$-dimensional manifold. We model $G_{*}(\cdot,\cdot):{\mathcal{Z}}\times{\mathcal{X}}\subset{\mathbb{R}}^{
\mathfrak{d}}\times{\mathbb{R}}^{\mathfrak{p}}\rightarrow\mathcal{Y}\subset{
\mathbb{R}}^{D}$ using a deep neural network, leading to a *conditional deep generative model* for (2).

In the next section, we present a more general result in terms of the entropy bound (variance) for the true function class of $G_{*}$ and the approximability (bias) of the search class. We then proceed to a simplified understanding in the context of conditional deep generative models in subsequent sections.

### 2.1 Convergence rates of the Sieve MLE

In light of equation (2), it is evident that the distribution of $Y|X=x$ results from the convolution of two distinct distributions: the pushforward of $Z$ through $G_{*}$ with $X=x$, and $\varepsilon$ following an independent $D$-dimensional normal distribution. The density corresponding to the true distribution $P_{*}(\cdot|X=x)$ can thus be expressed as:

$$ $p_{*}(y|x)=\int\phi_{\sigma_{*}}(y-G_{*}(z,x))\,dP_{Z},$ $$

where $\phi_{\sigma_{*}}$ is the density of $\mathsf{N}({0},\sigma^{2}_{*}I_{d})$.
We define the class of conditional distributions ${\mathcal{P}}$ as

$$ ${\mathcal{P}}=\Big{\{}P_{g,\sigma}:g(\cdot,x)\in\mathcal{F},\sigma\in[\sigma_{ \min},\sigma_{\max}]\Big{\}},$ (3) $$

where $P_{g,\sigma}$ represents the distribution with density $p_{g,\sigma}=\int\phi_{\sigma}(y-g(z,x))dP_{Z}$. In this notation, $P_{*}=P_{G_{*},\sigma_{*}}$ and $p_{*}=p_{G_{*},\sigma_{*}}$. The elements of ${\mathcal{P}}$ comprise two components: $g$ originating from the underlying function class $\mathcal{F}$, and $\sigma$, which characterizes the noise component.
This class enables us to obtain separate estimates for $G_{*}$ and $\sigma_{*}$, furnishing us with both the canonical estimator for the distribution of $Y|X=x$ and enhancing our comprehension of the singular distribution of $G_{*}(Z,x)$, supported on a low-dimensional manifold.

Given a data set $\{(X_{i},Y_{i})\}_{i=1}^{n}$, the log-likelihood function is defined as $\ell_{n}(g,\sigma)=n^{-1}\sum_{i=1}^{n}\log p_{g,\sigma}(Y_{i}|X_{i})$. For a sequence $\eta_{n}\downarrow 0$ as $n\to\infty$, a *sieve* maximum likelihood estimator (MLE) (Geman and Hwang,, ) is any estimator $(\widehat{g},\widehat{\sigma})\in\mathcal{F}\times[\sigma_{\min},\sigma_{\max}]$ that satisfies

$$ $\ell_{n}\left(\widehat{g},\widehat{\sigma}\right)\geq\sup_{\overset{\sigma\in[ \sigma_{\min},\sigma_{\max}]}{g\in\mathcal{F}}}\ell_{n}(g,\sigma)-\eta_{n}.$ (4) $$

Here $\widehat{g}\in\mathcal{F}$ and $\widehat{\sigma}\in[\sigma_{\min},\sigma_{\max}]$ are the estimators, and $\eta_{n}$ represents the optimization error. The dependence of $\widehat{g}$ and $\widehat{\sigma}$ on $n$ illustrates the sieve’s role in approximating the true distribution when optimization is performed over the class ${\mathcal{P}}$. The estimated density $\widehat{p}=p_{\widehat{g},\widehat{\sigma}}$ provides an estimator for $p_{*}(\cdot|\cdot)$, and $Q_{\widehat{g}}(\cdot|X=x)$ serve as the estimator for $Q_{*}(\cdot|X=x)$.

In this section, we formulate the main results, which provide convergence rates in the Hellinger distance for our sieve MLE estimator. The convergence rate was derived for any search functional class $\mathcal{F}$, with a brief emphasis on their entropy and approximation capabilities.

###### Assumption 1 (True distribution) .

Denote $\mu_{X}^{\ast}(x)$ as the distribution of $X$.
We denote the true conditional densities as $p_{*}=\{p_{*}(\cdot|x),x\in\mathbb{R}^{\mathfrak{p}}\}$. It is natural to assume that
the data is generated from $p_{*}$ from model (2) with some true generator $G_{*}$ and $\sigma_{*}$.
We denote $Q_{*}(\cdot|X=x)$ (or $Q_{G_{*}}$) as the distribution of $G_{*}(Z,x)$ for some distribution $P_{Z}$.

A function $g$ is said to have a composite structure (Schmidt-Hieber,, ; Kohler and Langer,, ) if it takes the form as

$$ $g=f_{q}\circ f_{q-1}\circ\cdots\circ f_{1}$ (5) $$

where $f_{j}:(a_{j},b_{j})^{d_{j}}\to(a_{j+1},b_{j+1})^{d_{j+1}}$, $d_{0}={\mathfrak{p}}+{\mathfrak{d}}$ and $d_{q+1}=D$. Denote $f_{j}=(f_{j}^{(1)},\ldots,f_{j}^{(d_{j+1})})$ as the components of $f_{j}$, let $t_{j}$ be the maximal number of variables on which each of the $f^{(i)}_{j}$ depends and let $f^{(i)}_{j}\in\mathcal{H}^{\beta_{j}}\left((a_{j},b_{j})^{t_{j}},K\right)$ (see Section 2.4.1 for the definition of the Hölder class $\mathcal{H}^{\beta}$). A composite structure is very general which includes smooth functions and additive structure as special cases. In addition, in the next section, we show the class of conditional distributions $\left\{Q_{G_{*}}(\cdot|X=x):x\in\mathbb{R}^{\mathfrak{p}},G_{*}\in\mathcal{G}\right\}$ induced by the composite structure is broad.

###### Assumption 2 (composite structure ) .

Denote $\mathcal{G}=\mathcal{G}\left(q,\bm{d},\bm{t},\bm{\beta},K\right)$ as a collection of functions of form (5), where $\bm{d}=(d_{0},\ldots,d_{q+1})$, $\bm{t}=(t_{0},\ldots,t_{q+1})$, and $\bm{\beta}=(\beta_{0},\ldots,\beta_{q+1})$. We regard $(q,\bm{d},\bm{t},\bm{\beta},K)$ as constants in our setup, and assume that the true generator $G_{*}(\cdot,x)$ as in (2) belongs to $\mathcal{G}$, for all $x\in{\mathcal{X}}$. Additionally, we assume $\||G_{*}|_{\infty}\|_{\infty}\leq K$.

$$ $\widetilde{\beta}_{j}=\beta_{j}\prod_{l=j+1}^{q}\left(\beta_{l}\wedge 1\right) ,\qquad j_{*}=\operatorname*{argmax}_{j\in\{0,\ldots,q\}}\frac{t_{j}}{ \widetilde{\beta}_{j}},\qquad\beta_{*}=\widetilde{\beta}_{j_{*}},\qquad t_{*}= t_{j_{*}}.$ $$

The quantities $t_{*}$ and $\beta_{*}$ are called intrinsic dimension and smoothness of $G_{*}$ (or of $\mathcal{G}$).

###### Remark 1 (Strength of the Composition Structure) .

The expression $(a_{j},b_{j})\subset[-K,K]$ can be intuitively visualized by setting $a_{j}=-K$ and $b_{j}=K$. To illustrate the impact of intrinsic dimensionality and smoothness, consider a function $f:\mathbb{R}^{d}\to\mathbb{R}$ defined as $f(x)=f_{1}(x_{1})+\ldots+f_{d}(x_{d})$, where $x=(x_{1},\ldots,x_{d})$ and $f_{j}\in\mathcal{H}^{\beta}((-K,K),K)$ for $j=1,\ldots,d$. While $f\in\mathcal{H}^{\beta}((-K,K)^{d},K)$, its intrinsic dimension is $t_{*}=1$ with intrinsic smoothness $\beta$. This mitigates the curse of dimensionality.

###### Assumption 3 .

Let ${\mathcal{M}}_{*}$ be the closure of $G_{*}({\mathcal{Z}},{\mathcal{X}})$. We assume that ${\mathcal{M}}_{*}$ does not have an interior point, and reach$({\mathcal{M}}_{*})=\mathsf{r}_{\ast}$ with $\mathsf{r}_{\ast}>0$.

2 permits low intrinsic dimensionality within the learnable function class. 3 imposes the strong identifiability condition necessary for efficient estimation, as seen in manifold literature (Aamari and Levrard,, ; Tang and Yang,, ).

Given two conditional densities $p_{1}(\cdot|x),p_{2}(\cdot|x)$ and $\mu_{X}^{\ast}$ denoting the density of $X$, we use integrated distances for a measure of evaluation. With a slight abuse of notation, we denote $d_{1}(p_{1},p_{2})=\mathbb{E}_{X}\left[d_{1}(p_{1}(\cdot|x),p_{2}(\cdot|x))\right]$ and $d_{H}(p_{1},p_{2})=\mathbb{E}_{X}\left[d_{H}(p_{1}(\cdot|x),p_{2}(\cdot|x)\right]$, where $d_{1}$ and $d_{H}$ represent the $L_{1}$ and the Hellinger distance as $d_{1}(p_{1}(\cdot|x),p_{2}(\cdot|x))=\int\left|p_{1}(y|x)-p_{2}(y|x)\right|dy$ and $d_{H}(p_{1},p_{2})=(\int\int[\sqrt{p_{1}(y|x)}-\sqrt{p_{2}(y|x)}]^{2}\,dy)^{1/2}$ respectively.
Denote ${\mathcal{N}}(\delta,\mathcal{F},d)$ and ${\mathcal{N}}_{[]}(\delta,\mathcal{F},d)$ as covering and bracketing numbers of the function class $\mathcal{F}$ with respect to the (pseudo)-metric $d$.

We first present Lemma 1, which establishes the bracketing entropy of the functional class $\mathcal{P}$ with respect to Hellinger
distance in terms of the covering entropy of the search class $\mathcal{F}$.
This enables us to transfer the entropy control of the individual components $\mathcal{F}$ and $\sigma$ to the entire $\mathcal{P}$.

###### Lemma 1 .

Let $\mathcal{F}$ be class of functions from ${\mathcal{Z}}\times{\mathcal{X}}$ to $\mathbb{R}^{D}$ such that $\||g|_{\infty}\|_{\infty}\leq K$ for every $g\in\mathcal{F}$. Let ${\mathcal{P}}=\left\{P_{g,\sigma}:g\in\mathcal{F},\sigma\in[\sigma_{\min},
\sigma_{\max}]\right\}$ with $\sigma_{\min}\leq 1$. Then, there exist constants $c=c(\sigma_{\max},K,D)$ and $C=C(\sigma_{\max},K,D)$ and $\delta_{*}=\delta_{*}(D)$ such that for every $\delta\in(0,\delta_{*}]$,

$$ $\log{\mathcal{N}}_{[]}(\delta,{\mathcal{P}},d_{H})\leq\log{\mathcal{N}}(c \sigma_{\min}^{D+3}\delta^{4},\mathcal{F},\||\cdot|_{\infty}\|_{\infty})+\log \left(\frac{C}{\sigma_{\min}^{D+2}\delta^{4}}\right),$ (6) $$

The proof of Lemma 1 is provided in the Appendix D.
Theorem 1 presents the convergence rate of the sieve-MLE to the true distribution (see Appendix E for the proof).

###### Theorem 1 .

Let $\mathcal{F},\mathcal{P},\sigma_{\min}$ and $\delta_{*}=\delta_{*}(D)$ be given as in Lemma 1, and $n\geq 1$. Suppose that $\log{\mathcal{N}}\left(\delta,\mathcal{F},\||\cdot|_{\infty}\|_{\infty}\right)
\leq\xi\left\{A+1\vee\log\delta^{-1}\right\}$ for every $\delta\in(0,\delta_{*}]$ and some $A,\xi>0$. Suppose that there exists a $G\in\mathcal{F}$ and some $\delta_{\operatorname*{approx}}\in(0,\delta_{*}]$ such that $\||G-G_{*}|_{\infty}\|_{\infty}\leq\delta_{\operatorname*{approx}}$. Furthermore, suppose that $s\geq 1$, $A\geq 1$, $\sigma_{min}\leq 1$, $\delta_{\operatorname*{approx}}\leq 1$ and $\sigma_{*}\in[\sigma_{\min},\sigma_{\max}]$. Then

$$ $P_{*}\left(d_{H}(\widehat{p},p_{*})>\varepsilon_{n}^{*}\right)\leq 5e^{-C_{1}n \varepsilon_{n}^{*2}}+C_{2}n^{-1}$ (7) $$

provided that $\eta_{n}\leq n\varepsilon^{*2}_{n}/6$ and $\varepsilon^{*}_{n}\leq\sqrt{2}\delta_{*}$, where

$$ $\varepsilon_{n}^{*}=C_{3}\left(\sqrt{\frac{\xi\left\{A+\log\left(n/\sigma_{ \min}\right)\right\}}{n}}\vee\frac{\delta_{\operatorname*{approx}}}{\sigma_{*} }\right),$ (8) $$

$C_{1}$ is an absolute constant, $C_{2}=C_{2}(D)$ and $C_{3}=C_{3}(D,K,\sigma_{\max})$.

The outlined rate has two components: the statistical component, expressed as an upper bound to the metric entropy of $\mathcal{F}$, and the approximation component, denoted as $\delta_{\operatorname*{approx}}$. The statistical error is quantified by measuring the complexity of the class $\mathcal{P}$, as formulated in Lemma 1. The approximation error is assessed through the ability of the provided function class to approximate the true distribution.

###### Assumption 1 (True distribution) .

###### Assumption 2 (composite structure ) .

###### Remark 1 (Strength of the Composition Structure) .

###### Assumption 3 .

###### Lemma 1 .

###### Theorem 1 .

### 2.2 Neural network class

We model $G_{*}(\cdot,\cdot)$ using a deep neural network.
More specifically,
we parameterize the true generator $G_{*}$ with a deep neural neural architecture $(L,{\bf r})$ of the form

$$ $f:\mathbb{R}^{r_{0}}\to\mathbb{R}^{r_{L+1}},\qquad z\mapsto f(z)=W_{L}\rho_{v_ {L}}W_{L-1}\rho_{v_{L-}}\ldots W_{1}\rho_{v_{1}}W_{0}z,$ (9) $$

where $W_{j}\in\mathbb{R}^{r_{j+1}\times r_{j}},v_{j}\in\mathbb{R}^{r_{j}}$, $\rho_{v_{j}}(\cdot)=\mathrm{ReLU}(\cdot-v_{j})$ and ${\bf r}=(r_{0},\ldots,r_{L+1})\in{\mathbb{N}}^{L+2}$. The constant $L$ is the number of hidden layers and $r=(r_{0},\ldots,r_{L+1})$ represents the number of nodes in each layer.

We define the sparse neural architecture class $\mathcal{F}_{s}(L,{\bf r},s,B,K)$ as set of functions of form (9) satisfying

$$ $\max_{0\leq j\leq L}|W_{j}|_{\infty}\vee|v_{j}|_{\infty}\leq B,\qquad\sum_{j=1 }^{L}|W_{j}|_{0}+|v_{j}|_{0}\leq s,\qquad\||f|_{\infty}\|_{\infty}\leq K,$ $$

with $r_{0}={\mathfrak{d}}+{\mathfrak{p}}$ and $r_{L+1}=D$, where $|\cdot|_{0}$ and $|\cdot|_{\infty}$ stand for the $L^{0}$ and $L^{\infty}$ vector norms, and $\||f|_{\infty}\|_{\infty}=\sup_{x\in{\mathbb{R}}^{r_{0}}}\max_{i=1,\ldots,D}|f
_{i}(x)|$, $s$ is sparsity parameter and $K$ is functional bound.

The fully connected neural architecture class $\mathcal{F}_{c}=\mathcal{F}_{c}\left(L,{\bf r},B,K\right)$ is set of functions of form (9) satisfying

$$ $\max_{0\leq j\leq L}|W_{j}|_{\infty}\vee|v_{j}|_{\infty}\leq B,\qquad\||f|_{ \infty}\|_{\infty}\leq K.$ $$

Both classes $\mathcal{F}_{s}$ and $\mathcal{F}_{c}$ for the deep generator will be considered in our analysis of the resulting sieve maximum likelihood estimator. We denote the corresponding sieve-MLE as $\widehat{p}_{s}$ and $\widehat{p}_{c}$, respectively. When we use $r$ instead of ${\bf r}$, it refers to $r_{1}=\ldots=r_{L}=r$ along with $r_{0}={\mathfrak{d}}+{\mathfrak{p}}$ and $r_{L+1}=D$.

We can simplify and visualize the result stated in Theorem 1 in both cases: when the sieve-MLE is obtained with optimization performed over the class $\mathcal{F}_{s}$ and $\mathcal{F}_{c}$.
To fulfill the conditions stated in the Theorem 1, we need to establish entropy bounds for these function classes, $\mathcal{F}_{s}$ and $\mathcal{F}_{c}$, and gain insight into their approximation capabilities for the composite structure class described in 2.

For the sparse neural architecture class $\mathcal{F}_{s}(L,r,s,K)$, the entropy, formally stated as Proposition 1 in Ohn and Kim, (), is bounded as follows.

$$ $\log{\mathcal{N}}(\delta,\mathcal{F}_{s},\||\cdot|_{\infty}\|_{\infty}) \lesssim sL\,\{\log(BLr)+\log\delta^{-1}\}.$ (10) $$

From an entropy perspective, the fully connected neural architecture class $\mathcal{F}_{c}(L,r,B,K)$ can be viewed as $\mathcal{F}_{s}$ without any sparsity constraint, meaning $s\asymp r^{2}L$. Therefore, we have

$$ $\log{\mathcal{N}}(\delta,\mathcal{F}_{c},\||\cdot|_{\infty}\|_{\infty}) \lesssim L^{2}r^{2}\{\log(BLr)+\log\delta^{-1}\}.$ (11) $$

The approximation properties of the sparse and fully connected network are provided in Lemma 4.1 and Lemma 4.2 of the Appendix J, respectively.

Having established the essential components for $\mathcal{F}_{c}$ in (11) and Lemma 4.2, and for $\mathcal{F}_{s}$ in (10) and Lemma 4.1, respectively, we can simplify Theorem 1 and state Corollary 1.

###### Corollary 1 .

Suppose that Assumptions 1 and 2 hold, and $\sigma_{*}\in\left[\sigma_{\min},\sigma_{\max}\right]$ with $\sigma_{\min}\leq 1$ and $\sigma_{\max}<\infty$. Moreover, assume that the noise $\sigma_{*}$ decays at rate $\alpha$, i.e., $\sigma_{*}\asymp n^{-\alpha}$, and $\sigma_{\min}=n^{-\gamma}$ for some $\gamma\geq\alpha\geq 0$.
Then, for every $\delta_{\operatorname*{approx}}\in[0,1]$, the following holds:

- 1.
Let $\mathcal{F}_{s}=\mathcal{F}_{s}\left(L,r,s,B,K\right)$ with $\delta_{*}=\delta_{*}(D)$ be as given in Lemma 1, and $L\asymp\log\delta^{-1}_{\operatorname*{approx}}$, $r\asymp\delta^{-t_{*}/\beta_{*}}_{\operatorname*{approx}}$, $s\asymp\delta^{-t_{*}/\beta_{*}}_{\operatorname*{approx}}\log\delta^{-1}_{
\operatorname*{approx}}$, $B\asymp\delta^{-1}_{\operatorname*{approx}}$. Then the sieve MLE $\widehat{p}_{s}$ satisfies (7) with $\varepsilon_{n}^{*}$ as in (8) with $\xi=\delta^{-t_{*}/\beta_{*}}_{\operatorname*{approx}}\log^{2}(\delta^{-1}_{
\operatorname*{approx}})$ and $A=\log^{2}(\delta^{-1}_{\operatorname*{approx}})$ provided that $\eta_{n}\leq n\varepsilon^{*2}_{n}/6$ and $\varepsilon^{*}_{n}\leq\sqrt{2}\delta_{*}$.
- 2.
Let $\mathcal{F}_{c}=\mathcal{F}_{c}\left(L,r,B,K\right)$ with $\delta_{*}=\delta_{*}(D)$ be as given in Lemma 1, and $L\asymp\log\delta^{-1}_{\operatorname*{approx}}$, $r\asymp\delta^{-t_{*}/2\beta_{*}}_{\operatorname*{approx}}$, $B\asymp\delta^{-1}_{\operatorname*{approx}}$. Then the sieve MLE $\widehat{p}_{c}$ satisfies (7) with $\varepsilon_{n}^{*}$ as in (8) with $\xi=\delta^{-t_{*}/\beta_{*}}_{\operatorname*{approx}}\log^{2}(\delta^{-1}_{
\operatorname*{approx}})$ and $A=\log^{2}(\delta^{-1}_{\operatorname*{approx}})$ provided that $\eta_{n}\leq n\varepsilon^{*2}_{n}/6$ and $\varepsilon^{*}_{n}\leq\sqrt{2}\delta_{*}$.

In particular, choosing
$\delta_{\operatorname*{approx}}:=\left(\sigma_{*}^{2}/n\right)^{\beta_{*}/(2
\beta_{*}+t_{*})}$
minimizes $\varepsilon_{n}^{*}\asymp\sqrt{\xi\left\{A+\log\left(n/\sigma_{\min}\right)
\right\}/n}\vee\delta_{\operatorname*{approx}}/\sigma_{*}$, and gives

$$ $\varepsilon_{n}^{*}\asymp n^{-\frac{\beta_{*}-t_{*}\alpha}{2\beta_{*}+t_{*}}} \log^{2}(n).$ (12) $$

###### Remark 2 .

The convergence rate in (12) illustrates the influence of intrinsic dimensionality, smoothness, and noise level on the estimation process. Note that $\alpha$ is upper bounded as $\varepsilon_{n}^{*}\leq\sqrt{2}\delta_{*}(D)$. For large values of $\alpha$, estimation of $G_{*}$ is inherent difficult as the data is very close on the singular support. To address this, a small noise injection, as described in Corollary 2, can smooth the estimation and ensure consistency.

The proof of Corollary 1 is provided in Appendix F. For the composite structural class $\mathcal{G}$, the effective smoothness is denoted by $\beta_{*}$, and the dimension is $t_{*}$. This effectively mitigates the curse of dimensionality. The convergence rate at (12) also recovers the optimal rate when $q=1$ and $\alpha=0$, and there is a small lag of polynomial factor $t_{*}\alpha/(2\beta_{*}+t_{*})$ when $\alpha>0$ (Norets and Pati,, ). This lag arises due to the presence of full-dimensional noise in the response observation $Y$.
Note that when the noise is small, that is $\alpha$ is large, achieving a sharp estimation of $p_{*}$ requires an equally accurate estimate of $G_{*}$. This can be quite challenging. Our practically tractable approach attempts to address this without initially estimating the singular support.

###### Corollary 1 .

###### Remark 2 .

### 2.3 Wasserstein convergence of the intrinsic (conditional) distributions

Using Wasserstein distance as a metric for distributions $Q_{g}$ is meaningful due to their singularity in ambient space: when ${\mathfrak{d}}<D$, the conditional distribution is singular with respect to the Lebesgue measure on ${\mathbb{R}}^{D}$.

The integrated Wasserstein distance, for $r\geq 1$, between $P_{1}(\cdot|X)$ and $P_{2}(\cdot|X)$ is defined as

$$ $W_{r}\left(P_{1},P_{2}\right)=\mathbb{E}_{X}\left[\inf_{\beta\in\Gamma(P_{1},P _{2})}\left(\mathbb{E}_{(U_{1},U_{2})\sim\beta}\big{[}|U_{1}-U_{2}|_{r}^{r} \big{]}\right)^{1/r}\right],$ $$

where $\Gamma(P_{1},P_{2})$ is the set of all couplings between $P_{1}$ and $P_{2}$ that preserves the two marginals. The (dual) representation of this norm, $W_{r}(P_{1},P_{2})=\mathbb{E}_{X}\left[\sup_{\|f\|_{{\mathrm{L}ip}_{r}}\leq 1}
\left\{\,\mathbb{E}_{P_{1}}[f]-\mathbb{E}_{P_{2}}[f]\,\right\}\right]$ (Villani et al.,, ) with $\|\cdot\|_{{\mathrm{L}ip}_{r}}$ denoting the $r$-Lipschitz norm, is particularly useful in our proofs.

###### Theorem 2 .

Suppose that Assumption 3 holds. If $d_{H}(p_{g,\sigma},p_{*})\leq\varepsilon$ holds for some $\varepsilon\in[0,1]$ and some $p_{g,\sigma}\in{\mathcal{P}}$, then we have

$$ $W_{1}(Q_{g},Q_{*})\leq C\left(\varepsilon+\sigma_{*}\sqrt{\log\varepsilon^{-1} }\right),$ $$

where $C=C(D,K,\mathsf{r}_{\ast})$ depends only on $(D,K,\mathsf{r}_{\ast})$.

The proof of Theorem 2 is provided in Appendix G. Theorem 2 guarantees that $W_{1}\left(\widehat{Q}_{\widehat{g}},Q_{*}\right)\lesssim_{\log}d_{H}(\widehat
{p},p_{*})+\sigma_{*}$, where $\lesssim_{\log}$ represents less than or equal up to a logarithmic factor of $n$. Following from Corollary 1, the Wasserstein convergence rate, $n^{-(\beta_{*}-t_{*}\alpha)/(2\beta_{*}+t_{*})}\log^{2}(n)\vee\sigma_{*}\log^{
1/2}(n)$, comprises two components: the convergence rate in the Hellinger distance and the standard deviation of the true noise sequence. It is noteworthy that the first expression is influenced by the variance of noise by the factor $\alpha$. When $\alpha$ is very small, indicating that the data $Y_{j}$ lies very close to the manifold, the second expression $n^{-\alpha}$ in the overall rate dominates. Intuitively, this phenomenon arises from the underlying structural challenges in related manifold estimation problems with noisy data, as discussed by Genovese et al., (). To address this issue, we propose a data perturbation strategy by transforming the data $\{(Y_{j},X_{j})\}_{j=1}^{n}$ into $\{(\widetilde{Y}_{j},X_{j})\}_{j=1}^{n}$, where $\widetilde{Y}_{j}=Y_{j}+{\bm{\epsilon}}_{j}$ and ${\bm{\epsilon}}_{j}\sim{\mathsf{N}}\left(0_{D},n^{-\beta_{*}/(\beta_{*}+t_{*})
}\,I_{D}\right)$. The resulting estimation error bound is summarized below, whose proof is provided in Appendix H.

###### Corollary 2 .

Suppose that Assumption 1, 2, and 3 hold, and $\sigma_{*}\in[\sigma_{\min},\sigma_{\max}]$ with $\sigma_{*}=n^{-\alpha}$ and $\sigma_{\min}=n^{-\gamma}$ for some $0\leq\alpha\leq\gamma$. Then for each of the network architecture classes (sparse and fully connected) with the network parameters specified in Corollary 1, the sieve MLE $\widehat{p}_{per}$ and $\widehat{Q}_{per}$ based on the perturbed data $\{(\widetilde{Y}_{j},X_{j})\}_{j=1}^{n}$ satisfies

$$ $P_{*}\left[W_{1}\left(\widehat{Q}_{per},Q_{*}\right)\geq\left(\varepsilon^{*}_ {n}+\sigma_{*}\sqrt{\log((\varepsilon^{*}_{n})^{-1})}\right)\right]\lesssim 5e ^{-C_{1}n{\varepsilon^{*}_{n}}^{2}}+\frac{C_{2}}{n}$ $$

where $\varepsilon^{*}_{n}$ can be chosen such that

$$ $\varepsilon^{*}_{n}+\sigma_{*}\sqrt{\log((\varepsilon^{*}_{n})^{-1})}\,\asymp \begin{cases}n^{-\frac{\beta_{*}-t_{*}\alpha}{2\beta_{*}+t_{*}}}\,\log^{2}(n), \qquad&\text{if}\;\alpha<\beta_{*}/\{2(\beta_{*}+t_{*})\},\\ n^{-\frac{\beta_{*}}{2(\beta_{*}+t_{*})}}\,\log^{2}(n),\qquad&\text{otherwise} .\end{cases}$ (13) $$

###### Theorem 2 .

###### Corollary 2 .

### 2.4 Characterization of the learnable distribution class

Section 2.2 focuses on the true generator $G_{*}$ within the class of functions with composition structures.
In this subsection, we show that such a conditional distribution class achieved by the push-forward map
$G_{*}$ is broad and includes many existing distribution classes for $Q_{*}$ as special cases.

#### 2.4.1 Smooth conditional density

For $\beta>0$, let $\mathcal{H}^{\beta}(D,M)$ be the class of all $\beta$-Hölder functions $f:D\subset{\mathbb{R}}^{\mathfrak{d}}\to{\mathbb{R}}$ with $\beta$-Hölder norm bounded by $M>0$. Let $\mathcal{H}^{\beta}(D)=\cup_{M>0}\mathcal{H}^{\beta}(D,M)$. See Appendix A for their formal definitions.

###### Lemma 2 .

Suppose that (i) ${\mathcal{Z}}\times{\mathcal{X}}$ and ${\mathcal{Y}}$ are uniformly convex and (ii) $p_{Z}\in\mathcal{H}^{\beta_{Z}}({\mathcal{Z}})$, $\mu^{*}_{X}\in\mathcal{H}^{\beta_{X}}({\mathcal{X}})$ and $q_{*}\in\mathcal{H}^{\beta_{Q}}({\mathcal{Y}})$ for some $\beta_{Z},\beta_{X},\beta_{Q}>0$ and are bounded above and below.
Then, there exists a map $g(\cdot,\cdot):{\mathcal{Z}}\times{\mathcal{X}}\to{\mathcal{Y}}$ such that $Q_{*}(\cdot|\cdot)=Q_{g}$ and $g\in\mathcal{H}^{\beta_{\min}+1}({\mathcal{Z}}\times{\mathcal{X}})$,
where $\beta_{\min}=\min\{\beta_{Z},\beta_{X},\beta_{Q}\}$.

Lemma 2 establishes that the learnable distribution class includes Hölder-smooth functions with smoothness parameter $\beta_{\min}$ and intrinsic dimension ${\mathfrak{d}}$. As a result, following Corollary 1, the convergence rate for density estimation is given by $\varepsilon^{*}_{n}\asymp n^{-(\beta_{\min}+1-{\mathfrak{d}}\alpha)/(2\beta_{
\min}+2+{\mathfrak{d}})}$.
A push-forward map is a transport map between two distributions.
The well-established regularity theory of transport map in optimal transport is directly applicable here [see Villani et al., () and Villani, ()].
The proof of Lemma 2 is based on Theorem 12.50 of Villani et al., () and Caffarelli, (), which establishes the regularity of this transport map and its existence follows from Brenier, ().
When $p_{Z}$ is selected as a well-behaved parametric distribution, the regularity of the transport map is determined by the smoothness of both $\mu_{X}^{*}$ and $Q_{*}$. For a more detailed discussion on this, please refer to Appendix B.

###### Lemma 2 .

#### 2.4.2 A broader conditional distribution class with smoothness disparity

In Appendix K, we present a novel approximation result for the function class exhibiting smoothness disparity in Theorem 5. This new result facilitates the study of theoretical properties of estimators when the generator $G_{*}\in\mathcal{H}_{{\mathfrak{d}},{\mathfrak{p}}}^{\beta_{Z},\beta_{X}}\left
({\mathcal{Z}},{\mathcal{X}},K\right)$.
Note that such a function class defined in (16) in Appendix K is much broader compared to the smoothness class in Section 2.4.1 as $Z$ and $X$ do not have to be jointly smooth and it allows for smoothness disparity among them. The subsequent Theorem 3 combines our approximation result with (11) and enables us to specialize Theorem 1 to this class (see Appendix I for the proof).

###### Theorem 3 .

Let $G_{*}\in\mathcal{H}_{{\mathfrak{d}},{\mathfrak{p}}}^{\beta_{Z},\beta_{X}}\left
({\mathcal{Z}},{\mathcal{X}},K\right)$. Suppose that Assumption 1 holds and $\sigma_{*}\in\left[\sigma_{\min},\sigma_{\max}\right]$ with $\sigma_{\min}\leq 1$ and $\sigma_{\max}<\infty$. Moreover, we assume $\sigma_{*}\asymp n^{-\alpha}$, and $\sigma_{\min}=n^{-\gamma}$ for some $0\leq\alpha\leq\gamma\leq(\beta_{Z}^{-1}{\mathfrak{d}}+\beta_{X}^{-1}{
\mathfrak{p}})^{-1}$. Then, for every $\delta_{\operatorname*{approx}}\in[0,1]$, we have: Let $\mathcal{F}_{s}=\mathcal{F}_{s}\left(L,r,s,1,K\right)$ with $L\asymp\log\delta^{-1}_{\operatorname*{approx}}$, $r\asymp\delta^{-(\beta_{Z}^{-1}{\mathfrak{d}}+\beta_{X}^{-1}{\mathfrak{p}})}_{
\operatorname*{approx}}$, $s\asymp\delta^{-(\beta_{Z}^{-1}{\mathfrak{d}}+\beta_{X}^{-1}{\mathfrak{p}})}_{
\operatorname*{approx}}\log\delta^{-1}_{\operatorname*{approx}}$. Then the sieve MLE $\widehat{p}_{s}$ satisfies (7) with the rate outlined in (8) with $\xi=\delta^{-(\beta_{Z}^{-1}{\mathfrak{d}}+\beta_{X}^{-1}{\mathfrak{p}})}_{
\operatorname*{approx}}\log^{2}\delta^{-1}_{\operatorname*{approx}}$ and $A=\log^{2}\delta^{-1}_{\operatorname*{approx}}$, provided that $\eta_{n}\leq n\varepsilon^{*2}_{n}/6$.
In particular, choosing
$\delta_{\operatorname*{approx}}:=\left(\sigma_{*}^{2}/n\right)^{1/\left(2+
\beta_{Z}^{-1}{\mathfrak{d}}+\beta_{X}^{-1}{\mathfrak{p}}\right)}\leq 1$
minimizes $\varepsilon_{n}^{*}\asymp\sqrt{\xi\left\{A+\log\left(n/\sigma_{\min}\right)
\right\}/n}\vee\delta_{\operatorname*{approx}}/\sigma_{*}$, and gives

$$ $\varepsilon_{n}^{*}\asymp n^{-\frac{1-\alpha(\beta_{Z}^{-1}{\mathfrak{d}}+ \beta_{X}^{-1}{\mathfrak{p}})}{2+\beta_{Z}^{-1}{\mathfrak{d}}+\beta_{X}^{-1}{ \mathfrak{p}}}}\log^{2}(n).$ (14) $$

The proof of Theorem 3 is provided in Appendix I.
In the special case when $\alpha=0$ and ${\mathfrak{d}}=D$, our convergence rate in (14) recovers the minimax optimal rate for conditional density estimation based on kernel smoothing, as established in Li et al., ().

###### Theorem 3 .

#### 2.4.3 Conditional distribution on manifolds

In this part, we extend Lemma 2 and provide the existence of the generator when the conditional distribution is supported on a compact manifold with dimension $\mathsf{d}_{\ast}\leq D$. Due to space constraints, we provide only a sketched proof here; the detailed proof can be found in Appendix C. Specifically, we first present arguments for the existence of the generator when ${\mathcal{Y}}$ is covered by a single chart. We then extend this to the multiple chart case using the technique of partition of unity.

In the simpler case when there exists a single $({\mathcal{Y}},\varphi)$ covering ${\mathcal{Y}}$, where $\varphi:{\mathcal{B}}_{1}({0}_{\mathsf{d}_{\ast}})\to{\mathcal{Y}}$ is a homeomorphism, we assume $\varphi\in\mathcal{H}^{\beta_{\min}+1}$. In this case, we use the change of variable formula to transfer the measure on ${\mathcal{B}}_{1}({0}_{\mathsf{d}_{\ast}})$ (unit ball in ${\mathbb{R}}^{\mathsf{d}_{\ast}}$) from ${\mathcal{Y}}$. Following Lemma 2, we can find a transport map $g\in\mathcal{H}^{\beta_{\min}}$ mapping from ${\mathcal{Z}}\times{\mathcal{X}}$ to ${\mathcal{B}}_{1}({0}_{\mathsf{d}_{\ast}})$. The map $g\circ\varphi$ then serves as our generator.

In the general case where the compact manifold $\mathcal{Y}$ needs to be covered by multiple charts, demonstrating the existence of a transport or push-forward map is challenging because $\mathcal{Y}$ is not uniformly convex. Suppose that $\{(U_{k},\varphi_{k})\}_{k=1}^{K}$ forms a cover of ${\mathcal{Y}}$. Due to the compactness of ${\mathcal{Y}}$, the number of charts $K$ is finite. Analogous to the single chart scenario, we first construct $g_{k}\circ\varphi_{k}$ to transport the measure on each chart. We then patch these local transport maps together to construct a global transport map; see Appendix C for full details. As a result, following Corollary 1, the convergence rate for density estimation shall be given by $\varepsilon^{*}_{n}\asymp n^{-(\beta_{\min}-{\mathfrak{d}}\alpha)/(2\beta_{
\min}+{\mathfrak{d}})}$.

## 3 Numerical Results

In this section, we present numerical experiments to validate and complement our theoretical findings using two synthetic dataset examples. These experiments cover a range of scenarios, including full-dimensional cases as well as benchmark examples involving manifold-based data. Additionally, we provide a real data example to further enrich our experimentation and validation process.
It is worth noting that, although not significant, the computational cost of fitting a conditional generative model is higher compared to fitting an unconditional one, as the input dimension of the deep neural network (DNN) is ${\mathfrak{p}}+{\mathfrak{d}}$ rather than just ${\mathfrak{d}}$.

##### Learning algorithm to compute sieve MLE.

For the computational algorithm, we adopt a common conditional variational auto-encoder (VAE) architecture to maximize the following log-likelihood term:$\sum_{j=1}^{n}{\mathcal{L}}_{\mathrm{V}AE}(g,\sigma,\phi;Y_{j},X_{j})$, where

$$ ${\mathcal{L}}_{\mathrm{V}AE}(g,\sigma,\phi;y,x)=\log\left(\frac{p_{g,\sigma}(y ,x,z)}{q_{\phi}(Z|y,x)}\right).$ $$

The variational distribution $q_{\phi}(Z|y,x)$ is chosen as the standard normal family ${\mathsf{N}}(\mu_{\phi}(y,x),\Sigma_{\phi}(y,x))$.

We examine two classes of datasets: (i) full-dimensional response and (ii) response residing on a low-dimensional manifold. The first highlights the generality of our proposed approach, while the second underscores its efficiency in terms of the Wasserstein metric and validates the small noise perturbation strategy outlined in Corollary 2.

Simulation from full dimension distribution. We use the following models for data generation.

- •
FD1 : $Y=\mathbb{I}_{\{U<0.5\}}\,{\mathsf{N}}\left(-X,0.25^{2}\right)+\mathbb{I}_{\{U
>0.5\}}\,{\mathsf{N}}\left(X,0.25^{2}\right)$, where $U\sim\mathrm{Unif}(0,1)$ and $X_{1}\sim{\mathsf{N}}(3,1)$.
- •
FD2 : $Y=X_{1}^{2}+e^{(X_{2}+X_{3}/3)}+\sin(X_{4}+X_{5})+\varepsilon$, where $\{X_{j}\}_{j=1}^{5}\stackrel{{\scriptstyle i.i.d}}{{\sim}}{\mathsf{N}}(0,1)$
- •
FD3 : $Y=X_{1}^{2}+e^{(X_{2}+X_{3}/3)}+X_{4}-X_{5}+0.5\,(1+X_{2}^{2}+X_{5}^{2})\times\varepsilon$, where $\{X_{j}\}_{j=1}^{5}\stackrel{{\scriptstyle i.i.d}}{{\sim}}{\mathsf{N}}(0,1)$.

These are examples of a mixture model, an additive noise model, and a multiplicative noise model, respectively. The neural architecture for both the encoder and decoder consists of two deep layers, i.e., $L=2$. The hyperparameters are as follows: $r_{\mathrm{enc}}=({\mathfrak{p}}+1,10,10)$ for $\mu_{\phi}$ and $\Sigma_{\phi}$, and $r_{\mathrm{dec}}=(10+{\mathfrak{p}},10,1)$ for $g$. The sample size used for simulation is 5000, with a training-to-testing ratio of $4:1$. We employ a batch size of $64$ with a learning rate of $10^{-3}$.

We compare the sieve MLE with CKDE (Hall et al.,, ) and FlexCode proposed by Izbicki and Lee, (). To evaluate their performance, we compute the mean squared error (MSE) for both the mean and the standard deviation.
We use Monte Carlo approximation to compute the mean and standard deviation for the sieve MLE, and numerical integration for CKDE and Flexcode. This evaluation strategy resembles that implemented by Zhou et al., (). Table 1 summarizes the findings.

**Table 1: MSE for the estimated conditional mean and the standard deviation.**
|  |  | Sieve MLE | CKDE | FlexCode |
| --- | --- | --- | --- | --- |
| FD1 | MEAN | $\mathbf{0.0379}\pm 0.0170$ | $1.0053\pm 0.1004$ | $1.1660\pm 0.1076$ |
| SD | $\mathbf{0.0280}\pm 0.0045$ | $0.9887\pm 0.0347$ | $1.2000\pm 0.0126$ |  |
| FD2 | MEAN | $\mathbf{0.1943}\pm 0.0427$ | $0.2640\pm 0.0515$ | $0.3954\pm 0.0571$ |
| SD | $\mathbf{0.2843}\pm 0.0093$ | $0.2853\pm 0.0213$ | $5.8278\pm 0.1607$ |  |
| FD3 | MEAN | $\mathbf{0.2337}\pm 0.0453$ | $0.2967\pm 0.0537$ | $1.3419\pm 0.1087$ |
| SD | $1.6394\pm 0.0861$ | $\mathbf{0.6334}\pm 0.0460$ | $11.4898\pm 0.1559$ |  |

Note that the sieve MLE outperforms all other methods in all scenarios except for the MSE(SD) for the FD3 dataset. However, for the FD3 dataset, we found that as the training sample size increases further, the MSE(SD) of the sieve MLE achieves performance increasingly comparable to CKDE.

Simulation from distributions on manifolds.
We consider two examples of manifolds with an intrinsic dimension ${\mathfrak{d}}=1$, while the ambient dimension is $D=2$.

- •
M1 : $Y=G_{*}(Z,U)+\varepsilon$, with $G_{*}=(G_{*}^{(1)},G_{*}^{(2)})$, $G_{*}^{(1)}=\mathbb{I}_{\{U<0.5\}}\left(1-\cos(Z)\right)+\mathbb{I}_{\{U>0.5\}
}\cos(Z)$, $G_{*}^{(2)}=\mathbb{I}_{\{U<0.5\}}\left(0.5-\sin(Z)\right)+\mathbb{I}_{\{U>0.5
\}}\sin(Z)$; where $Z\sim\mathrm{Unif}(0,\pi)$ and $U\sim\mathrm{Unif}(0,1)$.
- •
M2 : $Y=G_{*}(Z,U)+\varepsilon$, with $G_{*}=\left(G_{*}^{(1)},G_{*}^{(2)}\right)$, $G_{*}^{(1)}=\mathbb{I}_{\{U<0.5\}}\cos(Z)+\mathbb{I}_{\{U>0.5\}}\,2\cos(Z)$, $G_{*}^{(2)}=\mathbb{I}_{\{U<0.5\}}\,0.5\sin(Z)+\mathbb{I}_{\{U>0.5\}}\sin(Z)$; where $Z\sim\mathrm{Unif}(0,2\pi)$ and $U\sim\mathrm{Unif}(0,1)$.

The manifold $M_{1}$ consists of two moons. The manifold $M_{2}$ comprises ellipses, with conditions distinguishing the inner and outer confocal ellipses. The noise sequence follows a two-dimensional centered Gaussian distribution, $\varepsilon\sim{\mathsf{N}}(0_{2},\sigma_{*}^{2}I_{2})$. We investigated this setup across various noise variances $\sigma_{*}^{2}$. Our neural architecture employed $r_{\mathrm{enc}}=({\mathfrak{p}}+2,100,100,2)$ for $\mu_{\phi}$ and $\Sigma_{\phi}$, and $r_{\mathrm{dec}}=(2+{\mathfrak{p}},100,100,2)$ for $g$. We utilized a sample size of 5000 for simulation, with a training-to-testing ratio of $4:1$. A batch size of $100$ was employed, with a learning rate of $10^{-3}$.

Figure: Figure 1: Generated samples from manifold $M_{1}$ and $M_{2}$ are displayed in the left panel. The right panel shows box plots for the empirical Wasserstein distance at different noise levels $\sigma_{*}$.
Refer to caption: extracted/5897099/images/twomoon_sample_0.01.png

We computed the empirical $W_{1}$ distance using the algorithm proposed by Cuturi, () to evaluate the performance. The right panel of Figure 1 presents the boxplots of $W_{1}$ between the true and learned distribution for $M_{1}$ and $M_{2}$ across 20 repetitions. The left panel highlights the following general behaviors:

- •
When $\alpha$ is small and close to zero, the noise variance is large, making estimation challenging due to the singularity of the true data distribution.
- •
When $\alpha$ is large, the noise variance is very small, and the perturbed data facilitates efficient estimation.

This observed pattern, as emphasized in Corollary 2, closely aligns with the results achieved in (13).

##### Numerical result for real data

We utilized the widely used MNIST dataset for two purposes: to demonstrate the generalizability of our approach to a benchmark image dataset where the intrinsic dimension ${\mathfrak{d}}$ is much lesser than the ambient dimension $D=784$ and to underscore the effectiveness of sparse networks as outlined in Lemma 4.1 and Corollary 1.1.

For the fully connected architecture, we set $r_{\mathrm{enc}}=(10+784,512,2)$ for $\mu_{\phi}$ and $\Sigma_{\phi}$, and $r_{\mathrm{dec}}=(10+2,512,784)$ for $g$. For the sparse architecture, we use $r_{\mathrm{enc}}=(10+784,608,432,256,2)$ for $\mu_{\phi}$ and $\Sigma_{\phi}$, and $r_{\mathrm{dec}}=(10+2,256,432,608,784)$ for $g$. The input dimension of $10$ for both the encoder and decoder corresponds to the one-hot encoding of the labels. We employ a batch size of $64$ with a learning rate of $10^{-3}$.

Figure 2 presents a visual comparison between real and generated images, organized according to their respective labels. The real images were randomly sampled from the training set along with their corresponding labels, while the generated images were produced using these labels (conditions) and random seeds.

Figure: Figure 2: MNIST images: real images (left panel), generated images with sparse architecture (central panel), and generated images with fully connected architecture (right panel)
Refer to caption: extracted/5897099/images/MNIST_original.png

This MNIST example highlights a case where the intrinsic dimension is significantly smaller than the ambient data dimension. This example serves to validate the proposed methodology in high-dimensional settings.

## 4 Discussion

We investigated statistical properties of a likelihood-based conditional deep generative model for distribution regression in a scenario where the response variable is situated in a high-dimensional ambient space but is centered around a potentially lower-dimensional intrinsic structure. Our analysis established favorable rates in both the Hellinger and Wasserstein metrics which are dependent on only the intrinsic dimension of the data. Our theoretical findings show that the conditional deep generative models can circumvent the curse of dimensionality for high-dimensional distribution regression. To the best of our knowledge, our work is the first of its kind.

Given the novelty of emerging statistical methodologies with intricate structural considerations in the study of deep generative models, there exist numerous paths for future exploration. Among these potential directions, we are particularly interested in investigating controllable generation via penalized optimization methods, studying statistical properties of deep generative models trained via matching flows, as well as delving into the hypothesis testing problem within the framework of deep generative models, among others. Another interesting direction is to explore residual neural network structure for modeling time series of distributions with interesting temporal dependence structures.

## Appendix A Notation

We denote $a\vee b$ and $a\wedge b$ as the maximum and minimum of two real numbers $a$ and $b$, respectively. The notation $\lceil a\rceil$ represents the smallest integer greater than or equal to $a$. The inequality $a\lesssim b$ indicates that $a$ is less than or equal to $b$ up to a multiplicative constant. When we write $a\lesssim_{\log}b$, it means that $a$ is less than or equal to $b$ up to a logarithmic factor, specifically $\log(n)$. We denote $a\asymp b$ when both $a\lesssim b$ and $b\lesssim a$ hold. For vector norms, $|\cdot|_{{\mathfrak{p}}}$ represents the $\ell^{{\mathfrak{p}}}$ norm, while $\|\cdot\|_{{\mathfrak{p}}}$ denotes the $L^{{\mathfrak{p}}}$-norm of a function for $1\leq{\mathfrak{p}}\leq\infty$. Lastly, ${\mathcal{B}}_{\epsilon}(u)$ signifies the Euclidean open ball with radius $\epsilon$ centered at $u$.

We use the multi-index notation through the main paper and the appendix. Denote $\mathbb{N}$ as the set of natural numbers and $\mathbb{N}_{0}$ as $\mathbb{N}\cup\{0\}$. For a vector $\mathbf{x}\in{\mathbb{R}}^{r}$, we denote the components as $\mathbf{x}=(x^{(1)},\ldots,x^{(r)})$. Given a function $f:D\subset{\mathbb{R}}^{r}\to{\mathbb{R}}$, the operator is defined as $\partial^{\bm{\alpha}}:=\partial^{\alpha^{(1)}}\ldots\partial^{\alpha^{(r)}}$ with $\bm{\alpha}\in\mathbb{N}_{0}^{r}$, where $\partial^{\alpha^{(j)}}f:=\partial^{\alpha^{(j)}}f(\mathbf{x})/\partial x^{(j)}$. For $\bm{\alpha}\in\mathbb{N}_{0}^{r}$, the expression $|\bm{\alpha}|=\sum_{j=1}^{r}|\alpha^{(j)}|$. Given a function $f(\cdot,\cdot):D\times D_{\prime}\subset{\mathbb{R}}^{r}\times{\mathbb{R}}^{r_
{\prime}}\to{\mathbb{R}}$, we denote the operator $\partial^{\bm{\alpha}+\bm{\alpha}_{\prime}}:=\partial^{\alpha^{(1)}}\ldots
\partial^{\alpha^{(r)}}\,\partial^{\alpha_{\prime}^{(1)}}\ldots\partial^{
\alpha_{\prime}^{(r_{\prime})}}$, with $\bm{\alpha}\in\mathbb{N}_{0}^{r}$ and $\bm{\alpha}_{\prime}\in\mathbb{N}_{0}^{r_{\prime}}$, where $\partial^{\alpha^{(j)}}f(\mathbf{x},\mathbf{y})=\partial^{\alpha^{(j)}}f(
\mathbf{x},\mathbf{y})/\partial^{\alpha^{(j)}}x^{(j)}$ and $\partial^{\alpha_{\prime}^{(j)}}f(\mathbf{x},\mathbf{y})=\partial^{\alpha_{
\prime}^{(j)}}f(\mathbf{x},\mathbf{y})/\partial y^{(j)}$, with $\mathbf{x}\in D$ and $\mathbf{y}\in D_{\prime}$. This notation allows us to represent the derivative with variable $\mathbf{x}$ and $\mathbf{y}$ separately through the vector $\bm{\alpha}$ and $\bm{\alpha}_{\prime}$ which is required to tackle the smoothness disparity along $x$ and $y$ variable. The $\beta-$Hölder class functions are defined as

$$ $\displaystyle\mathcal{H}_{{r}}^{\beta}$ $\displaystyle\left(D,M\right)=\Big{\{}f:D\subset{\mathbb{R}}^{r}\rightarrow \mathbb{R}:$ (15) $\displaystyle\sum_{\begin{subarray}{c}{\bm{\alpha}}:|{\bm{\alpha}}|<\beta\end{ subarray}}\|\partial^{{\bm{\alpha}}}f\|_{\infty}+\sum_{\begin{subarray}{c}{\bm {\alpha}}:|{\bm{\alpha}}|=\lfloor\beta\rfloor\end{subarray}}\,\sup_{\begin{ subarray}{c}\mathbf{u}_{1},\mathbf{u}_{2}\in D\\ \mathbf{u}_{1}\neq\mathbf{u}_{2}\end{subarray}}\frac{|\partial^{{\bm{\alpha}}} f(\mathbf{u}_{1})-\partial^{{\bm{\alpha}}}f(\mathbf{u}_{2})|}{|\mathbf{u}_{1}- \mathbf{u}_{2}|_{\infty}^{\beta-\lfloor\beta\rfloor}}\leq M\Big{\}},$ $$

We extend this definition to include the Hölder class of functions with differences in smoothness (smoothness disparity) along two variables. This class is defined as

$$ $\displaystyle\mathcal{H}_{{r},{r_{\prime}}}^{\beta,\beta_{\prime}}$ $\displaystyle\left(D,D_{\prime},M\right)=\Big{\{}f(\cdot,\cdot):D\times D_{ \prime}\subset\mathbb{R}^{{r}}\times\mathbb{R}^{{r_{\prime}}}\rightarrow \mathbb{R}:$ (16) $\displaystyle\sum_{\begin{subarray}{c}{\bm{\alpha}}:|{\bm{\alpha}}|<\beta\\ {\bm{\alpha}_{\prime}}:|{\bm{\alpha}_{\prime}}|<\beta_{\prime}\end{subarray}} \|\partial^{{\bm{\alpha}}+{\bm{\alpha}_{\prime}}}f\|_{\infty}+\sum_{\begin{ subarray}{c}{\bm{\alpha}}:|{\bm{\alpha}}|=\lfloor\beta\rfloor\\ {\bm{\alpha}_{\prime}}:|{\bm{\alpha}_{\prime}}|=\lfloor\beta_{\prime}\rfloor \end{subarray}}\,\sup_{\begin{subarray}{c}\mathbf{u}_{1},\mathbf{u}_{2}\in D_{ X}\\ \mathbf{v}_{1},\mathbf{v}_{2}\in D_{Y}\\ \mathbf{u}_{1}\neq\mathbf{u}_{2}\\ \mathbf{v}_{1}\neq\mathbf{v}_{2}\end{subarray}}\frac{|\partial^{{\bm{\alpha}}+ {\bm{\alpha}_{\prime}}}f(\mathbf{v}_{1},\mathbf{u}_{1})-\partial^{{\bm{\alpha} }+{\bm{\alpha}_{\prime}}}f(\mathbf{v}_{2},\mathbf{u}_{2})|}{|\mathbf{u}_{1}- \mathbf{u}_{2}|_{\infty}^{\beta-\lfloor\beta\rfloor}\vee|\mathbf{v}_{1}- \mathbf{v}_{2}|_{\infty}^{\beta_{\prime}-\lfloor\beta_{\prime}\rfloor}}\leq M \Big{\}}.$ $$

We denote $\mathcal{H}_{{r}}^{\beta}(D)=\cup_{M>0}\mathcal{H}_{{r}}^{\beta}\left(D,M\right)$ and $\mathcal{H}_{{r},{r_{\prime}}}^{\beta,\beta_{\prime}}\left(D,D_{\prime}\right)
=\cup_{M>0}\mathcal{H}_{{r},{r_{\prime}}}^{\beta,\beta_{\prime}}\left(D,D_{
\prime},M\right)$.

## Appendix B More on Smooth conditional density

###### Theorem 4 ( Villani et al., ( 2009 ) Theorem 12.50) .

Suppose that

- (i)
$\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ are uniformly convex, bounded, open subsets of ${\mathbb{R}}^{\mathfrak{d}}$ with $\mathcal{C}^{\lfloor\beta\rfloor+2}$ (continuously differentiable up to order ${\lfloor\beta\rfloor+2}$) boundaries,
- (ii)
$h_{1}\in\mathcal{H}^{\beta}(\mathcal{A}_{1})$ and $h_{2}\in\mathcal{H}^{\beta}(\mathcal{A}_{2})$ for some $\beta>0$, are probability densities bounded above and below.

Then, there exists a unique map (up to an additive constant) $g:\mathcal{A}_{1}\to\mathcal{A}_{2}$ with $g\in\mathcal{H}^{\beta+1}(\mathcal{A}_{1})$, such that if $U\sim h_{1}$ then $g(U)\sim h_{2}$.

###### Proof of Lemma 2 .

Given that $Z$ and $X$ is independent, the product measure on ${\mathcal{Z}}\times{\mathcal{X}}$ is $p_{Z}\mu_{X}^{*}$. Following the smoothness from $p_{Z}$ and $\mu_{X}^{*}$, the map $p_{Z}(\cdot)\mu_{X}^{*}(\cdot)\in\mathcal{H}^{\min\{\beta_{Z},\beta_{X}\}}({
\mathcal{Z}}\times{\mathcal{X}})$. This implies that $p_{Z}(\cdot)\mu_{X}^{*}(\cdot)\in\mathcal{H}^{\min\{\beta_{Z},\beta_{X},\beta_
{Q}\}}({\mathcal{Z}}\times{\mathcal{X}})$. Again $q_{*}\in\mathcal{H}^{\beta_{Q}}({\mathcal{Y}})$ implies $q_{*}\in\mathcal{H}^{\min\{\beta_{Z},\beta_{X},\beta_{Q}\}}({\mathcal{Y}})$. The result now follows directly from Theorem 4.
∎

Many of the problems in the conditional setting have an analog in the joint setup. Our proposed approach has a direct statistical extension to this setup. The sufficiency of such extension follows from the observation in the subsequent Lemma 3 which is based on Lemma 2.1 and Lemma 2.2 of Zhou et al., () (see also Theorem 5.10 of Kallenberg, ()).

###### Lemma 3 (Noise Outsourcing Lemma) .

Let $(Y,X)\in{\mathcal{Y}}\times{\mathcal{X}}$ with joint distribution $P_{X,Y}$. Suppose $Y$ is standard Borel space, then there exists $Z\sim{\mathsf{N}}(0,I_{m})$ for any given $m\geq 1$, independent of $X$ and Borel measurable function $G:\mathbb{R}^{m}\times{\mathcal{X}}\to{\mathcal{Y}}$ such that

$$ $\left(X,G(Z,X)\right)\sim(Y,X).$ (17) $$

Moreover, the condition (17) is equivalent of

$$ $G(Z,x)\sim P_{Y|X=x}.$ $$

###### Theorem 4 ( Villani et al., ( 2009 ) Theorem 12.50) .

###### Proof of Lemma 2 .

###### Lemma 3 (Noise Outsourcing Lemma) .

## Appendix C More on Conditional distribution on manifolds

Suppose $({\mathcal{Y}},\varphi)$ is the single chart covering ${\mathcal{Y}}$, where $\varphi:{\mathcal{B}}_{1}({0}_{\mathsf{d}_{\ast}})\to{\mathcal{Y}}$ is a homeomorphism. We assume that $\varphi\in\mathcal{H}^{\beta_{\min}+1}$, and that $\inf_{\mathbf{u}\in{\mathcal{B}}_{1}({0}_{\mathsf{d}_{\ast}})}|J_{\varphi}(
\mathbf{u})|$ is bounded below by a positive constant, where

$$ $|J_{\varphi}(\mathbf{u})|=\sqrt{\mathrm{det}\left(\frac{\partial\varphi}{ \partial\mathbf{u}^{\top}}\frac{\partial\varphi}{\partial\mathbf{u}}\right)}$ $$

is the Jacobian determinant of $\varphi$.

Note that when $\mathsf{d}_{\ast}<D$, the distribution $Q_{*}$ cannot possess a Lebesgue density because of the singularity of ${\mathcal{Y}}$. We, therefore consider a density with respect to the $\mathsf{d}_{\ast}-$dimensional Hausdorff measure in ${\mathbb{R}}^{D}$, denoted by $\mathsf{H}_{\mathsf{d}_{\ast}}$. Suppose that $Q$ allows the Radon-Nikodym derivative $q$ with respect to $\mathsf{H}_{\mathsf{d}_{\ast}}$. We further assume that $q$ is bounded from above and below and that $q\circ\varphi\in\mathcal{H}^{\beta_{\min}}$. Then by change of variable formula, the Lebesgue density of $\widetilde{Q}$, the push-forward measure on ${\mathcal{B}}_{1}({0}_{\mathsf{d}_{\ast}})$ through the map $\varphi^{-1}$, is given as

$$ $\widetilde{q}(\mathbf{u})=q(\varphi(\mathbf{u}))|J_{\varphi}(\mathbf{u})|.$ $$

Following the assumptions on the Jacobian determinant and $\varphi\in\mathcal{H}^{\beta_{\min}+1}$, it follows that $|J_{\varphi}(\mathbf{u})|$ is bounded from above and below, and the map $\mathbf{u}\mapsto|J_{\varphi}(\mathbf{u})|$ belongs to $\mathcal{H}^{\beta_{\min}}$. Therefore, $\widetilde{q}$ is bounded above and below, belongs to $\mathcal{H}^{\beta_{\min}}({\mathcal{B}}_{1}({0}_{\mathsf{d}_{\ast}}))$. By Lemma 2, assuming $\beta_{\min}\leq\beta_{Z}\wedge\beta_{X}$, there exists $g\in\mathcal{H}^{\beta_{{\min}}+1}$ such that $\widetilde{Q}=Q_{g}$. Thus, we have $Q=Q_{g\circ\varphi}$, where $g\circ\varphi:{\mathcal{Z}}\times{\mathcal{X}}\to{\mathcal{Y}}$. Following Lemma 4, it is possible to find the appropriate neural network approximating them.

Suppose ${\mathcal{Y}}$ is covered by the charts $\left\{\left(U_{k},\varphi_{k}\right)\right\}_{k=1}^{K}$, with $1<K<\infty$, where $\varphi_{k}:{\mathcal{B}}_{1}(0_{\mathsf{d}_{\ast}})\to U_{k}$ is a homeomorphism. As before, we assume $\varphi_{k}\in\mathcal{H}^{\beta_{\min}+1}$, $|J_{\varphi_{k}}(\mathbf{u})|$ is bounded below by a positive constant, $Q$ possesses density $q$ with respect to $\mathsf{H}_{\mathsf{d}_{\ast}}$ that is bounded above and below, and that $q\circ\phi_{k}\in\mathcal{H}^{\beta_{\min}}$. Let $Q_{k}(\cdot)=Q(\cdot)/Q(U_{k})$ be the normalized measure of $Q$ over $U_{k}$.

We denote $q_{k}$ as the corresponding density with respect to $\mathsf{H}_{\mathsf{d}_{\ast}}$. For $\mathbf{u}\in U_{k}\cap U_{\ell}$, $q_{k}(\mathbf{u})Q(U_{k})=q_{\ell}(\mathbf{u})Q(U_{\ell})=q(\mathbf{u})$ holds due to the measure $Q(\cdot)$ being compatible with the charts. This is ensured because the densities $Q(U_{k})q_{k}(\cdot)$ and $Q(U_{\ell})q_{\ell}(\cdot)$ are consistent and align with the measure $Q$ over the overlapping regions of the charts. This compatibility is essential for constructing a coherent global measure from local chart densities.

A compact manifold ${\mathcal{Y}}$ admits a finite partition of unity $\{\tau_{k},k=1,\ldots,K\}$, each of which is sufficiently smooth (Lee,, ). Recall by definition of partition of unity we have, $\tau_{k}(\mathbf{u})=0$ for $\mathbf{u}\notin U_{k}$ and $\sum_{k=1}^{K}\tau_{k}(\mathbf{u})=1$, for all $\mathbf{u}\in{\mathcal{Y}}$. Since, $q(\mathbf{u})=Q(U_{k})q_{k}(\mathbf{u})$ for each $k$ and $\mathbf{u}\in U_{k}$, one may write $q(\mathbf{u})=\sum_{k=1}^{K}Q(U_{k})\tau_{k}(\mathbf{u})q_{k}(\mathbf{u})$. Denote $c_{k}=\int\tau_{k}(\mathbf{u})dQ_{k}(\mathbf{u})$ be normalizing constant and $q^{\prime}_{k}(\mathbf{u})=\tau_{k}(\mathbf{u})q_{k}(\mathbf{u})/c_{k}$. Then we have $q(\mathbf{u})=\sum_{k=1}^{K}=\pi_{k}q^{\prime}_{k}(\mathbf{u})$, where $\pi_{k}=c_{k}Q(u_{k})$. That is, $q$ is a mixture of $q\prime_{k}^{\prime}$s.

A compact manifold ${\mathcal{Y}}$ can be covered by a finite partition of unity $\{\tau_{k},k=1,\ldots,K\}$, each sufficiently smooth (Lee,, ). By definition, each function in this partition satisfies $\tau_{k}(\mathbf{u})=0$ for $\mathbf{u}\notin U_{k}$ and $\sum_{k=1}^{K}\tau_{k}(\mathbf{u})=1$ for all $\mathbf{u}\in{\mathcal{Y}}$. Given that $q(\mathbf{u})=Q(U_{k})q_{k}(\mathbf{u})$ for each $k$ and $\mathbf{u}\in U_{k}$, we can express $q(\mathbf{u})$ as:

$$ $q(\mathbf{u})=\sum_{k=1}^{K}Q(U_{k})\tau_{k}(\mathbf{u})q_{k}(\mathbf{u}).$ $$

To normalize, let $c_{k}=\int\tau_{k}(\mathbf{u})dQ_{k}(\mathbf{u})$ and define $q^{\prime}_{k}(\mathbf{u})=\tau_{k}(\mathbf{u})q_{k}(\mathbf{u})/c_{k}$. Thus, we can rewrite $q(\mathbf{u})$ as:

$$ $q(\mathbf{u})=\sum_{k=1}^{K}\pi_{k}q^{\prime}_{k}(\mathbf{u}),$ $$

where $\pi_{k}=c_{k}Q(U_{k})$. This formulation reveals that $q$ is a mixture of the component densities $q^{\prime}_{k}(\mathbf{u})$, weighted by $\pi_{k}$. This mixture approach ensures compatibility across different charts, providing a unified density representation over the entire manifold ${\mathcal{Y}}$.

With all the necessary tools at our disposal, we are ready to begin the construction of the desired map. We start by partitioning the measure on ${\mathcal{Z}}\times{\mathcal{X}}$ into subsets with weights $\{\pi_{k}\}_{k=1}^{K}$, denoted as $\{\mathcal{V}_{k}\}_{k=1}^{K}$. By following a similar approach to the single chart scenario and utilizing Lemma 2, we obtain the push-forward maps $g_{k}\circ\varphi_{k}:\mathcal{V}_{k}\subset{\mathcal{Z}}\times{\mathcal{X}}
\to U_{k}\subset{\mathcal{Y}}$, which transport the normalized measure from $\mathcal{V}_{k}$ to the measure $Q^{\prime}_{k}$, corresponding to the density $q^{\prime}_{k}$ for each $1\leq k\leq K$.

The final step is to combine these local maps into a single coherent global map. This is achieved by summing the local maps weighted by their corresponding measures:

$$ $G=\sum_{k=1}^{K}\pi_{k}g^{\prime}_{k},$ $$

where $g^{\prime}_{k}$ represents $g_{k}\circ\varphi_{k}$ on $\mathcal{V}_{k}$ and is defined as 0 outside $\mathcal{V}_{k}$. This ensures that the push-forward map accurately transports the measure from ${\mathcal{Z}}\times{\mathcal{X}}$ to ${\mathcal{Y}}$ while maintaining the necessary smoothness and compatibility across different charts.

## Appendix D Proof of Lemma 1

###### Proof.

For $g_{1}(\cdot|x),g_{2}(\cdot|x)\in\mathcal{F}$ with $\||g_{1}-g_{2}|_{\infty}\|_{\infty}\leq\eta_{1}$. Then

$$ $\displaystyle p_{g_{1},\sigma}(y|x)-p_{g_{2},\sigma}(y|x)$ $\displaystyle=$ $\displaystyle\int\phi_{\sigma}(y-g_{1}(x,z))\left(1-\frac{\phi_{\sigma}(y-g_{2 }(x,z))}{\phi_{\sigma}(y-g_{1}(x,z))}\right)dP_{Z}(z)$ $\displaystyle=$ $\displaystyle\int\phi_{\sigma}(y-g_{1}(x,z))\left(1-\exp{\left\{-\frac{|y-g_{2 }(x,z)|^{2}_{2}-|y-g_{1}(x,z)|^{2}_{2}}{2\sigma^{2}}\right\}}\right)dP_{Z}(z)$ $\displaystyle\leq$ $\displaystyle\int\phi_{\sigma}(y-g_{1}(x,z))\left(\frac{|y-g_{2}(x,z)|^{2}_{2} -|y-g_{1}(x,z)|^{2}_{2}}{2\sigma^{2}}\right)dP_{Z}(z)$ (18) $\displaystyle=$ $\displaystyle\int\phi_{\sigma}(y-g_{1}(x,z))\left(\frac{|g_{2}(x,z)-g_{1}(x,z) |^{2}_{2}-2(y-g_{1}(x,z))^{T}(g_{2}(x,z)-g_{1}(x,z))}{2\sigma^{2}}\right)dP_{Z }(z)$ $\displaystyle\leq$ $\displaystyle\int\phi_{\sigma}(y-g_{1}(x,z))\left(\frac{|g_{2}(x,z)-g_{1}(x,z) |^{2}_{2}}{2\sigma^{2}}+\frac{2|y-g_{1}(x,z)|_{1}|g_{2}(x,z)-g_{1}(x,z)|_{ \infty}}{2\sigma^{2}}\right)dP_{Z}(z)$ $\displaystyle\leq$ $\displaystyle\int\phi_{\sigma}(y-g_{1}(x,z))\frac{2KD\eta_{1}}{2\sigma^{2}}dP_ {Z}(z)+\frac{2\eta_{1}}{2\sigma^{2}}\int|y-g_{1}(x,z)|_{1}\phi_{\sigma}(y-g_{1 }(x,z))dP_{Z}(z)$ (19) $\displaystyle\leq$ $\displaystyle\frac{2KD\eta_{1}}{2\sigma^{2}}\frac{1}{\left(\sqrt{2\pi\sigma^{2 }}\right)^{D}}+\frac{\eta_{1}}{\sigma^{2}}\int\sqrt{\frac{D}{2\pi e}}\frac{1}{ (\sqrt{2\pi\sigma^{2}})^{D-1}}dP_{Z}(z)$ (20) $\displaystyle\leq$ $\displaystyle c_{1}(K,D)\sigma_{\min}^{-(D+2)}\eta_{1}.$ (21) $$

For the last line, we use the fact that $\sigma_{\min}\leq 1$. The inequality at (18) follows from $e^{-x}\geq(1-x)$. The ones at (19) follows using

$$ $\displaystyle|g_{2}(x,z)-g_{1}(x,z)|_{2}^{2}\leq 2K|g_{2}(x,z)-g_{1}(x,z)|_{1}$ $\displaystyle\leq 2KD|g_{2}(x,z)-g_{1}(x,z)|_{\infty}$ $\displaystyle\leq 2KD\||g_{1}-g_{2}|_{\infty}\|_{\infty}\leq 2KD\eta_{1}$ $$

and $|g_{2}(x,z)-g_{1}(x,z)|_{\infty}\leq\eta_{1}$. The change at (20) follows from $\phi_{\sigma}(y-g_{1}(x,z))\leq\left(\sqrt{2\pi\sigma^{2}}\right)^{-D}$ and the bound

$$ $|v|_{1}\phi_{\sigma}(v)\leq\sqrt{\frac{D}{2\pi e}}\frac{1}{(\sqrt{2\pi\sigma^{ 2}})^{D-1}}.$ $$

Now for $\sigma_{1},\sigma_{2}\in\left[\sigma_{\min},\sigma_{\max}\right]$ with $|\sigma_{1}-\sigma_{2}|\leq\eta_{2}$. It holds that $\left|\sigma_{1}^{-2}-\sigma_{2}^{-2}\right|\leq\sigma_{1}^{-2}\sigma_{2}^{-2}
\left(\sigma_{1}+\sigma_{2}\right)\eta_{2}$ and $\left|\log\left(\frac{\sigma_{2}}{\sigma_{1}}\right)\right|\leq\frac{\eta_{2}}
{\min\{\sigma_{1},\sigma_{2}\}}$. We have

$$ $\displaystyle p_{g,\sigma_{1}}(y|x)-p_{g_{2},\sigma_{2}}(y|x)$ $\displaystyle=$ $\displaystyle\int\phi_{\sigma_{1}}(y-g(x,z)\left(1-\left(\frac{\sigma_{1}}{ \sigma_{2}}\right)^{D}\exp\left\{\frac{|y-g(x,z)|_{2}^{2}}{2}\left(\frac{1}{ \sigma_{1}^{2}}-\frac{1}{\sigma_{2}^{2}}\right)\right\}\right)dP_{Z}(z)$ $\displaystyle\leq$ $\displaystyle\int\phi_{\sigma_{1}}(y-g(x,z)\left[\frac{|y-g(x,z)|_{2}^{2}}{2} \left(\frac{1}{\sigma_{2}^{2}}-\frac{1}{\sigma_{1}^{2}}\right)-D\log\left( \frac{\sigma_{1}}{\sigma_{2}}\right)\right]dP_{Z}(z)$ (22) $\displaystyle\leq$ $\displaystyle\int\phi_{\sigma_{1}}(y-g(x,z)\left[\frac{|y-g(x,z)|_{2}^{2}}{2} \left(\frac{\sigma_{1}+\sigma_{2}}{\sigma_{1}^{2}\sigma_{2}^{2}}\right)\eta_{2 }+\frac{D\eta_{2}}{\min\{\sigma_{1},\sigma_{2}\}}\right]dP_{Z}(z)$ $\displaystyle\leq$ $\displaystyle\frac{1}{(\sqrt{2\pi\sigma_{1}^{2}})^{D}}\frac{\sigma_{1}+\sigma_ {2}}{e\sigma_{2}^{2}}\eta_{2}+\frac{1}{\left(\sqrt{2\pi\sigma_{1}^{2}}\right)^ {D}}\frac{D\eta_{2}}{\min\{\sigma_{1},\sigma_{2}\}}$ (23) $\displaystyle\leq$ $\displaystyle c_{2}(D)\sigma_{\min}^{-(D+1)}\eta_{2}.$ (24) $$

The (22) follows from $1-e^{-\alpha}\leq\alpha$. The change at (23) follows from $\phi_{\sigma_{1}}(y-g(x,z))\leq\left(\sqrt{2\pi\sigma_{1}^{2}}\right)^{-D}$ and

$$ $|v|_{2}^{2}\phi_{\sigma}(v)\leq\frac{\sigma^{2}}{(\sqrt{2\pi\sigma^{2}})^{D}} \frac{2}{e}.$ $$

Let $\varepsilon>0$. Let $\{g_{1},\ldots,g_{N_{1}}\}$ be $\eta_{1}-$covering of $\mathcal{F}$ and $\{\sigma_{1},\ldots,\sigma_{N_{2}}\}$ be $\eta_{2}-$covering of $\left[\sigma_{\min},\sigma_{\max}\right]$ with respect to $\||\cdot|_{\infty}\|_{\infty}$ and $|\cdot|_{\infty}$. By (21) and (24), $\eta_{1}=c_{1}^{-1}\sigma_{\min}^{D+2}\varepsilon/4$ and $\eta_{2}=c_{2}^{-2}\sigma_{\min}^{D+1}\varepsilon/4$ implies

$$ $\left\{P_{g_{i},\sigma_{j}}(\cdot|\cdot):i=1,\ldots,N_{1},j=1,\ldots,N_{2}\right\}$ $$

forms an $\varepsilon/2-$covering for ${\mathcal{P}}$ with respect to $\|\cdot\|_{\infty}$. Denote the envelope function of $\mathcal{F}$

$$ $\displaystyle H(y,x)=\sup_{p\in{\mathcal{P}}}p(y|x)$ $\displaystyle\leq\frac{1}{\left(2\pi\sigma^{2}_{\min}\right)^{-D/2}}\exp\left \{-\frac{|y|_{2}^{2}-4K^{2}D}{4\sigma^{2}_{\max}}\right\}$ $\displaystyle=e^{K^{2}D/2\sigma_{\max}^{2}}2^{D/2}\left(\frac{\sigma_{\max}}{ \sigma_{\min}}\right)^{D}\phi_{\sqrt{2}\sigma_{\max}}(y).$ $$

Following from $\int_{|y|_{\infty}>t}\phi_{\sigma}(y)dy\leq 2De^{-t^{2}/2\sigma^{2}}$, we have

$$ $\int\int_{|y|_{\infty}>B}H(y,x)\mu(y,x)dydx=\int\left(\int_{|y|_{\infty}>B}H(y ,x)\mu(y|x)dy\right)\mu_{X}^{\ast}(x)dx<\varepsilon,$ $$

where

$$ $B=2\sigma_{\max}\left(\log\frac{1}{\varepsilon}+D\log\frac{\sigma_{\max}}{ \sigma_{\min}}+\frac{K^{2}D}{2\sigma_{\max}^{2}}+\log 2D\right)^{1/2}.$ $$

For each $(i,j)$ define

$$ $l_{ij}(y,x)=\max\left\{p_{g_{i},\sigma_{j}}(y,x)-\varepsilon/2,0\right\}\qquad \text{and}\qquad u_{ij}(y,x)=\min\left\{p_{g_{i},\sigma_{j}}(y,x)+\varepsilon/ 2,H(y,x)\right\}.$ $$

It follows that

$$ $\begin{split}&\int\int\left\{u_{ij}(y,x)-l_{ij}(y,x)\right\}\mu_{X}^{\ast}(x) dydx\\ \leq&\int\int_{|y|_{\infty}\leq B}\varepsilon\mu_{X}^{\ast}(x)dydx+\int\int_{| y|_{\infty}>B}H(y,x)\mu_{X}^{\ast}(x)dydx\\ \leq&\left\{(2B)^{D}+1\right\}\varepsilon.\end{split}$ (25) $$

Denote $\delta^{2}:=\left\{(2B)^{D}+1\right\}$. With $d_{H}^{2}(u_{ij},l_{ij})\leq d_{1}(u_{ij},l_{ij})$, we have

$$ ${\mathcal{N}}_{[]}(\delta,{\mathcal{P}},d_{H})\leq{\mathcal{N}}_{[]}(\delta^{2 },{\mathcal{P}},d_{1})\leq N_{1}N_{2}\leq\frac{\sigma_{\max}-\sigma_{\min}}{ \eta_{2}}{\mathcal{N}}(\eta_{1},\mathcal{F},\||\cdot|_{\infty}\|_{\infty}).$ (26) $$

It is possible to write

$$ $\delta^{2}=\varepsilon\leq C_{1}(\sigma_{\max},D)\left[\varepsilon(\log \varepsilon^{-1})^{D/2}+\varepsilon C_{2}(K)+\varepsilon\left(\log\frac{\sigma _{\max}}{\sigma_{\min}}\right)^{D/2}\right],$ $$

where $C_{1}(\sigma_{\max},D)$ and $C_{2}(K)$ is a constant. There exists small enough $\varepsilon_{*}(D)$ such that for all $\varepsilon\in(0,\varepsilon_{*}]$

$$ $\delta^{2}\leq C_{3}(\sigma_{\max},D,K)\sqrt{\varepsilon}\left(\log\frac{ \sigma_{\max}}{\sigma_{\min}}\right)^{D/2}.$ $$

Consequently, there exists $\delta_{*}=\delta_{*}(D)$, such that for all $\delta\leq\delta_{*}$, we have

$$ $C_{3}^{2}(\sigma_{\max},K,D)\delta^{4}\left(\log\frac{\sigma_{\max}}{\sigma_{ \min}}\right)^{-D}\leq\varepsilon.$ $$

It lead us to, for all $\delta\leq\delta_{*}$

$$ $\eta_{1}\geq\frac{c_{1}^{-1}C_{3}^{2}\sigma_{\min}^{D+3}\delta^{4}}{\sigma_{ \min}\{\log(\sigma_{\max}/\sigma_{\min})\}^{D}}\geq c\sigma_{\min}^{D+3}\delta ^{4},$ (27) $$

where $c(\sigma_{\max},K,D)$ is a constant. We use the fact that $\sigma_{\min}\{\log(\sigma_{\max}/\sigma_{\min})\}^{D}$ is bounded above by some constant depending only upon $\sigma_{\max}$ as $\sigma_{\min}\leq 1$. Similar to (27), it is possible to write for all $\delta>\delta_{*}$

$$ $\eta_{2}\geq c^{\prime}\sigma_{\min}^{D+2}\delta^{4},\qquad\qquad\text{ for all }\delta\leq\delta_{*},$ (28) $$

where $c^{\prime}(\sigma_{\max},K,D)$ is some constant.

The result now follows directly (28) and (27) with (26).
∎

###### Proof.

## Appendix E Proof of Theorem 1

###### Proof.

Choose four absolute constants $c_{1},\ldots,c_{4}$ as in Theorem 1 of Wong and Shen, (). Define $c$ and $C$ in the statement of Lemma 1. The proof closely follows Chae et al., (). We have therein the proof of Theorem 3 that

$$ $\begin{split}&\int_{\varepsilon^{2}/2^{8}}^{\sqrt{2}\varepsilon}\sqrt{\log{ \mathcal{N}}_{[]}(\delta/c_{3},\mathcal{P},d_{H})}d\delta\\ \leq&\,\sqrt{2}\varepsilon\sqrt{\xi A+(D+3)(s+1)\log\sigma_{min}^{-1}+c_{5}\xi }+\sqrt{2}\varepsilon\sqrt{4(\xi+1)}\sqrt{\log(2^{8}/\varepsilon^{2})},\end{split}$ (29) $$

for every $\varepsilon\leq\sqrt{2}\leq c_{3}\delta_{*}/\sqrt{2}$, where $c_{5}=c_{5}(c,C,c_{3})$. Observe that $c_{4}\sqrt{n}\varepsilon_{n}^{2}$ is upper bound to (29) and Eq. (3.1) of Wong and Shen, () is satisfied.

Using B.12 of Ghosal and van der Vaart, (), we have

$$ $\displaystyle K(p_{G_{*},\sigma_{*}},p_{g,\sigma_{*}})$ $\displaystyle\leq\int\int K\Big{(}N\left(G_{*}(z,x),\sigma_{*}^{2}\right),N \left(g(z,x),\sigma_{*}^{2}\right)\Big{)}\mu_{X}^{\ast}(x)\thinspace dx \thinspace dP_{Z}(z)$ $\displaystyle=\int\int\frac{|G_{*}(z,x)-g(z,x)|^{2}_{2}}{2\sigma_{*}^{2}}\mu_{ X}^{\ast}(x)\thinspace dx\thinspace dP_{Z}(z)\leq\frac{D\delta_{\operatorname* {approx}}^{2}}{2\sigma_{*}^{2}}=:\delta_{n}.$ $$

One may easily see that

$$ $\int\left(\log\frac{\phi_{\sigma}(x)}{\phi_{\sigma}(x-y)}\right)^{2}\phi_{ \sigma}(x)dx=\int\frac{|y|^{4}_{2}+4|x^{T}y|^{2}}{4\sigma^{2}}\phi_{\sigma}(x) dx\leq\frac{|y|_{2}^{4}}{4\sigma^{2}}+|y|_{2}^{2}\int\frac{|x|_{2}^{2}}{\sigma ^{2}}\phi_{\sigma}(x)dx.$ $$

Combining this with Example B.12, (B.17) and Exercise B.8 of Ghosal and van der Vaart, (), we have

$$ $\displaystyle\int\int\left(\log\frac{p_{G_{*},\sigma_{*}}(y|x)}{p_{g,\sigma_{* }}(y|x)}\right)^{2}\,dP_{*}(y|x)\,\mu_{X}^{\ast}(x)dx$ $\displaystyle\leq$ $\displaystyle\int\int\int\left(\log\frac{\phi_{\sigma}(y-G_{*}(z,x)}{\phi_{ \sigma}(y-G(z,x)}\right)^{2}\phi_{\sigma}(y-G_{*}(z,x))\,dy\thinspace dP_{Z}(z )\thinspace\mu_{X}^{\ast}(x)dx$ $\displaystyle\leq$ $\displaystyle\frac{D^{2}\delta^{4}_{\operatorname*{approx}}}{4\sigma_{*}^{2}}+ D\delta^{2}_{\operatorname*{approx}}\int\frac{|x|^{2}_{2}}{\sigma_{*}^{2}}\phi _{\sigma_{*}}(y)dy+\frac{2D\delta^{2}_{\operatorname*{approx}}}{\sigma_{*}^{2} }\leq c_{7}\frac{\delta_{\operatorname*{approx}}^{2}}{\sigma_{*}^{2}}=:\tau_{n},$ $$

where $c_{7}=c_{7}(D)$. We are using $\delta_{n}$ and $\tau_{n}$, although they are independent of $n$, for notational consistency with Theorem 4 of Wong and Shen, (). Let $\varepsilon_{n}^{*}=\varepsilon_{n}\vee\sqrt{12\delta_{n}}$. Then, using Theorem 4 of Wong and Shen, (), we have

$$ $P_{*}\left(d_{H}(\widehat{p},p_{*})>\varepsilon_{n}\right)\leq 5e^{-c_{2}n \varepsilon_{n}^{*2}}+\frac{\tau_{n}}{n\delta_{n}}=5e^{-c_{2}n\varepsilon_{n}^ {*2}}+\frac{2c_{7}^{2}}{Dn}.$ $$

The proof is complete after redefining constants.
∎

###### Proof.

## Appendix F Proofs of Corollary 1

###### Proof.

For the sparse case in 1.1, utilizing the entropy bound from (10), we observe that

$$ $\xi\{A+\log(n/\sigma_{\min})\}\asymp\delta_{\operatorname*{approx}}^{-t_{*}/ \beta_{*}}\log^{3}(\delta^{-1}_{\operatorname*{approx}}),$ $$

which naturally leads to the required convergence rate.

Similarly for the fully connected case 1.2, utilizing the entropy bound from (11) , we observe that

$$ $\xi\{A+\log(n/\sigma_{\min})\}\asymp\delta_{\operatorname*{approx}}^{-t_{*}/ \beta_{*}}\log^{3}(\delta^{-1}_{\operatorname*{approx}}),$ $$

which naturally leads to the required convergence rate.
∎

###### Proof.

## Appendix G Proof of Theorem 2

###### Proof.

It is suffice to assume that $\varepsilon$ and $\sigma_{*}\sqrt{\log\varepsilon^{-1}}$ are sufficiently small. If not, let $\varepsilon+\sigma_{*}\sqrt{\log\varepsilon^{-1}}\geq c_{0}$, where $c_{0}(K,D,\mathsf{r}_{\ast})$. Then Theorem 2 holds trivially by taking a large enough constant depending just on $D$, $K$, and $\mathsf{r}_{\ast}$.

Let $V\sim Q(\cdot|X=x)$, $V_{*}\sim Q(\cdot|X=x)$, ${\bm{\epsilon}}\sim{\mathsf{N}}(0_{D},\sigma^{2}\mathbb{I}_{d})$ and ${\bm{\epsilon}}_{*}\sim{\mathsf{N}}(0_{D},\sigma_{*}^{2}\mathbb{I}_{d})$ be independent with underlying probability density $\nu$. We truncate the random variable ${\bm{\epsilon}}$ and ${\bm{\epsilon}}_{*}$ componentwise as $\left({\bm{\epsilon}}_{K}\right)_{j}=\max\{-K,\min\{K,{\bm{\epsilon}}_{j}\}\}$ and $\left({\bm{\epsilon}}_{*K}\right)_{j}=\max\{-K,\min\{K,\left({\bm{\epsilon}}_{
*}\right)_{j}\}\}$ respectively. We denote $P_{g,\sigma}$ as $P$, $Q_{g}$ as $Q$, $\widetilde{P}$ as distribution of $V+{\bm{\epsilon}}_{K}$ and $\widetilde{P}_{*}$ as the distribution of $V_{*}+{\bm{\epsilon}}_{*K}$. One may note that $W_{1}(\widetilde{P}_{*},Q_{*})\leq W_{2}(\widetilde{P}_{*},Q_{*})\leq\sqrt{
\mathbb{E}\big{[}|{\bm{\epsilon}}_{*K}|_{2}^{2}\big{]}}\leq\sqrt{\mathbb{E}
\big{[}|{\bm{\epsilon}}_{*}|_{2}^{2}\big{]}}\leq\sigma_{*}\sqrt{D}$. Similarly, $W_{1}(\widetilde{P},Q)\leq\sigma\sqrt{D}$. The $\ell_{1}$ diameter of $[-2K,2K]^{D}$, where the support of $\widetilde{P}$ and $\widetilde{P}_{*}$, is 4KD. Observe that

$$ $W_{1}\left(\widetilde{P}_{*},\widetilde{P}\right)\leq 4\,KD\,d_{1}\left( \widetilde{P}_{*},\widetilde{P}\right)\leq 4\,KD\,d_{1}(P_{*},P)\leq 8\,KD\,d_ {H}(P_{*},P),$ $$

where the first inequality follows from Theorem 4 of Gibbs and Su, (), the second inequality follows from the fact the distance between two truncated distributions is always lesser than the original distributions and the last inequality follows from $d_{1}\leq 2d_{H}$. Hence,

$$ $W_{1}\left(Q_{*},Q\right)\leq W_{2}\left(Q_{*},\widetilde{P}_{*}\right)+W_{1} \left(\widetilde{P}_{*},\widetilde{P}\right)+W_{2}\left(\widetilde{P},Q\right) \leq\sigma_{*}\sqrt{D}+8\,KD\,\varepsilon+\sigma\sqrt{D}.$ $$

Now it is suffice to show that $\sigma\leq c\;\sigma_{*}\sqrt{\log\varepsilon^{-1}}$, where $c=c(D,K,r*)$ is a constant, because we have assumed that $\varepsilon$ is small enough. We establish this in the rest of the proof. Let $t_{*}=\left[2\,\sigma_{*}^{2}\,D\log\left(\frac{2D}{\varepsilon}\right)\right]
^{1/2}$. Observe that

$$ $\int_{|x|_{2}>t_{*}}\phi_{\sigma_{*}}(x)dx\leq\int_{|x|_{\infty}>t_{*}/\sqrt{D }}\phi_{\sigma_{*}}(x)dx\leq 2\,De^{-t_{*}^{2}/2D\sigma^{2}}\leq\varepsilon.$ $$

Let $\mathcal{M}_{*}^{t_{*}}=\mathcal{M}_{*}\oplus{\mathcal{B}}_{t_{*}}(0_{D})$. We may write

$$ $\displaystyle 1-P_{*}\left(\mathcal{M}_{*}^{t_{*}}\right)=\nu\left(Y_{*}+{\bm{ \epsilon}}_{*}\notin\mathcal{M}_{*}^{t_{*}}\right)\leq\nu\left(|{\bm{\epsilon} }_{*}|_{2}>t_{*}\right)$ (30) $\displaystyle\implies$ $\displaystyle P\left(\mathcal{M}_{*}^{t_{*}}\right)\geq 1-2\varepsilon,$ $$

the implication in the last line follows from $\sup_{B}\left|P(B)-P_{*}(B)\right|\leq d_{H}(P,P_{*})\leq\varepsilon$. For the sake of contradiction, let $\sigma\in\left[2t_{*},r^{*}/2\right]\cup\left(\mathsf{r}_{\ast}/2,\infty\right)$ ($t_{*}$ is sufficiently small, from the assumption we made at the beginning of this proof). If $\sigma>\mathsf{r}_{\ast}/2$, then

$$ $2\varepsilon\geq 1-P\left(\mathcal{M}_{*}^{t_{*}}\right)\geq 1-P\left([-K,K]^{ D}\right)\geq c_{2}(K,D,r*)$ $$

where $c_{2}$ is some positive constant. It is a contradiction following from the smallness of $\varepsilon$. Lets make a claim that if $\sigma\in\left[2t_{*},\mathsf{r}_{\ast}/2\right]$, then for every $y\in\mathbb{R}^{D}$, there is some $z\in\mathbb{R}^{D}$ such that $|z-y|_{2}\leq\sigma$ and ${\mathcal{B}}_{\sigma/2}(z)\cap\mathcal{M}_{*}^{t_{*}}=\emptyset$.

Following from the claim, we have

$$ $\nu\left(Y+{\bm{\epsilon}}\notin\mathcal{M}_{*}^{t_{*}}\big{|}Y=y\right)\geq \nu\left({\bm{\epsilon}}\in{\mathcal{B}}_{\sigma/2}(z-y)\right).$ $$

Since $|z-y|_{2}\leq\sigma$, the right hand side is bounded below by a positive constant depending just on $D$ which is again a contradiction to (30). This proves the assertion made in the theorem.

The proof of the claim is divided into three cases. Let $\rho\left(y,\mathcal{M}_{*}\right)=\inf\{|y-y^{\prime}|_{2}:y^{\prime}\in
\mathcal{M}_{*}\}$ be the $\ell_{2}$ set distance.

Case 1. $\rho(y,\mathcal{M}_{*})\geq\sigma:$ We may choose $z=y$.

Case 2. $\rho(y,\mathcal{M}_{*})\in(0,\sigma):$ Let $y_{0}$ be the unique Euclidean projection of $y$ onto $\mathcal{M}_{*}$. Such a unique projection exists because $\sigma<\mathsf{r}_{\ast}$ is within the reach and $y\in\mathcal{M}_{*}$, since $\mathcal{M}_{*}$ is closed. Suppose $y_{t}=y_{0}+t(y-y_{0})$. We shall define two continuous functions $d_{0}(t)=|y_{t}-y_{0}|_{2}$ and $d(t)=\rho(y_{t},\mathcal{M}_{*})$. It is obvious that $d(t)\leq d_{0}(t)$. For $t\in\big{[}0,1+\sigma/|y-y_{0}|_{2}\big{]}$, $d_{0}(t)\leq d(t)$ because $y_{0}$ is the unique projection for all the points that lie on the line segment including the farthest point with $t=1+\sigma/|y-y_{0}|_{2}$. Otherwise, say $d(t)=\rho(y_{t},z)$ and

$$ $|y-y_{0}|_{2}=|y-y_{t}|_{2}+|y_{t}-y_{0}|_{2}>|y-y_{t}|+|y_{t}-z|\geq|y-z|_{2}$ $$

which contradicts $y_{0}$ being a unique projection. The claim holds for the point $z=y_{1+\sigma/|y-y_{0}|_{2}}$. To see this, observe $|z-y|=\sigma$ and ${\mathcal{B}}_{\sigma/2}(z)\cap\mathcal{M}_{*}^{t_{*}}=\emptyset$ because $t_{*}\leq\sigma/2$ and the ball ${\mathcal{B}}_{\sigma/2}(z)\subset\mathcal{M}_{*}^{\mathsf{r}_{\ast}}$ is within the reach of the manifold.

Case 3. $\rho(y,\mathcal{M}_{*})=0:$ Because $\mathcal{M}_{*}$ has empty interior, for all $\gamma>0$, we always find a point $y_{\gamma}$, which in ${\mathcal{B}}_{\gamma}(y)$ which away from $\mathcal{M}_{*}$. For small enough $\gamma$, we reduce to case 2 by taking $\gamma\to 0$, the limit point of $y_{\gamma}$ has the required behavior.

∎

###### Proof.

## Appendix H Proof of Corollary 2

###### Proof.

The effective noise variance after the perturbation would be

$$ $\displaystyle\widetilde{\sigma}_{*}=n^{-\alpha}+n^{-\beta_{*}/2(\beta_{*}+t_{* })}\asymp\begin{cases}n^{-\alpha},\qquad&\alpha<\beta_{*}/\{2(\beta_{*}+t_{*}) \}\\ n^{\beta_{*}/2(\beta_{*}+t_{*})},\qquad&\text{otherwise}.\end{cases}$ $$

Following this and the Theorem 2, for the rate we have

$$ $\displaystyle\varepsilon^{*}_{n}+\sigma_{*}\sqrt{\log((\varepsilon^{*}_{n})^{- 1})}$ $\displaystyle\asymp\left(n^{-\frac{\beta_{*}-t_{*}\alpha}{2\beta_{*}+t_{*}}}+n ^{-\alpha}\right)\log^{2}(n)$ $\displaystyle\asymp\begin{cases}n^{-\frac{\beta_{*}-t_{*}\alpha}{2\beta_{*}+t_ {*}}}\,\log^{2}(n),\quad&\text{if}\;\alpha<\beta_{*}/\{2(\beta_{*}+t_{*})\},\\ n^{-\frac{\beta_{*}}{2(\beta_{*}+t_{*})}}\,\log^{2}(n),\quad&\text{otherwise}. \end{cases}$ $$

∎

###### Proof.

## Appendix I Proof of Theorem 3

###### Proof.

With $m=\lceil\log_{2}(n)\rceil$ and $N=\left(n^{(\beta_{Z}^{-1}d+\beta_{X}^{-1}p)\left[1+\alpha(\beta_{Z}^{-1}d+
\beta_{X}^{-1}p)\right]/\left[2+\beta_{Z}^{-1}d+\beta_{X}^{-1}p\right]}\right)$ in Theorem 5, we can find a network $G$ with the mentioned architecture such that

$$ $\||G-G_{*}|_{\infty}\|_{\infty}\leq\delta_{\operatorname*{approx}}.$ $$

Following the entropy bound from (10), we have

$$ $\displaystyle\log{\mathcal{N}}(\delta,\mathcal{F}_{s},\||\cdot|_{\infty}\|_{ \infty})$ $\displaystyle\lesssim sL\,\{\log(rL)+\log\delta^{-1}\}$ $\displaystyle\lesssim\delta^{-(\beta_{Z}^{-1}d+\beta_{X}^{-1}p)}_{ \operatorname*{approx}}\log^{2}\delta^{-1}_{\operatorname*{approx}}\left\{\ \log\left(\delta^{-1}_{\operatorname*{approx}}\log\left(\delta^{-1}_{ \operatorname*{approx}}\right)\right)+\log\left(\delta^{-1}_{\operatorname*{ approx}}\right)\right\}.$ $$

The rest directly follows from the Theorem 1
∎

###### Proof.

## Appendix J Approximation properties of the sparse and fully connected DNNs

The approximability of the sparse network is detailed in Lemma 4.1, which restates Lemma 5 from Chae et al., (). For the fully connected network, Lemma 4.2 demonstrates its approximation capabilities, derived directly from Theorem 2 and the proof of Theorem 1 in Kohler and Langer, (). Additionally, the inclusion of the class $\mathcal{G}$ in the fully connected setup is supported by the discussion in Section 1 of Kohler and Langer, ().

###### Lemma 4 .

Suppose that $G_{*}\in\mathcal{G}$. Then, for every small enough $\delta\in(0,1)$,

- 1.
there exists a sparse network $G\in\mathcal{F}_{s}=\mathcal{F}_{s}\left(L,r,s,K\vee 1\right)$ with $L\lesssim\log\delta^{-1}$, $r\lesssim\delta^{-t_{*}/\beta_{*}}$, $s\lesssim\delta^{-t_{*}/\beta_{*}}\log\delta^{-1}$ satisfying $\||G-G_{*}|_{\infty}\|_{\infty}\leq\delta$.
- 2.
there exists a fully connected network $G\in\mathcal{F}_{c}$ with $L\lesssim\log\delta^{-1}$, $r\lesssim\delta^{-t_{*}/2\beta_{*}}$, $B\lesssim\delta^{-1}$ satisfying $\left\|\left|G-G_{*}\right|_{\infty}\right\|_{\infty}\leq\delta$.

###### Lemma 4 .

## Appendix K A new approximation result for functions with smoothness disparity

In this section, we prove the approximability of the sparse neural network for the Hölder class of function $f\in{\mathcal{H}}_{{r},{r_{\prime}}}^{\beta,\beta_{\prime}}\left(D,D_{\prime},
K\right)$.

###### Theorem 5 .

Let $f\in{\mathcal{H}}_{{r},{r_{\prime}}}^{\beta,\beta_{\prime}}\left([0,1]^{{r}},[
0,1]^{{r_{\prime}}},K\right)$. Denote ${\mathsf{r}_{\mathrm{sum}}}=r+r_{\prime}$ and ${\beta_{\mathrm{sum}}}=\beta+\beta_{\prime}$. Then for any integers $m\geq 1$ and $N\geq({\beta_{\mathrm{sum}}}+1)^{{\mathsf{r}_{\mathrm{sum}}}}\vee(K+1)e^{{
\mathsf{r}_{\mathrm{sum}}}},$ there exists a network

$$ $\widetilde{f}\in\mathcal{F}_{s}\big{(}L,\big{(}{\mathsf{r}_{\mathrm{sum}}},6({ \mathsf{r}_{\mathrm{sum}}}+\lceil{\beta_{\mathrm{sum}}}\rceil)N,\ldots,6({ \mathsf{r}_{\mathrm{sum}}}+\lceil{\beta_{\mathrm{sum}}}\rceil)N,1\big{)},s, \infty\big{)}$ $$

with depth

$$ $L=8+(m+5)\left(1+\left\lceil\log_{2}\big{(}{\mathsf{r}_{\mathrm{sum}}}\vee{ \beta_{\mathrm{sum}}}\big{)}\right\rceil\right)$ $$

and the number of parameters

$$ $\displaystyle s\leq 109\big{(}{\mathsf{r}_{\mathrm{sum}}}+{\beta_{\mathrm{sum} }}+1\big{)}^{3+{\mathsf{r}_{\mathrm{sum}}}}N(m+6),$ $$

such that

$$ $\displaystyle\|\widetilde{f}-f\|_{L^{\infty}([0,1]^{{\mathsf{r}_{\mathrm{sum}} }})}\leq(2K+1)\left(1+{\mathsf{r}}_{\mathrm{sum}}^{2}+\beta_{\mathrm{sum}}^{2} \right)6^{{\mathsf{r}_{\mathrm{sum}}}}\,N\,2^{-m}+K\,3^{{\mathsf{r}_{\mathrm{ sum}}}/(\beta^{-1}r+\beta_{\prime}^{-1}r_{\prime})}\,N^{-1/(\beta^{-1}r+\beta_ {\prime}^{-1}r_{\prime})}.$ $$

We denote $\widetilde{\beta}=(\beta+\beta_{\prime})^{-1}\beta\beta_{\prime}$ and $\widetilde{r}=(\beta+\beta_{\prime})^{-1}(r\beta+r_{\prime}\beta_{\prime})$.
Before presenting the proof of Theorem 5, we formulate some required results.

We follow the classical idea of function approximation by local Taylor approximations that have previously been used for network approximations in Yarotsky, () and Schmidt-Hieber, (). For a vector $\mathbf{a}\in[0,1]^{r}$ define

$$ $\displaystyle P_{\mathbf{a},\mathbf{b}}^{\beta,\beta_{\prime}}f(\mathbf{u}, \mathbf{v})=\sum_{\begin{subarray}{c}0\leq|\bm{\alpha}|<\beta\\ 0\leq|\bm{\alpha}_{\prime}|<\beta_{\prime}\end{subarray}}(\partial^{\bm{\alpha }+\bm{\alpha}_{\prime}}f)(\mathbf{a},\mathbf{b})\frac{(\mathbf{u}-\mathbf{a})^ {\bm{\alpha}}(\mathbf{v}-\mathbf{b})^{\bm{\alpha}_{\prime}}}{\bm{\alpha}!\,\bm {\alpha}_{\prime}!}.$ (31) $$

We use the notation the $\mathbf{u}=(u^{(j)})_{j}$ to represent the component of the vector when the index $j$ is well understood. Accordingly we have $\mathbf{v}=(v^{(j)})_{j}$, $\mathbf{a}=(a^{(j)})_{j}$ and $\mathbf{b}=(b^{(j)})_{j}$ . By Taylor’s theorem for multivariate functions, we have for a suitable $\xi\in[0,1],$

$$ $\displaystyle f(\mathbf{u},\mathbf{v})=$ $\displaystyle\sum_{\begin{subarray}{c}\bm{\alpha}:|\bm{\alpha}|<\beta-1\\ \bm{\alpha}_{\prime}:|\bm{\alpha}_{\prime}|<\beta_{\prime}-1\end{subarray}}( \partial^{\bm{\alpha}+\bm{\alpha}_{\prime}}f)(\mathbf{a},\mathbf{b})\frac{( \mathbf{u}-\mathbf{a})^{\bm{\alpha}}(\mathbf{v}-\mathbf{b})^{\bm{\alpha}_{ \prime}}}{\bm{\alpha}!\,\bm{\alpha}_{\prime}!}$ $\displaystyle\,+$ $\displaystyle\sum_{\begin{subarray}{c}\beta-1\leq|\bm{\alpha}|<\beta\\ \beta_{\prime}-1\leq|\bm{\alpha}_{\prime}|<\beta_{\prime}\end{subarray}}( \partial^{\bm{\alpha}+\bm{\alpha}_{\prime}}f)(\mathbf{a}+\xi(\mathbf{u}- \mathbf{a}),\mathbf{b}+\xi(\mathbf{v}-\mathbf{b}))\frac{(\mathbf{u}-\mathbf{a} )^{\bm{\alpha}}(\mathbf{v}-\mathbf{b})^{\bm{\alpha}_{\prime}}}{\bm{\alpha}!\, \bm{\alpha}_{\prime}!}.$ $$

We have $|(\mathbf{u}-\mathbf{a})^{\bm{\alpha}}|=\prod_{j=1}^{r}|u_{j}-a_{j}|^{\alpha^{
(j)}}\leq|\mathbf{u}-\mathbf{a}|_{\infty}^{|\bm{\alpha}|}$ and $|(\mathbf{v}-\mathbf{b})^{\bm{\alpha}_{\prime}}|=\prod_{j=1}^{r_{\prime}}|v_{j
}-b_{j}|^{\alpha_{\prime}^{(j)}}\leq|\mathbf{v}-\mathbf{b}|_{\infty}^{|\bm{
\alpha}_{\prime}|}$. Consequently, for $f\in{\mathcal{H}}_{{r},{r_{\prime}}}^{\beta,\beta_{\prime}}\left([0,1]^{r},[0,
1]^{r_{\prime}},K\right)$,

$$ $\displaystyle\big{|}f(\mathbf{u},\mathbf{v})-P_{\mathbf{a},\mathbf{b}}^{\beta, \beta_{\prime}}f(\mathbf{u},\mathbf{v})\big{|}$ $\displaystyle\leq$ $\displaystyle\sum_{\begin{subarray}{c}\beta-1\leq|\bm{\alpha}|<\beta\\ \beta_{\prime}-1\leq|\bm{\alpha}_{\prime}|<\beta_{\prime}\end{subarray}}\left( \partial^{\bm{\alpha}+\bm{\alpha}_{\prime}}f(\mathbf{a}+\xi(\mathbf{u}-\mathbf {a}),\mathbf{b}+\xi(\mathbf{v}-\mathbf{b}))-\partial^{\bm{\alpha}+\bm{\alpha}_ {\prime}}f(\mathbf{a},\mathbf{b})\right)\frac{(\mathbf{u}-\mathbf{a})^{\bm{ \alpha}}(\mathbf{v}-\mathbf{b})^{\bm{\alpha}_{\prime}}}{\bm{\alpha}!\,\bm{ \alpha}_{\prime}!}$ (32) $\displaystyle\leq$ $\displaystyle K\left(|\mathbf{u}-\mathbf{a}|_{\infty}^{\beta}\vee|\mathbf{v}- \mathbf{b}|_{\infty}^{\beta_{\prime}}\right)$ $$

We may also write (31) as a linear combination of monomials

$$ $\displaystyle P_{\mathbf{a},\mathbf{b}}^{\beta,\beta_{\prime}}f(\mathbf{u}, \mathbf{v})=\sum_{\begin{subarray}{c}0\leq|\bm{\gamma}|<\beta\\ 0\leq|\bm{\gamma}_{\prime}|<\beta_{\prime}\end{subarray}}c_{\bm{\gamma},\bm{ \gamma}_{\prime}}\mathbf{u}^{\bm{\gamma}}\mathbf{v}^{\bm{\gamma}_{\prime}},$ (33) $$

for suitable coefficients $c_{\bm{\gamma},\bm{\gamma}_{\prime}}.$ For convenience, we omit the dependency on $\mathbf{a}$ and $\mathbf{b}$ in $c_{\bm{\gamma},\bm{\gamma}_{\prime}}.$ Since $\partial^{\bm{\gamma},\bm{\gamma}_{\prime}}P_{\mathbf{a},\mathbf{b}}^{\beta,
\beta_{\prime}}f(\mathbf{u},\mathbf{v})\,|_{(\mathbf{u}=0,\mathbf{v}=0)}=\bm{
\gamma}!\,\bm{\gamma}_{\prime}!\,c_{\bm{\gamma},\bm{\gamma}_{\prime}},$ we must have

$$ $\displaystyle c_{\bm{\gamma},\bm{\gamma}_{\prime}}=\sum_{\begin{subarray}{c} \bm{\gamma}\leq\bm{\alpha}\&|\bm{\alpha}|<\beta\\ \bm{\gamma}_{\prime}\leq\bm{\alpha}_{\prime}\&|\bm{\alpha}_{\prime}|<\beta_{ \prime}\end{subarray}}(\partial^{\bm{\alpha}+\bm{\alpha}_{\prime}}f)(\mathbf{a },\mathbf{b})\frac{(-\mathbf{a})^{\bm{\alpha}-\bm{\gamma}}\,(-\mathbf{b})^{\bm {\alpha}_{\prime}-\bm{\gamma}_{\prime}}}{\bm{\gamma}!\,\bm{\gamma}_{\prime}!\, (\bm{\alpha}-\bm{\gamma})!\,(\bm{\alpha}_{\prime}-\bm{\gamma}_{\prime})!}.$ $$

Notice that since $\mathbf{a}\in[0,1]^{r}$, $\mathbf{b}\in[0,1]^{r_{\prime}}$, and $f\in{\mathcal{H}}_{{r},{r_{\prime}}}^{\beta,\beta_{\prime}}\left([0,1]^{r},[0,
1]^{r_{\prime}},K\right)$,

$$ $\displaystyle|c_{\bm{\gamma}\bm{\gamma}_{\prime}}|\leq K/(\bm{\gamma}!\,\bm{ \gamma}_{\prime}!)\quad\text{and}\ \ \sum_{\begin{subarray}{c}\bm{\gamma}\geq{ 0}\\ \bm{\gamma}_{\prime}\geq{0}\end{subarray}}|c_{\bm{\gamma},\,\bm{\gamma}_{ \prime}}|\leq K\prod_{i=1}^{r}\prod_{j=1}^{r_{\prime}}\sum_{\gamma^{(i)}\geq 0 }\sum_{\gamma_{\prime}^{(j)}\geq 0}\frac{1}{\gamma^{(i)}!}\,\frac{1}{\gamma_{ \prime}^{(j)}!}=Ke^{r+r_{\prime}},$ (34) $$

where $\gamma=(\gamma^{(1)},\dots,\gamma^{(r)})$ and $\gamma_{\prime}=(\gamma_{\prime}^{(1)},\dots,\gamma_{\prime}^{(r_{\prime})})$.

Consider the set of grid points

$$ $\displaystyle\mathbf{D}(M):=\{\mathbf{u}_{{\bm{\ell}}^{(1)}}=$ $\displaystyle(\ell^{(1)}_{j}/M_{1})_{j=1,\ldots,r}\text{ and }\mathbf{v}_{{\bm {\ell}}^{(2)}}=(\ell^{(2)}_{j}/M_{2})_{j=1,\ldots,r_{\prime}}$ $\displaystyle:\,$ $\displaystyle{\bm{\ell}}^{(1)}=(\ell^{(1)}_{1},\ldots,\ell^{(1)}_{r})\in\{0,1, \ldots,M_{1}\}^{r},$ $\displaystyle{\bm{\ell}}^{(2)}=(\ell^{(2)}_{1},\ldots,\ell^{(2)}_{r})\in\{0,1, \ldots,M_{2}\}^{r_{\prime}},M_{1}=M^{\widetilde{\beta}/\beta},M_{2}=M^{ \widetilde{\beta}/\beta_{\prime}}\}.$ $$

The cardinality of this set is $(M_{1}+1)^{r}\cdot(M_{2}+1)^{r_{\prime}}.$ We write $\mathbf{u}_{{\bm{\ell}}^{(1)}}=(u^{(j)}_{{\bm{\ell}}^{(1)}})_{j=1,\ldots,r}$ and $\mathbf{v}_{{\bm{\ell}}^{(2)}}=(v^{(j)}_{{\bm{\ell}}^{(2)}})_{j=1,\ldots,r_{
\prime}}$ to denote the components of $\mathbf{u}_{{\bm{\ell}}^{(1)}}$ and $\mathbf{v}_{{\bm{\ell}}^{(2)}}$ respectively. With slight abuse of notation we denote $\mathbf{w}=(\mathbf{u},\mathbf{v})=(u^{(1)},\ldots,u^{(r)},v^{(1)},\ldots,v^{(
r_{\prime})})$, ${\bm{\ell}}=({\bm{\ell}}^{(1)},{\bm{\ell}}^{(2)})=(\ell^{(1)}_{1},\ldots,\ell^
{(1)}_{r},\ell^{(2)}_{1},\ldots,\ell^{(2)}_{r_{\prime}})$ and $\mathbf{w}_{\bm{\ell}}=(w^{(j)}_{\bm{\ell}})_{j=1,\ldots,r+r_{\prime}}=(
\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{v}_{{\bm{\ell}}^{(2)}})=(u^{(1)}_{{\bm{
\ell}}^{(1)}},\ldots,u^{(r)}_{{\bm{\ell}}^{(1)}},v^{(1)}_{{\bm{\ell}}^{(2)}},
\ldots,u^{(r_{\prime})}_{{\bm{\ell}}^{(2)}})$.
Define

$$ $\displaystyle P^{\beta,\beta_{\prime}}f(\mathbf{u},\mathbf{v})$ $\displaystyle=\,$ $\displaystyle P^{\beta,\beta_{\prime}}f(\mathbf{w})$ $\displaystyle:=$ $\displaystyle\sum_{\mathbf{w}_{{\bm{\ell}}}\in\mathbf{D}(M)}P_{\mathbf{w}_{{ \bm{\ell}}}}^{\beta,\beta_{\prime}}f(\mathbf{w})\prod_{j=1}^{r+r_{\prime}}(1-M _{j}|w^{(j)}-w_{{\bm{\ell}}}^{(j)}|)_{+}$ $\displaystyle=$ $\displaystyle\sum_{\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{v}_{{\bm{\ell}}^{(2) }}\in\mathbf{D}(M)}P_{\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{v}_{{\bm{\ell}}^{ (2)}}}^{\beta,\beta_{\prime}}f(\mathbf{u},\mathbf{v})\left(\prod_{j=1}^{r}(1-M _{1}|u^{(j)}-u_{{\bm{\ell}}^{(1)}}^{(j)}|)_{+}\right)\left(\prod_{j=1}^{r_{ \prime}}(1-M_{2}|v^{(j)}-v_{{\bm{\ell}}^{(2)}}^{(j)}|)_{+}\right),$ $$

where $M_{j}=M_{1}$ for $j=1,\ldots,r$ and $M_{j}=M_{2}$ for $j=r+1,\ldots,r+r_{\prime}$.

###### Lemma 5 .

If $f\in{\mathcal{H}}_{{r},{r_{\prime}}}^{\beta,\beta_{\prime}}\left([0,1]^{r},[0,
1]^{r_{\prime}},K\right)$, then $\|P^{\beta,\beta_{\prime}}f-f\|_{L^{\infty}[0,1]^{r+r_{\prime}}}\leq KM^{-
\widetilde{\beta}}.$

###### Proof.

Since for all $\mathbf{w}=(w^{(1)},\ldots,w^{(r+r_{\prime})})\in[0,1]^{r+r_{\prime}},$

$$ $\displaystyle\sum_{\mathbf{w}_{{\bm{\ell}}}\in\mathbf{D}(M)}\prod_{j=1}^{r+r_{ \prime}}(1-M_{j}|w^{(j)}-w_{{\bm{\ell}}}^{(j)}|)_{+}=\prod_{j=1}^{r+r_{\prime} }\sum_{\ell=0}^{M_{j}}(1-M_{j}|w^{(j)}-\ell/M_{j}|)_{+}=1,$ (35) $$

we have

$$ $\displaystyle f(\mathbf{w})$ $\displaystyle=f(\mathbf{u},\mathbf{v})$ $\displaystyle=\sum_{\begin{subarray}{c}\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{ v}_{{\bm{\ell}}^{(2)}}\in\mathbf{D}(M):\\ \|\mathbf{u}-\mathbf{u}_{{\bm{\ell}}^{(1)}}\|_{\infty}\leq 1/M_{1}\\ \|\mathbf{v}-\mathbf{v}_{{\bm{\ell}}^{(2)}}\|_{\infty}\leq 1/M_{2}\end{ subarray}}\,f(\mathbf{u},\mathbf{v})\left(\prod_{j=1}^{r}(1-M_{1}|u^{(j)}-u_{{ \bm{\ell}}^{(1)}}^{(j)}|)_{+}\right)\left(\prod_{j=1}^{r_{\prime}}(1-M_{2}|v^{ (j)}-v_{{\bm{\ell}}^{(2)}}^{(j)}|)_{+}\right)$ $$

and with (32),

$$ $\displaystyle\big{|}P^{\beta,\beta_{\prime}}f(\mathbf{u},\mathbf{v})-f(\mathbf {u},\mathbf{v})\big{|}$ $\displaystyle\leq\max_{\begin{subarray}{c}\mathbf{u}_{{\bm{\ell}}^{(1)}}, \mathbf{v}_{{\bm{\ell}}^{(2)}}\in\mathbf{D}(M):\\ \|\mathbf{u}-\mathbf{u}_{{\bm{\ell}}^{(1)}}\|_{\infty}\leq 1/M_{1}\\ \|\mathbf{v}-\mathbf{v}_{{\bm{\ell}}^{(2)}}\|_{\infty}\leq 1/M_{2}\end{ subarray}}\big{|}P_{\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{v}_{{\bm{\ell}}^{(2 )}}}^{\beta,\beta_{\prime}}f(\mathbf{u},\mathbf{v})-f(\mathbf{u},\mathbf{v}) \big{|}$ $\displaystyle\leq K\left(M_{1}^{-\beta}\vee M_{2}^{-\beta_{\prime}}\right)=K\, M^{-\widetilde{\beta}}.$ $$

∎

In the next few steps, we describe how to build a network that approximates $P^{\beta,\beta_{\prime}}f$.

###### Lemma 6 .

Let $M,m,$ be any positive integer. Denote $M_{1}=M^{\widetilde{\beta}/\beta}$, $M_{2}=M^{\widetilde{\beta}/\beta_{\prime}}$, $\mathsf{M}=(M_{1}+1)^{r}(M_{2}+1)^{r_{\prime}}$ and ${\mathsf{r}_{\mathrm{sum}}}=r+r_{\prime}$. Then there exists a network

$$ $\operatorname{Hat}^{{\mathsf{r}_{\mathrm{sum}}}}\in\mathcal{F}\left(2+(m+5) \lceil\log_{2}({\mathsf{r}_{\mathrm{sum}}})\rceil,{\mathsf{r}_{\mathrm{sum}}}, 2{\mathsf{r}_{\mathrm{sum}}}\mathsf{M},{\mathsf{r}_{\mathrm{sum}}}\mathsf{M},6 {\mathsf{r}_{\mathrm{sum}}}\mathsf{M},\ldots,6{\mathsf{r}_{\mathrm{sum}}} \mathsf{M},\mathsf{M}),s,1\right)$ $$

with $s\leq 37{\mathsf{r}_{\mathrm{sum}}}^{2}\mathsf{M}(m+5)\lceil\log_{2}({\mathsf{
r}_{\mathrm{sum}}})\rceil,$ such that $\operatorname{Hat}^{r}\in[0,1]^{\mathsf{M}}$ and for any $\mathbf{u}=(u^{(1)},\ldots,u^{(j)})\in[0,1]^{r}$ and for any $\mathbf{v}=(v^{(1)},\ldots,v^{(j)})\in[0,1]^{r_{\prime}}$

$$ $\displaystyle\Bigg{|}\operatorname{Hat}^{{\mathsf{r}_{\mathrm{sum}}}}(\mathbf{ u},\mathbf{v})-\Bigg{\{}$ $\displaystyle\Big{(}\prod_{j=1}^{r}(1/M_{1}-|u^{(j)}-u^{(j)}_{{\bm{\ell}}^{(1) }}|)_{+}\Big{)}\times$ $\displaystyle\Big{(}\prod_{j=1}^{r_{\prime}}(1/M_{2}-|v^{(j)}-v^{(j)}_{{\bm{ \ell}}^{(2)}}|)_{+}\Big{)}\Bigg{\}}_{\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{v} _{{\bm{\ell}}^{(2)}}\in\mathbf{D}(M)}\Bigg{|}_{\infty}\leq{\mathsf{r}_{\mathrm {sum}}}^{2}2^{-m}.$ $$

For any $\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{v}_{{\bm{\ell}}^{(2)}}\in\mathbf{D}(M),$ the support of the function $(\mathbf{u},\mathbf{v})\mapsto(\operatorname{Hat}^{r+r_{\prime}}(\mathbf{u},
\mathbf{v}))_{\mathbf{u}_{\ell^{(1)}},\mathbf{v}_{\ell^{(2)}}}$ is moreover contained in the support of the function

$$ $(\mathbf{u},\mathbf{v})\mapsto\left\{\bigg{(}\prod_{j=1}^{r}(1/M-|u^{(j)}-u^{( j)}_{{\bm{\ell}}^{(1)}}|)_{+}\bigg{)}\bigg{(}\prod_{j=1}^{r_{\prime}}(1/M-|v^{ (j)}-v^{(j)}_{{\bm{\ell}}^{(2)}}|)_{+}\bigg{)}\right\}.$ $$

###### Proof.

Step $1$: (For $r+r_{\prime}=1$) Without loss of generality we consider the case when $r=1$ and $r_{\prime}=0$.
We compute the functions $\{(u^{(j)}-\ell/M_{1})_{+}\}_{j=1,\ell=0}^{r,M_{1}}$ and $\{(\ell/M_{1}-u^{(j)})_{+}\}_{j=1,\ell=0}^{r,M_{1}}$ for the first hidden layer of the network. This requires $2r(M_{1}+1)$ units (nodes) and $2r(M_{1}+1)$ non-zero parameters.

For the second hidden layer we compute the functions $(1/M_{1}-|u^{(j)}-\ell/M_{1}|)_{+}=(1/M_{1}-(u^{(j)}-\ell/M_{1})_{+}-(\ell/M_{
1}-u^{(j)})_{+})_{+}$ using the output $(u^{(j)}-\ell/M_{1})_{+}$ and $(\ell/M_{1}-u^{(j)})_{+}$ from the output of the first hidden layer. This requires $r(M_{1}+1)+r_{\prime}(M_{2}+1)$ units (nodes) and $2r(M_{1}+1)$ non-zero parameters. This proves the result for the base case when $r+r_{\prime}=1$.

Step $2$:
For $r+r_{\prime}>1,$ we compose the obtained network with networks that approximately compute the following

$$ $\left\{\left(\prod_{j=1}^{r}(1/M_{1}-|u^{(j)}-u^{(j)}_{{\bm{\ell}}^{(1)}}|)_{+ }\right)\left(\prod_{j=1}^{r_{\prime}}(1/M_{2}-|v^{(j)}-v^{(j)}_{{\bm{\ell}}^{ (2)}}|)_{+}\right)\right\}_{\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{v}_{{\bm{ \ell}}^{(2)}}\in\mathbf{D}(M)}.$ $$

For fixed $\mathbf{u}_{{\bm{\ell}}^{(1)}}$ and $\mathbf{v}_{{\bm{\ell}}^{(2)}}$, and from the use of Lemma 8 there exist $\operatorname{Mult}_{m}^{r+r_{\prime}}$ networks in the class

$$ $\mathcal{F}\left(2+(m+5)\lceil\log_{2}(r+r_{\prime})\rceil,(r+r_{\prime},2(r+r _{\prime}),r+r_{\prime},6(r+r_{\prime}),6(r+r_{\prime}),\ldots,6(r+r_{\prime}) ,1)\right)$ $$

computing $(\prod_{j=1}^{r}(1/M_{1}-|u^{(j)}-u_{{\bm{\ell}}^{(1)}}|)_{+})\times(\prod_{j=
1}^{r_{\prime}}(1/M_{2}-|v^{(j)}-v_{{\bm{\ell}}^{(2)}}|)_{+})$ up to an error that is bounded by $(r+r_{\prime})^{2}\,2^{-m}$. Observe that we have two extra hidden layers to compute $(1/M_{1}-|u^{(j)}-u_{{\bm{\ell}}^{(1)}}|)_{+})$ and $(1/M_{2}-|v^{(j)}-v_{{\bm{\ell}}^{(2)}}|)_{+})$ for fixed $\mathbf{u}_{{\bm{\ell}}^{(1)}}$ and $\mathbf{v}_{{\bm{\ell}}^{(2)}}$ respectively, before we enter into the multinomial computation by regime invoking Lemma 8. Observe that the number of parameters in this network is upper bounded by $37(r+r_{\prime})^{2}(m+5)\lceil\log_{2}(r+r_{\prime})\rceil$.

Now we use the parallelization technique to have $(M_{1}+1)^{r}\cdot(M_{1}+1)^{r}$ parallel architecture for all elements of $\mathbf{D}(M)$. This provides the existence of the network with the number of non-zero parameters bounded by $37(r+r_{\prime})^{2}(M_{1}+1)^{r}(M_{2}+1)^{r_{\prime}}(m+5)\lceil\log_{2}(r+r
_{\prime})\rceil$

By Lemma 8, for any $\mathbf{x}\in{\mathbb{R}}^{r}$, $\operatorname{Mult}_{m}^{r}(\mathbf{x})=0$ if one of the components of $\mathbf{x}$ is zero. This shows that for any $\mathbf{u}_{\ell^{(1)}},\mathbf{v}_{\ell^{(2)}}\in\mathbf{D}(M)$, the support of the function $(\mathbf{u},\mathbf{v})\mapsto(\operatorname{Hat}^{r+r_{\prime}}(\mathbf{u},
\mathbf{v}))_{\mathbf{u}_{\ell^{(1)}},\mathbf{v}_{\ell^{(2)}}}$ is contained in the support of the function $(\mathbf{u},\mathbf{v})\mapsto\left(\prod_{j=1}^{r}(1/M-|u^{(j)}-u^{(j)}_{{\bm
{\ell}}^{(1)}}|)_{+}\prod_{j=1}^{r_{\prime}}(1/M-|v^{(j)}-v^{(j)}_{{\bm{\ell}}
^{(2)}}|)_{+}\right)$.

∎

###### Proof of Theorem 5 .

All the constructed networks in this proof are of the form $\mathcal{F}(L,\mathbf{p},s)=\mathcal{F}(L,\mathbf{p},s,\infty)$ with $F=\infty$. Denote $M_{1}=M^{\widetilde{\beta}/\beta}$, $M_{2}=M^{\widetilde{\beta}/\beta_{\prime}}$, ${\beta_{\mathrm{sum}}}=\beta+\beta_{\prime}$, and ${\mathsf{r}_{\mathrm{sum}}}=r+r_{\prime}$. Let $M$ be the largest integer such that $\mathsf{M}=(M_{1}+1)^{r}(M_{2}+1)^{r_{\prime}}\leq N$ and define $L^{*}:=(m+5)\lceil\log_{2}({\beta_{\mathrm{sum}}}\vee{\mathsf{r}_{\mathrm{sum}
}})\rceil$. Thanks to (34), (33) and Lemma 9, we can add one hidden layer to the network $\operatorname{Mon}_{m,{\beta_{\mathrm{sum}}}}^{\mathsf{r}_{\mathrm{sum}}}$ to obtain a network

$$ $Q_{1}\in\mathcal{F}\big{(}2+L^{*},(r,6\lceil\beta\rceil C_{{\mathsf{r}_{ \mathrm{sum}}},{\beta_{\mathrm{sum}}}},\ldots,6\lceil\beta\rceil C_{{\mathsf{r }_{\mathrm{sum}}},{\beta_{\mathrm{sum}}}},C_{{\mathsf{r}_{\mathrm{sum}}},{ \beta_{\mathrm{sum}}}},\mathsf{M})\big{)},$ $$

such that $Q_{1}(\mathbf{u},\mathbf{v})\in[0,1]^{\mathsf{M}}$ and for any $\mathbf{u}\in[0,1]^{r}$ and for any $\mathbf{v}\in[0,1]^{r_{\prime}}$

$$ $\displaystyle\Big{|}Q_{1}(\mathbf{u},\mathbf{v})-\Big{(}\frac{P^{\beta,\beta_{ \prime}}f(\mathbf{u},\mathbf{v})}{B}+\frac{1}{2}\Big{)}_{\mathbf{u}_{{\bm{\ell }}^{(1)}},\mathbf{v}_{{\bm{\ell}}^{(2)}}\in\mathbf{D}(M)}\Big{|}_{\infty}\leq \beta^{2}_{\mathrm{sum}}2^{-m}$ (36) $$

with $B:=\lceil 2Ke^{\mathsf{r}_{\mathrm{sum}}}\rceil$. The total number of non-zero parameters in the $Q_{1}$ network is $6{\mathsf{r}_{\mathrm{sum}}}({\beta_{\mathrm{sum}}}+1)C_{{\mathsf{r}_{\mathrm{
sum}}},{\beta_{\mathrm{sum}}}}+42({\beta_{\mathrm{sum}}}+1)^{2}C_{{\mathsf{r}_
{\mathrm{sum}}},{\beta_{\mathrm{sum}}}}^{2}(L^{*}+1)+C_{{\mathsf{r}_{\mathrm{
sum}}},{\beta_{\mathrm{sum}}}}\mathsf{M}$.

Recall that the network $\operatorname{Hat}^{\mathsf{r}_{\mathrm{sum}}}$ computes the products of hat functions (splines) $(\prod_{j=1}^{r}(1/M_{1}-|u^{(j)}-u_{{\bm{\ell}}^{(1)}}|)_{+})(\prod_{j=1}^{r_
{\prime}}(1/M_{2}-|v^{(j)}-v_{{\bm{\ell}}^{(2)}}|)_{+})$ up to an error that is bounded by $\mathsf{r}^{2}_{\mathrm{sum}}2^{-m}.$ It requires at most $37\mathsf{r}^{2}_{\mathrm{sum}}NL^{\ast}$ active parameters. Observe that $C_{{\mathsf{r}_{\mathrm{sum}}},{\beta_{\mathrm{sum}}}}\leq({\beta_{\mathrm{sum
}}}+1)^{\mathsf{r}_{\mathrm{sum}}}\leq N$ by the definition of $C_{r,\beta}$ and the assumptions on $N.$ By Lemma 6, the networks $Q_{1}$ and $\operatorname{Hat}^{\mathsf{r}_{\mathrm{sum}}}$ can be embedded into a joint parallel network $(Q_{1},\operatorname{Hat}^{\mathsf{r}_{\mathrm{sum}}})$ with $2+L^{*}$ hidden layers of size $({\mathsf{r}_{\mathrm{sum}}},6({\mathsf{r}_{\mathrm{sum}}}+\lceil{\beta_{
\mathrm{sum}}}\rceil)N,\ldots,6({\mathsf{r}_{\mathrm{sum}}}+\lceil{\beta_{
\mathrm{sum}}}\rceil)N,2\mathsf{M})$. Using $C_{r,\beta}\vee(M+1)^{r}\leq N$ again, the number of non-zero parameters in the combined network $(Q_{1},\operatorname{Hat}^{r})$ is bounded by

$$ $\displaystyle\begin{split}&6{\mathsf{r}_{\mathrm{sum}}}({\beta_{\mathrm{sum}}} +1)C_{{\mathsf{r}_{\mathrm{sum}}},{\beta_{\mathrm{sum}}}}+42({\beta_{\mathrm{ sum}}}+1)^{2}C_{{\mathsf{r}_{\mathrm{sum}}},{\beta_{\mathrm{sum}}}}^{2}(L^{*}+ 1)+C_{{\mathsf{r}_{\mathrm{sum}}},{\beta_{\mathrm{sum}}}}\mathsf{M}+37\mathsf{ r}^{2}_{\mathrm{sum}}NL^{\ast}\\ &\leq 42({\mathsf{r}_{\mathrm{sum}}}+{\beta_{\mathrm{sum}}}+1)^{2}C_{{\mathsf{ r}_{\mathrm{sum}}},{\beta_{\mathrm{sum}}}}N(1+L^{*})\\ &\leq 84({\mathsf{r}_{\mathrm{sum}}}+{\beta_{\mathrm{sum}}}+1)^{3+{\mathsf{r}_ {\mathrm{sum}}}}N(m+5),\end{split}$ (37) $$

where for the last inequality, we used $C_{{\mathsf{r}_{\mathrm{sum}}},{\beta_{\mathrm{sum}}}}\leq({\beta_{\mathrm{sum
}}}+1)^{\mathsf{r}_{\mathrm{sum}}},$ the definition of $L^{*}$ and that for any $x\geq 1,$ $1+\lceil\log_{2}(x)\rceil\leq 2+\log_{2}(x)\leq 2(1+\log(x))\leq 2x.$

Next, we pair the $(\mathbf{u}_{{\bm{\ell}}^{(1)}},\mathbf{v}_{{\bm{\ell}}^{(2)}})$-th entry of the output of $Q_{1}$ and $\operatorname{Hat}^{r}$ and apply to each of the $\mathsf{M}$ pairs the $\operatorname{Mult}_{m}$ network described in Lemma 7. In the last layer, we add all entries. By Lemma 7 this requires at most $24(m+5)\mathsf{M}+\mathsf{M}\leq 25(m+5)N$ active parameters for the $\mathsf{M}$ multiplications and the sum. Using Lemma 7, Lemma 6, (36) and triangle inequality, there exists a network $Q_{2}\in\mathcal{F}(2+L^{*}+m+6,({\mathsf{r}_{\mathrm{sum}}},6({\mathsf{r}_{
\mathrm{sum}}}+\lceil{\beta_{\mathrm{sum}}}\rceil)N,\ldots,6({\mathsf{r}_{
\mathrm{sum}}}+\lceil{\beta_{\mathrm{sum}}}\rceil)N,1))$ such that for any $\mathbf{u}\in[0,1]^{r}$ and for any $\mathbf{v}\in[0,1]^{r_{\prime}}$

$$ $\displaystyle\Bigg{|}Q_{2}(\mathbf{u},\mathbf{v})-\sum_{\mathbf{u}_{{\bm{\ell} }^{(1)}},\mathbf{v}_{{\bm{\ell}}^{(2)}}\in\mathbf{D}(M)}\Big{(}\frac{P^{\beta, \beta_{\prime}}f(\mathbf{u},\mathbf{v})}{B}+\frac{1}{2}\Big{)}$ $\displaystyle\Big{(}\prod_{j=1}^{r}(1/M_{1}-|u^{(j)}-u^{(j)}_{{\bm{\ell}}^{(1) }}|)_{+}\Big{)}$ $\displaystyle\Big{(}\prod_{j=1}^{r_{\prime}}(1/M_{2}-|v^{(j)}-v^{(j)}_{{\bm{ \ell}}^{(2)}}|)_{+}\Big{)}\Bigg{|}$ $$

$$ $\displaystyle\leq\sum_{\begin{subarray}{c}\mathbf{u}_{{\bm{\ell}}^{(1)}}, \mathbf{v}_{{\bm{\ell}}^{(2)}}\in\mathbf{D}(M):\\ \|\mathbf{u}-\mathbf{u}_{{\bm{\ell}}^{(1)}}\|_{\infty}\leq 1/M_{1}\\ \|\mathbf{v}-\mathbf{v}_{{\bm{\ell}}^{(2)}}\|_{\infty}\leq 1/M_{2}\end{ subarray}}(1+\mathsf{r}^{2}_{\mathrm{sum}}+\beta^{2}_{\mathrm{sum}})2^{-m}$ $\displaystyle\leq(1+\mathsf{r}^{2}_{\mathrm{sum}}+\beta^{2}_{\mathrm{sum}})2^{ r-m}.$ (38) $$

Here, the first inequality follows from the fact that the support of $(\operatorname{Hat}^{r+r_{\prime}}(\mathbf{u},\mathbf{v}))_{\mathbf{u}_{\ell^{
(1)}},\mathbf{v}_{\ell^{(2)}}}$ is contained in the support of $\left(\prod_{j=1}^{r}(1/M-|u^{(j)}-u^{(j)}_{{\bm{\ell}}^{(1)}}|)_{+}\prod_{j=1
}^{r_{\prime}}(1/M-|v^{(j)}-v^{(j)}_{{\bm{\ell}}^{(2)}}|)_{+}\right)$ (see Lemma 6). Because of (37), the network $Q_{2}$ has at most

$$ $\displaystyle\begin{split}&109({\mathsf{r}_{\mathrm{sum}}}+{\beta_{\mathrm{sum }}}+1)^{3+{\mathsf{r}_{\mathrm{sum}}}}N(m+5)\end{split}$ (39) $$

non-zero parameters.

To obtain a network reconstruction of the function $f$, it remains to scale and shift the output entries. This is not entirely trivial because of the bounded parameter weights in the network. Recall that $B=\lceil 2Ke^{r}\rceil.$ The network $x\mapsto BM_{1}^{r}M_{2}^{r_{\prime}}x$ is in the class $\mathcal{F}(3,(1,M_{1}^{r}M_{2}^{r_{\prime}},1,\lceil 2Ke^{r}\rceil,1))$ with shift vectors $\mathbf{v}_{j}$ are all equal to zero and weight matrices $W_{j}$ with all entries equal to one. Because of $N\geq(K+1)e^{\mathsf{r}_{\mathrm{sum}}},$ the number of parameters of this network is bounded by $2M_{1}^{r}M_{2}^{r_{\prime}}+2\lceil 2Ke^{r}\rceil\leq 6N$. This shows existence of a network in the class $\mathcal{F}(4,(1,2,2M_{1}^{r}M_{2}^{r_{\prime}},2,2\lceil 2Ke^{r}\rceil,1))$ computing $a\mapsto BM_{1}^{r}M_{2}^{r_{\prime}}(a-c)$ with $c:=1/(2M_{1}^{r}M_{2}^{r_{\prime}}).$ This network computes in the first hidden layer $(a-c)_{+}$ and $(c-a)_{+}$ and then applies the network $x\mapsto BM_{1}^{r}M_{2}^{r_{\prime}}x$ to both units. In the output layer, the second value is subtracted from the first one. This requires at most $6+12N$ active parameters.

Because of (38) and (35), there exists a network $Q_{3}$ in

$$ $\displaystyle\mathcal{F}\big{(}(m+13)+L^{*},({\mathsf{r}_{\mathrm{sum}}},6({ \mathsf{r}_{\mathrm{sum}}}+\lceil{\beta_{\mathrm{sum}}}\rceil)N,\ldots,6({ \mathsf{r}_{\mathrm{sum}}}+\lceil{\beta_{\mathrm{sum}}}\rceil)N,1)\big{)}$ $$

such that

$$ $\displaystyle\Bigg{|}Q_{3}(\mathbf{u},\mathbf{v})-\sum_{\mathbf{u}_{{\bm{\ell} }^{(1)}},\mathbf{v}_{{\bm{\ell}}^{(2)}}\in\mathbf{D}(M)}P^{\beta,\beta_{\prime }}f(\mathbf{u},\mathbf{v})$ $\displaystyle\Big{(}\prod_{j=1}^{r}(1/M_{1}-|u^{(j)}-u^{(j)}_{{\bm{\ell}}^{(1) }}|)_{+}\Big{)}$ $\displaystyle\Big{(}\prod_{j=1}^{r_{\prime}}(1/M_{2}-|v^{(j)}-v^{(j)}_{{\bm{ \ell}}^{(2)}}|)_{+}\Big{)}\Bigg{|}$ $$

$$ $\displaystyle\leq(2K+1)M_{1}^{r}M_{2}^{r_{\prime}}(1+\mathsf{r}^{2}_{\mathrm{ sum}}+\beta^{2}_{\mathrm{sum}})(2e)^{\mathsf{r}_{\mathrm{sum}}}2^{-m},\ \ \text{for all}\ (\mathbf{u},\mathbf{v})\in[0,1]^{\mathsf{r}_{\mathrm{sum}}}.$ $$

With (39), the number of non-zero parameters of $Q_{3}$ is bounded by

$$ $\displaystyle 109({\mathsf{r}_{\mathrm{sum}}}+{\beta_{\mathrm{sum}}}+1)^{3+{ \mathsf{r}_{\mathrm{sum}}}}N(m+6).$ $$

Observe that by construction $\mathsf{M}=(M_{1}+1)^{r}(M_{2}+1)^{r_{\prime}}\leq N\leq(3M_{1})^{r}(3M_{2})^{
r_{\prime}}=3^{\mathsf{r}_{\mathrm{sum}}}M^{\widetilde{r}}$ and hence $M^{-\widetilde{\beta}}\leq N^{-\widetilde{\beta}/\widetilde{r}}3^{{\mathsf{r}_
{\mathrm{sum}}}\widetilde{\beta}/\widetilde{r}}.$ Together with Lemma 5, the result follows.
∎

###### Theorem 5 .

###### Lemma 5 .

###### Proof.

###### Lemma 6 .

###### Proof.

###### Proof of Theorem 5 .

### K.1 Embedding properties of neural network function classes

We denote $\mathcal{F}(L,\bm{p})$ as the class of neural networks with $L$ hidden layers and $\bm{p}\in\mathbb{N}^{L+2}$ nodes per layer. The class $\mathcal{F}(L,\bm{p})$ is subset of $\mathcal{F}(L,\bm{p})$ with the sparsity parameter $s$.

For the approximation of a function by a network, we first construct smaller networks computing simpler objects. Let $\mathbf{p}=(p_{0},\ldots,p_{L+1})$ and $\mathbf{p}^{\prime}=(p_{0}^{\prime},\ldots,p_{L+1}^{\prime}).$ To combine networks, we make frequent use of the following rules.

Enlarging: $\mathcal{F}(L,\mathbf{p},s)\subseteq\mathcal{F}(L,\mathbf{q},s^{\prime})$ whenever $\mathbf{p}\leq\mathbf{q}$ componentwise and $s\leq s^{\prime}.$

Composition: Suppose that $f\in\mathcal{F}(L,\mathbf{p})$ and $g\in\mathcal{F}(L^{\prime},\mathbf{p}^{\prime})$ with $p_{L+1}=p_{0}^{\prime}.$ For a vector $\mathbf{v}\in\mathbb{R}^{p_{L+1}}$ we define the composed network $g\circ\sigma_{\mathbf{v}}(f)$ which is in the space $\mathcal{F}(L+L^{\prime}+1,(\mathbf{p},p_{1}^{\prime},\ldots,p_{L^{\prime}+1}^
{\prime})).$ In most of the cases that we consider, the output of the first network is non-negative and the shift vector $\mathbf{v}$ will be taken to be zero.

Additional layers/depth synchronization: To synchronize the number of hidden layers for two networks, we can add additional layers with an identity weight matrix, such that

$$ $\displaystyle\mathcal{F}(L,\mathbf{p},s)\subset\mathcal{F}(L+q,(\underbrace{p_ {0},\ldots,p_{0}}_{q\text{ times}},\mathbf{p}),s+qp_{0}).$ (40) $$

Parallelization: Suppose that $f,g$ are two networks with the same number of hidden layers and the same input dimension, that is, $f\in\mathcal{F}(L,\mathbf{p})$ and $g\in\mathcal{F}(L,\mathbf{p}^{\prime})$ with $p_{0}=p_{0}^{\prime}.$ The parallelized network $(f,g)$ computes $f$ and $g$ simultaneously in a joint network in the class $\mathcal{F}(L,(p_{0},p_{1}+p_{1}^{\prime},\ldots,p_{L+1}+p_{L+1}^{\prime})).$

### K.2 Technical lemmas for the proof of Theorem 5

We use $\mathcal{F}(L,{\bf{r}})$ to denote a fully connected network with $L$ deep layers and ${\bf r}\in\mathbb{N}_{0}^{L+2}$ representing the nodes in each layer.

The following technical lemmas are required for the proof of Theorem 5. Lemma 7, Lemma 8, and Lemma 9 restate Lemma A.2, Lemma A.3, and Lemma A.4 from Schmidt-Hieber, (), respectively.

###### Lemma 7 .

For any positive integer $m,$ there exists a network $\operatorname{Mult}_{m}\in\mathcal{F}(m+4,(2,6,6,\ldots,6,1)),$ such that $\operatorname{Mult}_{m}(x,y)\in[0,1],$

$$ $\displaystyle\big{|}\operatorname{Mult}_{m}(x,y)-xy\big{|}\leq 2^{-m},\quad \text{for all}\ x,y\in[0,1],$ $$

and $\operatorname{Mult}_{m}(0,y)=\operatorname{Mult}_{m}(x,0)=0.$

###### Lemma 8 .

For any positive integer $m,$ there exists a network

$$ $\operatorname{Mult}_{m}^{r}\in\mathcal{F}((m+5)\lceil\log_{2}r\rceil,(r,6r,6r, \ldots,6r,1))$ $$

such that $\operatorname{Mult}_{m}^{r}\in[0,1]$ and

$$ $\displaystyle\Big{|}\operatorname{Mult}_{m}^{r}(\mathbf{x})-\prod_{i=1}^{r}x_{ i}\Big{|}\leq r^{2}2^{-m},\quad\text{for all}\ \ \mathbf{x}=(x_{1},\ldots,x_{r })\in[0,1]^{r}.$ $$

Moreover, $\operatorname{Mult}_{m}^{r}(\mathbf{x})=0$ if one of the components of $\mathbf{x}$ is zero.

The number of monomials with degree $|\bm{\alpha}|<\gamma$ is denoted by $C_{r,\gamma}.$ Obviously, $C_{r,\gamma}\leq(\gamma+1)^{r}$ since each $\alpha_{i}$ has to take values in $\{0,1,\ldots,\lfloor\gamma\rfloor\}.$

###### Lemma 9 .

For $\gamma>0$ and any positive integer $m,$ there exists a network

$$ $\displaystyle\operatorname{Mon}_{m,\gamma}^{r}\in\mathcal{F}\big{(}1+(m+5) \lceil\log_{2}(\gamma\vee 1)\rceil,(r,6\lceil\gamma\rceil C_{r,\gamma},\ldots, 6\lceil\gamma\rceil C_{r,\gamma},C_{r,\gamma})\big{)},$ $$

such that $\operatorname{Mon}_{m,\gamma}^{r}\in[0,1]^{C_{r,\gamma}}$ and

$$ $\displaystyle\Big{|}\operatorname{Mon}_{m,\gamma}^{r}(\mathbf{x})-(\mathbf{x}^ {\bm{\alpha}})_{|\bm{\alpha}|<\gamma}\Big{|}_{\infty}\leq\gamma^{2}2^{-m}, \quad\text{for all}\ \mathbf{x}\in[0,1]^{r}.$ $$

###### Lemma 7 .

###### Lemma 8 .

###### Lemma 9 .

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 4
questions: 1. In regards to the discussion in 2.4.3 and Appendix C with the multiple manifold case: If a partition of unity is used to define local transport maps, we require that the local distributions on each chart need to have density lower bounded. However, a partition of unity will have its weights vanish to 0 on the boundary, so the induced local distribution would have density that vanishes to 0 on the boundary. How can you apply optimal transport regularity theory in that case, when your target has density that is not lower bounded?
soundness: 4
strengths: 1. Very well written paper
2. Method is clearly explained and proof are followable
summary: A statistical analysis of conditional deep generative models using a likelihood approach (specifically a sieve MLE) is proposed and analyzed. A convergence rate for the estimator is established using neural network approximation theory for the (conditional) generator. Numerical examples are provided showing the efficacy of the method.
weaknesses: 1. The method seems to be extending the techniques of previous work and filling in some details with approximation.
2. Question on the applicability to multiple chart manifolds.
3. Numerical examples are very simple.

---

confidence: 4
contribution: 3
initial_rating: 3
presentation: 1
questions: Do numerical results support the rates proved in theory? Maybe we can take logarithm on the sample size and inspect the slope of the error.

Assumption 3 excludes interior points. Does this assumption hold (approximately) in practical applications?
soundness: 3
strengths: The theoretical results appear to be correct and sound. Moreover, there are detailed analyses for different setups and matching lower bound.

There are numerical results to support the performance of neural network generative models for distribution regression.
summary: This paper studies using neural network generative models for distribution regression. Statistical rates are established, and a significant effort is devoted to resolving the curse of dimensionality issue. Numerical results are provided to support the theory.
weaknesses: The paper is not very easy to follow in places, mostly due to the technical presentations and pure statistical discussions. I would suggest authors to prepare a preliminary section before mixing background introduction with main results in current Section 2.

It is relatively difficult for me to judge the technical novelty of the paper. On the one hand, the distribution regression problem and low-dimensional generative neural networks seem to pose technical difficulties. On the other hand, the results as well as most of the assumptions are replication of many existing works, such as Schmidt-Hieber, 2020 and Tang and Yang, 2023.

Practical impact of this work may be limited, due to many technical assumptions.

Numerical results loosely connect to the theoretical results.

The paper is not prepared in the standard ICLR style --- the left and right margins are altered significantly. Under the standard style template, this submission will go beyond the page limit.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 2
questions: How easy are some of the conditions to verify for the convergence theorems to hold?

How do you interpret the constants in the various theoretical statements? For example, some are dependent on the domain or other parameters, in what sense do they grow for perturbations of their dependencies? Can these be estimated in practice?
soundness: 2
strengths: The proposed analysis is helpful in understanding why generative models work as they do. Additionally the proposed error analysis is useful for anticipating the number of data points needed to achieve a specific targeted error with respect to the Wasserstein or Hellinger distances. 

The empirical results support the main claims of the work to some extent.
summary: The authors describe a method for estimating high dimensional conditional distributions that concentrate around a manifold. Such distributions appear in common settings such as conditional generative modeling. The authors consider the convergence of an estimator that improves based on the number of data points to the true distribution in the Hellinger distance sense. The authors then consider applying this estimator to some synthetic and real datasets to understand its empirical performance. The analysis considers two types of neural networks: a sparse neural network and fully connected one that parameterize the push forward of the latent variable to the conditional outcome. The analysis is heavily dependent on empirical process theory to establish different properties of the architectures.
weaknesses: Some of the theoretical results seem hard to verify in practice. It would be interesting for some contrived examples if the authors can show the convergence of the distances as a function of the number of points. 

The theorems contain some constants that do not have a clear way of being estimated making some of the errors hard to use in practice. 

Note that the paper does not appear to be in the right ICLR format since the margins are considerably smaller.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: see the weaknesses, they are a bit intertwined.

POST REBUTTAL: Concerns have been mostly addressed. raising score from 6 to 8.
soundness: 3
strengths: The paper tackles a very interesting topic. There has not been much theory of conditional generative models, in particular not on their statistical rates with respect to the true smoothness and intrinsic dimension. 

Theorem 1 presents a very nice approach on how to tackle these problems, which allows for splitting these problems up into the covering number and the accuracy of the approximation. This gives a nice guideline: Pick a function class that has a good covering scaling and good approximation capabilities. In particular the noise injection strategy seems very neat in corollary 2, if the data distribution is close to a singular one. Generally, the writing and the maths is very rigoruous, with clearly defined objects up to a few exceptions (see weaknesses). It is also nice that the qualitiative behaviour of the theorems is verified in understandable toy examples.
summary: This paper aims at providing statistical bounds for conditional generative models under various smoothness assumptions. The overall goal is to model conditional distributions $Y|X$ via conditional generative models, i.e., $y = G(x)_{*}P_Z$ of some latent distribution $P_Z$. First empirical convergence rates of so-called sieve MLE estimators are provided (Theorem 1), showing that as paired data tends to infinity the approximation of the measure tends to zero in probability.  Then also the Wasserstein distance between the pushforward distribution of the MLE estimate and the true generated distribution can be controlled via a pertubation strategy that depends on the data dimension and smoothness. Then with some ansatz function class spaces for the true generating process, a convergence theorem for sieve MLE with generating functions in some sparse NN class is presented. 

Finally, the sieve MLE is validated on some toy examples which is learned via a VAE and on the MNIST dataset.
weaknesses: First, let me say, that I come from a background in conditional generative models, not so much in statistical theory. Therefore I had difficulties following some claims in the paper. 

1) First, when I think of where conditional generative models are actually used, I am confused about the whole modelling assumption in section 2.1. Usually the latent space Z should take care of the uncertainty in Y, so that T(x,Z) models this perfectly. of course for trained models there will still be a mismatch, but treating this as an independent Gaussian seems very unrealistic. Can you elaborate this assumption maybe taking the example of inverse problems?

Can you maybe find for lets say a normalizing flow always a similar regular normalizing flow with a small gaussian noise added to it such that it becomes the "real" generative model? This would increase my confidence in this assumption.

2) It would be good to have a running example, something super simple like taking the pushforward of a Gaussian under a linear operator (with non trivial rank) to explain the rates and corollaries. 

3) I am worried about the scaling in corollary 1 for the $\sigma^{*}$, basically this means that if the noise level tends to zero, that the approximation accuracy has to increase by a lot. This might be due to the choice of the Hellinger distance. Why not try to control the Wasserstein distance in equation 7 in the first space? Since $W_1(p, p \star \mathcal{N}(0,\sigma))$ are close for small $\sigma$ this might alleviate this problem.

4) I feel like section 2.4 is not rigorous in the sense that it is very unclear if the article looks at joint measures of $(X,Y)$ or conditional measures $Y|X=x$. Looking at the proofs I think it holds for the condiitonals and arbitrary $X=x$? ´

5) For the experiments part, it would be very enlightening to consider manifolds with a known dimension (and maybe make them conditional) to validate how the intrinsic dimension comes into play [1]. Further, for the MNIST example, one does not need to take the MNIST dataset as the ground truth, but rather one can take a pretrained generator as the ground truth and then validate how smoothness and approximation error impact the Wasserstein/FID/... between the pretrained generator and a learned model. This would be much more in spirit with validating the bounds, since the point of this paper is not to come up with a practical algorithm I presume. 

6) Can the in probability bounds be strengethened to expectation bounds over x between the conditional measures under stronger assumptions, see also [4].Conditional Wasserstein Distances with Applications in Bayesian OT Flow Matching

7) Also, one algorithm where the composition approach might come into play is flow matching/diffusion with odes [2]. Here, if one takes a simple Euler scheme to sample one indeed has such a simple compositional structure. I would find it very interesting to briefly state the relation to those, as they are bit more SOTA and practical interest than VAEs at the moment. 

8) Can you make a practical example where the smoothness disparity is obvious?

9) This paper is missing some citations on generative models and manifolds [1,3], and the relation between conditional measures and joint measures [4] and also in what sense conditional generative models optimize the MLE estimator (somewhat clear for NFs and VAEs, less clear for diffusion/...). Here I am thinking along the lines of [5]

Overall, I am mainly missing practical takeaways and a running example where I can see that the assumptions made are indeed practical. Indeed, also the modelling assumptions remain unclear, and would greatly benefit from more explanations. I feel like the whole paper can be a bit streamlined to make it more digestible for the learning community as opposed to the statistics community, but overall I appreciate the deep theoretical analysis if my points can be addressed. 

[1] Diffusion Models Encode the Intrinsic Dimension of Data Manifolds, ICML 

[2] Flow Matching for Generative Modeling, ICLR 

[3] Score-Based Generative Models Detect Manifolds, NeurIPS

[4] Conditional Wasserstein Distances with Applications in Bayesian OT Flow Matching, arXiv

[5] Maximum Likelihood Training of Score-Based Diffusion Models, NeurIPS

Meta-review: This paper presents a theoretical analysis of conditional deep generative models, focusing on estimating high-dimensional conditional distributions that concentrate around manifolds. The authors provide statistical guarantees and convergence rates for different network architectures and present empirical results to support their theoretical claims.

The main strength of the paper lies in its attempt to provide a theoretical foundation for conditional generative models, offering convergence rates and error bounds for different network architectures. However, the paper suffers from several weaknesses. The practical implications and immediate impact on deep learning applications are limited. Many of the assumptions and constants in the theoretical framework are difficult to verify or estimate in practice. Moreover, the theoretical contribution appears to primarily combine existing work without offering significant new insights.

The decision to reject is primarily based on the limited novelty and practical impact of the work. While the paper provides a theoretical analysis, it does not offer substantial new insights or immediate practical applications in deep learning. The theoretical results largely align with expected outcomes and do not demonstrate significant deviations from existing work. Additionally, the initial submission's violation of formatting guidelines raises concerns about adherence to conference standards.

Decision: Reject