Learning Theory of Distribution Regression with Neural Networks
In this paper, we aim at establishing an approximation theory and a learning
theory of distribution regression via a fully connected neural network (FNN).
In contrast to the classical regression methods, the input variables of
distribution regression are probability measures. Then we often need to perform
a second-stage sampling process to approximate the actual information of the
distribution. On the other hand, the classical neural network structure
requires the input variable to be a vector. When the input samples are
probability distributions, the traditional deep neural network method cannot be
directly used and the difficulty arises for distribution regression. A
well-defined neural network structure for distribution inputs is intensively
desirable. There is no mathematical model and theoretical analysis on neural
network realization of distribution regression. To overcome technical
difficulties and address this issue, we establish a novel fully connected
neural network framework to realize an approximation theory of functionals
defined on the space of Borel probability measures. Furthermore, based on the
established functional approximation results, in the hypothesis space induced
by the novel FNN structure with distribution inputs, almost optimal learning
rates for the proposed distribution regression model up to logarithmic terms
are derived via a novel two-stage error decomposition technique.
Zhongjie Shi, Zhan Yu, Ding-Xuan Zhou
# Learning Theory of Distribution Regression with Neural Networks

###### Abstract

In this paper, we aim at establishing an approximation theory and a learning theory of distribution regression via a fully connected neural network (FNN). In contrast to the classical regression methods, the input variables of distribution regression are probability measures. Then we often need to perform a second-stage sampling process to approximate the actual information of the distribution. On the other hand, the classical neural network structure requires the input variable to be a vector. When the input samples are probability distributions, the traditional deep neural network method cannot be directly used and the difficulty arises for distribution regression. A well-defined neural network structure for distribution inputs is intensively desirable. There is no mathematical model and theoretical analysis on neural network realization of distribution regression. To overcome technical difficulties and address this issue, we establish a novel fully connected neural network framework to realize an approximation theory of functionals defined on the space of Borel probability measures. Furthermore, based on the established functional approximation results, in the hypothesis space induced by the novel FNN structure with distribution inputs, almost optimal learning rates for the proposed distribution regression model up to logarithmic terms are derived via a novel two-stage error decomposition technique.

_Keywords_: Learning theory, distribution regression, neural networks, approximation rates, learning rates.

## 1 Introduction

Recent years have witnessed a vast number of applications of deep learning in various fields of science and engineering. Deep learning based on deep neural networks has become a
powerful tool to provide various models and structures for many complex learning tasks, taking advantages of the huge improvement of computing power in the era of big data [15]. The classical fully connected neural network for learning functions of input variable vector \(x=(x_{i})_{i=1}^{d}\in\mathbb{R}^{d}\) with \(J\) layers of neurons \(\{H^{(k)}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d_{k}}\}_{k=1}^{J}\) with width \(\{d_{k}\}_{k=1}^{J}\) is defined iteratively by

\[H^{(k)}(x)=\sigma\left(F^{(k)}H^{(k-1)}(x)-b^{(k)}\right),\quad k=1,2,...,J, \tag{1.1}\]

where \(\sigma:\mathbb{R}\rightarrow\mathbb{R}\) is an activation function acting on vectors component-wise, \(F^{(k)}\) is a \(d_{k}\times d_{k-1}\) matrix, \(b^{(k)}\) is a bias vector, \(H^{(0)}(x)=x\) with width \(d_{0}=d\). Starting from 1980s, approximation theory of neural networks has been well developed in lots of studies, such as [8, 18, 1, 24, 4, 36, 2, 41, 21, 9].

The classical theory of deep neural networks circling around the generalization with FNNs has been walking towards mature within the last few past decades, and the regression problem by using neural networks is always a popular topic. There have been several studies on learning rates of classical least square regression with deep neural networks (DNNs), and satisfactory learning rates have been derived [5, 28]. However, in the existing work related to the regression analysis with DNNs, the input data are still limited to the Euclidean vectors. A natural question appears when the data we need to handle are probability measures and the neural network structure defined in aforementioned works is not suitable for handling the distribution data directly. That is, how can we utilize deep neural networks to handle such distribution data? Solving this problem is desirable, since in recent years there are a lot of applications dealing with distribution data [27, 33]. One of the most popular objects in statistical learning theory is the distribution regression. However, it is still difficult to perform distribution regression in FNNs since an initial model to realize distribution regression with neural networks is still lacking. One of the main difficulties to realize distribution regression by utilizing neural networks is to establish an appropriate network structure induced by FNNs that is fit for distribution inputs. To overcome the limitation of existing neural network techniques and answer the above question, in this paper we provide a novel FNN structure to realize distribution regression and establish a meaningful regression model, finally we derive the satisfactory learning rates.

Recently, there has been some progress in distribution regression which aims at regressing from probability measures to real-valued responses [26, 32, 33, 13, 37]. In such problems, in contrast to classical regression problems with vector inputs, we need to handle the inputs consisting of probability measures. Hence the sampling process often involves two stages. In the first stage, the distribution samples are selected from some meta distribution. To further get access to the actual information of the probability measure samples, we need to draw samples from the distribution samples, that is the so called second stage process. In the literature of distribution regression, most of the existing works and their related regression models are based on a kernel mean embedding technique to transform probability measure samples to the space of their mean embeddings which is compact. In these works, they choose the hypothesis space as the reproducing kernel Hilbert space induced by some Mercer kernel \(K\), and use some regularization and integral operator techniques (e.g.[29, 30]), some nice learning rates are derived. However, to realize the learning
theory of distribution regression with FNNs, the aforementioned traditional kernel-based techniques break down. One of the main reasons is that the architecture of the neural network is essentially different from the kernel based hypothesis space. Furthermore, as we need to handle the distribution data instead of the classical Euclidean vectors in neural networks, the classical neural network structure as in (1.1) cannot be utilized directly. Hence, some new neural network structures must be established for handling distributions. To this end, we will first define a novel neural network structure with the variable of probability measures that will be given in Definition 1 below. In the hypothesis space induced by the novel network structure, we establish a distribution regression framework with probability samples drawn from a class of Wasserstein space \((\mathcal{P}(\Omega),W_{p})\), where \(\mathcal{P}(\Omega)\) denotes a collection of Borel probability measures defined on a compact space \(\Omega\subset\mathbb{R}^{d}\), \(\mathcal{P}(\Omega)\) is naturally a subset of the space \(\mathcal{B}(\Omega)\) consisting of all Borel measures on \(\Omega\), and \(W_{p}\) denotes the well-known Wasserstein metric that will be defined later.

We first introduce our two-stage distribution regression model where the inputs are distribution samples on the Wasserstein space \((\mathcal{P}(\Omega),W_{p})\). In our distribution regression model, we assume that the distribution samples cannot be observed directly, the data are generated through a two-stage sampling process. In the first stage of sampling, the dataset \(D=\{(\mu_{i},y_{i})\}_{i=1}^{m}\) is i.i.d. sampled from a meta Borel distribution \(\rho\) on \(\mathcal{U}\times\mathcal{Y}\), where \(\mathcal{U}=\mathcal{P}(\Omega)\), and \(\mathcal{Y}=\mathbb{R}\) is the output space. In the second stage of sampling, the dataset is \(\hat{D}=\{(\{x_{i,j}\}_{j=1}^{n_{i}},y_{i})\}_{i=1}^{m}\), where \(\{x_{i,j}\in\Omega\}_{j=1}^{n_{i}}\) are i.i.d. sampled from the probability distribution \(\mu_{i}\), which is a first-stage sample. If we use the empirical terminology \(\hat{\mu}_{i}^{n_{i}}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\delta_{x_{i,j}}\), then we can also denote the second-stage sample set by

\[\hat{D}=\{(\hat{\mu}_{i}^{n_{i}},y_{i})\}_{j=1}^{m}.\]

The classical fully connected neural network takes a vector as input, which is not suitable for distribution regression. Since we may have samples of different size \(n_{i}\) for different sample distribution \(\mu_{i}\) in the second stage of sampling, if we just vectorize all the \(n_{i}\) samples in the second stage to a vector as the input of a FNN, when the second-stage sample size \(n_{i}\) is different for different first-stage sample \((\mu_{i},y_{i})\), the input dimension of the vector for the FNN is also different, thus making the traditional FNN with vector inputs fails in the distribution regression problem. Therefore, we propose a novel general FNN structure which takes a distribution rather than a vector as the input for the learning of distribution regression. In practical settings, this network structure is able to utilize the empirical distribution

\[\hat{\mu}_{i}^{n_{i}}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\delta_{x_{i,j}} \tag{1.2}\]

of the second-stage sample rather than the vector as the input. For a Borel measure \(\mu\) and an integrable function vector \(h:\Omega\rightarrow\mathbb{R}^{D}\), if we use the vector integral notation

\[\int_{\Omega}h(x)d\mu(x)=\left[\begin{array}{c}\int_{\Omega}(h(x))_{1}d\mu( x)\\ \int_{\Omega}(h(x))_{2}d\mu(x)\\ \vdots\\ \int_{\Omega}(h(x))_{D}d\mu(x)\end{array}\right],\]
then the FNN scheme for distribution regression is given in the following definition.

**Definition 1** (FNN for distribution regression).: _Let \(\mu\in\mathcal{P}(\Omega)\) be an input distribution, the FNN for distribution regression \(\{h^{(j)}:\mathcal{P}(\Omega)\rightarrow\mathbb{R}^{d_{j}}\}_{j=J_{1}}^{J}\) of type \((J_{1},J)\), with depth \(J\in\mathbb{N}\) and realizing level \(J_{1}\in\{0,1,...,J\}\), and width \(\{d_{j}\}_{j=1}^{J}\) is defined iteratively by:_

\[h^{(j)}(\mu)=\left\{\begin{array}{ll}\int_{\Omega}\sigma\left(F^{(J_{1})} \sigma\left(\cdots\sigma\left(F^{(1)}x-b^{(1)}\right)\cdots\right)-b^{(J_{1}) }\right)d\mu(x),&\mbox{if }j=J_{1},\\ \sigma\left(F^{(j)}h^{(j-1)}(\mu)-b^{(j)}\right),&\mbox{if }j=J_{1}+1,\ldots,J. \end{array}\right. \tag{1.3}\]

_where \(\sigma\) is the ReLU activation function given by \(\sigma(u)=\max\{u,0\}\), \(F^{(j)}\in\mathbb{R}^{d_{j}\times d_{j-1}}\) is a full connection matrix with \(d_{0}=d\), and \(b^{(j)}\in\mathbb{R}^{d_{j}}\) is a bias vector._

An output of the above FNN is just a linear combination of the last layer

\[c\cdot h^{(J)}(\mu),\]

where \(c\in\mathbb{R}^{d_{J}}\) is a coefficient vector. In practical settings, an input probability measure always appears as an empirical distribution. For example, in practice, the actual input always has the form \(\hat{\mu}_{i}^{n_{i}}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\delta_{x_{i,j}}\), \(i=1,2,...,m\), where \(x_{i,j}\in\Omega\subset\mathbb{R}^{d}\). Then the explicit form of \(h^{(J)}(\hat{\mu}_{i}^{n_{i}})\) can be represented layer by layer from

\[h^{(j)}(\hat{\mu}_{i}^{n_{i}})=\left\{\begin{array}{ll}\frac{1}{n_{i}}\sum_{ j=1}^{n_{i}}\sigma\left(F^{(J_{1})}\sigma\left(\cdots\sigma\left(F^{(1)}x_{i,j}-b^{( 1)}\right)\cdots\right)-b^{(J_{1})}\right),&\mbox{if }j=J_{1},\\ \sigma\left(F^{(j)}h^{(j-1)}(\hat{\mu}_{i}^{n_{i}})-b^{(j)}\right),&\mbox{if }j=J_{1}+1, \ldots,J.\end{array}\right.\] (1.3 \[*\] )

then the FNN structure becomes more flexible to handle such empirical variables.

If we use \(\mathcal{H}_{\mathrm{FNN}}\) to denote some hypothesis space induced by the above FNN structure for probability measures in Definition 1 (the rigorous definition will be given in Section 2.3), then a distribution regression scheme can be described as

\[f_{\hat{D},\mathcal{H}_{\mathrm{FNN}}}=\arg\min_{f\in\mathcal{H}_{\mathrm{FNN}}} \frac{1}{m}\sum_{i=1}^{m}\left(f(\hat{\mu}_{i}^{n_{i}})-y_{i}\right)^{2}. \tag{1.4}\]

In this paper, we first derive rates of approximating by FNNs two classes of composite functionals with distribution variables. The functional induced by ridge functions which we call a ridge functional will be considered first. Then the functional induced by a composite function associated with a polynomial \(Q\) will be considered. We will construct the approximation for the former functional with a FNN of \(J=2\), \(J_{1}=1\) and then for the latter one with a FNN of \(J=3\), \(J_{1}=2\) that are both involved in the FNN structure given in Definition 1. Then with the help of the functional approximation results and a covering number bound developed for the hypothesis space induced by the FNN structure defined above, we derive learning rates for the FNN distribution regression scheme. The covering number argument is in fact based on a crucial fact of the compactness of the hypothesis space, which is also the foundation of the learning theory framework. In this paper, we have rigorously shown that the new-defined hypothesis space with distribution inputs from the Wasserstein space is compact.
In fact, for a general distribution regression scheme without kernel regularization, developing appropriate theoretical results on the learning ability in some hypothesis space is still an open question itself. On the way of investigating the new distribution regression scheme and deriving learning rates, one novel two-stage error decomposition technique is proposed in this paper. To the best of our knowledge, such a two-stage error decomposition technique appears for the first time in the literature of the learning theory of distribution regression. It overcomes the difficulty of deriving learning rates for the non-kernel based distribution regression scheme for which the existing literature has not considered yet. Almost optimal learning rates for distribution regression up to logarithmic terms are derived by utilizing neural networks. Our method of deriving almost optimal learning rates also highlights the importance of the FNN structure we construct in the two-stage sampling distribution regression scheme.

**Notations**: we fix some notations and terminologies related to the \(\|\cdot\|_{\infty}\) norm of function or functional in different settings of this paper. In additional to the notations of \(\|\cdot\|_{\infty}\) for common scale vectors or matrices, without specification, for a continuous function (or say a "functional") \(f\) defined on the compact space \((\mathcal{P}(\Omega),W_{p})\), we use the notation \(\|h\|_{\infty}=\sup_{\mu\in\mathcal{P}(\Omega)}|h(\mu)|\). For a vector of functions \(h:\mathcal{P}(\Omega)\to\mathbb{R}^{D}\) defined on \(\mathcal{P}(\Omega)\), we use \(\|h\|_{\infty}=\max_{1\leq i\leq D}\sup_{\mu\in\mathcal{P}(\Omega)}|(h(\mu)) _{i}|\). For a vector of functions \(h:\Omega\to\mathbb{R}^{D}\), we also use \(\|h\|_{\infty}=\max_{1\leq i\leq D}\sup_{x\in\Omega}|(h(x))_{i}|\).

## 2 Main Results

We show the capacity of our proposed FNN for distribution regression by deriving the approximation rate of a class of composite nonlinear functional \(f\circ L_{G}\) with the variable of probability measures and a univariate function \(f:\mathbb{R}\to\mathbb{R}\)

\[f\circ L_{G}:\mu\mapsto f\left(L_{G}(\mu)\right)=f\left(\int_{\Omega}G(x)d\mu( x)\right). \tag{2.1}\]

The input \(\mu\) is a probability distribution defined on a compact set \(\Omega\subset\mathbb{B}\), where \(\mathbb{B}:=\{x\in\mathbb{R}^{d}:\|x\|_{2}\leq 1\}\) is the unit ball of the Euclidean space \(\mathbb{R}^{d}\), \(L_{G}\) is the inner functional of the composite nonlinear functional \(f\circ L_{G}\) in the form of

\[L_{G}(\mu)=\int_{\Omega}G(x)d\mu(x)\]

defined on the space \(\mathcal{P}(\Omega)\) of all Borel probability measures on \(\Omega\), and \(G\) is a function that induces the functional \(L_{G}\). In our setting, since \(\Omega\) is compact, the Riesz-Markov-Kakutani representation theorem asserts that the space \(\mathcal{B}(\Omega)\) of all finite Borel measures on \(\Omega\) is the dual of the space \(C(\Omega)\) of continuous functions on \(\Omega\). So, \(C(\Omega)\) is isometrically embedded into the closed space of the dual space of \(\mathcal{B}(\Omega)\). Hence, by considering a continuous function \(G\) on \(\Omega\) as the inducing function of the functional \(L_{G}\), the formulation of \(L_{G}(\mu)\) is naturally meaningful as a continuous functional defined on \(\mathcal{P}(\Omega)\subset\mathcal{B}(\Omega)\). Furthermore, since \(\Omega\) is a compact metric space equipped with the standard Euclidean metric, the space \(\mathcal{P}(\Omega)\) is also a compact metric space under the Wasserstein metric, which ensures that the
learning theory architecture of the proposed neural networks for distribution regression is well-defined. In the following, we will consider two classes of functionals, one is induced by the ridge function in the form of \(G(x)=g(\xi\cdot x)\) which is extremely popular in the realm of applications related to neural networks. The other one is induced by a more general composite function in the form of \(G(x)=g(Q(x))\) with a \(q\)-degree polynomial \(Q\).

### Rates of approximating ridge functionals

We first consider the case where \(G(x)=g(\xi\cdot x)\) is actually a ridge function, with a feature vector \(\xi\in\mathbb{R}^{d}\). We assume \(B_{\xi}=\|\xi\cdot x\|_{C(\Omega)}\) and \(g\in C^{0,1}[-B_{\xi},B_{\xi}]\), the space of Lipschitz functions on \([-B_{\xi},B_{\xi}]\), with semi-norm \(|g|_{C^{0,1}}:=\sup\limits_{x_{1}\neq x_{2}\in[-B_{\xi},B_{\xi}]}\frac{|g(x_{1 })-g(x_{2})|}{\|x_{1}-x_{2}\|_{2}}\), and \(B_{g}=\|g\|_{C([-B_{\xi},B_{\xi}])}\), and \(f\in C^{0,\beta}[-B_{G},B_{G}]\), the space of Lipschitz-\(\beta\) functions on \([-B_{G},B_{G}]\) with \(0<\beta\leq 1\), with semi-norm \(|f|_{C^{0,\beta}}:=\sup\limits_{x_{1}\neq x_{2}\in[-B_{G},B_{G}]}\frac{|f(x_{ 1})-f(x_{2})|}{\|x_{1}-x_{2}\|_{2}^{\beta}}\), in which \(B_{G}\) is defined as

\[B_{G}=B_{g}+2B_{\xi}|g|_{C^{0,1}} \tag{2.2}\]

and \(\|f\|_{\infty}=\|f\|_{C[-B_{G},B_{G}]}\). Then the functional to be approximated is of the form

\[f\circ L_{G}^{\xi}:\mu\mapsto f\left(L_{G}^{\xi}(\mu)\right)=f\left(\int_{ \Omega}g\left(\xi\cdot x\right)d\mu(x)\right),\ \mu\in\mathcal{P}(\Omega). \tag{2.3}\]

During the construction of a FNN to realize an approximation of the functional \(f\left(L_{G}(\mu)\right)\), we need to first approximate the inner functional \(L_{G}^{\xi}\) (we call it a ridge functional) induced by the feature vector \(\xi\in\mathbb{R}^{d}\) defined by

\[L_{G}^{\xi}:\mu\mapsto\int_{\Omega}g\left(\xi\cdot x\right)d\mu(x), \tag{2.4}\]

Our first result gives rates of approximating the composite functional (2.3) by a two-layer FNN for distribution regression defined in Definition 1. The _free parameters_ are defined as the implicit training parameters in our approximator construction, which is different from the total training parameters in the classical FNN structure due to the special structure of our construction.

**Theorem 1**.: _Let the functional \(f\circ L_{G}^{\xi}\) be in the ridge form (2.3), where \(\xi\in\mathbb{R}^{d}\) is a feature vector, \(g\in C^{0,1}[-B_{\xi},B_{\xi}]\), and \(f\in C^{0,\beta}[-B_{G},B_{G}]\) for some \(0<\beta\leq 1\). Then for any \(N\in\mathbb{N}\), there exists a FNN of type \((1,2)\) and width \(2N+3\) having the structure of Definition 1 with \(F^{(j)}\), \(b^{(j)}\), \(j=1,2\), explicitly constructed such that the following approximation rates hold,_

\[\inf_{\|c\|_{\infty}\leq\frac{4\|f\|_{\infty}N}{B_{G}}}\left\{\sup_{\mu\in \mathcal{P}(\Omega)}\left|c\cdot h^{(2)}(\mu)-f(L_{G}^{\xi}(\mu))\right|\right\} \leq\frac{2B_{G}^{\beta}\left|f\right|_{C^{0,\beta}}+\left(2B_{\xi}|g|_{C^{0, 1}}\right)^{\beta}\left|f\right|_{C^{0,\beta}}}{N^{\beta}}. \tag{2.5}\]

_The total number of free parameters in this network is \(\mathcal{N}=8N+d+12\)._
The architecture of the FNN of type \((1,2)\) we construct actually depends upon the formulation of the ridge functional \(f\left(L_{G}^{\xi}(\mu)\right)\). The functional \(L_{G}^{\xi}(\mu)\) includes the procedure of taking the integration w.r.t. the input distribution \(\mu\) on the ridge function \(g(\xi\cdot x)\), which can be approximated by a shallow net with one hidden layer. Therefore, the FNN we construct also takes the integration after the first hidden layer, corresponding to the number of hidden layer to be used for approximating ridge functions. Then we utilize another hidden layer outside the integration layer because the function \(f\) composed with the functional \(L_{G}^{\xi}\) can just be approximated by one hidden layer.

We end this subsection by presenting an interesting example. Recall the Laplace transform of a probability measure \(\mu\in\mathcal{P}(\Omega)\) is defined as

\[\mathcal{L}(\mu)(\xi)=\int_{\Omega}e^{-x\cdot\xi}d\mu(x),\quad\xi\in\mathbb{R }^{d}. \tag{2.6}\]

Then the FNN structure can be used to approximate a class of functionals \(\mathcal{L}_{\xi}\) induced by the Laplace transform of the elements in \(\mathcal{P}(\Omega)\) at \(\xi\):

\[\mathcal{L}_{\xi}:\mu\mapsto\mathcal{L}(\mu)(\xi). \tag{2.7}\]

The following corollary describes rates of this approximation.

**Corollary 1**.: _Let the functional \(\mathcal{L}_{\xi}\) be defined by (2.6) and (2.7), \(\xi\in\mathbb{R}^{d}\) is a feature vector. \(B_{\xi}=\|\xi\cdot x\|_{C(\Omega)}\). Then for any \(N\in\mathbb{N}\), there exists a FNN of type \((1,2)\) and width \(2N+3\) following the structure of Definition 1 with \(F^{(j)}\), \(b^{(j)}\), \(j=1,2\), explicitly constructed such that the following approximation rates hold,_

\[\inf_{\|c\|_{\infty}\leq 4N}\left\{\sup_{\mu\in\mathcal{P}(\Omega)}\Big{|}c \cdot h^{(2)}(\mu)-\mathcal{L}_{\xi}(\mu)\Big{|}\right\}\leq\frac{2e^{B_{\xi} }+6B_{\xi}e^{B_{\xi}}}{N}. \tag{2.8}\]

_The total number of free parameters in this network is \(\mathcal{N}=8N+d+12\). Furthermore, for any \(\xi\in\Omega\), the above bound can be further improved to_

\[\inf_{\|c\|_{\infty}\leq 4N}\left\{\sup_{\mu\in\mathcal{P}(\Omega)}\Big{|}c \cdot h^{(2)}(\mu)-\mathcal{L}_{\xi}(\mu)\Big{|}\right\}\leq\frac{8e}{N}. \tag{2.9}\]

### Rates of approximating composite functionals with polynomial features

Recall an important fact on the space \(\mathcal{P}_{q}^{h}(\mathbb{R}^{d})\) of homogeneous polynomials on \(\mathbb{R}^{d}\) of degree \(q\) from [20] and [22] that \(\mathcal{P}_{q}^{h}(\mathbb{R}^{d})\) has a basis \(\left\{(\xi_{k}\cdot x)^{q}\right\}_{k=1}^{n_{q}}\) for some vector set \(\left\{\xi_{k}\right\}_{k=1}^{n_{q}}\subset\mathbb{R}^{d}\setminus\{0\}\) and this vector set can even be chosen in such a way that the homogeneous polynomial set \(\left\{(\xi_{k}\cdot x)^{\ell}\right\}_{k=1}^{n_{q}}\) spans the space \(\mathcal{P}_{\ell}^{h}(\mathbb{R}^{d})\) for every \(\ell\in\{1,\ldots,q-1\}\), where \(n_{q}=\left(\begin{array}{c}d-1+q\\ q\end{array}\right)\) is the dimension of \(\mathcal{P}_{q}^{h}(\mathbb{R}^{d})\). Applying this fact to a polynomial \(Q\) of degree \(q\) yields the following lemma stated in [40].
**Lemma 1**.: _Let \(d\in\mathbb{N}\) and \(q\in\mathbb{N}\). Then there exists a set \(\left\{\xi_{k}\right\}_{k=1}^{n_{q}}\subset\left\{\xi\in\mathbb{R}^{d}:|\xi|=1\right\}\) of vectors with \(\ell_{2}\)-norm \(1\) such that for any \(Q\in\mathcal{P}_{q}(\mathbb{R}^{d})\) we can find a set of coefficients \(\left\{\gamma_{k,\ell}^{Q}:k=1,\ldots,n_{q},\ell=1,\ldots,q\right\}\subset \mathbb{R}\) such that_

\[Q(x)=Q(0)+\sum_{k=1}^{n_{q}}\sum_{\ell=1}^{q}\gamma_{k,\ell}^{Q}(\xi_{k}\cdot x )^{\ell},\qquad x\in\mathbb{R}^{d}. \tag{2.10}\]

Now we consider a general class of functional \(f\left(L_{G}^{Q}(\mu)\right)\) to be approximated, where \(G\) is of a composite form \(G(x)=g\left(Q(x)\right)\) and \(Q\) is a polynomial of degree \(q\) on \(\Omega\) with the constant \(\widehat{B}_{Q}=\|Q\|_{C(\Omega)}\). For the polynomial \(Q\), Lemma 1 indicates that there is a vector

\[\gamma_{Q}=\left[\gamma_{1,1}^{Q}\ \gamma_{1,2}^{Q}\ \ldots\ \gamma_{1,q}^{Q}\ \gamma_{2,1}^{Q}\ \gamma_{2,2}^{Q}\ \ldots\ \gamma_{2,q}^{Q}\ \ldots\ \gamma_{n_{q},1}^{Q}\ \gamma_{n_{q},2}^{Q}\ \ldots\ \gamma_{n_{q},q}^{Q}\right],\]

such that

\[Q(x)=Q(0)+\sum_{k=1}^{n_{q}}\sum_{\ell=1}^{q}\gamma_{k,\ell}^{Q}(\xi_{k}\cdot x )^{\ell},\qquad x\in\mathbb{R}^{d}. \tag{2.11}\]

We are now ready to define the following constant

\[B_{Q}=\widehat{B}_{Q}+2q\|\gamma_{Q}\|_{1}. \tag{2.12}\]

For the function \(g\), we assume \(g\in C^{0,1}[-B_{Q},B_{Q}]\), the space of Lipschitz-1 functions on \([-B_{Q},B_{Q}]\), with semi-norm \(|g|_{C^{0,1}}\), and \(B_{g}=\|g\|_{C[-B_{Q},B_{Q}]}\). We can then define the second constant

\[B_{G}=B_{g}+3B_{Q}|g|_{C^{0,1}} \tag{2.13}\]

for introducing the following Theorem 2.

For function \(f\), we assume \(f\in C^{0,\beta}[-B_{G},B_{G}]\), with semi-norm \(\left|f\right|_{C^{0,\beta}}\), and \(\|f\|_{\infty}=\|f\|_{C[-B_{G},B_{G}]}\). In this subsection, we aim at approximating the functional

\[f\circ L_{G}^{Q}:\ \mu\mapsto f\left(L_{G}^{Q}(\mu)\right)=f\left(\int_{\Omega}g \left(Q(x)\right)d\mu(x)\right),\ \mu\in\mathcal{P}(\Omega). \tag{2.14}\]

As in subsection 2.1, we need to first approximate the inner functional \(L_{G}^{Q}(\mu)\) on \(\mathcal{P}(\Omega)\) defined as

\[L_{G}^{Q}:\mu\mapsto\int_{\Omega}g\left(Q(x)\right)d\mu(x), \tag{2.15}\]

and derive the corresponding approximation rates of the composite nonlinear functional \(\mu\mapsto f\left(L_{G}^{Q}(\mu)\right)\) defined on \(\mathcal{P}(\Omega)\), with one more hidden layer to be used in the integral part of the FNN comparing with that of the ridge functional in subsection 2.1.

**Theorem 2**.: _Let the functional \(f\circ L_{G}^{Q}\) be in the composite form (2.14), where \(Q\) is a polynomial of degree \(q\) on \(\Omega\), \(g\in C^{0,1}[-B_{Q},B_{Q}]\), and \(f\in C^{0,\beta}[-B_{G},B_{G}]\) for some \(0<\beta\leq 1\). Then for any \(N\in\mathbb{N}\), there exists a FNN of type \((2,3)\) and widths \(d_{1}=n_{q}(2N+3)\)
\(d_{2}=d_{3}=2N+3\) following the structure of Definition 1 with \(F^{(j)}\), \(b^{(j)}\), \(j=1,2,3\) explicitly constructed, such that the following approximation rates hold,_

\[\inf_{\|c\|_{\infty}\leq\frac{4\|f\|_{\infty}N}{B_{G}}}\left\{\sup_{\mu\in \mathcal{P}(\Omega)}\left|c\cdot h^{(3)}(\mu)-f(L_{G}^{Q}(\mu))\right|\right\} \leq\frac{2B_{G}^{\beta}\left|f\right|_{C^{0,\beta}}+(3B_{Q}|g|_{C^{0,1}})^{ \beta}\left|f\right|_{C^{0,\beta}}}{N^{\beta}}. \tag{2.16}\]

_The total number of free parameters in this network is \(\mathcal{N}=(2q+10)N+(d+q)n_{q}+3q+15\)._

The reason that we construct a FNN of type \((2,3)\) also depends upon the formulation of the composite functional \(f\left(L_{G}^{Q}(\mu)\right)\). The number of hidden layers in the integration part is chosen to be two since we actually need a net with two hidden layers to approximate the composite function \(g\left(Q(x)\right)\). In fact, the position of the integration w.r.t. the input distribution \(\mu\) for the FNN we construct should be consistent with that in the true functional, otherwise the derivation of an approximation rate might be unsuccessful.

Specifically, when the polynomial \(Q\) in the functional-inducing function \(g(Q(x))\) takes the special form of \(Q(x)=\|x\|_{2}^{2}=x_{1}^{2}+x_{2}^{2}+\cdots+x_{d}^{2}\), then the composite functional-inducing function \(g\left(Q(x)\right)\) becomes the radial function \(g\left(\|x\|_{2}^{2}\right)\) and the functional has the form

\[f\circ L_{G}^{\|\cdot\|_{2}}:\mu\mapsto f\left(L_{G}^{\|\cdot\|_{2}}(\mu) \right)=f\left(\int_{\Omega}g\left(\|x\|_{2}^{2}\right)d\mu(x)\right). \tag{2.17}\]

Since we can write

\[Q(x)=\|x\|_{2}^{2}=\sum_{k=1}^{d}\left(e_{k}\cdot x\right)^{2}\]

in terms of the standard basis \(\{e_{k}\}_{k=1}^{d}\) of \(\mathbb{R}^{d}\), also \(\widehat{B}_{Q}=\|Q\|_{C(\Omega)}\leq 1\), we have \(B_{Q}=\widehat{B}_{Q}+2q\|\gamma_{Q}\|_{1}\leq 1+4d\).

We are then able to derive the following rates of approximating \(f\circ L_{G}^{\|\cdot\|_{2}}\) by utilizing Theorem 2 directly.

**Corollary 2**.: _Let the functional \(f\circ L_{G}^{\|\cdot\|_{2}}\) be in the composite form (2.17), where \(g\in C^{0,1}[-B_{Q},B_{Q}]\), and \(f\in C^{0,\beta}[-B_{G},B_{G}]\), for some \(0<\beta\leq 1\). Then for any \(N\in\mathbb{N}\), there exists a FNN of type \((2,3)\) and widths \(d_{1}=(2N+3)d\), \(d_{2}=d_{3}=2N+3\) following the structure of Definition 1, such that the following approximation rates hold,_

\[\inf_{\|c\|_{\infty}\leq\frac{4\|f\|_{\infty}N}{B_{G}}}\left\{\sup_{\mu\in \mathcal{P}(\Omega)}\left|c\cdot h^{(3)}(\mu)-f\left(L_{G}^{\|\cdot\|_{2}}( \mu)\right)\right|\right\}\leq\frac{2B_{G}^{\beta}\left|f\right|_{C^{0,\beta} }+(12d+3)^{\beta}|g|_{C^{0,1}}^{\beta}\left|f\right|_{C^{0,\beta}}}{N^{\beta}}. \tag{2.18}\]

_The total number of free parameters in this network is \(\mathcal{N}=12N+d^{2}+d+18\)._

### Distribution regression with FNN

In this subsection, we conduct generalization analysis of the commonly used empirical risk minimization (ERM) algorithm for distribution regression based on FNNs, which measures the generalization ability of the truncated empirical target function of the ERM algorithm.
We start with a formulation of the distribution regression model in this paper. First we recall some concepts on Wasserstein metric and the corresponding Wasserstein space. Such a space would serve as the basic underlying space for our regression analysis. Now let \((\Omega,d)\) be a complete, separable metric (Polish) space. Let \(p\in[1,\infty)\), for two probability measures \(\mu\) and \(\nu\) on \(\Omega\), the Wasserstein metric of order \(p\) is defined as

\[W_{p}(\mu,\nu) = \left(\inf_{\pi\in\Pi(\mu,\nu)}\int_{\Omega\times\Omega}d(x,y)^{p }d\pi(x,y)\right)^{1/p}\] \[= \inf\left\{\left[\mathbb{E}d(X,Y)^{p}\right]^{\frac{1}{p}},\quad \mathrm{law}(X)=\mu,\quad\mathrm{law}(Y)=\nu\right\},\]

where the infimum is taken over all pairs of random vectors \(X\) and \(Y\) marginally distributed as \(\mu\) and \(\nu\) (\(X\sim\mu\), \(Y\sim\nu\)). In above definition, \(\pi\) is also called a transference plan from \(\mu\) to \(\nu\) satisfying \(\int_{\Omega}d\pi(x,\cdot)=d\mu(x),\ \int_{\Omega}d\pi(\cdot,y)=d\nu(y)\). \(\Pi(\mu,\nu)\) is used to denote all such transference plans. For \(W_{p}\) defined above, earlier literatures have shown that it satisfies the basic axioms of a metric [34]. As an classical example, the distance \(W_{1}\), that is often called Kantorovich-Rubinstein metric, has played an important role in recent statistical learning problems [39, 38].

Now we consider a working space that is convenient for the study of our distribution regression model with FNNs. The space is the so-called Wasserstein space. Denote the collection of all Borel probability measures on \(\Omega\) by \(\mathcal{P}(\Omega)\). Let \(x_{0}\) be a fixed but arbitrarily chosen point in \(\Omega\), the Wasserstein space of order \(p\) is defined as

\[\mathcal{P}_{p}(\Omega)=\left\{\mu\in\mathcal{P}(\Omega):\quad\int_{\Omega}d \left(x_{0},x\right)^{p}d\mu(x)<+\infty\right\}.\]

This space does not depend on the choice of the point \(x_{0}\). Then \(W_{p}\) defines a (finite) metric on \(\mathcal{P}_{p}(\Omega)\).

It can be found in [34] that, for any fixed \(p\in[1,\infty)\), the metric space \((\mathcal{P}_{p}(\Omega),W_{p})\) is a compact metric space when \((\Omega,d)\) is a compact subset of \(\mathbb{R}^{n}\) equipped with the standard Euclidean metric \(d(x,y)=\|x-y\|_{2}\) as in aforementioned subsection. On the other hand, it is easy to see that when \(d(x,y)=\|x-y\|_{2}\), for any \(\mu\in\mathcal{P}(\Omega)\), there holds

\[\int_{\Omega}\|x_{0}-x\|_{2}^{p}\,d\mu(x)<+\infty\]

when \(\Omega\) is compact. Hence, under this setting, we have \(\mathcal{P}_{p}(\Omega)=\mathcal{P}(\Omega)\) and the metric space \((\mathcal{P}(\Omega),W_{p})\) is compact. In the following of the paper, we will use the notation that \(\Omega\) is a compact space with standard Euclidean metric. For fixed \(p\in[1,\infty)\) which will be selected depending on different learning situations, the distribution regression framework will be considered in such a compact metric space \((\mathcal{P}(\Omega),W_{p})\).

Inspired by the two-stage regression scheme for distribution regression, we propose our distribution regression model in Wasserstein space \((\mathcal{P}(\Omega),W_{p})\). The first-stage data set \(D=\{(\mu_{i},y_{i})\}_{i=1}^{m}\) are i.i.d. sampled from an unknown meta distribution that is a Borel probability measure \(\rho\) on \(\mathcal{Z}=\mathcal{U}\times\mathcal{Y}\), where \(\mathcal{U}=(\mathcal{P}(\Omega),W_{p})\) is the input metric space of
Borel probability measures on \(\Omega\) with Wasserstein metric, and \(\mathcal{Y}=[-M,M]\) is the output space with some \(M>0\). The regression function \(f_{\rho}\) on \(\mathcal{U}\) is defined as

\[f_{\rho}(\mu)=\int_{\mathcal{Y}}yd\rho(y|\mu), \tag{2.19}\]

where \(\rho(y|\mu)\) is the conditional distribution at \(\mu\) induced by \(\rho\), and it minimizes the mean squared error

\[\mathcal{E}(f)=\int_{\mathcal{Z}}(f(\mu)-y)^{2}d\rho.\]

We denote \(\rho_{\mathcal{U}}\) as the marginal distribution of \(\rho\) on \(\mathcal{U}\), and \(\left(L^{2}_{\rho_{\mathcal{U}}},\|\cdot\|_{\rho}\right)\) as the space of square integrable functions with respect to \(\rho_{\mathcal{U}}\). For convenience of later representations in some places, for a continuous function \(f\) defined on the compact space \((\mathcal{P}(\Omega),W_{p})\), we also use the notation

\[\|f\|_{\infty}=\sup_{\mu\in\mathcal{P}(\Omega)}|f(\mu)|. \tag{2.20}\]

The hypothesis space we use for the ERM algorithm follows the 3-layer FNN structure constructed in the proof of Theorem 2 with a positive constant \(R\) depending on \(d,Q,f,g\) associated with the fixed functional \(f\circ L^{Q}_{G}\) to bound the parameters in the FNN given by

\[\mathcal{H}_{(2,3),R,N}=\left\{c\cdot h^{(3)}(\mu):\|F^{(j)}\|_{\infty}\leq RN ^{2},\|b^{(j)}\|_{\infty}\leq R,\text{ for }j=1,2,3,\|c\|_{\infty}\leq RN\right\}, \tag{2.21}\]

where the \(\|\cdot\|_{\infty}\) norm for matrices is the maximum value of \(\ell_{1}\)-norms of its rows. For the hypothesis space \(\mathcal{H}_{(2,3),R,N}\), one of the most fundamental problems is whether it is well-defined as a compact hypothesis space for the learning theory of our distribution regression model? In the next theorem, we give a positive answer by rigorously proving the compactness of space \(\mathcal{H}_{(2,3),R,N}\), we show that \(\mathcal{H}_{(2,3),R,N}\) can be actually used as the hypothesis space of the distribution regression model (1.4). We use \((C(\mathcal{P}(\Omega)),\|\cdot\|_{\infty})\) to denote space of continuous functions on \((\mathcal{P}(\Omega),W_{p})\) equipped with the infinity norm defined in (2.20).

**Theorem 3**.: _The space \(\mathcal{H}_{(2,3),R,N}\) is a compact metric subspace of the space \((C(\mathcal{P}(\Omega)),\|\cdot\|_{\infty})\)._

Recall that for a subset \(\mathcal{H}\) of a normed space equipped with a norm \(\|\cdot\|\), the \(\epsilon\)-covering number \(\mathcal{N}(\mathcal{H},\epsilon,\|\cdot\|)\) is the minimum number of balls with radius \(\epsilon>0\) that cover \(\mathcal{H}\). After showing the compactness of \(\mathcal{H}_{(2,3),R,N}\), the covering number

\[\mathcal{N}(\mathcal{H}_{(2,3),R,N},\epsilon,\|\cdot\|_{\infty})\]

makes sense, is finite and would act as one of the main objects to establish other main results. Estimates of \(\mathcal{N}(\mathcal{H}_{(2,3),R,N},\epsilon,\|\cdot\|_{\infty})\) will be given in subsection 4.4.

Since for the distribution regression model, the concrete information of the first stage distribution samples are unavailable, we can only observe the second-stage data
\(\{(\{x_{i,j}\}_{j=1}^{n_{i}},y_{i})\}_{i=1}^{m}\), where \(\{x_{i,j}\in\Omega\}_{j=1}^{n_{i}}\) are i.i.d. sampled from the probability distribution \(\{\mu_{i}\}_{i=1}^{m}\), so we have to use the empirical distribution

\[\hat{\mu}_{i}^{n_{i}}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\delta_{x_{i,j}}\]

as the input of our FNN structure. Denote the empirical error with respect to the second-stage data \(\hat{D}\) as

\[\mathcal{E}_{\hat{D}}(f):=\frac{1}{m}\sum_{i=1}^{m}\left(f(\hat{\mu}_{i}^{n_{i }})-y_{i}\right)^{2},\]

then the empirical target function from the ERM algorithm using the hypothesis space (2.21) is the function in \(\mathcal{H}_{(2,3),R,N}\) that minimizes the empirical error:

\[f_{\hat{D},R,N}:=\arg\min_{f\in\mathcal{H}_{(2,3),R,N}}\mathcal{E}_{\hat{D}}(f). \tag{2.22}\]

The compactness of the hypothesis space \(\mathcal{H}_{(2,3),R,N}\) proved in Theorem 3 ensures the existence of a minimizer \(f_{\hat{D},R,N}\) of the above variational problem, and then the learning theory of the above distribution regression scheme makes sense in such a FNN-inducing space. We now define the projection operator \(\pi_{M}\) on the space of functional \(f:\mathcal{P}(\Omega)\to\mathbb{R}\) as

\[\pi_{M}(f)(\mu)=\begin{cases}M,&\text{ if }f(\mu)>M,\\ -M,&\text{ if }f(\mu)<-M,\\ f(\mu),&\text{ if }-M\leq f(\mu)\leq M.\end{cases}\]

Since the regression function \(f_{\rho}\) is bounded by \(M\), we use the truncated empirical target function

\[\pi_{M}f_{\hat{D},R,N}\]

as the final estimator.

We first derive an oracle inequality for the distribution regression scheme (2.22), which is different from those for the traditional one-stage statistical regression or kernel-based two-stage distribution regression. The oracle inequality for the traditional regression framework in [28, lemma 4] fails in the distribution regression framework, since the empirical target function of the ERM algorithm is learned from the second-stage data rather than the first-stage data. We utilize a novel two-stage error decomposition method for distribution regression by including the empirical error of the first-stage sample

\[\mathcal{E}_{D}(f):=\frac{1}{m}\sum_{i=1}^{m}\left(f(\mu_{i})-y_{i}\right)^{2}\]

as an intermediate term of the error decomposition. Such a two-stage error decomposition is new in the literature of learning theory. We present a novel error decomposition in the following proposition.
**Proposition 1**.: _Denote \(\mathcal{H}=\mathcal{H}_{(2,3),R,N}\), then for any \(h\in\mathcal{H}\) and \(f_{\hat{D},R,N}\) defined in (2.22),_

\[\mathcal{E}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}\left(f_ {\rho}\right)\leq\mathcal{E}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}_{D }\left(\pi_{M}f_{\hat{D},R,N}\right)+\mathcal{E}_{D}\left(\pi_{M}f_{\hat{D},R,N }\right)\] \[-\mathcal{E}_{\hat{D}}\left(\pi_{M}f_{\hat{D},R,N}\right)+ \mathcal{E}_{\hat{D}}\left(h\right)-\mathcal{E}_{D}\left(h\right)+\mathcal{E}_ {D}\left(h\right)-\mathcal{E}(h)+\mathcal{E}(h)-\mathcal{E}(f_{\rho}),\]

_which can be further bounded by \(I_{1}(D,\mathcal{H})+I_{2}(D,\mathcal{H})+\left|I_{3}(\hat{D},\mathcal{H}) \right|+\left|I_{4}(\hat{D},\mathcal{H})\right|+R(\mathcal{H})\), in which_

\[I_{1}(D,\mathcal{H}) =\left\{\mathcal{E}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E }\left(f_{\rho}\right)\right\}-\left\{\mathcal{E}_{D}\left(\pi_{M}f_{\hat{D}, R,N}\right)-\mathcal{E}_{D}\left(f_{\rho}\right)\right\},\] \[I_{2}(D,\mathcal{H}) =\left\{\mathcal{E}_{D}\left(h\right)-\mathcal{E}_{D}\left(f_{ \rho}\right)\right\}-\left\{\mathcal{E}(h)-\mathcal{E}\left(f_{\rho}\right) \right\},\] \[I_{3}(\hat{D},\mathcal{H}) =\mathcal{E}_{D}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}_{ \hat{D}}\left(\pi_{M}f_{\hat{D},R,N}\right),\] \[I_{4}(\hat{D},\mathcal{H}) =\mathcal{E}_{\hat{D}}\left(h\right)-\mathcal{E}_{D}\left(h\right),\qquad R(\mathcal{H})=\mathcal{E}(h)-\mathcal{E}(f_{\rho}).\]

Based on the two-stage error decomposition, we derive a new oracle inequality for the distribution regression scheme using our hypothesis space \(\mathcal{H}_{(2,3),R,N}\) in the following theorem.

**Theorem 4**.: _Consider the distribution regression framework described above with the first stage sample size \(m\), and the second stage sample size \(n_{1}=n_{2}=\cdots=n_{m}=n\). Then for any \(h\in\mathcal{H}_{(2,3),R,N}\) and \(\epsilon>0\), we have_

\[\text{Prob}\left\{\left\|\pi_{M}f_{\hat{D},R,N}-f_{\rho}\right\|_ {\rho}^{2}>2\left\|h-f_{\rho}\right\|_{\rho}^{2}+8\epsilon\right\}\] \[\leq \exp\left\{T_{1}N\log\frac{16M\widehat{R}}{\epsilon}+T_{2}N\log N -\frac{3m\epsilon}{2048M^{2}}\right\}\] \[+ \exp\left\{-\frac{m\epsilon^{2}}{2\left(3M+\left\|h\right\|_{ \infty}\right)^{2}\left(\left\|h-f_{\rho}\right\|_{\rho}^{2}+\frac{2}{3} \epsilon\right)}\right\}\] \[+ \exp\left\{\log 4m+T_{1}N\log\frac{80M\widehat{R}R^{2}N^{4}}{ \epsilon}+T_{2}N\log N-\frac{n\epsilon^{2}}{115200\max\{\left\|h\right\|_{ \infty}^{2},M^{2}\}R^{8}N^{16}\right\},\]

_where \(R\), \(\widehat{R}\) are constants depending on \(d,Q,g,f\), and \(T_{1}\), \(T_{2}\) are constants depending on \(d,Q\) given explicitly in the proof._

Utilizing the above oracle inequality for distribution regression, we are able to achieve almost optimal learning rates for the proposed distribution regression scheme on the hypothesis space \(\mathcal{H}_{(2,3),R,N}\) that includes the FNN structure we construct.

**Theorem 5**.: _If \(f_{\rho}=f\circ L_{G}^{Q}\) in the composite form (2.14), where \(Q\) is a polynomial of degree \(q\) on \(\Omega\), \(g\in C^{0,1}[-B_{Q},B_{Q}]\), and \(f\in C^{0,\beta}[-B_{G},B_{G}]\) for some \(0<\beta\leq 1\). If the first stage sample size \(m\) satisfies the restriction that_

\[\log\left(4m\right)\leq A_{6}m^{\frac{1}{2\beta+1}},\]
_and the neural network parameter \(N\) and the second stage sample size \(n\) are chosen by_

\[N=\left[A_{4}m^{\frac{1}{2\beta+1}}\right],\quad n\geq\left\lceil A_{5}m^{\frac{4 \beta+17}{2\beta+1}}\right\rceil,\]

_then there exists a constant \(A_{7}\) such that_

\[\mathbb{E}\left\{\mathcal{E}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E} \left(f_{\rho}\right)\right\}\leq A_{7}m^{-\frac{2\beta}{2\beta+1}}\log m,\]

_where \(A_{4},A_{5},A_{6}\) and \(A_{7}\) are constants depending on \(d,Q,g,f\) given explicitly in the proof._

Notice that in the proof of Theorem 4 and Theorem 5, we assume that all \(n_{i}=n\) for simplicity, but the proof still works for \(n_{i}\) chosen to be different, in which we just need to replace \(n\) by \(\min_{i}\{n_{i}\}\) in the proof, this also stresses the advantage of using distribution rather than vectors as inputs of FNNs for the distribution regression problem. We end this subsection by providing an intrinsic fact that the composite functional \(f\left(L_{G}(\mu)\right)\) should be continuous w.r.t. the probability variable \(\mu\). We put this fact in the following proposition and the proof will be given in the appendix.

**Proposition 2**.: _For any fixed \(p\in[1,\infty)\), the composite functional \(f\left(L_{G}(\mu)\right)\) with distribution variable \(\mu\) is continuous on the metric space \(\left(\mathcal{P}(\Omega),W_{p}\right)\)._

## 3 Related work and discussion

In this section, we review some related works on approximation and generalization performance of fully connected neural networks for regression with vector inputs, and some researches on distribution regression, to demonstrate the novelty and superiority of our results.

Dating back to thirty years ago, there have been many works analysing the approximation and generalization ability of fully connected neural networks on learning various classes of functions, wherein the Holder spaces are one of the most important classes of functions to be studied. At the beginning, most results about the approximation ability of the fully connected neural networks are derived using the \(C^{\infty}\) activation functions that satisfy two assumptions, that is for some \(b\in\mathbb{R}\), \(\sigma^{k}(b)\neq 0\) for any non-negative integer \(k\), and for some integer \(q\neq 1\), \(\lim_{u\rightarrow-\infty}\sigma(u)/|u|^{q}=0\) and \(\lim_{u\rightarrow\infty}\sigma(u)/u^{q}=1\). Such approximation rates can be found in [24] through a localized Taylor expansion method. It was shown that for any \(f\in W_{\infty}^{\beta}\left([-1,1]^{d}\right)\), there exists a net \(f_{N}\) such that \(\left\|f_{N}-f\right\|_{C([-1,1]^{d})}\leq c_{f,d,\beta}N^{-\beta/d}\), with a constant \(c_{f,d,\beta}\) independent of \(N\). ReLU was more efficient and effective in deep learning [15] and thus became more and more popular in practice, therefore it is desirable to find approximation results for deep ReLU neural networks. Such approximation results are derived in [36] that \(f\in W_{\infty}^{\beta}\left([-1,1]^{d}\right)\) can be approximated with an error \(\epsilon\), by using at most \(c\left(\log 1/\epsilon+1\right)\) layers and at most \(c\epsilon^{-d/\beta}\left(\log 1/\epsilon+1\right)\) weights and computation units. Based on these approximation results, generalization results for ERM algorithms using fully connected neural networks are then conducted in [5] and [28] to achieve optimal convergence rate in the traditional regression
framework. However, these research results are all for the traditional regression problem, the vector-input FNN does not work for the distribution regression framework with two stages of sampling due to the fact that the different sizes of second-stage sample will result in different dimensions of input vectors if we just vectorize all the second-stage samples. Furthermore, the generalization analysis for the traditional regression problem also fails in the distribution regression since the ERM algorithm is based on the second-stage data \(\hat{D}\) rather than the first-stage data \(D\), resulting in the desire of a new method to bound the concentration error between the generalization error and the second-stage empirical error for the two-stage sampling distribution regression framework.

To solve the limitation of the usage of vector-input neural network structure for the distribution regression framework, we propose a novel FNN structure that is able to take the empirical distribution as the input for the learning of distribution regression in practice. The idea of constructing the FNN structure with distribution input is slightly related to the work of [42], in which they transform the problem of learning symmetric functions to the functional learning problem with distribution input in a Wasserstein space, and construct a shallow net

\[f(\mu)=\frac{1}{m^{\prime}}\sum_{j^{\prime}=1}^{m^{\prime}}b_{j^{\prime}} \tilde{\sigma}\left(\frac{1}{m}\sum_{j=1}^{m}c_{j^{\prime},j}\int\sigma_{ \alpha}\left(\left\langle w_{j^{\prime},j},\tilde{x}\right\rangle\right)d \left(\mu(x)\right)\right),\quad\tilde{x}=[x,R]^{T}. \tag{3.1}\]

The architecture of the FNN for functional approximation and distribution regression we construct (Definition 1) is more general than (3.1) and can be deeper, which enables us to learn more complex functional. Utilizing our novel FNN structure, we construct a three-layer FNN to derive the approximation rate \(O\left(N^{-\beta}\right)\) for learning the functional \(f\left(L_{G}^{Q}(\mu)\right)\) induced by the polynomial \(Q\) and Holder functions \(f,g\), which is the same as the rate of learning the composite function \(f\left(g\left(Q(x)\right)\right)\) itself without taking into consideration the input of probability measures [28]. Actually, our FNN structure has the potential to learn functionals that are induced by more composite functions, where we just need to construct a deep FNN with the number of layers determined by the number of composite functions that induce the functionals. Furthermore, we utilize a novel error decomposition method to derive an oracle inequality for the distribution regression framework with two-stage sampling, and then combine it with the approximation rate to derive an almost optimal learning rate \(O\left(m^{-\frac{2\beta}{2\beta+1}}\log m\right)\) as [5] and [28].

Let us review some related works on distribution regression in the literature of learning theory [6, 31]. The recent main works on the two-stage distribution regression mainly include [33], [13], [11], [12], [25] and [37]. These works mainly rely on a mean embedding technique to transform probability measures to a space consisting of all mean embeddings via some Mercer kernel. In this theoretical model, the learning theory framework is in fact performed on the space of all mean embeddings which form a compact subspace of the continuous function space defined on the underlying space of the probability measures. To derive nice learning rates, the works of [33], [17], [13], [11], [12] and [37] deeply rely on the kernel regularization and the corresponding integral operator technique. In fact, to the best of our knowledge, learning rate for distribution regression are unexplored in
non-kernel or non-regularized settings where the integral operator approach cannot always be applied. Moveover, most of the classical kernel-based learning theory approaches such as [30], [3], [17], and [16] fail in the deep learning setting, and such classical techniques are often inappropriate for modern fully connected neural network structures. One of the main reasons is that the structure of FNN-based hypothesis spaces is totally different from that of the kernel-based hypothesis space. Hence, we are mainly in face of three aspects of difficulties in this work. The first one is how to provide an appropriate FNN structure with distribution inputs and a meaningful distribution regression model for the hypothesis space induced by the FNN. The second one is how to establish approximation results without utilizing the traditional kernel methods and its related integral operator approaches that the existing works on distribution regression used. The third one is how to realize distribution regression without using regularization in a novel hypothesis space that relies on FNN structures and derive optimal learning rates, since there is currently no theoretical result on distribution regression derived in such a setting.

In this paper, we overcome aforementioned difficulties and establish a learning theory framework of distribution regression with neural networks. At first, some new functional approximation results are proposed. Then we provide the proof of the crucial fact that the space \(\mathcal{H}_{(2,3),R,N}\) is compact. In contrast to the existing works on distribution regression, the underlying hypothesis space defined by us possesses nice FNN learning structures. Furthermore, though the classical effective dimension based analysis related to kernel methods in [33], [13], [25] and [37] fails in our setting, new estimates based on the covering number of the newly defined hypothesis space \(\mathcal{H}_{(2,3),R,N}\) are obtained. In contrast to the previous works utilizing the classical error decomposition for kernel based regression such as some earlier works [6, 35, 7] and recent developments [14, 19], another technical novelty is the crucial two-stage error decomposition which appears first in the literature of learning theory. Based on the functional approximation results, covering number estimates and the two-stage error decomposition, almost optimal learning rates of the proposed FNN distribution regression scheme are obtained. To the best of our knowledge, the work provides a beginning of the approximation theory of functionals with variable of probability measures using neural networks. It is also the first work to establish the learning theory of distribution regression and obtain almost optimal learning rates with neural networks. The methods in this work have the potential to be further applied to more general deep learning settings.

## 4 Proof of Main Results

### Proof of Theorem 1

In the construction of the FNN to approximate the functional, we need to use the following lemma to approximate Lipschitz continuous univariate functions using continuous piecewise linear functions (splines) \(\{\sigma(\cdot-t_{i})\}_{i=1}^{2N+3}\), with \(t_{i}=-1+\frac{i-2}{N}\), which can be found in [40] and [23].

**Lemma 2**.: _For \(N\in\mathbb{N}\), let \(\textbf{t}=\left\{t_{i}:=-1+\frac{i-2}{N}\right\}_{i=1}^{2N+3}\) be the uniform mesh on \(\left[\,-1-\frac{1}{N}\right.\)_
\(1+\frac{1}{N}\)], \(L_{t}\) be a linear operator on \(C[-B,B]\) given by_

\[L_{t}(g)=\frac{N}{B}\sum_{i=1}^{2N+3}\left(\mathcal{L}_{N}\left(\{g(Bt_{k})\}_{k =2}^{2N+2}\right)\right)_{i}\sigma\left(\cdot-Bt_{i}\right),\qquad g\in C[-B,B]. \tag{4.1}\]

_where \(\mathcal{L}_{N}:\mathbb{R}^{2N+1}\to\mathbb{R}^{2N+3}\) is a linear operator that for \(\zeta=(\zeta_{i})_{i=2}^{2N+2}\in\mathbb{R}^{2N+1}\)_

\[\left(\mathcal{L}_{N}(\zeta)\right)_{i}=\begin{cases}\zeta_{2},&\text{for }i=1, \\ \zeta_{3}-2\zeta_{2},&\text{for }i=2,\\ \zeta_{i-1}-2\zeta_{i}+\zeta_{i+1},&\text{for }3\leq i\leq 2N+1,\\ \zeta_{2N+1}-2\zeta_{2N+2},&\text{for }i=2N+2,\\ \zeta_{2N+2},&\text{for }i=2N+3.\end{cases} \tag{4.2}\]

_Then for \(g\in C^{0,\alpha}[-B,B]\) with \(0<\alpha\leq 1\), we have_

\[\left\|L_{\mathbf{t}}(g)-g\right\|_{C[-B,B]}\leq\frac{2B^{\alpha}|g|_{0,\alpha }}{N^{\alpha}},\]

_where \(|g|_{0,\alpha}\) is the semi-norm of the Lipschitz-\(\alpha\) continuous function \(g\)._

_Proof of Theorem 1._ For the first hidden layer of our FNN, we aim to realize the ridge functionals \(\left\{\int_{\Omega}\sigma(\xi\cdot x-B_{\xi}t_{j})d\mu\right\}_{j}\) leading to approximating the inner functional \(L_{G}^{\xi}(\mu)\). Take the connection matrix and bias vector as

\[F^{(1)}=\left[\begin{array}{c}\xi^{T}\\ \xi^{T}\\ \vdots\\ \xi^{T}\end{array}\right]\in\mathbb{R}^{(2N+3)\times d},\quad b^{(1)}=B_{\xi} \left[\begin{array}{c}t_{1}\\ t_{2}\\ \vdots\\ t_{2N+3}\end{array}\right]\in\mathbb{R}^{2N+3},\]

then the output of the first hidden layer \(h^{(1)}\in\mathbb{R}^{2N+3}\) is

\[\left(h^{(1)}(\mu)\right)_{j}=\int_{\Omega}\sigma\left(\xi\cdot x-B_{\xi}t_{j }\right)d\mu(x),\quad j=1,2,\ldots,2N+3. \tag{4.3}\]

The free parameters in the first hidden layer come only from \(\xi\) and \(B_{\xi}t_{j}\), thus the number of free parameters in this layer is

\[\mathcal{N}_{1}=d+2N+3.\]

The second hidden layer of our FNN aims to realize the ridge functions \(\left\{\sigma(\widetilde{\mathcal{S}}(\mu)-B_{G}t_{j})\right\}_{j}\) leading to approximating \(f\left(\widetilde{\mathcal{S}}(\mu)\right)\), where \(\widetilde{\mathcal{S}}(\mu)\) is the approximation of \(L_{G}^{\xi}(\mu)\) by the linear combination of \(h^{(1)}(\mu)\), and \(B_{G}\) is the upper bound of \(\left|\widetilde{\mathcal{S}}(\mu)\right|\). From Lemma 2, we know that there exists a linear combination of \(h^{(1)}(\mu)\) that approximates \(L_{G}^{\xi}(\mu)\), and denote this approximation as

\[\widetilde{\mathcal{S}}(\mu):=\frac{N}{B_{\xi}}\sum_{j=1}^{2N+3}\left(\mathcal{ L}_{N}\left(\left\{g\left(B_{\xi}t_{k}\right)\right\}_{k=2}^{2N+2}\right) \right)_{j}\int_{\Omega}\sigma\left(\xi\cdot x-B_{\xi}t_{j}\right)d\mu(x),\]
then according to Lemma 2, we have

\[\sup_{\mu\in\mathcal{P}(\Omega)}\left|\widetilde{\mathcal{S}}(\mu)-L _{G}^{\xi}(\mu)\right|\] \[\leq \sup_{\mu\in\mathcal{P}(\Omega)}\left|\frac{N}{B_{\xi}}\sum_{j=1}^{ 2N+3}\left(\mathcal{L}_{N}\left(\left\{g\left(B_{\xi}t_{k}\right)\right\}_{k=2 }^{2N+2}\right)\right)_{j}\int_{\Omega}\sigma\left(\xi\cdot x-B_{\xi}t_{j} \right)d\mu(x)-\int_{\Omega}g(\xi\cdot x)d\mu(x)\right|\] \[\leq \sup_{\mu\in\mathcal{P}(\Omega)}\int_{\Omega}\left\|\frac{N}{B_{ \xi}}\sum_{j=1}^{2N+3}\left(\mathcal{L}_{N}\left(\left\{g\left(B_{\xi}t_{k} \right)\right\}_{k=2}^{2N+2}\right)\right)_{j}\sigma\left(\xi\cdot x-B_{\xi}t _{j}\right)-g(\xi\cdot x)\right\|_{C(\Omega)}d\mu(x)\] \[\leq \frac{2B_{\xi}|g|_{0,1}}{N},\]

which also implies the bound of \(\left|\widetilde{\mathcal{S}}(\mu)\right|\),

\[\sup_{\mu\in\mathcal{P}(\Omega)}\left|\widetilde{\mathcal{S}}(\mu)\right|\leq \sup_{\mu\in\mathcal{P}(\Omega)}\left|L_{G}^{\xi}(\mu)\right|+\sup_{\mu\in \mathcal{P}(\Omega)}\left|\widetilde{\mathcal{S}}(\mu)-L_{G}^{\xi}(\mu) \right|\leq B_{g}+2B_{\xi}|g|_{0,1}=B_{G}.\]

Then we take the connection matrix and bias vector of the second hidden layer as

\[F^{(2)}=\frac{N}{B_{\xi}}\left[\begin{array}{c}\mathcal{L}_{N}^{T}\left( \left\{g\left(B_{\xi}t_{k}\right)\right\}_{k=2}^{2N+2}\right)\\ \mathcal{L}_{N}^{T}\left(\left\{g\left(B_{\xi}t_{k}\right)\right\}_{k=2}^{2N+ 2}\right)\\ \vdots\\ \mathcal{L}_{N}^{T}\left(\left\{g\left(B_{\xi}t_{k}\right)\right\}_{k=2}^{2N+ 2}\right)\end{array}\right]\in\mathbb{R}^{(2N+3)\times(2N+3)},\quad b^{(2)}=B_ {G}\left[\begin{array}{c}t_{1}\\ t_{2}\\ \vdots\\ t_{2N+3}\end{array}\right]\in\mathbb{R}^{2N+3}.\]

Therefore the output of the second hidden layer \(h^{(2)}(\mu)\in\mathbb{R}^{2N+3}\) is

\[\left(h^{(2)}(\mu)\right)_{j}=\sigma\left(\widetilde{\mathcal{S}}(\mu)-B_{G} t_{j}\right)_{j},\quad j=1,2,\ldots,2N+3. \tag{4.4}\]

The free parameters in the second hidden layer come only from \(\frac{N}{B_{\xi}}\mathcal{L}_{N}\left(\left\{g\left(B_{\xi}t_{k}\right)\right\} _{k=2}^{2N+2}\right)\) and \(B_{G}t_{j}\), thus the number of free parameters in this layer is

\[\mathcal{N}_{2}=2N+3+2N+3=4N+6.\]

Finally, choose \(c=\frac{N}{B_{G}}\mathcal{L}_{N}\left(\left\{f\left(B_{G}t_{k}\right)\right\}_{ k=2}^{2N+2}\right)\in\mathbb{R}^{2N+3}\), according to Lemma 2, we have

\[\sup_{\mu\in\mathcal{P}(\Omega)}\left|c\cdot h^{(2)}(\mu)-f(\widetilde{ \mathcal{S}}(\mu))\right|\leq\frac{2B_{G}^{\beta}|f|_{C^{0,\beta}}}{N^{\beta}}. \tag{4.5}\]

Then we can derive the approximation rate of \(f(L_{G}^{\xi}(\mu))\) with our constructed FNN by combining (4.4) and (4.5), that is

\[\sup_{\mu\in\mathcal{P}(\Omega)}\left|c\cdot h^{(2)}(\mu)-f(L_{G }^{\xi}(\mu))\right|\] \[\leq\sup_{\mu\in\mathcal{P}(\Omega)}\left|c\cdot h^{(2)}(\mu)-f( \widetilde{\mathcal{S}}(\mu))\right|+\left|f\right|_{C^{0,\beta}}\sup_{\mu\in \mathcal{P}(\Omega)}\left(\left|L_{G}^{\xi}(\mu)-\widetilde{\mathcal{S}}(\mu )\right|\right)^{\beta} \tag{4.6}\] \[\leq\frac{2B_{G}^{\beta}\left|f\right|_{C^{0,\beta}}}{N^{\beta}}+ \left|f\right|_{C^{0,\beta}}\left(\frac{2B_{\xi}|g|_{C^{0,1}}}{N}\right)^{ \beta}=\frac{2B_{G}^{\beta}\left|f\right|_{C^{0,\beta}}+\left|f\right|_{C^{0, \beta}}(2B_{\xi}|g|_{C^{0,1}})^{\beta}}{N^{\beta}}.\]
From the expression of the coefficient vector \(c=\frac{N}{B_{G}}\mathcal{L}_{N}\left(\{f\left(B_{G}t_{k}\right)\}_{k=2}^{2N+2}\right)\), we have known \(\|c\|_{\infty}\leq\frac{4\|f\|_{\infty}N}{B_{G}}\) and hence (2.5) holds.

The free parameters of our constructed FNN are from the first, second hidden layer and the coefficients \(c\), therefore, the total number of free parameters is

\[\mathcal{N}=\mathcal{N}_{1}+\mathcal{N}_{2}+2N+3=8N+d+12. \tag{4.7}\]

This completes the proof of Theorem 1.

Proof of Corollary 1.: First we define a function

\[g:[-B_{\xi},B_{\xi}]\rightarrow\mathbb{R},x\mapsto e^{-x}.\]

Then we know that \(\left|g\right|_{C^{0,1}}=e^{B_{\xi}}\) and \(B_{g}=\|g\|_{C[-B_{\xi},B_{\xi}]}=e^{B_{\xi}}\). Then the corresponding \(B_{G}\) in (2.2) can be explicitly written as

\[B_{G}=B_{g}+2B_{\xi}|g|_{C^{0,1}}=e^{B_{\xi}}+2B_{\xi}e^{B_{\xi}}.\]

Define a function \(f\) as the identity on the interval

\[[-(e^{B_{\xi}}+2B_{\xi}e^{B_{\xi}}),e^{B_{\xi}}+2B_{\xi}e^{B_{\xi}}]\]

which satisfies \(\left|f\right|_{C^{0,1}}=1\) and \(\|f\|_{\infty}=e^{B_{\xi}}+2B_{\xi}e^{B_{\xi}}\). We know that \(\|f\|_{\infty}=B_{G}\) and the functional \(\mathcal{L}_{\xi}\) defined in (2.7) can be exactly represented by

\[\mathcal{L}_{\xi}(\mu)=f\left(\int_{\Omega}g(\xi\cdot x)d\mu(x)\right)\]

with \(f\), \(g\) defined as above. Then the first result follows directly from Theorem 1 by substituting \(\beta=1\) and the above bounds.

The second inequality follows directly from Schwarz inequality \(|\xi\cdot x|\leq\|\xi\|_{2}\|x\|_{2}\leq 1\), which implies that \(B_{\xi}=\|\xi\cdot x\|_{C(\Omega)}\leq 1\). 

### Proof of Theorem 2

Proof of Theorem 2.: We first construct the two-hidden layer structures in the integral part, which aims to realize an approximation of \(g(Q(x))\) for \(x\in\Omega\). The first layer of this two-layer structure in the integral part aims to realize the ridge functions \(\{\xi_{i}\cdot x-t_{j}\}_{i,j}\) leading to approximating the polynomials \(\{(\xi_{i}\cdot x)^{\ell}\}_{i,\ell}\). Take the connection matrix as

\[F^{(1)}=\left[\begin{array}{c}\mathbf{1}_{2N+3}\xi_{1}^{T}\\ \mathbf{1}_{2N+3}\xi_{2}^{T}\\ \vdots\\ \mathbf{1}_{2N+3}\xi_{n_{q}}^{T}\end{array}\right]\in\mathbb{R}^{n_{q}(2N+3) \times d},\]

where \(\mathbf{1}_{2N+3}\) is the constant \(1\) vector in \(\mathbb{R}^{2N+3}\), and the bias vector \(b^{(1)}\in\mathbb{R}^{n_{q}(2N+3)}\) as

\[b^{(1)}_{(k-1)(2N+3)+j}=t_{j},\quad k=1,2,\ldots,n_{q},\ j=1,2,\ldots,2N+3,\]
then we have that

\[\left(\sigma(F^{(1)}x-b^{(1)})\right)_{(k-1)(2N+3)+j}=\sigma\left(\xi_{k}\cdot x -t_{j}\right),\quad k=1,2,\ldots,n_{q},\ j=1,2,\ldots,2N+3.\]

The free parameters in this layer are only \(\xi_{k}\) and \(t_{j}\), so the number of free parameters in the first hidden layer is

\[\mathcal{N}_{1}=n_{q}d+2N+3.\]

The second layer of this two-layer structure in the integral part aims to realize the ridge functions \(\left\{\sigma\left(\widetilde{Q}(x)-B_{Q}t_{j}\right)\right\}_{j}\), where \(\widetilde{Q}(x)\) is an approximation of \(Q(x)\), and \(B_{Q}\) is the upper bound of \(\left|\widetilde{Q}(x)\right|\) defined explicitly in (2.12). Then after the integral with respect to the input distribution \(\mu\), we get the output of the integral-part layer, which is then used to realize an approximation of the inner functional \(L_{G}^{Q}\).

Denote \(v^{[l]}=\mathcal{L}_{N}\left(\left\{t_{j}^{l}\right\}_{j=2}^{2N+2}\right)\in \mathbb{R}^{2N+3}\), \(\mathbf{V}=\left[v^{[1]}\ v^{[2]}\ \ldots\ v^{[q]}\right]^{T}\in\mathbb{R}^{q\times(2N+3)}\), \(O\) the zero matrix of size \(q\times(2N+3)\) and

\[F^{[N]}=N\left[\begin{array}{cccc}\mathbf{V}&O&\ldots&O\\ O&\mathbf{V}&\ldots&O\\ \vdots&\vdots&\vdots&\vdots\\ O&O&\ldots&\mathbf{V}\end{array}\right]\in\mathbb{R}^{qn_{q}\times n_{q}(2N+3)}.\]

Then we have

\[F^{[N]}\sigma\left(F^{(1)}x-b^{(1)}\right)=\left[P_{1}^{[1]}(x)\ \cdots\ P_{1}^{[q]}(x)\ P_{2}^{[1]}(x)\ \cdots\ P_{2}^{[q]}(x)\cdots P_{n_{q}}^{[1]}(x)\ \cdots\ P_{n_{q}}^{[q]}(x)\right]^{T},\]

where for \(k=1,2,\ldots,n_{q}\) and \(\ell=1,2,\ldots,q\),

\[P_{k}^{[\ell]}(x)=N\sum_{j=1}^{2N+3}v_{j}^{[\ell]}\sigma\left(\xi_{k}\cdot x-t _{j}\right).\]

If we denote

\[\gamma_{Q}=\left[\gamma_{1,1}^{Q}\ \gamma_{1,2}^{Q}\ \ldots\ \gamma_{1,q}^{Q}\ \gamma_{2,1}^{Q}\ \gamma_{2,2}^{Q}\ \ldots\ \gamma_{2,q}^{Q}\ \ldots\ \gamma_{n_{q},1}^{Q}\ \gamma_{n_{q},2}^{Q}\ \ldots\ \gamma_{n_{q},q}^{Q}\right],\]

then,

\[\gamma_{Q}F^{[N]}\sigma\left(F^{(1)}x-b^{(1)}\right)=\sum_{k=1}^{n_{q}}\sum_{ \ell=1}^{q}\gamma_{k,\ell}^{Q}N\sum_{j=1}^{2N+3}v_{j}^{[\ell]}\sigma\left(\xi _{k}\cdot x-t_{j}\right)\]

According to Lemma 2, for \(k=1,2,\ldots,n_{q}\), and \(\ell=1,2,\ldots,q\),

\[\sup_{x\in\Omega}\left|N\sum_{j=1}^{2N+3}v_{j}^{[\ell]}\sigma\left(\xi_{k} \cdot x-t_{j}\right)-\left(\xi_{k}\cdot x\right)^{\ell}\right|\leq\frac{2\ell }{N}, \tag{4.8}\]
denote

\[\widetilde{Q}(x)=Q(0)+\sum_{k=1}^{n_{q}}\sum_{\ell=1}^{q}\gamma_{k,\ell}^{Q}N \sum_{j=1}^{2N+3}v_{j}^{[\ell]}\sigma\left(\xi_{k}\cdot x-t_{j}\right),\]

then we have

\[\sup_{x\in\Omega}\left|\widetilde{Q}(x)-Q(x)\right|\leq\frac{2q\|\gamma_{Q}\|_{ 1}}{N}, \tag{4.9}\]

which also implies the upper bound of \(\left|\widetilde{Q}(x)\right|\) as

\[\sup_{x\in\Omega}\left|\widetilde{Q}(x)\right|\leq\sup_{x\in\Omega}|Q(x)|+\sup _{x\in\Omega}\left|\widetilde{Q}(x)-Q(x)\right|\leq\widehat{B}_{Q}+2q\|\gamma _{Q}\|_{1}=B_{Q}.\]

Take the connection matrix of the second hidden layer as

\[F^{(2)}=\mathbf{1}_{2N+3}\gamma_{Q}F^{[N]}\in\mathbb{R}^{(2N+3)\times n_{q}(2N +3)},\]

and the bias vector \(b^{(2)}\in\mathbb{R}^{2N+3}\) as

\[b_{j}^{(2)}=-Q(0)+B_{Q}t_{j},\quad\text{for }j=1,2,\ldots,2N+3.\]

Then the output of the integral-part layer \(h^{(2)}\in\mathbb{R}^{2N+3}\) is

\[\left(h^{(2)}(\mu)\right)_{j}=\int_{\Omega}\sigma\left(\widetilde{Q}(x)-B_{Q} t_{j}\right)d\mu(x),\quad\text{for }j=1,2,\ldots,2N+3. \tag{4.10}\]

The free parameters in this layer are only \(\gamma_{Q},Nv^{[\ell]}\) and \(-Q(0)+B_{Q}t_{j}\), so the number of free parameters in the second hidden layer is

\[\mathcal{N}_{2}=qn_{q}+q(2N+3)+2N+3=2(q+1)N+qn_{q}+3q+3.\]

After the realization of the polynomial \(Q\), the rest part of the proof is similar as that in the proof of Theorem 1. The third hidden layer of our FNN aims to realize the ridge functions \(\left\{\sigma(\widetilde{\mathcal{S}}(\mu)-B_{G}t_{j})\right\}_{j}\) leading to approximating \(f\left(\widetilde{\mathcal{S}}(\mu)\right)\), where \(\widetilde{\mathcal{S}}(\mu)\) is the approximation of \(L_{G}^{Q}(\mu)\) by the linear combination of \(h^{(2)}(\mu)\) that is defined in following equation, and \(B_{G}\) is the upper bound of \(\left|\widetilde{\mathcal{S}}(\mu)\right|\) defined explicitly in (2.13). Denote

\[\widetilde{\mathcal{S}}(\mu):=\frac{N}{B_{Q}}\sum_{j=1}^{2N+3}\left(\mathcal{ L}_{N}\left(\left\{g\left(B_{Q}t_{k}\right)\}_{k=2}^{2N+2}\right)\right)_{j} \int_{\Omega}\sigma\left(\widetilde{Q}(x)-B_{Q}t_{j}\right)d\mu(x),\]

then utilizing the error decomposition,

\[\sup_{\mu\in\mathcal{P}(\Omega)}\left|\widetilde{\mathcal{S}}(\mu )-L_{G}^{Q}(\mu)\right|\] \[\leq\sup_{\mu\in\mathcal{P}(\Omega)}\left|\frac{N}{B_{Q}}\sum_{j =1}^{2N+3}\left(\mathcal{L}_{N}\left(\left\{g\left(B_{Q}t_{k}\right)\}_{k=2}^{ 2N+2}\right)\right)_{j}\int_{\Omega}\sigma\left(\widetilde{Q}(x)-B_{Q}t_{j} \right)d\mu-\int_{\Omega}g(\widetilde{Q}(x))d\mu\right|\]
\[\left(h^{(3)}(\mu)\right)_{j}=\sigma\left(\widetilde{\mathcal{S}}(\mu)-B_{G}t_{j} \right)_{j},\quad j=1,2,\ldots,2N+3. \tag{4.12}\]

The free parameters in the third hidden layer are only \(\frac{N}{B_{Q}}\mathcal{L}_{N}\left(\left\{g\left(B_{\xi}t_{k}\right)\right\}_{ k=2}^{2N+2}\right)\) and \(B_{G}t_{j}\), thus the number of free parameters in this layer is

\[\mathcal{N}_{3}=2N+3+2N+3=4N+6.\]

Finally, choose \(c=\frac{N}{B_{G}}\mathcal{L}_{N}\left(\left\{f\left(B_{G}t_{k}\right)\right\}_{ k=2}^{2N+2}\right)\in\mathbb{R}^{2N+3}\), according to Lemma 2, we have

\[\sup_{\mu\in\mathcal{P}(\Omega)}\left|c\cdot h^{(3)}(\mu)-f(\widetilde{ \mathcal{S}}(\mu))\right|\leq\frac{2B_{G}^{\beta}\left|f\right|_{C^{0,\beta}}} {N^{\beta}}. \tag{4.13}\]
Then we can derive the approximation rate of our constructed FNN by combining (4.11) and (4.13), and using the Lipschitz-\(\beta\) condition of \(f\), that is

\[\begin{split}&\sup_{\mu\in\mathcal{P}(\Omega)}\Big{|}c\cdot h^{(3)}( \mu)-f(L_{G}^{Q}(\mu))\Big{|}\\ &\leq\sup_{\mu\in\mathcal{P}(\Omega)}\Big{|}c\cdot h^{(3)}(\mu)- f(\widetilde{\mathcal{S}}(\mu))\Big{|}+|f|_{C^{0,\beta}}\sup_{\mu\in\mathcal{P}( \Omega)}\Big{(}\Big{|}L_{G}^{Q}(\mu)-\widetilde{\mathcal{S}}(\mu)\Big{|} \Big{)}^{\beta}\\ &\leq\frac{2B_{G}^{\beta}\,|f|_{C^{0,\beta}}}{N^{\beta}}+|f|_{C^{ 0,\beta}}\left(\frac{3B_{Q}|g|_{C^{0,1}}}{N}\right)^{\beta}=\frac{2B_{G}^{ \beta}\,|f|_{C^{0,\beta}}+|f|_{C^{0,\beta}}\,(3B_{Q}|g|_{C^{0,1}})^{\beta}}{N ^{\beta}}.\end{split} \tag{4.14}\]

From the expression \(c=\frac{N}{B_{G}}\mathcal{L}_{N}\left(\{f\left(B_{G}t_{k}\right)\}_{k=2}^{2N+2 }\right)\), we know \(\|c\|_{\infty}\leq\frac{4\|f\|_{\infty}N}{B_{G}}\). Then we have (2.16). The free parameters of the whole network are from the first, second, third hidden layer and the coefficient vector \(c\), therefore the total number of free parameters is

\[\mathcal{N}=\mathcal{N}_{1}+\mathcal{N}_{2}+\mathcal{N}_{3}+2N+3=(2q+10)N+(d+ q)n_{q}+3q+15. \tag{4.15}\]

This proves Theorem 2. 

### Compactness of the hypothesis space \(\mathcal{H}_{(2,3),R,N}\)

In this section, we rigorously prove the compactness of \(\mathcal{H}_{(2,3),R,N}\). The argument is based on a general version of the Ascoli-Arzela theorem for an abstract framework [10, Theorem 19.1c].

We first recall some notations on general equi-continuity and equi-boundedness. Let \(\{f_{n}\}\) be a countable collection of continuous functions from a separable topological space \((X,\mathcal{T})\) into a metric space \((Y,d_{Y})\). The functions \(\{f_{n}\}\) are equibounded at \(x\) if the closure in \((Y,d_{Y})\) of the set \(\{f_{n}(x)\}\) is compact. The functions \(\{f_{n}\}\) are equi-continuous at a point \(x\in X\) if for every \(\varepsilon>0\), there exists an open set \(\mathcal{O}\in\mathcal{T}\) containing \(x\) and such that \(d_{Y}\left(f_{n}(x),f_{n}(y)\right)\leq\varepsilon\quad\) for all \(y\in\mathcal{O}\quad\) and all \(n\in\mathbb{N}\). Then the general Ascoli-Arzela theorem can be stated as follow.

**Lemma 3**.: _Let \(\{f_{n}\}\) be a sequence of continuous functions from a separable topological space \((X,\mathcal{T})\) into a metric space \((Y,d_{Y})\). Assume that the functions \(\{f_{n}\}\) are equibounded and equi-continuous at each \(x\in X\). Then, there exists a subsequence \(\{f_{n^{\prime}}\}\subset\{f_{n}\}\) and a continuous function \(f:X\to Y\) such that \(\{f_{n^{\prime}}\}\to f\) pointwise in \(X.\) Moreover the convergence is uniform on compact subsets of \(X.\)_

Proof of Theorem 3.: The idea is to show \(\mathcal{H}_{(2,3),R,N}\) is a sequentially compact subset of the metric space \((C(\mathcal{P}(\Omega)),\|\cdot\|_{\infty})\). Then \(\mathcal{H}_{(2,3),R,N}\) is compact. We have already known from aforementioned argument that the Wasserstein space \((\mathcal{P}(\Omega),W_{p})\) is a compact separable metric space. Hence, it is naturally a separable topological space with the topology induced by the Wasserstein metric \(\mathcal{T}=W_{p}\).

For any countable collection of functions \(\{f_{n}\}\subset\mathcal{H}_{(2,3),R,N}\) with

\[f_{n}:(\mathcal{P}(\Omega),W_{p})\rightarrow(\mathbb{R},|\cdot|). \tag{4.16}\]
We first show that \(\{f_{n}\}\) is equi-continuous at any point \(\mu\in{\cal P}(\Omega)\). From the structure of space \({\cal H}_{(2,3),R,N}\), we know that for any function \(f_{n}\) in the collection, there exist \(c_{f_{n}}\), \(F_{f_{n}}^{(j)}\), \(b_{f_{n}}^{(j)}\) with \(\|F_{f_{n}}^{(j)}\|_{\infty}\leq RN^{2}\), \(\|b_{f_{n}}^{(j)}\|_{\infty}\leq R\), \(\|c_{f_{n}}\|_{\infty}\leq RN\), for \(j=1,2,3\) and any \(n\in\mathbb{N}\) such that

\[f_{n}(\mu)=c_{f_{n}}\cdot\sigma\left(F_{f_{n}}^{(3)}\int_{\Omega}H_{f_{n}}^{(2 )}(x)d\mu(x)-b_{f_{n}}^{(3)}\right) \tag{4.17}\]

with

\[H_{f_{n}}^{(2)}(x)=\sigma\left(F_{f_{n}}^{(2)}\sigma(F_{f_{n}}^{(1)}x-b_{f_{n} }^{(1)})-b_{f_{n}}^{(2)}\right).\]

Then for any \(\mu,\nu\in{\cal P}(\Omega)\) and \(n\in\mathbb{N}\),

\[|f_{n}(\mu)-f_{n}(\nu)|\] \[=\left|c_{f_{n}}\cdot\left\{\sigma\left(F_{f_{n}}^{(3)}\int_{ \Omega}H_{f_{n}}^{(2)}(x)d\mu(x)-b_{f_{n}}^{(3)}\right)-\sigma\left(F_{f_{n}} ^{(3)}\int_{\Omega}H_{f_{n}}^{(2)}(x)d\nu(x)-b_{f_{n}}^{(3)}\right)\right\}\right|\] \[\leq\|c_{f_{n}}\|_{1}\|F_{f_{n}}^{(3)}\|_{\infty}\left\|\int_{ \Omega}H_{f_{n}}^{(2)}(x)d\mu(x)-\int_{\Omega}H_{f_{n}}^{(2)}(x)d\nu(x)\right\| _{\infty}. \tag{4.18}\]

The definition of the infinity norm for the vector implies that the above bound equals

\[\|c_{f_{n}}\|_{1}\|F_{f_{n}}^{(3)}\|_{\infty}\max_{1\leq i\leq 2N+3}\left| \int_{\Omega}\left(H_{f_{n}}^{(2)}(x)\right)_{i}d\mu(x)-\int_{\Omega}\left(H_ {f_{n}}^{(2)}(x)\right)_{i}d\nu(x)\right|\]

Recall the well-known duality formula for the Kantorovich-Rubinstein metric

\[W_{1}(\mu,\nu)=\sup_{\psi:\|\psi\|_{C^{0,1}}\leq 1}\left\{\int_{\Omega}\psi d \mu-\int_{\Omega}\psi d\nu\right\}.\]

Note that \(\|c_{f_{n}}\|_{1}\leq(2N+3)\|c_{f_{n}}\|_{\infty}\), we have (4.18) is bounded by

\[(2N+3)\|c_{f_{n}}\|_{\infty}\|F_{f_{n}}^{(3)}\|_{\infty}\max_{1\leq i\leq 2N +3}\left\|\left(H_{f_{n}}^{(2)}\right)_{i}\right\|_{C^{0,1}}W_{1}(\mu,\nu).\]

Also, for any \(x,y\in\Omega\) and \(n\in\mathbb{N}\), note that for any \(1\leq i\leq 2N+3\) there holds

\[\left|\left(H_{f_{n}}^{(2)}(x)\right)_{i}-\left(H_{f_{n}}^{(2)}(y )\right)_{i}\right|\leq\left\|H_{f_{n}}^{(2)}(x)-H_{f_{n}}^{(2)}(y)\right\|_{\infty}\] \[=\left\|\sigma\left(F_{f_{n}}^{(2)}\sigma(F_{f_{n}}^{(1)}x-b_{f_{ n}}^{(1)})-b_{f_{n}}^{(2)}\right)-\sigma\left(F_{f_{n}}^{(2)}\sigma(F_{f_{n}}^{(1)}y -b_{f_{n}}^{(1)})-b_{f_{n}}^{(2)}\right)\right\|_{\infty}\] \[\leq\|F_{f_{n}}^{(2)}\|_{\infty}\left\|F_{f_{n}}^{(1)}(x-y)\right\| _{\infty}\leq\|F_{f_{n}}^{(2)}\|_{\infty}\|F_{f_{n}}^{(1)}\|_{\infty}\|x-y\|_{ 2}.\]

Then it follows that

\[\max_{1\leq i\leq 2N+3}\left\|\left(H_{f_{n}}^{(2)}\right)_{i}\right\|_{C^{0,1} }\leq\|F_{f_{n}}^{(2)}\|_{\infty}\|F_{f_{n}}^{(1)}\|_{\infty}\]

and

\[|f_{n}(\mu)-f_{n}(\nu)|\leq(2N+3)\|c_{f_{n}}\|_{\infty}\prod_{j=1}^{3}\|F_{f_{ n}}^{(j)}\|_{\infty}W_{1}(\mu,\nu)\leq(2N+3)R^{4}N^{7}W_{1}(\mu,\nu),\ \forall\ n\in\mathbb{N}.\]
We know from [34] that there holds \(W_{1}(\mu,\nu)\leq W_{p}(\mu,\nu)\) for any \(p\in[1,\infty)\) which is in fact as a result of Holder's inequality, then we arrive at

\[|f_{n}(\mu)-f_{n}(\nu)|\leq(2N+3)R^{4}N^{7}W_{p}(\mu,\nu),\text{ for any }\mu,\nu\in\mathcal{P}(\Omega),n\in\mathbb{N}. \tag{4.19}\]

Then we know that, at any point \(\mu\in\mathcal{P}(\Omega)\), for every \(\epsilon>0\), there always exists an open ball in the space \((\mathcal{P}(\Omega),W_{p})\) centered at \(\mu\) in the form of

\[\mathcal{O}_{\mu,\epsilon}=\left\{\nu\in\mathcal{P}(\Omega):W_{p}(\nu,\mu)< \frac{\epsilon}{(2N+3)R^{4}N^{7}}\right\}\]

such that \(|f_{n}(\mu)-f_{n}(\nu)|\leq\epsilon\) for all \(\nu\in\mathcal{O}_{\mu,\epsilon}\) and all \(n\in\mathbb{N}\). That is to say, for any countable collection \(\{f_{n}\}\subset\mathcal{H}_{(2,3),R,N}\), the equi-continuity of \(\{f_{n}\}\) holds at any point of \(\mathcal{P}(\Omega)\).

For the equi-boundedness of \(\{f_{n}\}\), we only need to show that for any \(\mu\in\mathcal{P}(\Omega)\) and any \(n\in\mathbb{N}\), the function \(f_{n}\) in (4.17) as a map of (4.16) is uniformly bounded in \(\mathbb{R}\). To this end, since

\[f_{n}(\mu)=c_{f_{n}}\cdot h_{f_{n}}^{(3)}(\mu)\]

for some \(h_{f_{n}}^{(3)}\) defined layer by layer as in definition 1. Then for \(h_{f_{n}}^{(1)}(x):=\sigma\left(F_{f_{n}}^{(1)}x-b_{f_{n}}^{(1)}\right)\),

\[\left\|h_{f_{n}}^{(1)}(x)\right\|_{\infty}\leq\|F_{f_{n}}^{(1)}\|_{\infty}\|x \|_{\infty}+\|b_{f_{n}}^{(1)}\|_{\infty}\leq RN^{2}+R\leq 2RN^{2},\quad\forall n \in\mathbb{N}.\]

For \(h_{f_{n}}^{(2)}(\mu)\), we have

\[\left\|h_{f_{n}}^{(2)}(\mu)\right\|_{\infty}=\left\|\int_{\Omega} \sigma\left(F_{f_{n}}^{(2)}h_{f_{n}}^{(1)}(x)-b_{f_{n}}^{(2)}\right)d\mu \right\|_{\infty}\leq\int_{\Omega}\left\|\sigma\left(F_{f_{n}}^{(2)}h_{f_{n} }^{(1)}(x)-b_{f_{n}}^{(2)}\right)\right\|_{\infty}d\mu\\ \leq\|F_{f_{n}}^{(2)}\|_{\infty}\left\|h_{f_{n}}^{(1)}(x)\right\| _{\infty}+\|b_{f_{n}}^{(2)}\|_{\infty}\leq(2R^{2}+R)N^{4},\ \forall n\in\mathbb{N}.\]

Then for the third hidden layer,

\[\left\|h_{f_{n}}^{(3)}(\mu)\right\|_{\infty}\leq\|F_{f_{n}}^{(3)}\|_{\infty} \left\|h_{f_{n}}^{(2)}(\mu)\right\|_{\infty}+\|b_{f_{n}}^{(3)}\|_{\infty}\leq (2R^{3}+R^{2}+R)N^{6},\forall n\in\mathbb{N}.\]

Finally,

\[|f_{n}(\mu)|=|c_{f_{n}}\cdot h_{f_{n}}^{(3)}(\mu)|\leq\|c_{f_{n}}\|_{1}\|h_{f _{n}}^{(3)}(\mu)\|_{\infty}\leq(2N+3)(2R^{4}+R^{3}+R^{2})N^{7},\forall n\in \mathbb{N}. \tag{4.20}\]

From the above process, we have in fact shown that \(\sup_{\mu\in\mathcal{P}(\Omega)}|f_{n}(\mu)|\leq(2N+3)(2R^{4}+R^{3}+R^{2})N^{7}\) holds for all \(n\in\mathbb{N}\). Hence the equi-boundedness of \(\{f_{n}\}\) holds at any point of \(\mathcal{P}(\Omega)\). Combining the above arguments and using the fact that a subset of a metric space is compact if and only if it is a sequentially compact set, we obtain that \(\mathcal{H}_{(2,3),R,N}\) is compact.
### Covering number of space \(\mathcal{H}_{(2,3),R,N}\)

The generalization (estimation) error bound can be derived by a bias and variance trade-off method. The bias corresponds to the approximation rate of our proposed hypothesis space for approximating the target function, which is shown in Theorem 2, while we also need to demonstrate the approximation ability of our proposed hypothesis space by showing that the parameters of our constructed FNN can actually be bounded as required in the hypothesis space. The variance term can be measured by the complexity of the hypothesis space, such as pseudo dimension, Rademacher complexity and covering number, and we utilize the covering number as the tool of estimating the complexity. We first bound the parameters in our constructed FNN structure in the following lemma.

**Lemma 4**.: _Let \(Q\) be a polynomial of degree \(q\) on \(\Omega\), \(g\in C^{0,1}[-B_{Q},B_{Q}]\), \(f\in C^{0,\beta}[-B_{G},B_{G}]\), for some \(0<\beta\leq 1\). Then for the FNN constructed in the proof of Theorem 2, there exists a constant \(R=R_{d,Q,f,g}\) depending on \(d,Q,f,g\) such that for \(j=1,2,3\),_

\[\|F^{(j)}\|_{\infty}\leq RN^{2},\quad\|b^{(j)}\|_{\infty}\leq R,\quad\|c\|_{ \infty}\leq RN.\]

Proof.: For the first hidden layer, since \(|\xi|=1\), and \(|t_{j}|=1\), for all \(j=1,2,\ldots,2N+3\), we have that \(\|F^{(1)}\|_{\infty}\leq\sqrt{d}\) and \(\|b^{(1)}\|_{\infty}\leq 2\).

For the second hidden layer, since \(\|v^{(\ell)}\|_{\infty}\leq 4\), for all \(\ell=1,2,\ldots,q\), and \(\|Q\|_{C(\Omega)}=\widehat{B}_{Q}\leq B_{Q}\), we have that \(\|F^{(2)}\|_{\infty}\leq 4\|\gamma_{Q}\|_{1}N(2N+3)\) and \(\|b^{(2)}\|_{\infty}\leq 3B_{Q}\).

For the third hidden layer, since \(\|g\|_{C([-B_{Q},B_{Q}])}=B_{g}\leq B_{G}\), we have that \(\|F^{(3)}\|_{\infty}\leq\frac{4B_{G}}{B_{Q}}N(2N+3)\) and \(\|b^{(3)}\|_{\infty}\leq 2B_{G}\).

For the coefficient vector \(c\), we have that \(\|c\|_{\infty}\leq\frac{4\|f\|_{\infty}}{B_{G}}N\).

Thus we finish the proof by choosing

\[R=\max\left\{2\sqrt{d},20\|\gamma_{Q}\|_{1},3B_{Q},\frac{20B_{G}}{B_{Q}},2B_{G },\frac{4\|f\|_{\infty}}{B_{G}}\right\}.\]

In the following lemma, we bound the covering number of our hypothesis space \(\mathcal{H}_{(2,3),R,N}\), which is then used to derive an upper bound of the estimation error for the ERM algorithm.

**Lemma 5**.: _For \(R\geq 1,N\in\mathbb{N}\), and \(\widehat{R}=3R^{4}(10n_{q}+d+18)\), there exist two constants \(T_{1},T_{2}\) depending on \(d,q\) such that_

\[\log\mathcal{N}\left(\mathcal{H}_{(2,3),R,N},\epsilon,\|\cdot\|_{\infty} \right)\leq T_{1}N\log\frac{\widehat{R}}{\epsilon}+T_{2}N\log N.\]

Proof.: Denote \(h^{(1)}(x):=\sigma\left(F^{(1)}x-b^{(1)}\right)\), then

\[\left\|h^{(1)}(x)\right\|_{\infty}\leq\|F^{(1)}\|_{\infty}\|x\|_{\infty}+\|b^ {(1)}\|_{\infty}\leq RN^{2}+R\leq 2RN^{2}.\]
For the second layer, we have

\[\left\|h^{(2)}(\mu)\right\|_{\infty}=\left\|\int_{\Omega}\sigma \left(F^{(2)}h^{(1)}(x)-b^{(2)}\right)d\mu\right\|_{\infty}\leq\sup_{\mu\in \mathcal{P}(\Omega)}\int_{\Omega}\left\|\sigma\left(F^{(2)}h^{(1)}(x)-b^{(2)} \right)\right\|_{\infty}d\mu\] \[\leq\|F^{(2)}\|_{\infty}\left\|h^{(1)}(x)\right\|_{\infty}+\|b^{( 2)}\|_{\infty}\leq 2R^{2}N^{4}+R\leq 3R^{2}N^{4}.\]

Then for the third hidden layer,

\[\left\|h^{(3)}(\mu)\right\|_{\infty}\leq\|F^{(3)}\|_{\infty}\left\|h^{(2)}(\mu) \right\|_{\infty}+\|b^{(3)}\|_{\infty}\leq 3R^{3}N^{6}+R\leq 4R^{3}N^{6}.\]

If we choose another functional \(\widehat{c}\cdot\widehat{h}^{(3)}(\mu)\) in the hypothesis space \(\mathcal{H}_{(2,3),R,N}\) induced by \(\widehat{F}^{(j)}\), \(\widehat{b}^{(j)}\) and \(\widehat{c}\), satisfying the restriction that,

\[\left|F_{ik}^{(j)}-\widehat{F}_{ik}^{(j)}\right|\leq\epsilon,\ \left\|b^{(j)}- \widehat{b}^{(j)}\right\|_{\infty}\leq\epsilon,\ \text{and}\left\|c-\widehat{c}\right\|_{\infty}\leq\epsilon,\]

then by the Lipschitz property of ReLU,

\[\left\|h^{(1)}(x)-\widehat{h}^{(1)}(x)\right\|_{\infty}\leq\left\|F^{(1)}- \widehat{F}^{(1)}\right\|_{\infty}\left\|x\right\|_{\infty}+\left\|b^{(1)}- \widehat{b}^{(1)}\right\|_{\infty}\leq(d+1)\epsilon.\]

Consider the random variable \(\xi\), since \(\left|\mathbb{E}(\xi)\right|\leq\mathbb{E}\left(\left|\xi\right|\right)\), and \(\mathbb{E}(\xi)\) is bounded by the maximum value of \(\xi\), utilize the Lipschitz property of ReLU,

\[\left\|h^{(2)}(\mu)-\widehat{h}^{(2)}(\mu)\right\|_{\infty}\leq \left\|\int_{\Omega}\sigma\left(F^{(2)}h^{(1)}(x)-b^{(2)}\right)d\mu-\int_{ \Omega}\sigma\left(\widehat{F}^{(2)}\widehat{h}^{(1)}(x)-\widehat{b}^{(2)} \right)d\mu\right\|_{\infty}\] \[\leq \sup_{\mu\in\mathcal{P}(\Omega)}\int_{\Omega}\left\|\sigma \left(F^{(2)}h^{(1)}(x)-b^{(2)}\right)-\sigma\left(\widehat{F}^{(2)}\widehat{h }^{(1)}(x)-\widehat{b}^{(2)}\right)\right\|_{\infty}d\mu\] \[\leq \left\|F^{(2)}\left(h^{(1)}(x)-\widehat{h}^{(1)}(x)\right) \right\|_{\infty}+\left\|\left(F^{(2)}-\widehat{F}^{(2)}\right)\widehat{h}^{( 1)}(x)\right\|_{\infty}+\left\|b^{(2)}-\widehat{b}^{(2)}\right\|_{\infty}\] \[\leq RN^{2}(d+1)\epsilon+2RN^{2}n_{q}(2N+3)\epsilon+\epsilon\leq(1 0n_{q}+d+2)RN^{3}\epsilon,\]

where \(\|\cdot\|_{\infty}\) in above inequalities obeys the notations at the end of the introduction, specifically, \(\|\cdot\|_{\infty}\) on the right hand side of the first inequality is taken w.r.t. vector of functions with variables in \(\mathcal{P}(\Omega)\). Similarly, for the third hidden layer,

\[\left\|h^{(3)}(\mu)-\widehat{h}^{(3)}(\mu)\right\|_{\infty}\leq \left\|F^{(3)}\left(h^{(2)}(x)-\widehat{h}^{(2)}(x)\right)\right\|_{\infty}+ \left\|\left(F^{(3)}-\widehat{F}^{(3)}\right)\widehat{h}^{(2)}(x)\right\|_{\infty}\] \[+\left\|b^{(3)}-\widehat{b}^{(3)}\right\|_{\infty}\leq RN^{2}(10n _{q}+d+2)RN^{3}\epsilon+3R^{2}N^{4}(2N+3)\epsilon+\epsilon\] \[\leq(10n_{q}+d+18)R^{2}N^{5}\epsilon.\]

Finally, the output of these two functionals satisfy

\[\left\|c\cdot h^{(3)}(\mu)-\widehat{c}\cdot\widehat{h}^{(3)}(\mu )\right\|_{\infty}\leq\left\|c\cdot\left(h^{(3)}(\mu)-\widehat{h}^{(3)}(\mu )\right)\right\|_{\infty}+\left\|(c-\widehat{c})\cdot\widehat{h}^{(3)}(\mu) \right\|_{\infty}\] \[\leq(2N+3)RN(10n_{q}+d+18)R^{2}N^{5}\epsilon+(2N+3)4R^{3}N^{6} \epsilon\leq 5(10n_{q}+d+22)R^{3}N^{7}\epsilon=:\widetilde{\epsilon}. \tag{4.21}\]
Therefore, by taking a \(\epsilon\)-net of each free parameter in \(F^{(j)}\), \(b^{(j)}\) and \(c\), the \(\widetilde{\epsilon}\)-covering number of \(\mathcal{H}_{(2,3),R,N}\) can be bounded by

\[\mathcal{N}\left(\mathcal{H}_{(2,3),R,N},\widetilde{\epsilon}, \left\|\cdot\right\|_{\infty}\right)\leq\left\lceil\frac{2RN^{2}}{\epsilon} \right\rceil^{n_{q}d+qn_{q}+q(2N+3)+2N+3}\left\lceil\frac{2R}{\epsilon}\right \rceil^{6N+9}\left\lceil\frac{2RN}{\epsilon}\right\rceil^{2N+3} \tag{4.22}\] \[\leq \left(\frac{3R}{\epsilon}\right)^{(2q+10)N+(d+q)n_{q}+3q+15}N^{(4 q+6)N+2(d+q)n_{q}+6q+9}\leq\left(\frac{\widehat{R}}{\widetilde{\epsilon}} \right)^{T_{1}N}N^{T_{2}N},\]

where \(\widehat{R}=15R^{4}(10n_{q}+d+18)\), \(T_{1}=(d+q)n_{q}+5q+5\), and \(T_{2}=9(d+q)n_{q}+45q+190\). Thus we finish the proof by taking the logarithm.

### Proof of Theorem 4

We first prove the two-stage error decomposition which is crucial for the proof of Theorem 4.

Proof of Proposition 1.: With simple computations, we have

\[\mathcal{E}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}\left(f _{\rho}\right)=\mathcal{E}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}_{D }\left(\pi_{M}f_{\hat{D},R,N}\right)+\mathcal{E}_{D}\left(\pi_{M}f_{\hat{D},R, N}\right)\] \[-\mathcal{E}_{\hat{D}}\left(\pi_{M}f_{\hat{D},R,N}\right)+ \mathcal{E}_{\hat{D}}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}_{\hat{D} }\left(h\right)+\mathcal{E}_{\hat{D}}\left(h\right)-\mathcal{E}_{D}\left(h \right)+\mathcal{E}_{D}\left(h\right)-\mathcal{E}(h)\] \[+\mathcal{E}(h)-\mathcal{E}(f_{\rho}).\]

Then the first decomposition follows from the fact \(\mathcal{E}_{\hat{D}}\left(\pi_{M}f_{\hat{D},R,N}\right)\leq\mathcal{E}_{\hat {D}}\left(h\right)\). The second further decomposition follows immediately by inserting \(\pm\mathcal{E}(f_{\rho})\), \(\pm\mathcal{E}_{D}(f_{\rho})\) to the above terms.

We introduce several concentration inequalities based on Bernstein's inequality and Hoeffding's inequality which can be found in or derived from [6], [29] and [7].

**Lemma 6** (A-inequality).: _Let \(\mathcal{G}\) be a set of continuous functions on a probability space \(\mathcal{Z}\) such that, for some \(B>0,c>0\), \(|G-\mathbb{E}(G)|\leq B\) almost surely and \(\mathbb{E}\left(G^{2}\right)\leq c\mathbb{E}(G)\) for all \(f\in\mathcal{G}\), then for any \(\epsilon>0\) and \(0<\alpha\leq 1\),_

\[\text{Prob}\left\{\sup_{G\in\mathcal{G}}\frac{\mathbb{E}(G)-\frac{1}{m}\sum_{i= 1}^{m}G(z_{i})}{\sqrt{\mathbb{E}(G)+\epsilon}}>4\alpha\sqrt{\epsilon}\right\} \leq\mathcal{N}\left(\mathcal{G},\alpha\epsilon,\|\cdot\|_{\infty}\right)\exp \left\{-\frac{\alpha^{2}m\epsilon}{2c+\frac{2B}{3}}\right\}.\]

**Lemma 7** (B-inequality).: _Let \(\xi\) be a random variable on a probability space \(\mathcal{Z}\) with mean \(\mathbb{E}(\xi)\) and variance \(\sigma^{2}(\xi)=\sigma^{2}\). If \(|\xi(z)-\mathbb{E}(\xi)|\leq M_{\xi}\) for almost all \(z\in\mathcal{Z}\), then for any \(\epsilon>0\),_

\[\text{Prob}\left\{\frac{1}{m}\sum_{i=1}^{m}\xi(z_{i})-\mathbb{E}(\xi)> \epsilon\right\}\leq\exp\left\{-\frac{m\epsilon^{2}}{2\left(\sigma^{2}+\frac{1 }{3}M_{\xi}\epsilon\right)}\right\}.\]
**Lemma 8** (C-inequality).: _Let \(\mathcal{H}_{2}\) be a set of continuous functions on a probability space \(\mathcal{X}\) such that, for some \(M_{\mathcal{H}_{2}}>0\), \(|f-\mathbb{E}(f)|\leq M_{\mathcal{H}_{2}}\) almost surely for all \(f\in\mathcal{H}_{2}\), then for any \(\epsilon>0\),_

\[\text{Prob}\left\{\sup_{f\in\mathcal{H}_{2}}\left|\mathbb{E}(f)-\frac{1}{n} \sum_{i=1}^{n}f(x_{i})\right|>\epsilon\right\}\leq 2\mathcal{N}\left(\mathcal{H}_{2 },\frac{\epsilon}{4},\|\cdot\|_{\infty}\right)\exp\left\{-\frac{n\epsilon^{2} }{8M_{\mathcal{H}_{2}}^{2}}\right\}.\]

Now we prove the oracle inequality for the distribution regression framework based on the two-stage error decomposition method shown in Proposition 1 and these concentration inequalities.

Proof of Theorem 4.: We denote \(\mathcal{H}=\mathcal{H}_{(2,3),R,N}\), the two-stage error decomposition and the basic fact \(\left\|\pi_{M}f_{\hat{D},R,N}-f_{\rho}\right\|_{\rho}^{2}=\mathcal{E}\left( \pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}\left(f_{\rho}\right)\) imply that \(\left\|\pi_{M}f_{\hat{D},R,N}-f_{\rho}\right\|_{\rho}^{2}\) can be bounded by the sum of \(I_{1}(D,\mathcal{H})\), \(I_{2}(D,\mathcal{H})\), \(\left|I_{3}(\hat{D},\mathcal{H})\right|\), \(\left|I_{4}(\hat{D},\mathcal{H})\right|\) and \(R(\mathcal{H})\).

We first derive a bound for \(I_{1}(D,\mathcal{H})\) in a probability form. Consider the functional class

\[\mathcal{G}:=\left\{G=\left(\pi_{M}f(\mu)-y\right)^{2}-\left(f_{\rho}(\mu)-y \right)^{2}:f\in\mathcal{H}_{(2,3),R,N}\right\}.\]

For any fixed \(G\in\mathcal{G}\), there exists a \(f\in\mathcal{H}_{(2,3),R,N}\) such that \(G(z)=\left(\pi_{M}f(\mu)-y\right)^{2}-\left(f_{\rho}(\mu)-y\right)^{2}\), and

\[\mathbb{E}(G)=\mathcal{E}\left(\pi_{M}f\right)-\mathcal{E}\left(f_{\rho} \right)=\left\|\pi_{M}f-f_{\rho}\right\|_{\rho}^{2},\]

\[\frac{1}{m}\sum_{i=1}^{m}G(z_{i})=\mathcal{E}_{D}\left(\pi_{M}f\right)- \mathcal{E}_{D}\left(f_{\rho}\right).\]

Furthermore, since \(\left|\pi_{M}f(\mu)\right|\leq M\), \(\left|f_{\rho}(\mu)\right|\leq M\) and \(\left|y\right|\leq M\) almost surely, we have

\[\left|G(z)\right|=\left|\left(\pi_{M}f(\mu)-f_{\rho}(\mu)\right)\left(\pi_{M} f(\mu)+f_{\rho}(\mu)-2y\right)\right|\leq 8M^{2}.\]

Thus \(\left|G(z)-\mathbb{E}(G)\right|\leq 16M^{2}\) and \(\mathbb{E}\left(G^{2}\right)\leq 16M^{2}\left\|\pi_{M}f-f_{\rho}\right\|_{\rho}^{2}=16M^{2} \mathbb{E}(G)\). Moreover, since for any \(f_{1},f_{2}\in\mathcal{H}_{(2,3),R,N}\),

\[\left|\left(\pi_{M}f_{1}(\mu)-y\right)^{2}-\left(\pi_{M}f_{2}(\mu)-y\right)^ {2}\right|\leq 4M\left|\pi_{M}f_{1}(\mu)-\pi_{M}f_{2}(\mu)\right|\leq 4M\left|f_{1 }(\mu)-f_{2}(\mu)\right|,\]

an \(\frac{\epsilon}{4M}\)-covering of \(\mathcal{H}_{(2,3),R,N}\) results in an \(\epsilon\)-covering of \(\mathcal{G}\), we have

\[\mathcal{N}\left(\mathcal{G},\epsilon,\|\cdot\|_{\infty}\right)\leq\mathcal{N }\left(\mathcal{H}_{(2,3),R,N},\frac{\epsilon}{4M},\|\cdot\|_{\infty}\right).\]

Therefore, by applying A-inequality to \(\mathcal{G}\) with \(B=c=16M^{2}\) and \(\alpha=\frac{1}{4}\), we have that

\[\text{Prob}\left\{\sup_{f\in\mathcal{H}_{(2,3),R,N}}\frac{\mathcal{ E}\left(\pi_{M}f\right)-\mathcal{E}\left(f_{\rho}\right)-\left(\mathcal{E}_{D} \left(\pi_{M}f\right)-\mathcal{E}_{D}\left(f_{\rho}\right)\right)}{\sqrt{ \mathcal{E}\left(\pi_{M}f\right)-\mathcal{E}\left(f_{\rho}\right)+\epsilon}}> \sqrt{\epsilon}\right\}\] \[\leq\mathcal{N}\left(\mathcal{G},\frac{\epsilon}{4},\|\cdot\|_{ \infty}\right)\exp\left\{-\frac{3m\epsilon}{2048M^{2}}\right\}\leq\mathcal{N} \left(\mathcal{H}_{(2,3),R,N},\frac{\epsilon}{16M},\|\cdot\|_{\infty}\right) \exp\left\{-\frac{3m\epsilon}{2048M^{2}}\right\}.\]
Since \(\sqrt{\epsilon\left(\mathcal{E}\left(\pi_{M}f\right)-\mathcal{E}\left(f_{\rho} \right)+\epsilon\right)}\leq\frac{1}{2}\left(\mathcal{E}\left(\pi_{M}f\right)- \mathcal{E}\left(f_{\rho}\right)\right)+\epsilon\), by choosing \(f=f_{\hat{D},R,N}\), we conclude that

\[\begin{array}{l}\mathit{Prob}\left\{I_{1}(D,\mathcal{H})>\frac{1}{2}\left( \mathcal{E}\left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}\left(f_{\rho} \right)\right)+\epsilon\right\}\\ \leq\mathcal{N}\left(\mathcal{H}_{(2,3),R,N},\frac{\epsilon}{16M},\| \cdot\|_{\infty}\right)\exp\left\{-\frac{3m\epsilon}{2048M^{2}}\right\}.\end{array} \tag{4.23}\]

Next we derive a bound for \(I_{2}(D,\mathcal{H})\) in a probability form. Consider the random variable \(\xi\) on \(\mathcal{Z}\) defined by

\[\xi(z)=\left(h(\mu)-y\right)^{2}-\left(f_{\rho}(\mu)-y\right)^{2}.\]

Since \(|f_{\rho}(\mu)|\leq M\) and \(|y|\leq M\) almost surely, we have

\[\begin{array}{l}|\xi(z)|\leq\left(3M+\|h\|_{\infty}\right)^{2},\quad|\xi- \mathbb{E}(\xi)|\leq 2\left(3M+\|h\|_{\infty}\right)^{2},\\ \sigma^{2}\leq\mathbb{E}(\xi^{2})\leq\left(3M+\|h\|_{\infty}\right)^{2}R( \mathcal{H}).\end{array}\]

Hence, by applying B-inequality with \(M_{\xi}=2\left(3M+\|h\|_{\infty}\right)^{2}\), we conclude that

\[\mathit{Prob}\left\{I_{2}(D,\mathcal{H})>\epsilon\right\}\leq\exp\left\{- \frac{m\epsilon^{2}}{2\left(\sigma^{2}+\frac{1}{3}M_{\xi}\epsilon\right)} \right\}\leq\exp\left\{-\frac{m\epsilon^{2}}{2\left(3M+\|h\|_{\infty}\right) ^{2}\left(R(\mathcal{H})+\frac{2}{3}\epsilon\right)}\right\}. \tag{4.24}\]

Then we derive a bound for \(\left|I_{3}(\hat{D},\mathcal{H})\right|\) in a probability form. Since \(\left\|\pi_{M}f_{\hat{D},R,N}\right\|_{\infty}\leq M\) and \(|y_{i}|\leq M\),

\[\begin{array}{l}\left|I_{3}(\hat{D},\mathcal{H})\right|=&\left|\frac{1}{m} \sum_{i=1}^{m}\left(\pi_{M}f_{\hat{D},R,N}\left(\hat{\mu}_{i}^{n}\right)-y_{i} \right)^{2}-\left(\pi_{M}f_{\hat{D},R,N}\left(\mu_{i}\right)-y_{i}\right)^{2 }\right|\\ \leq&\frac{1}{m}\sum_{i=1}^{m}4M\left|f_{\hat{D},R,N}\left(\hat{\mu }_{i}^{n}\right)-f_{\hat{D},R,N}\left(\mu_{i}\right)\right|,\end{array}\]

Then we only need to bound \(\left|f_{\hat{D},R,N}\left(\hat{\mu}_{i}^{n}\right)-f_{\hat{D},R,N}\left(\mu_ {i}\right)\right|\). Recall that the structure of the FNN in our hypothesis space \(\mathcal{H}_{(2,3),R,N}\) includes the integration w.r.t. the input distribution \(\mu\) after the second layer, that is, for any \(f\in\mathcal{H}_{(2,3),R,N}\), there always exists a vector of functions

\[H_{f}^{(2)}(x)=\sigma\left(F_{f}^{(2)}\sigma(F_{f}^{(1)}x-b_{f}^{(1)})-b_{f}^ {(2)}\right) \tag{4.25}\]

and \(c_{f}\), \(F_{f}^{(3)}\), \(b_{f}^{(3)}\) satisfying \(\|F_{f}^{(j)}\|_{\infty}\leq RN^{2}\), \(\|b_{f}^{(j)}\|_{\infty}\leq R\), \(\|c_{f}\|_{\infty}\leq RN\), \(j=1,2,3\) such that

\[f(\mu)=c_{f}\cdot\sigma\left(F_{f}^{(3)}h_{f}^{(2)}(\mu)-b_{f}^{(3)}\right)=c_ {f}\cdot\sigma\left(F_{f}^{(3)}\int_{\Omega}H_{f}^{(2)}(x)d\mu(x)-b_{f}^{(3)} \right).\]

If we denote a function class \(\mathcal{H}_{2}\) by

\[\begin{array}{l}\mathcal{H}_{2}=\Big{\{}H_{2}(x)=\left(\sigma\left(F^{(2)} \sigma\left(F^{(1)}x-b^{(1)}\right)-b^{(2)}\right)\right)_{1}:\\ \|F^{(j)}\|_{\infty}\leq RN^{2},\|b^{(j)}\|_{\infty}\leq R,j=1,2\Big{\}},\end{array} \tag{4.26}\]
which is just the first element of the output a second-layer network in the part of the hypothesis space \(\mathcal{H}_{(2,3),R,N}\). Therefore, for any \(f\in\mathcal{H}_{(2,3),R,N}\), there exists a function \(H_{f}\in\mathcal{H}_{2}\), with

\[\left(h_{f}^{(2)}(\mu_{i})\right)_{1}=\mathbb{E}_{x_{i}\sim\mu_{i}}\left(H_{f} \left(x_{i}\right)\right),\quad i=1,2,...,m,\]

and

\[\left(h_{f}^{(2)}(\hat{\mu}_{i}^{n})\right)_{1}=\frac{1}{n}\sum_{j=1}^{n}H_{f} (x_{i,j}),\quad i=1,2,...,m.\]

From the proof of Lemma 5 or Lemma 9 (in Appendix), we know that \(\sup_{H_{2}\in\mathcal{H}_{2}}|H_{2}|\leq 3R^{2}N^{4}\), thereby \(\sup_{H_{2}\in\mathcal{H}_{2}}|H_{2}-\mathbb{E}(H_{2})|\leq 6R^{2}N^{4}\). Hence, by applying C-inequality to \(\mathcal{H}_{2}\) with \(M_{\mathcal{H}_{2}}=6R^{2}N^{4}\) and using the above notations, we conclude that

\[Prob\left\{\sup_{f\in\mathcal{H}_{(2,3),R,N}}\left|\left(h_{f}^{ (2)}(\mu_{i})\right)_{1}-\left(h_{f}^{(2)}(\hat{\mu}_{i}^{n})\right)_{1}\right| >\epsilon\right\}\] \[= Prob\left\{\sup_{f\in\mathcal{H}_{(2,3),R,N}}\left|\mathbb{E}_{x _{i}\sim\mu_{i}}\left(H_{f}\left(x_{i}\right)\right)-\frac{1}{n}\sum_{j=1}^{n }H_{f}(x_{i,j})\right|>\epsilon\right\}\] \[\leq Prob\left\{\sup_{H_{2}\in\mathcal{H}_{2}}\left|\mathbb{E}_{x_{i} \sim\mu_{i}}\left(H_{2}\left(x_{i}\right)\right)-\frac{1}{n}\sum_{j=1}^{n}H_{ 2}(x_{i,j})\right|>\epsilon\right\}\] \[\leq 2\mathcal{N}\left(\mathcal{H}_{2},\frac{\epsilon}{4},\|\cdot\|_ {C(\Omega)}\right)\exp\left\{-\frac{n\epsilon^{2}}{288R^{4}N^{8}}\right\}.\]

From the structure of the hypothesis space \(\mathcal{H}_{(2,3),R,N}\), since each row of the connection matrix \(F^{(j)}\) and bias \(b^{(j)}\) satisfies the same symmetric restriction, we observe from the following explicit form

\[\sup_{\begin{subarray}{c}\left\|F^{(j)}\right\|_{\infty}\leq RN^ {2}\\ \left\|\phi^{(j)}\right\|_{\infty}\leq R,j=1,2\end{subarray}} \bigg{|}\int_{\Omega}\left(\sigma\left(F^{(2)}\sigma(F^{(1)}x-b^{( 1)})-b^{(2)}\right)\right)_{s}d\mu_{i}(x)\] \[-\frac{1}{n}\sum_{j=1}^{n}\left(\sigma\left(F^{(2)}\sigma(F^{(1) }x_{ij}-b^{(1)})-b^{(2)}\right)\right)_{s}\bigg{|}\]

that for any fixed \(\mu_{i}\) and \(\hat{\mu}_{i}=\frac{1}{n}\sum_{j=1}^{n}\delta_{x_{ij}}\), the sup operation w.r.t. \(F^{(j)}\) and \(b^{(j)}\) contribute equally to any \(s\)-th row of the \(2N+3\) rows. Then we know that

\[\sup_{f\in\mathcal{H}_{(2,3),R,N}}\left|\left(h_{f}^{(2)}(\mu_{i})\right)_{1}- \left(h_{f}^{(2)}(\hat{\mu}_{i}^{n})\right)_{1}\right|\leq\epsilon\]

implies that

\[\sup_{f\in\mathcal{H}_{(2,3),R,N}}\left|\left(h_{f}^{(2)}(\mu_{i})\right)_{s}- \left(h_{f}^{(2)}(\hat{\mu}_{i}^{n})\right)_{s}\right|\leq\epsilon,\quad\forall s =1,2,\ldots,2N+3,\]
which then implies that \(\sup_{f\in\mathcal{H}_{(2,3),R,N}}\left\|h_{f}^{(2)}(\mu_{i})-h_{f}^{(2)}(\hat{ \mu}_{i}^{n})\right\|_{\infty}\leq\epsilon\). Then it follows that

\[\sup_{f\in\mathcal{H}_{(2,3),R,N}}|f(\mu_{i})-f(\hat{\mu}_{i}^{n})|\] \[\leq \sup_{f\in\mathcal{H}_{(2,3),R,N}}\|c_{f}\|_{1}\left\|\sigma\left(F _{f}^{(3)}h_{f}^{(2)}(\mu_{i})-b_{f}^{(3)}\right)-\sigma\left(F_{f}^{(3)}h_{f} ^{(2)}(\hat{\mu}_{i}^{n})-b_{f}^{(3)}\right)\right\|_{\infty}\] \[\leq (2N+3)\sup_{f\in\mathcal{H}_{(2,3),R,N}}\|c_{f}\|_{\infty}\left\| \sigma\left(F_{f}^{(3)}h_{f}^{(2)}(\mu_{i})-b_{f}^{(3)}\right)-\sigma\left(F_ {f}^{(3)}h_{f}^{(2)}(\hat{\mu}_{i}^{n})-b_{f}^{(3)}\right)\right\|_{\infty}\] \[\leq (2N+3)RN\|F_{f}^{(3)}\|_{\infty}\sup_{f\in\mathcal{H}_{(2,3),R,N }}\left\|h_{f}^{(2)}(\mu_{i})-h_{f}^{(2)}(\hat{\mu}_{i}^{n})\right\|_{\infty} \leq 5R^{2}N^{4}\epsilon.\]

Thus we have that, for \(i=1,2,...,m\),

\[\text{\it Prob}\left\{\sup_{f\in\mathcal{H}_{(2,3),R,N}}|f(\mu_{i})-f(\hat{ \mu}_{i}^{n})|>5R^{2}N^{4}\epsilon\right\}\leq 2\mathcal{N}\left(\mathcal{H}_{2}, \frac{\epsilon}{4},\|\cdot\|_{C(\Omega)}\right)\exp\left\{-\frac{n\epsilon^{2 }}{288R^{4}N^{8}}\right\}.\]

Now combine all the \(m\) terms with \(i=1,2,\ldots,m\), we have

\[\text{\it Prob}\left\{\left|I_{3}(\hat{D},\mathcal{H})\right|>20 MR^{2}N^{4}\epsilon\right\}\] \[\leq \text{\it Prob}\left\{\frac{1}{m}\sum_{i=1}^{m}4M\left|f_{\hat{D},R,N}\left(\mu_{i}\right)-f_{\hat{D},R,N}\left(\hat{\mu}_{i}^{n}\right)\right| >20MR^{2}N^{4}\epsilon\right\}\] \[\leq \sum_{i=1}^{m}\text{\it Prob}\left\{\left|f_{\hat{D},R,N}(\mu_{i })-f_{\hat{D},R,N}(\hat{\mu}_{i}^{n})\right|>5R^{2}N^{4}\epsilon\right\}\] \[\leq \sum_{i=1}^{m}\text{\it Prob}\left\{\sup_{f\in\mathcal{H}_{(2,3),R,N}}|f(\mu_{i})-f(\hat{\mu}_{i}^{n})|>5R^{2}N^{4}\epsilon\right\}\] \[\leq 2m\mathcal{N}\left(\mathcal{H}_{2},\frac{\epsilon}{4},\|\cdot\|_ {C(\Omega)}\right)\exp\left\{-\frac{n\epsilon^{2}}{288R^{4}N^{8}}\right\}.\]

Finally, replacing \(\epsilon\) by \(\frac{\epsilon}{20MR^{2}N^{4}}\), we get

\[\text{\it Prob}\left\{\left|I_{3}(\hat{D},\mathcal{H})\right|>\epsilon\right\} \leq 2m\mathcal{N}\left(\mathcal{H}_{2},\frac{\epsilon}{80MR^{2}N^{4}},\|\cdot\|_ {C(\Omega)}\right)\exp\left\{-\frac{n\epsilon^{2}}{115200M^{2}R^{8}N^{16}} \right\}. \tag{4.27}\]

Similarly, we have

\[\text{\it Prob}\left\{\left|I_{4}(\hat{D},\mathcal{H})\right|> \epsilon\right\}\leq 2m\mathcal{N}\left(\mathcal{H}_{2},\frac{\epsilon}{80MR^{2}N^{4}}, \|\cdot\|_{C(\Omega)}\right) \tag{4.28}\] \[\exp\left\{-\frac{n\epsilon^{2}}{28800\left(\|h\|_{\infty}+M\right) ^{2}R^{8}N^{16}}\right\}.\]
Finally, combining (4.23), (4.24), (4.27) and (4.28), we have

\[\mathit{Prob}\left\{\left\|\pi_{M}f_{\hat{D},R,N}-f_{\rho}\right\|_ {\rho}^{2}>2\left\|h-f_{\rho}\right\|_{\rho}^{2}+8\epsilon\right\}\] \[\leq \mathcal{N}\left(\mathcal{H}_{(2,3),R,N},\frac{\epsilon}{16M}, \left\|\cdot\right\|_{\infty}\right)\exp\left\{-\frac{3m\epsilon}{2048M^{2}}\right\}\] \[+\exp\left\{-\frac{m\epsilon^{2}}{2\left(3M+\left\|h\right\|_{ \infty}\right)^{2}\left(\left\|h-f_{\rho}\right\|_{\rho}^{2}+\frac{2}{3} \epsilon\right)}\right\}\] \[+4m\mathcal{N}\left(\mathcal{H}_{2},\frac{\epsilon}{80MR^{2}N^{4 }},\left\|\cdot\right\|_{C(\Omega)}\right)\exp\left\{-\frac{n\epsilon^{2}}{11 5200\max\{\left\|h\right\|_{\infty}^{2},M^{2}\}R^{8}N^{16}}\right\}.\]

Thus we complete the proof by utilizing the covering number bound of the set \(\mathcal{H}_{(2,3),R,N}\) and \(\mathcal{H}_{2}\) from Lemma 5 and Lemma 9. 

### Proof of Theorem 5

Based on the oracle inequality for distribution regression in Theorem 4, we are now ready to give the proof of Theorem 5.

Proof of Theorem 5.: According to Theorem 2, there exists \(h^{*}\in\mathcal{H}_{(2,3),R,N}\), such that

\[\left\|h^{*}-f_{\rho}\right\|_{\rho}\leq\sup_{\mu\in\mathcal{P}(\Omega)}\left| h^{*}(\mu)-f_{\rho}(\mu)\right|\leq C_{1}N^{-\beta},\]

where \(C_{1}=2B_{G}^{\beta}\left|f\right|_{C^{0,\beta}}+\left(3B_{Q}|g|_{C^{0,1}} \right)^{\beta}\left|f\right|_{C^{0,\beta}}\), which also implies that

\[\left\|h^{*}\right\|_{\infty}\leq\sup_{\mu\in\mathcal{P}(\Omega)}\left|h^{*}( \mu)-f_{\rho}(\mu)\right|+\sup_{\mu\in\mathcal{P}(\Omega)}\left|f_{\rho}(\mu) \right|\leq C_{1}+M.\]

Then utilizing Theorem 4 by taking \(h=h^{*}\), it follows that

\[\mathit{Prob}\left\{\left\|\pi_{M}f_{\hat{D},R,N}-f_{\rho}\right\| _{\rho}^{2}>2C_{1}^{2}N^{-2\beta}+8\epsilon\right\}\] \[\leq \exp\left\{T_{1}N\log\frac{16M\widehat{R}}{\epsilon}+T_{2}N\log N -\frac{3m\epsilon}{2048M^{2}}\right\}\] \[+\exp\left\{-\frac{m\epsilon^{2}}{2\left(4M+C_{1}\right)^{2} \left(C_{1}^{2}N^{-2\beta}+\frac{2}{3}\epsilon\right)}\right\}\] \[+\exp\left\{\log 4m+T_{1}N\log\frac{80M\widehat{R}R^{2}N^{4}}{ \epsilon}+T_{2}N\log N-\frac{n\epsilon^{2}}{115200\left(M+C_{1}\right)^{2}R^{ 8}N^{16}}\right\}.\]
If we restrict \(\epsilon\geq 2C_{1}^{2}N^{-2\beta}\log N\), we have

\[\begin{split}&\mathit{Prob}\left\{\|\pi_{M}f_{\hat{D},R,N}-f_{ \rho}\|_{\rho}^{2}>9\epsilon\right\}\leq\mathit{Prob}\left\{\|\pi_{M}f_{\hat{D},R,N}-f_{\rho}\|_{\rho}^{2}>2C_{1}^{2}N^{-2\beta}+8\epsilon\right\}\\ \leq&\exp\left\{T_{1}N\log\frac{8M\widehat{R}N^{2 \beta}}{C_{1}^{2}}+T_{2}N\log N-\frac{3m\epsilon}{2048M^{2}}\right\}+\exp \left\{-\frac{3m\epsilon}{8\left(4M+C_{1}\right)^{2}}\right\}\\ +&\exp\left\{\log 4m+T_{1}N\log\frac{40M\widehat{R}R^{ 2}N^{2\beta+4}}{C_{1}^{2}}+T_{2}N\log N-\frac{n\epsilon^{2}}{115200\left(M+C_ {1}\right)^{2}R^{8}N^{16}}\right\}\\ \leq&\exp\left\{A_{1}N\log N-\frac{3m\epsilon}{2048M ^{2}}\right\}+\exp\left\{-\frac{3m\epsilon}{8\left(4M+C_{1}\right)^{2}} \right\}\\ +&\exp\left\{\log 4m+A_{2}N\log N-\frac{n\epsilon^{2}}{A_{ 3}N^{16}}\right\},\end{split}\]

where

\[A_{1} =T_{1}\left(\log\frac{8M\widehat{R}}{C_{1}^{2}}+2\beta\right)+T_{ 2},\] \[A_{2} =T_{1}\left(\log\frac{40M\widehat{R}R^{2}}{C_{1}^{2}}+2\beta+4 \right)+T_{2},\] \[A_{3} =115200\left(M+C_{1}\right)^{2}R^{8}.\]

If we choose the neural network parameter

\[N=\left[A_{4}m^{\frac{1}{2\beta+1}}\right],\]

with \(A_{4}=\left(\min\left\{\frac{3C_{1}^{2}}{2048M^{2}A_{1}},\frac{3C_{1}^{2}}{40 96M^{2}A_{2}}\right\}\right)^{\frac{1}{2\beta+1}}\), and \([x]\) denotes the largest integer that is smaller than or equal to \(x\). Moreover, we choose the second stage sample size

\[n\geq\left\lceil A_{5}m^{\frac{4\beta+17}{2\beta+1}}\right\rceil,\]

with \(A_{5}=\frac{3A_{3}A_{4}^{2\beta+16}}{4096M^{2}C_{1}^{2}}\), and \(\lceil x\rceil\) denotes the smallest integer that is greater than \(x\). Then when \(m\) satisfying the restriction that

\[\log 4m\leq\frac{3C_{1}^{2}A_{4}^{-2\beta}}{4096M^{2}}m^{\frac{1}{2\beta+1}}=A_ {6}m^{\frac{1}{2\beta+1}}, \tag{4.29}\]

we conclude that

\[\begin{split}&\mathit{Prob}\left\{\|\pi_{M}f_{\hat{D},R,N}-f_{ \rho}\|_{\rho}^{2}>9\epsilon\right\}\leq\exp\left\{\frac{3m\epsilon}{4096M^{2} }-\frac{3m\epsilon}{2048M^{2}}\right\}+\exp\left\{-\frac{3m\epsilon}{8\left(4 M+C_{1}\right)^{2}}\right\}\\ &+\exp\left\{\frac{3m\epsilon}{8192M^{2}}+\frac{3m\epsilon}{8192 M^{2}}-\frac{3m\epsilon}{2048M^{2}}\right\}\leq 3\exp\left\{-\frac{3m\epsilon}{2 56\left(4M+C_{1}\right)^{2}}\right\},\end{split}\]
take \(t=9\epsilon\), then when \(t\geq 18C_{1}^{2}N^{-2\beta}\log N\), we have

\[\mathit{Prob}\left\{\left\|\pi_{M}f_{\hat{D},R,N}-f_{\rho}\right\|_{\rho}^{2}>t \right\}\leq 3\exp\left\{-\frac{m^{\frac{2\beta}{2\beta+1}}t}{768\left(4M+C_{1} \right)^{2}}\right\}. \tag{4.30}\]

Then using the property that for the non-negative random variable \(\xi\), \(\mathbb{E}(\xi)=\int_{0}^{\infty}P(\xi>t)dt\), with \(\xi=\left\|\pi_{M}f_{\hat{D},R,N}-f_{\rho}\right\|_{\rho}^{2}=\mathcal{E} \left(\pi_{M}f_{\hat{D},R,N}\right)-\mathcal{E}\left(f_{\rho}\right)\), from (4.30), we yield

\[\mathbb{E}\left\{\mathcal{E}\left(\pi_{M}f_{\hat{D},R,N}\right)- \mathcal{E}\left(f_{\rho}\right)\right\}=\int_{0}^{\infty}Prob\left\{\left\| \pi_{M}f_{\hat{D},R,N}-f_{\rho}\right\|_{\rho}^{2}>t\right\}dt\] \[= \left(\int_{0}^{18C_{1}^{2}N^{-2\beta}\log N}+\int_{18C_{1}^{2}N ^{-2\beta}\log N}\right)Prob\left\{\left\|\pi_{M}f_{\hat{D},R,N}-f_{\rho} \right\|_{\rho}^{2}>t\right\}dt\] \[\leq 18C_{1}^{2}N^{-2\beta}\log N+\int_{18C_{1}^{2}N^{-2\beta}\log N }^{\infty}3\exp\left\{-\frac{m^{\frac{2\beta}{2\beta+1}}t}{768\left(4M+C_{1} \right)^{2}}\right\}\] \[\leq\] \[\leq A_{7}m^{-\frac{2\beta}{2\beta+1}}\log m,\]

where \(A_{7}=18C_{1}^{2}2^{2\beta}A_{4}^{-2\beta}\left(\log A_{4}+\frac{1}{2\beta+1} \right)+2304\left(4M+C_{1}\right)^{2}\). Thus we complete the proof of Theorem 5. 

## Acknowledgments

The work described in this paper is supported partially by the NSFC/RGC Joint Research Scheme [RGC Project No. N-CityU102/20 and NSFC Project No. 12061160462], Germany/Hong Kong Joint Research Scheme [Project No. G-CityU101/20], Laboratory for AI-Powered Financial Technologies, Hong Kong Institute for Data Science, and the CityU Strategic Interdisciplinary Research Grant [Project No. 7020010].

## Appendix

This appendix provides the proof of the continuity of \(f\circ L_{G}^{Q}\) and the bound for the covering number of \(\mathcal{H}_{2}\).

Proof of Proposition 2.: It suffices to prove that \(L_{G}^{Q}(\cdot)\) is continuous. We use again the Kantorovich-Rubinstein distance

\[W_{1}(\mu,\nu)=\sup_{\psi:\left\|\psi\right\|_{C^{0,1}}\leq 1}\left\{\int_{ \Omega}\psi d\mu-\int_{\Omega}\psi d\nu\right\}.\]
Now return to the definition of the functional \(L_{G}^{Q}(\cdot)\) which is defined as

\[L_{G}^{Q}(\mu)=\int_{\Omega}g(Q(x))d\mu.\]

For any polynomial \(Q\) on \(\Omega\), it is easy to see that for any \(x,y\in\Omega\), \(|Q(x)-Q(y)|\leq\|\nabla Q\|_{C(\Omega)}|x-y|\). Then it follows that \(\|g\circ Q\|_{C^{0,1}}\leq|g|_{C^{0,1}}\|\nabla Q\|_{C(\Omega)}\). Combining the above equations yields that

\[\left|L_{G}^{Q}(\mu)-L_{G}^{Q}(\nu)\right|\leq|g|_{C^{0,1}}\|\nabla Q\|_{C( \Omega)}W_{1}(\mu,\nu),\forall\mu,\nu\in\mathcal{P}(\Omega).\]

Holder's inequality shows that, when \(p\leq q\), there holds, \(W_{p}(\mu,\nu)\leq W_{q}(\mu,\nu),\forall\mu,\nu\in\mathcal{P}(\Omega)\). Finally there holds

\[\left|L_{G}^{Q}(\mu)-L_{G}^{Q}(\nu)\right|\leq|g|_{C^{0,1}}\|\nabla Q\|_{C( \Omega)}W_{p}(\mu,\nu),\forall\mu,\nu\in\mathcal{P}(\Omega),\]

which shows the continuity of \(L_{G}(\cdot)\).

Now we give a bound for the covering number \(\mathcal{N}\left(\mathcal{H}_{2},\epsilon,\left\|\cdot\right\|_{C(\Omega)}\right)\)

**Lemma 9**.: _For \(R\geq 1,N\in\mathbb{N}\), and \(\widehat{R}=3R^{4}(10n_{q}+d+18)\), there exists two constants \(T_{1},T_{2}\) depending on \(d,q\) such that_

\[\log\mathcal{N}\left(\mathcal{H}_{2},\epsilon,\left\|\cdot\right\|_{C(\Omega) }\right)\leq T_{1}N\log\frac{\widehat{R}}{\epsilon}+T_{2}N\log N.\]

Proof.: Since \(\mathcal{H}_{2}\) is just the component of a classical two-layer neural network, for the first layer,

\[\left\|\sigma\left(F^{(1)}x-b^{(1)}\right)\right\|_{\infty}\leq\|F^{(1)}\|_{ \infty}\|x\|_{\infty}+\|b^{(1)}\|_{\infty}\leq RN^{2}+R\leq 2RN^{2},\]

for the second layer, we have

\[\left\|H_{2}(x)\right\|_{C(\Omega)}\leq\|F^{(2)}\|_{\infty}\left\|\sigma\left( F^{(1)}x-b^{(1)}\right)\right\|_{\infty}+\|b^{(2)}\|_{\infty}\leq 2R^{2}N^{4}+R \leq 3R^{2}N^{4}.\]

If we choose another function \(\widehat{H}_{2}\) in the hypothesis space \(\mathcal{H}_{2}\) induced by \(\widehat{F}^{(j)}\) and \(\widehat{b}^{(j)}\), satisfying the restriction that,

\[\left|F^{(j)}_{ik}-\widehat{F}^{(j)}_{ik}\right|\leq\epsilon,\ \left\|b^{(j)}- \widehat{b}^{(j)}\right\|_{\infty}\leq\epsilon,\]

then by the Lipschitz property of ReLU,

\[\left\|\sigma\left(F^{(1)}x-b^{(1)}\right)-\sigma\left(\widehat{ F}^{(1)}x-\widehat{b}^{(1)}\right)\right\|_{\infty}\leq \left\|F^{(1)}-\widehat{F}^{(1)}\right\|_{\infty}\left\|x\right\|_{ \infty}+\left\|b^{(1)}-\widehat{b}^{(1)}\right\|_{\infty}\] \[\leq (d+1)\epsilon.\]
It follows that

\[\left\|H_{2}(x)-\widehat{H}_{2}(x)\right\|_{C(\Omega)}\leq \left\|F^{(2)}\left(\sigma\left(F^{(1)}x-b^{(1)}\right)-\sigma \left(\widehat{F}^{(1)}x-\widehat{b}^{(1)}\right)\right)\right\|_{\infty}\] \[+ \left\|\left(F^{(2)}-\widehat{F}^{(2)}\right)\sigma\left(\widehat{ F}^{(1)}x-\widehat{b}^{(1)}\right)\right\|_{\infty}+\left\|b^{(2)}-\widehat{b}^{(2)} \right\|_{\infty}\] \[\leq RN^{2}(d+1)\epsilon+2RN^{2}n_{q}(2N+3)\epsilon+\epsilon\leq(10n _{q}+d+2)RN^{3}\epsilon=:\tilde{\epsilon},\]

Therefore, by taking a \(\epsilon\)-net of each free parameter in \(F^{(j)}\) and \(b^{(j)}\), the \(\widetilde{\epsilon}\)-covering number of \(\mathcal{H}_{2}\) can be bounded by

\[\mathcal{N}\left(\mathcal{H}_{2},\widetilde{\epsilon},\left\|\cdot\right\|_{ C(\Omega)}\right)\leq\left\lceil\frac{2RN^{2}}{\epsilon}\right\rceil^{n_{q}d+ qn_{q}+q(2N+3)}\left\lceil\frac{2R}{\epsilon}\right\rceil^{4N+6}\leq\left( \frac{\widehat{R}}{\widetilde{\epsilon}}\right)^{T_{1}N}N^{T_{2}N}, \tag{31}\]

where \(\widehat{R},T_{1}\) and \(T_{2}\) are defined in Lemma 5. Thus we finish the proof by taking the logarithm.
