# Wasserstein Generative Regression

Shanshan Song Tong Wang∗ Guohao Shen Yuanyuan Lin ∗∗ Jian Huang Shangshan Song
and Tong Wang contributed equally to this work.Department of Statistics, The
Chinese University of Hong Kong, Hong Kong SAR, ChinaDepartment of Statistics,
The Chinese University of Hong Kong, Hong Kong SAR, ChinaDepartment of Applied
Mathematics, The Hong Kong Polytechnic University, Hong Kong SAR,
ChinaDepartment of Statistics, The Chinese University of Hong Kong, Hong Kong
SAR, ChinaDepartment of Applied Mathematics, The Hong Kong Polytechnic
University, Hong Kong SAR, ChinaCorresponding authors, Yuanyuan Lin, email:
<EMAIL_ADDRESS>Jian Huang, email<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a new and unified approach for nonparametric
regression and conditional distribution learning. Our approach simultaneously
estimates a regression function and a conditional generator using a generative
learning framework, where a conditional generator is a function that can
generate samples from a conditional distribution. The main idea is to estimate
a conditional generator that satisfies the constraint that it produces a good
regression function estimator. We use deep neural networks to model the
conditional generator. Our approach can handle problems with multivariate
outcomes and covariates, and can be used to construct prediction intervals. We
provide theoretical guarantees by deriving non-asymptotic error bounds and the
distributional consistency of our approach under suitable assumptions. We also
perform numerical experiments with simulated and real data to demonstrate the
effectiveness and superiority of our approach over some existing approaches in
various scenarios.

keywords: Conditional distribution, deep neural networks, generative learning,
nonparametric regression, non-asymptotic error bounds.

##  1 Introduction

Regression models and conditional distributions play a key role in a variety
of prediction and inference problems in statistics. There is a vast literature
on nonparametric methods for regression analysis and conditional density
estimation. Most existing methods use smoothing and basis expansion
techniques, including kernel smoothing, local polynomials, and splines
(Silverman,, 1986; Scott,, 1992; Fan and Gijbels,, 1996; Györfi et al.,, 2002;
Wasserman,, 2006; Tsybakov,, 2008). However, the existing nonparametric
regression and conditional density estimation methods suffer from the “curse
of dimensionality”, that is, their performance deteriorates dramatically as
the dimensionality of data increases. Indeed, most existing methods can only
effectively handle up to a few predictors. Moreover, most existing methods
only consider the case when the response is a scalar, but are not applicable
to the settings with a high-dimensional response vector.

To circumvent the curse of dimensionality, many researchers have proposed and
studied non- and semi-parametric models that impose certain structural
constraints that reduce the model dimensionality. Some notable examples
include the single index model (Ichimura,, 1993; Hardle et al.,, 1993), the
generalized additive model (Hastie and Tibshirani,, 1986; Stone,, 1986), and
the projection pursuit model (Friedman and Stuetzle,, 1981), among others.
However, these methods make strong assumptions about the model structure,
which may not hold in reality. Moreover, these methods aim at estimating the
regression function, but do not learn the conditional distribution. Therefore,
they can only provide point prediction, not interval prediction with a measure
of uncertainty.

In recent years, there have been many important developments in deep
generative learning (Salakhutdinov,, 2015), in which deep neural networks are
used to approximate high-dimensional functions, such as generator and
discriminator functions. In particular, for learning distributions of high-
dimensional data arising in image analysis and natural language processing,
the generative adversarial networks (GANs) (Goodfellow et al.,, 2014; Arjovsky
et al.,, 2017) have proven to be effective and achieved impressive success
(Reed et al.,, 2016; Zhu et al.,, 2017). Instead of estimating the functional
form of a density function, GANs start from a known reference distribution and
learn a map that pushes the reference distribution to the data distribution.
GANs have also been extended to learn conditional distributions (Mirza and
Osindero,, 2014; Kovachki et al.,, 2021; Zhou et al.,, 2022; Liu et al.,,
2021).

One of the main challenges in nonparametric regression is to estimate a
function that can accurately capture the relationship between covariate and
response variables. GANs can learn complex distributions. However, to the best
of our knowledge, there have not been systematic studies on how GANs can be
used for nonparametric regression, despite their successes in distribution
learning. Furthermore, conditional GANs, which are a natural extension of GANs
for learning conditional distributions, do not automatically guarantee a good
estimation of a regression function.

We propose a new and unified approach for nonparametric regression and
conditional distribution estimation. Our approach estimates the regression
function and a conditional generator at the same time using a generative
learning method. A conditional generator is a function that transforms a
random vector from a known reference distribution to the response variable
space, which can be used to sample from a conditional distribution. Thus, when
a conditional generator is estimated, it can be used to explore the target
conditional distribution. Theoretically, the regression function is the
expectation of the conditional generator with respect to the reference
distribution. However, empirically such an expectation may not produce a good
estimator of the regression function.

Our main idea is to constrain the conditional generator to produce samples
that minimize the quadratic loss of the regression function, which is computed
as the expectation of the conditional generator. Specifically, in the
objective function for estimating the conditional generator based on
distribution matching using the Wasserstein distance, we incorporate a
quadratic loss term to control the error of the estimated regression function.
We use deep neural networks to approximate the conditional generator, which
can capture the complex structure of the data distribution. In principle,
other approximation methods such as splines can also be used. However, deep
neural network approximation has the important advantage of being able to
adapting to the latent structure of the data distribution. For simplicity, we
call our method Wasserstein generative regression (WGR).

The proposed method has several attractive properties. First, it is applicable
to problems with a high-dimensional response variable, while the existing
methods typically only consider the case of a scalar response. Second, the
proposed method allows continuous, discrete and mixed types of predictors and
responses, while the smoothing and basis expansion methods are mainly
applicable to continuous-type variables. Third, since the proposed method
learns a conditional distribution generator, it can be used for constructing
prediction intervals. In comparison, the existing nonparametric regression can
only give point prediction. Finally, the proposed method is able to adapt to
the latent data structure in a data-driven manner and thus can mitigate the
curse of dimensionality, under the assumption the data distribution is
supported on an approximate low-dimensional set.

The rest of the paper is organized as follows. In Section 2 we describe the
proposed WGR method. We present the implementation details in Section 3. In
Section 4 we establish non-asymptotic error bounds for the proposed estimator
and show that it is consistent. In Section 5 we conduct numerical experiments,
including simulation studies and real data analysis, to evaluate the
performance of the proposed method. Technical proofs and additional numerical
experiments are given in Appendix.

##  2 Method

Consider a pair of random vectors $(X,Y)\in\mathcal{X}\times\mathcal{Y}$,
where $X$ is a vector of predictors and $Y$ is a vector of response variables.
Suppose $\mathcal{X}\subseteq\mathbb{R}^{d}$ and
$\mathcal{Y}\subset\mathbb{R}^{q}$ with $d,q\geq 1$. We allow either or both
of $X$ and $Y$ to be high-dimensional. The predictor $X$ or the response $Y$
can contain both continuous and categorical components. Our goal is to learn
the conditional distribution of $Y$ given $X=x$ and estimate the regression
function $\mathbb{E}(Y|X=x)$ in a unified framework.

We describe the proposed WGR method in detail below, which has three main
ingredients, a conditional distribution generator, a quadratic loss for
regression, and the Wasserstein metric for distribution matching.

###  2.1 Conditional generator

The theoretical foundation of WGR is the noise outsourcing lemma (Kallenberg,,
2002). This lemma states that, if $\mathcal{Y}$ is a standard Borel space
(Preston,, 2009), there exist a Borel-measurable function
$g^{*}:\mathcal{X}\times\mathbb{R}^{m}\to\mathcal{Y}$ and a random variable
$\eta\sim\text{Uniform}[0,1]$ such that $\eta$ is independent of $X$ and

$(X,Y)=(X,g^{*}(X,\eta))\ \ \ \ \text{almost surely.}$ (1)

We note that the condition that $\mathcal{Y}$ being a standard Borel space is
satisfied in all the applications we are interested in. In (1), for
simplicity, $\eta$ is taken to be a uniform random variable. In general, we
can take $\eta$ to be a random vector from a given reference distribution
$P_{\eta}$ that is easy to sample from. For example, we can take $P_{\eta}$ to
be the standard multivariate normal $N(\mathbf{0},\mbox{\bf I}_{m})$ with
$m\geq 1,$ which allows us to control the noise level more easily in practice.

We call the function $g^{*}$ in (1) a conditional generator, since if $g^{*}$
satisfies (1), it also satisfies

$g^{*}(\eta,x)\sim P_{Y|X=x},\eta\sim P_{\eta},x\in\mathcal{X}.$ (2)

So for a given $x$, to sample from the conditional distribution $P_{Y|X=x}$,
we can first generate $\eta\sim P_{\eta}$, then calculate $g^{*}(\eta,x)$,
which gives a sample from $P_{Y|X=x}.$ In addition, we can calculate any
moments of $P_{Y|X=x}$ via $g^{*}(\cdot,x)$. In particular, we have

$\mathbb{E}(Y|X=x)=\mathbb{E}_{\eta}[g^{*}(\eta,x)],x\in\mathcal{X}.$

In summary, we can determine the usual regression function (the conditional
mean) and sample from the conditional distribution as follows:

  * •

Regression function or conditional mean:
$\mathbb{E}(Y|X=x)=\mathbb{E}_{\eta}g^{*}(x;\eta),x\in\mathcal{X},$

  * •

Conditional distribution: $g^{*}(x;\eta)\sim
P_{Y|X=x},x\in\mathcal{X},\eta\sim P_{\eta}.$

Hence, the conditional generator provides a basis for a unified framework for
nonparametric regression and conditional distribution learning.

###  2.2 Objective function

Let $P_{X,g}$ denote the joint distribution of $(X,g(X,\eta))$, which is the
generated distribution based on a conditional generator $g(x,\eta),\eta\sim
P_{\eta},x\in\mathcal{X}.$ One possible way to measure the quality of a
conditional generator $g$ is to compare the generated distribution $P_{X,g}$
with the data distribution $P_{X,Y}$. A good conditional generator $g$ should
ensure that the generated distribution is close to the data distribution in
some sense. For example, one can use a distance metric such as the Wasserstein
distance or a divergence measure such as the Kullback-Leibler divergence to
quantify the discrepancy between the two distributions.

Let $\mathbb{D}$ be a divergence measure for the difference between $P_{X,g}$
and $P_{X,Y}$. Then, we formulate an objective function that combines this
divergence measure with the least squares loss to minimize the distribution
mismatch and the prediction error simultaneously. The objective function is

$\lambda_{w}{\mathbb{D}(P_{X,g}\|P_{X,Y})}+\lambda_{\ell}\mathbb{E}\|Y-\mathbb{E}_{\eta}g(X,\eta)\|^{2}.$
(3)

Here, both $\lambda_{\ell}$ and $\lambda_{w}$ are tuning parameters weighing
two losses, which are assumed to be nonnegative and
$\lambda_{\ell}+\lambda_{w}=1$. The objective function (3) combines two types
of losses: the first one evaluates how closely the generated distribution
$P_{X,g}$ resembles the data distribution $P_{X,Y};$ the second one is a
criterion quantifying how well the regression function fits the data.
Intuitively, the objective function (3) tries to learn the conditional
distribution of $Y$ given $X$ with the regularization that the conditional
mean is well estimated.

We take $\mathbb{D}$ to be the 1-Wasserstein distance. A computationally
convenient form of the 1-Wasserstein distance between $P_{X,g}$ and $P_{X,Y}$
is the Monge-Rubinstein dual (Villani,, 2009):

$\displaystyle\mathbb{D}_{W}(P_{X,g},P_{X,Y})=\sup_{f\in\mathcal{F}^{1}_{\text{Lip}}}\big{\\{}\mathbb{E}_{(X,\eta)}f(X,g(X,\eta))-\mathbb{E}_{(X,Y)}f(X,Y)\big{\\}},$
(4)

where
$\mathcal{F}^{1}_{\text{Lip}}=\\{f:\mathcal{X}\times\mathcal{Y}\to\mathbb{R},|f(u)-f(v)|\leq\|u-v\|_{2},\forall
u,v\in\mathcal{X}\times\mathcal{Y}\\}$ is a 1-Lipschitz class of functions on
$\mathcal{X}\times\mathcal{Y}.$ The Lipschitz function $f$ in (4) is often
called a critic or a discriminator.

Then, based on (3), the population objective function for the proposed
Wasserstein generative regression (WGR) is:

$\displaystyle
L(g)=\lambda_{w}\sup_{f\in\mathcal{F}^{1}_{\text{Lip}}}L_{\text{W}}(g,f)+\lambda_{\ell}L_{\text{LS}}(g),$
(5)

where

$\displaystyle L_{\text{W}}(g,f)$
$\displaystyle=\mathbb{E}_{(X,\eta)}f(X,g(X,\eta))-\mathbb{E}_{(X,Y)}f(X,Y),$
$\displaystyle L_{\text{LS}}(g)$
$\displaystyle=\mathbb{E}_{(X,Y)}\|Y-\mathbb{E}_{\eta}g(X,\eta)\|^{2}.$

Suppose we have a random sample $\\{(X_{i},Y_{i}),i=1,2,\ldots,n\\}$ from
$P_{X,Y}$, where $n\geq 1$ is the sample size. Let
$\\{\eta_{i},i=1,\ldots,n\\}$ and
$\\{\eta_{ij},i=1,2,\ldots,n,j=1,\ldots,J\\}$ with $J\geq 1$ be random
variables generated independently from $P_{\eta}$. We parameterize the
generator function $g$ and the discriminator $f$ by neural network functions
$g_{\bm{\theta}}$ and $f_{\bm{\phi}}$ with parameters (weights and biases)
$\bm{\theta}$ and $\bm{\phi}$, respectively. That is, we use neural network
functions to approximate the generator and critic functions and optimize the
objective function given below over the neural networks to obtain an estimator
of $g$. In addition, since $\mathbb{E}_{\eta}g_{\bm{\theta}}(X_{i},\eta)$
generally does not have a close form expression, we approximate it by the
sample average $J^{-1}\sum_{j=1}^{J}g_{\bm{\theta}}(X_{i},\eta_{ij}).$ Then,
the empirical objective function for estimating $(\bm{\theta},\bm{\phi})$ is

$\displaystyle\widehat{L}(g_{\bm{\theta}},f_{\bm{\phi}})=\lambda_{w}\widehat{L}_{\text{W}}(g_{\bm{\theta}},f_{\bm{\phi}})+\lambda_{\ell}\widehat{L}_{\text{LS}}(g_{\bm{\theta}}),$
(6)

where

$\displaystyle\widehat{L}_{\text{W}}(g_{\bm{\theta}},f_{\bm{\phi}})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left\\{f_{\bm{\phi}}(X_{i},g_{\bm{\theta}}(X_{i},\eta_{i}))-f_{\bm{\phi}}(X_{i},Y_{i})\right\\},$
$\displaystyle\widehat{L}_{\text{LS}}(g_{\bm{\theta}})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\big{\|}Y_{i}-\frac{1}{J}\sum_{j=1}^{J}g_{\bm{\theta}}(X_{i},\eta_{ij})\big{\|}^{2}.$

Let $(\hat{\bm{\theta}},\hat{\bm{\phi}})$ be a solution to the minimax problem

$\displaystyle(\hat{\bm{\theta}},\hat{\bm{\phi}})=\arg\min_{\bm{\theta}}\max_{\bm{\phi}}\widehat{L}(g_{\bm{\theta}},f_{\bm{\phi}}).$
(7)

Then, the estimated conditional generator is
$\hat{g}(x,\eta)=g_{\hat{\bm{\theta}}}(x,\eta)$ and the estimated regression
function is obtained by taking the expectation of $\hat{g}(x,\eta)$ with
respect to $\eta$, that is, $\hat{g}(x)=\mathbb{E}_{\eta}\hat{g}(x,\eta).$
Since there is no analytical expression for the expectation
$\mathbb{E}_{\eta}\hat{g}(x,\eta),$ we approximate it using an empirical
average based on a random sample
$\\{\eta_{1}^{\prime},\ldots,\eta_{K}^{\prime}\\},K\geq 1,$ from $P_{\eta}$,

$\hat{g}(x)\approx\frac{1}{K}\sum_{k=1}^{K}\hat{g}(x,\eta_{k}^{\prime}),$

which gives the estimated regression function.

We note that for a given $x\in\mathcal{X}$,
$\\{\hat{g}(x,\eta_{k}^{\prime}),k=1,\ldots,K\\}$ are approximately
distributed as $P_{Y|X=x}.$ We can use
$\\{\hat{g}(x,\eta_{k}^{\prime}),k=1,\ldots,K\\}$ to explore any aspects of
$P_{Y|X=x}$ that we are interested in such as its higher moments and
quantiles.

##  3 Implementation

In this section, we present the details for implementing WGR. We first
describe the neural networks used in the approximation of $g$ and $f$. We then
present the computational algorithm we implemented in detail.

###  3.1 ReLU Feedforward Neural Networks

We first give a brief description of feedforward neural networks (FNN) with
rectified linear unit (ReLU) activation function. The ReLU function is denoted
by $\sigma(x)\coloneqq\max(x,0)$, and it is defined for each component of $x$
if $x$ is a vector. A neural network can be expressed as a composite function
$\zeta(x)=\mathcal{L}_{H}\circ\sigma\circ\mathcal{L}_{H-1}\circ\sigma\circ\cdots\circ\sigma\circ\mathcal{L}_{1}\circ\sigma\circ\mathcal{L}_{0}(x),x\in\mathbb{R}^{p_{0}},$
where $\mathcal{L}_{i}(x)=W_{i}x+b_{i}$ with a weight matrix
$W_{i}\in\mathbb{R}^{p_{i+1}\times p_{i}}$ and bias vector
$b_{i}\in\mathbb{R}^{p_{i+1}}$ in the $i$-th linear transformation, and
$p_{i}$ is the width of the $i$-th layer, $i=0,1,\ldots,H$. The width and
depth of the network are described by $W=\max\\{p_{1},\ldots,p_{H}\\}$ and
$H$, respectively. To ease the presentation, we use
$\mathcal{NN}(p_{0},p_{H+1},W,H)$ to denote the neural networks with input
dimension $p_{0}$, output dimension $p_{H+1}$, width at most $W$ and depth at
most H.

We now specify the function classes below:

  * •

For the generator network class $\mathcal{G}$: Let
$\mathcal{G}\equiv\mathcal{NN}(d+m,q,W_{\mathcal{G}},H_{\mathcal{G}})$ be a
class of ReLU-activated FNNs,
$g_{\bm{\theta}}:\mathbb{R}^{d+m}\to\mathbb{R}^{q},$ with parameter
$\bm{\theta}$, width $W_{\mathcal{G}}$, and depth $H_{\mathcal{G}}$.

  * •

For the discriminator network class $\mathcal{D}$: Let
$\mathcal{D}\equiv\mathcal{NN}(d+q,1,W_{\mathcal{D}},H_{\mathcal{D}})\cap\text{Lip}(\Omega;K_{\mathcal{D}})$
be a class of ReLU-activated FNNs, $f_{\bm{\phi}}:\Omega\to\mathbb{R},$ with
parameter $\bm{\phi}$, width $W_{\mathcal{D}}$, and depth $H_{\mathcal{D}}$,
where for some $K_{\mathcal{D}}>0$, $\text{Lip}(\Omega;K_{\mathcal{D}})$ is a
class of Lipschitz functions defined below.

For any function $f:\Omega\to\mathbb{R}$, the Lipschitz constant of $f$ is
denoted by

$\text{Lip}(f)=\sup_{x,y\in\Omega,x\neq y}\frac{|f(x)-f(y)|}{\|x-y\|}.$

For a given $0<K<\infty$, denote $\text{Lip}(\Omega;K)$ as the set of all
functions $f:\Omega\mapsto\mathbb{R}$ with $\text{Lip}(f)\leq K$. And let
$\text{Lip}(\Omega;K,C):=\\{f\in\text{Lip}(\Omega;K):\|f\|_{\infty}\leq C\\}$,
where $\|f\|_{\infty}=\sup_{x\in\Omega}\|f(x)\|_{\infty}$ and $C$ is a
positive constant. Hence, $\mathcal{F}^{1}_{\text{Lip}}$ defined in (4) is
$\text{Lip}(\mathbb{R}^{d+q};1)$.

Algorithm 1 WGR Algorithm

1: Require: (a) Labeled data $\\{(X_{i},Y_{i})\\}_{i=1}^{n}$; (b) Minibatch
size $v\leq\min(n,N)$; (c) $J$, the size of noise vector $\eta$ to compute the
conditional expectation; $\lambda$, the gradient penalty parameter.

2: for number of training iterations do

3: Sample i.i.d. $\\{\eta_{ij},i=1,2,\ldots,n,j=1,\ldots,J\\}$ from the
standard multivariate normal distribution.

4: With fixed $\bm{\theta}^{(v)}$, update the discriminator $f_{\bm{\phi}}$ by
ascending its stochastic gradient:
$\displaystyle\nabla_{\phi}\frac{\lambda_{w}}{n}\sum_{i=1}^{n}\left\\{f_{\bm{\phi}}\left(X_{i},g_{\bm{\theta}}\left(X_{i},\eta_{i0}\right)\right)-f_{\bm{\phi}}\left(X_{i},Y_{i}\right)-\lambda\left(\left\|\nabla_{(x,y)}f_{\bm{\phi}}\left(X_{i},Y_{i}\right)\right\|_{2}-1\right)^{2}\right\\}$
(8) where the resulting parameter is denoted by $\bm{\phi}^{(v+1)}$.

5: With fixed $\bm{\phi}^{(v+1)}$, update the generator $g_{\bm{\theta}}$ by
descending its stochastic gradient:
$\displaystyle\nabla_{\theta}\left[\frac{\lambda_{l}}{n}\sum_{i=1}^{\tilde{n}}\left\\{Y_{i}-\frac{1}{J}\sum_{j=1}^{J}g_{\bm{\theta}}(X_{i},\eta_{ij})\right\\}^{2}+\frac{\lambda_{w}}{n}\sum_{i=1}^{n}f_{\bm{\phi}}\left(X_{i},g_{\bm{\theta}}\left(X_{i},\eta_{i0}\right)\right)\right],$
where the resulting parameter is denoted by $\bm{\theta}^{(v+1)}$.

6: end for

###  3.2 Computation

We now describe the implementation of WGR. For training the conditional
distribution generator $g_{\bm{\theta}}$ and the discriminator
$f_{\bm{\phi}}$, we use the leaky rectified linear unit (leaky ReLU) as the
activation function in $g_{\bm{\theta}}$ and $f_{\bm{\phi}}$. The training
algorithm is presented in Algorithm 1. We have implemented it in Pytorch.

To constrain the discriminator $f_{\bm{\phi}}$ to the class of 1-Lipschitz
functions, a gradient penalty is used in (8), which is a slightly modified
version of the algorithm proposed by (Gulrajani et al.,, 2017). The difference
is that we evaluate the gradients at the sample points in the penalty, instead
of using generated intermediate points. Another approach that we have tried to
enforce the Lipschitz constraint is the clipping method (Arjovsky et al.,,
2017), which also produces acceptable results, but seems to be less stable
than the penalty method described in Algorithm 1. We use traversal to select
the tuning parameters $\lambda_{l}$ and $\lambda_{w}$ that control the trade-
off between label and word embeddings. The constraints are that $\lambda_{l}$
and $\lambda_{w}$ must add up to 1 and have one decimal place each. In Section
5, we demonstrate the effectiveness of this algorithm in various numerical
experiments. However, we do not have a theoretical analysis of its convergence
behavior and we leave this as an open problem for future research.

##  4 Error analysis and convergence

In this section, we first develop an error decomposition, which decomposes the
estimation errors into approximation errors and stochastic errors of the
generator and discriminator. We then derive non-asymptotic error bounds for
WGR based on this error decomposition.

###  4.1 Error decomposition

We present a high-level description of the error decomposition for WGR. For
the estimator $\hat{g}$ define in (7), the estimation error consists of two
parts: the $L_{2}$-based excess risk
$\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}$,
and the integral probability metric (Müller,, 1997)
$d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})$ defined as

$d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})=\sup_{f\in\mathcal{F}^{1}_{B}}\\{\mathbb{E}_{(X,\eta)}f(X,\hat{g}(X,\eta))-\mathbb{E}_{X,Y}f(X,Y)\\},$

where
$\mathcal{F}^{1}_{B}=\\{f:\mathbb{R}^{d+q}\mapsto\mathbb{R},|f(z_{1})-f(z_{2})|\leq\|z_{1}-z_{2}\|,z_{1},z_{2}\in\mathbb{R}^{d+q},\|f\|_{\infty}\leq
B\\}$ is the bounded 1-Lipschitz function class. Clearly, if $P_{X,Y}$ has a
bounded support, $d_{\mathcal{F}^{1}_{B}}$ is the 1-Wasserstein distance. A
function class $\mathcal{F}$ is called symmetric if $f\in\mathcal{F}$ implies
$-f\in\mathcal{F}$.

We introduce a new error decomposition method in Lemma 1, which decomposes the
estimation error into approximation error and stochastic error of the
generator and discriminator.

######  Lemma 1.

Assume that the discriminator network class $\mathcal{D}$ is symmetric and the
probability measures of $(X,Y)$ and $(X,g(X,\eta))$ are supported on a compact
set $\Omega\subseteq\mathbb{R}^{d+q}$ for any $g\in\mathcal{G}$. Then, for the
WGR estimator defined in (7),

$\displaystyle\mathbb{E}_{\mathcal{S}}\\!\left\\{\\!\lambda_{\ell}\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)\\!-\\!\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}\\!+\\!\lambda_{w}d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\\!\right\\}$
$\displaystyle\quad\leq\lambda_{l}\mathcal{E}_{1}+4\lambda_{l}\mathcal{E}_{2}+2\mathcal{E}_{3}+2\lambda_{w}\mathcal{E}_{4}+3\lambda_{w}\mathcal{E}_{5}+3\lambda_{w}\mathcal{E}_{6},$
(9)

where $\mathcal{S}=\\{(X_{i},Y_{i},\eta_{i})\\}_{1\leq i\leq
n}\cup\\{\eta_{ij}\\}_{1\leq i\leq n,1\leq j\leq J}$ and

$\displaystyle\mathcal{E}_{1}$
$\displaystyle:=\mathbb{E}_{\mathcal{S}}\Big{\\{}\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\mathbb{E}\|Y-\mathbb{E}_{\eta}\hat{g}(X,\eta)\|^{2}-\frac{2}{n}\sum_{i=1}^{n}\|Y_{i}-\mathbb{E}_{\eta}\hat{g}(X_{i},\eta)\|^{2}\Big{\\}},$
$\displaystyle\mathcal{E}_{2}$
$\displaystyle:=\mathbb{E}_{\mathcal{S}}\Big{\\{}\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}|\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}g(X_{i},\eta_{ij})\|^{2}-\|Y_{i}-\mathbb{E}_{\eta}g(X_{i},\eta)\|^{2}|\Big{\\}},$
$\displaystyle\mathcal{E}_{3}$
$\displaystyle:=\inf_{g\in\mathcal{G}}\Big{[}\lambda_{l}\mathbb{E}\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}\sup_{f\in\mathcal{D}}\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\\}\Big{]},$
$\displaystyle\mathcal{E}_{4}$
$\displaystyle:=\sup_{h\in\mathcal{F}^{1}_{B}}\inf_{f\in\mathcal{D}}\|h-f\|_{\infty},$
$\displaystyle\mathcal{E}_{5}$
$\displaystyle:=\mathbb{E}_{\mathcal{S}}\Big{[}\sup_{f\in\mathcal{D}}\Big{\\{}\mathbb{E}f(X,Y)-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},Y_{i})\Big{\\}}\Big{]},$
$\displaystyle\mathcal{E}_{6}$
$\displaystyle:=\mathbb{E}_{\mathcal{S}}\Big{[}\sup_{f\in\mathcal{D},g\in\mathcal{G}}\Big{\\{}\mathbb{E}f(X,g(X,\eta))-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},g(X_{i},\eta_{i}))\Big{\\}}\Big{]}.$

According to their definitions,
$\mathcal{E}_{1},\mathcal{E}_{2},\mathcal{E}_{5}$ and $\mathcal{E}_{6}$ are
stochastic errors; $\mathcal{E}_{3}$ and $\mathcal{E}_{4}$ are approximation
errors. Lemma 1 provides a general error decomposition method, which covers
the error decomposition inequality in Jiao et al., (2023) for the traditional
nonparametric regression as a special case (corresponding to the case that
$\lambda_{\ell}=1,\lambda_{w}=0$ and $J=\infty$). It can also be utilized for
the error analysis for the conditional WGAN (corresponding to the case that
$\lambda_{\ell}=0,\lambda_{w}=1$). Moreover, when $\lambda_{\ell}=0$ and
$\lambda_{w}=1$, our error decomposition result in (9) is in line with that in
Lemma 9 in Huang et al., (2022) for the general GANs. The main difference is
that the upper bound of $\mathcal{E}_{6}$ in (9) depends on the sample size
$n$, while it is determined by the sample size of the generated noise $\eta$
in Huang et al., (2022). This is the key difference in the theoretical
analysis between conditional WGAN and general WGAN.

###  4.2 Non-asymptotic error bounds

More notations are needed. For $x\in\mathbb{R}^{d}$, its
$\ell_{1},\ell_{2},\ell_{\infty}$-norm is defined as
$\|x\|_{1}=\sum_{k=1}^{d}|x_{k}|$,
$\|x\|_{2}=(\sum_{k=1}^{d}|x_{k}|^{2})^{1/2}$ and $\|x\|_{\infty}=\max_{1\leq
k\leq d}|x_{k}|$, respectively. Let $\mathbb{N}$ be the set of positive
integers and $\mathbb{N}_{0}\coloneqq\mathbb{N}\cup\\{0\\}$. The maximum and
minimum of $A$ and $B$ are denoted by $A\vee B$ and $A\wedge B$. Let $\lfloor
A\rfloor$ be the largest integer strictly smaller than $A$ and $\lceil
A\rceil$ be the smallest integer strictly larger than $A$. For any $\beta>0$
and a set $\Gamma\subseteq\mathbb{R}^{m+d}$, the Hölder class of functions
$\mathcal{H}^{\beta}(\Gamma,B_{1})$ with a constant $0<B_{1}<\infty$ is
defined as

$\displaystyle\mathcal{H}^{\beta}(\Gamma,B_{1})=\left\\{f:\Gamma\rightarrow\mathbb{R},\max_{\|\alpha\|_{1}\leq\lfloor\beta\rfloor}\|D^{\alpha}f\|_{\infty}\leq
B_{1},\max_{\|\alpha\|_{1}=\lfloor\beta\rfloor}\sup_{x,y\in\Gamma,x\neq
y}\frac{|\partial^{\alpha}f(x)-\partial^{\alpha}f(y)|}{\|x-y\|^{r}}\leq
B_{1}\right\\},$

where $\partial^{\alpha}=\partial^{\alpha_{1}}\cdots\partial^{\alpha_{m+d}}$
with $\alpha=(\alpha_{1},\ldots,\alpha_{m+d})^{\top}\in\mathbb{N}_{0}^{m+d}$.

The following assumptions are needed.

######  Condition 1.

The probability measures of $(X,Y)$ and $(X,g(X,\eta))$ are supported on a
compact set
$\Omega_{X}\times\Omega_{Y}\subseteq[-B_{0},B_{0}]^{d+q}\subseteq\mathbb{R}^{d+q}$
for any $g\in\mathcal{G}$, where $0<B_{0}<\infty$ is a constant.

######  Condition 2.

The probability measure of $\eta$ is supported on
$\Omega_{\eta}\subseteq[-B_{0},B_{0}]^{m}$.

######  Condition 3.

For $g^{*}=(g^{*}_{1},\ldots,g^{*}_{q})^{\top}$,
$g^{*}_{k}\in\mathcal{H}^{\beta}(\Omega_{\eta}\times\Omega_{X},B_{1}),k=1,\ldots,q$,
where $\beta>0$ and $0<B_{1}<\infty$.

######  Condition 4.

For any $x\in\Omega_{X}$, there exists a vector $\eta_{x}\in\Omega_{\eta}$
such that for any $g,\tilde{g}\in\mathcal{G}$,

$\|\mathbb{E}_{\eta}g(x,\eta)-\mathbb{E}_{\eta}\tilde{g}(x,\eta)\|_{1}\leq\|g(x,\eta_{x})-\tilde{g}(x,\eta_{x})\|_{1}.$

Let $W,\bar{W},\bar{H}\in\mathbb{N}$, which may depend on $n$. We also make
the following assumptions on the network classes $\mathcal{D}$ and
$\mathcal{G}$.

######  ND 1.

The discriminator ReLU network class
$\mathcal{D}=\mathcal{NN}(d+q,1,W_{\mathcal{D}},H_{\mathcal{D}})\cap\text{Lip}([-B_{0},$
$B_{0}]^{d+q};K_{\mathcal{D}})$ has width
$W_{\mathcal{D}}=W^{d+q}\\{9(W+1)+5(d+q)-1\\}$, depth
$H_{\mathcal{D}}=3+14(d+q)(d+q-1)$ and Lipschitz constant $K_{\mathcal{D}}\leq
54B2^{d+q}(d+q)^{1/2}W^{2}$.

######  NG 1.

The generator ReLU network class
$\mathcal{G}=\mathcal{NN}(m+d,q,W_{\mathcal{G}},H_{\mathcal{G}})$ has width
$W_{\mathcal{G}}=38q(\lfloor\beta\rfloor+1)^{2}3^{(m+d)}(m+d)^{\lfloor\beta\rfloor+1}\bar{W}\lceil\log_{2}(8\bar{W})\rceil$
and depth $H_{\mathcal{G}}=21(\lfloor\beta\rfloor+1)^{2}\bar{H}$
$\lceil\log_{2}(8\bar{H})\rceil+2(m+d)$.

Conditions 1-3 require that $P_{X,Y}$, $P_{X,g}$ and $P_{\eta}$ have a bounded
support. Condition 3 is a smoothness condition for $g^{*}$ in (1). Condition 4
is a technical condition. The upper bound of the Lipschitz constant in ND 1 is
needed to achieve small approximation error. More details can be found in
Lemma 6 in the Appendix. To lighten the notations, we define the following two
quantities:

$\displaystyle a$
$\displaystyle:=\frac{\beta}{2\beta+\\{3(m+d)\\}\vee\\{2\beta(d+q+1)\\}},$
$\displaystyle b$
$\displaystyle:=\frac{3(m+d)}{2[2\beta+\\{3(m+d)\\}\vee\\{2\beta(d+q+1)\\}]}.$

######  Theorem 1.

Suppose that Conditions 1 \- 4 hold and the network parameters of
$\mathcal{D}$ and $\mathcal{G}$ satisfy ${\bf ND}$ 1 and ${\bf NG}$ 1 with
$W=\lceil n^{a}\rceil$, $\bar{W}=\lceil n^{b}/\log^{2}n\rceil$ and
$\bar{H}=\lceil\log n\rceil$. Then, for $J\gtrsim n$ and given weights
$\lambda_{l}$ and $\lambda_{w}$ satisfying $0<\lambda_{l},\lambda_{w}<1$ and
$\lambda_{l}+\lambda_{w}=1$, we have

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}_{X}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}\right\\}\leq
C_{1}n^{-\frac{\beta}{2\beta+\\{3(m+d)\\}\vee\\{2\beta(d+q+1)\\}}}(\log
n)^{\frac{2\beta}{m+d}\vee 1},$

where $C_{1}$ is a positive constant independent of $n$ and $J$.

Theorem 1 establishes a non-asymptotic upper bound for the excess risk of WGR
using deep neural networks. The convergence rates in Theorem 1 is slightly
slower than $O(n^{-\frac{2\beta}{2\beta+d}}\log^{c}n)$, the rate in deep least
squares nonparametric regression as in Jiao et al., (2023). This is because in
nonparametric regression, there is no distributional matching constraint, thus
the noise vector $\eta$ is not involved and a faster convergence rate can be
achieved. In our proposed framework, we are not only interested in estimating
the mean regression function, but also the conditional generator, which
involves the noise vector from a reference distribution. This increases the
dimensionality of the problem and results in a slower convergence rate.

We next establish a non-asymptotic error bound for the integral probability
metric $d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})$.

######  Theorem 2.

Suppose that those conditions of Theorem 1 hold. Then, for $J\gtrsim n$ and
given weights $\lambda_{l}$ and $\lambda_{w}$ satisfying
$0\leq\lambda_{l}<1,0<\lambda_{w}\leq 1$ and $\lambda_{l}+\lambda_{w}=1$, we
have

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\right\\}\leq
C_{2}n^{-\frac{\beta}{2\beta+\\{3(m+d)\\}\vee\\{2\beta(d+q+1)\\}}}(\log
n)^{\frac{2\beta}{m+d}\vee 1},$

where $C_{2}$ is a positive constant independent of $n$ and $J$.

The non-asymptotic error bounds in Theorem 1 and Theorem 2 are established for
fixed positive weights $\lambda_{l}$ and $\lambda_{w}$. In this case, we
obtain the same convergence rate of the excess risk
$\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}$
and the integral probability metric
$d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})$. Note that the conditional
GAN is involved in our proposed procedure, thus the joint stochastic error of
the discriminator and generator is affected by the dimension of the noise
vector $\eta,$ leading to a slower convergence rate compared with the one in
Theorem 5 in Huang et al., (2022) for GAN estimators.

Next, we present a non-asymptotic error bound for varying weights
$\lambda_{l}$ and $\lambda_{w}$, which can diverge with the sample size $n$.

######  Theorem 3.

Suppose that Conditions 1 \- 4 hold and the network parameters of
$\mathcal{D}$ and $\mathcal{G}$ satisfy ${\bf ND}$ 1 and ${\bf NG}$ 1 with
$W=\lceil n^{a}\rceil$, $\bar{W}=\lceil n^{b}/\log^{2}n\rceil$ and
$\bar{H}=\lceil\log n\rceil$. Then, for $\lambda_{l}>0,\lambda_{w}>0$
satisfying $\lambda_{l}+\lambda_{w}=1$ and $\lambda_{w}=O(n^{-1/(d+q+2)})$,
when $2\beta(d+q+1)\geq 3(m+d)+\beta$ and $J\gtrsim
n^{\\{3(m+d)+6\beta\\}/\\{4\beta(d+q+2)\\}}$, we have

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}\right\\}\leq
C_{3}n^{-\frac{3}{2(d+q+2)}}(\log n)^{\frac{2\beta}{m+d}\vee 2},$

where $C_{3}$ is a positive constant independent of $n$ and $J$. Moreover, as
$n\rightarrow\infty$,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\right\\}\rightarrow
0.$

When $\lambda_{w}=O(n^{-1/(d+q+2)})$, Theorem 3 gives an improved non-
asymptotic error bound of the proposed estimator, and it also implies that our
estimated distribution $P_{X,\hat{g}}$ converges weakly to $P_{X,Y}$ as
$n\to\infty$.

##  5 Numerical studies

In this section, a number of experiments including simulation studies and real
data examples are conducted to assess the performance of the proposed method.
We implement WGR in Pytorch, and use the stochastic gradient descent algorithm
RMSprop in the training process.

For comparison, we also compute the nonparametric least squares regression
using neural networks (NLS) and conditional Wasserstein GAN (cWGAN) (Arjovsky
et al.,, 2017; Liu et al.,, 2021). Aside from the results presented in this
section, additional numerical results are provided in the supplementary
materials, including the experiments with data generated from other models,
experiments investigating the effects of the noise dimension $m$ and the size
$J$, and neural networks with different architectures.

###  5.1 Simulation Studies

We conduct simulation studies to evaluate the performance of WGR with
univariate or multi-dimensional response $Y$. We use five different models to
generate data for our analysis. Each model has its own parameters and
assumptions, which are summarized in Table 1. The table also shows the sample
size and the architectures of the neural networks used in the analysis. We
compare the performance of WGR with other existing methods under these models.

######  Model 1.

A nonlinear regression model with an additive error term:

$Y=X_{1}^{2}+\exp\left(X_{2}+X_{3}/3\right)+\sin\left(X_{4}+X_{5}\right)+\varepsilon,$

where $\varepsilon\sim N(0,1).$

######  Model 2.

A nonlinear regression model with additive heteroscedastic error:

$Y=X_{1}^{2}+\exp(X_{2}+X_{3}/3)+X_{4}-X_{5}+(0.5+X_{2}^{2}/2+X_{5}^{2}/2)\varepsilon,$

where $\varepsilon\sim N(0,1)$.

######  Model 3.

A Gaussian mixture model:

$Y_{1}=X+\varepsilon_{1},\ Y_{2}=X+\varepsilon_{2},$

where $\varepsilon_{i}\sim
I_{\\{U_{i}<1/3\\}}N(-2,0.25^{2})+I_{\\{1/3<U_{i}<2/3\\}}N(0,0.25^{2})+I_{\\{U_{i}>2/3\\}}N(2,0.25^{2})$,
$i=1,2,$ and $U_{1}\sim\text{Uniform}(0,1)$, $U_{2}\sim\text{Uniform}(0,1)$.

######  Model 4.

Involute model:

$Y_{1}=2X+U\sin(2U)+\varepsilon_{1},\ Y_{2}=2X+U\cos(2U)+\varepsilon_{2},$

where $U\sim\text{Uniform}(0,2\pi)$, $\varepsilon_{1}\sim N(0,0.4^{2})$,
$\varepsilon_{2}\sim N(0,0.4^{2})$.

######  Model 5.

Octagon Gaussian mixture:

$Y_{1}=X+\varepsilon_{1},\ Y_{2}=X+\varepsilon_{2},$

where
$\varepsilon_{1}\\!\sim\\!I_{\\{U_{1}\in(i-1,i)\\}}N(\bm{\mu}_{i},\bm{\Sigma}_{i}),$
$\varepsilon_{2}\\!\sim\\!I_{\\{U_{2}\in(i-1,i)\\}}N(\bm{\mu}_{i},\bm{\Sigma}_{i}),$
$\bm{\mu}_{i}=(3\cos\frac{\pi_{i}}{4},3\sin\frac{\pi_{i}}{4})$,
$\bm{\Sigma}_{i}=\left(\\!\begin{array}[]{cc}\cos^{2}\frac{\pi
i}{4}+0.16^{2}\sin^{2}\frac{\pi
i}{4}\\!&\\!\left(1-0.16^{2}\right)\sin\frac{\pi i}{4}\cos\frac{\pi i}{4}\\\
\left(1-0.16^{2}\right)\sin\frac{\pi i}{4}\cos\frac{\pi
i}{4}\\!&\\!\sin^{2}\frac{\pi i}{4}+0.16^{2}\cos^{2}\frac{\pi
i}{4}\end{array}\\!\right)$, $U_{1}\sim\text{Uniform}(0,8)$,
$U_{2}\sim\text{Uniform}(0,8)$, $i=1,\ldots,8$.

Table 1: Description of simulation settings.

Data structure
---
Model | M1, M2 | M3, M4, M5
Response $(Y)$ | $Y\in\mathbb{R}$ | $Y\in\mathbb{Y}^{2}$
Covariate $(X)$ |  $N(0,I_{5})$; $\ N(0,I_{100})$ | $N(0,1)$
Noise $(\eta)$ | $N(0,I_{3})$ | $N(0,I_{10})$
Sample size
$J$ | 200 | 50
Training | 5000 | 40000
Validation | 1000 | 2000
Testing | 1000 | 10000
Network architecture
Generator network | (32, 16) | (512, 512, 512)
Discriminator network $f_{\bm{\phi}}$ | (32, 16) | (512, 512, 512)

Notes: In Model 1 and Model 2, $X\sim N(0,I_{5})$ is the intrinsic dimensional
case, and $X\sim N(0,I_{100})$ is the high dimensional case under sparsity
assumption. $J$ is the size of the random noise sample. The networks used in
the simulation are fully-connected feedforward neural networks with widths
specified above.

In the validation and testing stages, for each realization of $X$, we generate
$J=500$ i.i.d. noise samples $\\{\eta_{j},j=1,\ldots,J\\}$ from the standard
normal distribution and compute
$\\{\hat{g}_{\bm{\theta}}(X,\eta_{j}),j=1,\ldots,J\\}$. To evaluate the
predictive power of the WGR estimator $\hat{g}_{\theta}$, we use the $L_{1}$
and $L_{2}$ errors defined as

$\displaystyle
L_{1}=\frac{1}{n}\sum_{i=1}^{n}\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}\hat{g}_{\bm{\theta}}(X_{i},\eta_{ij})\|,\
L_{2}=\frac{1}{n}\sum_{i=1}^{n}\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}\hat{g}_{\bm{\theta}}(X_{i},\eta_{ij})\|^{2}.$
(10)

In addition, we can used the estimated conditional generator to obtain the
estimated $\tau$-th conditional quantile for a given $X=X_{i}$, denoted by
$\hat{F}_{Y|X}^{-1}(\tau|X=X_{i})$, via Monte Carlo. Also, we can calculate
the estimated conditional mean by
$\hat{E}(Y|X=X_{i})=(1/J)\sum_{j=1}^{J}g(X_{i},\eta_{ij})$ and the estimated
conditional standard deviation by
$\hat{SD}(Y|X=X_{i})=[(1/J)\sum_{j=1}^{J}\\{g(X_{i},\eta_{ij})-\hat{E}(Y|X=X_{i})\\}^{2}]^{1/2}$.
Hence, we can compute the mean squared error (MSE) of the estimated
conditional mean, standard deviation and the estimated $\tau$-th quantile,
defined by

$\displaystyle\text{MSE(mean)}=\frac{1}{K}\sum_{i=1}^{K}\\{\hat{E}(Y|X=X_{i})-E(Y|X=X_{i})\\}^{2},$
$\displaystyle\text{MSE(sd)}=\frac{1}{K}\sum_{i=1}^{K}\\{\hat{SD}(Y|X=X_{i})-SD(Y|X=X_{i})\\}^{2},$
$\displaystyle\text{MSE}(\tau)=\frac{1}{K}\sum_{i=1}^{K}\\{\hat{F}_{Y|X}^{-1}(\tau|X=X_{i})-F_{Y|X}^{-1}(\tau|X=X_{i})\\}^{2},$

where $K$ is the size of the validation or testing set. We consider
$\tau=0.05,0.25,0.50,$ $0.75,0.95$.

Truth cWGAN WGR

Model 3

Model 4

Model 5

Figure 1: Comparison of conditional density estimation. The abbreviations
represent the same as in Table 2. The conditional density functions are
estimated using 5000 samples, which are generated by the conditional samplers
from the methods given the randomly selected value of $X$.  Table 2:
Comparison of WGR with NLS and cWGAN for Models 1 and 2

Model | $d$ | Method | $L_{1}$ | $L_{2}$ | Mean | Sd
---|---|---|---|---|---|---
M1 | $5$ | NLS | 0.83(0.02) | 1.08(0.05) | - | -
cWGAN | 1.00(0.02) | 1.97(0.25) | 0.98(0.25) | 0.09(0.03)
WGR |  0.82(0.02) |  1.07(0.05) |  0.06(0.01) |  0.04(0.01)
$100$ | NLS | 1.17(0.03) | 2.52(0.15) | - | -
cWGAN | 1.17(0.04) | 2.65(0.41) | 1.67(0.39) | 0.81(0.02)
WGR |  1.15(0.04) |  2.40(0.24) |  1.64(0.25) |  0.16(0.04)
M2 | $5$ | NLS | 1.24(0.04) |  2.45(0.38) | - | -
cWGAN | 1.32(0.05) | 4.11(0.44) | 0.85(0.24) | 0.37(0.10)
WGR |  1.23(0.04) | 3.41(0.38) |  0.19(0.10) |  0.22(0.04)
$100$ | NLS | 1.63(0.05) | 5.37(0.47) | - | -
cWGAN | 1.64(0.06) | 5.34(0.46) | 2.20(0.22) | 1.11(0.08)
WGR |  1.61(0.06) |  5.27(0.45) |  2.06(0.24) |  0.30(0.03)

Notes: $d$ is the dimension of covariate $X$. The corresponding standard
errors are given in parentheses. The smallest $L_{1}$ and $L_{2}$ error, and
the smallest MSEs are in boldface. NLS represents the nonparametric least
squares regression, cWGAN represents the conditional WGAN, WGR is our proposed
Wasserstein generative regression method.

Table 3: MSE of the estimated conditional quantile at different quantile
levels

|  | $X\sim N(0,I_{5})$ | $X\sim N(0,I_{100})$
---|---|---|---
Model | $\tau$ | cWGAN | WGR | cWGAN | WGR
M1 | 0.05 | 1.22(0.23) |  0.29(0.07) | 3.18(0.71) |  1.84(0.23)
0.25 | 1.04(0.25) |  0.10(0.01) | 1.83(0.22) |  1.69(0.20)
0.50 | 0.99(0.26) |  0.09(0.02) | 1.75(0.16) |  1.66(0.16)
0.75 | 1.03(0.24) |  0.10(0.03) | 1.89(0.21) |  1.88(0.13)
0.95 | 1.34(0.21) |  0.23(0.06) | 3.61(0.39) |  2.41(0.14)
M2 | 0.05 | 1.86(0.21) |  0.77(0.09) | 4.99(0.51) |  3.42(1.07)
0.25 | 0.94(0.26) |  0.31(0.06) | 2.63(0.24) |  2.26(0.28)
0.50 | 0.85(0.26) |  0.19(0.04) | 2.21(0.22) |  2.19(0.22)
0.75 | 1.00(0.21) |  0.27(0.05) | 2.79(0.36) |  2.57(0.27)
0.95 | 2.59(0.52) |  0.81(0.15) | 5.41(0.65) |  3.49(0.47)

Notes: The notations are the same as in Table 2. The corresponding standard
errors are given in parentheses. The smallest MSEs are in boldface.

We repeat the simulations 10 times. For each evaluation criterion, the
simulation standard errors are computed and provided in parentheses. Table 2
summaries the average $L_{1}$ error, $L_{2}$ error, MSE(mean) and MSE(sd).
Table 3 reports the average MSE($\tau$) for different $\tau$. In Figure 1, we
visualize the quality of the conditional samples and the conditional density
estimation given a random realization of $X$.

It can be seen that, for Models 1 and 2, the three methods are comparable in
terms of $L_{1}$ and $L_{2}$ errors. But for the conditional mean, conditional
standard deviation, and conditional quantile estimation, WGR has smaller MSEs
values compared with cWGAN, indicating that WGR works better in distributional
matching.

For Models 3-5, Figure 1 shows the kernel-smoothing conditional density
estimates for a randomly selected value of $X$ based on 5,000 samples
generated using the estimated conditional generator. It can be seen that WGR
can better estimate the underlying conditional distributions for these models.

###  5.2 Real data examples

We demonstrate the effectiveness of WGR on four datasets: CT slides (Graf et
al.,, 2011), UJIndoorLoc (Torres-Sospedra et al.,, 2014), MNIST (LeCun et
al.,, 2010), and STL10 (Coates et al.,, 2011). The results from the STL10
dataset are given in the Appendix. Table 4 gives a summary of the dimensions
and training sizes of these datasets and the noise vectors used in the
analysis.

Table 4:  Summary of datasets and noise dimension and size

| CT slides | UJIndoorLoc | MNIST | STL10
---|---|---|---|---
Dimension of $X$ | 383 | 520 | 588 | 36477
Dimension of $Y$ | 1 | 6 | 144 | 12675
Training size | 40000 | 14948 | 20000 | 10000
Validation size | 3500 | 1100 | 1000 | 1000
Testing size | 10000 | 5000 | 10000 | 2000
Dimension of $\eta$ | 50 | 50 | 100 | 12675
Size of $\eta$ $(J)$ | 200 | 200 | 1 | 1

Notes: The noise vector $\eta$ is sampled from the multivariate standard
normal distribution. $J$ is the size of the noise vector generated in each
iteration.

####  5.2.1 The CT slices dataset

We evaluate the methods on the CT slices dataset (Graf et al.,, 2011) and
compare their prediction accuracy. This dataset can be found at the UCI
machine learning repository
(https://archive.ics.uci.edu/ml/datasets/Relative+location+of+CT+slices+on+axial+axis).
The dataset contains 53,500 CT images from 74 patients (43 male, 31 female)
with different anatomical landmarks annotated on the axial axis of the human
body. Each CT image is represented by two histograms in polar space: one for
the bone structures and another for the air inclusions inside the body. The
covariate vector consists of 383 variables: 239 for the bone histogram and 145
for the air histogram. The response variable is the relative location of the
image on the axial axis, which ranges from $0$ to $180$, where $0$ indicates
the top of the head and $180$ indicates the soles of the feet.

The sample size of this dataset is 53,500. We use 40,000 observations for
training, 3,500 observations for validation, and 10,000 observations for
testing. Both the generator network and the critic network have two hidden
layers with widths 128 and 64, respectively. The LeakyReLU activation function
is used in both networks. The noise vector $\eta$ is generated from
$N(\bm{0},\bm{I}_{50}).$ The number of the noise vectors $J$ is set to be
$200.$

For evaluation, besides the $L_{1}$ and $L_{2}$ errors in (10), we also
compute the average length of the estimated 95$\%$ prediction interval (PI)
and the corresponding coverage probability (CP), defined as

$\displaystyle\text{PI}=\\!\frac{1}{K}\sum_{i=1}^{K}\\{\hat{F}_{Y|X}^{-1}(0.975|X_{i})-\hat{F}_{Y|X}^{-1}(0.025|X_{i})\\},$
$\displaystyle\text{CP}=\\!\frac{1}{K}\sum_{i=1}^{K}I\\{Y_{i}\in\\![\hat{F}_{Y|X}^{-1}(0.025|X_{i}),\hat{F}_{Y|X}^{-1}(0.975|X_{i})]\\},$

where $\hat{F}_{Y|X}^{-1}(\tau|X_{i})$ is the estimated $\tau$-th conditional
quantile for a given $X=X_{i}$, and $K$ is the sample size of the validation
or testing set. The numerical results are summarized in Table 5.

Figure 2 shows the prediction intervals for 200 test samples, which are sorted
in ascending order according to the value of $Y$. We randomlyselect 200
samples from the test dataset and estimate the conditional prediction interval
based on 10,000 observations. The prediction intervals are sorted in acsending
order according to the value of $Y$ and are shown in Figure 2. In addition, we
display the estimated conditional density functions for 10 test samples in
Figure 2. The conditional density function is estimated using kernel smoothing
based on 10,000 values calculated from the conditional generator.

Table 5:  Summary statistics for CT test data Method | L1 | L2 | PI | CP
---|---|---|---|---
NLS | 0.40 | 0.51 | - | -
cWGAN | 0.95 | 2.30 | 1.54 | 0.48
WGR | 0.36 | 0.48 | 2.80 | 0.96

Prediction interval Estimated conditional density

cWGAN

WGR

Figure 2: The prediction intervals for the testing dataset and the estimated
conditional density functions of 10 randomly selected samples. In the left two
panels, the blue line represents the truth, the orange and the green lines
represent the 2.5$\%$ and the 97.5$\%$ quantiles, respectively. In the right
two panels, each color represents an observation from the test dataset.

As shown in Table 5, WGR and NLS have similar performance in terms of $L_{1}$
and $L_{2}$ errors, and both methods are superior to cWGAN. Furthermore, the
CP of WGR is close to the nominal level of 95$\%$ and much higher than cWGAN.
In addition, Figure 2 illustrates that the conditional distributions estimated
by cWGAN are more peaked than those of our proposed method, and this accounts
for why the prediction interval obtained by cWGAN covers fewer points than
WGR.

####  5.2.2 The UJIndoorLoc dataset

We present an analysis of the UJIndoor dataset (Torres-Sospedra et al.,,
2014), a multi-building multi-floor indoor localization database that relies
on WLAN/WiFi fingerprinting. The dataset can be downloaded from the UCI
machine learning repository
(https://archive.ics.uci.edu/ml/datasets/UJIIndoorLoc). This dataset contains
21,048 observations, which are divided into three parts: 14948 for training,
1100 for validation, and 5000 for testing. Each observation has 529
attributes. The attributes include the WiFi fingerprint, which is composed of
520 intensity values of different detected Wireless Access Points (WAPs),
ranging from -104dBm (very weak signal) to 0dBm (strong signal), and 100 for
non-detected WAPs. The attributes also include the location information, which
consists of six variables: longitude, latitude, floor, building ID, space ID,
and relative position. The first two variables are continuous, while the
others are categorical with at least two levels. We apply standardization to
the data before training. Our goal is to predict the location information from
the WiFi fingerprint using different machine learning methods and compare
their performance.

The neural networks used are two-layer fully connected feedforward networks
with 256 and 128 nodes, respectively. The LeakyReLU activation function is
used in both the conditional generator and critic networks. The noise vector
$\eta$ is generated from $N(\bm{0},\bm{I}_{50})$. And we use $J=200$.

Table 6 presents the analysis results and Figure 3 shows the prediction
intervals for building ID, space ID, and relative position in the response
vector, based on 200 samples randomly selected from the test dataset. Compared
with cWGAN, the prediction intervals of WGR have a higher coverage probability
with a comparable length.

Table 6:  Analysis results of the UJIndoorLoc testing dataset

Method |  | LNG | LAT | Floor | B-ID | S-ID | RP
---|---|---|---|---|---|---|---
NLS | $L_{1}$ | 0.07 | 0.09 | 0.12 | 0.05 | 0.16 | 0.34
$L_{2}$ | 0.06 | 0.10 | 0.05 | 0.12 | 0.21 | 0.58
PI | - | - | - | - | - | -
CP | - | - | - | - | - | -
cWGAN | $L_{1}$ | 0.14 | 0.17 | 0.23 | 0.12 | 0.29 | 0.23
$L_{2}$ | 0.04 | 0.07 | 0.11 | 0.04 | 0.36 | 0.11
PI | 0.24 | 0.26 | 0.31 | 0.22 | 0.61 | 1.06
CP | 0.55 | 0.50 | 0.44 | 0.53 | 0.68 | 0.49
WGR | $L_{1}$ | 0.08 | 0.10 | 0.18 | 0.05 | 0.20 | 0.41
$L_{2}$ | 0.01 | 0.03 | 0.06 | 0.01 | 0.16 | 0.51
PI | 0.24 | 0.27 | 0.37 | 0.24 | 0.74 | 2.29
CP | 0.75 | 0.71 | 0.60 | 0.89 | 0.84 | 0.58

Notes: WGR is the proposed method. LNG is Longitude, LAT is latitude, B-ID is
building ID, S-ID is space ID, RP is relative position.

building ID space ID relative position

cWGAN

WGR

Figure 3: Prediction intervals for building ID, space ID, and relative
position in the UJIndoorLoc test dataset. The blue line represents the truth,
the orange line represents the 2.5$\%$ quantile, and the green line represents
the 97.5$\%$ quantile.

We remark that the prediction intervals in these two data examples do not
account for the uncertainty in estimating the conditional generator. To
achieve theoretically valid coverage probability, one could use the conformal
prediction framework (Papadopoulos et al.,, 2002; Vovk et al.,, 2005) to
adjust the prediction intervals accordingly. However, this is problem is
beyond the scope of the current paper and we leave it for future work.

####  5.2.3 MNIST handwritten digits dataset

We now demonstrate the performance of WGR on a high-dimensional problem, where
both $X$ and $Y$ are high-dimensional. We use the MNIST dataset (LeCun et
al.,, 2010) of handwritten digits that can be downloaded from
http://yann.lecun.com/exdb/mnist/. The MNIST dataset consists of $28\times 28$
matrices of gray-scale images with values ranging from 0 to 1. Each image has
a corresponding label in $\\{0,1,\ldots,9\\}$. We apply WGR to the task of
reconstructing the central part of an image that is masked. We assume that the
masked part is the response $Y\in\mathbb{R}^{14\times 14}$ and the remaining
part is the covariate $X$, which has a dimension of $28\times 28-14\times
14=588$.

To evaluate the quality of the reconstructed images, we randomly sampled two
images per digit from the test set and compared the results of three different
methods in Figure 4. The figure shows that WGR produces sharper and more
faithful images than the other methods, as it preserves more details and
reduces artifacts.

$X$ Truth WGR NLS cWGAN

Figure 4: Reconstructed images in MNIST test data.

####  5.2.4 MNIST dataset: effects of sample sizes and network architectures.

We conduct experiments to investigate how the training sample sizes and the
network architectures affect the quality of the generated images using the
MNIST dataset. We apply WGR with training sample sizes $n=2,000$ and
$n=20,000$. The validation sample size is 1,000 and the test sample size is
10,000. For the network architectures, we use a network with 2 CNN layers and
a network with 3 fully-connected layers, respectively. Figure 5 displays the
reconstructed images with the two different network architectures. WGR with
CNN layers is more stable than NLS and cWGAN when the training sample size
varies. Moreover, it can be observed that WGR tends to generate images with
higher quality.

WGR NLS cWGAN

$X$ Truth $2000$ $20000$ $2000$ $20000$ $2000$ $20000$

CNN

Fully-Connected

Figure 5: Reconstructed images in MNIST test data.

##  6 Conclusions

In this paper, we have proposed a generative regression approach, Wasserstein
generative regression (WGR), for simultaneously estimating a regression
function and a conditional generator. We have provided theoretical support for
WGR by establishing its non-asymptotic error bounds and convergence
properties. Our numerical experiments demonstrate that it works well in
various situations from the standard generalized nonparametric regression
problems to more complex image reconstruction tasks.

WGR can be viewed as a way of estimating a conditional generator with a data-
dependent regularization on the first conditional moment. However, our
framework is not limited to this problem and can be adapted to other
estimation tasks by choosing different loss functions. For instance, we can
estimate the conditional median function or the conditional quantile function
by using other losses that are more suitable for these objectives. We can also
impose regularization on higher conditional moments or other properties of the
conditional distribution, depending on the research question.

Although we have established non-asymptotic error bounds and convergence
properties of WGR, our analysis is only a first attempt to deal with a
challenging technical problem that involves empirical processes on complex
functional spaces and approximation properties of deep neural networks.
Further work is needed to better understand the properties of generative
regression methods, including the proposed WGR. For instance, it would be
interesting to know if the error bounds we derived are optimal or if they can
be improved. WGR is a nonparametric method. For statistical inference and
model interpretation, it is desirable to incorporate a semiparametric
structure (Bickel et al.,, 1998) or a variable selection and dimension
reduction component in WGR (Chen et al.,, 2022; Huang et al.,, 2012).

Generative regression leverages the power of deep neural networks to model
complex and high-dimensional conditional distributions. Unlike traditional
regression methods that only output point estimates, generative regression can
capture the uncertainty and variability of the data by generating samples from
the learned distribution. This allows for more interpretable results in
various statistical applications. Therefore, we expect generative learning to
be a useful addition to the existing methods for prediction and inference in
statistics.

## References

  * Anthony and Bartlett, (1999) Anthony, M. and Bartlett, P. L. (1999).  Neural Network Learning: Theoretical Foundations.  Cambridge University Press.
  * Arjovsky et al., (2017) Arjovsky, M., Chintala, S., and Bottou, L. (2017).  Wasserstein generative adversarial networks.  In Precup, D. and Teh, Y. W., editors, Proceedings of the 34th International Conference on Machine Learning, pages 214–223.
  * Bartlett et al., (2019) Bartlett, P. L., Harvey, N., Liaw, C., and Mehrabian, A. (2019).  Nearly-tight vc-dimension and pseudodimension bounds for piecewise linear neural networks.  The Journal of Machine Learning Research, 20(1):2285–2301.
  * Bickel et al., (1998) Bickel, P. J., Klaassen, C. A. J., Yaácov, R., and Wellner, J. A. (1998).  Efficient and Adaptive Estimation for Semiparametric Models.  Springer, New York.
  * Chen et al., (2022) Chen, Y., Gao, Q., and Wang, X. (2022).  Inferential Wasserstein generative adversarial networks.  Journal of the Royal Statistical Society Series B, 84(1):83–113.
  * Coates et al., (2011) Coates, A., Ng, A., and Lee, H. (2011).  An analysis of single-layer networks in unsupervised feature learning.  In Proceedings of the 14th International Conference on Artificial Intelligence and Statistics, pages 215–223.
  * Fan and Gijbels, (1996) Fan, J. and Gijbels, I. (1996).  Local polynomial modelling and its applications.  Monographs on statistics and applied probability series. Chapman & Hall, London.
  * Friedman and Stuetzle, (1981) Friedman, J. H. and Stuetzle, W. (1981).  Projection pursuit regression.  J. Amer. Statist. Assoc., 76(376):817–823.
  * Goodfellow et al., (2014) Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. (2014).  Generative adversarial nets.  In Proceedings of the 27th International Conference on Neural Information Processing Systems, pages 2672–2680.
  * Graf et al., (2011) Graf, F., Kriegel, H.-P., ólsterl, S., and Schubert, M. (2011).  Position prediction in ct volume scans.  In Proceedings of the 28th International Conference on Machine Learning. Workshop on Learning for Global Challenges.
  * Gulrajani et al., (2017) Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V., and Courville, A. C. (2017).  Improved training of Wasserstein GANs.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, page 5769–5779.
  * Györfi et al., (2002) Györfi, L., Kohler, M., Krzyżak, A., and Walk, H. (2002).  A Distribution-Free Theory of Nonparametric Regression.  Springer-Verlag, New York.
  * Hardle et al., (1993) Hardle, W., Hall, P., and Ichimura, H. (1993).  Optimal smoothing in single-index models.  Annals of Statistics, 21(1):157–178.
  * Hastie and Tibshirani, (1986) Hastie, T. and Tibshirani, R. (1986).  Generalized additive models.  Statistical Science, 1(3):297 – 310.
  * Heusel et al., (2017) Heusel, M., Ramsauer, H., Unterthiner, T., Nessler, B., and Hochreiter, S. (2017).  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.
  * Hon and Yang, (2022) Hon, S. and Yang, H. (2022).  Simultaneous neural network approximation for smooth functions.  Neural Networks, 154:152–164.
  * Huang et al., (2012) Huang, J., Breheny, P., and Ma, S. (2012).  A selective review of group selection in high-dimensional models.  Statistical Science, 27(4):481 – 499.
  * Huang et al., (2022) Huang, J., Jiao, Y., Li, Z., Liu, S., Wang, Y., and Yang, Y. (2022).  An error analysis of generative adversarial networks for learning distributions.  Journal of Machine Learning Research, 23(116):1–43.
  * Ichimura, (1993) Ichimura, H. (1993).  Semiparametric least squares (sls) and weighted sls estimation of single-index models.  Journal of Econometrics, 58(1):71–120.
  * Jiao et al., (2023) Jiao, Y., Shen, G., Lin, Y., and Huang, J. (2023).  Deep nonparametric regression on approximate manifolds: Nonasymptotic error bounds and polynomial prefactors.  Annals of Statistics, 51(2):691–716.
  * Kallenberg, (2002) Kallenberg, O. (2002).  Foundations of Modern Probability.  Springer, New York.
  * Kovachki et al., (2021) Kovachki, N., Baptista, R., Hosseini, B., and Marzouk, Y. (2021).  Conditional sampling with monotone gans.  arXiv 2006.06755.
  * LeCun et al., (2010) LeCun, Y., Cortes, C., and Burges, C. (2010).  Mnist handwritten digit database.  AT&T Labs [Online]. Available: http://yann. lecun. com/exdb/mnist, 2.
  * Liu et al., (2021) Liu, S., Zhou, X., Jiao, Y., and Huang, J. (2021).  Wasserstein generative learning of conditional distribution.  arXiv preprint arXiv:2112.10039.
  * Mirza and Osindero, (2014) Mirza, M. and Osindero, S. (2014).  Conditional generative adversarial nets.  arxiv:1411.1784.
  * Müller, (1997) Müller, A. (1997).  Integral probability metrics and their generating classes of functions.  Advances in Applied Probability, pages 429–443.
  * Nakada and Imaizumi, (2020) Nakada, R. and Imaizumi, M. (2020).  Adaptive approximation and generalization of deep neural network with intrinsic dimensionality.  Journal of Machine Learning Research, 21(174):1–38.
  * Papadopoulos et al., (2002) Papadopoulos, H., Proedrou, K., Vovk, V., and Gammerman, A. (2002).  Inductive confidence machines for regression.  In Proceedings of 13th European Conference on Machine Learning, pages 345–356. Springer.
  * Preston, (2009) Preston, C. (2009).  A note on standard Borel and related spaces.  Journal of Contemporary Mathematical Analysis, 44:63–71.
  * Reed et al., (2016) Reed, S., Akata, Z., Yan, X., Logeswaran, L., Schiele, B., and Lee, H. (2016).  Generative adversarial text to image synthesis.  In ICML.
  * Salakhutdinov, (2015) Salakhutdinov, R. (2015).  Learning deep generative models.  Annual Review of Statistics and Its Application, 2:361–385.
  * Scott, (1992) Scott, D. W. (1992).  Multivariate Density Estimation: Theory, Practice and Visualization.  Wiley, New York.
  * Silverman, (1986) Silverman, B. W. (1986).  Density Estimation for Statistics and Data Analysis.  Chapman & Hall, London.
  * Stone, (1986) Stone, C. J. (1986).  The dimensionality reduction principle for generalized additive models.  Ann. Statist., 14(2):590–606.
  * Torres-Sospedra et al., (2014) Torres-Sospedra, J., Montoliu, R., Martínez-Usó, A., Avariento, J. P., Arnau, T. J., Benedito-Bordonau, M., and Huerta, J. (2014).  Ujiindoorloc: A new multi-building and multi-floor database for wlan fingerprint-based indoor localization problems.  In International Conference on Indoor Positioning and Indoor Navigation, pages 261–270.
  * Tsybakov, (2008) Tsybakov, A. (2008).  Introduction to Nonparametric Estimation.  Springer Science & Business Media.
  * Villani, (2009) Villani, C. (2009).  Optimal Transport: Old and New.  Springer.
  * Vovk et al., (2005) Vovk, V., Gammerman, A., and Shafer, G. (2005).  Algorithmic Learning in a Random World.
  * Wasserman, (2006) Wasserman, L. (2006).  All of Nonparametric Statistics (Springer Texts in Statistics).  Springer-Verlag, Berlin, Heidelberg.
  * Zhou et al., (2022) Zhou, X., Jiao, Y., Liu, J., and Huang, J. (2022).  A deep generative approach to conditional sampling.  Journal of the American Statistical Association, 0:1–12.
  * Zhu et al., (2017) Zhu, J.-Y., Park, T., Isola, P., and Efros, A. A. (2017).  Unpaired image-to-image translation using cycle-consistent adversarial networks.  In ICCV.

Appendices

In this appendix, we provide detailed proofs of the main theorems and
additional numerical experiments, including simulation studies and real data
examples.

##  Appendix A Proof of Lemma 1

We first recall the error decomposition. For the proposed estimator $\hat{g}$,
the estimation errors can be decomposed as:

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\lambda_{\ell}\mathbb{E}\left\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\right\|^{2}+\lambda_{w}d_{\mathcal{F}_{B}^{1}}\left(P_{X,\hat{g}},P_{X,Y}\right)\right\\}$
$\displaystyle\leq\lambda_{l}\mathcal{E}_{1}+4\lambda_{l}\mathcal{E}_{2}+2\mathcal{E}_{3}+2\lambda_{w}\mathcal{E}_{4}+3\lambda_{w}\mathcal{E}_{5}+3\lambda_{w}\mathcal{E}_{6},$

where $\mathcal{S}=\left\\{\left(X_{i},Y_{i},\eta_{i}\right)\right\\}_{1\leq
i\leq n}\cup\left\\{\eta_{ij}\right\\}_{1\leq i\leq n,1\leq j\leq J}$ and

$\displaystyle\mathcal{E}_{1}:=\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}\left\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\right\|^{2}+\mathbb{E}\left\|Y-\mathbb{E}_{\eta}\hat{g}(X,\eta)\right\|^{2}-\frac{2}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}\left(X_{i},\eta\right)\right\|^{2}\right\\},$
(A.1)
$\displaystyle\mathcal{E}_{2}:=\mathbb{E}_{\mathcal{S}}\left\\{\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\left|\left\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}g\left(X_{i},\eta_{ij}\right)\right\|^{2}-\left\|Y_{i}-\mathbb{E}_{\eta}g\left(X_{i},\eta\right)\right\|^{2}\right|\right\\},$
(A.2)
$\displaystyle\mathcal{E}_{3}:=\inf_{g\in\mathcal{G}}\left[\lambda_{l}\mathbb{E}\left\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\right\|^{2}+\lambda_{w}\sup_{f\in\mathcal{D}}\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\\}\right],$
(A.3)
$\displaystyle\mathcal{E}_{4}:=\sup_{h\in\mathcal{F}_{B}^{1}}\inf_{f\in\mathcal{D}}\|h-f\|_{\infty},$
(A.4)
$\displaystyle\mathcal{E}_{5}:=\mathbb{E}_{\mathcal{S}}\left[\sup_{f\in\mathcal{D}}\left\\{\mathbb{E}f(X,Y)-\frac{1}{n}\sum_{i=1}^{n}f\left(X_{i},Y_{i}\right)\right\\}\right],$
(A.5)
$\displaystyle\mathcal{E}_{6}:=\mathbb{E}_{\mathcal{S}}\left[\sup_{f\in\mathcal{D},g\in\mathcal{G}}\left\\{\mathbb{E}f(X,g(X,\eta))-\frac{1}{n}\sum_{i=1}^{n}f\left(X_{i},g\left(X_{i},\eta_{i}\right)\right)\right\\}\right].$
(A.6)

We will prove Lemma 1 in several steps.

Step 1. To rewrite the estimation error. For the generalized nonparametric
regression model $Y=g^{*}(X;\eta),$ since the proposed objective function is

$\displaystyle
L(g,f)=\lambda_{\ell}\mathbb{E}\left\|Y-\mathbb{E}_{\eta}g(X,\eta)\right\|^{2}+\lambda_{w}\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\\},$

we have
$\sup_{f\in\mathcal{F}^{1}_{B}}L(g^{*},f)=\lambda_{\ell}\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}.$
Then, we write the estimation error

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\lambda_{\ell}\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\right\\}$
$\displaystyle\quad=\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{F}^{1}_{B}}L(\hat{g},f)-\sup_{f\in\mathcal{F}^{1}_{B}}L(g^{*},f)\right\\}$
$\displaystyle\quad=\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{F}^{1}_{B}}L(\hat{g},f)-\sup_{f\in\mathcal{D}}L(\hat{g},f)\right\\}+\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{D}}L(\hat{g},f)-\sup_{f\in\mathcal{F}^{1}_{B}}L(g^{*},f)\right\\}.$
(A.7)

Step 2. To bound the first term in (A.7). By Lemma 24 in Huang et al., (2022),

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{F}^{1}_{B}}L(\hat{g},f)-\sup_{f\in\mathcal{D}}L(\hat{g},f)\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})-d_{\mathcal{D}}(P_{X,\hat{g}},P_{X,Y})\right\\}$
$\displaystyle\leq$ $\displaystyle
2\sup_{h\in\mathcal{F}^{1}_{B}}\inf_{f\in\mathcal{D}}\|h-f\|_{\infty}=2\mathcal{E}_{4},$

where $\mathcal{E}_{4}$ is defined in (A.4).

Step 3. To bound the second term in (A.7). By the definition of $L(g,f)$ and
$g^{*}$,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{D}}L(\hat{g},f)-\sup_{f\in\mathcal{F}^{1}_{B}}L(g^{*},f)\right\\}$
$\displaystyle=$
$\displaystyle\lambda_{\ell}\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}\|Y-\mathbb{E}_{\eta}\hat{g}(X,\eta)\|^{2}-\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}\right\\}+\lambda_{w}\mathbb{E}_{\mathcal{S}}\left[d_{\mathcal{D}}(P_{X,\hat{g}},P_{X,Y})\right].$
(A.8)

We decompose the second term in (A) as

$\displaystyle\mathbb{E}_{\mathcal{S}}[d_{\mathcal{D}}(P_{X,\hat{g}},P_{X,Y})]$
$\displaystyle\leq\mathbb{E}_{\mathcal{S}}\left[\sup_{f\in\mathcal{D}}\left\\{\mathbb{E}f(X,Y)-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},Y_{i})\right\\}\right]$
$\displaystyle\quad+\mathbb{E}_{\mathcal{S}}\left[\sup_{f\in\mathcal{D},g\in\mathcal{G}}\left\\{\mathbb{E}f(X,g(X,\eta))-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},g(X_{i},\eta_{i}))\right\\}\right]$
$\displaystyle\quad+\mathbb{E}_{\mathcal{S}}\left[\sup_{f\in\mathcal{D}}\left\\{\frac{1}{n}\sum_{i=1}^{n}f(X_{i},\hat{g}(X_{i},\eta_{i}))-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},Y_{i})\right\\}\right]$
$\displaystyle\coloneqq\mathcal{E}_{5}+\mathcal{E}_{6}+\Delta,$ (A.9)

where $\mathcal{E}_{5}$ and $\mathcal{E}_{6}$ are defined in (A.5) and (A.6).
Note that $\Delta$ in (A) is non-negative as $\mathcal{D}$ is symmetric. To
deal with $\Delta$, we introduce a new estimator $\bar{g}$ defined as

$\bar{g}=\arg\inf_{g\in\mathcal{G}}\left[\lambda_{l}\mathbb{E}\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}\sup_{f\in\mathcal{D}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\right\\}\right].$
(A.10)

It then follows from the defintion of $\hat{g}$ that

$\displaystyle\Delta\leq\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{D}}\hat{L}(\hat{g},f)\right\\}\leq\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{D}}\hat{L}(\bar{g},f)\right\\}.$
(A.11)

Then, similar to (A), we have

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{D}}\hat{L}(\bar{g},f)\right\\}\leq$
$\displaystyle\lambda_{w}\mathcal{E}_{5}+\lambda_{w}\mathcal{E}_{6}+\lambda_{w}\sup_{f\in\mathcal{D}}\left\\{\mathbb{E}f(X,\bar{g}(X,\eta))-\mathbb{E}f(X,Y)\right\\}$
$\displaystyle\quad+\lambda_{l}\mathbb{E}_{\mathcal{S}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}\bar{g}(X_{i},\eta_{ij})\right\|^{2}\right\\}.$
(A.12)

To bound the last term in (A), we subtract
$\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}$ on both sides of (A).
Then,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}\bar{g}(X_{i},\eta_{ij})\right\|^{2}\right\\}-\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}\bar{g}(X_{i},\eta_{ij})\right\|^{2}\right\\}-\mathbb{E}_{\mathcal{S}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\|Y_{i}-\mathbb{E}_{\eta}\bar{g}(X_{i},\eta)\|^{2}\right\\}$
$\displaystyle+\mathbb{E}\|\mathbb{E}_{\eta}\bar{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}$
$\displaystyle\leq$
$\displaystyle\mathcal{E}_{2}+\mathbb{E}\|\mathbb{E}_{\eta}\bar{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2},$

where $\mathcal{E}_{2}$ is defined in (A.2). Then, by the definition of
$\bar{g}$,

$\displaystyle\mathbb{E}\|\mathbb{E}_{\eta}\bar{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}$
$\displaystyle\leq$
$\displaystyle\inf_{g\in\mathcal{G}}\left[\lambda_{l}\mathbb{E}\left\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\right\|^{2}+\lambda_{w}\sup_{f\in\mathcal{D}}\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\\}\right]\coloneqq\mathcal{E}_{3},$

where $\mathcal{E}_{3}$ is defined in (A.3). Thus,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{D}}\hat{L}(\hat{g},f)\right\\}-\lambda_{l}\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}\leq\lambda_{w}\mathcal{E}_{5}+\lambda_{w}\mathcal{E}_{6}+\lambda_{\ell}\mathcal{E}_{2}+\mathcal{E}_{3}.$
(A.13)

Further, by the defintion of $\hat{L}(g,f)$,

$\displaystyle-\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{D}}\hat{L}(\hat{g},f)\right\\}+\lambda_{l}\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}$
$\displaystyle=$
$\displaystyle-\lambda_{l}\mathbb{E}_{\mathcal{S}}\left[\sup_{f\in\mathcal{D}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}\hat{g}(X_{i},\eta_{ij})\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}\left(X_{i},\eta\right)\right\|^{2}\right\\}\right.$
$\displaystyle-\left.\lambda_{w}\\{\frac{1}{n}\sum_{i=1}^{n}f\left(X_{i},\hat{g}\left(X_{i},\eta_{i}\right)\right)-\frac{1}{n}\sum_{i=1}^{n}f\left(X_{i},Y_{i}\right)\right\\}$
$\displaystyle+\left.\\{\lambda_{l}\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}\left(X_{i},\eta\right)\right\|^{2}\\}\right]$
$\displaystyle\leq$
$\displaystyle\lambda_{l}\mathcal{E}_{2}-\lambda_{w}\Delta+\lambda_{l}\mathbb{E}_{\mathcal{S}}\\{\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}\left(X_{i},\eta\right)\right\|^{2}\\}.$
(A.14)

Combining (A.13) and (A), we have

$\displaystyle\ \lambda_{w}\Delta\leq$
$\displaystyle\lambda_{w}\mathcal{E}_{5}+\lambda_{w}\mathcal{E}_{6}+2\lambda_{\ell}\mathcal{E}_{2}+\mathcal{E}_{3}+\lambda_{l}\mathbb{E}_{\mathcal{S}}\\{\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}\left(X_{i},\eta\right)\right\|^{2}\\}.$

As a result,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\sup_{f\in\mathcal{D}}L(\hat{g},f)-\sup_{f\in\mathcal{F}_{B}^{1}}L\left(g^{*},f\right)\right\\}$
$\displaystyle\leq$
$\displaystyle\lambda_{\ell}\mathbb{E}_{S}\left\\{\mathbb{E}\left\|Y-\mathbb{E}_{\eta}\hat{g}(X,\eta)\right\|^{2}-\mathbb{E}\left\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\right\|^{2}\right\\}+2\lambda_{w}\mathcal{E}_{5}+2\lambda_{w}\mathcal{E}_{6}+2\lambda_{l}\mathcal{E}_{2}+\mathcal{E}_{3}$
$\displaystyle+\lambda_{l}\mathbb{E}_{\mathcal{S}}\\{\mathbb{E}\|Y-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}\left(X_{i},\eta\right)\right\|^{2}\\}$
$\displaystyle\leq$ $\displaystyle
2\lambda_{w}\mathcal{E}_{5}+2\lambda_{w}\mathcal{E}_{6}+2\lambda_{l}\mathcal{E}_{2}+\mathcal{E}_{3}+\mathbb{E}_{S}\left\\{\mathbb{E}\left\|Y-\mathbb{E}_{\eta}\hat{g}(X,\eta)\right\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}\left(X_{i},\eta\right)\right\|^{2}\right\\}$
$\displaystyle\coloneqq$ $\displaystyle
2\lambda_{w}\mathcal{E}_{5}+2\lambda_{w}\mathcal{E}_{6}+2\lambda_{l}\mathcal{E}_{2}+\mathcal{E}_{3}+\lambda_{l}\mathcal{E}_{1},$

where $\mathcal{E}_{1}$ is given in (A.1).

Consequently,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\lambda_{\ell}\mathbb{E}\left\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\right\|^{2}+\lambda_{w}d_{\mathcal{F}_{B}^{1}}\left(P_{X,\hat{g}},P_{X,Y}\right)\right\\}$
$\displaystyle\leq\lambda_{l}\mathcal{E}_{1}+4\lambda_{l}\mathcal{E}_{2}+2\mathcal{E}_{3}+2\lambda_{w}\mathcal{E}_{4}+3\lambda_{w}\mathcal{E}_{5}+3\lambda_{w}\mathcal{E}_{6}.$

The proof of Lemma 1 is complete.

##  Appendix B Proofs of the main theorems

### Proof of Theorem 1.

For any fixed $\lambda_{l}$ and $\lambda_{w}$ satisfying
$0<\lambda_{l},\lambda_{w}<1$, it follows from Lemma 1 that

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}\right\\}$
$\displaystyle\leq$
$\displaystyle\frac{1}{\lambda_{l}}\mathbb{E}_{\mathcal{S}}\\{\lambda_{l}\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\\}$
$\displaystyle\leq$
$\displaystyle\frac{1}{\lambda_{l}}(\lambda_{l}\mathcal{E}_{1}+4\lambda_{l}\mathcal{E}_{2}+2\mathcal{E}_{3}+2\lambda_{w}\mathcal{E}_{4}+3\lambda_{w}\mathcal{E}_{5}+3\lambda_{w}\mathcal{E}_{6}).$
(A.15)

The upper bounds of the error terms in (B) are given in Lemmas 3 to 9 in
Section C. Thus, under Conditions 1-4, we have

$\displaystyle\mathbb{E}_{\mathcal{S}}\\{\lambda_{l}\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\\}$
$\displaystyle\leq$ $\displaystyle\lambda_{l}C_{1}(B_{0}\vee
1)^{2}q^{2}\frac{H_{\mathcal{G}}S_{\mathcal{G}}\log S_{\mathcal{G}}\log n}{n}$
$\displaystyle+16\lambda_{l}C_{2}qB_{0}^{2}\sqrt{\frac{H_{\mathcal{G}}S_{\mathcal{G}}\log
S_{\mathcal{G}}\log J}{J}}$
$\displaystyle+2({K_{\mathcal{D}}\lambda_{w}}+2B_{0}q^{2}\lambda_{l})19B_{1}(\lfloor\beta\rfloor+1)^{2}(m+d)^{\lfloor\beta\rfloor+1/2+(\beta\vee
1)/2}(\bar{W}\bar{H})^{-2\beta/(m+d)}$
$\displaystyle+2\lambda_{w}2^{d+q}(d+q)BW^{-1}$
$\displaystyle+3\lambda_{w}C_{3}\\{B+2^{d+q}(d+q)B\\}\sqrt{\frac{H_{\mathcal{D}}S_{\mathcal{D}}\log(S_{\mathcal{D}})\log
n}{n}}$
$\displaystyle+3\lambda_{w}C_{4}B\\{B+2^{d+q}(d+q)B\\}\sqrt{\frac{(H_{\mathcal{G}}+H_{\mathcal{D}}+1)S_{\mathcal{G},\mathcal{D}}\log(S_{\mathcal{G},\mathcal{D}})\log
n}{n}},$ (A.16)

where $C_{1},C_{2},C_{3},C_{4}$ are positive constants independent of
$n,J,\beta,B_{0},B_{1},B$, $S_{\mathcal{G},\mathcal{D}}$ is the size of the
network in
$\mathcal{NN}(m+d,1,2W_{\mathcal{G}}+2W_{\mathcal{D}}+2d,H_{\mathcal{G}}+H_{\mathcal{D}}+1)$,
and the Lipschitz constant for the discriminator network class
$K_{\mathcal{D}}\leq 54B2^{d+q}(d+q)^{1/2}W^{2}$.

To simplify (A.16), we let the network parameters $W$, $\bar{W}$ and $\hat{H}$
in the discriminator network class ND 1 and generator network class NG 1 be

$\displaystyle W=\lceil n^{a}\rceil,\bar{W}=\lceil
n^{b}/\log^{2}n\rceil,\bar{H}=\lceil\log n\rceil,$

where

$\displaystyle a$
$\displaystyle=\frac{\beta}{2\beta+\\{3(m+d)\\}\vee\\{2\beta(d+q+1)\\}},$
$\displaystyle b$
$\displaystyle=\frac{3(m+d)}{2[2\beta+\\{3(m+d)\\}\vee\\{2\beta(d+q+1)\\}]}.$

Finally,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\lambda_{l}\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\right\\}$
$\displaystyle\quad\leq
Cn^{-\frac{\beta}{2\beta+\\{3(m+d)\\}\vee\\{2\beta(d+q+1)\\}}}(\log
n)^{\frac{2\beta}{m+d}\vee 1},$

where $C$ is independent of $(n,J)$. The proof of Theorem 1 is complete.

### Proof of Theorem 2.

For any fixed $\lambda_{l}$ and $\lambda_{w}$ satisfying
$0<\lambda_{l},\lambda_{w}<1$,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\right\\}\leq\frac{1}{\lambda_{w}}\mathbb{E}_{\mathcal{S}}\\{\lambda_{l}\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\\}.$

Thus, Theorem 2 can be proved using the results in the proof of Theorem 1. We
omit the details here.

### Proof of Theorem 3:

Following similar arguments as in the proof of Theorem 1, we know that the
Lipschitz constant $K_{\mathcal{D}}$ in (A.16) satisfies $K_{\mathcal{D}}\leq
54B2^{d+q}(d+q)^{1/2}\lceil n^{2a}\rceil.$ Thus, for any
$\lambda_{l},\lambda_{w}>0$ satisfying $\lambda_{l}+\lambda_{w}=1$ and
$\lambda_{w}=O(n^{-\frac{1}{d+q+2}}),$

$K_{\mathcal{D}}\lambda_{w}\leq C_{1}54B2^{d+q}(d+q)^{1/2},$

where $C_{1}$ is a positive constant independent of $(n,J)$. Therefore, for
$d,q,m$ satisfying $2\beta(d+q+1)\geq 3(m+d)+\beta$, $J\gtrsim
n^{\\{3(m+d)+6\beta\\}/\\{4\beta(d+q+2)\\}}$, by (A.16),

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}\right\\}$
$\displaystyle\leq$
$\displaystyle\frac{1}{\lambda_{l}}\mathbb{E}_{\mathcal{S}}\\{\lambda_{l}\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\\}$
$\displaystyle\leq$ $\displaystyle C_{3}n^{-\frac{3}{2(d+q+2)}}(\log
n)^{\frac{2\beta}{m+d}\vee 2},$

and

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\right\\}$
$\displaystyle\leq\frac{1}{\lambda_{w}}\mathbb{E}_{\mathcal{S}}\\{\lambda_{l}\mathbb{E}\|\mathbb{E}_{\eta}\hat{g}(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}+\lambda_{w}d_{\mathcal{F}^{1}_{B}}(P_{X,\hat{g}},P_{X,Y})\\}$
$\displaystyle\leq C_{3}n^{-\frac{1}{2(d+q+2)}}(\log
n)^{\frac{2\beta}{m+d}\vee 2}\rightarrow 0,$

as $n\rightarrow\infty$, where $C_{3}$ is a positive constant independent of
$(n,J)$.

##  Appendix C Supporting lemmas

In this section, we give some supporting lemmas that are used to establish the
upper bound for the terms in the error decomposition and are needed in the
proof of Theorems 1 \- 3.

We bound $\mathcal{E}_{1}$ in Lemma 2.

######  Lemma 2.

Suppose Conditions 1 and 4 hold. Then,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}\|Y-\mathbb{E}_{\eta}\hat{g}(X,\eta)\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}(X_{i},\eta)\right\|^{2}\right\\}$
$\displaystyle\leq$ $\displaystyle
CB_{0}q\sqrt{\frac{H_{\mathcal{G}}S_{\mathcal{G}}\log S_{\mathcal{G}}\\{\log
B_{0}+\log n\\}}{n}},$

where $C$ is a positive constant, and $H_{\mathcal{G}},S_{\mathcal{G}}$ are
the width and depth of the network class $\mathcal{G}$ respectively.

###### Proof.

The proof will be done in two steps.

Step 1. Symmetrization. Let $\\{X^{\prime}_{i},Y^{\prime}_{j}\\}_{i=1}^{n}$ be
$n$ i.i.d. copies of $(X,Y)$, independent of $\\{X_{i},Y_{i}\\}_{i=1}^{n}$.
And let $h_{g}(X,Y)\coloneqq\|Y-\mathbb{E}_{\eta}g(X,\eta)\|^{2}$ for any
$g\in\mathcal{G}$. Then, by standard symmetrization technique,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}\|Y-\mathbb{E}_{\eta}\hat{g}(X,\eta)\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}(X_{i},\eta)\right\|^{2}\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\mathcal{S}}\left[\sup_{g\in\mathcal{G}}\left\\{\mathbb{E}\|Y-\mathbb{E}_{\eta}g(X,\eta)\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}g(X_{i},\eta)\right\|^{2}\right\\}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathcal{S}}\left[\sup_{g\in\mathcal{G}}\left\\{\mathbb{E}h_{g}(X,Y)-\frac{1}{n}\sum_{i=1}^{n}h_{g}(X_{i},Y_{i})\right\\}\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\\{X_{i},Y_{i},X^{\prime}_{i},Y^{\prime}_{i}\\}_{i=1}^{n}}\left[\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\\{h_{g}(X^{\prime}_{i},Y^{\prime}_{i})-h_{g}(X_{i},Y_{i})\\}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\\{X_{i},Y_{i},X^{\prime}_{i},Y^{\prime}_{i},\epsilon_{i}\\}_{i=1}^{n}}\left[\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}\\{h_{g}(X^{\prime}_{i},Y^{\prime}_{i})-h_{g}(X_{i},Y_{i})\\}\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\\{X_{i},Y_{i},X^{\prime}_{i},Y^{\prime}_{i},\epsilon_{i}\\}_{i=1}^{n}}\left[\sup_{g\in\mathcal{G}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}\\{h_{g}(X^{\prime}_{i},Y^{\prime}_{i})-h_{g}(X_{i},Y_{i})\\}\right\\}\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\\{X^{\prime}_{i},Y^{\prime}_{i},\epsilon_{i}\\}_{i=1}^{n}}\left\\{\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}h_{g}(X^{\prime}_{i},Y^{\prime}_{i})\right\\}+\mathbb{E}_{\\{X_{i},Y_{i},\epsilon_{i}\\}_{i=1}^{n}}\left\\{\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}h_{g}(X_{i},Y_{i})\right\\}$
$\displaystyle=$ $\displaystyle
2\mathbb{E}_{\\{X_{i},Y_{i},\epsilon_{i}\\}_{i=1}^{n}}\left\\{\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}h_{g}(X_{i},Y_{i})\right\\},$

where $\\{\epsilon_{j}\\}_{j=1}^{J}$ are independent uniform $\\{\pm
1\\}$-valued Rademacher random variables. We use $\mathcal{R}(\mathcal{H})$ to
denote the Rademacher complexity for $\mathcal{H}$, that is

$\mathcal{R}(\mathcal{H})\coloneqq\mathbb{E}_{\\{X_{i},Y_{i},\epsilon_{i}\\}_{i=1}^{n}}\left[\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}h_{g}(X_{i},Y_{i})\right].$

Step 2. To bound the Rademacher complexity $\mathcal{R}(\mathcal{H})$. We
define a function class
$\mathcal{G}_{|\\{X_{i}\\}_{i=1}^{n}}=\\{\mathbb{E}_{\eta}g(X_{1},\eta),\ldots,\mathbb{E}_{\eta}g(X_{n},\eta):g\in\mathcal{G}\\}$.
Then, under Condition 1,

$\displaystyle\mathcal{R}(\mathcal{H})\leq$ $\displaystyle
4B_{0}\mathbb{E}_{\\{X_{i},Y_{i},\epsilon_{i}\\}_{i=1}^{n}}\left[\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}\\{Y_{i}-\mathbb{E}_{\eta}g(X_{i},\eta)\\}\right]$
$\displaystyle\leq$ $\displaystyle
4B_{0}\mathbb{E}_{\\{Y_{i},\epsilon_{i}\\}_{i=1}^{n}}\left[\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}Y_{i}\right]+4B_{0}\mathbb{E}_{\\{X_{i},\epsilon_{i}\\}_{i=1}^{n}}\left[\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}\mathbb{E}_{\eta}g(X_{i},\eta)\right]$
$\displaystyle\leq$ $\displaystyle
32\mathbb{E}_{\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\left[\inf_{0<\delta<B_{0}/2}\left\\{\delta+\frac{3}{\sqrt{n}}\int_{\delta}^{B_{0}/2}\sqrt{\mathcal{N}(\epsilon,\mathcal{G}_{|\\{X_{i}\\}_{i=1}^{n}},\|\cdot\|_{\infty})}d\epsilon\right\\}\right],$

where the last inequality is by Lemma 12 in Huang et al. (2022). By Condition
4, there exists a vector $\eta_{x}\in\Omega_{\eta}$ such that for any $g$,
$\tilde{g}\in\mathcal{G}$,

$\displaystyle\|\mathbb{E}_{\eta}g(x,\eta)-\mathbb{E}_{\eta}\tilde{g}(x,\eta)\|_{1}\leq\|g(x,\eta_{x})-\tilde{g}(x,\eta)\|_{1}.$

Thus,

$\displaystyle\mathcal{N}(\epsilon,\mathcal{G}_{|\\{X_{i}\\}_{i=1}^{n}},\|\cdot\|_{\infty})\leq\mathcal{N}(\epsilon,\mathcal{G}_{|\\{X_{i},\eta_{xi}\\}_{i=1}^{n}},\|\cdot\|_{\infty}).$

Next, we consider to use $\text{Pdim}(\mathcal{G})$, the pseudo-dimension of
$\mathcal{G}$, to bound
$\mathcal{N}(\epsilon,\mathcal{G}_{|\\{(X_{i},\eta_{i})\\}_{i=1}^{n}},\|\cdot\|_{\infty})$.
When $n\geq\text{Pdim}(\mathcal{G})$, its upper bound can be directly obtained
by Theorem 12.2 in Anthony and Bartlett, (1999). When
$n<\text{Pdim}(\mathcal{G})$,
$\mathcal{G}_{|\\{(X_{i},\eta_{i})\\}_{i=1}^{n}}\subseteq\\{x\in\mathbb{R}^{n}:\|x\|_{\infty}\leq
B_{2}\\}$ can be covered by at most $\lceil 2B_{0}/\epsilon\rceil^{n}$ balls
with radius $\epsilon$ in the distance $\|\cdot\|_{\infty}$. Thus, for any $n$
and $\epsilon>0$,

$\log\mathcal{N}(\epsilon,\mathcal{G}_{|\\{(X_{i},\eta_{i})\\}_{i=1}^{n}},\|\cdot\|_{\infty})\leq\text{Pdim}(\mathcal{G})\log\left[\frac{2eB_{0}n}{\epsilon}\right].$

Since $\mathcal{G}$ is a ReLU network class, according to Bartlett et al.,
(2019), its pseudo-dimension can be bounded by

$H_{\mathcal{G}}S_{\mathcal{G}}\log(S_{\mathcal{G}}/H_{\mathcal{G}})\lesssim\text{Pdim}(\mathcal{G})\lesssim
H_{\mathcal{G}}S_{\mathcal{G}}\log S_{\mathcal{G}},$ (A.17)

where $H_{\mathcal{G}},S_{\mathcal{G}}$ are the width and depth of the network
class $\mathcal{G}$, respectively. As a result, $\mathcal{E}_{1}$ can be
bounded by

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\mathbb{E}\|Y-\mathbb{E}_{\eta}\hat{g}(X,\eta)\|^{2}-\frac{1}{n}\sum_{i=1}^{n}\left\|Y_{i}-\mathbb{E}_{\eta}\hat{g}(X_{i},\eta)\right\|^{2}\right\\}$
$\displaystyle\leq$ $\displaystyle
32\inf_{0<\delta<B_{2}/2}\left(\delta+\frac{3\sqrt{H_{\mathcal{G}}S_{\mathcal{G}}\log
S_{\mathcal{G}}(\mathcal{G})}}{\sqrt{n}}\int_{\delta}^{B_{0}/2}\sqrt{\log[2eB_{2}n/\epsilon]}d\epsilon\right)$
$\displaystyle\leq$ $\displaystyle
32\inf_{0<\delta<B_{0}/2}\left(\delta+\frac{3B_{2}}{2}\sqrt{\frac{H_{\mathcal{G}}S_{\mathcal{G}}\log
S_{\mathcal{G}}(\mathcal{G})\log[2eB_{0}n/\delta]}{n}}\right)$
$\displaystyle\leq$ $\displaystyle
CB_{0}\sqrt{\frac{H_{\mathcal{G}}S_{\mathcal{G}}\log
S_{\mathcal{G}}(\mathcal{G})[\log n+\log B_{0}]}{n}},$

∎

Next, we intend to bound $\mathcal{E}_{2}$, the stochastic error for the
generator. Under Condition 1, for any $g\in\mathcal{G}$ and
$i\in\\{1,\ldots,n\\}$,

$\displaystyle\mathbb{E}_{\mathcal{S}}\left\\{\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i=1}^{n}\left|\left\|Y_{i}-\frac{1}{J}\sum_{j=1}^{J}g(X_{i},\eta_{ij})\right\|^{2}-\|Y_{i}-\mathbb{E}_{\eta}g(X_{i},\eta)\|^{2}\right|\right\\}$
$\displaystyle\leq$ $\displaystyle
4B_{0}\mathbb{E}_{\mathcal{S}}\left[\sup_{g\in\mathcal{G}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left\|\mathbb{E}_{\eta}g(X_{i},\eta)-\frac{1}{J}\sum_{j=1}^{J}g(X_{i},\eta_{ij})\right\|_{1}\right\\}\right].$
(A.18)

Then, we can establish the upper bound of (A.18) in the following lemma.

######  Lemma 3.

Assume Condition 1 hold. Then,

$\displaystyle\mathbb{E}_{\\{X_{i}\\}_{i=1}^{n}\cup\\{\eta_{ij}\\}_{j=1}^{J}}\left[\sup_{g\in\mathcal{G}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left\|\mathbb{E}_{\eta}g(X_{i},\eta)-\frac{1}{J}\sum_{j=1}^{J}g(X_{i},\eta_{ij})\right\|_{1}\right\\}\right]$
$\displaystyle\quad\leq CqB_{0}\sqrt{\frac{H_{\mathcal{G}}S_{\mathcal{G}}\log
S_{\mathcal{G}}\\{\log J+\log(B_{0})\\}}{J}},$

where $C>0$ is a universal constant.

###### Proof.

Recall that $\eta\sim P_{\eta}$ and
$\\{\eta_{ij},i=1,\ldots,n,j=1,\ldots,J\\}$ is $n\times J$ i.i.d. copies of
$\eta$, independent of $\\{X_{i}\\}_{i=1}^{n}$. The proof will be done in
three steps.

Step 1: To decompose the network class $\mathcal{G}$. The network class
$\mathcal{G}$ can be decomposed into the product
$\mathcal{G}_{1}\otimes\mathcal{G}_{2}\otimes\ldots\otimes\mathcal{G}_{q}$,
where $\mathcal{G}_{k}$ satisfies that for
$g=(g_{1},\ldots,g_{q})^{\top}\in\mathcal{G}$,
$g_{k}\in\mathcal{G}_{k}\subseteq\mathbb{R}$ for $k=1,\ldots,q$. Note that

$\displaystyle\mathbb{E}_{\\{X_{i}\\}_{i=1}^{n}\cup\\{\eta_{ij}\\}_{j=1}^{J}}\left[\sup_{g\in\mathcal{G}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\left\|\mathbb{E}_{\eta}g(X_{i},\eta)-\frac{1}{J}\sum_{j=1}^{J}g(X_{i},\eta_{ij})\right\|_{1}\right\\}\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\\{X_{i}\\}_{i=1}^{n}\cup\\{\eta_{ij}\\}_{j=1}^{J}}\left[\frac{1}{n}\sum_{i=1}^{n}\left\\{\sup_{g\in\mathcal{G}}\left\|\mathbb{E}_{\eta}g(X_{i},\eta)-\frac{1}{J}\sum_{j=1}^{J}g(X_{i},\eta_{ij})\right\|_{1}\right\\}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\\{X\\}\cup\\{\eta_{j}\\}_{j=1}^{J}}\left\\{\sup_{g\in\mathcal{G}}\left\|\mathbb{E}_{\eta}g(X,\eta)-\frac{1}{J}\sum_{j=1}^{J}g(X,\eta_{j})\right\|_{1}\right\\}$
$\displaystyle\leq$
$\displaystyle\sum_{k=1}^{q}\mathbb{E}_{\\{X\\}\cup\\{\eta_{j}\\}_{j=1}^{J}}\left\\{\sup_{g_{k}\in\mathcal{G}_{k}}\left|\mathbb{E}_{\eta}g_{k}(X,\eta)-\frac{1}{J}\sum_{j=1}^{J}g_{k}(X,\eta_{j})\right|\right\\}$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{q}\mathbb{E}_{X}\left[\mathbb{E}_{\\{\eta_{j}\\}_{j=1}^{J}}\left\\{\sup_{g_{k}\in\mathcal{G}_{k}}\left|\mathbb{E}_{\eta}g_{k}(X,\eta)-\frac{1}{J}\sum_{j=1}^{J}g_{k}(X,\eta_{j})\right|\right\\}\right].$
(A.19)

Step 2: To bound
$\mathbb{E}_{\\{\eta_{j}\\}_{j=1}^{J}}\\{\sup_{g_{k}\in\mathcal{G}_{k}}|\mathbb{E}_{\eta}g_{k}(X,\eta)-\frac{1}{J}\sum_{j=1}^{J}g_{k}(X,\eta_{j})|\\}$.
Let $\\{\eta^{\prime}_{j}\\}_{j=1}^{J}$ be $J$ i.i.d. copies of $\eta$,
independent of $\\{\eta_{j}\\}_{j=1}^{J}$ and $X$. It then follows from the
standard symmetrization technique that given $X$, for $k=1,\ldots,q$,

$\displaystyle\mathbb{E}_{\\{\eta_{j}\\}_{j=1}^{J}}\left\\{\sup_{g_{k}\in\mathcal{G}_{k}}\left|\mathbb{E}_{\eta}g_{k}(X,\eta)-\frac{1}{J}\sum_{j=1}^{J}g_{k}(X,\eta_{j})\right|\right\\}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\\{\eta_{j}\\}_{j=1}^{J}}\left[\sup_{g_{k}\in\mathcal{G}_{k}}\mathbb{E}_{\\{\eta^{\prime}_{j}\\}_{j=1}^{J}}\left|\frac{1}{J}\sum_{j=1}^{J}\\{g_{k}(X,\eta_{j})-g_{k}(X,\eta^{\prime}_{j})\\}\right|\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\\{\eta_{j},\eta^{\prime}_{j},\epsilon_{j}\\}_{j=1}^{J}}\left[\sup_{g_{k}\in\mathcal{G}_{k}}\left|\frac{1}{J}\sum_{j=1}^{J}\epsilon_{j}\\{g_{k}(X,\eta_{j})-g_{k}(X,\eta^{\prime}_{j})\\}\right|\right]$
$\displaystyle\leq$ $\displaystyle
2\mathbb{E}_{\\{\eta_{j},\epsilon_{j}\\}_{j=1}^{J}}\left\\{\sup_{g_{k}\in\mathcal{G}_{k}}\left|\frac{1}{J}\sum_{j=1}^{J}\epsilon_{j}g_{k}(X,\eta_{j})\right|\right\\}$
$\displaystyle=$ $\displaystyle 2\mathcal{R}(\mathcal{G}_{k}|X),$

where $\mathcal{R}(\mathcal{G}_{k}|X)$ is the conditional Rademacher
complexity of $\mathcal{G}_{k}$ given $X$.

Step 3: To bound the conditional Rademacher complexity
$\mathcal{R}(\mathcal{G}_{k}|X)$. For any
$g_{k},\tilde{g}_{k}\in\mathcal{G}_{k}$,

$\displaystyle\left|\frac{1}{J}\sum_{j=1}^{J}\epsilon_{j}g_{k}(X,\eta_{j})\right|-\left|\frac{1}{J}\sum_{j=1}^{J}\epsilon_{j}\tilde{g}_{k}(X,\eta_{j})\right|$
$\displaystyle\leq\left|\frac{1}{J}\sum_{j=1}^{J}\epsilon_{j}\\{g_{k}(X,\eta_{j})-\tilde{g}_{k}(X,\eta_{j})\\}\right|$
$\displaystyle\leq\max_{1\leq j\leq
J}|g_{k}(X,\eta_{j})-\tilde{g}_{k}(X,\eta_{j})|.$

We use $e_{X,J}(g_{k},\tilde{g}_{k})$ to denote the distance between $g_{k}$
and $\tilde{g}_{k}$ with respect to $\\{\eta_{j}\\}_{j=1}^{J}$ and $X$ defined
as

$\displaystyle
e_{X,J}(g_{k},\tilde{g}_{k})=\|g_{k}-\tilde{g}_{k}\|_{L^{\infty}(\\{(X,\eta_{j})\\}_{j=1}^{J})}.$

For any $\delta>0$, we define $\mathcal{G}_{k,\delta}$ to be a covering set of
$\mathcal{G}_{k}$ with radius $\delta$ with respect to $e_{X,J}(\cdot)$ and
$\mathcal{N}(\delta,\mathcal{G}_{k},e_{X,J})$ to be the $\delta$-covering
number of $\mathcal{G}_{k}$ with respect to the distance $e_{X,J}(\cdot)$.
The, by the triangle inequality and Lemma B.4 in Zhou et al., (2022), we have

$\displaystyle\mathcal{R}(\mathcal{G}_{k}|X)$
$\displaystyle\leq\delta+\mathbb{E}_{\\{\eta_{j},\epsilon_{j}\\}_{j=1}^{J}}\left\\{\sup_{g_{k}\in\mathcal{G}_{k,\delta}}\left|\frac{1}{J}\sum_{j=1}^{J}\epsilon_{j}g_{k}(X,\eta_{j})\right|\right\\}$
$\displaystyle\leq\delta+\frac{C_{1}}{J}\mathbb{E}_{\\{\eta_{j}\\}_{j=1}^{J}}\left[\left\\{\log\mathcal{N}(\delta,\mathcal{G}_{k},e_{X,J})\right\\}^{1/2}\left\\{\max_{g_{k}\in\mathcal{G}_{k,\delta}}\sum_{j=1}^{J}g_{k}^{2}(X,\eta_{j})\right\\}^{1/2}\right]$
$\displaystyle\leq\delta+\frac{C_{1}B_{0}}{\sqrt{J}}\mathbb{E}_{\\{\eta_{j}\\}_{j=1}^{J}}\left[\left\\{\log\mathcal{N}(\delta,\mathcal{G}_{k},e_{X,J})\right\\}^{1/2}\right],$

where $C_{1}>0$ is a constant and the third inequality holds since the
probability measure of $(X,g(X,\eta))$ is supported on
$\Omega\subseteq[-B_{0},B_{0}]^{d+q}\subseteq\mathbb{R}^{d+q}$ for any
$g\in\mathcal{G}$. It follows from Theorem 12.2 in Anthony and Bartlett,
(1999) that for $J\geq\text{Pdim}(\mathcal{G}_{k})$,

$\log\mathcal{N}(\delta,\mathcal{G}_{k},e_{X,J})=\log\mathcal{N}(\delta,\mathcal{G}_{k|\\{(X,\eta_{j})\\}_{j=1}^{J}},\|\cdot\|_{\infty})\leq\text{Pdim}(\mathcal{G}_{k})\log\left\\{\frac{2eJB_{0}}{\delta\text{Pdim}(\mathcal{G}_{k})}\right\\},$

where
$\mathcal{G}_{k|\\{(X,\eta_{j})\\}_{j=1}^{J}}=\\{(g_{k}(\eta_{1},X),\ldots,g_{k}(X,\eta_{j})):g_{k}\in\mathcal{G}_{k}\\}$
and
$\mathcal{N}(\delta,\mathcal{G}_{k|\\{(X,\eta_{j})\\}_{j=1}^{J}},\|\cdot\|_{\infty})$
is the $\delta$-covering number of
$\mathcal{G}_{k|\\{(X,\eta_{j})\\}_{j=1}^{J}}\subseteq\mathbb{R}^{J}$ with
respect to the distance $\|\cdot\|_{\infty}$. When
$J<\text{Pdim}(\mathcal{G}_{k})$,
$\mathcal{G}_{k|\\{(X,\eta_{j})\\}_{j=1}^{J}}\subseteq\\{x\in\mathbb{R}^{J}:\|x\|_{\infty}\leq
B_{0}\\}$ can be covered by at most $\lceil 2B_{0}/\delta\rceil^{J}$ balls
with radius $\delta$ with respect to the distance $\|\cdot\|_{\infty}$. This
means that for any $J$,

$\log\mathcal{N}(\delta,\mathcal{G}_{k},e_{X,J})\leq\text{Pdim}(\mathcal{G}_{k})\log\left(\frac{2eJB_{0}}{\delta}\right).$

Therefore,

$\mathcal{R}(\mathcal{G}_{k}|X)\leq\delta+\frac{C_{1}B_{0}}{\sqrt{J}}\log\left(\frac{2eJB_{0}}{\delta}\right)\text{Pdim}(\mathcal{G}_{k}).$

Moreover, Bartlett et al., (2019) gives the following result for the pseudo-
dimension of a ReLU network class:

$\text{Pdim}(\mathcal{G}_{k})\lesssim H_{\mathcal{G}}S_{\mathcal{G}}\log
S_{\mathcal{G}},$

where $H_{\mathcal{G}},S_{\mathcal{G}}$ are the width and depth of the network
class $\mathcal{G}$, respectively. By letting $\delta=1/J$, we have

$\displaystyle\mathbb{E}_{\\{\eta_{j}\\}_{j=1}^{J}}\left\\{\sup_{g_{k}\in\mathcal{G}_{k}}\left|\mathbb{E}_{\eta}g_{k}(X,\eta)-\frac{1}{J}\sum_{j=1}^{J}g_{k}(X,\eta_{j})\right|\right\\}\lesssim
B_{0}\sqrt{\frac{H_{\mathcal{G}}S_{\mathcal{G}}\log S_{\mathcal{G}}\\{\log
J+\log(B_{0})\\}}{J}}.$

Thus, we have proved Lemma 3 in view of (C).

∎

We next can bound $\mathcal{E}_{3}$ by combining Lemmas 4 and 5.

######  Lemma 4.

Suppose Condition 1 holds. Then,

$\mathcal{E}_{3}\leq(K_{\mathcal{D}}\lambda_{w}+2B_{0}q^{2}\lambda_{l})\inf_{g\in\mathcal{G}}\sup_{f\in\mathcal{F}^{1}_{B}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\right\\},$
(A.20)

where $K_{\mathcal{D}}$ is the Lipschitz constant defined in ND 1.

###### Proof.

Let $Y=(Y^{(1)},\ldots,Y^{(q)})^{\top}\in\mathbb{R}^{q}$. Under Condition 1,
for any $g=(g^{(1)},\ldots,g^{(q)})^{\top}\in\mathcal{G}$,

$\displaystyle\mathbb{E}\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}$
$\displaystyle\leq\mathbb{E}\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}(Y|X)\|_{1}^{2}$
$\displaystyle\leq
2B_{0}q\mathbb{E}\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}(Y|X)\|_{1}$
$\displaystyle\leq
2B_{0}q^{2}\max_{k=1,\ldots,q}\mathbb{E}\left|\mathbb{E}_{\eta}g^{(k)}(X,\eta)-\mathbb{E}(Y^{(k)}|X)\right|.$

Let $F_{g^{(k)}(X,\eta)|X}(\cdot)$ and $F_{Y^{(k)}|X}(\cdot)$ denote the
conditional cumulative distribution of $g^{(k)}(X,\eta)$ and $Y^{(k)}$,
respectively. For any $k\in\\{1,\ldots,q\\}$ and
$g=(g^{(1)},\ldots,g^{(q)})^{\top}\in\mathcal{G}$, by the definition of
Wasserstein distance,

$\displaystyle\mathbb{E}\left|\mathbb{E}_{\eta}g^{(k)}(X,\eta)-\mathbb{E}(Y^{(k)}|X)\right|$
$\displaystyle\leq$
$\displaystyle\mathbb{E}\mathcal{W}_{1}(P_{g^{(k)}(X,\eta)|X},P_{Y^{(k)}|X})\leq\mathcal{W}_{1}(P_{X,g(X,\eta)},P_{X,Y})=\sup_{f\in\mathcal{F}^{1}_{B}}\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\\},$

where the last equality holds for any $g\in\mathcal{G}$, and
$\mathcal{F}^{1}_{B}=\\{f:\mathbb{R}^{d+q}\mapsto\mathbb{R},|f(z_{1})-f(z_{2})|\leq\|z_{1}-z_{2}\|,z_{1},z_{2}\in\mathbb{R}^{d+q},\|f\|_{\infty}\leq
B\\}$. Hence, for any $g\in\mathcal{G}$,

$\displaystyle\mathbb{E}\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\|^{2}\leq
2B_{0}q^{2}\sup_{f\in\mathcal{F}^{1}_{B}}\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\\}.$
(A.21)

Under Conditions 1, ND 1 and NG 1,

$\displaystyle\sup_{f\in\mathcal{D}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\right\\}\leq
K_{\mathcal{D}}\sup_{f\in\mathcal{F}^{1}_{B}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\right\\}.$

Consequently, together with (A.21), we have

$\displaystyle\mathcal{E}_{3}=$
$\displaystyle\inf_{g\in\mathcal{G}}\left[\lambda_{l}\mathbb{E}\left\|\mathbb{E}_{\eta}g(X,\eta)-\mathbb{E}_{\eta}g^{*}(X,\eta)\right\|^{2}+\lambda_{w}\sup_{f\in\mathcal{D}}\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\\}\right]$
$\displaystyle\leq$
$\displaystyle(K_{\mathcal{D}}\lambda_{w}+2B_{0}q^{2}\lambda_{l})\inf_{g\in\mathcal{G}}\sup_{f\in\mathcal{F}^{1}_{B}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\right\\}.$

∎

######  Lemma 5.

For $g^{*}=(g^{*}_{1},\ldots,g^{*}_{q})^{\top}$, assume that
$g^{*}_{k}\in\mathcal{H}^{\beta}([-B_{0},B_{0}]^{m+d},B_{1}),k=1,\ldots,q$,
and the network class $\mathcal{G}$ satisfies NG 1. Then,

$\displaystyle\inf_{g\in\mathcal{G}}\sup_{f\in\mathcal{F}^{1}_{B}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\right\\}$
$\displaystyle\leq$ $\displaystyle
19B_{1}(\lfloor\beta\rfloor+1)^{2}(m+d)^{\lfloor\beta\rfloor+1/2+(\beta\vee
1)/2}(\bar{W}\bar{H})^{-2\beta/(m+d)}.$

###### Proof.

By the definition of $g^{*}$ in (1) and the Lipschitz continuity of
$f\in\mathcal{F}^{1}_{B}$, we have

$\displaystyle\inf_{g\in\mathcal{G}}\sup_{f\in\mathcal{F}^{1}_{B}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\right\\}$
$\displaystyle=\inf_{g\in\mathcal{G}}\sup_{f\in\mathcal{F}^{1}_{B}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,g^{*}(X,\eta))\right\\}$
$\displaystyle\leq\inf_{g\in\mathcal{G}}\left\\{\mathbb{E}\|g(X,\eta)-g^{*}(X,\eta)\|\right\\}$
$\displaystyle\leq\sqrt{m+d}\inf_{g\in\mathcal{G}}\|g-g^{*}\|_{L^{\infty}([-B_{0},B_{0}]^{m+d})}.$

For any $k\in\\{1,\ldots,q\\}$ and
$g^{*}_{k}\in\mathcal{H}^{\beta}([-B_{0},B_{0}]^{m+d},B_{1})$, we define a new
function

$\displaystyle\tilde{g}^{0}_{k}(x):=g^{*}_{k}\left(2B_{0}x-B_{0}\right)\in\mathcal{H}^{\beta}([0,1]^{m+d},B_{1}).$

Corollary 3.1 in Jiao et al., (2023) tells that there exists a ReLU network
function $\tilde{\phi}_{k}$, satisfying NG 1, such that for any $k$,

$\displaystyle\|\tilde{\phi}_{k}-\tilde{g}^{0}_{k}\|_{L^{\infty}([0,1]^{m+d})}\leq
19B_{1}(\lfloor\beta\rfloor+1)^{2}(m+d)^{\lfloor\beta\rfloor+(\beta\vee
1)/2}(WH)^{-2\beta/(m+d)}.$

Then, we define

$\displaystyle\phi(x):=\left(\tilde{\phi}_{1}\left(\frac{1}{2B_{0}}x-\frac{1}{2}\right),\ldots,\tilde{\phi}_{q}\left(\frac{1}{2B_{0}}x-\frac{1}{2}\right)\right)^{\top}.$

Hence, by Remark 14 in Nakada and Imaizumi, (2020), $\phi\in\mathcal{G}$ with
width $W_{\mathcal{G}}$ and $H_{\mathcal{G}}$ satisfying NG 1. Then,

$\displaystyle\|\phi-g^{*}\|_{L^{\infty}([-B_{0},B_{0}]^{m+d})}$
$\displaystyle\leq\max_{1\leq k\leq
q}\left\|\tilde{\phi}_{k}\left(\frac{1}{2B_{0}}x-\frac{1}{2}\right)-\tilde{g}^{0}_{k}\left(\frac{1}{2B_{0}}x-\frac{1}{2}\right)\right\|_{L^{\infty}([-B_{0},B_{0}]^{m+d})}$
$\displaystyle=\max_{1\leq k\leq
q}\|\tilde{\phi}_{k}-\tilde{g}^{0}_{k}\|_{L^{\infty}([0,1]^{m+d})}$
$\displaystyle\leq
19B_{1}(\lfloor\beta\rfloor+1)^{2}(m+d)^{\lfloor\beta\rfloor+(\beta\vee
1)/2}(WH)^{-2\beta/(m+d)}.$

Therefore,

$\displaystyle\inf_{g\in\mathcal{G}}\sup_{f\in\mathcal{F}^{1}_{B}}\left\\{\mathbb{E}f(X,g(X,\eta))-\mathbb{E}f(X,Y)\right\\}\leq\sqrt{m+d}\|\phi-g^{*}\|_{L^{\infty}([-B_{0},B_{0}]^{m+d})}$
$\displaystyle\leq$ $\displaystyle
19B_{1}(\lfloor\beta\rfloor+1)^{2}(m+d)^{\lfloor\beta\rfloor+1/2+(\beta\vee
1)/2}(WH)^{-2\beta/(m+d)}.$

∎

Next, we can bound $\mathcal{E}_{4}$ by Lemmas 6.

######  Lemma 6.

Suppose Condition 1 holds. Then, there exists a network function $\tilde{f}$
satisfying the width and depth in ND 1 such that on the domain
$[-B_{0},B_{0}]^{d+q}$,

$\sup_{z\in[-B_{0},B_{0}]^{d+q}}|\phi(z)-f(z)|\leq 2^{d+q}(d+q)BW^{-1}.$

###### Proof.

The proof is carried out in two steps. We begin with approximating $f$ by a
sum-product combination of univariate piecewise-linear functions. And then we
use a neural network to approximate it.

_Step 1:_ To approximate $f$ by a piecewise-linear function. Let $N$ be a
positive integer. First, we find a grid of $(N+1)^{d+q}$ functions
$\phi_{\mathbf{m}}(\mathbf{x})$ such that

$\displaystyle\sum_{\mathbf{m}\in\\{0,1,\cdots,N\\}^{d+q}}\phi_{\mathbf{m}}(\mathbf{x})\equiv
1,\quad\mathbf{x}\in[-B_{0},B_{0}]^{d+q}.$

To address this issue, for any $\mathbf{x}=(x_{1},\cdots,x_{d+q})^{\top}$, we
define

$\displaystyle\phi_{\mathbf{m}}(\mathbf{x})=\prod_{k=1}^{d+q}\psi\left(3N\left(x_{k}-\frac{m_{k}}{N}\right)\right),$
(A.22)

where $\mathbf{m}=(m_{1},\cdots,m_{d+q})\in\\{0,1,\cdots,N\\}^{d+q}$ and

$\displaystyle\psi(x)=\begin{cases}0,&|x|>2,\\\ 2-|x|,&1\geq|x|\leq 2,\\\
1,&|x|<1.\end{cases}$ (A.23)

Note that for arbitrary $\mathbf{m}$,
$\text{supp}\,\phi_{\mathbf{m}}\subset\\{\mathbf{x}:|x_{k}-m_{k}/N|\leq
1/N\\}$. Then, we construct the function $f_{1}$ to approximate $f$ defined as

$\displaystyle
f_{1}(\mathbf{x})=\sum_{\mathbf{m}\in\\{0,1,\cdots,N\\}^{d+q}}a_{\mathbf{m}}\phi_{\mathbf{m}}(\mathbf{x}),$

where $a_{\mathbf{m}}=f(\mathbf{m}/N)$. Then, by the properties of
$\phi_{\mathbf{m}}$ and the Lipschitz property of $f$, the approximation error
can be bounded by

$\displaystyle\|f-f_{1}\|_{L^{\infty}([-B_{0},B_{0}]^{d+q})}$
$\displaystyle=\|\sum_{\mathbf{m}\in\\{0,1,\cdots,N\\}^{d+q}}\phi_{\mathbf{m}}(\mathbf{x})(f(\mathbf{x})-a_{\mathbf{m}})\|_{L^{\infty}([-B_{0},B_{0}]^{d+q})}$
$\displaystyle\leq\sum_{\mathbf{m}:|x_{k}-m_{k}/N|\leq 1/N,\,\forall
k}|f(\mathbf{x})-a_{\mathbf{m}}|$ $\displaystyle\leq
2^{d+q}\max_{\mathbf{m}:|x_{k}-m_{k}/N|\leq 1/N,\,\forall
k}|f(\mathbf{x})-a_{\mathbf{m}}|$ $\displaystyle\leq
2^{d+q}\max_{\mathbf{m}:|x_{k}-m_{k}/N|\leq 1/N,\,\forall
k}\|\mathbf{x}-\mathbf{m}/N\|_{2}$
$\displaystyle\leq\frac{2^{d+q}\sqrt{d+q}}{N}.$

_Step 2:_ To approximate $f_{1}$ by a neural network. First, we construct a
ReLU activated $\psi$ satisfying (A.23):

$\displaystyle\psi(x)=\sigma(\sigma(x+2)-\sigma(x+1)+\sigma(2-x)-\sigma(1-x)-1).$
(A.24)

According to Hon and Yang, (2022), for an almost everywhere differentiable
function $g:[0,1]^{d+q}\rightarrow\mathbb{R}$, we can define its sobolev norm
as

$\|g\|_{\mathcal{W}^{1,\infty}([0,1]^{d+q})}=\max_{\mathbf{n}:|\mathbf{n}|_{1}=1}\sup_{\mathbf{x}\in[0,1]^{d+q}}|D^{\mathbf{n}}g(\mathbf{x})|,$

where $\mathbf{n}=(n_{1},\cdots,n_{d+q})\in\\{0,1\\}^{d+q}$ and
$D^{\mathbf{n}}g$ is the respective weak derivative. Then, we define

$\displaystyle\tilde{f}_{\mathbf{m}}(\mathbf{x})=$
$\displaystyle\varphi\left(\psi(3Nx_{1}-3m_{1}),\cdots,\psi(3Nx_{d+q}-3m_{d+q})\right),$
$\displaystyle\tilde{f}(\mathbf{x})=$
$\displaystyle\sum_{\mathbf{m}\in\\{0,1,\cdots,N\\}^{d+q}}a_{\mathbf{m}}\tilde{f}_{\mathbf{m}}(\mathbf{x}),$
(A.25)

where $\varphi$ is a ReLU neural network with width $9(N_{1}+1)+d+q-1$ and
depth $14(d+q)(d+q-1)L$, and
$\|\varphi\|_{\mathcal{W}^{1,\infty}([0,1]^{d+q})}\leq 18$. By Lemma 7, we can
implement $\tilde{f}$ by a neural network satisfying ND 1. Then, the
approximation error of $f_{1}$ can be bounded by

$\displaystyle\mathop{\sup}_{\mathbf{x}\in[-B_{0},B_{0}]^{d+q}}|\tilde{f}(\mathbf{x})-f_{1}(\mathbf{x})|$
$\displaystyle=\mathop{\sup}_{\mathbf{x}\in[-B_{0},B_{0}]^{d+q}}|\sum_{\mathbf{m}}a_{\mathbf{m}}\\{\tilde{f}_{\mathbf{m}}(\mathbf{x})-\phi_{\mathbf{m}}(\mathbf{x})\\}|$
$\displaystyle\leq
B(N+1)^{d+q}\mathop{\sup}_{\mathbf{x}\in[-B_{0},B_{0}]^{d+q}}\max_{\mathbf{m}}|\tilde{f}_{\mathbf{m}}(\mathbf{x})-\phi_{\mathbf{m}}(\mathbf{x})|$
$\displaystyle\leq B(N+1)^{d+q}\|\varphi-t_{1}t_{2}\cdots
t_{d+q}\|_{\mathcal{W}^{1,\infty}([0,1]^{d+q})}$ $\displaystyle\leq
B(N+1)^{d+q}10(d+q-1)(N+1)^{-7(d+q)}$
$\displaystyle=10B(d+q-1)(N+1)^{-6(d+q)},$

where the third inequality is by Lemma 3.5 in Hon and Yang, (2022). Therefore,

$\displaystyle\mathop{\sup}_{\mathbf{x}\in[-B_{0},B_{0}]^{d+q}}|\tilde{f}(\mathbf{x})-f(\mathbf{x})|$
$\displaystyle\leq\mathop{\sup}_{\mathbf{x}\in[-B_{0},B_{0}]^{d+q}}|\tilde{f}(\mathbf{x})-f_{1}(\mathbf{x})|+\|f-f_{1}\|_{L^{\infty}([-B_{0},B_{0}]^{d+q})}$
$\displaystyle\leq 2^{d+q}(d+q)^{1/2}N^{-1}+10B(d+q-1)(N+1)^{-6(d+q)}$
$\displaystyle\leq 2^{d+q}(d+q)BN^{-1}.$

∎

######  Lemma 7.

Suppose Condition 1 holds. Then, for the network function $\tilde{f}$ definied
in (C), its Lipschitz constant satisfies $Lip(\tilde{f})\leq
54B2^{d+q}(d+q)^{1/2}N^{2}$.

###### Proof.

To bound the Lipschitz constant of $\tilde{f}$, we need to find a positive
constant $C$ such that

$\max_{\mathbf{n}:|\mathbf{n}|_{1}=1}\sup_{\mathbf{x}\in[0,1]^{d+q}}|D^{\mathbf{n}}\tilde{f}(\mathbf{x})|\leq
C,$

as it is a sufficient condition for $\text{Lip}(\tilde{f})\leq\sqrt{d+q}C$. By
the definitions of $\tilde{f}$ and $\tilde{f}_{\mathbf{m}}$, we can decompose
the first-order partial derivative of $\tilde{f}$ over $x_{k_{0}}$ as

$\displaystyle\left|\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}(\mathbf{x})\right|=$
$\displaystyle\left|\sum_{\mathbf{m}\in\\{0,1,\cdots,N\\}^{d+q}}a_{\mathbf{m}}\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|$ $\displaystyle\leq$
$\displaystyle\left|\sum_{\mathbf{m}:\forall k,\atop
x_{k}\in\\{t:|t-m_{k}/N|\leq 1/N\\}}a_{\mathbf{m}}\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|+\left|\sum_{\mathbf{m}:\exists
k_{1},\atop x_{k_{1}}\notin\\{t:|t-m_{k_{1}}/N|\leq
1/N\\}}a_{\mathbf{m}}\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|$ $\displaystyle\leq$
$\displaystyle\left|\sum_{\mathbf{m}:\forall k,\atop
x_{k}\in\\{t:|t-m_{k}/N|\leq 1/N\\}}a_{\mathbf{m}}\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|+\left|\sum_{\mathbf{m}:\exists
k_{1}\neq k_{0},\atop x_{k_{1}}\notin\\{t:|t-m_{k_{1}}/N|\leq
1/N\\}}a_{\mathbf{m}}\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|$
$\displaystyle\quad+\left|\sum_{\mathbf{m}:x_{k_{0}}\notin\\{t:|t-m_{k_{0}}/N|\leq
1/N\\},\atop\forall k_{1}\neq k_{0},x_{k_{1}}\in\\{t:|t-m_{k_{1}}/N|\leq
1/N\\}}a_{\mathbf{m}}\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|$ $\displaystyle\coloneqq$
$\displaystyle E_{1}+E_{2}+E_{3}.$

We shall next bound $E_{1}$, $E_{2}$ and $E_{3}$, respectively. For $E_{1}$,
by the definition of $\tilde{f}_{\mathbf{m}}$ and the properties of $\varphi$,

$\displaystyle E_{1}$ $\displaystyle\leq 2^{d+q}B\max_{\mathbf{m}:\forall
k,\atop x_{k}\in\\{t:|t-m_{k}/N|\leq 1/N\\}}\left|\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|\leq 2^{d+q}B\cdot
3N\|\varphi\|_{\mathcal{W}^{1,\infty}([0,1]^{d+q})}\leq 54B2^{d+q}N,$

where the second inequality holds according to the chain derivative rule.

Next, we consider $E_{2}$. By the definition of $\phi_{\mathbf{m}}$, for any
$\mathbf{m}\in\\{\mathbf{m}:\exists k_{1}\neq
k_{0},x_{k_{1}}\notin\\{t:|t-m_{k_{1}}/N|\leq 1/N\\}\\}$, we have
$\partial\phi_{\mathbf{m}}/\partial x_{k_{0}}=0$. Thus, similar to $E_{1}$, we
have

$\displaystyle E_{2}$ $\displaystyle\leq B\sum_{\mathbf{m}:\exists k_{1}\neq
k_{0},x_{k_{1}}\notin\\{t:|t-m_{k_{1}}/N|\leq
1/N\\}}\left|\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|$ $\displaystyle\leq
B(N+1)^{d+q}\max_{\mathbf{m}:\exists k_{1}\neq
k_{0},x_{k_{1}}\notin\\{t:|t-m_{k_{1}}/N|\leq
1/N\\}}\left|\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|$
$\displaystyle=B(N+1)^{d+q}\cdot\max_{\mathbf{m}:\exists k_{1}\neq
k_{0},x_{k_{1}}\notin\\{t:|t-m_{k_{1}}/N|\leq
1/N\\}}\left|\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})-\frac{\partial}{\partial
x_{k_{0}}}\phi_{\mathbf{m}}(\mathbf{x})\right|$ $\displaystyle\leq
B(N+1)^{d+q}\cdot 3N\|\varphi-t_{1}t_{2}\cdots
t_{d+q}\|_{\mathcal{W}^{1,\infty}([0,1]^{d+q})}$ $\displaystyle\leq
30B(d+q-1)(N+1)^{-5(d+q)}.$

Third, $E_{3}$ satisfies

$\displaystyle E_{3}$ $\displaystyle\leq
B(N-1)2^{d+q-1}\max_{\mathbf{m}:x_{k_{0}}\notin\\{t:|t-m_{k_{0}}/N|\leq
1/N\\},\forall k_{1}\neq k_{0},x_{k_{1}}\in\\{t:|t-m_{k_{1}}/N|\leq
1/N\\}}\left|\frac{\partial}{\partial
x_{k_{0}}}\tilde{f}_{\mathbf{m}}(\mathbf{x})\right|$ $\displaystyle\leq
B(N-1)2^{d+q-1}\cdot 3N\|\varphi\|_{\mathcal{W}^{1,\infty}([0,1]^{d+q})}$
$\displaystyle\leq 542^{d+q-1}BN^{2}.$

As a result, we can bound the sobolev norm of $\tilde{f}$ by

$\displaystyle\max_{\mathbf{n}:|\mathbf{n}|_{1}=1}\sup_{\mathbf{x}\in[-B_{0},B_{0}]^{d+q}}|D^{\mathbf{n}}\tilde{f}(\mathbf{x})|\leq
54B2^{d+q}N^{2}.$

Therefore, on the domain $[-B_{0},B_{0}]^{d+q}$, $\text{Lip}(\tilde{f})\leq
54B2^{d+q}(d+q)^{1/2}N^{2}$.  ∎

In the following, we shall establish the bound of $\mathcal{E}_{5}$ in Lemma
8.

######  Lemma 8.

Suppose $\sup_{f\in\mathcal{D}}\|f\|_{\infty}\leq B_{2}$ for a constant
$B_{2}$ and the pseudo-dimension of $\mathcal{D}$ satisfies
$\text{Pdim}(\mathcal{D})<\infty$, then

$\displaystyle\mathbb{E}_{\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\left[\sup_{f\in\mathcal{D}}\left\\{\mathbb{E}f(X,Y)-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},Y_{i})\right\\}\right]\leq
CB_{2}\sqrt{\frac{H_{\mathcal{D}}S_{\mathcal{D}}\log(S_{\mathcal{D}})[\log
n+\log B_{2}]}{n}},$

where $C>0$ is a universal constant, $H_{\mathcal{D}},S_{\mathcal{D}}$ are the
width and depth of the network class $\mathcal{D}$, respectively.

###### Proof.

Define
$\mathcal{D}_{|\\{(X_{i},Y_{i})\\}_{i=1}^{n}}=\\{(f(X_{1},Y_{1}),\ldots,f(X_{n},Y_{n})):f\in\mathcal{D}\\}$,
where $\\{(X_{i},Y_{i})\\}_{i=1}^{n}$ are $n$ i.i.d. samples from $P_{X,Y}$.
Let
$\mathcal{N}(\epsilon,\mathcal{D}_{|\\{(X_{i},Y_{i})\\}_{i=1}^{n}},\|\cdot\|_{\infty})$
be the $\epsilon$-covering number of
$\mathcal{D}_{|\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\subseteq\mathbb{R}^{n}$ with
respect to the distance $\|\cdot\|_{\infty}$.

It follows from Lemma 12 in Huang et al., (2022) that

$\displaystyle\mathbb{E}_{\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\left[\sup_{f\in\mathcal{D}}\left\\{\mathbb{E}f(X,Y)-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},Y_{i})\right\\}\right]$
$\displaystyle\quad\leq
8\mathbb{E}_{\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\left[\inf_{0<\delta<B_{2}/2}\left\\{\delta+\frac{3}{\sqrt{n}}\int_{\delta}^{B_{2}/2}\sqrt{\mathcal{N}(\epsilon,\mathcal{D}_{|\\{(X_{i},Y_{i})\\}_{i=1}^{n}},\|\cdot\|_{\infty})}d\epsilon\right\\}\right].$

It follows from similar arguments as in the proof of Lemma 2 that
$\mathcal{E}_{5}$ can be bounded by

$\displaystyle\mathbb{E}_{\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\left[\sup_{f\in\mathcal{D}}\left\\{\mathbb{E}f(X,Y)-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},Y_{i})\right\\}\right]$
$\displaystyle\leq$ $\displaystyle
CB_{2}\sqrt{\frac{H_{\mathcal{D}}S_{\mathcal{D}}\log
S_{\mathcal{D}}(\mathcal{D})[\log n+\log B_{2}]}{n}},$ (A.26)

where $C$ is a universal positive constant.  ∎

We will bound $\mathcal{E}_{6}$ in Lemma 9.

######  Lemma 9.

Suppose $\sup_{f\in\mathcal{D}}\|f\|_{\infty}\leq B_{2}$ and the pseudo-
dimension of $\mathcal{G}$ and $\mathcal{D}$ satisfies
$\text{Pdim}(\mathcal{G})<\infty,\text{Pdim}(\mathcal{D})<\infty$. Let
$(W_{\mathcal{G}},H_{\mathcal{G}})$ be the width and depth of the network in
$\mathcal{G}$ and $(W_{\mathcal{D}},H_{\mathcal{D}})$ be the width and depth
of the network in $\mathcal{D}$. Let $S_{\mathcal{G},\mathcal{D}}$ denote the
size of the network in
$\mathcal{NN}(m+d,1,2W_{\mathcal{G}}+2W_{\mathcal{D}}+2d,H_{\mathcal{G}}+H_{\mathcal{D}}+1)$.
Then,

$\displaystyle\mathbb{E}_{\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\left[\sup_{f\in\mathcal{D},g\in\mathcal{G}}\left\\{\mathbb{E}f(X,g(X,\eta))-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},g(X_{i},\eta_{i}))\right\\}\right]$
$\displaystyle\quad\leq
CB_{2}\sqrt{\frac{(H_{\mathcal{G}}+H_{\mathcal{D}}+1)S_{\mathcal{G},\mathcal{D}}\log(S_{\mathcal{G},\mathcal{D}})[\log
n+\log B_{2}]}{n}},$

where $C>0$ is a universal constant.

###### Proof.

Let $h(X,\eta)=f(X,g(X,\eta))$ and
$\mathcal{F}=\mathcal{NN}(m+d,1,2W_{\mathcal{G}}+2W_{\mathcal{D}}+2d,H_{\mathcal{G}}+H_{\mathcal{D}}+1)+1$.
Then, $\sup_{h\in\mathcal{F}}\|h\|_{\infty}\leq B_{2}$. According to Remarks
13-14 in Nakada and Imaizumi, (2020), we can show by a similar argument to the
proof of Lemma 8 that

$\displaystyle\mathbb{E}_{\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\left[\sup_{f\in\mathcal{D},g\in\mathcal{G}}\left\\{\mathbb{E}f(X,g(X,\eta))-\frac{1}{n}\sum_{i=1}^{n}f(X_{i},g(X_{i},\eta_{i0}))\right\\}\right]$
$\displaystyle\leq$
$\displaystyle\mathbb{E}_{\\{(X_{i},Y_{i})\\}_{i=1}^{n}}\left[\sup_{h\in\mathcal{F}}\left\\{\mathbb{E}h(X,\eta)-\frac{1}{n}\sum_{i=1}^{n}h(X_{i},\eta_{i0})\right\\}\right]$
$\displaystyle\leq$ $\displaystyle
CB_{2}\sqrt{\frac{(H_{\mathcal{G}}+H_{\mathcal{D}}+1)S_{\mathcal{G},\mathcal{D}}\log(S_{\mathcal{G},\mathcal{D}})[\log
n+\log B_{2}]}{n}},$

where $C>0$ is a universal constant.  ∎

##  Appendix D Additional simulation results

In this section, we provide additional simulation results.  Other than Models
1 and 2 considered in the main text, three additional models are used to show
the performance of the proposed method.

Model 1. A nonlinear model with an additive error term:

$\displaystyle
Y=X_{1}^{2}+\exp\left(X_{2}+X_{3}/3\right)+\sin\left(X_{4}+X_{5}\right)+\varepsilon,\quad\varepsilon\sim
N(0,1).$

Model 2. A model with an additive error term whose variance depends on the
predictors:

$\displaystyle
Y=X_{1}^{2}+\exp\left(X_{2}+X_{3}/3\right)+X_{4}-X_{5}+\left(0.5+X_{2}^{2}/2+X_{5}^{2}/2\right)\times\varepsilon,\
\varepsilon\sim N(0,1)$

Model 6. A nonlinear model with a heavy tail additive error term:

$\displaystyle
Y=X_{1}^{2}+\exp\left(X_{2}+X_{3}/3\right)+\sin\left(X_{4}+X_{5}\right)+\varepsilon,\quad\varepsilon\sim
t(3).$

Model 7. A model with a multiplicative non-Gassisan error term:

$\displaystyle
Y=(5+X_{1}^{2}/3+X_{2}^{2}+X_{3}^{2}+X_{4}+X_{5})\times\exp(0.5\epsilon),$

where $\epsilon\sim 0.5N(-2,1)+0.5N(2,1).$

Model 8. A mixture of two normal distributions:

$\displaystyle
Y=I_{\\{U<0.5\\}}N(-X_{1},0.25^{2})+I_{\\{U>0.5\\}}N(X_{1},0.25^{2}),$

where $U\sim$Uniform$(0,1)$ and is independent of $X$.

We use the same evaluation criteria in Section 5.1.

###  D.1 Simulation results for different methods for Models 1-2, 6-8

We compare the performance of the proposed method to the nonparametric least
squares regression (NLS) and conditional WGAN (cWGAN). The numerical results
are summarized in Tables A.1 and A.2. It can be seen that the proposed method
has smaller MSEs for most cases, indicating that the proposed method can
improve the distribution matching to some extent.

Table A.1: The $L_{1}$ and $L_{2}$ errors, MSE of the estimated conditional
mean and the estimated standard deviation by different methods.

$X$ | Model | Method | $L_{1}$ | $L_{2}$ | Mean | SD
---|---|---|---|---|---|---
$N(0,I_{d})$ | M1 | NLS | 0.83(0.02) | 1.08(0.05) | - | -
cWGAN | 1.00(0.02) | 1.97(0.25) | 0.98(0.25) | 0.09(0.03)
WGR |  0.82(0.02) |  1.07(0.05) |  0.06(0.01) |  0.04(0.01)
M2 | NLS | 1.24(0.04) |  2.45(0.38) | - | -
cWGAN | 1.62(0.05) | 4.11(0.44) | 0.85(0.24) | 0.37(0.10)
WGR |  1.23(0.04) | 3.41(0.38) |  0.19(0.10) |  0.22(0.04)
M6 | NLS | 1.17(0.04) | 3.05(0.70) | - | -
cWGAN | 1.31(0.03) | 3.95(0.53) | 2.92(0.46) | 0.39(0.15)
WGR |  1.16(0.05) |  3.00(0.71) |  2.49(0.14) |  0.28(0.14)
M7 | NLS |  2.23(0.09) |  9.43(1.09) | - | -
cWGAN | 2.31(0.09) | 10.13(1.21) | 0.91(0.12) | 1.00(0.08)
WGR | 2.24(0.09) | 9.47(1.08) |  0.20(0.03) |  0.85(0.09)
M8 | NLS |  0.83(0.02) | 1.09(0.07) | - | -
cWGAN |  0.83(0.03) | 1.09(0.07) |  0.01(0.00) |  0.53(0.14)
WGR |  0.83(0.02) |  1.09(0.07) |  0.01(0.00) | 0.62(0.09)
$N(0,I_{100})$ | M1 | NLS | 1.17(0.03) | 2.52(0.15) | - | -
cWGAN | 1.17(0.04) | 2.65(0.41) | 1.67(0.39) | 0.81(0.02)
WGR |  1.15(0.04) |  2.40(0.24) |  1.64(0.25) |  0.16(0.04)
M2 | NLS | 1.63(0.05) | 5.37(0.47) | - | -
cWGAN | 1.64(0.06) | 5.34(0.46) | 2.20(0.22) | 1.11(0.08)
WGR |  1.61(0.06) |  5.27(0.45) |  2.06(0.24) |  0.23(0.03)
M6 | NLS |  1.58(0.03) |  5.09(0.62) | - | -
cWGAN | 1.71(0.08) | 6.01(0.80) | 5.16(0.74) | 1.19(0.45)
WGR | 1.62(0.06) | 5.39(0.66) |  4.20(0.53) |  0.47(0.19)
M7 | NLS | 2.51(0.09) | 11.52(1.17) | - | -
cWGAN | 2.51(0.09) | 12.02(1.49) | 2.88(0.29) | 0.96(0.25)
WGR |  2.45(0.11) |  11.37(1.30) |  2.22(0.17) |  0.73(0.06)
M8 | NLS | 0.85(0.02) | 1.13(0.07) | - | -
cWGAN | 0.89(0.03) | 1.29(0.10) |  0.00(0.00) |  0.06(0.00)
WGR |  0.83(0.02) |  1.09(0.07) | 0.19(0.05) |  0.06(0.30)

NOTE: The corresponding simulation standard errors are given in parentheses.
WGR is the propoed method.

Table A.2: MSE of the estimated conditional quantiles by different methods.

$X$ | Model | Method | $\tau=0.05$ | $\tau=0.25$ | $\tau=0.50$ | $\tau=0.75$ | $\tau=0.95$
---|---|---|---|---|---|---|---
$N(0,I_{d})$ | M1 | cWGAN | 1.22(0.23) | 1.04(0.25) | 0.99(0.26) | 1.03(0.24) | 1.34(0.21)
WGR |  0.29(0.07) |  0.10(0.01) |  0.09(0.02) |  0.10(0.03) |  0.23(0.06)
M2 | cWGAN | 1.86(0.21) | 0.94(0.26) | 0.85(0.26) | 1.00(0.21) | 2.59(0.52)
WGR |  0.77(0.09) |  0.31(0.06) |  0.19(0.04) |  0.27(0.05) |  0.81(0.15)
M6 | cWGAN | 3.57(0.73) | 3.29(0.58) | 2.95(0.50) | 3.01(0.46) | 6.14(1.44)
WGR |  3.39(0.48) |  2.70(0.24) |  2.47(0.12) |  2.50(0.21) |  2.76(0.50)
M7 | cWGAN |  0.17(0.03) |  0.28(0.05) | 0.60(0.09) | 1.99(0.21) | 8.40(0.79)
WGR | 2.01(0.47) | 0.41(0.10) |  0.39(0.06) |  0.74(0.11) |  4.04(0.40)
M8 | cWGAN |  0.03(0.02) | 0.04(0.03) |  0.36(0.04) | 0.04(0.03) | 0.03(0.02)
WGR |  0.03(0.01) |  0.02(0.01) | 0.40(0.04) |  0.02(0.01) |  0.03(0.02)
$N(0,I_{100})$ | M1 | cWGAN | 3.18(0.71) | 1.83(0.22) | 1.75(0.16) | 1.89(0.21) | 3.61(0.39)
WGR |  1.84(0.23) |  1.69(0.20) |  1.66(0.16) |  1.88(0.13) |  2.41(0.14)
M2 | cWGAN | 4.99(0.51) | 2.63(0.24) | 2.21(0.22) | 2.79(0.36) | 5.41(0.65)
WGR |  3.42(1.07) |  2.26(0.28) |  2.19(0.22) |  2.57(0.27) |  3.49(0.47)
M6 | cWGAN | 7.79(1.79) | 5.95(1.09) | 5.04(0.80) | 6.15(0.88) | 12.86(2.92)
WGR |  5.47(0.76) |  4.89(0.67) |  4.72(0.79) |  4.91(0.78) |  4.65(0.71)
M7 | cWGAN |  1.43(0.24) | 1.92(0.11) | 2.65(0.17) | 4.21(0.49) | 10.38(1.89)
WGR | 1.88(0.186) |  1.82(0.13) |  2.17(0.17) |  3.11(0.23) |  7.15(0.59)
M8 | cWGAN | 1.69(0.10) | 1.06(0.08) |  0.34(0.04) | 0.99(0.08) | 1.59(0.10)
WGR |  0.64(0.13) |  0.64(0.13) | 0.56(0.10) |  0.66(0.11) |  0.62(0.10)

NOTE: The corresponding simulation standard errors are given in parentheses.
WGR is the propoed method.

###  D.2 Simulation results with varying noise dimension

We check the performance of the proposed method with varying dimension of the
noise $\eta$. Other than $m=3$ considered in Section 5.1, we also try $m=10$
and $m=25$. The results are given in Tables A.3 and A.4. For larger $m$, there
tends to be a smaller MSE of the estimated conditional quantiles, but it takes
longer time to train. There is trade-off between performance and training
time. We also observe that the performance of the proposed method is somewhat
robust to the noise dimension $m$, as long as $m$ takes reasonable values.

Table A.3: The $L_{1}$ and $L_{2}$ error, MSE of the estimated conditional
mean and the estimated standard deviation for different $m$.

$X$ | Model | $m$ | $L_{1}$ | $L_{2}$ | Mean | SD
---|---|---|---|---|---|---
$N(0,I_{d})$ | M1 | 3 | 0.83(0.02) | 1.07(0.05) | 0.06(0.01) | 0.04(0.01)
10 | 0.83(0.02) | 1.08(0.05) | 0.06(0.01) | 0.04(0.2)
25 | 0.83(0.03) | 1.08(0.08) | 0.06(0.01) | 0.07(0.02)
M2 | 3 | 1.24(0.04) | 3.42(0.38) | 0.18(0.13) | 0.25(0.08)
10 | 1.23(0.05) | 3.38(0.39) | 0.14(0.03) | 0.14(0.03)
25 | 1.21(0.04) | 3.21(0.28) | 0.16(0.04) | 0.20(0.05)
M6 | 3 | 1.16(0.05) | 3.00(0.72) | 2.49(0.14) | 0.28(0.14)
10 | 1.16(0.05) | 3.02(0.72) | 2.36(0.16) | 0.20(0.06)
25 | 1.18(0.07) | 3.10(0.72) | 2.37(0.17) | 0.18(0.06)
M7 | 3 | 2.24(0.09) | 9.47(1.08) | 0.20(0.03) | 0.85(0.09)
10 | 2.24(0.09) | 9.43(1.09) | 0.20(0.02) | 0.80(0.15)
25 | 2.21(0.10) | 9.21(1.22) | 0.22(0.04) | 0.77(0.07)
M8 | 3 | 0.83(0.02) | 1.09(0.07) | 0.01(0.01) | 0.62(0.09)
10 | 0.83(0.02) | 1.09(0.07) | 0.00(0.00) | 0.63(0.10)
25 | 0.83(0.02) | 1.08(0.07) | 0.00(0.00) | 0.56(0.13)
$N(0,I_{100})$ | M1 | 3 | 1.15(0.04) | 2.40(0.24) | 1.64(0.25) | 0.16(0.04)
10 | 1.25(0.03) | 2.83(0.21) | 1.86(0.22) | 0.17(0.05)
25 | 1.18(0.04) | 2.47(0.17) | 1.44(0.13) | 0.16(0.05)
M2 | 3 | 1.62(0.05) | 5.36(0.46) | 2.11(0.22) | 0.36(0.12)
10 | 1.63(0.08) | 5.40(0.57) | 2.20(0.23) | 0.38(0.12)
25 | 1.62(0.06) | 5.31(0.43) | 2.17(0.19) | 0.45(0.19)
M6 | 3 | 1.62(0.06) | 5.39(0.66) | 4.20(0.53) | 0.47(0.19)
10 | 1.63(0.05) | 5.34(0.56) | 4.47(0.62) | 0.65(0.18)
25 | 1.58(0.03) | 4.92(0.33) | 4.20(0.39) | 0.50(0.13)
M7 | 3 | 2.45(0.11) | 11.37(1.30) | 2.22(0.17) | 0.73(0.06)
10 | 2.44(0.10) | 11.14(1.08) | 2.15(0.24) | 0.80(0.15)
25 | 2.46(0.09) | 11.26(1.20) | 2.10(0.25) | 0.73(0.06)
M8 | 3 | 0.83(0.02) | 1.09(0.07) | 0.19(0.05) | 0.30(0.16)
10 | 0.90(0.02) | 1.34(0.06) | 0.24(0.08) | 0.34(0.12)
25 | 0.89(0.04) | 1.27(0.11) | 0.15(0.08) | 0.40(0.13)

NOTE: The corresponding standard errors are given in parentheses.

Table A.4: Mean squared prediction error (MSE) of the estimated conditional
quantile for different values of $m$.

$X$ | Model | $m$ | $\tau=0.05$ | $\tau=0.25$ | $\tau=0.50$ | $\tau=0.75$ | $\tau=0.95$
---|---|---|---|---|---|---|---
$N(0,I_{d})$ | M1 | 3 | 0.29(0.07) | 0.10(0.01) | 0.09(0.02) | 0.10(0.03) | 0.23(0.06)
10 | 0.18(0.04) | 0.09(0.02) | 0.08(0.01) | 0.09(0.01) | 0.16(0.05)
25 | 0.26(0.05) | 0.10(0.02) | 0.07(0.02) | 0.11(0.10) | 0.29(0.10)
M2 | 3 | 0.77(0.09) | 0.31(0.06) | 0.19(0.04) | 0.27(0.05) | 0.81(0.15)
10 | 0.49(0.11) | 0.21(0.05) | 0.15(0.03) | 0.20(0.04) | 0.54(0.12)
25 | 0.62(0.18) | 0.24(0.07) | 0.16(0.04) | 0.23(0.04) | 0.76(0.16)
M6 | 3 | 3.39(0.48) | 2.70(0.24) | 2.47(0.12) | 2.50(0.21) | 2.76(0.50)
10 | 2.74(0.26) | 2.39(0.19) | 2.39(0.20) | 2.59(0.21) | 3.07(0.45)
25 | 2.97(0.30) | 2.46(0.20) | 2.37(0.20) | 2.55(0.21) | 3.02(0.35)
M7 | 3 | 2.01(0.47) | 0.40(0.10) | 0.39(0.06) | 0.74(0.11) | 4.04(0.40)
10 | 1.51(0.66) | 0.56(0.05) | 0.31(0.05) | 0.63(0.19) | 4.73(0.68)
25 | 1.07(0.14) | 0.50(0.11) | 0.32(0.09) | 0.57(0.10) | 3.88(0.51)
M8 | 3 | 0.03(0.01) | 0.02(0.01) | 0.40(0.04) | 0.02(0.01) | 0.03(0.02)
10 | 0.02(0.01) | 0.02(0.01) | 0.41(0.05) | 0.02(0.01) | 0.02(0.01)
25 | 0.03(0.02) | 0.02(0.02) | 0.41(0.04) | 0.02(0.02) | 0.03(0.01)
$N(0,I_{100})$ | M1 | 3 | 1.84(0.23) | 1.69(0.20) | 1.75(0.16) | 1.88(0.13) | 2.41(0.14)
10 | 2.05(0.30) | 1.91(0.25) | 1.92(0.22) | 2.02(0.23) | 2.52(0.36)
25 | 1.69(0.17) | 1.42(0.15) | 1.47(0.13) | 1.65(0.12) | 2.02(0.22)
M2 | 3 | 3.42(1.07) | 2.26(0.28) | 2.19(0.22) | 2.57(0.27) | 3.49(0.47)
10 | 3.07(0.67) | 2.29(0.21) | 2.33(0.26) | 2.64(0.37) | 3.55(0.44)
25 | 3.58(1.15) | 2.25(0.17) | 2.34(0.17) | 2.79(0.12) | 3.54(0.26)
M6 | 3 | 4.47(0.76) | 4.89(0.67) | 4.72(0.79) | 4.91(0.78) | 5.65(0.71)
10 | 5.91(0.55) | 4.68(0.55) | 4.55(0.61) | 4.94(0.68) | 6.17(0.74)
25 | 4.28(0.49) | 4.08(0.36) | 4.36(0.38) | 4.71(0.57) | 5.35(0.76)
M7 | 3 | 1.88(0.19) | 1.82(0.13) | 2.17(0.17) | 3.11(0.23) | 7.15(0.59)
10 | 1.91(0.29) | 1.87(0.17) | 2.16(0.21) | 3.06(0.30) | 7.22(0.91)
25 | 1.89(0.23) | 1.79(0.18) | 2.10(0.25) | 3.00(0.37) | 6.77(0.69)
M8 | 3 | 0.64(0.13) | 0.64(0.13) | 0.56(0.10) | 0.66(0.11) | 0.62(0.10)
10 | 0.67(0.12) | 0.62(0.10) | 0.63(0.08) | 0.64(0.12) | 0.59(0.08)
25 | 0.69(0.17) | 0.60(0.09) | 0.52(0.10) | 0.65(0.06) | 0.67(0.09)

NOTE: The corresponding simulation standard errors are given in parentheses.

###  D.3 Simulation results for different values of $J$

We conduct simulation studies to check the performance of the proposed method
for different $J$, the sample size of the generated noise vector $\eta$. The
generator and discriminator networks have two fully-connected hidden layers.
The LeakyReLU function is used as the active function. The noise vector $\eta$
is generated from $N(\bm{0},\bm{I}_{3})$. We take $J=10,50,200,500$, and
report the results in Tables A.5 and A.6. It can be seen that the proposed
method works comparably for different $J$.

Table A.5: The $L_{1}$ and $L_{2}$ error, MSE of the estimated conditional
mean and the estimated standard deviation for different $J$.

$X$ | Model | $J$ | $L_{1}$ | $L_{2}$ | Mean | SD
---|---|---|---|---|---|---
$N(0,I_{d})$ | M1 | 10 | 0.83(0.02) | 1.09(0.04) | 0.07(0.01) | 0.10(0.03)
50 | 0.82(0.02) | 1.07(0.04) | 0.06(0.02) | 0.05(0.02)
200 | 0.83(0.02) | 1.07(0.05) | 0.06(0.01) | 0.04(0.01)
500 | 0.82(0.02) | 1.06(0.06) | 0.05(0.01) | 0.05(0.02)
M2 | 10 | 1.24(0.04) | 3.42(0.38) | 0.18(0.13) | 0.25(0.08)
50 | 1.23(0.05) | 3.42(0.40) | 0.19(0.14) | 0.26(0.09)
200 | 1.23(0.04) | 3.41(0.38) | 0.19(0.10) | 0.22(0.04)
500 | 1.22(0.04) | 3.40(0.38) | 0.16(0.13) | 0.22(0.03)
M6 | 10 | 1.16(0.03) | 2.86(0.30) | 2.60(0.22) | 0.35(0.13)
50 | 1.16(0.04) | 3.02(0.66) | 2.51(0.15) | 0.24(0.14)
200 | 1.16(0.05) | 3.00(0.72) | 2.49(0.14) | 0.28(0.14)
500 | 1.16(0.05) | 3.00(0.71) | 2.49(0.14) | 0.28(0.14)
M7 | 10 | 2.23(0.09) | 9.44(1.08) | 0.21(0.03) | 0.95(0.12)
50 | 2.24(0.09) | 9.42(1.10) | 0.21(0.04) | 0.77(0.15)
200 | 2.24(0.09) | 9.47(1.08) | 0.20(0.03) | 0.85(0.09)
500 | 2.23(0.09) | 9.40(1.07) | 0.18(0.02) | 0.86(0.09)
M8 | 10 | 0.83(0.03) | 1.09(0.07) | 0.00(0.00) | 0.59(0.15)
50 | 0.83(0.02) | 1.08(0.07) | 0.00(0.00) | 0.62(0.14)
200 | 0.83(0.02) | 1.09(0.07) | 0.01(0.01) | 0.62(0.09)
500 | 0.83(0.02) | 1.09(0.07) | 0.00(0.00) | 0.60(0.11)
$N(0,I_{100})$ | M1 | 10 | 1.17(0.05) | 2.40(0.21) | 1.71(0.20) | 0.20(0.13)
50 | 1.15(0.03) | 2.50(022) | 1.65(0.22) | 0.14(0.04)
200 | 1.15(0.04) | 2.40(0.24) | 1.64(0.25) | 0.16(0.04)
500 | 1.20(0.03) | 2.45(0.20) | 1.74(0.18) | 0.13(0.03)
M2 | 10 | 1.62(0.05) | 5.36(0.46) | 2.11(0.22) | 0.36(0.12)
50 | 1.59(0.04) | 5.23(0.54) | 1.98(0.27) | 0.37(0.10)
200 | 1.61(0.06) | 5.27(0.45) | 2.06(0.24) | 0.30(0.03)
500 | 1.62(0.05) | 5.39(0.48) | 2.08(0.24) | 0.30(0.05)
M6 | 10 | 1.58(0.05) | 5.14(0.54) | 4.32(0.51) | 0.43(0.17)
50 | 1.57(0.05) | 5.06(0.52) | 4.18(0.59) | 0.53(0.17)
200 | 1.62(0.06) | 5.39(0.66) | 4.20(0.53) | 0.47(0.19)
500 | 1.57(0.05) | 5.05(0.56) | 4.14(0.50) | 0.44(0.13)