Improved Convergence Rate of Nested Simulation with LSE on Sieve]Improved Convergence Rate of Nested Simulation with LSE on Sieve [
Corresponding author<EMAIL_ADDRESS>
1]Ruoxue Liu[Data Science and Analytics Thrust, Information Hub, The Hong Kong
University of Science and Technology
(Guangzhou), 1 Duxue Rd, Nansha District, Guangzhou, Guangdong, 51453, China.  ]

2]Liang Ding[Schoold of Data Science, Fudan University, 220 Handan Rd, Yangpu District, Shanghai, 200437, China. ]

3]Wenjia Wangfoot:hkust

4]Lu  Zoufoot:hkust

Nested simulation encompasses the estimation of functionals linked to conditional expectations through simulation techniques. In this paper, we treat conditional expectation as a function of the multidimensional conditioning variable and provide asymptotic analyses of general Least Squared Estimators  on sieve, without imposing specific assumptions on the function's form. Our study explores scenarios in which the convergence rate surpasses that of the standard Monte Carlo method and the one recently proposed based on  kernel ridge regression. We also delve into the conditions that allow for achieving the best possible square root convergence rate among all methods. Numerical experiments are conducted to support our statements.

[Mathematics Subject Classification]62G08, 62G20, 90-10

§ INTRODUCTION

Simulation often aims to estimate a function linked to conditional expectations, which can take the form of either an expectation or a quantile. Typically, there is no analytical expression readily available for either the conditional expectation or the function itself. Therefore, the estimation process involves two simulation layers.  First, in the “outer” level, we simulate the random variable being conditioned on. Then, in the “inner” level, we simulate another random entity, conditioning on each previously generated sample (referred to as a “scenario”). This problem category is known as nested simulation. In this paper, our primary focus centers on using simulation to estimate the following quantity:

\begin{align}\label{eqgoal}
\theta = \mathcal{T}(\mathbb{E}[Y \vert  X]),
\end{align}

where $X$ is a random variable in $\mathbb{R}^d$, $Y$ is a real-valued random variable, and $\mathcal{T}$ is a functional mapping a probability distribution to a real number.
Nested simulation is closely related to conditional Monte Carlo, a versatile technique for reducing variance in simulations <cit.>. Any expectation, denoted as $\mathbb{E}[Y]$, can be expressed as $\mathbb{E}[\mathbb{E}[Y \vert X]]$. By the law of total variance, the variance of $\mathbb{E}[Y \vert  X]$ is less than or equal to the variance of $Y$. Consequently, $\mathbb{E}[\mathbb{E}[Y \vert X]]$ serves as an unbiased estimator with lower variance than $Y$. It's especially valuable when $X$ and $Y$ are strongly correlated, and when the conditional expectation can be computed or estimated efficiently. Additionally, conditional Monte Carlo can be used as a smoothing technique for gradient estimation [2, 3].

An important research topic in nested simulation is the allocation of simulation resources, specifically the determination of outer-level scenarios and inner-level samples. In this study, we use equal inner-level samples per outer-level scenario. For a given budget $\Gamma$, standard nested simulation has outer-level sample size $n$ growing as $\Gamma^{2/3}$ and inner-level sample size $m$ for each out-level sample as $\Gamma^{1/3}$. It has been shown that the standard nested simulation achieves a convergence rate at $\Gamma^{-1/3}$ [4, 5], which is slower than the optimal rate $\Gamma^{-1/2}$ over all Monte Carlo methods.

Recent studies [6, 7] have shed light on an opportunity to improve the optimal convergence rate from $\Gamma^{-1/3}$ to faster convergence rates. In particular, if the conditional expectation function $f(\BFx) = \mathbb{E}[Y\vert X=\BFx]$ is a Sobolev function with respect to $\BFx$, the convergence rate $\Gamma^{-1/2}$ can be attained in certain scenarios [7]. This condition allows for the exchange of information among different outer-level scenarios ${\BFx}$ and can be effectively estimated using kernel ridge regression (KRR) [8]. In this study, we extend the settings in [7] to more general function spaces and estimators. Specifically, we unify the traditional nested simulation and the recent KRR method under the so-called least squared estimation (LSE) on sieve. Within this framework, we have the flexibility to incorporate various advanced estimation techniques, such as KRR with inducing points and neural networks, to enhance the overall effectiveness of nested simulation. We then investigate the circumstances under which the best possible convergence rate $\Gamma^{-1/2}$ is attainable.

§.§ Notation and Organization

Throughout the paper, we use the following notation. For independent and identically distributed (i.i.d.) outer scenario $\{\BFx_i\}_{i=1}^n$, let $P_n=\frac{1}{n}\sum_{i=1}^n\delta_{\BFx_i}$ be the empirical density of $X$ with $\delta_{\BFx_i}$ a point mass at $\BFx_i$. For any functions $f$, $g$ defined on $\{\BFx_i\}_{i=1}^n$, define
the empirical inner product of $f$ and $g$ as $\langle f, g\rangle_n=\int f(\BFx)g(\BFx)dP_n$ and the empirical $L_2$ norm $\|f\|_n^2=\langle f,f\rangle_n.$ We also define the $L_2$ norm induced by the distribution of $X$, denoted by $P_X$, as $\|f\|_{L_2(X)}=\int f^2dP_X$. Let $\Omega$ be the support of the distribution of $X$. Then the $L_\infty$ norm of a function $f$ is defined as $\|f\|_{L_\infty}=\sup_{\BFx\in\Omega}|f(\BFx)|$. For any integer $K$, let $[K]$ denote the set $\{1,\cdots,K\}$

For two positive sequences $a_n$ and $b_n$, $a_n = \CalO_{\pr}(b_n)$ means  that for any $\varepsilon>0$, there exists $C>0$ such that
$\pr(a_n >C b_n) <\varepsilon$
for all $n$ large enough.  We also write
$a_n=\CalO(b_n)$ and $a_n\asymp b_n$
if there exist some constants $C,C^\prime >0$ such that
$a_n\leq C b_n$ and
$C^\prime \leq a_n/b_n \leq C $, respectively, for all $n$ large enough.

The remainder of the paper is organized as follows.
In Section <ref>, we introduce the background of nested simulation and formulate the research question.
In Section <ref>, we study the general framework LSE on sieve.
In Section <ref>, we analyze its asymptotic properties.
In Section <ref>, we use KRR with inducing points and neural network as examples.
In Section <ref>, we conduct numerical experiments to verify our theorems.
In Section <ref>, we conclude the paper.
Technical results and proofs are presented in the Appendix.

§ PROBLEM FORMULATION

In this paper, we consider the following forms of functional $\CalT$ in (<ref>):

* $\mathcal{T}(Z) = \E[\eta(Z)]$ for some function $\eta:\Real^d\mapsto \Real$;

* $\mathcal{T}(Z) =  \inf\{z\in\Real:\pr(Z\leq z)\geq \tau\}$, i.e., the VaR of $Z$ at level $\tau\in(0,1)$.

Functionals in the Form <ref> are known as nested expectations and are employed in supervised learning, including stochastic kriging and classification. On the other hand, functionals represented by Form <ref> are linked to risk measures and find applications in risk management.

The nested simulation usually involves the following three steps to generate data.
(i) Generate $n$ i.i.d.  outer-level scenarios $\{\BFx_i:i=1,\ldots,n\}$.
(ii) For each $\BFx_i$, generate $m$ i.i.d. inner-level samples from the conditional distribution of $Y$ given $X=\BFx_i$, denoted by $\{y_{ij}:j=1,\ldots,m\}$.
(iii) An estimator $\hat{f}(X)$ of $f(X)=\E[Y\vert X]$, adapted to the samples $(X, Y)$, is constructed to estimate the  target $\theta$.

Let $\hat{f}_{(1)} \leq \cdots \leq \hat{f}_{(n)}$  denote the order statistics of $\{\hat{f}(\BFx_1),\ldots, \hat{f}(\BFx_n)\}$. Then, $\hat{f}_{(\lceil \tau n\rceil)}$ is the sample quantile at level $\tau$,
where $\lceil z\rceil$ denotes the least integer greater than or equal to $z$.
Nested simulation then estimates $\theta$ via

\begin{equation}\label{eq:standard-estimator}
\hat\theta_{n,m} =
\left\{
\begin{array}{ll}
n^{-1}\sum_{i=1}^n \eta(\hat{f}(\BFx_i) ),     &\quad\mbox{if } \mathcal{T}(\cdot)= \E[\eta(\cdot)],  \\[0.5ex]
\hat{f}_{(\lceil \tau n\rceil)},     &\quad  \mbox{if }\mathcal{T}(\cdot) = \mbox{quantile}.
\end{array}
\right.
\end{equation}

A central problem in nested simulation is budget allocation. Given a simulation budget $\Gamma$ and an estimator $\hat\theta_{n,m}$, the optimal values for $n$ and $m$ to minimize the error in estimating $\theta$ can be determined by analyzing the asymptotic behavior of the statistics $|\hat\theta_{n,m} - \theta|$. [4] and [5] demonstrate that, without any assumption on $f$, the simulation budget should be allocated in such a way that $n\asymp \Gamma^{2/3}$ and $m\asymp \Gamma^{1/3}$ as $\Gamma\to\infty$, resulting in

\[
\left(\mathbb{E}\left[(\hat\theta_{n,m}-\theta)^2\right]\right)^{1/2} \asymp \Gamma^{-1/3}.
\]

[7] improves upon these results by imposing a smoothness condition on $f$ and shows that, for $f$ smooth enough, the simulation budget can be reduced to $n\asymp \Gamma$ and $m\asymp 1$ as $\Gamma\to\infty$, leading to

\[
|\hat\theta_{n,m}-\theta| = \mathcal{O}_p(\Gamma^{-1/2}).
\]

§ LSE ON SIEVE

In standard nested simulations, the estimator $\hat{f}(\BFx_i)$ can be calculated as the sample averages by $\bar{y}(\BFx_i) = \frac{1}{m}\sum_{j=1}^m y_{ij}$. Alternatively, it can be expressed as the following least squared estimator:

\begin{equation}\label{eq:sample_avg}
\hat{f} = \argmin_{h\in\mathcal{F}_n}\|h-\bar{y}\|_n^2
\end{equation}

where $\mathcal{F}_n$ is defined as $\{\ind_{\{\BFx_i=\BFx\}}:i=1,\cdots,n\}$. On the other hand, [7] constructs $\hat{f}$ using KRR, which can be seen as the following least squared estimator:

\begin{equation}\label{eq:KRR}
\hat{f} = \argmin_{h\in\mathcal{H}(\delta_n)}\|h-\bar{y}\|_n^2
\end{equation}

where $\mathcal{H}(\delta)$ is defined as $\{h\in\mathcal{H}: \|h\|_{\mathcal{H}}\leq \delta\}$, $\|\cdot\|_{\mathcal{H}}$ denotes the norm of a Sobolev space $\mathcal{H}$, and $\delta_n>0$ is adaptive to sample sizes $n,m$; see  9, 10 for details.

Both (<ref>) and (<ref>) are evidently instances of LSE on sieves. To see this, we begin by introducing the definition of LSE on sieve. Suppose the true function $\E[Y|X=\BFx]$ resides in a function space $\CalF$. Consider another space of functions, denoted as $\mathcal{F}_{n,m}$, which is adaptable to the available data. This space doesn't necessarily need to be a subspace of $\mathcal{F}$, but it should not be more complex than $\mathcal{F}$. Define

\begin{equation}\label{eq:sieve}
\hat{f}_{n,m}=\argmin_{h\in\mathcal{F}_{n,m}}\|h-\bar{y}\|_n^2.
\end{equation}

We say that $\CalF_{n,m}$ is a sieve and $\hat{f}_{n,m}$ is an LSE on sieve. The sieve $\CalF_{n,m}$ can be a function space, a vector space, or a neural network accordingly.

Similar to [7], we impose the following assumptions on the data.

The support of $X$ is a bounded convex set  $\Omega\subset \Real^d$.
Moreover, $X$ has a probability density function that is bounded above and below away from zero.

$y_{ij}$ conditioned on $\BFx_i$ is of the form

$$y_{ij} = f({\BFx}_i) + \varepsilon_{ij}$$

where $\varepsilon_{ij}$ represents independent zero-mean  random variables, which may not have identical distributions but satisfy the following $\sigma^2$-sub-Gaussian assumption:

\[\max_{i,j} \E e^{t\varepsilon_{ij}}\leq e^{\sigma^2t^2/2},\quad \forall t\in \Real.\]

In contrast to [7], which relies on specific assumptions about the underlying function $f$ and employs the KRR technique, our study adopts a more generalized framework. We investigate the necessary conditions governing the relationship between the “complexity” of $\mathcal{F}_{n,m}$ and the underlying function $f$, where the complexity of a space is quantified by the following measures.

Let $\mathcal{M}$ be a space equipped with metric $\|\cdot\|$. The $u$-covering number $N(u,\mathcal{M},\|\cdot\|)$ is defined as the minimum number of balls with radius $u$ necessary to cover $\mathcal{M}$. The $u$-entropy of $\mathcal{M}$ is then $H(u,\mathcal{M},\|\cdot\|)=\log N(u,\mathcal{M},\|\cdot\|) \vee 0$.

Following the above definition, we impose an additional assumption regarding the “regularity”  of  the spaces $\CalF_{n,m}$ and $\CalF$. This assumption is mild and can be satisfied by both  $\CalF_n$ in (<ref>) and $\CalH(\delta_n)$ in (<ref>).

(i) $H(u,\CalF_{n,m},\|\cdot\|_n)\leq H(u,\CalF,\|\cdot\|_n)<\infty$ for any $u>0$; (ii) there exists a constant $b<\infty$ such that $\sup_{h\in\CalF\cup\CalF_{n,m}}\|h\|_{L_\infty}<\infty$; (iii) $\CalF$ and $\CalF_{n,m}$ are Donsker with respect to the distribution of $X$ (see <cit.> for details).

Given the entropy of $\CalF_{n,m}$, we are ready to present our first theorem regarding the mean squared error (MSE) of LSE on sieve.

\[f^*_{n,m}=\argmin_{h\in\CalF_{n,m}}\|h-f\|_n,\quad\text{and}\quad \hat{f}_{n,m}=\argmin_{h\in\CalF_{n,m}}\|h-\bar{y}\|_n.\]

Define for $0<{\delta}<2^7\sigma/\sqrt{m}$,

\begin{align*}
&\CalF_n({\delta})=\{h\in\CalF_{n,m}: \|h-f^*_{n,m}\|_n\leq {\delta}\},\\
&  \mathcal{J}({\delta})=\int_{\frac{\delta^2\sqrt{m}}{2^7{\sigma}}}^{{\delta}}\sqrt{H(u,\CalF_n({\delta}),\|\cdot\|_n)}du \vee {\delta}.
\end{align*}

Suppose Assumptions <ref>– <ref> hold. Take any upper bound $\Psi({\delta})\geq \mathcal{J}({\delta})$ in  such a way that $\Psi({\delta})/{\delta}^2$ is a non-decreasing function for any $0<{\delta}<2^7{\sigma}/\sqrt{m}$ . Then for a constant $C$ only depending on $\sigma$, and for

\begin{equation}\label{eq:thm1_condition}
\sqrt{nm}\delta^2_{n,m}\geq c\Psi(\delta_{n,m}),\quad \forall n,m,
\end{equation}

we have

\begin{align}\label{eq:LSE_convergence}
&\pr\left( \|\hat{f}_{n,m}-f^*_{n,m}\|_n\geq \delta_{n,m}\right) \leq {\pr}\left(\|\bar{\varepsilon}\|_n\geq \frac{\sigma}{\sqrt{m}}\right)+C\exp\left[-\frac{nm\delta_{n,m}^2}{C^2}\right],\nonumber\\
&\|\hat{f}_{n,m}-f^*_{n,m}\|_n\leq \CalO_{\pr}(\delta_{n,m}+\|f^*_{n,m}-f\|_n).
\end{align}

Theorem <ref> considers both the sample sizes $n$ of outer scenario $\{\BFx_i\}$ and $m$ of the inner scenarios $\{y_{ij}\}$. From this theorem, the MSE of the LSE can be decomposed into two parts: the estimation error $\mathcal{O}_p(\delta_{n,m})$ and the approximation error $\|f_{n,m}^*-f\|_n$. In simple terms, the estimation error increases with the complexity of $\mathcal{F}_{n,m}$, while the approximation error tends to decrease when $\mathcal{F}_{n,m}$ effectively approximates the underlying function $f$.

§ ASYMPTOTIC ANALYSIS

If we establish the asymptotic behavior of $|\hat{\theta}_{n,m} - \theta|$ concerning the sample sizes $n$ and $m$, we can then determine the optimal allocation of the budget $\Gamma=nm$.

For the nested expectation functional $\CalT$, we have the following theorem:

Let $\mathcal{T}(\cdot) = \E[\eta(\cdot)]$ and $\eta$ be a twice-differentiable function with bounded first- and second-order derivatives.
Suppose Assumptions <ref>–<ref> hold. Moreover, $m$ satisfies $\sqrt{m}(\delta_{n,m}+\|f_{n,m}^*-f\|_n)=\CalO(1)$.
\[|\hat \theta_{n,m}-\theta| = O_{\pr}\left(n^{-\frac{1}{2}}+\delta_{n,m}^2+\|f^*_{n,m}-f\|_n+\hat{\CalL}\right),\]

where $\delta_{n,m}$ and $f^*_{n,m}$ are as defined in Theorem <ref> and

\begin{align*}
\hat{\CalL}=\left|\E\left[\eta'(f(X))\right]\E_{X,\{\varepsilon_i\}} [f^*(X)-\hat{f}(X)]\right|.
%&\Delta_n=\inf\{\delta: \sqrt{n}\delta^2\geq C_1\int^\delta_0\sqrt{H(u,\CalF_n(\delta),\|\cdot\|_n)}du\vee \delta\}.
\end{align*}

It is worth noticing that $\hat{\CalL}$ depends on the linearity $\hat{f}_{n,m}$. For linear estimators, i.e., $\hat{f}_{n,m}=\sum_{i=1}^n\phi_i(\cdot;X)\bar{y}_i$,  $\hat{\CalL}=0$ because

\begin{align*}
&\E_{X,\{\varepsilon_i\}} [f^*(X)-\hat{f}(X)]=\E_{X,\{\varepsilon_i\}}\sum_{i=1}^n\phi_i \left(f(\BFx_i)+\varepsilon_i\right)-\phi_if(\BFx_i)=0.
\end{align*}

Linear estimators include many commonly used estimators, such as KRR and spline regression. For some non-linear estimator, such as  wavelet-threshold estimators [12], it can be bounded by the approximation error $\|f^*_{n,m}-f_{n,m}\|_{L_2}$. For highly non-linear estimators such as neural networks, the analysis of $\hat{\CalL}$ becomes more involved.

If $\CalT$ is the risk measure, we need to impose two more assumptions.

(i) $\CalF$ and $\CalF_{n,m}$ are star-shaped, i.e., for any $h$ in $\CalF$ or $\CalF_{n,m}$, $ah$ is also in   $\CalF$ or $\CalF_{n,m}$ for all $0\leq a \leq 1$. (ii) for any $\delta >0$ and $h\in\CalF\cup\CalF_{n,m}$ with $\|h\|_{L_2(X)}\leq \delta$,  there exists universal constant $\alpha \in(0,1]$ such that $\|h\|_{L_\infty}\leq \CalO(\delta^\alpha)$.

There exist positive constants $C_1$, $C_2$, $t_0$, and $\beta\leq \gamma$ such that

\begin{align*}
C_2 t^{\gamma} \leq \pr(|f(X)-z|\leq t)\leq C_1 t^{\beta},\quad
\forall t\in(0,t_0],\;\forall z\in \{f(\BFx):\BFx\in\Omega \}.
\end{align*}

Assumption <ref> is for handling the error between the order statistics $\hat{f}_{(\lceil \tau n\rceil)}$ and the true quantile of $f(X)$, where we require point-wise asymptotic property. It is worth noting that this assumption is  mild and can be met when $\mathcal{F}$ is chosen from commonly encountered spaces. For example, when $\CalF$ is a order-$m$ Sobolev space, we have $\alpha=1-1/(2m)$ <cit.>, <cit.>; when $\CalF$ is the reproducing kernel Hilbert space (RKHS) generated by  Gaussian kernel, we have $\alpha=1$ <cit.>; when $\CalF$ is a  ReLU neural network, we have $\alpha=d/(d+2)$ as we will show in Section <ref>. These relationships are often referred as function space embeddings in general.

Assumption <ref> is Assumption 5 of [7]. For more discussion of Assumption <ref>, we refer to Remark 4 of [7]. In particular, Assumption <ref> implicitly implies $\beta \leq 1 \leq \gamma $.

Equipped with theses tools, we are now ready to present the following theorem.

Let $\mathcal{T}(\cdot) = \inf\{z\in\Real:\pr(\cdot \leq z)\geq \tau\}$ for some $\tau\in(0,1)$. Suppose Assumptions <ref>–<ref> hold. Then

\[|\hat \theta_{n,m}-\theta|=\CalO_{\pr}\left(\left|\frac{\Delta_n}{\sqrt{m}}+\|f^*_{n,m}-f\|_n+\delta_{n,m}\right|^{\kappa}+n^{-\frac{1}{2\gamma}}\right), \]

where $\kappa=\alpha\beta/\gamma$, $\delta_{n,m}$ and $f^*_{n,m}$ are as defined in Theorem <ref> and

\begin{equation*}
\Delta_n=\inf\left\{\delta: \frac{\delta^2}{b}\geq \inf_{\varepsilon\geq 0}\left\{4\varepsilon+\frac{C}{\sqrt{n}}\int_\varepsilon^{\delta} \sqrt{H(u,\CalH\cap B(\delta,\|\cdot\|_n),\|\cdot\|_n)}du\right\}\right\}.
\end{equation*}

§ EXAMPLES

In this section, we show some examples to illustrate our theorems. These examples make it clear that the Least Squared Estimation (LSE) on sieve offers a higher degree of flexibility and, in certain specific scenarios, outperforms the KRR estimator.

§.§ Kernel Ridge Regression

In this subsection, we show that KRR can be formulated as a LSE on a sieve that is independent of the data. We further illustrate that by selecting a sieve that adapts to the data, we can obtain an LSE on sieve equivalent to KRR with inducing points, and this estimator exhibits an improved convergence rate.

§.§.§ Nested Expectation

We first consider the following KRR in constrained form with sieve $\CalF_{n,m}$  a fixed subset of the RKHS  generated by the Matérn-$\nu$ kernel. This space is  norm-equivalent to the $s=\nu+d/2$-order Sobolev space  $\mathscr{H}^s(\Omega)$ equipped with norm $\|\cdot\|_{\mathscr{H}^s(\Omega)}$. Under our framework, we can analyze the asymptotic properties of KRR estimator regardless the specific form of Matérn kernel. More precisely, the KRR problem is equivalent to the following constrained optimization problem for $\nu>d/2$:

\begin{equation}\label{eq:constrained_KRR}
\min_{h}\|h-\bar{y}\|_n^2,\quad\text{s.t.}\ \|h\|_{\mathscr{H}^s(\Omega)}^2\leq \CalR
\end{equation}

where  $\CalR\geq \|f\|^2_{\mathscr{H}^s(\Omega)}$ is independent of $n,m$. In this case, the approximation error $\|f^*-f\|_n=0$ because $f\in\CalF_{n,m}=\{h:\|h\|_{\mathscr{H}^s(\Omega)}^2\leq \CalR\}$. Furthermore, according to [15], the $u$-entropy of $\mathscr{H}^s(\Omega)$ is $u^{-\frac{d}{s}}$, and according to <cit.>, $\CalF=\CalF_{n,m}$ is Donsker. Therefore, it is straightforward to check that $\CalF=\CalF_{n,m}$ satisfy Assumptions <ref> and <ref>. Applying Theorem <ref> yields

\[\delta_{n,m}=\inf\{\delta>0: \sqrt{nm}\delta^2\geq c \delta ^{1-\frac{d}{2m}}\}=\CalO\left(|nm|^{-\frac{s}{2s+d}}\right).\]

When $\eta$ has bounded twice derivative, we cau use Theorem <ref> to derive that

\[|\hat \theta_{n,m}-\theta|=|\frac{1}{n}
\sum_{i=1}^n\hat{f}\circ\eta(\BFx_i)-\E[\eta\circ f(X)]| = \CalO_{\pr}\left(n^{-\frac{1}{2}}+|nm|^{-\frac{2s}{2s+d}}\right).\]

By setting $m=1$ and $n=\Gamma$, it is straightforward to derive that we can have the improved convergence rate $|\hat \theta_{n,m}-\theta|=\CalO(\Gamma^{-\frac{1}{2}})$. This result is exactly the same as the improved convergence rate in  [7] for $\nu>d/2$.

However, for $\nu\leq d/2$, [7] has sub-optimal result. We now show that by adaptively choosing the sieve $\CalF_{n,m}$, we can still have the optimal convergence rate $\Gamma^{-\frac{1}{2}}$. Without loss of generality, suppose $X$ is uniformly distributed on $\Omega$. We then select $\CalS_{n}=\CalO(n^{\frac{d}{2s+d}})$ points from $\{\BFx_i\}$, denoted as $\{\tilde{\BFx}_i\}$ such that the filled distance of $\tilde{\BFx}_j$ satisfy

\[h_{\Omega,\tilde{\BFX}}=\sup_{\BFx\in\Omega}\inf_{j\in[\CalS_n]}\|\tilde{\BFx}_j-\BFx\|=\CalO(\CalS_n^{-\frac{1}{d}}).\]

Let $\CalF_{n,m}=\text{span}\{k(\tilde{\BFx}_i,\cdot)\}$ where $k$ is the Matérn-$\nu$ kernel. According to [16],

\[\inf_{h\in\CalF_{n,m}}\|h-f\|_n=\CalO_{\pr}(\CalS_n^{-s/d})=\CalO_{\pr}(n^{-\frac{s}{2s+d}}).\]

Because $\CalF_{n,m}$ is a finite dimensional space, we can use <cit.> to derive that

\[H(u,\CalF_n(\delta),\|\cdot\|_n)\leq Cn^{\frac{d}{2s+d}}\log(\frac{4\delta+u}{u}).\]

From straightforward calculations, it can be shown that $\delta_{n,m}=\CalO(n^{-\frac{s}{2s+d}}$). By setting $n=\Gamma$ and $m=1$, together with the fact that $s=\nu+d/2> d/2$, we have the same convergence rate $|\hat \theta_{n,m}-\theta| =\CalO_{\pr}(\Gamma^{-\frac{1}{2}})$.

The above LSE on sieve is also known as KRR with inducing points. We can notice that KRR with inducing points has two advantages compared with KRR. Firstly, it can achieve the optimal convergence rate for all choices of $\nu>0$; secondly, it has lower  computational complexity – the computational time complexity of KRR is $\CalO(n^3)$ while the one of ours is $\CalO(n\CalS_n^2)=\CalO(n^{\frac{2s+3d}{2s+d}})$. For analysis of computational time complexity of KRR with inducing points, please refer to [17, 18] for details.

§.§.§ Risk Measures

In [7], it has been shown that the convergence rate of $|\hat{f}_{(\lceil\frac{\tau}{n}\rceil)}-\text{VaR}_\tau(f(X))|$ cannot achieve the optimal rate $\Gamma^{-\frac{1}{2}}$  for any Matérn-$\nu$ kernel with $\nu<\infty$. In this study, we discuss the case that the underlying $f$ resides in the RKHS $\CalF=\CalH_k$ induced by the Gaussian kernel $k(\BFx,\BFx')=\exp(-\omega\|\BFx-\BFx'\|^2)$, which can be treated as the limit $\nu\to\infty$. Similar to Section <ref>,  we rewrite the KRR in the constrained optimization form (<ref>) with sieve $\CalF_{n,m}=\{h: \|h\|_{\CalH_k}^2\leq\CalR\}$ where $\|\cdot\|_{\CalH_k}$ is the RKHS norm induced by the Gaussian kernel. According to [19], the $u$-entropy of $\CalF_{n,m}$ can be bounded as follows:

\[H(u,\CalF_{n,m},\|\cdot\|_n)=\CalO(\frac{(\log \frac{1}{u})^{d+1}}{(\log \log  \frac{1}{u})^d}).\]

Direct calculations show that the $\delta_{n,m}$ in Theorem <ref> and $\Delta_n$ in Theorem <ref> are as follows:

\begin{align*}
&\delta_{n,m}=\CalO\left((nm)^{-\frac{1}{2}}[\log nm]^{\frac{d+1}{2}}\right),\quad
\Delta_n= \CalO\left(n^{-\frac{1}{2}}[\log n]^{\frac{d+1}{2}}\right).
\end{align*}

Then, we can use <cit.> to get  embedding constant $\alpha=1$ for the Gaussian RKHS:

\[\|f-\hat{f}\|_{L_\infty}\leq\CalO(\|f-\hat{f}\|_{L_2}),\quad f-\hat{f}\in\CalF_{n,m}.\]

Therefore, if the parameters $\beta=\gamma=1$, we can have the almost optimal convergence rate by setting $m=1$:

\[|\hat{f}_{(\lceil\frac{\tau}{n}\rceil)}-\text{VaR}_\tau(f(X))|\leq\CalO_{\pr}(\Gamma^{-\frac{1}{2}}[\log\Gamma]^{\frac{d+1}{2}}).\]

To lower the computational time complexity, we choose $\CalS_n=[\log n]^{\frac{d}{2}}$ inducing points $\tilde{\BFX}=\{\tilde{\BFx}_j\}$ from $\{\BFx_i\}$ such that the filled distance of $\tilde{\BFx}_j$ satisfy

\[h_{\Omega,\tilde{\BFX}}=\sup_{\BFx\in\Omega}\inf_{j\in[\CalS_n]}\|\tilde{\BFx}_j-\BFx\|=\CalO(\CalS_n^{-\frac{1}{d}})=\CalO([\log n]^{-\frac{1}{2}}).\]

Let $\CalF_{n,m}=\text{span}\{k(\tilde{\BFx}_j,\cdot)\}$ where $k$ is the Gaussian kernel. According to <cit.>,

\[\inf_{h\in\CalF_{n,m}}\|h-f\|_n\leq\CalO(\exp[-1/h_{\Omega,\tilde{\BFX}}])=\CalO(n^{-\frac{1}{2}}).\]

Because $\CalF_{n,m}$ is a finite dimensional space, we can use <cit.> to derive that

\[H(u,\CalF_n(\delta),\|\cdot\|_n)=[\log n]^{\frac{d}{2}}\log(\frac{4\delta+u}{u}).\]

Similarly, it can be shown that $\delta_{n,m}=\CalO(n^{-\frac{1}{2}}[\log n]^{\frac{d+1}{2}}$). By setting $n=\Gamma$ and $m=1$,  we have the same convergence rate $|\hat \theta_{n,m}-\theta| =\CalO_{\pr}(\Gamma^{-\frac{1}{2}}[\log \Gamma]^{\frac{d+1}{2}})$. Compared with the $\CalO(n^3)$ computational time complexity of KRR, the LSE on sieve only requires $\CalO(n\CalS_n^2)=\CalO(n[\log n]^{d+1})$ time to achieve the same convergence rate.

§.§ Deep Neural Networks

We now apply our theories to obtain bounds for a specific type of neural networks called feed-forward multi-layer networks with ReLU activation. Here, we consider Hölder-$s$  underlying $f$ with smoothness parameter $s> d/2$, meaning that there exists some constant $M$ such that

\[\|f\|_{C^{s}(\Omega)}=\max_{|\BFv|<\lfloor s\rfloor}\|D^{\BFv}f\|_{L_\infty}+\max_{|\BFv|=\lfloor s\rfloor}\sup_{\BFx,\BFx'\in\Omega}\frac{|D^{\BFv}f(\BFx)-D^{\BFv}f(\BFx')|}{\|\BFx-\BFx'\|^{s-\lfloor s\rfloor}}\leq L,\]

where $\BFv=(v_1,\cdots,v_d)\in\NatInt^d$ consists of $d$ integers and

\[D^{\BFv}f=\frac{\partial^{|\BFv|}f}{\partial x_1^{v_!}\cdots\partial x_d^{v_d}}.\]

Denote $A(x)=\max(x,0)$ as the ReLU activation function. An $H$-layer ReLU net $\Phi(H,W,S,B)$ is a parametrized function space consists of functions of the form:

Here, $A(\cdot)$ is entry-wise positive part operation, $\BFW^{(1)}\in\Real^{W\times d}$,   $\BFW^{(h)}\in\Real^{W\times W}$ and $\BFb^{(h)}\in\Real^{W}$ for $h\in[H]$ and they satisfy

\begin{align*}
&\sum_{h=1}^H\|\BFW^{(h)}\|_0+\|\BFb^{(h)}\|_0\leq S,\quad \max\{\max_{i,j,h}|\BFW^{(h)}_{i,j}|,\max_{i,h}|\BFb^{(h)}_i|\}\leq B,\\
\end{align*}

where $\|\BFM\|_{0}$ counts the number of non-zero elements of matrix $\BFM$. Intuitively, $H$ is the number of hidden layers, $W$ is the maximum size of weight matrices, $S$ is the sparsity, and $B$ is the maximum values allowed for the outputs of all activations. In summary, the target space $\CalF$ consists of all Hölder-$s$ functions
and the sieve $\CalF_{n,m}$ is a ReLU neureal network:

\[\CalF=\{f: \|f\|_{C^s(\Omega)}\leq L\},\quad \CalF_{n,m}=\Phi(H,W,S,B).\]

According to [21], $\Phi(H,W,S,B)$ is a subset of Hölder-1 space. Because it is also a parametric class, we can use the same reasoning in <cit.> to show that $\Phi(H,W,S,B)$ is Donsker. Moreover, we  apply  Gagliardo–Nirenberg interpolation inequality on Hölder-1 class, which leads to the following inequality for any function $f\in\CalF_{n,m}$:

\begin{equation}\label{eq:GNS_ineql}
\|f\|_{L_\infty}\leq C_1\|f\|_{L_2}^{\frac{d}{d+2}}\|D^{\BFv}f\|_{L_\infty}^{\frac{2}{d+2}}\leq C_2\|f\|_{C^1},
\end{equation}

where $C_1$ and $C_2$ are some universal constants. In summary,  $\Phi(H,W,S,B)$ satisfies Assumptions <ref> and <ref>. For target function space $\CalF$, it is well known that Hölder-$s$ space can be embedded on Sobolev space $\mathscr{H}^s(\Omega)$ for compact $\Omega$. Therefore, $\CalF=C^s(\Omega)$ also satisfies Assumptions <ref> and <ref> for $s>d/2$.

By <cit.>, the $u$-entropy of $\CalF_{n,m}$ can be bounded as follows:

\begin{equation*}
H(u,\Phi(H,W,S,B),\|\cdot\|_n)\leq 2SH\log [(B\vee 1)(W+1)] +S\log\frac{H}{u}.
\end{equation*}

Therefore, it can be shown that the $\delta_{n,m}$ in Theorem <ref> and the $\Delta_n$ in Theorem <ref> can be derived by

\begin{align*}
&\delta_{n,m}=\inf\{0<{\delta}<\frac{2^7\sigma}{\sqrt{m}}: \sqrt{nm}\delta^2\geq c\left( \sqrt{S}\int ^{\delta}_{\frac{\delta^2\sqrt{m}}{2^7\sigma}}\sqrt{\log\frac{H}{u} }du\vee \delta\right)\},\\
&\Delta_n=\inf\left\{\delta: \frac{\delta^2}{b}\geq \inf_{\varepsilon\geq 0}\left\{4\varepsilon+\frac{C\sqrt{S}}{\sqrt{n}}\int_\varepsilon^{\delta} \sqrt{\log\frac{H}{u} }du\right\}\right\},
\end{align*}

respectively. From direct calculations, we can check that

\begin{equation}\label{eq:stat_err_relu}
\delta_{n,m}=\CalO\left(\sqrt{\frac{S}{mn}\log[\frac{mnH^2}{S}]}\right),\quad \Delta_n=\CalO\left(\sqrt{\frac{b^2S}{n}}|\log[\sqrt{\frac{nH^2}{b^2S}}]|\right).
\end{equation}

For the approximation error of $\Phi(H,W,S,B)$, we have the following Proposition, which is adapted from <cit.>

Let $C_0$ and $W_0$ be some constants depending only on dimension $d$ and smoothness $s$ with $s\geq 1$. Let

\begin{align*}
&H=3+2\lceil\log_2\frac{3^{d\vee\lfloor s+2\rfloor}}{\delta C_0}+5\rceil\lceil\log_2 d\vee\lfloor s+2\rfloor\rceil,\\
&W=W_0 N,\quad
S=[(H-1)W_0^2+1]N,\quad B=N^{\frac{1}{d}},
\end{align*}

where $N=\lceil |\delta\log \delta|^{-\frac{d}{s}}\rceil $. Then it holds that

\[\sup_{f\in\CalF}\inf_{f^*\in\CalF_{n,m}}\|f-f^*\|_{n}\leq \sup_{f\in\CalF}\inf_{f^*\in\CalF_{n,m}}\|f-f^*\|_{L_\infty}=\CalO(\delta).\]

From Proposition <ref>, it is clear that $H,W,S,B$ can all be treated as functions of the approximation error $\delta$. We now can apply Theorem <ref> to get the convergence rate. We first let $m=1$ and
$\delta=n^{-\frac{s}{2s+d}}\log n$
to get the associate $H,W,S,B$ accordingly. Then we substitute $H,W,S,B$ into (<ref>) to get

\[\delta_{n,m}=\CalO(n^{-\frac{s}{2s+d}}\log n)=\Delta_n.\]

Therefore, the empirical error between $\hat{f}$ and $f$ is $\CalO_{\pr}(n^{-\frac{s}{2s+1}}\log n).$ However, the non-linearity term $\hat{\CalL}$ for neural net $\Phi(H,W,S,B)$ cannot be zero. This is because,  in contrast to the previous KRR examples, the set $\Phi(H,W,S,B)$ is not convex. Consequently, during the training process, there can be multiple resulting estimators ${\hat{f}_j}$ within $\mathcal{F}_{n,m}$ that closely approximate  $f$ on the dataset $\{\BFx_i\}$ but may deviate significantly from  $f$ outside of the dataset. However, the discussion of $\hat{\CalL}$ is beyond the scope of this study. We make the assumption that the training algorithm performs effectively, resulting in $\hat{\CalL}=\CalO(n^{-\frac{1}{2}})$. In this case, the nested expectation can achieve the optimal rate for $s>d/2$:

\[|\hat \theta_{n,m}-\theta|=|\frac{1}{n}
\sum_{i=1}^n\hat{f}\circ\eta(\BFx_i)-\E[\eta\circ f(X)]| = O_{\pr}\left(n^{-\frac{1}{2}}\right)=\CalO_{\pr}(\Gamma^{-\frac{1}{2}}).\]

For the rate under risk measure, we use (<ref>) to get $\alpha=d/(d+2)$. Substitute $\alpha$, $\delta_{n,m}$, and $\Delta_n$ into Theorem <ref>, we have

\[|\hat{f}_{(\lceil\frac{\tau}{n}\rceil)}-\text{VaR}_\tau(f(X))| \leq\CalO(\Gamma^{-\frac{s}{2s+d}\frac{d}{d+2}\frac{\beta}{\gamma}}[\log \Gamma]^{-\frac{d}{d+2}\frac{\beta}{\gamma}}+n^{-\frac{1}{2\gamma}}).\]

§ NUMERICAL EXPERIMENTS

We now conduct numerical comparisons between the KRR estimator and two LSE's on sieve as we introduced in Section <ref>. The first LSE utilizes inducing points $\{\tilde{\BFx}_i\}$ to construct sieve $\mathcal{F}_{n,m} = \text{span}\{k(\tilde{\BFx}_i,\cdot)\}$, while the second one employs sieve $\mathcal{F}_{n,m}$ as a three-layer ReLU deep neural network $\Phi$ with 256 activations in the first hidden layer and 128 activations in the second hidden layer. Furthermore, it's worth noting that throughout our experiments, we have set $m=1$ because generating duplicate samples at the same realization $\BFx_i$ does not lead to an improvement in the convergence rate, as we have demonstrated in Section <ref>. The network is trained with ADMM algorithm.  All of the experiments are implemented in Matlab R2023a on a MacBook with Apple M2 Max Chip and 32GB of RAM.

Let $f(\BFx) = \E[Y|X=\BFx]$ be the conditional expectation in the inner level.
We assume

\begin{equation} \label{eq:test-fun}
f(\BFx) = \frac{1}{N}\sum_{i=1}^N c_i k(\BFx, {\BFU}_i),\quad \BFx\in\Omega,
\end{equation}

$N= 1000$, $\Omega = [0,1]^d$, and
$k$ is set  as the Laplace kernel $k_{\CalL}$ or the Gaussian kernel $k_{\CalG}$
$$k_{\CalL}(\BFx,\BFx')=e^{-\frac{\|\BFx-\BFx'\|}{d}},\quad k_{\CalG}(\BFx,\BFx')=e^{-\frac{\|\BFx-\BFx'\|^2}{d}}.$$

Both $c_i$ and $\BFU_i$ are randomly generated—with the former from $\mathsf{Normal}(0, 1)$  (i.e., the standard normal distribution) and the latter from  $\mathsf{Uniform}[0,1]^d$  (i.e., the uniform distribution on $[0,1]^d$)—and then fixed.
According to the definition of RKHSs, $f$ must be a function in the RKHS generated by the kernel $k$ because

\[\|f\|_{\CalH_k}^2=\frac{1}{N^2}\sum_{i,j=1}^Nc_ic_jk(\BFU_i,\BFU_j)<\infty.\]

All experiments can be separated into four classes as shown in the following table so each experiment can be labeled by the vector $(\CalT,d)$.

10 dimension     50 dimension

$\CalT(Z)=Z^2$    $k_{\CalL}$, $1000-5000$ data     $k_{\CalL}$, $4000-20000$ data

$\CalT(Z)=\text{VaR}(Z)$     $k_{\CalG}$, $1000-5000$ data     $k_{\CalG}$, $4000-20000$ data.

We let outer-level scenario $X$ be uniformly distributed on $[0,1]^d$. Given a realization $\BFx_i$ of $X$,
the inner-level samples are simulated via $y_{ij} = f(\BFx_i) + \varepsilon_{ij}$, where $\varepsilon_{ij}$'s are i.i.d. $\mathsf{Normal}(0, 1)$ random variables. Once an estimator $\hat{\theta}=\CalT\circ \hat{f} (X)$ of $\theta=\CalT\circ f(X)$ is constructed, we sample 10000 points for $X$ to estimate the absolute value $|\hat{\theta}-\theta|$ for comparisons.

Upper row: numerical experiments on 10-dimensional data; bottom row: numerical experiments on 50-dimensional data. The lines represent the averaged errors and the shaded area represent the STD for nested expectation and one-tenth of the STD for VaR.

For experiments in dimension $d=10$, data sizes of $n=1000, 2000, \cdots, 5000$ are used. In the case of dimension $d=50$, data sizes of $n=4000,\cdots,20000$ are used. Each experiment are repeated 1000 times and we  report the averaged absolute error, its standard deviation
(STD) when $\CalT(Z)=Z^2$, and one-tenth of its STD when $\CalT(Z)=\text{VaR}(Z)$  as our result.

For experiments with functional $\CalT(Z)=Z^2$,  the underlying function $f$ in (<ref>) is linear combinarion of Laplace kernel $k_\CalL$. In this case, the RKHS generated by $k_\CalL$ is  the $(\frac{1}{2}+\frac{d}{2})$-order Sobolev space with $d=10$ or $50$. This aligns with the case where the smoothness parameter $\nu = \frac{1}{2} < \frac{d}{2}$, and consequently, KRR with inducing points and ReLU neural network both exhibit a faster convergence rates (given that the training algorithm of neural net performs well enough), as discussed in Sections <ref> and <ref>. The KRR with inducing point has inducing point set randomly selected from data $\{\BFx_i\}$ with size $\CalS_n=n^{\frac{1}{2}}$. The ReLU neural net is as introduced at the beginning of this section. As demonstrated in Figure <ref> (a) and (d), KRR with inducing points consistently outperforms its competitors. The ReLU network performs similarly to KRR when the dimension is $d=10$ and surpasses KRR when $d=50$. According to our results in Section <ref>, both KRR with inducing points and the ReLU neural network exhibit a faster convergence rate than KRR in this case. For KRR with inducing points, our experiments successfully validate this assertion. However, for the ReLU network, we believe that its suboptimal performance is due to its high $u$-entropy when dealing with small datasets $n=1000,\cdots,5000$. When the data size is sufficiently large, i.e. $n=4000,\cdots,20000$, the ReLU network outperforms KRR and shows similar performance to KRR with inducing points.

For experiments with functional $\CalT(Z)=\text{VaR}(Z)$,  the underlying function $f$ in (<ref>) is a linear combination of Gaussian kernels $k_\mathcal{G}$. In this case, the underlying function $f$ exhibits smoothness. As discussed in Section <ref>, both KRR and KRR with inducing points can achieve the optimal convergence rate of $\mathcal{O}(\Gamma^{-\frac{1}{2}})$. In our experiments, KRR with inducing points selects the inducing point set randomly from the data ${\mathbf{x}_i}$, with a size of $\mathcal{S}_n = |\log n|^3$. From Figure <ref> (b) and (e), it is evident that KRR with inducing points consistently outperforms KRR. We attribute this to the fact that KRR with inducing points requires significantly lower computational resources due to $\mathcal{S}_n = |\log n|^3$, resulting in much lower numerical error compared to KRR. As for the ReLU neural network, it is not a smooth estimator, and thus, it can never attain the optimal convergence rate, regardless of the smoothness of the underlying function $f$ as discussed in Section <ref>. Consequently, it consistently exhibits the poorest performance among all estimators.

We have also plotted the time required for training each estimator, as shown in Figure <ref> (c) and (f). The training times for the KRR and ReLU neural network are similar for both nested expectation and VaR. Therefore, we report the average training times for KRR and the neural network over all experiments with dimension $d=10$. Regarding KRR with inducing points, in the nested expectation case, Laplace kernel is used so there are $n^{\frac{1}{2}}$ inducing points, while in the VaR case, Gaussian kernel is used so there are $|\log n|^3$ inducing points. Therefore, we  separately plot these different scenarios. The running times are plotted on a logarithmic scale. It's evident that Gaussian kernel KRR with inducing points has the shortest running time among all cases. This is because it only requires $|\log n|^3$ inducing points, making it a highly computational efficient estimator.

§ CONCLUSION

In this paper, we generalize the KRR for nested simulation to LSE on sieve. We demonstrate that, for various scenarios of nested simulation, more accurate LSEs on sieve exist, surpassing KRR and the traditional Monte Carlo method in terms of performance. We relax the conditions required for KRR to achieve the square root convergence rate, which is the standard rate for conventional Monte Carlo simulations. Our theoretical framework for convergence rate analysis is more versatile and can accommodate different forms of nested simulation.

Our work can be extended in several ways. Firstly, we have shown that the convergence rate of neural network also depends the training algorithm. However, how to select training algorithm to achieve optimal convergence rate is unknown. Furthermore, there may be opportunities for enhancements in our upper bounds on the convergence rates for the VaR functional. This is primarily because the convergence rate under the VaR metric is closely related to the structure, or more precisely, the topology of the class $\CalF$ and the sieve $\CalF_{n,m}$. It is evident that the concept “entropy” in our study may not be the ideal tool for capturing the topology information, so achieving better results than those presented in [7] for the convergence rate under VaR may be challenging.

§ PROOF OF THEOREM <REF>

Let $\bar{\sigma}_i$ be the sub-Gaussianilty of the average $\bar{\varepsilon}_i=\frac{1}{m}\sum_{j=1}^m\varepsilon_{ij}$ where $\varepsilon_{ij}$ are the zero-mean sub-Gaussian noise satisfying assumption <ref>. Then
$$\max_{i=1,\cdots,n}\bar{\sigma_{i}}^2\leq C{\sigma}^2/m,$$

where $C$ is some universal constant.

Lemma <ref>  is a direct result of the Bernstein inequality.

For any $\{h_i\in\Real\}_{i=1}^n$ and $\delta>0$,

\[\pr\left(|\sum_{i=1}^n\bar{\varepsilon}_ih_i|\geq \delta\right)\leq 2\exp\left[-\frac{\delta^2 m}{C \sum_{i=1}^nh_i^2}\right]\]

where $\bar{\varepsilon}_i=1/m\sum_{j=1}^m\varepsilon_{ij}$ are the averaged noise in Lemma <ref> and $C$ is some universal constant.

From the sub-Gaussianilty of $\bar{\varepsilon}$ and Lemma <ref>, we can show that for any $\omega$

\[\E \exp\left[\omega \sum_{i=1}^nh_i\bar{\varepsilon}_i\right] =\prod_{i=1}^n\E \exp\left[\omega\bar{\varepsilon}_ih_i\right]\leq \exp\left[C_1\frac{\sigma^2\omega^2}{m}\sum_{i=1}^nh_i^2\right],\]

where $C_1>0$ is the universal constant in Lemma <ref>. We then use Chebyshev's inequality to show that for any $\omega>0$

\[\pr\left(\sum_{i=1}^n\bar{\varepsilon}_ih_i\geq \delta\right)\leq \exp\left[C_1\frac{\sigma^2\omega^2}{m}\sum_{i=1}^nh_i^2-\delta\omega\right].\]

Let $\omega=\delta m/(2C_1\sigma^2\sum_{i=1}^nh_i^2)$ to minimize the quantity $C_1\frac{\sigma^2\omega^2}{m}\sum_{i=1}^nh_i^2-\delta\omega$, we then have the following inequality independent of $\omega$

\[\pr\left(\sum_{i=1}^n\bar{\varepsilon}_ih_i\geq \delta\right)\leq \exp\left[-\frac{\delta^2m}{4C_1\sigma^2\sum_{i=1}^nh_i^2}\right].\]

Because the same inequality holds for $-\sum_{i=1}^n\bar{\varepsilon}_ih_i$, the lemma follows.

The following theorem is adapted from Lemma 3.2 in  [11]:
Suppose $W_i$ are i.i.d. $\sigma$-sub-Gaussian r.v.s. Assume $\sup_{h\in\CalH}\|h\|_n\leq R$. There exist constants $C$ and $C_1$  depending only on $\sigma$ such that for all $\delta>0$ and $K>0$ satisfying

\[\sqrt{n}\delta\geq 2C\int^{R}_{\delta/(2^3K)}\sqrt{H(u,\CalH,\|\cdot\|_n)}du\vee R,\]

we have

\begin{align*}
\pr\left(\{\sup_{h\in\CalH}|\langle W, h \rangle_n|\geq\delta\}\cap\{\|W\|_n\leq K\}\right)\leq C_1\exp\left[-\frac{n\delta^2}{1152\sigma^2R^2}\right].
\end{align*}

Suppose $\sup_{h\in\CalH}\|h\|_n\leq R$. Then for some constants $C_1$, $C_2$ depending only on $\sigma$, and  for $\bar{\delta},\underline{\delta}>0$ satisfying $\bar{\delta}/\underline{\delta}<R$ and

\[\sqrt{nm}\bar{\delta}\geq  2C_1\int^R_{\bar{\delta}/(2^3\underline{\delta})}\sqrt{H(u,\CalH,\|\cdot\|_n)}du\vee R,\]

we have

\[\pr\left(\{\sup_{h\in\CalH}|\langle \bar{\varepsilon},h\rangle_n|\geq \bar{\delta}\}\cap\{\|\bar{\varepsilon}\|_n^2\leq\underline{\delta}\} \right)\leq C_1 \exp\left[-\frac{mn{\bar{\delta}}^2}{C_2R^2}\right]\]

where $\bar{\varepsilon}(\BFx_i)=1/m\sum_{j=1}^m\varepsilon_{ij}$ are the averaged noise in Lemma <ref>.

Substitute Lemma <ref> for the sub-Gaussianality $\sigma$ in Lemma <ref>, we can immediately have the result.

For simplicity, write $\hat{f}_{n,m}$ as $\hat{f}$, and write $f^*_{n,m}$ as $f^*$. Because $\hat{f}$ is the miniimizer of the problem $\min_{h\in\CalF_{n,m}}\|h-\bar{y} \|_n^2=\min_{h\in\CalF_{n,m}}\|h-f-\bar{\varepsilon} \|_n^2$, we can have $\|\hat{f}-f-\bar{\varepsilon}\|_n^2\leq \|{f}^*-f-\bar{\varepsilon}\|_n^2$. This inequality can be rewritten as

\begin{equation}\label{eq:thm1_pf_1}
\|\hat{f}-f\|_n^2\leq 2\langle \bar{\varepsilon},\hat{f}-f^*\rangle_n+\|f^*-f\|_n^2.
\end{equation}

By triangle inequality, we can further have

\begin{align}
\|\hat{f}-f^*\|_n^2&\leq 2\|\hat{f}-f\|^2_n+2\|f^*-f\|^2_n\nonumber\\
&\leq 4\langle \bar{\varepsilon},\hat{f}-f^*\rangle_n+4\|f^*-f\|_n^2\label{eq:thm1_pf_2}
\end{align}

where the second line is from (<ref>).

If $\|f^*-f\|_n^2\geq \langle \bar{\varepsilon},\hat{f}-f^*\rangle_n $, (<ref>) can be derived directly from (<ref>).

If $\|f^*-f\|_n^2\leq \langle \bar{\varepsilon},\hat{f}-f^*\rangle_n $, we have

\begin{equation}\label{eq:thm1_pf_3}
\|\hat{f}-f^*\|_n^2\leq 8  \langle \bar{\varepsilon},\hat{f}-f^*\rangle_n.
\end{equation}

By the Cauchy-Schwarz inequality,  we can deduce from (<ref>) that $\|\hat{f}-f^*\|_n\leq 8\underline{\delta}$ on the event $\|\bar{\varepsilon}\|_n\leq \underline{\delta}$. On this event, it suffices to discuss the probability of the event $\|\hat{f}-f^*\|_n\leq \underline{\delta}$ for $\bar{\delta}\leq 8\underline{\delta}$. Let $S=\min\{s=0,1,\cdots:2^s\bar{\delta}\geq 8\underline{\delta}\}$. Then using the peeling technique, we have

\begin{align}
&\pr\left( \{\|\hat{f}-f^*\|_n^2\geq \bar{\delta}^2\}\cap\{\|\bar{\varepsilon}\|_n\leq\underline{\delta}\}\right)\\
\leq{}&  \pr\left( \{\langle \bar{\varepsilon},\hat{f}-f^*\rangle_n\geq \bar{\delta}^2/8\}\cap\{\|\bar{\varepsilon}\|_n\leq\underline{\delta}\}\right)\nonumber\\
%&\leq \pr\left(\bigcup_{s=0}^\infty\{\sup_{h\in\CalF_n(2^{s+1}\delta)}\langle \bar{\varepsilon},h-f^*\rangle_n\geq 2^{2s-1}\delta^2\}\right)\nonumber\\
={}& \pr\left( \bigcup_{s= 0}^S\{2^{2s+2}\bar{\delta}^2/8\geq \langle \bar{\varepsilon},\hat{f}-f^*\rangle_n\geq 2^{2s}\bar{\delta}^2/8\}\cap\{\|\bar{\varepsilon}\|_n\leq\underline{\delta}\}\right)\nonumber\\
%\leq{}& \sum_{s\geq 0} \pr\left( 2^{2s+2}\delta^2\geq \langle \bar{\varepsilon},\hat{f}-f^*\rangle_n\geq 2^{2s}\delta^2\right)\nonumber\\
\leq{}& \sum_{s=0}^S\pr\left(\sup_{h\in\CalF_n(2^{s+1}\bar{\delta}/\sqrt{8})}\langle \bar{\varepsilon},h-f^*\rangle_n\geq 2^{2s-3}\bar{\delta^2}\cap\{\|\bar{\varepsilon}\|_n\leq\underline{\delta}\}\right)\label{eq:thm1_pf_4}.
\end{align}

Notice that if we set $\underline{\delta}=\sigma/\sqrt{m}$ and choose $\sqrt{nm}\delta_{n,m}^2\geq c\Psi(\delta_{n,m}) $, i.e., $\delta_{n,m}$ satisfies condition (<ref>), then for any $2^{s+1}{\delta}\geq \delta_{n,m}$ with $2^{s+1}\delta<2^7\sigma/\sqrt{m}$, condition (<ref>) is also satisfied:

\[\sqrt{nm}2^{2s+2}\delta^2\geq c\Psi(2^{s+1}\delta)\]

for $\Psi(\delta)/\delta^2$ is a non-decreasing function of $\delta$ for $0<\delta<2^7\sigma/\sqrt{m}$. We then can apply Proposition <ref> to each term of the summation (<ref>):

\begin{align*}
\pr\left(\{\sup_{h\in\CalF_n(\frac{2^{s+1}\bar{\delta}}{\sqrt{8}})}\langle \bar{\varepsilon},h-f^*\rangle_n\geq 2^{2s-3}\bar{\delta}^2\}\cap\{\|\bar{\varepsilon}\|_n\leq\underline{\delta}\}\right)\leq C_1\exp\left[-\frac{mn2^{4s}\bar\delta^4}{C_22^{2s+2}\bar\delta^2}\right]
\end{align*}

where $C_1,C_2$ only depend on $\sigma$. So (<ref>) becomes

\begin{align}
\pr\left( \{\|\hat{f}-f^*\|_n\geq \delta_{n,m}\}\cap\{\|\bar{\varepsilon}\|_n\leq \frac{\sigma}{\sqrt{m}} \}\right)&\leq \sum_{s= 0}^S C\exp\left[-\frac{mn2^{2s}\delta_{n,m}^2}{4C}\right]\nonumber\\
&\leq C\exp\left[-\frac{mn\delta_{n,m}^2}{C}\right] \label{eq:thm1_pf_5}.
\end{align}

where $C$  only depends on $\sigma$.  So

\[ \pr\left( \{\|\hat{f}-f^*\|_n\geq \delta_{n,m}\right)\leq C\exp\left[-\frac{mn\delta_{n,m}^2}{C}\right]  +\pr\left(\|\bar{\varepsilon}\|_n\geq\frac{\sigma}{\sqrt{m}}\right). \]

§ PROOF OF THEOREM <REF>

If $m$ satisfies
$\sqrt{m}(\delta_{n,m}+\|{f}^*-f\|_n)=\CalO(1)$, $\sqrt{m}\|\hat{f}-f^*\|_n=\CalO(1)$.

From triangle inequality:

\begin{align*}
\sqrt{m}\|\hat{f}-f^*\|_n&\leq \sqrt{m}\left(\|\hat{f}-f\|_n+\|{f}^*-f\|_n\right)\\
&\leq \sqrt{m}\left(\delta_{n,m}+2\|{f}^*-f\|_n\right)=\CalO(1).
\end{align*}

By the triangle inequality, we have

\begin{align}\label{eq:thm2_1}
|\hat \theta_{n,m}-\theta| ={}& \left|\E [\eta(f(X))] -\frac{1}{n}\sum_{i=1}^n \eta(f(\BFx_i))+ \frac{1}{n}\sum_{i=1}^n \eta(f(\BFx_i))- \frac{1}{n}\sum_{i=1}^n \eta(\hat{f}(\BFx_i))\right|\nonumber\\
\leq{}& \underbrace{\left|\E [\eta(f(X))] -\frac{1}{n}\sum_{i=1}^n \eta(f(\BFx_i))\right|}_{I_1}
+ \underbrace{ \left|\frac{1}{n}\sum_{i=1}^n \bigl[\eta(f(\BFx_i))- \eta(\hat{f}(\BFx_i)) \bigr] \right| }_{I_2}.
% \leq{}& \left|\E [\eta(f(X))] -\frac{1}{n}\sum_{i=1}^n \eta(f(\BFx_i))\right|+  \left|\frac{1}{n}\sum_{i=1}^n \eta(f(\BFx_i))- \frac{1}{n}\sum_{i=1}^n \eta(\hat{f}(\BFx_i))\right|\nonumber\\
% \coloneqq {}& I_1 + I_2.
\end{align}

Assumption <ref> and conditions on $\eta$ imply that both $\|f\|_{L_\infty} $ and $\|\eta(f(\cdot))\|_{L_\infty}$ are finite.
Therefore, $\eta(f(\BFx_i))$'s are bounded random variables, thereby being sub-Gaussian.
Then, the central limit theorem  implies that

\begin{align}\label{eq:thm2_2}
I_1 = O_{\pr}(n^{-1/2}).
\end{align}

It remains to bound $I_2$.  It follows from Taylor's expansion and the triangle inequality that

\begin{align*}
I_2 ={}& \left|\frac{1}{n}\sum_{i=1}^n \eta'(f(\BFx_i))(f(\BFx_i)-\hat{f}(\BFx_i)) + \frac{1}{2n}\sum_{i=1}^n \eta''(\tilde{z}_i)(f(\BFx_i)-\hat{f}(\BFx_i))^2 \right|\\
\leq{}& \underbrace{\left|\frac{1}{n}\sum_{i=1}^n \eta'(f(\BFx_i))(f(\BFx_i)-\hat{f}(\BFx_i))\right|}_{I_{21}}
+ \underbrace{\left|\frac{1}{2n}\sum_{i=1}^n \eta''(\tilde{z}_i)(f(\BFx_i)-\hat{f}(\BFx_i))^2 \right|}_{I_{22}},
\end{align*}

where $\tilde{z}_i$ is a value between $f(\BFx_i)$.

For $I_{21}$, we use the Triangle inequality again to have

\begin{align}
&  {\left|\frac{1}{n}\sum_{i=1}^n \eta'(f(\BFx_i))(f(\BFx_i)-\hat{f}(\BFx_i))\right|}\nonumber\\
\leq{}&\left|\frac{1}{n}\sum_{i=1}^n \eta'(f(\BFx_i))(f(\BFx_i)-{f}^*(\BFx_i))\right| +\left|\frac{1}{n}\sum_{i=1}^n \eta'(f(\BFx_i))(f^*(\BFx_i)-\hat{f}(\BFx_i))\right| \nonumber\\
\leq{}& \|\eta'\circ f\|_n\|f-f^*\|_n+\left|\frac{1}{n}\sum_{i=1}^n \eta'(f(\BFx_i))(f^*(\BFx_i)-\hat{f}(\BFx_i))\right|, \label{eq:thm2_pf_I21}
\end{align}

where the first term in the  last line is from Cauchy–Schwarz inequality. For the second term in the last line, we further decompose it as follows:

\begin{align}
& \left|\frac{1}{n}\sum_{i=1}^n \eta'(f(\BFx_i))(f^*(\BFx_i)-\hat{f}(\BFx_i))\right| \nonumber\\
\leq {}& \left|\frac{1}{n}\sum_{i=1}^n \left(\eta'(f(\BFx_i))-\E \eta'(f(\BFx_i)) \right)(f^*(\BFx_i)-\hat{f}(\BFx_i))\right|\label{eq:thm2_3}\\
&\quad +\left|\frac{1}{n}\sum_{i=1}^n \E\left[\eta'(f(\BFx_i))\right](f^*(\BFx_i)-\hat{f}(\BFx_i))\right| \label{eq:thm2_4}.
\end{align}

For (<ref>), because we have assumed that any function in $\CalF$ and $\CalF_{n,m}$ is bounded by $b$ and $P_X$-Donsker and, in the later proofs for Theorem <ref>, we will show that $|f^*(\BFx_i)-\hat{f}(\BFx_i)|\leq\|\hat{f}-f^*\|_{L_\infty}=o_{\pr}(1)$. Therefore, it is straightforward to use central limit theorem to get the following bound directly:

\begin{align}
&\left|\frac{1}{n}\sum_{i=1}^n \E\left[\eta'(f(\BFx_i))\right](f^*(\BFx_i)-\hat{f}(\BFx_i))\right|\nonumber\\
\leq& \left|\E\left[\eta'(f(X))\right]\E_{X,\{\varepsilon_i\}} [f^*(X)-\hat{f}(X)]\right|+\CalO_{\pr}(n^{-\frac{1}{2}})\label{eq:thm2_5}.
\end{align}

For (<ref>), because $\eta'$ is bounded so $W_i=\eta'(f(\BFx_i))-\E \eta'(f(\BFx_i))$ is a zero-mean sub-Gaussian r.v. for each $i$.  Also, from Lemma <ref>, we can notice that for any $m$ satisfying condition $\sqrt{m}(\delta_{n,m}+\|f-f^*\|_n)=\CalO(1)$, then $\|f^*-\hat{f}\|_n\leq R/\sqrt{m}$ for some constant $R$.  Let $\sigma_\eta$ denote the sub-Gaussianality of $W_i$. Then from Lemma <ref>, we have for any $\delta$ satisfying $R/\sqrt{m}>\delta/\sigma_\eta$ and

\begin{equation}\label{eq:thm2_condition_W}
\sqrt{n}\delta\geq 2C\left(\int^{R/\sqrt{m}}_{\delta/(8\sigma_\eta)}\sqrt{H(u,\CalF_{n,m},\|\cdot\|_n)}du\vee \frac{R}{\sqrt{m}}\right)
\end{equation}

we have

\begin{align*}
\pr\left(\langle W,f^*-\hat{f}\rangle_n>\delta\}\right)
%=&\pr\left(\bigcup_{s=1}^S\{2^{s} \delta\geq\langle W,f^*-\hat{f}\rangle_n\geq 2^{s-1}\delta\}\cap\{\|W\|_n\leq \sigma_\eta\}\right)\\
%\leq& \sum_{s=0}^S\pr\left(\sup_{h\in\CalF_n(2^{s}\delta)}\langle W,h\rangle_n\geq 2^{s-1}\delta\right)\\
\leq   C_3\exp\left[-\frac{mn\delta^2}{C_4R^2}
\right]+\pr\left(\|W\|_n\leq \sigma_\eta\right).
\end{align*}

Recall that $\delta_{n,m}$ satisfies $\sqrt{m}\delta_{n,m}=\CalO(1)$ and

\[\sqrt{mn}\delta^2_{n,m}\geq c\left(\int_{\frac{\delta_{n,m}^2\sqrt{m}}{2^7{\sigma}}}^{{\delta_{n,m}}}\sqrt{H(u,\CalF_n({\delta}_{n,m}),\|\cdot\|_n)}du \vee {\delta}_{n,m}\right).\]

Let $\delta=\delta_{n,m}^2$, then $\delta$ satisfies condition (<ref>). So

\begin{equation}\label{eq:thm2_6}
\left|\frac{1}{n}\sum_{i=1}^n \left(\eta'(f(\BFx_i))-\E \eta'(f(\BFx_i)) \right)(f^*(\BFx_i)-\hat{f}(\BFx_i))\right|\leq \CalO_{\pr}(\delta_{n,m}^2).
\end{equation}

For $I_{22}$, from the boundedness of $\eta''$, it is straightforward to derive that

\begin{equation}
I_{22}\leq \frac{\|\eta''\|_{L_\infty}}{2}\|f-\hat{f}\|_n^2\leq\CalO_{\pr}(\delta^2_{n,m}+\|f-f^*\|_n^2)\label{eq:thm2_7}.
\end{equation}

Substitute (<ref>), (<ref>), (<ref>), and (<ref>) into $I_2$, together with (<ref>) for $I_1$, we can have the final result.

§ PROOF OF THEOREM <REF>

We first define two quantities known as the empirical Rademacher complexity and $L_2$-Rademacher complexity of a space $\CalH$, respectively:

\begin{align*}
&\hat{\CalR}_n(\delta,\CalH)=\E_{r}\left[\sup_{h\in\CalH\cap B(\delta,\|\cdot\|_n)} \langle r,h\rangle_n\right],\\
&{\CalR}_n(\delta,\CalH)=\E_{\{\BFx_i\},r}\left[\sup_{h\in\CalH\cap B(\delta,\|\cdot\|_{L_2})} \langle r,h\rangle_n\right]
\end{align*}

where $B(\delta,\|\cdot\|)=\{h: \|h\|\leq \delta\}$ is the ball with radius $\delta$ under norm $\|\cdot\|$ and $r(\BFx_i)=r_i$ are i.i.d. Rademacher r.v.s, i.e., $\pr(r_i=1)=\pr(r_i=-1)=1/2$. According to the Dudley's entropy integral bound,

\begin{equation}\label{eq:Dudley}
\hat{\CalR}_n(\delta,\CalH)\leq \inf_{\varepsilon\geq 0}\left\{4\varepsilon+\frac{C}{\sqrt{n}}\int_\varepsilon^{\delta} \sqrt{H(u,\CalH\cap B(\delta,\|\cdot\|_n),\|\cdot\|_n)}du\right\},
\end{equation}

where $C>0$ is some universal constant. Please refer to <cit.> for details.

We now can state the following theorem regarding the relationships between the empirical and $L_2$ norms/Rademacher complexities. This theorem is adapted from Theorem 14.1 and Proposition 14.25 of [24]:

Suppose $\sup_{h\in\CalH}\|h\|_{\infty}\leq b$ and $\CalH$ is a star-shaped. Let

\begin{align}\label{eq:thm2_wainwright}
&\hat{\delta}_n=\inf\{\delta: \hat{\CalR}_n(\delta,\CalH)\leq \delta^2/b\},\ {\delta}_n=\inf\{\delta: {\CalR}_n(\delta,\CalH)\leq \delta^2/b\}.
\end{align}

Assume $\log\log 1/\delta=\CalO(n\delta_n^2)$. Then there exist constants $m_1$,$m_2$, $c_1$ and $a$ such that with probability at least $1-ae^{-n\delta_n^2/(c_0b)}$:

\begin{align*}
&m_1\delta_n\leq \hat{\delta}_n\leq m_2\delta_n,\\
&\sup_{h\in\CalH}\left|\|h\|_n-\|h\|_{L_2}\right|\leq c_1\delta_n.
\end{align*}

Remind that $\CalF_{n,m}$ and $\CalF$ satisfy Assumptions <ref> and <ref>. Using Theorem <ref>, we can immediately derive the following bound regarding the $L_2$ and empirical norms for any function in $\CalF\cup\CalF_{n,m}$:

There exist a universal constant $C$ such that for $\delta>0$ satisfying:

\begin{equation}\label{eq:thm3_coro_condition}
\frac{\delta^2}{b}\geq \inf_{\varepsilon\geq 0}\left\{4\varepsilon+\frac{C}{\sqrt{n}}\int_\varepsilon^{\delta} \sqrt{H(u,\CalH\cap B(\delta,\|\cdot\|_n),\|\cdot\|_n)}du\right\},
\end{equation}

\[\sup_{h\in\CalF\cup\CalF_{n,m}}\|h\|_{L_2}\leq \sup_{h\in\CalF\cup\CalF_{n,m}}\|h\|_n+\CalO_{\pr}(\delta).\]

According to Assumption <ref>, $H(u,\CalF_{n,m},\|\cdot\|_n)\leq H(u,\CalF,\|\cdot\|_n)$. So it is straightforward to derive that $H(u,\CalF\cup\CalF_{n,m},\|\cdot\|_n) \leq 2H(u,\CalF,\|\cdot\|_n)$ according to the definition of $u$-entropy.  Using Dudley's integral (<ref>), we have

\begin{align*}
\hat{\CalR}_n(\delta,\CalF\cup\CalF_{n,m})&\leq \inf_{\varepsilon\geq 0}\left\{4\varepsilon+\frac{C}{\sqrt{n}}\int_\varepsilon^{\delta} \sqrt{H(u,(\CalF\cup\CalF_{n,m})\cap B(\delta,\|\cdot\|_n),\|\cdot\|_n)}du\right\}
\end{align*}

where $C$ is some universal constant. Therefore, for any $\delta$ satisfying the condition (<ref>), $\delta$ must be the upper bound of $\hat{\delta}_n$ in (<ref>) up to a multiplicative constant. Then according to Assumption <ref> and Theorem <ref>, we have

\begin{align*}
\sup_{h\in\CalF\cup\CalF_{n,m}}\|h\|_{L_2}&= \sup_{h\in\CalF\cup\CalF_{n,m}}\|h\|_{n}+\|h\|_{L_2}-\|h\|_{n}\\
&\leq \sup_{h\in\CalF\cup\CalF_{n,m}}\|h\|_{n}+ \sup_{h\in\CalF\cup\CalF_{n,m}}\left||h\|_{L_2}-\|h\|_{n}\right|\\
&\leq \sup_{h\in\CalF\cup\CalF_{n,m}}\|h\|_n+\CalO_{\pr}(\delta).
\end{align*}

Suppose Assumption <ref>–<ref> holds. Let $\delta_{n,m}$, $f^*$, and $\hat{f}$ be as defined in Theorem <ref>. Then

\[\rho_n:=\|\hat{f}-f^*\|_{L_\infty}\leq \CalO_{\pr}\left(\left|\delta_{n,m}+\|f^*-f\|_n+{\delta}/{\sqrt{m}}\right|^\alpha\right)\]

where $\delta$ satisfies condition (<ref>) and $m$ satisfies $\sqrt{m}(\delta_{n,m}+\|f^*-f\|_n)=\CalO(1)$.

Lemma <ref> is a direct result of Theorem <ref>, Assumption <ref>, and Lemma <ref>.

The proof of Theorem <ref> is merely repeating the proof of Theorem 4 of [7]. The only difference is that the term $l_n$ in the proof of Theorem 4 of [7] is

\begin{align*}
l_n = \left|\delta_{n,m}+\|f^*-f\|_n+{\delta}/{\sqrt{m}}\right|^\alpha
\end{align*}

in our case.

#1ISBN #1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1<https://doi.org/#1>et al.#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1<><#>1#1#1#1#1#1#1#1#1#1#1#1#1#1PreBibitemsHook

[1]

Asmussen, S.,
Glynn, P.W.:
Stochastic Simulation: Algorithm and Analysis.

[2]

Fu, M.C.,
Hu, J.-Q.:
Conditional Monte Carlo: Gradient Estimation and Optimization

[3]

Fu, M.C.,
Hong, L.J.,
Hu, J.-Q.:
Conditional Monte Carlo estimation of quantile sensitivities.
Manag. Sci.

[4]

Gordy, M.B.,
Juneja, S.:
Nested simulation in portfolio risk measurement.
Manag. Sci.

[5]

Zhang, K.,
Liu, G.,
Wang, S.:
Technical note—Bootstrap-based budget allocation for nested
Oper. Res.

[6]

Hong, L.J.,
Juneja, S.,
Liu, G.:
Kernel smoothing for nested estimation with application to portfolio
risk measurement.
Oper. Res.

[7]

Wang, W.,
Wang, Y.,
Zhang, X.:
Smooth Nested Simulation: Bridging Cubic and Square Root Convergence Rates in
High Dimensions

[8]

Hastie, T.:
Ridge regularization: An essential concept in data science.

[9]

Kanagawa, M.,
Hennig, P.,
Sejdinovic, D.,
Sriperumbudur, B.K.:
Gaussian Processes and Kernel Methods: A Review on Connections and
Preprint available at <https://arxiv.org/abs/1807.02582>

[10]

Adams, R.A.,
Fournier, J.J.:
Sobolev Spaces,
2nd edn.
Academic Press,
Cambridge, MA

[11]

Geer, S.:
Empirical Processes in M-Estimation.
Cambridge University Press,

[12]

Donoho, D.L.:
De-noising by soft-thresholding.
IEEE transactions on information theory

[13]

Geer, S.:
On the uniform convergence of empirical norms and inner products,
with application to causal inference.
Electronic Journal of Statistics

[14]

Ding, L.,
Hu, T.,
Jiang, J.,
Li, D.,
Wang, W.,
Yao, Y.:
Random Smoothing Regularization in Kernel Gradient Descent Learning.
Preprint available at <https://arxiv.org/abs/2305.03531>

[15]

Birman, M.S.,
Solomjak, M.Z.:
Piecewise-polynomial approximations of functions of the classes.
Mathematics of the USSR-Sbornik

[16]

Wu, Z.-m.,
Schaback, R.:
Local error estimates for radial basis function interpolation of
scattered data.
IMA journal of Numerical Analysis

[17]

Davies, A.J.:
Effective implementation of gaussian process regression for machine learning.
PhD thesis,
University of Cambridge

[18]

Ding, L.,
Tuo, R.,
Shahrampour, S.:
Generalization guarantees for sparse kernel approximation with
entropic optimal features.
In: International Conference on Machine Learning,
pp. 2545–2555

[19]

Kühn, T.:
Covering numbers of gaussian reproducing kernel hilbert spaces.
Journal of Complexity

[20]

Wendland, H.:
Scattered Data Approximation.
Cambridge University Press,

[21]

Gouk, H.,
Frank, E.,
Pfahringer, B.,
Cree, M.J.:
Regularisation of neural networks by enforcing lipschitz continuity.
Machine Learning

[22]

Vaart, A.W.:
Asymptotic Statistics.
Cambridge University Press,

[23]

Suzuki, T.:
Adaptivity of deep relu network for learning in besov and mixed smooth besov
spaces: optimal rate and curse of dimensionality.
arXiv preprint arXiv:1810.08033

[24]

Wainwright, M.J.:
High-Dimensional Statistics: A Non-Asymptotic Viewpoint.
Cambridge University Press,


Further author information: (Send correspondence to G.M. and A.S.C.)
G.M.: E-mail<EMAIL_ADDRESS>
A.S.C.: E-mail<EMAIL_ADDRESS>
# Deep and shallow data science for multi-scale optical neuroscience

Gal Mishne  Halıcıoğlu Data Science Institute, Department of Electrical and
Computer Engineering and the Neurosciences Graduate Program, UC San Diego,
9500 Gilman Drive, La Jolla, CA 92093 USA  Adam Charles  Department of
Biomedical Engineering, Kavli Neuroscience Discovery Institute, Center for
Imaging Science, Department of Neuroscience, and Mathematical Institute for
Data Science, Johns Hopkins University, Baltimore, MD 21287 USA

###### Abstract

Optical imaging of the brain has expanded dramatically in the past two
decades. New optics, indicators, and experimental paradigms are now enabling
in-vivo imaging from the synaptic to the cortex-wide scales. To match the
resulting flood of data across scales, computational methods are continuously
being developed to meet the need of extracting biologically relevant
information. In this pursuit, challenges arise in some domains (e.g., SNR and
resolution limits in micron-scale data) that require specialized algorithms.
These algorithms can, for example, make use of state-of-the-art machine
learning to maximally learn the details of a given scale to optimize the
processing pipeline. In contrast, other methods, however, such as graph signal
processing, seek to abstract away from some of the details that are scale-
specific to provide solutions to specific sub-problems common across scales of
neuroimaging. Here we discuss limitations and tradeoffs in algorithmic design
with the goal of identifying how data quality and variability can hamper
algorithm use and dissemination.

###### keywords:

fluorescence microscopy, calcium imaging, functional imaging, data analysis

##  1 INTRODUCTION

Data science, in the form of statistical modeling and image processing, has
become a central topic in extracting the most out of complex imaging data,
such as taken by the latest in optical designs [1]. Computational approaches
in particular have blossomed as imaging methods have extended across the
micro- and macro- scales that aim to glean information about how the brain
processes information at all scales [2]. These methods fundamentally take into
account—implicitly or explicitly—properties of the data imbued by the imaging
target and the imaging device. Colloquially we refer to these as the
likelihood and the prior, after the successful terminology from the Bayesian
inference and estimation literature.

Due to the enormous range of targets (different brain areas across species) at
different scales (cell bodies, widefield, dendrites, axons, vasculature) that
exhibit different temporal dynamics (different indicators of calcium and
voltage, hemodynamics, long-term AMPA receptor changes, etc.), fundamentally
different tools have been developed by the community. Some of these tools have
aimed to solve core processing challenges, for example, denoising or
identifying functional components in videos, across many different imaging
contexts. Other tools instead have focused on concentrating effort toward
solving a very specific set of challenges in a single context. As an extreme
example of the latter, computational imaging techniques use specialized
algorithms to process and interpret data from one specific instance of an
imaging device.

Inherent to these methods are thus a set of tradeoffs surrounding the ability
of these algorithms to generalize across contexts. These tradeoffs can impact
how new methods can be robustly assessed and shared. We discuss here a number
of considerations one should consider both in the development of such
algorithms and in the deployment and dissemination thereof.

##  2 The promise of deep learning?

One increasingly popular group of approaches in the analysis of functional
optical imaging is to leverage advances in deep learning, where deep learning
has been integrated into various points of the processing pipeline for
functional imaging, from preprocessing, through inference to postprocessing,
which we briefly review below.

One challenge deep learning methods have sought to address is denoising of the
acquired video. Functional imaging suffers from imaging-dependent noise, where
signal-to-noise (SNR) can vary across different neuronal components due to
differences in the expression of the indicator. Identifying neural components,
e.g., cell bodies, is enabled because the location of cells is fixed in the
field-of-view and their activity results in fluorescence persisting over
multiple frames. Extraction of spatial components and their corresponding time
traces can be improved by incorporating a denoising step as part of
preprocessing (following motion correction). Recently, general advances from
deep learning [3, 4, 5] have been proposed for denoising calcium imaging. One
class of methods, including DeepInterpolation [3], SUPPORT [6] and DeepVID [7,
8], leverage the spatio-temporal consistency of the fluorescence signal to
separate it from the noise. In essence, thse methods train a neural network to
predict a movie frame (for DeepInterpolation) or pixel (for SUPPORT) based on
previous and past frames or the surrounding pixels in space and time,
respectively. As independent noise cannot be predicted from the spatio-
temporal neighborhood, it is thus effectively removed. A broader class of
image-restoration methods are based on more traditional static architectures.
CARE [9], and the more recetnt XTC that specializes for sypatic imaging [10],
are based on a U-net architecture trained on high- and low-resolution images.
These images can be obtained by more expensive measurements in vivo, or as in
XTC, by creating training data ex-vivo using super-resolution microscopy.

The fundamental step for the analysis of functional imaging data is
identifying spatial regions of interest (ROIs) and their corresponding time
traces. While earlier methods mostly relied on regularized matrix
decomposition [11, 12, 13, 14, 15], recently, general deep learning methods
for image segmentation have been adapted to the ROI extraction problem [16,
17, 18, 19, 20, 21, 22]. These are supervised approaches that require labeled
training data in order to train a network to identify ROIs in new datasets.

Some networks tackle only the problem of identifying spatial components and
ignore the temporal dimension. UNet2DS [16] is based on the fully
convolutional UNet [17] model, which was developed for biomedical image
segmentation. UNet2DS takes the mean image as an input and outputs the
probability of a pixel belonging to a cell or background. UNet2DS, however,
ignores important temporal information, leading to a difficulty in separating
overlapping neurons. Furthermore, it is dependent on the summary image that
may be biased towards cells that fire more often, or have high expression
despite exhibiting no activity. Similar to UNet2DS, Cellpose [22] is similarly
based on a UNet architecture, with an added global “style” designator to help
differentiate different image classes, and so suffers the same drawbacks.

Other networks have added temporal statistics or activity to their design [18,
19, 20]. A (2+1)D convolutional neural network [18] was trained on
spatiotemporal sliding windows and the output represented the probability of a
pixel belonging to an ROI centroid. This approach better suppressed noisy
detections compared to a 2D network that only took as input the time-averaged
image, yet still only provided spatial components as outputs.

STNeuroNett [19] is a 3D convolutional neural network for neuron segmentation,
trained on overlapping spatiotemporal blocks and outputs a binarized
probability map. Shallow U-Net Neuron Segmentation (SUNS) [20] simplifies
STNeuroNet’s architecture for speed and incorporates temporal information by
applying a temporal matched filter and whitening to the input video. Both
methods post-process the output probability maps by clustering to detect
individual cells and aggregating detected cells across the maps. The design of
these networks incorporates temporal information while efficiently reducing it
in order to obtain fast processing speeds and low-latency inference in closed-
loop experiments. Finally, DISCo (Deep learning, Instance Segmentation, and
Correlations) [21] uses a combination of time-trace correlation-based pixel
segmentation on a graph and a convolutional neural network to identify
individual spatial profiles. While these these methods incorporate temporal
information in different ways, they only identify cells and do not address the
issue of estimating time traces.

Finally, deep learning has also served in postprocessing to classify
identified regions-of-interest (ROIs) as either true or false, based on human-
annotated training data, to improve the precision of calcium imaging analysis
methods, e.g., in CaImAn [12].

On the practical side, training these networks might involve a lot of training
data and time. The gain is the fast processing speed of these networks at test
time. Critical to their success is 1) the availability of training data, 2)
sufficient computing resources, and 3) the generalization of the trained
system to new datasets. These impact the performance of deep learning models
when deployed to unseen data in a new lab not encountered in training. As with
general deep learning models, the exact expected biases of the black box
system are unknown.

##  3 10 labs, 11 opinions: the out-of-distribution problem.

One of the most frustrating situations is the inability of an algorithm that
worked fantastically on one dataset but failed utterly on another. In
neuroimaging—wherein surgical protocols, custom microscopes, and variability
in environment can impact signal statistics from lab-to-lab—this effect is
further frustrated when a method thoroughly validated on many datasets from
one lab do not perform to the same degree in another lab. One reason can be
traced back to an insidious type of overfitting to the distribution of one
lab’s data, while even a similar lab might have data that may as well be _out
of distribution (OOD)_.

OOD is a central challenge in machine learning more generally [23, 24]. In
particular, the fact that much of modern-day machine learning focuses on
deriving much of what can be derived from the data itself, a central axiom is
that the full dataset follows a static distribution. When, of the full
distribution (say all videos of calcium imaging of GCaMP-6f transgenic mice in
Hipocampus area CA1), only a biased subset is used to train the ML algorithm
(say all videos of calcium imaging of GCaMP-6f transgenic mice in Hipocampus
area CA1 taken at Johns Hopkins University) then the extra assumption is
needed that the subset of the data follows the same distribution of the full
data. If it does not then a different subset (say all videos of calcium
imaging of GCaMP-6f transgenic mice in Hippocampus area CA1 taken at UC San
Diego) will fail to leverage what the ML algorithm learned from the first
subset. This effect is being continuously studied in more complex, modern, ML
systems, however will continue to have an impact in how we can expect ML
solutions to generalize across users.

One direct consequence of the “every lab is OOD” problem is that algorithms
that learn too much from data can suffer drawbacks in dissemination. An
obvious case is the use of deep learning. Deep learning fits highly over-
parameterized models and is known to be able to interpolate within a
distribution extremely well, while extrapolate very poorly to new domains.
Thus, despite the fact that pre-trained networks are fast to run, OOD create a
high level of uncertainty in using pre-trained networks accurately. In
developing methods for functional imaging analysis, and especially for
training deep learning models during method development, most use NeuroFinder
[25] and/or the Allen Institute Mouse Brain Observatory [26]. However, both
datasets suffer from both mislabeled data and missing labels, or erroneous
time-trace estimates [27]. Thus limited and mislabeled training data creates a
limited and potentially distorted _in distribution_ that the deep learning
model learns, which can affect performance in unexpected ways. Not all is
lost, however, as fine-tuning models on new datasets can be a work-around for
new data sets, albeit at the cost of further annotation or data collection. In
fact, many labs have manually labeled data acquired from their own optical
systems and in the brain regions they typically study, this barrier can be
somewhat mitigated provided easy processes by which to use this data to fine-
tune the model.

Thus, even though pre-trained networks are fast to run, these challenges
create a high level of uncertainty in using pre-trained networks accurately.
Individual labs can instead choose to annotate their own data to train
networks from scratch, such that they are optimized for the imaging used
locally. While this solves some of the aforementioned challenges (assuming
care is taken in annotation), the training procedure can be computationally
intensive relative to a lab’s computational capabilities.

Some methods explicitly seek to transcend such challenges by providing models
that learn from data in an instance-specific setting. For example, matrix
factorization and other such methods for source identification use more
general assumptions of low-rank and non-negativity rather than learning full
feature spaces as in deep learning. Danger, however, lies here as well. These
methods also have multiple (numerous) parameters that are often fit in certain
domains, or have regularizations that are fixed to work well for certain
imaging contexts. For example, algorithms that constrain components to be
contained in a small area in the image work excellently for somatic imaging
across brain areas and species but fail for different imaging scales and non-
compact or non-continuous morphologies (e.g., dendrites or widefield).

##  4 One pipeline does not fit all.

The generalization constraints discussed create a fundamental foil in the
field’s goal of creating highly disseminable and reproducible data pipelines
that can work “out of the box”. Such a goal is highly appealing, as we should
seek as much as possible to leverage the advances of others to optimize our
data analysis. To realize this goal, however, two challenges must be met:
robust assessment, and adaptive retraining.

The first of these challenges must face the fact that OOD means that the
assessment performed of the pipeline and its constituent parts might not be
correct for all labs/users. To enable new users to correctly choose an optimal
computational approach, new approaches to validation must be developed. These
validation approaches must promote data diversity and scale, rather than
dissecting singular datasets. Promising avenues in this direction include the
growing repositories of data in, e.g., the Dandi archive, and new simulation-
based data generators [28].

The second challenge lies in having the pipeline meet the users. Currently,
the goal of robust dissemination has been to make the code as immutable as
possible [29]. In some dimensions, this is the correct requirement, as one
does not want the code itself to change. In others, the pipeline must be
adaptive to identify when the approach is failing at processing a dataset and
which modification must be taken to improve the approach. Simple versions of
such methods would be simple parameter re-fitting (e.g., using BayesOpt) or
fine-tuning of deep learning models [30]. More complex solutions should be
able to swap and reorder stages in the pipeline. None of these is currently
“out of the box” for users and require significant manual intervention to
adapt pipelines.

##  5 Conclusions

The pace of neuroimaging analysis is rapid and unrelenting. Progress is often
achieved within the context of specific labs tackling problems which, while
critical to their work, would be invaluable to the community at large. Robust
dissemination is thus key and must account for the OOD problem.

New efforts offer promise. Efforts at the Allen Institute, the International
Brain Lab, DataJoint, and NeuroCAAS, amongst others, are striving to create
stable pipelines may lead to easier access. These efforts can form central
points for tools that assess and correct OOD problems as new users build up
their neuroimaging needs.

###### Acknowledgements.

G.M. was partially funded by a Kavli Institute for Brain and Mind (KIBM)
Innovative Research Grant award. We would like to thank Michael Xie for
helpful feedback.

## References

  * [1] Benisty, H., Song, A., Mishne, G., and Charles, A. S., “Review of data processing of functional optical microscopy for neuroscience,” Neurophotonics 9(4), 041402–041402 (2022).
  * [2] Charles, A. S., Falk, B., Turner, N., Pereira, T. D., Tward, D., Pedigo, B. D., Chung, J., Burns, R., Ghosh, S. S., Kebschull, J. M., et al., “Toward community-driven big open brain science: open big data and tools for structure, function, and genetics,” Annual review of neuroscience 43, 441–464 (2020).
  * [3] Lecoq, J., Oliver, M., Siegle, J. H., Orlova, N., Ledochowitsch, P., and Koch, C., “Removing independent noise in systems neuroscience data using deepinterpolation,” Nature Methods 18(11), 1401–1408 (2021).
  * [4] Li, X., Zhang, G., Wu, J., Zhang, Y., Zhao, Z., Lin, X., Qiao, H., Xie, H., Wang, H., Fang, L., et al., “Reinforcing neuron extraction and spike inference in calcium imaging using deep self-supervised denoising,” Nature Methods 18(11), 1395–1400 (2021).
  * [5] Weigert, M., Schmidt, U., Boothe, T., Müller, A., Dibrov, A., Jain, A., Wilhelm, B., Schmidt, D., Broaddus, C., Culley, S., et al., “Content-aware image restoration: pushing the limits of fluorescence microscopy,” Nature methods 15(12), 1090–1097 (2018).
  * [6] Eom, M., Han, S., Park, P., Kim, G., Cho, E.-S., Sim, J., Lee, K.-H., Kim, S., Tian, H., Böhm, U. L., et al., “Statistically unbiased prediction enables accurate denoising of voltage imaging data,” Nature Methods 20(10), 1581–1592 (2023).
  * [7] Liu, C., Platisa, J., Ye, X., Ahrens, A. M., Chen, I. A., Davison, I. G., Pieribone, V. A., Chen, J. L., and Tian, L., “Two-photon voltage imaging denoising by self-supervised learning,” in [Neural Imaging and Sensing 2023 ], 12365, 13–14, SPIE (2023).
  * [8] Platisa, J., Ye, X., Ahrens, A. M., Liu, C., Chen, I. A., Davison, I. G., Tian, L., Pieribone, V. A., and Chen, J. L., “High-speed low-light in vivo two-photon voltage imaging of large neuronal populations,” Nature Methods , 1–9 (2023).
  * [9] Wei, Z., Lin, B.-J., Chen, T.-W., Daie, K., Svoboda, K., and Druckmann, S., “A comparison of neuronal population dynamics measured with calcium imaging and electrophysiology,” PLoS computational biology 16(9), e1008198 (2020).
  * [10] Xu, Y. K. T., Graves, A. R., Coste, G. I., Huganir, R. L., Bergles, D. E., Charles, A. S., and Sulam, J., “Cross-modality supervised image restoration enables nanoscale tracking of synaptic plasticity in living mice,” Nature Methods , 1–10 (2023).
  * [11] Pachitariu, M., Packer, A. M., Pettit, N., Dalgleish, H., Hausser, M., and Sahani, M., “Suite2p: beyond 10,000 neurons with standard two-photon microscopy,” bioRxiv (2016).
  * [12] Giovannucci, A., Friedrich, J., Gunn, P., Kalfon, J., Brown, B. L., Koay, S. A., Taxidis, J., Najafi, F., Gauthier, J. L., Zhou, P., Khakh, B. S., Tank, D. W., Chklovskii1, D. B., and Pnevmatikakis, E. A., “CaImAn an open source tool for scalable calcium imaging data analysis,” Elife 8, e38173 (2019).
  * [13] Mishne, G. and Charles, A. S., “Learning spatially-correlated temporal dictionaries for calcium imaging,” in [IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP-2019) ], 1065–1069, IEEE (2019).
  * [14] Charles, A. S., Cermak, N., Affan, R. O., Scott, B. B., Schiller, J., and Mishne, G., “Graft: Graph filtered temporal dictionary learning for functional neural imaging,” IEEE Transactions on Image Processing 31, 3509–3524 (2022).
  * [15] Mukamel, E. A., Nimmerjahn, A., and Schnitzer, M. J., “Automated analysis of cellular signals from large-scale calcium imaging data,” Neuron 63(6), 747–760 (2009).
  * [16] Klibisz, A., Rose, D., Eicholtz, M., Blundon, J., and Zakharenko, S., “Fast, simple calcium imaging segmentation with fully convolutional networks,” in [Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support ], 285–293, Springer (2017).
  * [17] Ronneberger, O., Fischer, P., and Brox, T., “U-net: Convolutional networks for biomedical image segmentation,” in [International Conference on Medical image computing and computer-assisted intervention ], 234–241, Springer (2015).
  * [18] Apthorpe, N., Riordan, A., Aguilar, R., Homann, J., Gu, Y., Tank, D., and Seung, H. S., “Automatic neuron detection in calcium imaging data using convolutional networks,” Advances in neural information processing systems 29, 3270–3278 (2016).
  * [19] Soltanian-Zadeh, S., Sahingur, K., Blau, S., Gong, Y., and Farsiu, S., “Fast and robust active neuron segmentation in two-photon calcium imaging using spatiotemporal deep learning,” Proceedings of the National Academy of Sciences 116(17), 8554–8563 (2019).
  * [20] Bao, Y., Soltanian-Zadeh, S., Farsiu, S., and Gong, Y., “Segmentation of neurons from fluorescence calcium recordings beyond real time,” Nature Machine Intelligence , 1–11 (2021).
  * [21] Kirschbaum, E., Bailoni, A., and Hamprecht, F. A., “Disco: Deep learning, instance segmentation, and correlations for cell segmentation in calcium imaging,” in [International Conference on Medical Image Computing and Computer-Assisted Intervention ], 151–162, Springer (2020).
  * [22] Stringer, C., Wang, T., Michaelos, M., and Pachitariu, M., “Cellpose: a generalist algorithm for cellular segmentation,” Nature methods 18(1), 100–106 (2021).
  * [23] Geisa, A., Mehta, R., Helm, H. S., Dey, J., Eaton, E., Dick, J., Priebe, C. E., and Vogelstein, J. T., “Towards a theory of out-of-distribution learning,” arXiv preprint arXiv:2109.14501 (2021).
  * [24] Fort, S., Ren, J., and Lakshminarayanan, B., “Exploring the limits of out-of-distribution detection,” Advances in Neural Information Processing Systems 34, 7068–7081 (2021).
  * [25] Berens, P., Theis, L., Stone, J. S., Sofroniew, N., Tolias, A., Bethge, M., and Freeman, J., “Standardizing and benchmarking data analysis for calcium imaging,” in [COSYNE ], (2017).
  * [26] “Allen institute brain observatory whitepaper.” https://help.brain-map.org/display/observatory/Documentation.  Accessed: 2019-09-27.
  * [27] Gauthier, J. L., Koay, S. A., Nieh, E. H., Tank, D. W., Pillow, J. W., and Charles, A. S., “Detecting and correcting false transients in calcium imaging,” bioRxiv , 473470 (2018).
  * [28] Song, A., Gauthier, J. L., Pillow, J. W., Tank, D. W., and Charles, A. S., “Neural anatomy and optical microscopy (naomi) simulation for evaluating calcium imaging methods,” Journal of Neuroscience Methods 358, 109173 (2021).
  * [29] Abe, T., Kinsella, I., Saxena, S., Buchanan, E. K., Couto, J., Briggs, J., Kitt, S. L., Glassman, R., Zhou, J., Paninski, L., et al., “Neuroscience cloud analysis as a service: An open-source platform for scalable, reproducible data analysis,” Neuron 110(17), 2771–2789 (2022).
  * [30] Pachitariu, M. and Stringer, C., “Cellpose 2.0: how to train your own model,” Nature methods 19(12), 1634–1641 (2022).