\section{Introduction}








Many computational problems from engineering and science can be cast as certain parametric approximation problems (cf., e.g., \cite{Cohen2015,Khoo2021,Kutyniok2021,Bhattacharya2020,hs99,berner2020numerically,Vidales18,Chkifa15,Villa22,MR4153855} and references mentioned therein). 
In particular, parametric PDEs are of fundamental importance in various applications, where one is not only interested in an approximation of the solution of the approximation problem at one fixed (space-time) point but where one is interested to evaluate the approximative solution again and again as, for instance, in financial engineering where prices of financial products are computed approximately many times each trading day with (slightly) different parameters in each calculation. 
Moreover, the problems appearing in such financial applications are often high-dimensional, as the dimension usually corresponds to the number of assets/financial contracts in the considered trading portfolio. 


It is a widespread issue that the majority of algorithms for such parametric approximation problems suffer from the \textit{curse of dimensionality} (cf., e.g., Bellman \cite{Bellman1957} and Novak \& Wo\'{z}niakowski \cite[Chapter 1]{Novak08}) in the sense that the computational effort of the approximation methods grows exponentially in the dimension of the approximation problem or in the required approximation accuracy, making them useless for high-dimensional problems.
In the \textit{information based complexity} literature there are fundamental lower bounds which generally reveal the impossibility to approximate the solutions of certain classes of high-dimensional approximation problems without the curse of dimensionality among general classes of approximation algorithms; see, e.g.,
Grohs \& Voigtlaender \cite{Grohs2021},
Heinrich \cite{Heinrich2006},
Heinrich \& Sindambiwe \cite{hs99}, and
Novak \& Wo\'{z}niakowski \cite{NovakWozniakowski08}.
Developing methods which produce good approximations for high-dimensional problems is thus an exceedingly hard task and, among the general classes of algorithms considered in the above named references, essentially impossible.



In this paper we present a new method to tackle high-dimensional parametric approximation problems.
Roughly speaking, our strategy is based on the idea to \textit{combine Monte Carlo} (MC) \textit{algorithms} (such as, e.g., standard MC methods, multilevel Monte Carlo (MLMC) methods, or multilevel Picard (MLP) methods) \textit{with stochastic gradient descent} (SGD) \textit{optimization methods} by viewing the employed realizations of random variables in the MC approximation as training parameters for the SGD optimization method. 
In other words, in this approach we intend to employ SGD optimization methods not to train standard artificial neural networks (ANNs) but to learn random variables appearing in MC approximations.

To make this idea more concrete, we now sketch this \textit{learning the random variables} (LRV) approximation strategy in the context of a basic example of a parametric approximation problem.
Let $\mathfrak{p} \in \mathbb{N}$, suppose that we intend to approximate a target function $u \colon [0,1]^\mathfrak{p} \to \mathbb{R}$, 
let $\phi \colon [0,1]^\mathfrak{p} \times \mathbb{R} \to \mathbb{R} $ be measurable and bounded, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let $ \mathcal{W} \colon \Omega \to \mathbb{R} $ be a random variable, 
and suppose that $u \colon [0,1]^\mathfrak{p} \to \mathbb{R}$ admits the probabilistic representation that for all 
	$p \in [0, 1]^\mathfrak{p}$
we have that
\begin{equation}
\label{intro:eq0}
\begin{split} 
  u( p ) = 
  \mathbb{E}\!\left[ 
    \phi( 
      p, \mathcal{W}
    )
  \right]
\end{split}
\end{equation}
(parametric integration problem; cf., e.g.,
Cohen \& DeVore \cite{Cohen2015} and
Heinrich \& Sindambiwe \cite{hs99}). 
We  note that we only chose the random variable $\mathcal{W}$ to be 1-dimensional for simplicity and refer to \cref{sect:MC} for the case when $\mathcal{W}$ is a possibly high-dimensional random variable.
Our first step to derive the proposed approximation algorithm is to recall standard MC approximations for the parametric expectation in \eqref{intro:eq0}. 
Let $ \mathfrak{M} \in \mathbb{N} $, 
let 
$ 
  W^{ m, \mathfrak{m} } \colon \Omega \to \mathbb{R} 
$, $ m, \mathfrak{m} \in \mathbb{N}_0 $, 
be i.i.d.\ random variables which satisfy for all $ A \in \mathcal{B}( \mathbb{R} ) $ that 
$
  \P( W^{ 0, 0 } \in A ) = \P( \mathcal{W} \in A ) 
$, and observe 
for all $ p \in [0,1]^\mathfrak{p} $ that
\begin{equation}
\label{intro:eq1}
  u( p ) = 
  \mathbb{E}\!\left[ 
    \phi( 
      p, \mathcal{W}
    )
  \right]
  =
  \mathbb{E}\!\left[ 
    \phi( 
      p, W^{ 0, 1 }
    )
  \right]
  \approx 
  \tfrac{ 1 }{ \mathfrak{M} }
  \big[
    {\textstyle \sum\nolimits_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }}
    \phi( 
      p, W^{ 0, \mathfrak{m} }
    )
  \big]
  .
\end{equation}
In the next step we introduce a parametric function on $[0,1]^\mathfrak{p}$ with parameter set $\mathbb{R}^{ \mathfrak{M}}$ to reformulate the standard MC approximation in \eqref{intro:eq1} in such a way that the random variables 
$
	W^{ 0, m } \colon \Omega \to \mathbb{R}
$, $m \in \{1, 2, \ldots, \mathfrak{M}\}$, 
correspond to the parameters of the parametric function.
More precisely, let $\mathscr{N} \colon [0,1]^\mathfrak{p} \times \mathbb{R}^{ \mathfrak{M}} \to \mathbb{R}$ satisfy for all
	$p \in [0,1]^\mathfrak{p}$, 
	$\theta = ( \theta_1, \dots, \theta_{ \mathfrak{M} } ) \in \mathbb{R}^{\mathfrak{M}}$
that
$
	\mathscr{N}(p, \theta)
=
    \frac{ 1 }{ \mathfrak{M} }
      \big[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \phi( 
          p, 
          \theta_{ \mathfrak{m} }
        )
      \big]
$
and note that \eqref{intro:eq1} suggests for all 
	$ p \in [0,1]^\mathfrak{p} $ 
that
\begin{equation}
\label{intro:eq3}
  u( p )
  \approx 
  \tfrac{ 1 }{ \mathfrak{M} }
  \big[
    {\textstyle \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
    \phi( 
      p, W^{ 0, \mathfrak{m} }
    )}
  \big]
 =
 \mathscr{N}\!\big(p,  (W^{ 0, 1 }, W^{ 0, 2 }, \ldots, W^{ 0,  \mathfrak{M}})\big)
  .
\end{equation}
Next we employ the SGD optimization method to train the right hand side of \eqref{intro:eq3} in search of “better random realizations” to approximate the target function $u$ 
than those provided by the random variables 
$
	W^{ 0, m } \colon \Omega \to \mathbb{R}
$, $m \in \{1, 2, \ldots, \mathfrak{M}\}$, 
in the standard MC approximation in \eqref{intro:eq3}.
The random variables 
$
	W^{ 0, m } \colon \Omega \to \mathbb{R}
$, $m \in \{1, 2, \ldots, \mathfrak{M}\}$, 
on the right hand side of \eqref{intro:eq3} then only supply the initial guess in the SGD training procedure. 
More specifically, 
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to [0,1]^\mathfrak{p}
$, $ m, {\bf m} \in \mathbb{N} $, 
be independent continuous uniformly distributed random variables, 
let $ {\bf M} \in \mathbb{N} $, 
for every 
	$m \in \mathbb{N}$
let 
$ 
  F_m 
  \colon \mathbb{R}^{ \mathfrak{M} } 
  \times \Omega \to \mathbb{R} 
$
satisfy 
for all 
$ 
  \theta = 
  ( \theta_1, \dots, \theta_{ \mathfrak{M} } ) 
  \in \mathbb{R}^{ \mathfrak{M} } 
$ 
that
\begin{equation}
\label{intro:eq4}
\begin{split}
  F_m( \theta ) 
  &=
  \tfrac{ 1 }{ {\bf M} }
  \big[
    {\textstyle \sum_{ {\bf m} = 1 }^{ {\bf M} }
    \left|
      \phi( 
        P_{ 
          m, {\bf m} 
        }, 
        W^{ m, {\bf m} } 
      )
      -
      \mathscr{N}(
        P_{ 
          m, {\bf m} 
        }, 
        \theta
      )
    \right|^2
    }
  \big],
\end{split}
\end{equation}
assume for all 
$p \in [0,1]^\mathfrak{p}$ that
$
	(\mathbb{R} \ni w \mapsto \phi(p, w) \in \mathbb{R}) \in C^1(\mathbb{R}, \mathbb{R})
$,
for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{ \mathfrak{M} } \times \Omega 
  \to \mathbb{R}^{ \mathfrak{M} } 
$
satisfy 
for all 
	$ \theta \in \mathbb{R}^{ \mathfrak{M} } $, $ \omega \in \Omega $
that
$
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega )
$,
let $ (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) $, and 
let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \mathfrak{M} } 
$ 
satisfy 
for all $ m \in \mathbb{N} $ that 
\begin{equation}
\label{intro:eq5}
  \Theta_0 = 
  ( W^{ 0, 1 }, W^{ 0, 2 }, \dots, W^{ 0, \mathfrak{M} } )
 \qquad\text{and}\qquad
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
\end{equation}
Note that the recursion in \eqref{intro:eq5} describes nothing else but the standard SGD optimization method with the learning rate schedule
$
\mathbb{N} \ni m \mapsto
\gamma_m \in \mathbb{R}
$.
For every sufficiently large 
$m \in \mathbb{N}$
we then propose to employ the random function 
\begin{equation}
\label{intro:eq6}
\begin{split} 
	[0,1]^\mathfrak{p} \times \Omega \ni (p, \omega) \mapsto \mathscr{N}(p, \Theta_m(\omega)) \in \mathbb{R}
\end{split}
\end{equation}
as an approximation for the target function
$
	[0,1]^\mathfrak{p} \ni p \mapsto u(p) \in \mathbb{R}
$. 
Note that the set 
$
	[0,1]^\mathfrak{p} \ni p \mapsto \mathscr{N}(p, \theta) \in \mathbb{R}
$, $\theta \in \mathbb{R}^{\mathfrak{M}}$,
of potential approximating functions for $u \colon [0,1]^\mathfrak{p} \to \mathbb{R}$
does not consist of standard fully-connected feedforward ANNs but
is a very problem specific class of approximating functions, which we refer to as \textit{MC neural networks},
determined by the MC method in \eqref{intro:eq1}.
Moreover, observe that through \eqref{intro:eq3} the MC method also naturally specifies a favorable initializing law for the SGD training procedure in \eqref{intro:eq5}.

In the special case of the LRV strategy illustrated in \eqref{intro:eq1}--\eqref{intro:eq6} above, the standard MC approximation method serves as a proposal algorithm for the SGD training procedure. 
More generally, the LRV strategy can, in principle, be used on any high-dimensional approximation problem as soon as there is a reasonable stochastic proposal algorithm for the considered approximation problem available.
The LRV strategy thereby naturally specifies the compositional
architecture of the involved approximating functions and also naturally specifies the initializing law in the SGD
training procedure for each specific approximation problem.
In particular, the LRV strategy can be used 
	with the standard MC method (see \cref{sect:MC}),
	the MC-Euler-Maruyama method (see \cref{sect:MCEUler}), or 
	the MLMC method (see \cref{sect:MLMC}) 
as the proposal algorithms to approximate solutions of stochastic differential equations (SDEs) and Kolmogorov PDEs, respectively, and 
the LRV strategy can be used with 
	the MLP method (see \cref{sect:MLP}) 
as the proposal algorithm to approximate solutions of semilinear PDEs.




\begin{table}[htb] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|,
	separator=semicolon,
     table head=
     \hline Approximation method &\thead{Number of \\ trainable \\ parameters} & \thead{Number\\ of \\ MC/QMC \\ samples} & \thead{$L^2$-approx. \\ error}& \thead{$L^\infty$-approx. \\ error}& \thead{Training \\ time in \\ seconds} & \thead{Evaluation \\ time in \\ seconds}  \\\hline,
    late after line=\\\hline]{Table_1.csv}
{
	method = \method, 
	num-params = \nparams,
	num-samples = \nsamp,
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	train-time=\train, 
	eval-time=\eval
}
{ \method  & \nparams & \nsamp & \llll & \linf &\train & \eval }}
\caption{\label{table:BS1_Summary}Pricing European call options in the Black-Scholes model}
\end{table}




We now illustrate the effectiveness of the LRV strategy on a simple but famous numerical example, where the LRV strategy leads to impressively precise approximations.
In \cref{table:BS1_Summary} we present results for the approximative computation of prices of European call options in the Black-Scholes model 
	by means of the deep learning method induced by Becker et al.\ \cite{BeckJafaari21} with $140000$ Adam (see Kingma \& Ba \cite{Kingma2014}) training steps
		(rows 2--5 in \cref{table:BS1_Summary}),
	by means of the standard MC method with 32768 MC samples 
		(row 6 in \cref{table:BS1_Summary}), 
	by means of the antithetic MC method with 32768 MC samples 
			(row 7 in \cref{table:BS1_Summary}), 
	by means of the Quasi-Monte Carlo (QMC) method using Sobol sequences with 32768 QMC samples 
		(row 8 in \cref{table:BS1_Summary}), and
	by means of the antithetic QMC method using Sobol sequences with 32768 QMC samples 
		(row 9 in \cref{table:BS1_Summary}), 
	by means of the LRV strategy with the standard MC method as the proposal algorithm with $140000$ Adam training steps
		(row 10 in \cref{table:BS1_Summary}), and
	by means of the LRV strategy with the antithetic MC method as the proposal algorithm with $140000$ Adam training steps
		(row 11 in \cref{table:BS1_Summary}).
In \cref{table:BS1_Summary} the $L^2$-error and the $L^\infty$-error have been computed approximately on the region
\begin{equation}
\begin{split} 
\label{intro:eq7}
	(\xi,T,r,\sigma,K) 
\in 
	[90,110]
	\times[ \tfrac{1}{100} ,1]
	\times[- \tfrac{1}{10} , \tfrac{1}{10} ]
	\times[ \tfrac{1}{100} , \tfrac{1}{2}]
	\times[90,100]
\end{split}
\end{equation}
using 8\,192\,000 evaluations of the considered approximation method. 
In \eqref{intro:eq7} we have 
	that $\xi$ stands for the initial price, 
	that $T$ stands for the time of maturity, 
	that $r$ stands for the drift rate, 
	that $\sigma$ stands for the volatility, and
	that $K$ stands for the strike price; see, e.g., \cite[Lemma 4.4]{Becker2019published} (with c = 0 in the notation of \cite[Lemma 4.4]{Becker2019published}). 
The reference solution values to compute the errors in \cref{table:BS1_Summary} have been computed with the famous Black-Scholes formula (see, e.g., \cite[Lemma 4.4]{Becker2019published} or \eqref{BS1:eq1} in \cref{simul:BS1}). 
The training time and evaluation time\footnote{The numerical experiments have been performed in {\sc Python} using {\sc TensorFlow} on a {\sc NVIDIA GeForce RTX 2080 Ti GPU}. The operating system was an AMD Ryzen 9 3950X CPU with 64 GB DDR4 memory running {\sc TensorFlow} 2.1 on {\sc Ubuntu} 20.04.}
are presented in seconds in \cref{table:BS1_Summary}. The evaluation time corresponds to the time required to compute 8\,192\,000 evaluations.
The numbers in \cref{table:BS1_Summary} are taken from \cref{table:BS_LRV,table:BS_MC,table:BS_ANNs,table:BS_QMC} in \cref{simul:BS1}. 
We refer to \cref{simul:BS1} for more details on the results in \cref{table:BS1_Summary}.


Note that \cref{table:BS1_Summary} indicates that the algorithm obtained by the LRV strategy not only produces very accurate prices in the $L^2$-norm, but even has a very high accuracy in the uniform $L^\infty$-norm.
Concretely, for this 5-dimensional approximation problem, this strongly suggests that for any choice of parameters in the region considered in \eqref{intro:eq7}, the LRV strategy offers an approximation with an error smaller than $\frac{1}{1000}$.
Based on 
	this, 
	on the other numerical results presented in this paper, as well as
	on preliminary analytic investigations (cf., e.g., Gonon et al.\ \cite[Lemma 2.16]{Gonon19Uniform})
	we conjecture that the LRV strategy can overcome the curse of dimensionality for certain classes of parametric PDE problems in the $L^\infty$-norm. 
We would like to emphasize that large classes of algorithms for such approximation problems have been shown not to be able to overcome the curse of dimensionality, see, e.g., Heinrich \cite[Theorem 2.4]{hs99}, Heinrich \& Sindambiwe \cite[Theorem 1]{Heinrich2006}, and Grohs \& Voigtlaender \cite{Grohs2021}.
However our conjecture does not contradict to the general lower bounds  established in the
above mentioned references, due to the fact that the LRV strategy does not belong to the class of algorithms considered in the
above mentioned references: roughly speaking, in the LRV strategy there are two stages of computational procedures, the \textit{main computational procedure} in which the “best random variables” are learned through SGD (corresponding to the 6th column in \cref{table:BS1_Summary}) and the \textit{evaluation procedure} where the computed approximation of the target function is evaluated (corresponding to the 7th column in \cref{table:BS1_Summary}). In the LRV strategy we consider the situation where it is allowed to perform function evaluations both during the {main computational} procedure \textit{and} the {evaluation procedure} while the lower bounds in \cite{hs99,Heinrich2006,Grohs2021} consider the situation where function evaluations are only allowed during the main computational procedure \textit{but not} during the evaluation procedure. 
The LRV strategy being outside of the classes of algorithms considered in the above named references is thus so not constrained by the established lower bounds and hence holds the potential to overcome the curse of dimensionality, even in the $L^\infty$-norm.
We note that in practically relevant situations, just as in the considered derivative pricing problem, it is often possible to perform function evaluations also in the evaluation procedure, thus making the LRV strategy an applicable method for practically relevant approximation problems.





We now compare the proposed LRV strategy to existing algorithms and computational methods in the scientific literature.
As the LRV strategy employs SGD-type method to "learn" parametric functions it is related to deep learning methods, where instead of random variables, optimal weights of ANNs are "learned".
From this point of view, the LRV strategy can be seen as a machine learning approach where instead of employing generic ANNs, very problem specific parametric functions are used, which contain a lot of human insight about the problem at hand.
There is a plethora of deep learning methods for the approximation of PDEs, which have been developed recently, and seem to be very effective for the approximation of high-dimensional PDEs:
cf., e.g., 
\cite{Weinan2017,han2018solving,sirignano2018dgm,E17Ritz,Fujii19Asymptotic,dissanayake1994neural,lagaris1998artificial,jianyu2003numerical,chen2020friedrichs,ChanMikaelWarin2019,Germain2020,Han20,han2020solving,henry2017deep,Hure20,Jacquier2019a,Shaolin20,Kremsner2020,WangEtal2019,pham2019neural,raissi2018forward,DeepSplitting,BeckJafaari21,berner2020numerically,HanEtal2020,Nusken21,samaniego2020energy,Becker2018,Becker2019published,Chen2019,Goudenege2019,Lyu2020,Lu2019,Zang2020weak,Han2019deep,Lye20,Zhu19,LongoEtal2020,Nuesken2021,Martin2020}. We refer to the survey articles \cite{beck2020overview,blechschmidt2021three,Germain2021,Weinan2020a} for a more detailed overview.
For methods which are specifically designed for parametric PDEs we refer to, e.g., \cite{Vidales18,Khoo2021,Bhattacharya2020,berner2020numerically,Villa22}.
In addition, we want to highlight a connection between the LRV strategy and the deep learning method for Kolmogorov PDEs developed in Becker et al.\ \cite{BeckJafaari21} and further specialized to parametric Kolmogorov PDEs in Berner et al.\ \cite{berner2020numerically}. Very roughly speaking, in case of parametric Kolmogorov PDEs the LRV strategy corresponds to the algorithm in \cite{berner2020numerically} if the ANNs appearing in the algorithm in \cite{berner2020numerically} are replaced by parametric functions induced by the proposal algorithm (typically given by the MC algorithm through the Feynman-Kac formula).


Furthermore, the LRV strategy can be associated with a broad subcategory of machine learning called \textit{algorithm learning} (cf., e.g., \cite{Towell1994,Chen2020,Chen2020RNA,gregor2010learning,Borgerding2016,Chen2018,Yoon2013,mensch18a,Wilder2019}).
Roughly speaking, algorithm learning refers to the idea to employ existing algorithms with known empirical or theoretical qualities as a basis to construct or extend ANNs or more general parametric function families.
In many cases the employed algorithm relies on certain hyper-parameters 
(such as, e.g., 
the learning rates in case of SGD-type methods: cf., e.g., Chen et al.\ \cite{Chen2020})
which typically are added to the set of trainable parameters of the ANN.
It is in this point that the LRV strategy differs from existing algorithm learning methods since the LRV strategy considers the random variables and not the hyper-parameters of the proposal algorithm as learnable parameters.
It thereby has the advantage that the initialization of all trainable parameters is implicitly given through the proposal algorithm.
However, the LRV strategy could very well be combined with ideas of algorithm learning. We leave this task open for future research.


Another tranche of literature connected to the LRV strategy are QMC methods (cf., e.g., \cite{drovandi2018improving,Morokoff1995,Caflisch1998}).
Roughly speaking, the idea of QMC methods is to replace uniformely distributed random variables of the MC method by a more suitable sequence of deterministic points to obtain a higher rate of convergence.
The idea to improve the choice of random variables in MC methods is a common feature of QMC methods with the LRV strategy. A key difference between the two methods is that QMC methods construct new integration points which have good properties for a wide class of integrands while the LRV strategy aims to "learn" new integration points which are specific to each considered integrand.



Next we illustrate a link between the LRV strategy and quantization methods (cf., e.g., \cite{Pages98,pages2015introduction,pages2004optimal,pages2003optimal,pages2005functional,dereich2013constructive,frikha2012quantization,pham2005approximation,rudd2017fast,Gronbach13Quantization,Altmayer14}). 
Quantization methods are concerned with approximating a continuously distributed random variable by a random variable with a finite image.
Quantization was first introduced for signal processing and information theory (cf., e.g., \cite{Graf20Foundations,pages2003optimal}) but can also be used for the numerical approximation of expectations involving the original random variable (cf., e.g., \cite{Pages98,pages2003optimal}). 
For the latter, the original random variable in the expectation is replaced by its quantization, resulting in an expected value which can easily be computed by a finite sum of weighted function evaluations.
A good quantization
is typically found with optimization methods such as, e.g., the Newton method in low dimensions or SGD-type optimization methods in high-dimensions (cf., e.g., \cite{pages2003optimal,pages2015introduction,pages2004optimal,pham2005approximation}).
To make this more concrete we now roughly illustrate this in the context of the setting described above.
For all 
	$\theta = (\theta_1, \ldots, \theta_\mathfrak{M}) \in \mathbb{R}^\mathfrak{M}$, 
	$C_1, C_2, \ldots, C_\mathfrak{M} \in \mathcal{B}(\mathbb{R})$ 
with 
	$\forall \, m,n \in \{1, 2, \ldots, \mathfrak{M}\}, m \neq n \colon C_m \cap C_n = \{\}$ 
	and 
	$\cup_{m = 1}^\mathfrak{M} C_m = \mathbb{R}$ 
we consider a quantization
	$Q_{\theta,C_1, C_2, \ldots, C_\mathfrak{M}} \colon \Omega \to \mathbb{R}^\mathfrak{M}$ 
of $\mathcal{W}$ given by
$
	Q_{\theta,C_1, C_2, \ldots, C_\mathfrak{M}}
=
	\sum_{\mathfrak{m} = 1}^\mathfrak{M} \theta_{\mathfrak{m}} \mathbbm{1}_{C_\mathfrak{m}}(\mathcal{W})
$.
A good quantization is found by minimizing an error between the quantization and the original random variable such as, e.g., the squared $L^2$-error given for all 
	$\theta = (\theta_1, \ldots, \theta_\mathfrak{M}) \in \mathbb{R}^\mathfrak{M}$, 
	$C_1, C_2, \ldots, C_\mathfrak{M} \in \mathcal{B}(\mathbb{R})$ 
by
\begin{equation}
\label{intro:eq8}
\begin{split} 
\Exp{\Vert\mathcal{W} - Q_{\theta,C_1, C_2, \ldots, C_\mathfrak{M}}\Vert^2},
\end{split}
\end{equation}
with an SGD-type optimization method.
Once appropriate 
	$\theta = (\theta_1, \ldots, \theta_\mathfrak{M}) \in \mathbb{R}^\mathfrak{M}$, 
	$C_1, C_2, \ldots, C_\mathfrak{M} \in \mathcal{B}(\mathbb{R})$ 
	with
		$\forall \, m,n \in \{1, 2, \ldots, \mathfrak{M}\}, m \neq n \colon C_m \cap C_n = \{\}$ 
		and 
		$\cup_{m = 1}^\mathfrak{M} C_m = \mathbb{R}$ 
have been found,
they can be employed to approximate the expectation in \eqref{intro:eq0} for every $p \in \mathfrak{P}$ by
\begin{equation}
\label{intro:eq9}
\begin{split} 
\textstyle
	u(p) 
=
	\Exp{\phi(p, \mathcal{W})}
\approx
 	\Exp{\phi(p, Q_{\theta,C_1, C_2, \ldots, C_\mathfrak{M}})}
=
	\sum_{\mathfrak{m} = 1}^\mathfrak{M} \phi(p, \theta_{\mathfrak{m}}) \P( \mathcal{W} \in C_\mathfrak{m}).
\end{split}
\end{equation}
Note that this approximation has a similar form as the proposed LRV approximation in \eqref{intro:eq6} and
in both cases the points at which $\phi$ is evaluated are found through an SGD-type optimization method.
The main difference is that the optimization problem in \eqref{intro:eq8}, which is used to determine the evaluation points in \eqref{intro:eq9}, only depends on the random variable $\mathcal{W}$ whereas the optimization problem (see \eqref{intro:eq4}) to determine the evaluation points in \eqref{intro:eq6} in the LRV strategy depends on the random variable $\mathcal{W}$ and the function $\phi$.


The concept of optimizing random MC samples, a central aspect of the LRV strategy, also appears in the Bayesian statistics literature in the context of inducing points in Bayesian learning (cf., e.g., \cite{snelson2005sparse,ranganath2014black,dellaportas2019gradient,Betancourt2017}).
Such inducing points are employed to find sparse representation for large data sets.
One way to find good inducing points is to first sample them randomly and then optimize them for example with gradient based methods (cf., e.g., Snelson \& Ghahramani \cite{snelson2005sparse}), which bears some resemblance to the LRV strategy.



The reminder of this article is structured as follows.
In \cref{sect:MC,sect:MCEUler,sect:MLMC,sect:MLP} we present the LRV strategy for increasingly complex approximation problems and proposal algorithms. 
As proposal algorithms we consider 
	the standard MC method in \cref{sect:MC},
	the MC-Euler method in \cref{sect:MCEUler},
	the MLMC method in \cref{sect:MLMC}, and
	the MLP method in \cref{sect:MLP}.
In \cref{sect:general} we present the most general case of a generic proposal algorithm, which includes all previous sections as special cases.
The results of numerical experiments for the LRV approximation strategy are presented in \cref{sect:numerics}.
Specifically, 
	we consider 1-dimensional Black-Scholes equations for European call options (resulting in a 5-dimensional parametric approximation problem) in \cref{simul:BS1},
	we consider 3-dimensional Black-Scholes equations for worst-of basket put options (resulting in a 15-dimensional parametric approximation problem)  in \cref{simul:BSworstPut},
	we consider 3-dimensional Black-Scholes equations for average basket put options with knock-in barriers (resulting in a 16-dimensional parametric approximation problem) in \cref{simul:BSaverageBarrier},
	and
	we consider stochastic Lorentz equations (resulting in a 10-dimensional parametric approximation problem) in \cref{simul:Lorentz}.
Finally, in \cref{sect:source} we present a representative {\sc Python} code, which has been used in this paper to implement the LRV strategy. 















\section{Learning the random variables (LRV) strategy in the case of Monte Carlo approximations}
\label{sect:MC}

In this section we employ the LRV strategy for the approximation of parametric expectations (see \eqref{MC:eq1} in  \cref{MC:target}).
This includes as a special case the approximative pricing of European options in the Black-Scholes model, as illustrated in \cref{MC:example}.
The LRV strategy 
based on standard MC averages as proposal algorithms 
for the general parametric approximation problem of \cref{MC:target}
is elaborated in \cref{MC:proposal,MC:replacing,MC:Loss,MC:SGD}.
The resulting method is summarized in a single framework in \cref{MC:description}.



\subsection{Parametric expectations involving vector valued random variables}
\label{MC:target}


Let $ \mathfrak{p},\mathbf{d}  \in \mathbb{N} $, 
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let $ u \colon \mathfrak{P} \to \mathbb{R} $ be a function, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let $ \mathcal{W} \colon \Omega \to \mathbb{R}^{ \mathbf{d} } $ be a random variable, 
let $ \phi \colon \mathfrak{P} \times \mathbb{R}^{ \mathbf{d} } \to \mathbb{R} $ be measurable, 
assume for all $ p \in \mathfrak{P} $ that
$
  \mathbb{R}^{ \mathbf{d} } \ni w \mapsto \phi( p, w ) \in \mathbb{R} 
$
is continuously differentiable, 
and assume for all $ p \in \mathfrak{P} $ that 
$
  \mathbb{E}\big[
    | \phi( p, \mathcal{W} ) |
  \big]
  < \infty
$
and
\begin{equation}
\label{MC:eq1}
  u( p ) = 
  \mathbb{E}\!\left[ 
    \phi( 
      p, \mathcal{W}
    )
  \right]
  .
\end{equation}
The goal of this section is to derive an algorithm to approximately compute the function
$
	u \colon \mathfrak{P} \to \mathbb{R}
$
given through the parametric expectation in \eqref{MC:eq1}.

\subsection{Approximative pricing of European options in the Black-Scholes model}
\label{MC:example}


\begin{lemma}\label{cor:black_scholes}
Let
$ \xi, T, \sigma \in ( 0, \infty ) $,
$ r \in \mathbb{R} $,
let $ \mathfrak{N} \colon \mathbb{R} \to \mathbb{R} $ satisfy for all $ x \in \mathbb{R} $ that
$
\mathfrak{N}( x ) = \int_{ - \infty }^x \frac{ 1 }{ \sqrt{ 2 \pi } } \, \exp({ - \frac{ y^2 }{ 2 }  }) \d y
$,
let $ ( \Omega, \mathcal{F}, ( \mathbb{F}_t )_{ t \in [0,T] }, \P ) $ be a filtered probability space
which satisfies the usual conditions,
let $ W \colon [0,T] \times \Omega \to \mathbb{R} $
be a standard $( \mathbb{F}_t )_{ t \in [0,T] }$-Brownian motion
with continuous sample paths,
and
let $ X \colon [0,T] \times \Omega \to \mathbb{R} $
be an $( \mathbb{F}_t )_{ t \in [0,T] }$-adapted stochastic process with continuous sample paths
which satisfies that
for all $ t \in [ 0, T ] $
it holds $ \P $-a.s.\ that
\begin{equation}
X_t
= \xi +
r \int_{0}^{t} X_s \, \d s +
\sigma \int_{0}^{t} X_s \, \d W_s
.
\end{equation}
Then 
\begin{enumerate}[label=(\roman *)]
\item \label{cor:black_scholes:item1}
it holds for all $ t \in [ 0, T ] $ that
$
\P\big(X_t
=
\exp\bigl(
\bigl[ r - \tfrac{ \sigma^2 }{2} \bigr] t
+
\sigma \, W_t \bigr)
\xi
\big)
=
1
$
and
\item \label{cor:black_scholes:item2}
it holds
for all $ K \in \mathbb{R} $ that
\begin{equation}
\begin{split}
& \mathbb{E} \bigl[
\exp({ -r T })
\max \{ X_T - K, 0 \}
\bigr]
\\ & =
\begin{cases}
\xi
\,
\mathfrak{N} \Bigl(
\tfrac{ ( r + \frac{ \sigma^2 }{2} ) T + \ln( {\xi}/{K} ) }{ \sigma \sqrt{T} }
\Bigr)
-
K \exp({ - r \, T})
\,
\mathfrak{N} \Bigl(
\tfrac{ ( r - \frac{ \sigma^2 }{2} ) T + \ln( {\xi}/{K} ) }{ \sigma \sqrt{T} }
\Bigr)
&
\colon
K > 0
\\[1ex]
 \xi
- K \exp({ - r \, T})
&
\colon
K \leq 0.
\end{cases}
\end{split}
\end{equation}

\end{enumerate}


\end{lemma}
\begin{proof}[Proof of \cref{cor:black_scholes}]
Observe that, e.g.,  \cite[item (i) in Proposition~4.3]{Becker2019published} implies that
for all $ t \in [ 0, T ] $
it holds $ \P $-a.s.\ that
\begin{equation}
\label{cor:black_scholes:eq1}
X_t
=
\exp\bigl(
\bigl[ r - \tfrac{ \sigma^2 }{2} \bigr] t
+
\sigma \, W_t \bigr)
\xi
=
\exp\big({
 ( r - \tfrac{ \sigma^2 }{2} ) t
+
\ln( \xi )
+
\sigma W_t}\big)
.
\end{equation}
This establishes item~\ref{cor:black_scholes:item1}.
Moreover, note that, e.g., \cite[Lemma~4.4]{Becker2019published} and \eqref{cor:black_scholes:eq1}  show
for all $ K \in \mathbb{R} $ that
\begin{equation}
\begin{split}
& \mathbb{E} \bigl[
\exp({ -r T })
\max \{ X_T - K, 0 \}
\bigr]
=
\exp({ -r T }) \,
\mathbb{E} \Bigl[
\max \Bigl\{ \exp\big({
( r - \tfrac{ \sigma^2 }{2} ) T
+
\ln( \xi )
+
\sigma W_T }\big) - K, 0 \Bigr\}
\Bigr]
\\ & =
\begin{cases}
 \xi
\,
\mathfrak{N} \Bigl(
\tfrac{ ( r + \frac{ \sigma^2 }{2} ) T + \ln( {\xi}/{K} ) }{ \sigma \sqrt{T} }
\Bigr)
-
K \exp({ - r \, T})
\,
\mathfrak{N} \Bigl(
\tfrac{ ( r - \frac{ \sigma^2 }{2} ) T + \ln( {\xi}/{K} ) }{ \sigma \sqrt{T} }
\Bigr)
&
\colon
K > 0
\\[1ex]
\xi
- K \exp({ - r \, T})
&
\colon
K \leq 0
\end{cases}
.
\end{split}
\end{equation}
The proof of \cref{cor:black_scholes} is thus complete.
\end{proof}



In the case 
where 
  $ \mathfrak{p} = 5 $,
  $ \mathbf{d} = 1 $, 
   and
  $
    \mathfrak{P} = [90, 110] \times [0.01,1] \times [-0.1,0.1] \times [0.01,0.5] \times  [90, 110]
  $, 
where 
for all 
$
  p = (\xi, T, r, \sigma, K ) \in \mathfrak{P}
$,
$
  w \in \mathbb{R}
$
it holds
that
$
  \phi( p, w ) = 
  \exp({-r T})
  \max\!\big\{ \!
    \exp\!\big(
      [ r - \frac{ \sigma^2 }{ 2 } ] T + T^{ 1 / 2 } \sigma w  
    \big) \xi
    - K
    , 0
  \big\}
  $,
  where $ \mathcal{W} $ is a standard normal random variable, 
  and
  where $ \Phi \colon \mathbb{R} \to \mathbb{R} $ satisfies for all $ x \in \mathbb{R} $ that
$
\mathfrak{N}( x ) = \int_{ - \infty }^x \frac{ 1 }{ \sqrt{ 2 \pi } } \, \exp({ - \frac{ y^2 }{ 2 } }) \, \d y
$
observe that 
\begin{enumerate}[label=(\roman *)]
\item 
it holds for all $ p = (\xi, T, r, \sigma, K )  \in \mathfrak{P} $
that
\begin{equation}
  \begin{split}
    u(p) &= \mathbb{E} \bigl[
     \exp({-r T})
      \max\!\big\{ \!
        \exp\!\big(
          [ r - \tfrac{ \sigma^2 }{ 2 } ] T + T^{ 1 / 2 } \sigma \mathcal{W}  
        \big) \xi
        - K
        , 0
      \big\}
\bigr]
\\
  &=
  \xi
\,
\mathfrak{N} \Bigl(
\tfrac{ ( r  + \frac{ \sigma^2 }{2} ) T + \ln( {\xi}/{K} ) }{ \sigma \sqrt{T} }
\Bigr)
-
K \exp({ - r \, T})
\,
\mathfrak{N} \Bigl(
\tfrac{ ( r - \frac{ \sigma^2 }{2} ) T + \ln( {\xi}/{K} ) }{ \sigma \sqrt{T} }
\Bigr)
  \end{split}
\end{equation}
(cf.\ \eqref{MC:eq1} and \cref{cor:black_scholes})
and
\item 
it holds for all 
	$p = (\xi, T, r, \sigma, K ) \in \mathfrak{P}$ 
that
$u(p) \in [0,\infty)$ is
the price of an European call option 
in the Black-Scholes model 
with 
initial underlying price $ \xi $, 
time of maturity $ T $, 
drift rate $ r $, 
volatility $ \sigma $, 
and strike price $ K $.
\end{enumerate}





\subsection{Monte Carlo approximations}
\label{MC:proposal}


Let $ \mathfrak{M} \in \mathbb{N} $
and let 
$ 
  W^{ m, \mathfrak{m} } \colon \Omega \to \mathbb{R}^{ \mathbf{d} } 
$, 
$ m, \mathfrak{m} \in \mathbb{N}_0 $, 
be i.i.d.\ random variables 
which satisfy for all $ A \in \mathcal{B}( \mathbb{R}^{ \mathbf{d} } ) $ that 
$
  \P( W^{ 0, 0 } \in A ) = \P( \mathcal{W} \in A ) 
$.
Observe that \eqref{MC:eq1} suggests that 
for all $ p \in \mathfrak{P} $ it holds that
\begin{equation}
\label{mc:eq1}
  u( p ) = 
  \mathbb{E}\!\left[ 
    \phi( 
      p, \mathcal{W}
    )
  \right]
  =
  \mathbb{E}\!\left[ 
    \phi( 
      p, W^{ 0, 1 }
    )
  \right]
  \approx 
  \frac{ 1 }{ \mathfrak{M} }
  \left[
    \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
    \phi( 
      p, W^{ 0, \mathfrak{m} }
    )
  \right]
  .
\end{equation}

\subsection{Replacing the random variables in Monte Carlo approximations}
\label{MC:replacing}

Let $\mathscr{N} \colon \mathfrak{P} \times \mathbb{R}^{ \mathfrak{M} \mathbf{d}} \to \mathbb{R}$ satisfy for all
	$p \in \mathfrak{P}$, 
	$\theta = ( \theta_1, \dots, \theta_{ \mathfrak{M} \mathbf{d} } ) \in \mathbb{R}^{\mathfrak{M} \mathbf{d}}$
that
\begin{equation}
\label{replacing:eq1}
\begin{split} 
	\mathscr{N}(p, \theta)
=
    \frac{ 1 }{ \mathfrak{M} }
      \left[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \phi\big( 
          p, 
          ( 
            \theta_{ ( \mathfrak{m} - 1 ) \mathbf{d} + k }
          )_{k \in \{1, 2, \ldots, \mathbf{d} \}}
        \big)
      \right].
\end{split}
\end{equation}
Note that \eqref{mc:eq1} and \eqref{replacing:eq1} suggest that 
for all $ p \in \mathfrak{P} $ it holds that
\begin{equation}
  u( p )
  \approx 
  \frac{ 1 }{ \mathfrak{M} }
  \left[
    \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
    \phi( 
      p, W^{ 0, \mathfrak{m} }
    )
  \right]
 =
 \mathscr{N}\!\big(p,  (W^{ 0, 1 }, W^{ 0, 2 }, \ldots, W^{ 0,  \mathfrak{M}})\big)
  .
\end{equation}

\subsection{Random loss functions for fixed random variables in Monte Carlo approximations} 
\label{MC:Loss}

Let $ {\bf M} \in \mathbb{N} $, 
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
for every 
	$m \in \mathbb{N}$
let 
$ 
  F_m 
  \colon \mathbb{R}^{ \mathfrak{M} \mathbf{d} } 
  \times \Omega \to \mathbb{R} 
$
satisfy 
for all 
$ 
  \theta = 
  ( \theta_1, \dots, \theta_{ \mathfrak{M} \mathbf{d} } ) 
  \in \mathbb{R}^{ \mathfrak{M} \mathbf{d} } 
$ 
that
\begin{equation}
\begin{split}
  F_m( \theta ) 
  &=
  \frac{ 1 }{ {\bf M} }
  \left[
    \sum_{ {\bf m} = 1 }^{ {\bf M} }
    \left|
      \phi( 
        P_{ 
          m, {\bf m} 
        }, 
        W^{ m, {\bf m} } 
      )
      -
      \mathscr{N}(
        P_{ 
          m, {\bf m} 
        }, 
        \theta
      )
    \right|^2
  \right] 
\\&= 
  \frac{ 1 }{ {\bf M} }
  \left[
    \sum_{ {\bf m} = 1 }^{ {\bf M} }
    \left|
      \phi( 
        P_{ 
          m, {\bf m} 
        }, 
        W^{ m, {\bf m} } 
      )
      -
      \frac{ 1 }{ \mathfrak{M} }
      \left[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \phi\big( 
          P_{ 
            m, {\bf m} 
          }, 
          ( 
            \theta_{ ( \mathfrak{m} - 1 ) \mathbf{d} + k }
          )_{k \in \{1, 2, \ldots, \mathbf{d} \}}
        \big)
      \right]
    \right|^2
  \right] ,
\end{split}
\end{equation}
and for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{ \mathfrak{M} \mathbf{d} } \times \Omega 
  \to \mathbb{R}^{ \mathfrak{M} \mathbf{d} } 
$
satisfy 
for all 
	$ \theta \in \mathbb{R}^{ \mathfrak{M} \mathbf{d} } $, 
	$ \omega \in \Omega $
that
\begin{equation}
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega ).
\end{equation}


\subsection{Learning the random variables with stochastic gradient descent}
\label{MC:SGD}

Let $ (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) $ and 
let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \mathfrak{M} \mathbf{d} } 
$ 
satisfy 
for all $ m \in \mathbb{N} $ that 
$ 
  \Theta_0 = 
  ( W^{ 0, 1 }, W^{ 0, 2 }, \dots, W^{ 0, \mathfrak{M} } )
$
and
\begin{equation}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
\end{equation}
\new{
For every sufficiently large 
$m \in \mathbb{N}$
we propose to employ the random function 
$
	\mathfrak{P} \times \Omega \ni (p, \omega) \mapsto \mathscr{N}(p, \Theta_m(\omega)) \in \mathbb{R}
$
as an approximation for the target function
$
	\mathfrak{P} \ni p \mapsto u(p) \in \mathbb{R}
$ in \eqref{MC:eq1}. 
}



\subsection{Description of the proposed approximation algorithm}
\label{MC:description}

\begin{algo}
Let $ \mathfrak{p},\mathbf{d}, \mathfrak{M}, {\bf M}\in \mathbb{N} $, 
$ (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) $,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let $ \phi \in C( \mathfrak{P} \times \mathbb{R}^{ \mathbf{d} } , \mathbb{R} ) $, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
let 
$ 
  W^{ m, \mathfrak{m} } \colon \Omega \to \mathbb{R}^{ \mathbf{d} } 
$, 
$ m, \mathfrak{m} \in \mathbb{N}_0 $, 
be i.i.d.\ random variables, 
assume that 
$
  (
    P_{ 
      m, {\bf m} 
    }
  )_{
    (m, {\bf m}) \in \mathbb{N}^2
  }
$
and 
$
  (
    W^{ m, \mathfrak{m} }
  )_{
    (m, \mathfrak{m}) \in \mathbb{N}^2
  }
$
are independent, 
assume for all $ p \in \mathfrak{P} $ that
$
  \mathbb{R}^{ \mathbf{d} } \ni w \mapsto \phi( p, w ) \in \mathbb{R} 
$
is continuously differentiable, 
for every 
	$m \in \mathbb{N}$
let 
$ 
  F_m 
  \colon \mathbb{R}^{ \mathfrak{M} \mathbf{d} } 
  \times \Omega \to \mathbb{R} 
$
satisfy 
for all 
$ 
  \theta = 
  ( \theta_1, \dots, \theta_{ \mathfrak{M} \mathbf{d} } ) 
  \in \mathbb{R}^{ \mathfrak{M} \mathbf{d} } 
$ 
that
\begin{equation}
  F_m( \theta ) 
  =
  \frac{ 1 }{ {\bf M} }
  \left[
    \sum_{ {\bf m} = 1 }^{ {\bf M} }
    \left|
      \phi( 
        P_{ 
          m, {\bf m} 
        }, 
        W^{ m, {\bf m} } 
      )
      -
      \frac{ 1 }{ \mathfrak{M} }
      \left[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \phi\big( 
          P_{ 
            m, {\bf m} 
          }, 
          ( 
            \theta_{ ( \mathfrak{m} - 1 ) \mathbf{d} + k }
          )_{k \in \{1, 2, \ldots, \mathbf{d} \}}
        \big)
      \right]
    \right|^2
  \right] ,
\end{equation}
for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{ \mathfrak{M} \mathbf{d} } \times \Omega 
  \to \mathbb{R}^{ \mathfrak{M} \mathbf{d} } 
$
satisfy 
for all 
	$ \theta \in \mathbb{R}^{ \mathfrak{M} \mathbf{d} } $, $ \omega \in \Omega $ 
that
$
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega )
$, 
and
let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \mathfrak{M} \mathbf{d} } 
$ 
satisfy 
for all $ m \in \mathbb{N} $ that 
$ 
  \Theta_0 = 
  ( W^{ 0, 1 }, W^{ 0, 2 }, \dots, W^{ 0, \mathfrak{M} } )
$
and
\begin{equation}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
\end{equation}
\end{algo}






\section{LRV strategy in the case of Monte Carlo-Euler-Maruyama approximations}
\label{sect:MCEUler}

In the previous section we derived the LRV strategy for the approximation of parametric expectations involving finite-dimensional random variables.
In this section we consider parametric expectations involving standard Brownian motions as random variables (see \eqref{MCEUler:eq0} in \cref{MCEUler:target}).
To derive an LRV algorithm for these approximation problems we first discretize the Brownian motions in \cref{MCEUler:discretization} to obtain parametric expectations involving only finite-dimensional random variables. 
In \cref{MCEUler:proposal,MCEUler:replacing,MCEUler:Loss,MCEUler:SGD} the LRV strategy is then subsequently applied to these discretized parametric expectations as in \cref{sect:MC} above (cf.\ \cref{MC:proposal,MC:replacing,MC:Loss,MC:SGD}).
In \cref{MCEUler:example} we illustrate a special instance
of the parametric expectations in \eqref{MCEUler:eq0} and of associated discretizations in the situation where the parametric expectations in \eqref{MCEUler:eq0} involve solutions of SDEs.
Finally, in \cref{MCEUler:description} we summarize the algorithm derived in this section in one single framework.


\subsection{Parametric expectations involving vector valued stochastic processes}
\label{MCEUler:target}

Let $ \mathfrak{p},d \in \mathbb{N} $, $T \in (0, \infty)$,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let $ u \colon \mathfrak{P} \to \mathbb{R} $ be a function, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let $ \mathcal{W} \colon [0,T] \times \Omega \to \mathbb{R}^{ d } $ be a 
standard Brownian motion with continuous sample paths, 
let $ \phi \colon \mathfrak{P} \times C( [0,T], \mathbb{R}^{ d } ) \to \mathbb{R} $ be measurable, 
and
assume for all $ p \in \mathfrak{P} $ that 
$
  \mathbb{E}\big[
    | \phi( p, \mathcal{W} ) |
  \big]
  < \infty
$
and
\begin{equation}
\label{MCEUler:eq0}
  u( p ) = 
  \mathbb{E}\!\left[ 
    \phi( 
      p, \mathcal{W}
    )
  \right].
\end{equation}
The goal of this section is to derive an  algorithm to approximately compute the function
$
	u \colon \mathfrak{P} \to \mathbb{R}
$
given through the parametric expectation in \eqref{MCEUler:eq0}.

\subsection{Temporal discretizations of the Brownian motion}
\label{MCEUler:discretization}

Let $N \in \mathbb{N}$,
let $ \mathfrak{W} \colon \Omega \to \mathbb{R}^{Nd}$ satisfy 
$
	\mathfrak{W}
=
	\sqrt{{N}/{T}}\,
	(
		\mathcal{W}_{T / N} - \mathcal{W}_0, 
		\mathcal{W}_{ 2 T / N } - \mathcal{W}_{ T / N }, 
		\dots, 
		\mathcal{W}_T -\mathcal{W}_{ (N-1) T / N } 
	)
$,
let
$
\Phi \colon
\mathfrak{P} \times \mathbb{R}^{N d} \to
\mathbb{R}
$
satisfy for all $ p \in \mathfrak{P} $ that
\(
  \mathbb{R}^{N d} \ni w \mapsto \Phi(p, w) \in \mathbb{R}
\)
is continuously differentiable,
and assume for all $p \in \mathfrak{P}$ that
$\mathbb{E}[|\Phi(p, \mathfrak{W})|]<\infty$.
We think of
$\Phi \colon
\mathfrak{P} \times \mathbb{R}^{N d} \to
\mathbb{R}$
as a suitable approximation of
$\phi \colon  \mathfrak{P} \times C( [0,T], \mathbb{R}^{ d } ) \to \mathbb{R}$
in the sense that for all $ p \in \mathfrak{P} $ it holds that
\begin{equation}
  \label{eq:u_SP_approx} 
    \phi( 
      p, \mathcal{W}
    )
  \approx 
    \Phi\big( 
    p, \sqrt{{N}/{T}}\,
    (\mathcal{W}_{T / N} - \mathcal{W}_0, \mathcal{W}_{ 2 T / N } - \mathcal{W}_{ T / N }, \dots, \mathcal{W}_T -\mathcal{W}_{ (N-1) T / N } )
    \big)
  =
 \Phi(p, \mathfrak{W}).
\end{equation}


\subsection{Euler-Maruyama approximations for parametric stochastic differential equations}
\label{MCEUler:example}



In the case 
where $ \mathfrak{p} = 1 + d $,
where
$
\mathfrak{P} = [0,T] \times [-10,10]^d
$,
where
$\mu \colon \mathbb{R}^d \to \mathbb{R}^d$ is globally Lipschitz continuous,
where
$
X = (X^{\xi, w}_t)_{(t, \xi, w) \in [0,T] \times \mathbb{R}^{d} \times C([0,T],\mathbb{R}^d)} \colon [0,T] \times \mathbb{R}^{d} \times C([0,T],\mathbb{R}^d) \to \mathbb{R}^{d}
$
satisfies for all
	$t \in [0,T]$,
	$\xi \in \mathbb{R}^{d}$,
	$w=(w_s)_{s\in [0,T]} \in C([0,T], \mathbb{R}^d)$
that 
\begin{equation}
\label{MCeuler:eq1}
  X_t^{\xi, w} = \xi + \int_0^t \mu(X_s^{\xi, w}) \d s + w_t,
\end{equation}
where $g\colon [0,T] \times \mathbb{R}^{d} \to \mathbb{R}$
satisfies for all $t \in [0,T]$ that 
$\mathbb{R}^{d} \ni x \mapsto g(t, x) \in \mathbb{R}$
is continuously differentiable,
where it holds for all
$p=(t,\xi)\in \mathfrak{P}$,
$
w \in C([0,T], \mathbb{R}^d)
$
that
$\phi(p,w) = g(t, X^{\xi,w}_t)$,
where
$
\mathcal{X}^{\xi, \theta} = (\mathcal{X}^{\xi, \theta}_t)_{t \in [0,T]}\colon [0,T] \to \mathbb{R}^d
$,
	$\theta \in \mathbb{R}^{N d}$,
	$\xi \in \mathbb{R}^{d}$,
satisfy for all 
	$\theta \in \mathbb{R}^{N d}$,
	$\xi \in \mathbb{R}^{d}$,
	$n\in \{1, 2, \dots N\}$,
	$t\in [\tfrac{(n-1)T}{N},\tfrac{nT}{N}]$
that $\mathcal{X}_0^{\xi,\theta} = \xi$
and
\begin{equation}
\label{MCeuler:eq2}
  \mathcal{X}_{t}^{\xi, \theta} =  \mathcal{X}_{(n-1)T/N}^{\xi, \theta} + \big(\tfrac{tN}{T} - n + 1\big)\bigl(\tfrac{T}{N} \mu(\mathcal{X}_{(n-1)T/N}^{\xi, \theta}) + \tfrac{\sqrt{T}}{\sqrt{N}} (\theta_{(n-1)d + k})_{k \in \{1, 2, \ldots, d\}}\bigr),
\end{equation}
and where it holds for all
$
p \in \mathfrak{P}
$,
$
\theta \in \mathbb{R}^{Nd}
$
that
$\Phi( p, \theta ) = g(t, \mathcal{X}^{\xi, \theta}_t)$
observe that 
\begin{enumerate}[label=(\roman *)]
\item 
it holds 
for all 
	$ p = (t,\xi)\in \mathfrak{P} $
that
\begin{equation}
\begin{split} 
	u(p) = u(t, \xi) = \mathbb{E}[g(t,X^{\xi,\mathcal{W}}_t) ]
\end{split}
\end{equation}
is the expectation of the test function $\mathbb{R}^d \ni x\mapsto g(t,x) \in \mathbb{R}$ evaluated at time $t$ of the solution process $(X^{\xi,\mathcal{W}}_s)_{s\in [0,T]}$ of the additive noise driven SDE in \eqref{MCeuler:eq1}
and
\item 
it holds for all
	$p = (t,\xi)\in \mathfrak{P} $
that
$\Phi(p, \mathfrak{W})$ is an approximation
\begin{equation}
\begin{split} 
	\Phi(p, \mathfrak{W}) = g(t, \mathcal{X}^{\xi, \mathfrak{W}}_t)
\approx 
	g(t, X^{\xi,\mathcal{W}}_t) = \phi(p,\mathcal{W})
\end{split}
\end{equation}
of $\phi(p,\mathcal{W})$ based on linearly interpolated Euler-Maruyama approximations $(\mathcal{X}_{s}^{\xi, \mathfrak{W}})_{s \in [0,T]}$ with $N$ timesteps of the solution $(X_{s}^{\xi, \mathcal{W}})_{s \in [0,T]}$ of the SDE in \eqref{MCeuler:eq1}.


\end{enumerate}








\subsection{Monte Carlo approximations}
\label{MCEUler:proposal}

Let $ \mathfrak{M} \in \mathbb{N} $
and 
let 
$ 
  W^{ m, \mathfrak{m} } \colon \Omega \to \mathbb{R}^{ N d } 
$, 
$ m, \mathfrak{m} \in \mathbb{N}_0 $, 
be i.i.d.\ standard normal random vectors.
Observe that \eqref{MCEUler:eq0} and \eqref{eq:u_SP_approx} suggest that 
for all $ p \in \mathfrak{P} $ it holds that
\begin{equation}
\label{mc:eq2}
    u( p ) = 
  \mathbb{E}\!\left[ 
    \phi( 
      p, \mathcal{W}
    )
  \right]
  \approx
  \mathbb{E}[
  	\Phi(p, \mathfrak{W})
  ]
  =
  \mathbb{E}[
    \Phi(p, W^{0, 1})
  ]
  \approx
  \frac{ 1 }{ \mathfrak{M} }
  \left[
    \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
    \Phi\big( 
    p, 
    W^{0,\mathfrak{m}}
    \big)
  \right]  
  .
\end{equation}


\subsection{Replacing the random variables in Monte Carlo approximations}
\label{MCEUler:replacing}

Let $\mathscr{N} \colon \mathfrak{P} \times\mathbb{R}^{ \mathfrak{M} N d }  \to \mathbb{R}$ satisfy for all
	$p \in \mathfrak{P}$, 
	$\theta = (\theta_1, \ldots, \theta_{\mathfrak{M} Nd}) \in \mathbb{R}^{ \mathfrak{M} N d } $
that
\begin{equation}
\label{replacing:eq2}
\begin{split} 
	\mathscr{N}(p, \theta)
=
   \frac{ 1 }{ \mathfrak{M} }
      \left[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \Phi\Big( 
          p, 
          (\theta_{( \mathfrak{m} - 1 ) N d + k})_{k \in\{1, 2, \dots, N d\}}          
        \Big)
      \right].
\end{split}
\end{equation}
Note that \eqref{mc:eq2} and \eqref{replacing:eq2} suggest that 
for all $ p \in \mathfrak{P} $ it holds that
\begin{equation}
  u( p )
  \approx
    \frac{ 1 }{ \mathfrak{M} }
    \left[
      \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
      \Phi\big( 
      p, 
      W^{0, \mathfrak{m}}
      \big)
    \right]  
 =
  \mathscr{N}\!\big(
    p,  
      (
        W^{ 0,  1  },
        W^{ 0,  2  },
        \ldots,
        W^{ 0,  \mathfrak{M}  }
       )
  \big)
  .
\end{equation}


\subsection{Random loss functions for fixed random variables in Monte Carlo approximations} 
\label{MCEUler:Loss}
Let $ {\bf M} \in \mathbb{N} $,  
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
for every 
	$m \in \mathbb{N}$
let
$
F_m
  \colon \allowbreak\mathbb{R}^{ \mathfrak{M} N d } 
  \times\Omega\to\mathbb{R}
$
satisfy 
for all 
$
\theta = (\theta_1, \ldots, \theta_{\mathfrak{M} Nd}) \in \mathbb{R}^{\mathfrak{M}Nd}
$
that
\begin{equation}
\begin{split}
  F_m( \theta ) 
  &=
  \frac{ 1 }{ {\bf M} }
  \left[
    \sum_{ {\bf m} = 1 }^{ {\bf M} }
    \left|
      \Phi( 
        P_{ 
          m, {\bf m} 
        }, 
        W^{ m, {\bf m} } 
      )
      -
      \mathscr{N}(
        P_{ 
          m, {\bf m} 
        }, 
        \theta
      )
    \right|^2
  \right] 
\\&= 
  \frac{ 1 }{ {\bf M} }
  \left[
    \sum_{ {\bf m} = 1 }^{ {\bf M} }
    \left|
      \Phi( 
        P_{ 
          m, {\bf m} 
        }, 
        W^{ m, {\bf m} } 
      )
      -
      \frac{ 1 }{ \mathfrak{M} }
      \left[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \Phi\big( 
          P_{ 
            m, {\bf m} 
          }, 
          ( 
            \theta_{ ( \mathfrak{m} - 1 ) Nd + k }
          )_{k \in \{1, 2, \ldots, Nd \}}
        \big)
      \right]
    \right|^2
  \right] ,
\end{split}
\end{equation}
and for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{\mathfrak{M} N d } \times \Omega 
  \to \mathbb{R}^{\mathfrak{M} N d } 
$
satisfy 
for all 
	$ \theta \in \mathbb{R}^{\mathfrak{M}Nd} $, $ \omega \in \Omega $ 
that
\begin{equation}
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega ).
\end{equation}


\subsection{Learning the random variables with stochastic gradient descent}
\label{MCEUler:SGD}

Let
$
  (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) 
$
and let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{\mathfrak{M}Nd} 
$ 
satisfy for all
$ m \in \mathbb{N} $ that 
$
\Theta_0
=
(W^{0,1}, W^{0,2}, \ldots, W^{0,\mathfrak{M}})
$
and
\begin{equation}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
\end{equation}
\new{
For every sufficiently large 
$m \in \mathbb{N}$
we propose to employ the random function 
$
	\mathfrak{P} \times \Omega \ni (p, \omega) \mapsto \mathscr{N}(p, \Theta_m(\omega)) \in \mathbb{R}
$
as an approximation for the target function
$
	\mathfrak{P} \ni p \mapsto u(p) \in \mathbb{R}
$ in \eqref{MCEUler:eq0}.}




\subsection{Description of the proposed approximation algorithm}
\label{MCEUler:description}

\begin{algo}
Let $\mathfrak{p}, d, N, \mathfrak{M}, {\bf M} \in \mathbb{N} $, 
$ (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) $,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
  $, $ m, {\bf m} \in \mathbb{N} $,
be i.i.d.\ random variables, 
let 
$ 
  W^{ m, \mathfrak{m} } \colon [0, T] \times \Omega \to \mathbb{R}^{ d } 
$, 
$ m, \mathfrak{m} \in \mathbb{N}_0 $, 
be i.i.d.\ standard Brownian motions,
assume that 
$
  (
    P_{ 
      m, {\bf m} 
    }
  )_{
    (m, {\bf m}) \in \mathbb{N}^2
  }
$
and 
$
  (
    W^{ m, \mathfrak{m} }
  )_{
    (m, \mathfrak{m}) \in \mathbb{N}^2
  }
$
are independent, 
let
$
\Phi \colon
\mathfrak{P} \times \mathbb{R}^{N d} \to
\mathbb{R}
$ satisfy for all $ p \in \mathfrak{P} $ that
\(
  \mathbb{R}^{N d} \ni w \mapsto \Phi(p, w) \in \mathbb{R}
\)
is continuously differentiable,
for every 
	$m \in \mathbb{N}$
let 
$ 
  F_m 
  \colon \mathbb{R}^{\mathfrak{M}Nd} 
  \times \Omega \to \mathbb{R} 
$
satisfy 
for all 
$
\theta = (\theta_1, \ldots, \theta_{\mathfrak{M}Nd})
  \in \mathbb{R}^{\mathfrak{M}Nd} 
$
that
\begin{multline}
  F_m( \theta ) 
  =
  \frac{ 1 }{ {\bf M} }
  \left[
    \sum_{ {\bf m} = 1 }^{ {\bf M} }
    \left|
    \Phi\Big( 
        P_{ 
          m, {\bf m} 
        }, 
        \sqrt{{N}/{T}}\,
        (W^{ m,  {\bf m}  }_{T / N} - W^{ m,  {\bf m}  }_{ 0 }, W^{ m,  {\bf m}  }_{ 2 T / N } - W^{ m,  {\bf m}  }_{ T / N }, \dots, W^{ m,  {\bf m}  }_T -W^{ m,  {\bf m}  }_{ (N-1) T / N } )    
       \Big)
       \vphantom{
       - \frac{ 1 }{ \mathfrak{M} }
      \left[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \Phi\Big( 
          P_{ 
            m, {\bf m} 
          }, 
          (\theta_{( \mathfrak{m} - 1 ) d + k, n })_{(k, n) \in\{1, 2, \dots, d\} \times \{1, 2, \dots, N\}}          
        \Big)
      \right]
       }       \right.
       \vphantom{\left.
       - \frac{ 1 }{ \mathfrak{M} }
      \left[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \Phi\Big( 
          P_{ 
            m, {\bf m} 
          }, 
          (\theta_{( \mathfrak{m} - 1 ) d + k, n })_{(k, n) \in\{1, 2, \dots, d\} \times \{1, 2, \dots, N\}}
        \Big)
      \right]
      \right|^2}
       \right.\\
      \left.\left.
       - \frac{ 1 }{ \mathfrak{M} }
      \left[ 
        \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
        \Phi\Big( 
          P_{ 
            m, {\bf m} 
          }, 
          ( 
            \theta_{ ( \mathfrak{m} - 1 ) Nd + k }
          )_{k \in \{1, 2, \ldots, Nd \}}         
        \Big)
      \right]
    \right|^2
  \right] ,
\end{multline}
for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{\mathfrak{M} Nd } \times \Omega 
  \to \mathbb{R}^{\mathfrak{M} Nd } 
$
satisfy 
for all 
	$ \theta \in \mathbb{R}^{\mathfrak{M}Nd } $, 
	$ \omega \in \Omega $ 
that
$
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega )
  $,
and
let 
$ 
  \Theta = (\Theta^{(1)}, \ldots, \Theta^{(\mathfrak{M}Nd)}) \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{\mathfrak{M}Nd } 
$ 
satisfy 
for all
	$ m \in \mathbb{N} $,
	$k \in \{1, 2, \ldots, \mathfrak{M}\}$,	
	$n \in \{1, 2, \ldots, N\}$
that 
$
	(
		\Theta_0^{([(k-1)N + n-1]d + 1)},
		\Theta_0^{([(k-1)N + n-1]d + 2)}, \allowbreak
		\ldots,
		\Theta_0^{([(k-1)N + n-1]d + d)}
	)
=
  \sqrt{{N}/{T}}
   (W^{0, k}_{nT/ N} - W^{0, k}_{(n-1)T/ N})
$
and
$
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
$
\end{algo}


\section{LRV strategy in the case of multilevel Monte Carlo approximations}
\label{sect:MLMC}

In this section we consider the problem of approximating a parametric expectation involving a general measure space-valued random variable for which a sequence of finite-dimensional approximations is available (see \cref{MLMC:target,MLMC:discretization}).
We illustrate such a situation in \cref{MLMC:example} in which the original random variable is a Brownian motion driving an SDE and the finite-dimensional approximations consist of Euler discretizations of the SDE with decreasing step sizes. 
The employed proposal algorithm for the LRV strategy in this section is a form of MLMC method 
(cf.\ Heinrich \cite{Heinrich98} and Giles \cite{Giles2008} and cf., e.g., \cite{Heinrich2001,Creutzig09,hs99,Giles15})
as described in \cref{MLMC:proposal}.
Based on this proposal algorithm, the LRV strategy is applied to the considered approximation problem in \cref{MLMC:replacing,MLMC:Loss,MLMC:SGD} and, thereafter, summarized in one single framework in \cref{MLMC:description}.

\subsection{Parametric expectations involving measure space valued random variables}
\label{MLMC:target}

Let $ \mathfrak{p} \in \mathbb{N} $,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let $(V, \mathcal{V})$ be a measurable space,
let $ \phi \colon \mathfrak{P} \times V \to \mathbb{R} $ be continuous, 
let $ u \colon \mathfrak{P} \to \mathbb{R} $ be a function, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let $ \mathcal{W} \colon \Omega \to V $ be a random variable,
and
assume for all $ p \in \mathfrak{P} $ that 
$
  \mathbb{E}\big[
    | \phi( p, \mathcal{W} ) |
  \big]
  < \infty
$
and
\begin{equation}
\label{MLMC:eq1}
  u( p ) = 
  \mathbb{E}\!\left[ 
    \phi( 
      p, \mathcal{W}
    )
  \right].
\end{equation}
The goal of this section is to derive an  algorithm to approximately compute the function
$
	u \colon \mathfrak{P} \to \mathbb{R}
$
given through the parametric expectation in \eqref{MLMC:eq1}.


\subsection{Finite-dimensional approximations}
\label{MLMC:discretization}

Let $\mathbf{d} = (\mathbf{d}_l)_{l \in \mathbb{N}_0} \colon \mathbb{N}_0 \to \mathbb{N}$ be non-decreasing,
let $\mathcal{P}_l \colon \mathbb{R}^{\mathbf{d}_l} \to \mathbb{R}^{\mathbf{d}_{l-1}}$, $l \in \mathbb{N}$, be continuously differentiable,
let
$\mathfrak{W}_l \colon \Omega \to \mathbb{R}^{\mathbf{d}_l}$, $l \in \mathbb{N}_0$, be random variables which satisfy for all $l \in \mathbb{N}$ that 
$\mathcal{P}_{l}(\mathfrak{W}_{l}) = \mathfrak{W}_{l-1}$,
let 
$
\Phi_l \colon
\mathfrak{P} \times \mathbb{R}^{\mathbf{d}_l} \to
\mathbb{R}
$,
$l \in \mathbb{N}_0$,
satisfy for all 
	$l \in \mathbb{N}_0$,
	$ p \in \mathfrak{P} $
that
\(
  \mathbb{R}^{\mathbf{d}_l} \ni w \mapsto \Phi_l(p, w) \in \mathbb{R}
\)
is continuously differentiable,
and assume for all $ p \in \mathfrak{P} $, $l \in \mathbb{N}_0$ that
$\mathbb{E}[|\Phi_l(p, \mathfrak{W}_l)|]<\infty$.
We think of
$\mathfrak{W}_l$, $l \in \mathbb{N}_0$, 
and 
$\Phi_l$,
$l \in \mathbb{N}_0$,
as suitable approximations of $\mathcal{W}$ and 
$\phi$
in the sense that for all $ p \in \mathfrak{P} $ and all sufficiently large $l \in \mathbb{N}_0$ it holds that
\begin{equation}
  \label{MLMC:eq2} 
    \phi( 
      p, \mathcal{W}
    )
  \approx 
    \Phi_l( 
    p,
    \mathfrak{W}_l
    ).
\end{equation}



\subsection{Multilevel Monte Carlo approximations}
\label{MLMC:proposal}


Let $L \in \mathbb{N}_0$, 
$ \mathfrak{M}_0, \mathfrak{M}_1, \ldots, \mathfrak{M}_L \in \mathbb{N} $
satisfy $\mathfrak{M}_0 \geq \mathfrak{M}_1 \geq \ldots \geq \mathfrak{M}_L$, 
let 
$ 
  W_l^{ m, \mathfrak{m} } \colon \Omega \to \mathbb{R}^{ \mathbf{d}_l } 
$, 
$ l, m, \mathfrak{m} \in \mathbb{N}_0 $, 
be independent random variables,
and
assume 
for all
$ l,m, \mathfrak{m} \in \mathbb{N}_0 $, $A \in \mathcal{B}( \mathbb{R}^{ \mathbf{d}_l } ) $
that
$
  \P( W_l^{ m, \mathfrak{m}} \in A ) = \P( \mathfrak{W}_l \in A ) 
$.
Observe that \eqref{MLMC:eq1} and \eqref{MLMC:eq2} suggest that 
for all $ p \in \mathfrak{P} $ it holds that
\begin{equation}
\label{MLMC:eq3}
\begin{split} 
    u( p ) = 
  \mathbb{E}\!\left[ 
    \phi( 
      p, \mathcal{W}
    )
  \right]
  &\approx
  \mathbb{E}[
  	\Phi_L(p, \mathfrak{W}_L)
  ]
  =
  \mathbb{E}[
  	\Phi_0(p, \mathfrak{W}_0)
  ]  
  +
  \sum_{l = 1}^L
  	\mathbb{E}[
  	  	\Phi_l(p, \mathfrak{W}_l)
  	  	-
  	  	\Phi_{l-1}(p, \mathcal{P}_l(\mathfrak{W}_l))
  	  ] 
\\
  &\approx
  \frac{1}{\mathfrak{M}_0}
  \left[
	  \sum_{m = 1}^{\mathfrak{M}_0}
		 	\Phi_0(p, W_0^{0,m})
  \right]
  +
  \sum_{l = 1}^L
    \frac{1}{\mathfrak{M}_l}
    \left[
      \sum_{m = 1}^{\mathfrak{M}_l}
      \Phi_l(p, W_l^{0,m})
      -
      \Phi_{l-1}(p, \mathcal{P}_l(W_l^{0,m}))
    \right].
\end{split}
\end{equation}




\subsection{Multilevel Monte Carlo approximations for parametric stochastic differential equations}
\label{MLMC:example}

In the case 
where $d \in \mathbb{N}$,
where 
	$ \mathfrak{p} = 1 + d $, 
	$
	\mathfrak{P} = [0,T] \times [-1000,1000]^d
	$, and
	$(V, \mathcal{V}) = \linebreak (C([0,T],\mathbb{R}^d), \mathcal{B}(C([0,T],\mathbb{R}^d)))$,
where
$\mu \colon \mathbb{R}^d \to \mathbb{R}^d$ is globally Lipschitz continuous,
where
$
X = \linebreak (X^{\xi, w}_t)_{(t, \xi, w) \in [0,T] \times \mathbb{R}^d \times C([0,T],\mathbb{R}^d)} \colon [0,T] \times \mathbb{R}^d \times C([0,T],\mathbb{R}^d) \to \mathbb{R}^{d}
$
satisfies for all
$t \in [0,T]$,
$\xi \in \mathbb{R}^d$,
$
w=(w_s)_{s\in [0,T]} \in C([0,T], \mathbb{R}^d)
$
that 
\begin{equation}
  X_t^{\xi, w} = \xi + \int_0^t \mu(X_s^{\xi, w}) \d s + w_t,
\end{equation}
where $g\colon [0,T] \times \mathbb{R}^{d} \to \mathbb{R}$
satisfies for all $t \in [0,T]$ that 
$\mathbb{R}^{d} \ni x \mapsto g(t, x) \in \mathbb{R}$
is continuously differentiable,
where it holds for all
	$p=(t,\xi)\in \mathfrak{P}$,
	$w \in C([0,T], \mathbb{R}^d)$
that
	$\phi(p,w) = g(t, X^{\xi,w}_t)$,
where 
$\mathcal{W} \colon \Omega \to C([0,T],\mathbb{R}^d) $ is a standard Brownian motion, 
where for all $l \in \mathbb{N}_0$ it holds that $\mathbf{d}_l = 2^l d$,
where for all 
	$l \in \mathbb{N}_0$
it holds that
\begin{equation}
\label{MLMC:eq5}
\begin{split} 
	\mathfrak{W}_{l}
=
	\sqrt{{2^l}/{T}}\,
    (\mathcal{W}_{T / 2^l} - \mathcal{W}_0, \mathcal{W}_{ 2 T / 2^l } - \mathcal{W}_{ T / 2^l }, \dots, \mathcal{W}_T -\mathcal{W}_{ (2^l-1) T / 2^l } ),
\end{split}
\end{equation}
where
$
\mathcal{X}^{\xi, \theta, N} = (\mathcal{X}^{\xi, \theta, N}_t)_{t \in [0,T]}\colon [0,T] \to \mathbb{R}^d
$,
	$N \in \mathbb{N}$,
	$\theta \in \mathbb{R}^{N d}$,
	$\xi \in \mathbb{R}^{d}$,
satisfy for all 
	$N \in \mathbb{N}$,
	$\theta \in \mathbb{R}^{N d}$,
	$\xi \in \mathbb{R}^{d}$,
	$n\in \{1, 2, \dots N\}$,
	$t\in [\tfrac{(n-1)T}{N},\tfrac{nT}{N}]$
that $\mathcal{X}_0^{\xi,\theta,N} = \xi$
and
\begin{equation}
\label{MLMC:eq6}
  \mathcal{X}_{t}^{\xi, \theta, N} =  \mathcal{X}_{(n-1)T/N}^{\xi, \theta, N} + \big(\tfrac{tN}{T} - n + 1\big)\big(\tfrac{T}{N} \mu(\mathcal{X}_{(n-1)T/N}^{\xi, \theta, N}) + \tfrac{\sqrt{T}}{\sqrt{N}} (\theta_{(n-1)d + k})_{k \in \{1, 2, \ldots, d\}}\big),
\end{equation}
where for all
$l \in \mathbb{N}_0$,
$
p = (t, \xi) \in \mathfrak{P}
$,
$
\theta \in \mathbb{R}^{2^l d}
$
it holds that
$\Phi_l( p, \theta ) = g(t, \mathcal{X}^{\xi, \theta, 2^l}_t)$,
and where for all
	$l \in \{0, 1, \ldots, L\}$
it holds that
$\mathfrak{M}_l = 2^{L-l}$
observe that
\begin{enumerate}[label=(\roman *)]
\item
it holds for all 
	$l \in \mathbb{N}$, 
	$\theta_1, \theta_2\, \ldots, \theta_{2^l}, \vartheta_1, \vartheta_2\, \ldots, \vartheta_{2^{l-1}} \in \mathbb{R}^{d}$ 
with 
$
	\forall \, n \in \{1, 2, \ldots, 2^{l-1}\}
	\colon
	\vartheta_n = \frac{1}{\sqrt{2}}(\theta_{2n} + \theta_{2n-1})
$
that
\begin{equation}
\begin{split} 
	\mathcal{P}_l(\theta_1, \theta_2\, \ldots, \theta_{2^l}) = (\vartheta_1, \vartheta_2\, \ldots, \vartheta_{2^{l-1}}),
\end{split}
\end{equation}

\item 
it holds for all 
	$ p = (t,\xi)\in \mathfrak{P} $
that
\begin{equation}
\begin{split} 
	u(p) = u(t, \xi) = \mathbb{E}[g(t,X^{\xi,\mathcal{W}}_t) ]
\end{split}
\end{equation}
is the expectation of the test function $\mathbb{R}^d \ni x\mapsto g(t,x) \in \mathbb{R}$ evaluated at time $t$ of the solution process $(X^{\xi,\mathcal{W}}_s)_{s\in [0,T]}$ of the additive noise driven SDE in \eqref{MLMC:eq5},

\item 
it holds for all 
	$l \in \mathbb{N}_0$, 
	$p = (t,\xi)\in \mathfrak{P} $
that
$\Phi_l(p, \mathfrak{W}_l)$ is an approximation
\begin{equation}
\begin{split} 
	\Phi_l(p, \mathfrak{W}_l) = g(t, \mathcal{X}^{\xi, \mathfrak{W}_l, 2^l}_t)
\approx 
	g(t, X^{\xi,\mathcal{W}}_t) = \phi(p,\mathcal{W})
\end{split}
\end{equation}
of $\phi(p,\mathcal{W})$ based on linearly interpolated Euler-Maruyama approximations $(\mathcal{X}_{s}^{\xi, \mathfrak{W}_l, 2^l})_{s \in [0,T]}$ with $2^l$ timesteps of the solution $(X^{\xi,\mathcal{W}}_s)_{s\in [0,T]}$ of the SDE in \eqref{MLMC:eq5},
and

\item 
it holds for all 
	$p = (t,\xi)\in \mathfrak{P} $
that
\begin{equation}
\begin{split} 
  &\frac{1}{\mathfrak{M}_0}
  \left[
	  \sum_{m = 1}^{\mathfrak{M}_0}
		 	\Phi_0(p, W_0^{0,m})
  \right]
  +
  \sum_{l = 1}^L
    \frac{1}{\mathfrak{M}_l}
    \left[
      \sum_{m = 1}^{\mathfrak{M}_l}
      \Phi_l(p, W_l^{0,m})
      -
      \Phi_{l-1}(p, \mathcal{P}_l(W_l^{0,m}))
    \right]\\
&=
  \frac{1}{2^L}
  \left[
	  \sum_{m = 1}^{2^L}
		 	g(t, \mathcal{X}^{\xi, W_0^{0,m}, 1}_t)
  \right]
  +
  \sum_{l = 1}^L
    \frac{1}{2^{L-l}}
    \left[
      \sum_{m = 1}^{2^{L-l}}
      g(t, \mathcal{X}^{\xi, W_l^{0,m}, 2^l}_t)
      -
      g(t, \mathcal{X}^{\xi, \mathcal{P}_l(W_l^{0,m}), 2^{l-1}}_t)
    \right]   
\end{split}
\end{equation}
is a MLMC approximation of $u(t,\xi)$ as proposed in Giles \cite{Giles2008}.
\end{enumerate}


\subsection{Replacing the random variables in multilevel Monte Carlo approximations}
\label{MLMC:replacing}

Let $\mathscr{N} \colon \mathfrak{P} \times \mathbb{R}^{  \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L}  \to \mathbb{R}$ satisfy for all
	$p \in \mathfrak{P}$, 
	$\theta = (\theta_1, \ldots, \theta_{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L})
	\in \mathbb{R}^{ \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } $
that
\begin{multline}  
\label{MLMC:eq4}
	\mathscr{N}(p, \theta)
=
     \frac{1}{\mathfrak{M}_0}
     \left[
   	  \sum_{m = 1}^{\mathfrak{M}_0}
   		 	\Phi_0\big(p, (\theta_{(m-1)\mathbf{d}_0 + k})_{k \in \{1, 2, \ldots, \mathbf{d}_0\}}\big)
     \right]\\
      +
     \sum_{l = 1}^L
       \frac{1}{\mathfrak{M}_l}
       \Bigg[
         \sum_{m = 1}^{\mathfrak{M}_l}
         \Phi_l\big(p, (\theta_{
         	\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_{l-1} \mathbf{d}_{l-1} + (m-1)\mathbf{d}_l + k
         })_{k \in \{1, 2, \ldots, \mathbf{d}_l\}}\big)\\
          -
         \Phi_{l-1}\big(p, \mathcal{P}_l(
         	(\theta_{
         	         	\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_{l-1} \mathbf{d}_{l-1} + (m-1)\mathbf{d}_l + k
         	         })_{k \in \{1, 2, \ldots, \mathbf{d}_l\}}
         )\big)
       \Bigg].
\end{multline}
Note that \eqref{MLMC:eq3} and \eqref{MLMC:eq4} suggest that 
for all $ p \in \mathfrak{P} $ it holds that
\begin{equation}
\begin{split}
  u( p )
  &\approx
  \frac{1}{\mathfrak{M}_0}
  \left[
	  \sum_{m = 1}^{\mathfrak{M}_0}
		 	\Phi_0(p, W_0^{0,m})
  \right]
  +
  \sum_{l = 1}^L
    \frac{1}{\mathfrak{M}_l}
    \left[
      \sum_{m = 1}^{\mathfrak{M}_l}
      \Phi_l(p, W_l^{0,m})
      -
      \Phi_{l-1}(p, \mathcal{P}_l(W_l^{0,m}))
    \right] \\
  &=
    \mathscr{N}\big(p, (
      W_0^{0, 1}, W_0^{0, 2}, \ldots, W_0^{0, \mathfrak{M}_1},
      W_1^{0, 1}, W_1^{0, 2}, \ldots, W_1^{0, \mathfrak{M}_2},
      \ldots,
      W_L^{0, 1}, W_L^{0, 2}, \ldots, W_L^{0, \mathfrak{M}_L}
    )\big).
\end{split}
\end{equation}


\subsection{Random loss functions for fixed random variables in multilevel Monte Carlo approximations}
\label{MLMC:Loss}
Let ${\bf M} \in \mathbb{N}$, $\mathcal{L} \in \mathbb{N}_0$, 
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
assume that 
$
  (
    P_{ 
      m, {\bf m} 
    }
  )_{
    (m, {\bf m}) \in \mathbb{N}^2
  }
$
and
$
  (
  W_l^{ m, \mathfrak{m} }
  )_{
    (l, m, \mathfrak{m}) \in (\mathbb{N}_0)^3
  }
$
are independent,
for every 
	$m \in \mathbb{N}$
let
$
  F_m 
  \colon \allowbreak\mathbb{R}^{  \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } 
  \times\Omega\to\mathbb{R}
$
satisfy 
for all 
$\theta = (\theta_1, \ldots, \theta_{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L})
\in \mathbb{R}^{ \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } $
that
\begin{equation}
\begin{split} 
	F_m( \theta ) 
&=
	\frac{ 1 }{ {\bf M} }
	\Bigg[
		\sum_{ {\bf m} = 1 }^{ {\bf M} }
		\left|
			\Phi_\mathcal{L}\big( 
				P_{ 
					m, {\bf m} 
				},
				W_\mathcal{L}^{m,  {\bf m}}
			\big)
			-
			\mathscr{N}(
				P_{ 
					m, {\bf m} 
				},
				\theta
			)
		\right|^2
	\Bigg],
\end{split}
\end{equation}
and for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L} \times \Omega 
  \to \mathbb{R}^{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } 
$
satisfy 
for all 
  $ \theta \in \mathbb{R}^{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } $, 
  $ \omega \in \Omega $ 
that
\begin{equation}
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega ).
\end{equation}


\subsection{Learning the random variables with stochastic gradient descent}
\label{MLMC:SGD}
Let $ (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) $ and 
let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } 
$ 
satisfy 
for all $ m \in \mathbb{N} $ that 
$ 
  \Theta_0 = 
  (
    W_0^{0, 1}, W_0^{0, 2}, \ldots, W_0^{0, \mathfrak{M}_1},
    W_1^{0, 1}, W_1^{0, 2}, \ldots, W_1^{0, \mathfrak{M}_2},
    \ldots,
    W_L^{0, 1}, W_L^{0, 2}, \ldots, W_L^{0, \mathfrak{M}_L}
  )
$
and
\begin{equation}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
\end{equation}
\new{
For every sufficiently large 
$m \in \mathbb{N}$
we propose to employ the random function 
$
	\mathfrak{P} \times \Omega \ni (p, \omega) \mapsto \mathscr{N}(p, \Theta_m(\omega)) \in \mathbb{R}
$
as an approximation for the target function
$
	\mathfrak{P} \ni p \mapsto u(p) \in \mathbb{R}
$ in \eqref{MLMC:eq1}.}




\subsection{Description of the proposed approximation algorithm}
\label{MLMC:description}
\begin{algo}
Let $ \mathfrak{p}, {\bf M} \in \mathbb{N} $, $L, \mathcal{L} \in \mathbb{N}_0$, 
$ \mathfrak{M}_0, \mathfrak{M}_1, \ldots, \mathfrak{M}_L \in \mathbb{N} $,
$(\mathbf{d}_l)_{l \in \mathbb{N}_0} \subseteq \mathbb{N}$, $ (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) $,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let 
$
\Phi_l \colon
\mathfrak{P} \times \mathbb{R}^{\mathbf{d}_l} \to
\mathbb{R}
$,
$l \in \mathbb{N}_0$,
satisfy for all 
	$l \in \mathbb{N}_0$,
	$ p \in \mathfrak{P} $ 
that
\(
  \mathbb{R}^{\mathbf{d}_l} \ni w \mapsto \Phi_l(p, w) \in \mathbb{R}
\)
is continuously differentiable,
let $\mathcal{P}_l \colon \mathbb{R}^{\mathbf{d}_l} \to \mathbb{R}^{\mathbf{d}_{l-1}}$, $l \in \{1, 2, \ldots, L \}$, be continuously differentiable,
let $\mathscr{N} \colon \mathfrak{P} \times \mathbb{R}^{  \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L}  \to \mathbb{R}$ satisfy for all
	$p \in \mathfrak{P}$, 
	$\theta = (\theta_1, \ldots, \theta_{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L})
	\in \mathbb{R}^{ \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } $
that
\begin{multline}  
	\mathscr{N}(p, \theta)
=
     \frac{1}{\mathfrak{M}_0}
     \left[
   	  \sum_{m = 1}^{\mathfrak{M}_0}
   		 	\Phi_0\big(p, (\theta_{(m-1)\mathbf{d}_0 + k})_{k \in \{1, 2, \ldots, \mathbf{d}_0\}}\big)
     \right]\\
      +
     \sum_{l = 1}^L
       \frac{1}{\mathfrak{M}_l}
       \Bigg[
         \sum_{m = 1}^{\mathfrak{M}_l}
         \Phi_l\big(p, (\theta_{
         	\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_{l-1} \mathbf{d}_{l-1} + (m-1)\mathbf{d}_l + k
         })_{k \in \{1, 2, \ldots, \mathbf{d}_l\}}\big)\\
          -
         \Phi_{l-1}\big(p, \mathcal{P}_l(
         	(\theta_{
         	         	\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_{l-1} \mathbf{d}_{l-1} + (m-1)\mathbf{d}_l + k
         	         })_{k \in \{1, 2, \ldots, \mathbf{d}_l\}}
         )\big)
       \Bigg],
\end{multline}
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let 
$ 
  W_l^{ m, \mathfrak{m} } \colon \Omega \to \mathbb{R}^{ d } 
$, 
$l, m, \mathfrak{m} \in \mathbb{N}_0 $, 
be independent random variables which
satisfy 
for all
$l, m, \mathfrak{m} \in \mathbb{N}_0 $, $A \in \mathcal{B}( \mathbb{R}^{ \mathbf{d}_l } ) $
that
$
  \P( W_l^{ m, \mathfrak{m}} \in A ) = \P( W_l^{ 0, 0} \in A ) 
$,
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
assume that 
$
  (
    P_{ 
      m, {\bf m} 
    }
  )_{
    (m, {\bf m}) \in \mathbb{N}^2
  }
$
and 
$
  (
  W_l^{ m, \mathfrak{m} }
  )_{
    (l, m, \mathfrak{m}) \in (\mathbb{N}_0)^3
  }
$
are independent,
for every 
	$m \in \mathbb{N}$
let
$
F_m 
  \colon \allowbreak\mathbb{R}^{  \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } 
  \times\Omega\to\mathbb{R}
$
satisfy 
for all 
$
\theta 
  \in \mathbb{R}^{ \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L} 
$
that
\begin{equation}
\label{T_B_D}
\begin{split} 
	F_m( \theta ) 
&=
	\frac{ 1 }{ {\bf M} }
	\Bigg[
		\sum_{ {\bf m} = 1 }^{ {\bf M} }
		\left|
			\Phi_\mathcal{L}\big( 
				P_{ 
					m, {\bf m} 
				},
				W_\mathcal{L}^{m,  {\bf m}}
			\big)
			-	
			\mathscr{N}(
				P_{ 
					m, {\bf m} 
				},
				\theta)
		\right|^2
	\Bigg],
\end{split}
\end{equation}
for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L} \times \Omega 
  \to \mathbb{R}^{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } 
$
satisfy 
for all 
$ \theta \in \linebreak \mathbb{R}^{\mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } $, 
$ \omega \in \Omega $
that
$
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega )
$,
and
let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \mathfrak{M}_0 \mathbf{d}_0 + \mathfrak{M}_1 \mathbf{d}_1 + \ldots + \mathfrak{M}_L \mathbf{d}_L } 
$ 
satisfy 
for all $ m \in \mathbb{N} $ that 
$ 
  \Theta_0 = 
  (
    W_0^{0, 1}, W_0^{0, 2}, \ldots, W_0^{0, \mathfrak{M}_1},
    W_1^{0, 1}, W_1^{0, 2}, \allowbreak \ldots, W_1^{0, \mathfrak{M}_2}, 
    \ldots,
    W_L^{0, 1}, W_L^{0, 2}, \ldots, W_L^{0, \mathfrak{M}_L}
  )
$
and
\begin{equation}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
\end{equation}
\end{algo}


\section{LRV strategy in the case of multilevel Picard approximations}
\label{sect:MLP}


In this section we consider the case where the target function to be approximated with the LRV strategy is given as the solution of a suitable stochastic fixed point equation (SFPE); see \eqref{MLP:eq1} in \cref{MLP:target}.
The setup in \cref{MLP:target} includes, as important special cases, the solutions of several semilinear parabolic PDEs. 
Two examples of such semilinear parabolic PDEs
(heat PDEs with Lipschitz nonlinearities and
Black-Scholes PDEs with Lipschitz nonlinearities)
are presented  in \cref{MLP:example}.
As proposal algorithms for the LRV strategy we present in \cref{MLP:proposal} a slight generalization of the MLP algorithm for semilinear PDEs in Hutzenthaler et al.\ \cite{Hutzenthaler_2020} (cf.\ also E et al.\ \cite{E2019,EHutzenthalerJentzenKruse16published}).
In analogy to the previous sections,  in \cref{MLP:replacing,MLP:Loss,MLP:SGD} an algorithm for the considered approximation problem is derived based on the LRV strategy with this MLP algorithm as proposal algorithm. 
Finally, the problem and the algorithm are summarized in one single framework in \cref{MLP:description}.


\subsection{Parametric solutions of stochastic fixed point equations}
\label{MLP:target}

Let $ \mathfrak{p}, \mathbf{d} \in \mathbb{N} $,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let 
$ \phi \colon \mathfrak{P} \times \mathbb{R}^\mathbf{d} \times \mathbb{R} \to \mathbb{R} $,
$\mathcal{P} \colon \mathfrak{P} \times \mathbb{R}^\mathbf{d} \to \mathfrak{P}$, and
$ u \colon \mathfrak{P} \to \mathbb{R} $ be measurable, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let $ \mathcal{W} \colon \Omega \to \mathbb{R}^\mathbf{d} $ be a random variable,
and
assume for all $ p \in \mathfrak{P} $ that 
$
  \mathbb{E}\big[
    | \phi( p, \mathcal{W}, u( \mathcal{P}(p, \mathcal{W}))) |
  \big]
  < \infty
$
and
\begin{equation}
\label{MLP:eq1}
  u( p ) = 
  \mathbb{E} \big[ 
    \phi\big( p, \mathcal{W}, u\big( \mathcal{P}(p, \mathcal{W})\big)\big)
  \big].
\end{equation}
The goal of this section is to derive an  algorithm to approximately compute the function
$
	u \colon \mathfrak{P} \to \mathbb{R}
$
given through the SFPE in \eqref{MLP:eq1}.

\subsection{Approximations for parametric semilinear partial differential equations}
\label{MLP:example}
\subsubsection{Heat partial differential equations with Lipschitz nonlinearities}
In the case 
where $T \in (0,\infty)$, $d \in \mathbb{N}$,
where 
	$\mathfrak{p} = \mathbf{d} = 1 + d$ and
	$\mathfrak{P} = [0,T] \times\mathbb{R}^d$,
where $u \in C([0,T] \times\mathbb{R}^d,\mathbb{R})$ is at most polynomially growing,
where $f,g \in C(\mathbb{R}^d, \mathbb{R})$ are at most polynomially growing,
where for all 
	$p = (t, x) \in \mathfrak{P}$,
	$w = (w_1, w_2) \in \mathbb{R}^{d} \times [0,1]$,
	$v \in \mathbb{R}$
it holds that
$
	\phi(p, w, v)
=
	g(x + (T-t)^{\nicefrac{1}{2}}w_1) + (T-t)f(v)
$
and
$
	\mathcal{P}(p, w)
=
	(t + w_2(T-t), x + (w_2(T-t))^{\nicefrac{1}{2}}w_1)
$,
where $N \colon \Omega \to \mathbb{R}^{d}$ is a standard normal random vector,
where $R \colon \Omega \to [0,1]$ is a continuous uniformly distributed random variable on $[0,1]$ (is a $\mathcal{U}_{[0,1]}$-distributed random variable),
where $N$ and $R$ are independent, and
where $\mathcal{W} = (N, R)$
observe that 
\begin{enumerate}[label=(\roman *)]
\item 
it holds for all
	$p = (t, x) \in \mathfrak{P}$
that
\begin{equation}
\begin{split} 
	u(p) = u(t, x) 
&= 
	\Exp{
		g(x + (T-t)^{\nicefrac{1}{2}}N) + (T-t) f\big(u(t + R(T-t), x + (R(T-t))^{\nicefrac{1}{2}}N)\big)
	}\\
&=
	\Exp{
		g(x + (T-t)^{\nicefrac{1}{2}}N) + \int_t^T f\big(u(s, x + (s-t)^{\nicefrac{1}{2}}N)\big) \d s 
	}
\end{split}
\end{equation}
and
\item 
it holds that $u$ is a viscosity solution of 
\begin{equation}
	(\tfrac{\partial u}{\partial t})(t,x) 
	+
	\tfrac{1}{2} 
	(\Delta_x u)(t, x) 
	+
	f(u(t, x))
=
	0
\end{equation}
with $u(T, x) = g(x)$ for $(t, x) \in (0,T) \times \mathbb{R}^d$
(cf.\ Beck et al.\ \cite[Theorem 1.1]{beck21}).
\end{enumerate}

\subsubsection{Black-Scholes partial differential equations with Lipschitz nonlinearities}

In the case 
where $T \in (0,\infty)$, $d \in \mathbb{N}$,
where 
	$\mathfrak{p} = \mathbf{d} = 1 + 3d $ and
	$\mathfrak{P} = [0,T] \times\mathbb{R}^d\times\mathbb{R}^d\times\mathbb{R}^d$,
where for all 
	$\alpha, \beta \in \mathbb{R}^d$
it holds that 
$
	[0,T] \times\mathbb{R}^d \ni (t, x) 
	\mapsto
	u(t,x,\alpha,\beta) \in \mathbb{R}
$ 
is continuous and at most polynomially growing,
where 
	$X = (X^{(i)})_{i \in \{1, 2, \ldots, d\}} \colon [0,T] \times \mathbb{R}^d\times\mathbb{R}^d\times\mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d$
satisfies for all 
	$i \in \{1, 2, \ldots, d\}$
	$t \in [0,T]$,
	$
		x = (x_1,\ldots,x_d), 
		\alpha = (\alpha_1,\ldots,\alpha_d), 
		\beta = (\beta_1,\ldots,\beta_d),
		n = (n_1, \ldots, n_d)
	\in
		\mathbb{R}^{d}
	$
that
\begin{equation}
\begin{split} 
	X^{(i)}(t, x, \alpha, \beta, n)
=
	x_i \exp\!\left(
		\big(\alpha_i - \tfrac{|\beta_i|^2}{2}\big)t + t^{1/2}\beta_i n_i 
	\right),
\end{split}
\end{equation}
where $f,g \in C(\mathbb{R}^d, \mathbb{R})$ are at most polynomially growing,
where for all 
	$p = (t, x, \alpha, \beta) \in \mathfrak{P}$,
	$w = (w_1, w_2) \in \mathbb{R}^{d} \times [0,1]$,
	$v \in \mathbb{R}$
it holds that
$
	\phi(p, w, v)
=
	g(X(T-t, x, \alpha, \beta, w_1)) + (T-t)f(v)
$
and
$
	\mathcal{P}(p, w)
=
	(t + w_2(T-t),  X(w_2(T-t), x, \alpha, \beta, w_1))
$,
where $N \colon \Omega \to \mathbb{R}^{d}$ is a standard normal random vector,
where $R \colon \Omega \to [0,1]$ is an $\mathcal{U}_{[0,1]}$-distributed random variable,
where $N$ and $R$ are independent, and
where $\mathcal{W} = (N, R)$
observe that 
\begin{enumerate}[label=(\roman *)]
\item 
it holds for all
	$p = (t, x, \alpha, \beta) \in \mathfrak{P}$
that
\begin{equation}
\begin{split} 
	u(p) 
&= 
	u(t, x, \alpha, \beta) \\
&= 
	\Exp{
		g\big(
			X(T-t, x, \alpha, \beta, N)
		\big) 
		+ 
		(T-t)
		f\Big(
			u\big(t + R(T-t),  X(R(T-t), x, \alpha, \beta, N)\big)
		\Big)
	}\\
&=
	\Exp{
		g\big(
			X(T-t, x, \alpha, \beta, N)
		\big) 
		+ 
		\int_t^T
			f\big(
				u\big(s,  X(s-t, x, \alpha, \beta, N)\big)
			\big)
		\d s
	}
\end{split}
\end{equation}
and
\item 
it holds for all
$
	\alpha = (\alpha_1,\ldots,\alpha_d), 
	\beta = (\beta_1,\ldots,\beta_d)
	\in \mathbb{R}^d
$
that 
$
	[0,T] \times\mathbb{R}^d \ni (t, x) 
	\mapsto
	u(t,x,\alpha,\beta) \in \mathbb{R}
$ 
is a viscosity solution of 
\begin{multline}
	\big(\tfrac{\partial u}{\partial t}\big)(t,x,\alpha,\beta) 
	+
	\left[
  	\sum_{i = 1}^d 
	    \tfrac{ | \beta_i |^2 |x_i|^2}{2} 
  	    \big( \tfrac{\partial^2 u}{\partial (x_i)^2 } \big)(t,x,\alpha,\beta)
    \right] \\
    +
	\left[
    \sum_{i = 1}^d 
	    \alpha_i
	     x_i
	    \big( \tfrac{\partial u}{\partial x_i } \big)(t,x,\alpha,\beta)
	\right]
	+
  f( u(t, x,\alpha,\beta) )
=
	0
\end{multline}
with $u(T, x) = g(x)$ for $(t, x) = (t, x_1, x_2, \ldots, x_d) \in (0,T) \times \mathbb{R}^d$
(cf., e.g., Beck et al.\ \cite[Theorem 1.1]{beck21}).
\end{enumerate}



\subsection{Multilevel Picard approximations}
\label{MLP:proposal}
Let $\mathfrak{M} \in \mathbb{N}$, $\mathbf{I} = \cup_{n = 1}^\infty \mathbb{Z}^n$, 
let $W^{\mathbf{i}} \colon \Omega \to \mathbb{R}^\mathbf{d}$, $\mathbf{i} \in \mathbf{I}$, be i.i.d.\ random variables which satisfy for all 
	$ A \in \mathcal{B}( \mathbb{R}^{ \mathbf{d} } ) $ 
that 
$
  \P( W^{ 0 } \in A ) = \P( \mathcal{W} \in A ) 
$, and
let $V_n^\mathbf{i} \colon \mathfrak{P} \times \Omega \to \mathbb{R}$, $n \in \mathbb{Z}$, $\mathbf{i} \in \mathbf{I}$, satisfy for all 
	$n \in \mathbb{Z}$, 
	$\mathbf{i} \in \mathbf{I}$,
	$p \in \mathfrak{P}$
that
\begin{multline}
\label{MLP:eq2}
	V_n^\mathbf{i}(p)
=
	\sum_{k = 0}^{n-1}
		\frac{1}{\mathfrak{M}^{n-k}}
		\sum_{\mathfrak{m} = 1}^{\mathfrak{M}^{n-k}} \Big[
			\phi\big(
				p, 
				W^{(\mathbf{i},k,\mathfrak{m})},
				V_k^{(\mathbf{i},k,\mathfrak{m})}\big(
					\mathcal{P}(p, W^{(\mathbf{i},k,\mathfrak{m})})
				\big)
			\big)\\
			-
			\mathbbm{1}_\mathbb{N}(k)
			\phi\big(
				p, 
				W^{(\mathbf{i},k,\mathfrak{m})},
				V_{k-1}^{(\mathbf{i},k-1,-\mathfrak{m})}\big(
					\mathcal{P}(p, W^{(\mathbf{i},k,\mathfrak{m})})
				\big)
			\big)
		\Big].
\end{multline}
Under suitable assumptions we expect for sufficiently large $n \in \mathbb{N}$ that
\begin{equation}
\label{MLP:eq3}
\begin{split} 
	u(p)
\approx
	V_n^0(p).
\end{split}
\end{equation}


\subsection{Replacing the random variables in multilevel Picard approximations}
\label{MLP:replacing}

\newcommand{c}{c}

Let $(C_n)_{n \in \mathbb{Z}} \subseteq \mathbb{N}_0$ satisfy for all
	$n \in \mathbb{Z}$
that
\begin{equation}
\begin{split} 
	C_n = \sum_{k = 0}^{n-1} \mathfrak{M}^{n-k}(1 + C_k + C_{k-1}),
\end{split}
\end{equation}
for every 
	$n, \mathfrak{m} \in \mathbb{N}$, $l \in \{1, 2, \ldots, n\}$
let $c_{n,l, \mathfrak{m}} \in \mathbb{N}$ satisfy
\begin{equation}
\begin{split} 
	c_{n,l, \mathfrak{m}} 
= 
	\left[
		\sum_{k = 0}^{l-1} \mathfrak{M}^{n-k}(1 + C_k +C_{k-1})
	\right]
	+
	(\mathfrak{m}-1) (1 + C_l + C_{l-1})
	+ 
	1,
\end{split}
\end{equation}
let $\mathscr{N}_n \colon \mathfrak{P} \times \mathbb{R}^{C_n \mathbf{d}}$, $n \in \mathbb{N}_0$, satisfy for all
	$n \in \mathbb{N}$,
	$p \in \mathfrak{P}$,
	$\theta_1,\theta_2,\ldots, \theta_{C_n} \in  \mathbb{R}^{\mathbf{d}}$
that
	$\mathscr{N}_0(p) = 0$
and
\begin{multline}
\label{MLP:eq5}
	\mathscr{N}_n(p, (\theta_r)_{r \in \{1, 2, \ldots, C_n \}})
=
	\left[
		\frac{1}{\mathfrak{M}^{n}}
		\sum_{\mathfrak{m} = 1}^{\mathfrak{M}^{n}}
			\phi\big(
				p, 
				\theta_{\mathfrak{m}},
				0
			\big)
	\right]
\\ 
	+
	\sum_{k = 1}^{n-1}
		\frac{1}{\mathfrak{M}^{n-k}}
		\sum_{\mathfrak{m} = 1}^{\mathfrak{M}^{n-k}} \bigg[
			\phi\Big(
				p, 
				\theta_{c_{n,k,\mathfrak{m}}}, 
				\mathscr{N}_k \big(
					\mathcal{P}(p, \theta_{c_{n,k,\mathfrak{m}}})
					,
					(\theta_{c_{n,k,\mathfrak{m}} + r})_{r \in \{1, 2, \ldots, C_k\}}					
				\big)
			\Big)
\\ 
			-
			\phi\Big(
				p, 
				\theta_{c_{n,k,\mathfrak{m}}},
				\mathscr{N}_{k-1}\big(
					\mathcal{P}(p, \theta_{c_{n,k,\mathfrak{m}}})
					,
					(\theta_{c_{n,k,\mathfrak{m}} + C_k + r})_{r \in \{1, 2, \ldots, C_{k-1}\}}
				\big)
			\Big)
		\bigg],
\end{multline}
let $N \in \mathbb{N}$,
assume for all $p \in \mathfrak{P}$ that
$\mathbb{R}^{ \MLPcount_\MLPlevel \dimWRV } \ni \theta \mapsto \mathscr{N}( p, \theta ) \in \mathbb{R}$ is continuously differentiable,  
and
let 
$\mathfrak{W}_k \colon \Omega \to \mathbb{R}^{\mathbf{d}}$, 
	$k \in \{1, 2, \ldots, C_N\}$,
be i.i.d.\ random variables which satisfy
for all
	$B \in \mathcal{B}(\mathbb{R}^\mathbf{d})$
that
$
	\P(\mathfrak{W}_1 \in B)
=
	\P(\mathcal{W} \in B)
$.
Note that \eqref{MLP:eq2} and \eqref{MLP:eq5} assure that for all 
	$p \in \mathfrak{P}$,
	$B \in \mathcal{B}(\mathbb{R})$
it holds that
\begin{equation}
\label{MLP:eq6}
\begin{split} 
	\P\big(
		\mathscr{N}_N\big(
			p, (\mathfrak{W}_1, \mathfrak{W}_2, \ldots, \mathfrak{W}_{C_N})
		\big) 
		\in B
	\big)
=
	\P(V_N^0(p) \in B).
\end{split}
\end{equation}
Combining this and \eqref{MLP:eq3} suggests that  for all
	$p \in \mathfrak{P}$
it holds that
\begin{equation}
\label{MLP:eq7}
\begin{split} 
	\mathscr{N}_N(p, \mathfrak{W})
\approx
	u(p).
\end{split}
\end{equation}

\subsection{Random loss functions for fixed random variables in multilevel Picard approximations}
\label{MLP:Loss}
Let ${\bf M}, \mathcal{N} \in \mathbb{N}$, 
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
assume that 
$
  (
    P_{ 
      m, {\bf m} 
    }
  )_{
    (m, {\bf m}) \in \mathbb{N}^2
  }
$,
$
  (
    W^\mathbf{i}
  )_{\mathbf{i} \in \mathbf{I}}
$,
and 
$(\mathfrak{W}_k)_{k \in \{1, 2, \ldots, C_N\} }$
are independent,
for every 
	$m \in \mathbb{N}$
let
$
F_m 
  \colon \allowbreak\mathbb{R}^{  \MLPcount_\MLPlevel \dimWRV } 
  \times\Omega\to\mathbb{R}
$
satisfy 
for all 
$
\theta 
  \in \mathbb{R}^{\MLPcount_\MLPlevel \dimWRV} 
$
that
\begin{equation}
\label{T_B_D}
\begin{split} 
	F_m( \theta ) 
&=
	\frac{ 1 }{ {\bf M} }
	\Bigg[
		\sum_{ {\bf m} = 1 }^{ {\bf M} }
		\left|
			V_{\mathcal{N}}^{(m, {\bf m})}( 
				P_{ 
					m, {\bf m} 
				}
			)
			-
			\mathscr{N}_{N}(
				P_{ 
					m, {\bf m} 
				},
				\theta
			)
		\right|^2
	\Bigg],
\end{split}
\end{equation}
and for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{\MLPcount_\MLPlevel \dimWRV} \times \Omega 
  \to \mathbb{R}^{\MLPcount_\MLPlevel \dimWRV} 
$
satisfy 
for all 
	$ \theta \in \mathbb{R}^{\MLPcount_\MLPlevel \dimWRV} $, $ \omega \in \Omega $ 
that
\begin{equation}
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega ).
\end{equation}

\subsection{Learning the random variables with stochastic gradient descent}
\label{MLP:SGD}
Let $ (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) $ and 
let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \MLPcount_\MLPlevel \dimWRV} 
$ 
satisfy 
for all $ m \in \mathbb{N} $ that 
$ 
	\Theta_0 
= 
	(\mathfrak{W}_1, \mathfrak{W}_2, \ldots, \mathfrak{W}_{C_N})
$
and
\begin{equation}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
\end{equation}
\new{
For every sufficiently large 
$m \in \mathbb{N}$
we propose to employ the random function 
$
	\mathfrak{P} \times \Omega \ni (p, \omega) \mapsto \mathscr{N}_N(p, \Theta_m(\omega)) \in \mathbb{R}
$
as an approximation for the target function
$
	\mathfrak{P} \ni p \mapsto u(p) \in \mathbb{R}
$ in \eqref{MLP:eq1}.}


\subsection{Description of the proposed  approximation algorithm}
\label{MLP:description}

\begin{algo}
Let 
	$ \mathfrak{p}, \mathbf{d}, \mathfrak{M}, N, \mathcal{N}, {\bf M}  \in \mathbb{N} $, 
	$ (\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty) $,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let 
$ \phi \colon \mathfrak{P} \times \mathbb{R}^\mathbf{d} \times \mathbb{R} \to \mathbb{R} $ and
$\mathcal{P} \colon \mathfrak{P} \times \mathbb{R}^\mathbf{d} \to \mathfrak{P}$ 
be measurable, 
let $(C_n)_{n \in \mathbb{Z}} \subseteq \mathbb{N}_0$ satisfy for all
	$n \in \mathbb{Z}$
that
\begin{equation}
\begin{split} 
	C_n = \sum_{k = 0}^{n-1} \mathfrak{M}^{n-k}(1 + C_k + C_{k-1}),
\end{split}
\end{equation}
for every
	$n, \mathfrak{m} \in \mathbb{N}$, $l \in \{1, 2, \ldots, n\}$
let $ c_{n,l, \mathfrak{m}}\in \mathbb{N} $ satisfy
\begin{equation}
\begin{split} 
	c_{n,l, \mathfrak{m}} 
= 
	\left[
		\sum_{k = 0}^{l-1} \mathfrak{M}^{n-k}(1 + C_k +C_{k-1})
	\right]
	+
	(\mathfrak{m}-1) (1 + C_l + C_{l-1})
	+ 
	1,
\end{split}
\end{equation}
let $\mathscr{N}_n \colon \mathfrak{P} \times \mathbb{R}^{C_n \mathbf{d}}$, $n \in \mathbb{N}_0$, satisfy for all
	$n \in \mathbb{N}$,
	$p \in \mathfrak{P}$,
	$\theta_1,\theta_2,\ldots, \theta_{C_n} \in  \mathbb{R}^{\mathbf{d}}$
that
	$\mathscr{N}_0(p) = 0$
and
\begin{multline}
	\mathscr{N}_n(p, (\theta_1,\theta_2,\ldots, \theta_{C_n}))
=
	\left[
		\frac{1}{\mathfrak{M}^{n}}
		\sum_{\mathfrak{m} = 1}^{\mathfrak{M}^{n}}
			\phi\big(
				p, 
				\theta_{\mathfrak{m}},
				0
			\big)
	\right]
\\ 
	+
	\sum_{k = 1}^{n-1}
		\frac{1}{\mathfrak{M}^{n-k}}
		\sum_{\mathfrak{m} = 1}^{\mathfrak{M}^{n-k}} \bigg[
			\phi\Big(
				p, 
				\theta_{c_{n,k,\mathfrak{m}}}, 
				\mathscr{N}_k \big(
					\mathcal{P}(p, \theta_{c_{n,k,\mathfrak{m}}})
					,
					(\theta_{c_{n,k,\mathfrak{m}} + r})_{r \in \{1, 2, \ldots, C_k\}}
			\Big)
\\ 
			-
			\phi\Big(
				p, 
				\theta_{c_{n,k,\mathfrak{m}}},
				\mathscr{N}_{k-1}\big(
					\mathcal{P}(p, \theta_{c_{n,k,\mathfrak{m}}})
					,
					(\theta_{c_{n,k,\mathfrak{m}} + C_k + r})_{r \in \{1, 2, \ldots, C_{k-1}\}}
				\big)
			\Big)
		\bigg],
\end{multline}
assume for all $p \in \mathfrak{P}$ that
$\mathbb{R}^{ \MLPcount_\MLPlevel \dimWRV } \ni \theta \mapsto \mathscr{N}( p, \theta ) \in \mathbb{R}$ is continuously differentiable,  
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let $\mathcal{W} \colon \Omega \to \mathbb{R}^\mathbf{d}$ be a random variable,
let $\mathfrak{W}_k \colon \Omega \to \mathbb{R}^{\mathbf{d}}$,
	$k \in \{1, 2, \ldots, C_N\}$,
be i.i.d.\ random variables,
let $W^{ m, {\bf m}}_k \colon \Omega \to \mathbb{R}^\mathbf{d}$, $m, {\bf m} \in \mathbb{N}$, $k \in \{1, 2, \ldots, C_\mathcal{N}\}$, be i.i.d.\ random variables,
assume for all
	$B \in \mathcal{B}(\mathbb{R}^\mathbf{d})$
that
$
	\P(\mathfrak{W}_1 \in B)
=
	\P(\mathcal{W} \in B)
=
	\P(W^{1,1}_1 \in B)
$,
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
assume that 
$
  (
    P_{ 
      m, {\bf m} 
    }
  )_{
    (m, {\bf m}) \in \mathbb{N}^2
  }
$, 
$
  (
    W^{ m, {\bf m}}_k
  )_{(m, {\bf m},k) \in \mathbb{N}^2 \times \{1, 2, \ldots, C_\mathcal{N}\}}
$,
and 
$(\mathfrak{W}_k)_{k \in \{1, 2, \ldots, C_N\} }$
are independent,
for every 
	$m \in \mathbb{N}$
let
$
F_m 
  \colon \allowbreak\mathbb{R}^{  \MLPcount_\MLPlevel \dimWRV} 
  \times\Omega\to\mathbb{R}
$
satisfy 
for all 
$
\theta 
  \in \mathbb{R}^{\MLPcount_\MLPlevel \dimWRV} 
$
that
\begin{equation}
\begin{split} 
	F_m( \theta ) 
&=
	\frac{ 1 }{ {\bf M} }
	\Bigg[
		\sum_{ {\bf m} = 1 }^{ {\bf M} }
		\left|
			\mathscr{N}_{\mathcal{N}}( 
				P_{ 
					m, {\bf m} 
				},
				( W^{ m, {\bf m}}_1,  W^{ m, {\bf m}}_2, \ldots,  W^{ m, {\bf m}}_{C_\mathcal{N}})
			)
			-
			\mathscr{N}_N(
				P_{ 
					m, {\bf m} 
				},
				\theta
			)
		\right|^2
	\Bigg],
\end{split}
\end{equation}
for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{\MLPcount_\MLPlevel \dimWRV} \times \Omega 
  \to \mathbb{R}^{\MLPcount_\MLPlevel \dimWRV} 
$
satisfy 
for all 
	$ \theta \in \mathbb{R}^{\MLPcount_\MLPlevel \dimWRV} $, $ \omega \in \Omega $ 
that
$
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega )
$,
and
let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \MLPcount_\MLPlevel \dimWRV} 
$ 
satisfy 
for all $ m \in \mathbb{N} $ that 
$ 
  \Theta_0 = 
  (\mathfrak{W}_1, \mathfrak{W}_2, \ldots, \mathfrak{W}_{C_N})
$
and
$
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
$.
\end{algo}



\section{LRV strategy in the case of a general proposal algorithm}
\label{sect:general}

In this section we derive and formulate the LRV strategy in its most general form, which contains all the algorithms derived in the previous sections as special cases.
Roughly speaking, we want to approximate a target function 
	(cf.\ $u \colon \mathfrak{P} \to \mathbb{R}^{k}$ in \cref{general:target}) 
for which we already have a generic stochastic approximation algorithm 
	(cf.\ $\Psi \colon \mathfrak{P} \times \mathbb{R}^{\mathfrak{d}} \to \mathbb{R}^{k}$ and $ \mathfrak{W} \colon \Omega \to \mathbb{R}^\mathfrak{d} $ in \eqref{general:eq1} in \cref{general:target}). 
We refer to this algorithm as proposal algorithm. 
Moreover, we assume that we are able to generate random variables which approximate the target function at every point in expectation 
	(cf.\ $\mathfrak{P} \times \Omega \ni (p, \omega) \mapsto \Xi(p, \mathbf{W}(\omega)) \in \mathbb{R}^{k}$ in \eqref{general:eq2} in \cref{general:target}).

In the next few sentences we briefly sketch in words the LRV strategy in this general case.
The first step of the LRV strategy is to consider the random variables in the stochastic approximation algorithm as parameters for a parametric family of functions (corresponding to 
	$\mathfrak{P} \ni p \mapsto \Psi(p, \theta) \in \mathbb{R}^k$, 
		$\theta \in \mathbb{R}^\mathfrak{d}$, 
in \cref{general:target}); 
see \cref{general:Loss}.
The goal of the LRV strategy is then to "learn" parameters whose corresponding function yields a good approximation of the target function 
$u \colon \mathfrak{P} \to \mathbb{R}^{k}$.
Taking this into account, the second step of the LRV strategy is to minimize a loss function (cf.\ \eqref{general:eq3} in \cref{general:Loss}) measuring the distance between the approximating function and the approximating random variables with an SGD-type optimization method; see \cref{general:SGD}.
As initial guess for the SGD-type learning procedure we suggest to randomly choose the parameters according to the distribution of the random variables appearing in the proposal algorithm, since we know that this already results in a passable approximation. 
This feature of the LRV strategy is an important advantage when compared to standard deep learning methods
in the sense that the LRV strategy has already in the beginning of the training procedure a relatively small loss function.
The entire approach is presented in one single framework in \cref{general:description}.

One of the differences between this section and the previous sections is that in the previous sections we only used, for simplicity, the plain vanilla SGD method, however in this section we allow for various more sophisticated SGD-type optimization methods (cf.\ \eqref{general:eq4} in \cref{general:SGD}).
Some of these more sophisticated SGD-type methods are presented in \cref{general:SGDtype} as special cases of the framework in \cref{general:description}.


\subsection{Stochastic approximations for general target functions related to parametric expectations}
\label{general:target}
Let $ \mathfrak{p}, \mathfrak{d}, k, \mathbf{d} \in \mathbb{N} $,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let 
	$\Psi \colon \mathfrak{P} \times \mathbb{R}^{\mathfrak{d}} \to \mathbb{R}^{k}$, 
	$\Xi \colon \mathfrak{P} \times \mathbb{R}^{\mathbf{d}} \to \mathbb{R}^{k}$,
	and
	$ u \colon \mathfrak{P} \to \mathbb{R}^{k} $
be measurable,
assume for all $p \in \mathfrak{P}$ that $\mathbb{R}^{ \mathfrak{d} } \ni w \mapsto \Psi( p, w ) \in \mathbb{R}^{k}$ is continuously differentiable, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let 
$ \mathfrak{W} \colon \Omega \to \mathbb{R}^\mathfrak{d} $ and 
$ \mathbf{W} \colon \Omega \to \mathbb{R}^\mathbf{d} $ be independent random variables,
and assume for all $p \in \mathfrak{P}$ that $\mathbb{E}[\|\Xi( p, \mathbf{W})\|]<\infty$.
We think of $(\Psi, \mathfrak{W})$ as a
stochastic approximation algorithm for $u \colon \mathfrak{P} \to \mathbb{R}^k$ in the sense that for all $p \in \mathfrak{P}$ it holds that
\begin{equation}
\label{general:eq1} 
	u(p)
\approx
	\Psi( 
		p, \mathfrak{W}
	)
\end{equation}
and for every $p \in \mathfrak{P}$ we think of $\mathbb{E}[\Xi(p, \mathbf{W})]$ as a suitable approximation
\begin{equation}
\label{general:eq2} 
	u(p)
\approx
	\mathbb{E}[
		\Xi( 
			p, \mathbf{W}
		)
	]
\end{equation}
of $u(p)$.
The goal of this section is to derive an  algorithm to approximately compute the function
$
	u \colon \mathfrak{P} \to \mathbb{R}
$.


\subsection{Random loss functions for fixed random variables in stochastic approximations}
\label{general:Loss}
Let ${\bf M} \in \mathbb{N}$, 
let $\mathfrak{C} \colon \mathbb{R}^{k}\times \mathbb{R}^{k} \to [0,\infty)$
continuously differentiable,
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
let 
$ 
  W^{ m, \mathfrak{m} } \colon \Omega \to \mathbb{R}^{ \mathbf{d} } 
$, 
$ m, \mathfrak{m} \in \mathbb{N} $, 
be i.i.d.\ random variables 
which satisfy for all $ A \in \mathcal{B}( \mathbb{R}^{ d } ) $ that 
$
  \P( W^{ 0, 0 } \in A ) = \P( \mathbf{W} \in A ) 
$,
assume that 
$
  (
    P_{ 
      m, {\bf m} 
    }
  )_{
    (m, {\bf m}) \in \mathbb{N}^2
  }
$
and 
$
	 (
	 	W^{ m, \mathfrak{m} }
	 )_{ (m, \mathfrak{m}) \in \mathbb{N}^2 }
$
are independent,
and for every 
	$m \in \mathbb{N}$
let
$
F_m 
  \colon \allowbreak\mathbb{R}^{ \mathfrak{d} } 
  \times\Omega\to\mathbb{R}
$
satisfy 
for all 
$
\theta 
  \in \mathbb{R}^{ \mathfrak{d}} 
$
that
\begin{equation}
\label{general:eq3}
\begin{split} 
	F_m( \theta ) 
&=
	\frac{ 1 }{ {\bf M} }
	\Bigg[
		\sum_{ {\bf m} = 1 }^{ {\bf M} }
		\mathfrak{C}\left(
			\Xi \big( 
				P_{ 
					m, {\bf m} 
				},
				W^{m,  {\bf m}}
			\big)
			,
			\Psi(
				P_{ 
					m, {\bf m} 
				},
				\theta
			)
		\right)
	\Bigg],
\end{split}
\end{equation}
and for every 
	$m \in \mathbb{N}$
let 
$ 
  G_m 
  \colon \mathbb{R}^{\mathfrak{d}} \times \Omega 
  \to \mathbb{R}^{\mathfrak{d}} 
$
satisfy 
for all 
	$ \theta \in \mathbb{R}^{\mathfrak{d} } $, $ \omega \in \Omega $ 
that
\begin{equation}
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega ).
\end{equation}


\subsection{Learning the random variables with stochastic gradient descent-type methods}
\label{general:SGD}
Let 
$
  \psi_m \colon \mathbb{R}^{m \mathfrak{d}} \to \mathbb{R}^{\mathfrak{d}}
$, $m \in \mathbb{N}$,
be functions and
let 
$ 
  \Theta \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \mathfrak{d} } 
$ 
satisfy 
for all $ m \in \mathbb{N} $ that 
$ 
  \Theta_0 = 
    \mathfrak{W}
$
and
\begin{equation}
\label{general:eq4}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
   - 
   \psi_m((G_1(\Theta_0), G_2(\Theta_1), \dots, G_{m}(\Theta_{m-1}))).
\end{equation}
\new{
For every sufficiently large 
$m \in \mathbb{N}$
we propose to employ the random function 
$
	\mathfrak{P} \times \Omega \ni (p, \omega) \mapsto \Psi(p, \Theta_m(\omega)) \in \mathbb{R}
$
as an approximation for the target function
$
	\mathfrak{P} \ni p \mapsto u(p) \in \mathbb{R}
$.}




\subsection{Description of the proposed approximation algorithm}
\label{general:description}

\begin{algo}\label{algo:general}
Let 
	$ \mathfrak{p}, \mathfrak{d}, \mathbf{d}, k, {\bf M} \in \mathbb{N} $,
let $ \mathfrak{P} \subseteq \mathbb{R}^{ \mathfrak{p} } $ be measurable, 
let
$\Psi \colon \mathfrak{P} \times \mathbb{R}^{\mathfrak{d}} \to \mathbb{R}^{k}$,
$\Xi \colon \mathfrak{P} \times \mathbb{R}^{\mathbf{d}} \to \mathbb{R}^{k}$,
and 
$\mathfrak{C} \colon \mathbb{R}^{k} \times \mathbb{R}^{k} \to [0,\infty)$
be functions,
let 
$
  \psi_m = (\psi_m^{(1)}, \ldots, \psi_m^{(\mathfrak{d})}) \colon \mathbb{R}^{m \mathfrak{d}} \to \mathbb{R}^{\mathfrak{d}}
$, $m \in \mathbb{N}$,
be functions,
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let 
$ \mathfrak{W} \colon \Omega \to \mathbb{R}^\mathfrak{d} $ be a random variable,
let 
$  
  P_{ 
    m, {\bf m} 
  }
  \colon \Omega \to \mathfrak{P}
$, $ m, {\bf m} \in \mathbb{N} $, 
be i.i.d.\ random variables, 
let 
$ 
  W^{ m, \mathfrak{m} } \colon \Omega \to \mathbb{R}^{ \mathbf{d} } 
$, 
$ m, \mathfrak{m} \in \mathbb{N} $, 
be i.i.d.\ random variables,
assume that 
$\mathfrak{W}$,
$
  (
    P_{ 
      m, {\bf m} 
    }
  )_{
    (m, {\bf m}) \in \mathbb{N}^2
  }
$, and
$
	 (
	 	W^{ m, \mathfrak{m} }
	 )_{ (m, \mathfrak{m}) \in \mathbb{N}^2 }
$
are independent,
for every 
	$m \in \mathbb{N}$
let
$
F_m 
  \colon \allowbreak\mathbb{R}^{ \mathfrak{d} } 
  \times\Omega\to\mathbb{R}
$
satisfy 
for all 
$
\theta 
  \in \mathbb{R}^{ \mathfrak{d}} 
$
that
\begin{equation}
\label{T_B_D}
\begin{split} 
	F_m( \theta ) 
&=
	\frac{ 1 }{ {\bf M} }
	\Bigg[
		\sum_{ {\bf m} = 1 }^{ {\bf M} }
		\mathfrak{C} \! \left(
			\Xi \big( 
				P_{ 
					m, {\bf m} 
				},
				W^{m,  {\bf m}}
			\big)
			,
			\Psi(
				P_{ 
					m, {\bf m} 
				},
				\theta
			)
		\right)
	\Bigg],
\end{split}
\end{equation}
for every 
	$ m \in \mathbb{N} $ 
let 
$ 
  G_m = (G_m^{(1)},\ldots, G_m^{(\mathfrak{d})})
  \colon \mathbb{R}^{\mathfrak{d}} \times \Omega 
  \to \mathbb{R}^{\mathfrak{d}} 
$
satisfy 
for all 
	$ \omega \in \Omega $, $ \theta \in \{ v \in \mathbb{R}^{\mathfrak{d} } \colon F_m( \cdot, \omega)  \text{ is} \linebreak \text{differentiable at } v\} $
that
$
  G_m( \theta, \omega )
  =
  ( \nabla_{ \theta } F_m )( \theta, \omega )
$,
and
let 
$ 
  \Theta = (\Theta^{(1)},\ldots, \Theta^{(\mathfrak{d})}) \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^{ \mathfrak{d} } 
$ 
be a stochastic process which satisfies  
for all $ m \in \mathbb{N} $ that 
$ 
  \Theta_0 = 
    \mathfrak{W}
$
and
\begin{equation}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
   - 
   \psi_m(G_1(\Theta_0), G_2(\Theta_1), \dots, G_{m}(\Theta_{m-1})).
\end{equation}

\end{algo}




\subsection{Explicit descriptions of some popular stochastic gradient descent-type methods}
\label{general:SGDtype}

\subsubsection{Standard stochastic gradient descent (SGD)}
\label{sect:vanillaSGD}

\begin{lemma}
\label{vanillaSGD}
Assume \cref{algo:general},
let $(\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty)$,
and assume for all 	
	$m \in \mathbb{N}$,
	$g_1, g_2, \ldots, g_m \in \mathbb{R}^{\mathfrak{d}}$
that
$
	\psi_m(g_1, g_2, \ldots, g_m) = \gamma_m g_m
$.
Then it holds for all
	$ m \in \mathbb{N} $ 
that
\begin{equation}
  \Theta_{ m }
  =
  \Theta_{ m - 1 }
  - 
  \gamma_m 
  G_m( \Theta_{ m - 1 } ) 
  .
\end{equation}
\end{lemma}


\subsubsection{Stochastic gradient descent with momentum (SGD with momentum)}
\label{sect:Momentum}


\begin{lemma}
\label{Momentum}
Assume \cref{algo:general},
let $(\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty)$, $\alpha \in (0,1)$,
and assume for all 	
	$m \in \mathbb{N}$,
	$g_1, g_2, \ldots, g_m \in \mathbb{R}^{\mathfrak{d}}$
that
$
	\psi_m(g_1, g_2, \ldots, g_m) = \gamma_m \sum_{k = 1}^m \alpha^{m-k}(1-\alpha) g_k
$.
Then 
there exists $\mathbf{m} \colon \mathbb{N}_0 \times \Omega \to \mathbb{R}^d$ such that 
for all 
	$m \in \mathbb{N}$ 
it holds that 
\begin{equation}
	\mathbf{m}_0 = 0,
\qquad
	\mathbf{m}_m = \alpha \mathbf{m}_{m-1} + (1-\alpha) G_m(\Theta_{m-1}), 
\qquad\text{and}\qquad
	\Theta_n = \Theta_{n-1} - \gamma_n  \mathbf{m}_n. 
\end{equation}
\end{lemma}






\subsubsection{Adaptive stochastic gradient descent (Adagrad)}
\label{sect:Adagrad}

\begin{lemma}
\label{Adagrad}
Assume \cref{algo:general},
let 
	$(\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty)$, 
	$\varepsilon \in (0,\infty)$, 
and assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$
		g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})}), 		g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})}),
		\ldots, 
		g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\left[\frac{\gamma_m}{
		\big(\varepsilon + \sum_{k = 1}^{m} |g_k^{(i)}|^2\big)^{\nicefrac{1}{2}}
	}\right]
	g_m^{(i)}.
\end{split}
\end{equation}
Then
it holds for all $m \in \mathbb{N}$, $i \in \{1,2,\ldots,\mathfrak{d}\}$ that 
\begin{equation}
	\Theta_m^{(i)}
= 
	\Theta_{m-1}^{(i)} - 
	\left[
	\frac{\gamma_m}{\big(\varepsilon +
	\sum_{k = 1}^m |G_k^{(i)}(\Theta_{k-1})|^2
	 \big)^{\nicefrac{1}{2}}} 
	 \right]G_m^{(i)}(\Theta_{m-1}).
\end{equation}
\end{lemma}


\subsubsection{Root mean square error propagation stochastic gradient descent (RMSprop)}
\label{sect:RMSprop}


\begin{lemma}
\label{RMSprop}
Assume \cref{algo:general},
let 
	$(\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty)$, 
	$\varepsilon \in (0,\infty)$, 
	$\beta \in (0,1)$,
and assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$
		g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})}), 		g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})}),
		\ldots, 
		g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\left[\frac{\gamma_m}{
		\big(\varepsilon + \sum_{k = 1}^{m} \beta^{m-k}(1-\beta)|g_k^{(i)}|^2\big)^{\nicefrac{1}{2}}
	}\right]
	g_m^{(i)}.
\end{split}
\end{equation}
Then
there exists  
$\mathbb{M} = (\mathbb{M}^{(1)},\ldots,\mathbb{M}^{(\mathfrak{d})}) \colon \mathbb{N}_0\times \Omega \to \mathbb{R}^\mathfrak{d}$ 
such that for all $m \in \mathbb{N}$, $i \in \{1,2,\ldots,\mathfrak{d}\}$ it holds that 
\begin{equation}
	\mathbb{M}_0 = 0, 
\qquad
	\mathbb{M}_m^{(i)} 
= 
	\beta\mathbb{M}_{m-1}^{(i)} + (1-\beta)|G_m^{(i)}(\Theta_{m-1})|^2,
\end{equation}
\begin{equation}
\text{and}\qquad
	\Theta_m^{(i)}
= 
	\Theta_{m-1}^{(i)} 
	- 
	\left[ 
		\frac{\gamma_m}{(\varepsilon +\mathbb{M}_m^{(i)} )^{\nicefrac{1}{2}}} 
	\right] 
	G_m^{(i)}(\Theta_{m-1}).
\end{equation}
\end{lemma}

\subsubsection{Adadelta stochastic gradient descent (Adadelta)}
\label{sect:Adadelta}



\begin{lemma}
\label{Adadelta}
Assume \cref{algo:general},
let 
	$\varepsilon \in (0,\infty)$, 
	$\beta, \delta \in (0,1)$,
and assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$
		g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})}), 		g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})}),
		\ldots, 
		g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\left[
		\frac{
			\varepsilon + \sum_{k = 1}^{m-1} \delta^{m-1-k}(1-\delta) |\psi_k^{(i)}(g_1, g_2, \ldots, g_k)|^2
		}{
			\varepsilon + \sum_{k = 1}^{m} \beta^{m-k}(1-\beta) |g_k^{(i)}|^2
		}
	\right]^{\nicefrac{1}{2}}
	g_m^{(i)}.
\end{split}
\end{equation}
Then
there exist 
	$\mathbb{M} = (\mathbb{M}^{(1)},\ldots,\mathbb{M}^{(d)})$
	and
	$\Delta = (\Delta^{(1)}, \ldots, \Delta^{(d)}) \colon \mathbb{N}_0\times \Omega \to \mathbb{R}^d$ 
	such that for all 
	$m \in \mathbb{N}$, 
	$i \in \{1,2,\ldots,d\}$ 
it holds that 
\begin{equation}
\mathbb{M}_0 = 0, \qquad \Delta_0 = 0, 
\qquad
	\mathbb{M}_m^{(i)}
= 
	\beta \,\mathbb{M}_{m-1}^{(i)} + (1-\beta)|G_m^{(i)}(\Theta_{m-1})|^2,
\end{equation}
\begin{equation}
	\Theta_m^{(i)} 
= 
	\Theta_{m-1}^{(i)} - \bigg[\frac{\varepsilon +\Delta_{m-1}^{(i)} }{\varepsilon +\mathbb{M}_m^{(i)} }\bigg]^{\nicefrac{1}{2}} 
	G_m^{(i)}(\Theta_{m-1}), 
\qquad\text{and}\qquad 
	\Delta_m^{(i)} = \delta\Delta_{m-1}^{(i)} + (1-\delta)|\Theta_m^{(i)}-\Theta_{m-1}^{(i)}|^2. 
\end{equation}

\end{lemma}


\subsubsection{Adamax stochastic gradient descent (Adamax)}
\label{sect:Adamax}

\begin{lemma}
Assume \cref{algo:general},
let 
	$\varepsilon \in (0,\infty)$, 
	$\alpha,\beta \in (0,1)$,
and assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$
		g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})}), 		g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})}),
		\ldots, 
		g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\gamma_m
	\left[
		\frac{
			\sum_{k = 1}^m \alpha^{m-k}(1-\alpha) g_k^{(i)}
		}{
			1-\alpha^m
		}
	\right]
	\left[
		\varepsilon +
		\max\big\{|g_m^{(i)}|, \beta |g_{m-1}^{(i)}|, \ldots, \beta^{m-1} |g_1^{(i)}| \big\}
	\right]^{-1}. 
\end{split}
\end{equation}
Then
there exist 
	$\mathbf{m}  = (\mathbf{m}^{(1)},\ldots,\mathbf{m}^{(d)})$ and 
	$\mathbb{M} = (\mathbb{M}^{(1)}, \ldots, \mathbb{M}^{(d)})  \colon \mathbb{N}_0\times \Omega \to \mathbb{R}^d$ 
such that for all 
	$m \in \mathbb{N}$, 
	$i \in \{1,2,\ldots,d\}$ 
it holds that 
\begin{equation}
	\mathbf{m}_0 = 0, 
\qquad
	\mathbf{m}_m
= 
	\alpha\mathbf{m}_{m-1} + (1-\alpha) G_m(\Theta_{m-1}),
\end{equation}
\begin{equation} 
	\mathbb{M}_0 = 0, 
\qquad
	\mathbb{M}_m^{(i)} = \max\{ \beta\,\mathbb{M}_{m-1}^{(i)}, |G_m^{(i)}(\Theta_{m-1})|^2\},
\end{equation}
\begin{equation}
\text{and}\qquad 
	\Theta_m^{(i)}
= 
	\Theta_{m-1}^{(i)} 
	-  
	\gamma_m
	\left[
		\tfrac{\mathbf{m}_m^{(i)}}{1-\alpha^m}
	\right]
	\left[
		\varepsilon + \mathbb{M}_m^{(i)}
	\right]^{-1}
	.
\end{equation}
\end{lemma}


\subsubsection{Adaptive moment estimation stochastic gradient desent (Adam)}
\label{sect:Adam}

\begin{lemma}
Assume \cref{algo:general},
let 
	$\varepsilon \in (0,\infty)$, 
	$\alpha,\beta \in (0,1)$,
and assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$
		g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})}), 		g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})}),
		\ldots, 
		g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\gamma_m
	\left[
		\frac{
			\sum_{k = 1}^m \alpha^{m-k}(1-\alpha) g_k^{(i)}
		}{
			1-\alpha^m
		}
	\right]
	\left[
		\varepsilon +\bigg[
			\frac{
				\sum_{k = 1}^m \beta^{m-k}(1-\beta) |g_k^{(i)}|^2
			}{
				1- \beta^m
			}
		\bigg]^{\nicefrac{1}{2}} 
	\right]^{-1}. 
\end{split}
\end{equation}
Then
there exist 
	$\mathbf{m}  = (\mathbf{m}^{(1)},\ldots,\mathbf{m}^{(d)})$ and 
	$\mathbb{M} = (\mathbb{M}^{(1)}, \ldots, \mathbb{M}^{(d)})  \colon \mathbb{N}_0\times \Omega \to \mathbb{R}^d$ 
such that for all 
	$m \in \mathbb{N}$, 
	$i \in \{1,2,\ldots,d\}$ 
it holds that 
\begin{equation}
	\mathbf{m}_0 = 0, 
\qquad
	\mathbf{m}_m
= 
	\alpha\mathbf{m}_{m-1} + (1-\alpha) G_m(\Theta_{m-1}),
\end{equation}
\begin{equation} 
	\mathbb{M}_0 = 0, 
\qquad
	\mathbb{M}_m^{(i)} = \beta\,\mathbb{M}_{m-1}^{(i)} + (1-\beta)  |G_m^{(i)}(\Theta_{m-1})|^2,
\end{equation}
\begin{equation}
\text{and}\qquad 
	\Theta_m^{(i)}
= 
	\Theta_{m-1}^{(i)} 
	-  
	\gamma_m
	\left[
		\tfrac{\mathbf{m}_m^{(i)}}{1-\alpha^m}
	\right]
	\left[
		\varepsilon +\big( \tfrac{\mathbb{M}_m^{(i)}}{1- \beta^m}\big)^{\!\!\nicefrac{1}{2}} 
	\right]^{-1}
	.
\end{equation}
\end{lemma}






\section{Numerical examples}
\label{sect:numerics}

\newcommand{\unclear}[1]{{\color{blue}#1}}

In this section we apply the LRV strategy to different parametric approximation problems from the literature.
Specifically,
	we consider the classical parametric Black-Scholes model for the pricing of European call options in \cref{simul:BS1},
	we consider a parametric Black-Scholes model for the pricing of worst-of basket put options on three underlying assets in \cref{simul:BSworstPut},
	we consider a parametric Black-Scholes model for the pricing of average basket put options on three underlying assets with knock-in barriers in \cref{simul:BSaverageBarrier}, and
	we consider a parametric stochastic Lorentz equation in \cref{simul:Lorentz}.
In the literature there are already a number of simulation results for SGD-based deep learning methods regarding the high-dimensional pricing of financial derivative contracts.
In particular 
	we refer to, e.g., 
	\cite{berner2020numerically,BeckJafaari21,Lokeshwar22,Ferguson2018,Germain2021,Biagini2021} 
	for parametric pricing results for European options,	
	we refer to, e.g.,  \cite{Andersson21,Salvador2020a,Becker2018,becker2019pricing,Lind2022,Chen21,Lapeyre21,Gaspar20,Ye2019,Becker2019published} 
	for the pricing of American options, and 
	we refer, e.g., to \cite{ruf2020neural} for further references.

In \cref{sect:anti} we briefly recall the antithetic MC method (cf., e.g., Glasserman \cite{Glasserman03}) and some well-known properties of it. 
This is a variance reduction technique for MC methods which we will employ in some of the proposal algorithms for the LRV strategy in case of some of the above mentioned approximation problems.
In each of the considered numerical examples we also compare the LRV strategy with existing approximation techniques from the literature such as
	the deep learning method induced by Becker et al.\ \cite{BeckJafaari21}, 
	MC methods, and
	QMC methods.
	
The numerical experiments have been performed in {\sc Python} using {\sc TensorFlow} in single precision on a {\sc NVIDIA GeForce RTX 2080 Ti GPU}. The operating system was an AMD Ryzen 9 3950X CPU with 64 GB DDR4 memory running {\sc TensorFlow} 2.1 on {\sc Ubuntu} 20.04.
The source codes for the numerical experiments are provided in \cref{sect:source}.


\subsection{Antithetic Monte Carlo approximations}
\label{sect:anti}

In this section we recall a special case of antithetic variates for MC methods (cf., e.g., Glasserman \cite[Section 4.2]{Glasserman03}) when the distribution of the sampled random variables is symmetric around the origin.
The following result, \cref{lem:anti1} below, shows that the resulting antithetic MC method achieves a higher or equal $L^2$-accuracy than the standard MC method when the same number of MC samples are used for both methods.
The subsequent result, \cref{lem:anti2} below, then provides a sufficient condition for the antithetic MC method to achieve a \emph{strictly} higher $L^2$-accuracy than the standard MC method when the same number of MC samples are used for both methods.


\begin{lemma}
\label{lem:anti1}
Let $d, \mathfrak{M} \in \mathbb{N}$, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let $X_m \colon \Omega \to \mathbb{R}^d$, $m \in \{1, 2, \ldots, \mathfrak{M}\}$, be i.i.d.\ random variables,
assume for all 
	$B \in \mathcal{B}(\mathbb{R}^d)$ 
that $\P(X_1 \in B) = \P(-X_1 \in B)$,
let $f \colon \mathbb{R}^d \to \mathbb{R}$ be measurable,  
assume
$\Exp{|f(X_1)|^2} < \infty$,
and let 
	$M \colon \Omega \to \mathbb{R}$
and
	$A \colon \Omega \to \mathbb{R}$
satisfy 
\begin{equation}
\label{anti1:ass1}
\begin{split} 
	M
=
	\left[ \frac{1}{\mathfrak{M}} \sum_{m = 1}^{\mathfrak{M}} f(X_m) \right]
\qquad\text{and}\qquad
	A
=
	\frac{1}{2\mathfrak{M}} \left[ \sum_{m = 1}^{\mathfrak{M}} \bigl( f(X_m) + f(-X_m) \bigr) \right].
\end{split}
\end{equation}
Then
\begin{enumerate}[label=(\roman *)]
\item \label{anti1:item1}
it holds that 
$
	\Exp{|M-\Exp{f(X_1)}|^2}
=
	\frac{\Var{f(X_1)}}{\mathfrak{M}}
$
and

\item \label{anti1:item2}
it holds that 
$
	\Exp{|A-\Exp{f(X_1)}|^2}
=
	\frac{\Var{f(X_1)} + \Cov{f(X_1)}{f(-X_1)}}{2\mathfrak{M}}
\leq
	\frac{\Var{f(X_1)}}{\mathfrak{M}}
$.
\end{enumerate}
\end{lemma}

\begin{proof}[Proof of \cref{lem:anti1}]
First, note that the fact that $X_m$, $m \in \{1, 2, \ldots, \mathfrak{M}\}$, are i.i.d.\ implies that for all
	$m,n \in \{1, 2, \ldots, \mathfrak{M}\}$
with $m \neq n$ it holds that
\begin{equation}
\label{anti1:eq1}
\begin{split} 
	\Cov{f(X_m)}{f(X_n)} = 0
\qquad\text{and}\qquad
	\Cov{f(X_m) + f(-X_m)}{f(X_n) + f(-X_n)} = 0
\end{split}
\end{equation}
(cf., e.g., Klenke \cite[Theorem 5.4]{Klenke14}).
This, 
the fact that $\Exp{M} = \Exp{f(X_1)}$, and 
the Bienaym\'e formula (cf, e.g., Klenke \cite[Theorem 5.7]{Klenke14}) assure that
\begin{equation}
\label{anti1:eq2}
\begin{split} 
	\Exp{|M-\Exp{f(X_1)}|^2}
&=
	\Var{M}
=
	\frac{1}{\mathfrak{M}^2}
	\Var{ {\textstyle \sum_{m = 1}^{\mathfrak{M}} f(X_m)} } \\
&=
	\frac{1}{\mathfrak{M}^2}
	{\textstyle \left(
		\left[ 
			\sum_{m = 1}^{\mathfrak{M}} \Var{f(X_m)} 
		\right]
		+
			\sum_{m,n \in \{1, 2, \ldots, \mathfrak{M}\}, m \neq n}
				\Cov{f(X_m)}{f(X_n)}
	\right)}\\
&=
	\frac{1}{\mathfrak{M}^2} \,
	{\textstyle
			\sum_{m = 1}^{\mathfrak{M}} \Var{f(X_1)} 
	}
=
	\frac{\Var{f(X_1)}}{\mathfrak{M}}.
\end{split}
\end{equation}
This proves \cref{anti1:item1}. 
Next observe that the fact that $X_1,X_2, \ldots, X_\mathfrak{M}, -X_1, -X_2, \ldots, -X_\mathfrak{M}$ are identically distributed implies that 
\begin{equation}
\label{anti1:eq3}
\begin{split} 
	\Exp{A}
=
	\frac{1}{2\mathfrak{M}} \sum_{m = 1}^{\mathfrak{M}} \left(\Exp{f(X_m)} + \Exp{f(-X_m)}\right)
=
	\frac{1}{2\mathfrak{M}} \sum_{m = 1}^{\mathfrak{M}} 2\Exp{f(X_1)}
=
	\Exp{f(X_1)}.
\end{split}
\end{equation}
This,
\eqref{anti1:eq1},
the fact that $X_1,X_2, \ldots, X_\mathfrak{M}, -X_1, -X_2, \ldots, -X_\mathfrak{M}$ are identically distributed, and
the Bienaym\'e formula (cf, e.g., Klenke \cite[Theorem 5.7]{Klenke14})
assure
that 
\begin{equation}
\label{anti1:eq4}
\begin{split} 
	\Exp{|A-\Exp{f(X_1)}|^2}
&=
	\Var{A}
=
	\frac{1}{(2\mathfrak{M})^2}
	\Var{ {\textstyle \sum_{m = 1}^{\mathfrak{M}} f(X_m) + f(-X_m)} } \\
&=
	\frac{1}{4\mathfrak{M}^2}
	{\textstyle \bigg(
		\left[ 
			\sum_{m = 1}^{\mathfrak{M}} \Var{f(X_m) + f(-X_m)} 
		\right]}\\
&\quad
		+
		{\textstyle
		\left[	
			\sum_{m,n \in \{1, 2, \ldots, \mathfrak{M}\}, m \neq n}
				\Cov{f(X_m) + f(-X_m)}{f(X_n) + f(-X_n)}
		\right]
	\bigg)}\\
&=
	\frac{1}{4\mathfrak{M}^2} \,
	{\textstyle 			\sum_{m = 1}^{\mathfrak{M}} \Var{f(X_1) + f(-X_1)} 
	} \\
&=
	\frac{\Var{f(X_1) + f(-X_1)} }{4\mathfrak{M}} \\
&=
	\frac{\Var{f(X_1)} + 2 \Cov{f(X_1)}{f(-X_1)} + \Var{f(-X_1)}}{4\mathfrak{M}} \\
&=
	\frac{\Var{f(X_1)} + \Cov{f(X_1)}{f(-X_1)}}{2\mathfrak{M}}.
\end{split}
\end{equation}
Moreover, note that 
the fact that $X_1$ and $-X_1$ are identically distributed and 
the Cauchy-Schwarz inequality (cf., e.g., Klenke \cite[Theorem 5.8]{Klenke14}) 
assure that
\begin{equation}
\label{anti1:eq5}
\begin{split} 
	(\Cov{f(X_1)}{f(-X_1)})^2
\leq
	\Var{f(X_1)}\Var{f(-X_1)}
=
	(\Var{f(X_1)})^2.
\end{split}
\end{equation}
Combining this with \eqref{anti1:eq4} demonstrates that 
\begin{equation}
\label{anti1:eq6}
\begin{split} 
	\Exp{|A-\Exp{f(X_1)}|^2}
&=
	\frac{\Var{f(X_1)} + \Cov{f(X_1)}{f(-X_1)}}{2\mathfrak{M}}\\
&\leq
	\frac{\Var{f(X_1)} + \Var{f(X_1)}}{2\mathfrak{M}}
=
	\frac{\Var{f(X_1)}}{\mathfrak{M}}.
\end{split}
\end{equation}
This establishes \cref{anti1:item2}.
The proof of \cref{lem:anti1} is thus complete.
\end{proof}




\begin{corollary}
\label{lem:anti2}
Let $d, \mathfrak{M} \in \mathbb{N}$, 
let $ ( \Omega, \mathcal{F}, \P ) $ be a probability space, 
let $X_m \colon \Omega \to \mathbb{R}^d$, $m \in \{1, 2, \ldots, \mathfrak{M}\}$, be i.i.d.\ random variables,
assume for all 
	$B \in \mathcal{B}(\mathbb{R}^d)$ 
that $\P(X_1 \in B) = \P(-X_1 \in B)$,
let $f \colon \mathbb{R}^d \to \mathbb{R}$ be measurable,  
assume
for all 
	$a,b,c \in \mathbb{R}$ 
that $\Exp{|f(X_1)|^2} < \infty$
and
$
	\P(af(X_1) + b f(-X_1) + c = 0) < 1
$,
and let 
	$M \colon \Omega \to \mathbb{R}$
and
	$A \colon \Omega \to \mathbb{R}$
satisfy 
\begin{equation}
\label{anti2:ass1}
\begin{split} 
	M
=
	\left[ \frac{1}{\mathfrak{M}} \sum_{m = 1}^{\mathfrak{M}} f(X_m) \right]
\qquad\text{and}\qquad
	A
=
	\frac{1}{2\mathfrak{M}} \left[ \sum_{m = 1}^{\mathfrak{M}} \bigl( f(X_m) + f(-X_m) \bigr) \right].
\end{split}
\end{equation}
Then
\begin{enumerate}[label=(\roman *)]
\item \label{anti2:item1}
it holds that 
$
	\Exp{|M-\Exp{f(X_1)}|^2}
=
	\frac{\Var{f(X_1)}}{\mathfrak{M}}
$
and

\item \label{anti2:item2}
it holds that 
$
	\Exp{|A-\Exp{f(X_1)}|^2}
=
	\frac{\Var{f(X_1)} + \Cov{f(X_1)}{f(-X_1)}}{2\mathfrak{M}}
<
	\frac{\Var{f(X_1)}}{\mathfrak{M}}
$.
\end{enumerate}
\end{corollary}

\begin{proof}[Proof of \cref{lem:anti2}]
First, note \cref{anti1:item1} in \cref{lem:anti1} establishes \cref{anti2:item1}.
Moreover, observe that 
the assumption that 
for all 
	$a,b,c \in \mathbb{R}$ 
it holds that
$
	\P(af(X_1) + b f(-X_1) + c = 0) < 1
$,
the fact that $X_1$ and $-X_1$ are identically distributed, and 
the strict Cauchy-Schwarz inequality (cf., e.g., Klenke \cite[Theorem 5.8]{Klenke14}) 
assure that
\begin{equation}
\label{anti2:eq5}
\begin{split} 
	(\Cov{f(X_1)}{f(-X_1)})^2
<
	\Var{f(X_1)}\Var{f(-X_1)}
=
	(\Var{f(X_1)})^2.
\end{split}
\end{equation}
Combining this with \cref{anti1:item2} in \cref{lem:anti1} demonstrates that 
\begin{equation}
\label{anti2:eq6}
\begin{split} 
	\Exp{|A-\Exp{f(X_1)}|^2}
&=
	\frac{\Var{f(X_1)} + \Cov{f(X_1)}{f(-X_1)}}{2\mathfrak{M}}\\
&<
	\frac{\Var{f(X_1)} + \Var{f(X_1)}}{2\mathfrak{M}}
=
	\frac{\Var{f(X_1)}}{\mathfrak{M}}.
\end{split}
\end{equation}
This establishes \cref{anti2:item2}.
The proof of \cref{lem:anti2} is thus complete.
\end{proof}






\subsection[Parametric Black-Scholes equations for European call options]{Parametric Black-Scholes partial differential equations for European call options}
\label{simul:BS1}

\newcommand{8\,192\,000}{8\,192\,000}
\newcommand{140000}{140000}




In this section we apply the LRV strategy to the problem of approximating the fair price of an European call option in the classical Black-Scholes model (cf.\ Black \& Scholes \cite{BlackScholes73} and Merton \cite{Merton73}).
A brief summary of the numerical results of this subsection can be found in \cref{table:BS1_Summary} in the introduction.
We start by introducing the Black-Scholes model in the context of \cref{algo:general}.

Assume \cref{algo:general},
assume
	$\mathfrak{p} = 5$
and
\begin{equation}
\begin{split} 
	\mathfrak{P} = {[90, 110] \times [0.01,1] \times [-0.1,0.1] \times [0.01,0.5] \times  [90, 110]},
\end{split}
\end{equation}
let $\mathfrak{N} \colon \mathbb{R}\to \mathbb{R}$ satisfy for all $z\in \mathbb{R}$ that
$
  \mathfrak{N}(z) =\tfrac{1}{\sqrt{2\pi}} \int_{-\infty}^z \exp({-\frac{y^2}{2}}) \d y
$,
and let $u \colon \mathfrak{P} \to \mathbb{R}$ satisfy for all
	$p = (\xi, T, r, \sigma, K) \in \mathfrak{P}$
that
\begin{equation}
\label{BS1:eq0}
\begin{split} 
	u(p)
=
	\xi \,
	\mathfrak{N} \Big(
		\tfrac{ (r + \frac{ \sigma^2 }{2}) T + \ln( {\xi}/{K} )}{ \sigma \sqrt{T} }
	\Big)
	-
	\exp({-r T})K
	\,
	\mathfrak{N} \Bigl(
		\tfrac{ (r - \frac{ \sigma^2 }{2}) T + \ln( {\xi}/{K} )}{ \sigma \sqrt{T} }
	\Bigr).
\end{split}
\end{equation}
Note that \cref{BS1:eq0} corresponds to the famous Black-Scholes formula for European call options.
In the economic interpretation of the Black-Scholes model, for every 
	$p = (\xi, T, r, \sigma, K) \in \mathfrak{P}$
the number 
	$u(p) \in \mathbb{R}$
thus corresponds to the fair price of a European call option with
	initial price $\xi$,	
	time of maturity $T$,
	drift rate $r$,
	volatility $\sigma$,
	and strike price $K$.

We now specify the mathematical objects in the LRV strategy appearing in \cref{algo:general} 
to approximately calculate the target function $u \colon \mathfrak{P} \to \mathbb{R}$ in \cref{BS1:eq0}.
Specifically, in addition to the assumptions above,
let 
	$\mathfrak{M} \in \mathbb{N}$,
	$a \in \{0,1\}$,
assume
	$\mathfrak{d} =  \mathfrak{M}$,
	$\mathbf{d} =  k = 1$,
	${\bf M} = {8192}$,
for every 
	$\xi, r, \sigma, w \in \mathbb{R}$,
	$T \in [0,\infty)$
let $X^{r, \sigma, \xi, w}_T \in \mathbb{R}$
satisfy
\begin{equation}
\begin{split} 
	X^{r, \sigma, \xi, w}_T
=
	\xi \exp \! \big(
		(r - \tfrac{\sigma^2}{2})T + \sigma \sqrt{T} w
	\big),
\end{split}
\end{equation}
let 
	$ \phi_\mathfrak{a} \colon \mathfrak{P} \times \mathbb{R}^{ \mathbf{d} } \to \mathbb{R} $, 
		$\mathfrak{a} \in \{0, 1\}$,
satisfy for all
	$p = (\xi, T, r, \sigma, K) \in \mathfrak{P}$,
	$w \in \mathbb{R}$
that
\begin{equation}
\label{BS1:eq0.1}
\begin{split} 
	\phi_0(p, w) 
=
	\exp({-r T})\max \{ X^{r, \sigma, \xi, w}_T - K,0\}
\qquad\text{and}
\end{split}
\end{equation}
\begin{equation}
\label{BS1:eq0.2}
\begin{split} 
\textstyle
	\phi_1(p, w) 
=
	\frac{\phi_0(p, w) + \phi_0(p, -w)}{2}
=
	\frac{\exp({-r T})}{2}(\max \{ X^{r, \sigma, \xi, w}_T - K,0\} + \max \{ X^{r, \sigma, \xi, -w}_T - K,0\}),
\end{split}
\end{equation}
assume for all 
	$p = (\xi, T, r, \sigma, K) \in \mathfrak{P}$,
	$w \in \mathbb{R}$,
	$ \theta = (\theta_1, \ldots, \theta_{\mathfrak{M}}) \in \mathbb{R}^{\mathfrak{M}}$
that
$
	\Xi(p, w)
=
	u(p)
$ 
and
\begin{equation}
\label{BS1:eq1}
\begin{split} 
	\Psi(p, \theta)
=
	\frac{ 1 }{ \mathfrak{M} }
	\textstyle
	\left[ 
		\sum\limits_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
			\phi_a( 
				p, 
				\theta_{ \mathfrak{m}}         
			)
	\right],
\end{split}
\end{equation}
assume for all
	$x, y \in \mathbb{R}$
that
	$\mathfrak{C}(x, y) = | x-y |^2$,
let
	$\mathbf{a} = {\frac{9}{10}}$,
	$\mathbf{b}  = {\frac{999}{1000}}$,
	$ \varepsilon \in (0,\infty)$,
let
	$(\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty)$
satisfy for all 
	$j \in \{1, 2, \ldots, 7\}$,
	$m \in \mathbb{N} \cap (20000 (j-1),  20000 j]$
that
$
	{\gamma_m = 10^{-j}}
$,
assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})})$, 	$g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})})$,
	$\dots$, 
	$g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\label{BS:eq2}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\gamma_m
	\left[
		\frac{
			\sum_{k = 1}^m \mathbf{a}^{m-k}(1-\mathbf{a}) g_k^{(i)}
		}{
			1-\mathbf{a}^m
		}
	\right]
	\left[
		\varepsilon +\bigg[
			\frac{
				\sum_{k = 1}^m \mathbf{b}^{m-k}(1-\mathbf{b}) |g_k^{(i)}|^2
			}{
				1- \mathbf{b}^m
			}
		\bigg]^{\nicefrac{1}{2}} 
	\right]^{-1},
\end{split}
\end{equation}	
assume that  
$ \mathfrak{W}$ is a standard normal random vector,
assume that  
$  
  P_{ 
    1, 1 
  }
$ is $\mathcal{U}_{\mathfrak{P}}$-distributed, and
assume that 
$ 
  W^{ 1, 1 }
$
is a standard normal random variable.

Let us add some comments regarding the setup introduced above.
Observe that 
	\eqref{BS1:eq0}, 
	\eqref{BS1:eq0.1}, 
	\eqref{BS1:eq0.2}, and 
	\cref{cor:black_scholes} 
assure that for all 
	$p = (\xi, T, r, \sigma, K) \in \mathfrak{P}$,
it holds that
\begin{equation}
\begin{split} 
	u(p)
&=
	\EXPP{
		\exp({-r T})\max \{ X^{r, \sigma, \xi, W^{1, 1}}_T - K,0\}
	}
=
	\Exp{
		\phi_0(p, W^{1, 1})
	}
=
	\Exp{
		\phi_1(p, W^{1, 1})
	}.
\end{split}
\end{equation}
Moreover, note that in the case $a = 0$ the proposal algorithm in \eqref{BS1:eq1} corresponds to the standard MC method and note that
in the case $a = 1$ the proposal algorithm in \eqref{BS1:eq1} corresponds to the antithetic MC method (cf.\ \cref{sect:anti}).
Furthermore, observe that \cref{BS:eq2} describes the Adam optimizer in the setup of \cref{algo:general} (cf.\ Kingma \& Ba \cite{Kingma2014} and \cref{sect:Adam}).

In \cref{table:BS_LRV} we approximately present for
	$\mathfrak{M} \in \{2^5, 2^6, \ldots, 2^{15}\}$,
	$a \in \{0,1\}$
one random realization of the $L^1(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BS_lrv:eq1}
\begin{split} 
\textstyle
	\int_\mathfrak{P} 
		|u(p) - \Psi(p, \Theta_{140000})|
	\d p
\end{split}
\end{equation}
(\nth{3} column in \cref{table:BS_LRV}),
one random realization of the $L^2(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BS_lrv:eq2}
\begin{split} 
\textstyle
	\bigl[
		\int_\mathfrak{P} 
			|u(p) - \Psi(p, \Theta_{140000})|^2
		\d p
	\bigr]^{\nicefrac{1}{2}}
\end{split}
\end{equation}
(\nth{4} column in \cref{table:BS_LRV}),
one random realization of the $L^\infty(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BS_lrv:eq3}
\begin{split} 
\textstyle
	\sup_{p \in \mathfrak{P}}
		|u(p) - \Psi(p, \Theta_{140000})|,
\end{split}
\end{equation}
(\nth{5} column in \cref{table:BS_LRV}),
the time to compute $\Theta_{140000}$
(\nth{6} column in \cref{table:BS_LRV}),
and the time to compute 
	$8\,192\,000$ 
evaluations of the function
$
	\mathfrak{P} \ni p \mapsto \Psi(p, \Theta_{140000}) \in \mathbb{R}
$
(\nth{7} column in \cref{table:BS_LRV}).
For every 	
	$\mathfrak{M} \in \{2^5, 2^6, \ldots, 2^{15}\}$,
	$a \in \{0,1\}$
we approximated the integrals in \eqref{BS_lrv:eq1} and \eqref{BS_lrv:eq2} with the MC method based on $8\,192\,000$ samples
and we approximated the supremum in \eqref{BS_lrv:eq3} based on $8\,192\,000$ random samples 
(cf., e.g., Beck et al.\ \cite[Lemma 3.5]{BeckJafaari21} and Beck et al.\ \cite[Section 3.3]{Beck2019}).

To compare the LRV strategy with existing approximation techniques from the literature, we also employ several other methods to approximate the function $u \colon \mathfrak{P} \to \mathbb{R}$ in \eqref{BS1:eq0}.
Specifically, 
	in \cref{table:BS_ANNs} we present numerical simulations for the deep learning method induced by Becker et al.\ \cite{BeckJafaari21}
	(with training values given by the exact solution,
	Adam $140000$ training steps, 
	batch size $8192$, and 
	learning rate schedule 
	$
		\mathbb{N} \ni j 
		\mapsto 
		\sum_{j = 1}^7
			\mathbbm{1}_{
				(20000 (j-1),  20000 j]
			}(m)
			10^{-j}
	$),
	in \cref{table:BS_MC} we present numerical simulations for the standard and the antithetic MC method, and
	in \cref{table:BS_QMC} we present numerical simulations for the standard and the antithetic QMC method with Sobol sequences.
In \cref{table:BS_ANNs,table:BS_MC,table:BS_QMC} 
	we have approximated the $L^1(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors 
		of the respective approximation methods with the MC method based on $8\,192\,000$ samples,
	we have approximated the $L^2(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the respective approximation methods with the MC method based on $8\,192\,000$ samples, and 
	we have approximated the $L^\infty(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the respective approximation methods based on $8\,192\,000$ random samples
		(cf., e.g., Beck et al.\ \cite[Lemma 3.5]{BeckJafaari21} and Beck et al.\ \cite[Section 3.3]{Beck2019}).
		







\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|c|,
     table head=
     \hline $\mathfrak{M}$ &  $a$ & \thead{Number \\of \\trainable \\parameters}
     & \thead{$L^1$-approx.\\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} & \thead{Training \\ time \\ in seconds} &\thead{Evaluation \\ time  for \\ $8\,192\,000$\\ evaluations \\ in seconds} \\\hline,
    late after line=\\\hline]{Table_2.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	antithetic=\anti, 
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	train-time=\train, 
	eval-time=\eval  ,
	dtype = \dtype
}
{\nsamp& \anti& \nsamp&  \lll & \llll & \linf &\train&\eval }}
\caption{\label{table:BS_LRV}Numerical simulations for the LRV strategy in case of the Black-Scholes model for European call options on one underlying described in \cref{simul:BS1} (5-dimensional approximation problem)}
\end{table}




\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|c|c|,
     table head=
     \hline \thead{Number \\ of \\ hidden \\ layers}& \thead{Number of \\ neurons \\ on  each \\ hidden \\ layer}& \thead{Number \\ of \\ trainable \\ parameters} & \thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} & \thead{Training \\ time \\ in \\ seconds} &\thead{Evaluation \\ time  for \\ $8\,192\,000$\\ evaluations \\ in seconds} \\\hline,
    late after line=\\\hline]{Table_3.csv}
{
	layers = \lyrs,
	inner-dim = \inner,
	num-weights = \weights,
	mc-samples = \mcsamp, 
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	train-time=\train, 
	eval-time=\eval  ,
	dtype = \dtype
}
{\pgfmathparse{int(\lyrs )}\pgfmathresult & \inner &\weights & \lll & \llll & \linf &\train & \eval }}
\caption{\label{table:BS_ANNs}Numerical simulations for the deep learning method induced by Becker et al.\ \cite{BeckJafaari21} in case of the Black-Scholes model European call options on one underlying described in \cref{simul:BS1} (5-dimensional approximation problem)}
\end{table}




\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|,
     table head=
     \hline \thead{Number \\ of \\MC \\samples}& \thead{MC Method \\ 0: standard \\ 1: antithetic} & \thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} &\thead{Evaluation \\ time  for \\ $8\,192\,000$\\ evaluations \\ in seconds} \\\hline,
    late after line=\\\hline]{Table_4.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	antithetic=\anti, 
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	time=\train, 
	dtype = \dtype
}
{\nsamp&\anti& \lll & \llll & \linf &\train }}
\caption{\label{table:BS_MC}Numerical simulations for the standard and the antithetic MC method in case of the Black-Scholes model for European call options on one underlying described in \cref{simul:BS1} (5-dimensional approximation problem)}
\end{table}








\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|,
     table head=
     \hline 
     	\thead{Number \\of \\MC \\samples}& \thead{QMC Method \\ 0: standard \\ 1: antithetic} & \thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} &\thead{Evaluation \\ time for \\  $8\,192\,000$\\ evaluations \\ in seconds} 
     \\\hline,
    late after line=\\\hline]{Table_5.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	antithetic=\anti, 
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	time=\train, 
	dtype = \dtype
}
{\nsamp&\anti& \lll & \llll & \linf &\train }}
\caption{\label{table:BS_QMC}Numerical simulations for the standard and the antithetic QMC method with Sobol sequences in case of the Black-Scholes model for European call options on one underlying described in \cref{simul:BS1} (5-dimensional approximation problem)}
\end{table}




\subsection[Parametric Black-Scholes equations for multi-asset worst-of basket put options]{Parametric Black-Scholes partial differential equations for multi-asset worst-of basket put options}
\label{simul:BSworstPut}

\newcommand{128000}{128000}
\newcommand{52428800}{52428800}
\newcommand{4096}{4096}
\newcommand{30000}{30000}



\newcommand{Q}{Q}
\newcommand{L}{L}
\newcommand{\mathcal{B}}{\mathcal{B}}
\newcommand{\mathbf{l}}{\mathbf{l}}




In this section we apply the LRV strategy for the approximation of the price of European worst-of basket put options in the Black-Scholes option pricing model (cf.\ Black \& Scholes \cite{BlackScholes73} and Merton \cite{Merton73}).
We start by introducing the Black-Scholes model for this pricing problem in the context of \cref{algo:general}.

Assume \cref{algo:general},
let 
	$d = 3$, 
	$e_1 = (1, 0, 0)$, 
	$e_2 = (0, 1, 0)$, 
	$e_3 = (0, 0, 1)$,
let $Q \colon \mathbb{R}^d \to \mathbb{R}^{d \times d}$ satisfy for all 
	$\rho = (\rho_1, \rho_2, \rho_3) \in \mathbb{R}^d$
that
\begin{equation}
\begin{split} 
	Q(\rho) 
= 
	\begin{pmatrix}
		1 &\rho_1 &\rho_2 \\
		\rho_1 &1 &\rho_3 \\
		\rho_2 &\rho_3 &1
	\end{pmatrix},
\end{split}
\end{equation}
let 
$
	R
\subseteq 
	\mathbb{R}^d
$
satisfy
$
	R 
= 
	\bigl\{
		\rho = (\rho_1, \rho_2, \rho_3) \in [-\frac{95}{100},\frac{95}{100}]^d 
		\colon  
		1 - |\rho_2|^2 - \frac{(\rho_3 - \rho_1\rho_2)^2}{(1-|\rho_1|^2)^{1/2}} \geq 0
	\bigr\}
$,
assume
	$\mathfrak{p} = 15$
and
\begin{equation}
\begin{split} 
	\mathfrak{P} 
= 
	[90, 110]^3 \times
	[\tfrac{1}{100}, 1] \times
	[-\tfrac{1}{20}, \tfrac{1}{20}] \times
	[0, \tfrac{1}{10}]^3 \times
	[\tfrac{1}{100}, \tfrac{1}{2}]^3 \times
	R \times 
	[90, 110],
\end{split}
\end{equation}
let $L \colon R \to \mathbb{R}^{d \times d}$ satisfy for all 
	$\rho = (\rho_1, \rho_2, \rho_3) \in R$
that
\begin{equation}
\begin{split} 
	L(\rho)
=
	\begin{pmatrix}
		1 & 0 & 0 \\
		\rho_1& (1 - |\rho_1|^2)^{{1}/{2}} &0 \\
		\rho_2& \frac{\rho_3 - \rho_1\rho_2}{(1 - |\rho_1|^2)^{{1}/{2}}} & \big(1 - |\rho_2|^2 - \frac{(\rho_3 - \rho_1\rho_2)^2}{1 - |\rho_1|^2}\big)^{{1}/{2}}
	\end{pmatrix},
\end{split}
\end{equation}
for every
	$\xi = (\xi_1, \xi_2, \xi_3) \in \mathbb{R}^d$,
	$T \in [0,\infty)$,
	$r \in \mathbb{R}$,
	$\delta = (\delta_1, \delta_2, \delta_3) \in \mathbb{R}^d$,
	$\sigma = (\sigma_1, \sigma_2, \sigma_3) \in \mathbb{R}^d$,
	$\rho \in R$,
	$w \in \mathbb{R}^d$
let 
$ 
	X^{\xi, r, \delta, \sigma, \rho, w}_T
= 
	(
	X^{r, \xi, \delta, \sigma, \rho, w, 1}_T, 
	X^{r, \xi, \delta, \sigma, \rho, w, 2}_T, 
	X^{r, \xi, \delta, \sigma, \rho, w, 3}_T
	) 
	\in \mathbb{R}^{ d }
$ 
satisfy\footnote{Note that for all
	$n \in \mathbb{N}$,
	$v= (v_1, \ldots, v_n)$, 
	$w = (w_1, \ldots, w_n) \in \mathbb{R}^n$
it holds that 
$
	\langle v, w \rangle = \sum_{i = 1}^n v_i w_i
$.
} for all
	$i \in \{1, 2, 3\}$
that
\begin{equation}
\label{BSworstPut:eq1}
\begin{split} 
	X^{\xi, r, \delta, \sigma, \rho, w, i}_T
=
	\xi_i 
	\exp\Bigl(
		\bigl[r - \delta_i -  \tfrac{|\sigma_i|^2}{2}\bigr]T
		+
		\sqrt{T}
		\sigma_i 
		\langle 
			L(\rho) w
			,
			e_i
		\rangle 
	\Bigr),
\end{split}
\end{equation}
let 
	$\mathbf{l} \colon \mathbb{R}^d \to \mathbb{R}$ 
satisfy for all 
	$x = (x_1, x_2, x_3) \in \mathbb{R}^d$
that
$\mathbf{l}(x) = \min\{ x_1, x_2, x_3 \}$,
let 
$
	\phi \colon \mathfrak{P} \times \mathbb{R}^d \to \mathbb{R}
$
satisfy for all
	$p = (\xi, T, r, \delta, \sigma, \rho, K) \in \mathfrak{P}$,
	$w \in \mathbb{R}^d$
that
\begin{equation}
\label{BSworstPut:eq2.2}
\begin{split} 
	\phi(p, w)
=
	\exp(-r T)
	\max\bigl\{
		K - 
		\mathbf{l}(
			X^{\xi, r, \delta, \sigma, \rho, w}_T
		)
		,
		0
	\bigr\},
\end{split}
\end{equation}
let $\mathcal{W} \colon \Omega \to \mathbb{R}^{ d } $ be a standard normal random vector,
and let $u \colon \mathfrak{P} \to \mathbb{R}$ satisfy for all
	$p \in \mathfrak{P}$
that
\begin{equation}
\label{BSworstPut:eq3}
\begin{split} 
	u(p)
=
	\Exp{
		\phi(p, \mathcal{W})
	}.
\end{split}
\end{equation}
In the economic interpretation of the Black-Scholes model for every 
	$p = (\xi, T, r, \delta, \sigma, \rho, K) \in \mathfrak{P}$
the number $u(p) \in \mathbb{R}$ corresponds to the fair price of a worst-of basekt put option on three underlying assets with 
	initial prices $\xi$,
	time of maturity $T$,
	risk free rate $r$,
	dividend yields of the respective underlying assets $\delta$,
	volatilities of the respective underlying assets $\sigma$,
	covariance matrix of the Brownian motions $Q(\rho)$,
	and
	strike price $K$.

We now specify the mathematical objects in the LRV strategy appearing in \cref{algo:general} 
to approximately calculate the target function $u \colon \mathfrak{P} \to \mathbb{R}$ in \cref{BSworstPut:eq3}.
Specifically, in addition to the assumptions above, 
let 
	$\mathfrak{M} \in \mathbb{N}$,
	$\mathscr{M} =  4096$,
assume
	$\mathfrak{d} =  \mathfrak{M} d$,
	$\mathbf{d} =  \mathscr{M} d$,
	$k = 1$,
	${\bf M} = {1024}$,
assume for all 
	$ p = (\xi, T, r, \delta, \sigma, \rho, K) \in \mathfrak{P} $,
	$ w = (w_1,\ldots, w_{\mathscr{M} d}) \in \mathbb{R}^{\mathscr{M} d}$,
	$ \theta = (\theta_1,\ldots, \theta_{\mathfrak{M} d}) \in \mathbb{R}^{\mathfrak{M} d}$
that
\begin{equation}
\label{BSworstPut:eq7}
\begin{split} 
	\Xi(p, w)
=
	\frac{ 1 }{ \mathscr{M} }
	\left[ 
		\textstyle \sum\limits_{ \mathfrak{m} = 1 }^{ \mathscr{M} }
			\phi\big( 
				p, 
				( 
				w_{ ( \mathfrak{m} - 1 )d + k }
				)_{k \in \{1, 2, 3 \}}         
			\big)
	\right]
\,\,\,\text{and}\,\,\,
	\Psi(p, \theta)
=
	\frac{ 1 }{ \mathfrak{M} }
	\left[ 
		\textstyle \sum\limits_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
			\phi\big( 
				p, 
				( 
				\theta_{ ( \mathfrak{m} - 1 )d + k }
				)_{k \in \{1, 2, 3 \}}         
			\big)
	\right],
\end{split}
\end{equation}
assume for all
	$x, y \in \mathbb{R}$
that
	$\mathfrak{C}(x, y) = | x-y |^2$,
let
	$\mathbf{a} = {\frac{9}{10}}$,
	$\mathbf{b}  = {\frac{999}{1000}}$,
	$ \varepsilon \in (0,\infty)$,
let
	$(\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty)$
satisfy for all 
	$m \in \{1, 2, \ldots, 30000 \}$
that
$
	\gamma_m 
=  
	{\mathbbm{1}_{(0,15000]}(m)}{10^{-3}}
	+
	{\mathbbm{1}_{(15000,25000]}(m)}{10^{-4}}
	+
	{\mathbbm{1}_{(25000,30000]}(m)}{10^{-5}}
$,
assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})})$, 	$g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})})$,
	$\dots$, 
	$g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\label{BSworstPut:eq7.2}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\gamma_m
	\left[
		\frac{
			\sum_{k = 1}^m \mathbf{a}^{m-k}(1-\mathbf{a}) g_k^{(i)}
		}{
			1-\mathbf{a}^m
		}
	\right]
	\left[
		\varepsilon +\bigg[
			\frac{
				\sum_{k = 1}^m \mathbf{b}^{m-k}(1-\mathbf{b}) |g_k^{(i)}|^2
			}{
				1- \mathbf{b}^m
			}
		\bigg]^{\nicefrac{1}{2}} 
	\right]^{-1},
\end{split}
\end{equation}	
assume that  
$ \mathfrak{W}$ is a standard normal random vector,
assume that  
$  
  P_{ 
    1, 1 
  }
$ is $\mathcal{U}_{\mathfrak{P}}$-distributed,
and
assume that 
$ 
  W^{ 1, 1 }
$
is a standard normal random vector.

Let us add some comments regarding the setup introduced above.
Note that the proposal algorithm in \cref{BSworstPut:eq3} corresponds to the MC method. 
Furthermore, observe that \cref{BSworstPut:eq7.2} describes the Adam optimizer in the setup of \cref{algo:general} (cf.\ Kingma \& Ba \cite{Kingma2014} and \cref{sect:Adam}).

In \cref{table:BSworstPut_LRV,table:BSworstPut_LRV_self} we approximately present for
	$\mathfrak{M} \in \{2^5, 2^6, \ldots, 2^{13}\}$
one random realization of the $L^1(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BSworstPut_LRV:eq1}
\begin{split} 
\textstyle
	\int_\mathfrak{P} 
		|u(p) - \Psi(p, \Theta_{30000})|
	\d p
\end{split}
\end{equation}
(\nth{3} column in \cref{table:BSworstPut_LRV,table:BSworstPut_LRV_self}),
one random realization of the $L^2(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BSworstPut_LRV:eq2}
\begin{split} 
\textstyle
	\bigl[
		\int_\mathfrak{P} 
			|u(p) - \Psi(p, \Theta_{30000})|^2
		\d p
	\bigr]^{\nicefrac{1}{2}}
\end{split}
\end{equation}
(\nth{4} column in \cref{table:BSworstPut_LRV,table:BSworstPut_LRV_self}),
one random realization of the $L^\infty(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BSworstPut_LRV:eq3}
\begin{split} 
\textstyle
	\sup_{p \in \mathfrak{P}}
		|u(p) - \Psi(p, \Theta_{30000})|,
\end{split}
\end{equation}
(\nth{5} column in \cref{table:BSworstPut_LRV,table:BSworstPut_LRV_self}), and
the time to compute $\Theta_{30000}$
(\nth{6} column in \cref{table:BSworstPut_LRV,table:BSworstPut_LRV_self}).
For every 	
	$\mathfrak{M} \in \{2^5, 2^6, \ldots, 2^{13}\}$
we approximated the integrals in \eqref{BSworstPut_LRV:eq1} and \eqref{BSworstPut_LRV:eq2} with the MC method based on $128000$ samples
and we approximated the supremum in \eqref{BSworstPut_LRV:eq3} based on $128000$ random samples 
(cf., e.g., Beck et al.\ \cite[Lemma 3.5]{BeckJafaari21} and Beck et al.\ \cite[Section 3.3]{Beck2019}).
In our approximations of \eqref{BSworstPut_LRV:eq1}, \eqref{BSworstPut_LRV:eq2}, and \eqref{BSworstPut_LRV:eq3} we have approximately computed for all required sample points 
	$p \in \mathfrak{P}$
the value $u(p)$ of the unknown exact solution 
by means of an MC approximation with $ 52428800$ MC samples in \cref{table:BSworstPut_LRV} and
by means of the LRV strategy itself (i.e., $\Psi(p, \Theta_{30000})$) with $\mathfrak{M} = 2^{15} = 32768$ in \cref{table:BSworstPut_LRV_self}. 


Besides the LRV strategy we also employed the standard MC method to approximate the function $u \colon \mathfrak{P} \to \mathbb{R}$ in \eqref{BSworstPut:eq3}. 
In \cref{table:BSworstPut_MC} we present the corresponding numerical simulation results.
In \cref{table:BSworstPut_MC}
	we have approximated the $L^1(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
		of the MC method with the MC method based on $128000$ samples,
	we have approximated the $L^2(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the MC method with the MC method based on $128000$ samples, and 
	we have approximated the $L^\infty(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the MC method based on $128000$ random samples
		(cf., e.g., Beck et al.\ \cite[Lemma 3.5]{BeckJafaari21} and Beck et al.\ \cite[Section 3.3]{Beck2019}).
In our approximations of the above mentioned approximation errors we have approximately computed for all required sample points 
	$p \in \mathfrak{P}$
the value $u(p)$ of the unknown exact solution by means of an MC approximation with $52428800$ MC samples.


		


\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|,
     table head=
     \hline $\mathfrak{M}$ & \thead{Number \\of \\trainable \\parameters}&\thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} 
     & \thead{Training \\ time \\ in seconds} 
     \\\hline,
    late after line=\\\hline,
    filter expr={
		test{\ifnumless{\thecsvinputline}{11}}
    }
    ]
{Table_6.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	num-params =\nparams,
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	train-time=\train,
	eval-time=\eval  ,
	dtype = \dtype
}
{ \nsamp& \nparams &  \lll & \llll & \linf &\train
}
}
\caption{\label{table:BSworstPut_LRV}Numerical simulations for the LRV strategy in case of the Black-Scholes model for European worst-of basket put options on three underlyings described in \cref{simul:BSworstPut} (15-dimensional approximation problem) based on reference solutions computed with the MC method}
\end{table}




\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|,
     table head=
     \hline $\mathfrak{M}$ & \thead{Number \\of \\trainable \\parameters}&\thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} 
     & \thead{Training \\ time \\ in seconds} 
     \\\hline,
    late after line=\\\hline,
    filter expr={
   		test{\ifnumless{\thecsvinputline}{11}}
    }
    ]
{Table_7.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	num-params =\nparams,
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	train-time=\train,
	eval-time=\eval  ,
	dtype = \dtype
}
{ \nsamp& \nparams &  \lll & \llll & \linf &\train
}
}
\caption{\label{table:BSworstPut_LRV_self}Numerical simulations for the LRV strategy in case of the Black-Scholes model for European worst-of basket put options on three underlyings described in \cref{simul:BSworstPut} (15-dimensional approximation problem) based on reference solutions computed with the LRV method itself}
\end{table}




\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|,
     table head=
     \hline 
     \thead{Number \\ of \\MC \\samples}&  
     \thead{Number of \\ scalar random \\ variables \\per evaluation} 
     &\thead{$L^1$-approx. \\ error} 
     &\thead{$L^2$-approx. \\error}
     &\thead{ $L^\infty$-approx. \\ error} 
     \\\hline,
    late after line=\\\hline]{Table_8.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	num-inputs = \ninputs,
	l1-error=\lll, 
	l1-error-ratio = \lratio,
	l2-error-ratio = \llratio,
	l-inf-error-ratio = \linfratio,
	l-2-error=\llll, 
	l-inf-error=\linf, 
	time=\train, 
	dtype = \dtype
}
{
	\nsamp
	& \ninputs 
	& \lll 
	& \llll 
	& \linf 
}}
\caption{\label{table:BSworstPut_MC}Numerical simulations for the standard MC method in case of the Black-Scholes model for European worst-of basket put options on three underlyings described in \cref{simul:BSworstPut} (15-dimensional approximation problem)}
\end{table}


\subsection[Parametric Black-Scholes equations for multi-asset average put options with barriers]{Parametric Black-Scholes partial differential equations for multi-asset average put options with knock-in barriers}
\label{simul:BSaverageBarrier}


\newcommand{16384}{16384}
\newcommand{52428800}{52428800}
\newcommand{1024}{1024}
\newcommand{40000}{40000}


\newcommand{\mathfrak{B}}{\mathfrak{B}}
\newcommand{\mathbf{h}}{\mathbf{h}}
\newcommand{\mathbf{m}}{\mathbf{m}}

\newcommand{\mathscr{T}}{\mathscr{T}}
\newcommand{\mathscr{U}}{\mathscr{U}}
\newcommand{\mathscr{L}}{\mathscr{L}}
\newcommand{\mathscr{P}}{\mathscr{P}}




In this section we apply the LRV strategy to a more complicated and practically relevant Black-Scholes option pricing problem. 
Specifically, we approximate the fair price of a European average basket put option on three underlyings with a knock-in barrier in the Black-Scholes model (cf.\ Black \& Scholes \cite{BlackScholes73} and Merton \cite{Merton73}).
We start by introducing the Black-Scholes model for this pricing problem in the context of \cref{algo:general}.

Assume \cref{algo:general},
let 
	$d = 3$, 
	$e_1 = (1, 0, 0)$, 
	$e_2 = (0, 1, 0)$, 
	$e_3 = (0, 0, 1)$,
let $Q \colon \mathbb{R}^d \to \mathbb{R}^{d \times d}$ satisfy for all 
	$\rho = (\rho_1, \rho_2, \rho_3) \in \mathbb{R}^d$
that
\begin{equation}
\begin{split} 
	Q(\rho) 
= 
	\begin{pmatrix}
		1 &\rho_1 &\rho_2 \\
		\rho_1 &1 &\rho_3 \\
		\rho_2 &\rho_3 &1
	\end{pmatrix},
\end{split}
\end{equation}
let 
$
	R
\subseteq 
	\mathbb{R}^d
$
satisfy
$
	R 
= 
	\bigl\{
		\rho = (\rho_1, \rho_2, \rho_3) \in [-\frac{95}{100},\frac{95}{100}]^d 
		\colon  
		1 - |\rho_2|^2 - \frac{(\rho_3 - \rho_1\rho_2)^2}{(1-|\rho_1|^2)^{1/2}} \geq 0
	\bigr\}
$,
assume
	$\mathfrak{p} = 16$
and
\begin{equation}
\begin{split} 
	\mathfrak{P} 
= 
	[90, 110]^3 \times
	[\tfrac{1}{2}, 1] \times
	[-\tfrac{1}{20}, \tfrac{1}{20}] \times
	[0, \tfrac{1}{10}]^3 \times
	[\tfrac{1}{100}, \tfrac{1}{2}]^3 \times 
	R \times 
	[90, 110] \times
	[70, 80],
\end{split}
\end{equation}
let $L \colon R \to \mathbb{R}^{d \times d}$ satisfy for all 
	$\rho = (\rho_1, \rho_2, \rho_3) \in R$
that
\begin{equation}
\begin{split} 
	L(\rho)
=
	\begin{pmatrix}
		1 & 0 & 0 \\
		\rho_1& (1 - |\rho_1|^2)^{{1}/{2}} &0 \\
		\rho_2& \frac{\rho_3 - \rho_1\rho_2}{(1 - |\rho_1|^2)^{{1}/{2}}} & \big(1 - |\rho_2|^2 - \frac{(\rho_3 - \rho_1\rho_2)^2}{1 - |\rho_1|^2}\big)^{{1}/{2}}
	\end{pmatrix},
\end{split}
\end{equation}
let $\mathcal{W} \colon \Omega \to C([0, \infty), \mathbb{R}^{ d }) $ be a standard Brownian motion with continuous sample paths,
for every
	$\xi = (\xi_1, \xi_2, \xi_3) \in \mathbb{R}^d$,
	$r \in \mathbb{R}$,
	$\delta = (\delta_1, \delta_2, \delta_3) \in \mathbb{R}^d$,
	$\sigma = (\sigma_1, \sigma_2, \sigma_3) \in \mathbb{R}^d$,
	$\rho \in \mathbb{R}^d$
let 
$ 
	X^{\xi, r, \delta, \sigma, \rho} 
= 
	(
	X^{r, \xi, \delta, \sigma, \rho, 1}, 
	X^{r, \xi, \delta, \sigma, \rho, 2}, 
	X^{r, \xi, \delta, \sigma, \rho, 3}
	) 
	\colon 
	C([0, \infty), \mathbb{R}^{ d })
	\to 
	C([0, \infty), \mathbb{R}^{ d })
$ 
satisfy for all
	$w \in C([0, \infty), \mathbb{R}^{ d })$,
	$t \in [0,\infty)$,
	$i \in \{1, 2, 3\}$
that
\begin{equation}
\label{BSaverageBarrier:eq1}
\begin{split} 
	X^{\xi, r, \delta, \sigma, \rho, i}_t(w)
=
	\xi_i 
	\exp\Bigl(
		\bigl[r - \delta_i -  \tfrac{|\sigma_i|^2}{2} \bigr]t 
		+
		\sigma_i 
		\langle 
			L(\rho) w_t
			,
			e_i
		\rangle
	\Bigr),
\end{split}
\end{equation}
let 
	$\mathbf{l} \colon \mathbb{R}^d \to \mathbb{R}$ and 
	$\mathbf{m} \colon \mathbb{R}^d \to \mathbb{R}$ 
satisfy for all 
	$x = (x_1, x_2, x_3) \in \mathbb{R}^d$
that
$\mathbf{l}(x) = \min\{ x_1, x_2, x_3 \}$ and
$\mathbf{m}(x) = \tfrac{ x_1 + x_2 + x_3 }{3}$,
for every 
	$T \in [0,\infty)$,
	$B \in \mathbb{R}$
let $\mathfrak{B}_{T, B} \subseteq C([0, \infty), \mathbb{R}^{ d })$ satisfy
\begin{equation}
	\label{BSaverageBarrier:eq2}
	\mathfrak{B}_{T, B}
=
	\bigl\{
		w \in C([0, \infty), \mathbb{R}^{ d })
		\colon 
		(
			\exists \, t \in [0,T] 
			\colon 
			\mathbf{l}(w_t)
			< B
		)  
	\bigr\},
\end{equation}
let 
$
	\phi \colon \mathfrak{P} \times C([0,\infty), \mathbb{R}^d) \to \mathbb{R}
$
satisfy for all
	$p = (\xi, T, r, \delta, \sigma, \rho, K, B) \in \mathfrak{P}$,
	$w \in C([0,\infty), \mathbb{R}^d)$
that
\begin{equation}
\label{BSaverageBarrier:eq2.2}
\begin{split} 
	\phi(p, w)
=
	\mathbbm{1}_{\mathfrak{B}_{T,B}}(X^{\xi, r, \delta, \sigma, \rho}(w))
	\exp(-r T)
	\max\bigl\{
		K - 
		\mathbf{m}(
			X^{\xi, r, \delta, \sigma, \rho}_T(w)
		)
		,
		0
	\bigr\},
\end{split}
\end{equation}
and let $u \colon \mathfrak{P} \to \mathbb{R}$ satisfy for all
	$p \in \mathfrak{P}$
that
\begin{equation}
\label{BSaverageBarrier:eq3}
\begin{split} 
	u(p)
=
	\Exp{
		\phi(p, \mathcal{W})
	}.
\end{split}
\end{equation}
In the economic interpretation of the Black-Scholes model for every 
	$p = (\xi, T, r, \delta, \sigma, \rho, K, B) \in \mathfrak{P}$
the number $u(p) \in \mathbb{R}$ corresponds to the fair price of a average basekt put option on three underlying assets with 
	initial prices $\xi$,
	time of maturity $T$,
	risk free rate $r$,
	dividend yields of the respective underlying assets $\delta$,
	volatilities of the respective underlying assets $\sigma$,
	covariance matrix of the Brownian motions $Q(\rho)$,
	strike price $K$, and
	knock-in barrier $B$.

We now specify the mathematical objects in the LRV strategy appearing in \cref{algo:general} 
to approximately calculate the target function $u \colon \mathfrak{P} \to \mathbb{R}$ in \cref{BSaverageBarrier:eq3}.
Specifically, in addition to the assumptions above, 
let 
	$N = 10$,
	$\mathscr{M} = 1024$,
	$\mathfrak{M} \in \mathbb{N}$,
assume
	$\mathfrak{d} =  \mathfrak{M} N d$,
	$\mathbf{d} = \mathscr{M} N d$,
	$k = 1$,
	${\bf M} = 512$,
for every 
	$p = (\xi, T, r, \delta, \sigma, \rho, K, B) \in \mathfrak{P}$,
	$w = (w_1, \ldots, w_N) \in \mathbb{R}^{N d}$
let 
$
	\mathcal{X}^{p, w} 
= 
	(\mathcal{X}^{p, w, 1}, \mathcal{X}^{p, w, 2}, \mathcal{X}^{p, w, 3} ) 
	\colon 
	\{0, 1, \ldots, N\} \to \mathbb{R}^d
$
satisfy for all 
	$n \in \{1, 2, \ldots, N \}$,
	$i \in \{1, 2, 3\}$
that
$
	\mathcal{X}^{p, w}_0 = \xi	
$
and
\begin{equation}
\label{BSaverageBarrier:eq4}
\begin{split} 
	\mathcal{X}^{p, w, i}_n
&= 
	\mathcal{X}^{p, w, i}_{n-1}
	\exp\Bigl(
		\tfrac{T}{N}
		\bigl[r - \delta_i -  \tfrac{|\sigma_i|^2}{2}\bigr]
		+
		\big[\tfrac{T}{N}\big]^{1/2}
		\sigma_i 
		\langle 
			L(\rho) w_n
			,
			e_i
		\rangle
	\Bigr),
\end{split}
\end{equation}
for every 
	$p = (\xi, T, r, \delta, \sigma_1, \sigma_2, \sigma_3, \rho, K, B) \in \mathfrak{P}$ 
let 
	$ \mathscr{T}_p = (\mathscr{T}_{p,1}, \mathscr{T}_{p,2}, \mathscr{T}_{p,3}) \colon (0,\infty)^d \times (0,\infty)^d \to \mathbb{R}^d$,
	$ \mathscr{U}_p \colon (0,\infty)^d \times (0,\infty)^d \to [0,1]$,
	$ \mathscr{L}_p \colon (0,\infty)^d \times (0,\infty)^d \to [0,1]$, and
	$ \mathscr{P}_p \colon ((0,\infty)^d)^{N+1} \to [0,1]$	
satisfy for all	
	$x = (x_1, x_2, x_3)$, 
	$y  = (y_1, y_2, y_3) \in (0,\infty)^d$,
	$\mathbf{x}_0, \mathbf{x}_1, \ldots, \mathbf{x}_N \in (0,\infty)^d$,
	$i \in \{1, 2, 3\}$
that
\begin{equation}
\label{BSaverageBarrier:eq5.1}
\begin{split} 
	\mathscr{T}_{p,i}(x, y) 
= 
	\begin{cases}
		1 &\colon \min\{ x, y \} < B \\
		\exp\!\left(
			-\frac{
				2 \ln(x_i/B)\ln(y_i/B)
			}{
				(\sigma_i)^2 T/N
			}
		\right)
		&\colon \min\{ x, y \} \geq B,
	\end{cases}
\end{split}
\end{equation}
\begin{equation}
\label{BSaverageBarrier:eq5.2}
\begin{split} 
	\mathscr{U}_{p}(x, y) 
= 
	1 - \max\bigl( \cup_{j \in \{1, 2, 3\}} \{\mathscr{T}_{p,j}(x, y) \} \bigr), \qquad
	\mathscr{L}_{p}(x, y) 
= 
	\max\bigl\{ 1 - \textstyle \sum_{j = 1}^3 \mathscr{T}_{p,j}(x, y), 0  \bigr\},
\end{split}
\end{equation}
\begin{equation}
\label{BSaverageBarrier:eq5.3}
\begin{split} \textstyle
\text{and}\qquad
	\mathscr{P}_{p}(\mathbf{x}_0, \mathbf{x}_1, \ldots, \mathbf{x}_N) 
= 
	\frac{1}{2}
	\left(
		2 
		- 
		\left[ 
			\prod\limits_{n = 1}^N
				\mathscr{U}_{p}(\mathbf{x}_{n-1}, \mathbf{x}_{n}) 
		\right]
		-
		\left[ 
			\prod\limits_{n = 1}^N
				\mathscr{L}_{p}(\mathbf{x}_{n-1}, \mathbf{x}_{n}) 
		\right]
	\right),
\end{split}
\end{equation}
let 
$
	\Phi \colon
	\mathfrak{P} \times \mathbb{R}^{N d} \to
	\mathbb{R}
$
satisfy for all 
	$ p = (\xi, T, r, \delta, \sigma, \rho, K, B) \in \mathfrak{P} $,
	$ w \in \mathbb{R}^{N d}$
that
\begin{equation}
\label{BSaverageBarrier:eq6}
\begin{split} 
	\Phi(p, w)
=
	\mathscr{P}_p(\mathcal{X}^{p, w}_{0}, \mathcal{X}^{p, w}_1, \ldots, \mathcal{X}^{p, w}_N)
	\exp(-r T)
	\max\bigl\{
		K - 
		\mathbf{m}(
			\mathcal{X}^{p, w}_N
		)
		,
		0
	\bigr\},
\end{split}
\end{equation}
assume for all 
	$ p = (\xi, T, r, \delta, \sigma, \rho, K, B) \in \mathfrak{P} $,
	$ w = (w_1,\ldots, w_{\mathscr{M} N d}) \in \mathbb{R}^{\mathscr{M} N d}$,
	$ \theta = (\theta_1,\ldots, \theta_{\mathfrak{M} N d}) \in \mathbb{R}^{\mathfrak{M} N d}$
that
\begin{equation}
\label{BSaverageBarrier:eq7.0}
\begin{split} 
	\Xi(p, w)
=
	\frac{ 1 }{ \mathscr{M} }
	\left[ 
		{\textstyle \sum\limits_{ \mathfrak{m} = 1 }^{ \mathscr{M} } }
			\Phi\big( 
				p, 
				( 
				w_{ ( \mathfrak{m} - 1 ) Nd + k }
				)_{k \in \{1, 2, \ldots, Nd \}}         
			\big)
	\right]
\end{split}
\end{equation}
\begin{equation}
\label{BSaverageBarrier:eq7}
\begin{split} 
\qquad\text{and}\qquad
	\Psi(p, \theta)
=
	\frac{ 1 }{ \mathfrak{M} }
	\left[ 
		{\textstyle \sum\limits_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }}
			\Phi\big( 
				p, 
				( 
				\theta_{ ( \mathfrak{m} - 1 ) Nd + k }
				)_{k \in \{1, 2, \ldots, Nd \}}         
			\big)
	\right],
\end{split}
\end{equation}
assume for all
	$x, y \in \mathbb{R}$
that
	$\mathfrak{C}(x, y) = | x-y |^2$,
let
	$\mathbf{a} = {\frac{9}{10}}$,
	$\mathbf{b}  = {\frac{999}{1000}}$,
	$ \varepsilon \in (0,\infty)$,
let
	$(\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty)$
satisfy for all 
	$m \in \{1, 2, \ldots, 40000 \}$
that
$
	\gamma_m 
=  
	{\mathbbm{1}_{(0,20000]}(m)}{10^{-3}}
	+
	{\mathbbm{1}_{(20000,30000]}(m)}{10^{-4}}
	+
	{\mathbbm{1}_{(30000,40000]}(m)}{10^{-5}}
$,
assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})})$, 	$g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})})$,
	$\dots$, 
	$g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\label{BSaverageBarrier:eq7.2}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\gamma_m
	\left[
		\frac{
			\sum_{k = 1}^m \mathbf{a}^{m-k}(1-\mathbf{a}) g_k^{(i)}
		}{
			1-\mathbf{a}^m
		}
	\right]
	\left[
		\varepsilon +\bigg[
			\frac{
				\sum_{k = 1}^m \mathbf{b}^{m-k}(1-\mathbf{b}) |g_k^{(i)}|^2
			}{
				1- \mathbf{b}^m
			}
		\bigg]^{\nicefrac{1}{2}} 
	\right]^{-1},
\end{split}
\end{equation}	
assume that  
$ \mathfrak{W}$ is a standard normal random vector,
assume that  
$  
  P_{ 
    1, 1 
  }
$ is $\mathcal{U}_{\mathfrak{P}}$-distributed,
and
assume that 
$ 
  W^{ 1, 1 }
$
is a standard normal random vector.

Let us add some comments regarding the setup introduced above.
The functions 
	$ \mathscr{P}_p \colon ((0,\infty)^d)^{N+1} \to [0,1]$,
	$p \in \mathfrak{P}$,
in \cref{BSaverageBarrier:eq4} are employed to estimate crossing probabilities of Brownian bridges as proposed in Shevchenko \cite{Shevchenko2003} (see also, e.g., Gobet \cite{Gobet2009}).
Specifically, note that Gobet \cite[Displays (12), (13)]{Gobet2009} suggests that for all
	$p = (\xi, T, r, \delta, \sigma, \rho, K, B) \in \mathfrak{P}$
we have $\P$-a.s.\ that
\begin{equation}
\label{T_B_D}
\begin{split} 
	&\mathscr{P}_p(
		X^{\xi, r, \delta, \sigma, \rho}_0(\mathcal{W}), 
		X^{\xi, r, \delta, \sigma, \rho}_{T/N}(\mathcal{W}), 
		\ldots, 
		X^{\xi, r, \delta, \sigma, \rho}_T(\mathcal{W})
	) \\
&\approx
	\P\bigl(
			\exists \, t \in [0, T] 
			\colon 
			\mathbf{l}(
				X^{\xi, r, \delta, \sigma, \rho}_t(\mathcal{W})
			)	
			<
			B
		\,\big|\,
			(
				\mathcal{W}_0, 
				\mathcal{W}_{T/N},
				\ldots, 
				\mathcal{W}_T 
			)
	\bigr).
\end{split}
\end{equation}
Combining this with 
	\cref{BSaverageBarrier:eq1},
	\cref{BSaverageBarrier:eq4},
	\cref{BSaverageBarrier:eq6}, 
	the tower property for conditional expectations, and
	the fact that Brownian motions have independent increments
suggests for all
	$p = (\xi, T, r, \delta, \sigma, \rho, K, B) \in \mathfrak{P}$
that
\begin{equation}
\label{T_B_D}
\begin{split} 
	&\Exp{\Phi(p, \mathfrak{W})}
=
	\Exp{
		\Phi\big( 
			p, \sqrt{{N}/{T}}\,
			(\mathcal{W}_{T / N} - \mathcal{W}_0, \mathcal{W}_{ 2 T / N } - \mathcal{W}_{ T / N }, \dots, \mathcal{W}_T -\mathcal{W}_{ (N-1) T / N } )
		\big)
	} \\
&=
	\mathbb{E} \Bigg[
		\mathscr{P}_p(
			X^{\xi, r, \delta, \sigma, \rho}_0(\mathcal{W}), 
			X^{\xi, r, \delta, \sigma, \rho}_{T/N}(\mathcal{W}), 
			\ldots, 
			X^{\xi, r, \delta, \sigma, \rho}_T(\mathcal{W})
		)\\
&\quad
		\exp(-r T)
		\max\bigl\{
			K - 
			\mathbf{m}(
				X^{\xi, r, \delta, \sigma, \rho}_{T}(\mathcal{W})
			)
			,
			0
		\bigr\}
	\Bigg]\\
&\approx
	\mathbb{E} \Bigg[
		\P\bigl(
				\exists \, t \in [0, T] 
				\colon 
				\mathbf{l}(
					X^{\xi, r, \delta, \sigma, \rho}_t(\mathcal{W}) 
				)	
				<
				B
			\,\big|\,
				(
					\mathcal{W}_0, 
					\mathcal{W}_{T/N},
					\ldots, 
					\mathcal{W}_T 
				)
		\bigr)\\
&\quad
		\exp(-r T)
		\max\bigl\{
			K - 
			\mathbf{m}(
				X^{\xi, r, \delta, \sigma, \rho}_{T}(\mathcal{W})
			)
			,
			0
		\bigr\}
	\Bigg]\\
&=
	\mathbb{E} \Bigg[
		\P\bigl(
			X^{\xi, r, \delta, \sigma, \rho}(\mathcal{W}) \in \mathfrak{B}_{T, B}
			\,\big|\,
				\mathcal{W}
		\bigr)
		\exp(-r T)
		\max\bigl\{
			K - 
			\mathbf{m}(
				X^{\xi, r, \delta, \sigma, \rho}_{T}(\mathcal{W})
			)
			,
			0
		\bigr\}
	\Bigg]\\
&=
	\mathbb{E} \Bigg[
		\mathbbm{1}_{\mathfrak{B}_{T, B}}(\mathcal{W})
		\exp(-r T)
		\max\bigl\{
			K - 
			\mathbf{m}(
				X^{\xi, r, \delta, \sigma, \rho}_{T}(\mathcal{W})
			)
			,
			0
		\bigr\}
	\Bigg]
=
	\Exp{
		\phi(p, \mathcal{W})
	}
=
	u(p).
\end{split}
\end{equation}
The proposal algorithm in \cref{BSaverageBarrier:eq3} thus corresponds to the MC method based on approximated MC samples.
Furthermore, observe that \cref{BSaverageBarrier:eq7.2} describes the Adam optimizer in the setup of \cref{algo:general} (cf.\ Kingma \& Ba \cite{Kingma2014} and \cref{sect:Adam}).



In \cref{table:BSaverageBarrier_LRV} we approximately present for
	$\mathfrak{M} \in \{2^5, 2^6, \ldots, 2^{13}\}$
one random realization of the $L^1(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BSaverageBarrier_LRV:eq1}
\begin{split} 
\textstyle
	\int_\mathfrak{P} 
		|u(p) - \Psi(p, \Theta_{40000})|
	\d p
\end{split}
\end{equation}
(\nth{3} column in \cref{table:BSaverageBarrier_LRV}),
one random realization of the $L^2(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BSaverageBarrier_LRV:eq2}
\begin{split} 
\textstyle
	\bigl[
		\int_\mathfrak{P} 
			|u(p) - \Psi(p, \Theta_{40000})|^2
		\d p
	\bigr]^{\nicefrac{1}{2}}
\end{split}
\end{equation}
(\nth{4} column in \cref{table:BSaverageBarrier_LRV}),
one random realization of the $L^\infty(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{BSaverageBarrier_LRV:eq3}
\begin{split} 
\textstyle
	\sup_{p \in \mathfrak{P}}
		|u(p) - \Psi(p, \Theta_{40000})|,
\end{split}
\end{equation}
(\nth{5} column in \cref{table:BSaverageBarrier_LRV}), and
the time to compute $\Theta_{40000}$
(\nth{6} column in \cref{table:BSaverageBarrier_LRV}).
For every 	
	$\mathfrak{M} \in \{2^5, 2^6, \ldots, 2^{13}\}$
we approximated the integrals in \eqref{BSaverageBarrier_LRV:eq1} and \eqref{BSaverageBarrier_LRV:eq2} with the MC method based on $16384$ samples
and we approximated the supremum in \eqref{BSaverageBarrier_LRV:eq3} based on $16384$ random samples 
(cf., e.g., Beck et al.\ \cite[Lemma 3.5]{BeckJafaari21} and Beck et al.\ \cite[Section 3.3]{Beck2019}).
In our approximations of \eqref{BSaverageBarrier_LRV:eq1}, \eqref{BSaverageBarrier_LRV:eq2}, and \eqref{BSaverageBarrier_LRV:eq3} we have approximately computed for all required sample points 
	$p \in \mathfrak{P}$
the value $u(p)$ of the unknown exact solution 
by means of an MC approximation with $52428800$ MC samples in \cref{table:BSaverageBarrier_LRV} and
by means of the LRV strategy itself (i.e., $\Psi(p, \Theta_{40000})$) with $\mathfrak{M} = 2^{13} = 8192$ in \cref{table:BSaverageBarrier_LRV_self}. 


Besides the LRV strategy we also employed the standard MC method to approximate the function $u \colon \mathfrak{P} \to \mathbb{R}$ in \eqref{BSaverageBarrier:eq3}. 
In \cref{table:BSaverageBarrier_MC} we present the corresponding numerical simulation results.
In \cref{table:BSaverageBarrier_MC}
	we have approximated the $L^1(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
		of the MC method with the MC method based on $16384$ samples,
	we have approximated the $L^2(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the MC method with the MC method based on $16384$ samples, and 
	we have approximated the $L^\infty(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the MC method based on $16384$ random samples
		(cf., e.g., Beck et al.\ \cite[Lemma 3.5]{BeckJafaari21} and Beck et al.\ \cite[Section 3.3]{Beck2019}).
In our approximations of the above mentioned approximation errors we have approximately computed for all required sample points 
	$p \in \mathfrak{P}$
the value $u(p)$ of the unknown exact solution by means of an MC approximation with $ 52428800$ MC samples.

		

\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|,
     table head=
     \hline $\mathfrak{M}$ & \thead{Number \\of \\trainable \\parameters}&\thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} 
     & \thead{Training \\ time \\ in seconds} 
     \\\hline,
    late after line=\\\hline]{Table_9.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	num-params =\nparams,
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	train-time=\train,
	eval-time=\eval  ,
	dtype = \dtype
}
{ \nsamp& \nparams &  \lll & \llll & \linf &\train }
}
\caption{\label{table:BSaverageBarrier_LRV}Numerical simulations for the LRV strategy in case of the Black-Scholes model for European average put options with knock-in barriers described in \cref{simul:BSaverageBarrier} (16-dimensional approximation problem) based on reference solutions computed with the MC method}
\end{table}


\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|,
     table head=
     \hline $\mathfrak{M}$ & \thead{Number \\of \\trainable \\parameters}&\thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} 
     & \thead{Training \\ time \\ in seconds} 
     \\\hline,
    late after line=\\\hline]{Table_10.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	num-params =\nparams,
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	train-time=\train,
	eval-time=\eval  ,
	dtype = \dtype
}
{ \nsamp& \nparams &  \lll & \llll & \linf &\train
}
}
\caption{\label{table:BSaverageBarrier_LRV_self}Numerical simulations for the LRV strategy in case of the Black-Scholes model for European average put options with knock-in barriers described in \cref{simul:BSaverageBarrier} (16-dimensional approximation problem) based on reference solutions computed with LRV method itself}
\end{table}



\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|c|c|,
     table head=
     \hline 
     \thead{Number \\ of \\MC \\samples}
     &\thead{Number of \\ scalar random \\ variables \\per evaluation} 
     &\thead{$L^1$-approx. \\ error} 
     & \thead{$L^2$-approx. \\error}
     &\thead{ $L^\infty$-approx. \\ error} 
     \\\hline,
    late after line=\\\hline]{Table_11.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	num-inputs = \ninputs,
	l1-error=\lll, 
	l1-error-ratio = \lratio,
	l2-error-ratio = \llratio,
	l-inf-error-ratio = \linfratio,
	l-2-error=\llll, 
	l-inf-error=\linf, 
	time=\train, 
	dtype = \dtype
}
{
	\nsamp
	& \ninputs 
	& \lll 
	& \llll 
	& \linf 
}}
\caption{\label{table:BSaverageBarrier_MC}Numerical simulations for the standard MC method in case of the Black-Scholes model for European average put options with knock-in barriers described in \cref{simul:BSaverageBarrier} (16-dimensional approximation problem)}
\end{table}








\subsection{Parametric stochastic Lorentz equations}
\label{simul:Lorentz}

\newcommand{10000}{10000}
\newcommand{6400} \newcommand{\numRefAlgSamplesLorentz}{{256}}{6400} \newcommand{\numRefAlgSamplesLorentz}{{256}}
\newcommand{{8\,196\,000}}{{8\,196\,000}}



In this section we apply the LRV strategy to the parametric stochastic Lorentz equation (cf., e.g, Schmallfuss \cite{Schmallfuss1997} and Hutzenthaler \& Jentzen \cite[Section 4.4]{Hutzenthaler12}).
A brief summary of the numerical results of this subsection can be found in \cref{table:Lorentz_summary} below.
We first introduce the parametric stochastic Lorentz equation in the context of \cref{algo:general}.

\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|,
	separator=semicolon,
     table head=
     \hline Approximation method& \thead{Number of \\ trainable \\ parameters} & \thead{Number\\ of \\ MC \\ samples} & \thead{$L^2$-approx. \\ error}& \thead{$L^\infty$-approx. \\ error}& \thead{Training \\ time in \\ seconds} 
     \\\hline,
    late after line=\\\hline]{Table_12.csv}
{
	method = \method, 
	num-samples = \nsamp, 
	num-params = \nparams,
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	train-time=\train, 
	eval-time=\eval
}
{ \method & \nparams &\nsamp& \llll & \linf &\train 
}
}
\caption{\label{table:Lorentz_summary}Parametric stochastic Lorentz equation}
\end{table}



Assume \cref{algo:general},
let $d = 3$, 
let $g \colon \mathbb{R}^d \to \mathbb{R}$ satisfy for all
	$x \in \mathbb{R}^d$
that
	$g(x) = \|x\|^2$,
for every let 
	$\alpha = (\alpha_1, \alpha_2, \alpha_3) \in \mathbb{R}^d$
let $\mu_\alpha \colon \mathbb{R}^d \to \mathbb{R}^d$ satisfy for all
	$x = (x_1, x_2, x_3) \in \mathbb{R}^d$
that
\begin{equation}
\begin{split} 
	\mu_\alpha(x) 
=
	(\alpha_1(x_2 - x_1), \alpha_2 x_1 - x_2 - x_1x_3, x_1x_2 - \alpha_3 x_3),
\end{split}
\end{equation}
and
for every 
	$\alpha = (\alpha_1, \alpha_2, \alpha_3), \beta = (\beta_1, \beta_2, \beta_3) \in \mathbb{R}^d$
let $u_{\alpha, \beta}=(u_{\alpha, \beta}(t,x))_{(t,x)\in [0,\infty)\times\mathbb{R}^d}\in C^{1,2}([0,\infty)\times\mathbb{R}^d,\mathbb{R})$
be an at most polynomially growing function which satisfies
for all 
	$t\in [0,\infty)$, 
	$x\in\mathbb{R}^d$ 
that 
\begin{multline}
\label{Lorentz:eq1}
(\tfrac{\partial u_{\alpha, \beta}}{\partial t})(t,x)
= 
\alpha_1(x_2-x_1) (\tfrac{\partial u_{\alpha, \beta}}{\partial x_1})(t,x)
+ 
(\alpha_2 x_1 - x_2 - x_1 x_3)(\tfrac{\partial u_{\alpha, \beta}}{\partial x_2})(t,x) \\
+ 
(x_1 x_2 - \alpha_3 x_3)(\tfrac{\partial u_{\alpha, \beta}}{\partial x_3})(t,x)
+ 
\sum_{i = 1}^d \tfrac{(\beta_i)^2}{2}(\tfrac{\partial u_{\alpha, \beta}}{\partial (x_i)^2})(t,x)
\end{multline}
and $u_{\alpha, \beta}(0,x) = g(x)$.
 
We now specify the mathematical objects in the LRV strategy appearing in \cref{algo:general} 
to approximately calculate $(u_{\alpha,\beta})_{(\alpha,\beta) \in \mathbb{R}^d \times \mathbb{R}^d}$.
Specifically, in addition to the assumptions above,
let $\odot \colon \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d$ satisfy for all
	$x = (x_1, x_2, x_3)$,
	$y = (y_1, y_2, y_3)$
that
$	
	x \odot y = (x_1 y_1, x_2 y_2, x_3 y_3)
$,
let 
	$N = {25}$,
	$\mathscr{M} = \numRefAlgSamplesLorentz$,
	$\mathfrak{M} \in \mathbb{N}$,
	$a \in \{0,1\}$,
assume
	$\mathfrak{p} = 10$,
	$\mathfrak{d} =  \mathfrak{M} N d$,
	$\mathbf{d} = \mathscr{M} N d$,
	$k = 1$,
	${\bf M} = {256}$, and
	\begin{equation}
	\begin{split} 
		\mathfrak{P} = [0.01,1] \times ([9, 11] \times [13, 15] \times [{1, 2}]) \times [0.05, 0.25]^3 \linebreak \times ([0.5, 2.5] \times [8, 10] \times [10, 12]),
	\end{split}
	\end{equation}
for every 
	$p = (T, \alpha, \beta, x) \in \mathfrak{P}$,
	$w = (w_1, \ldots, w_N) \in \mathbb{R}^{N d}$
let 
	$\mathcal{X}^{p, w} \colon \{0, 1, \ldots, N\} \to \mathbb{R}^d  $
satisfy for all 
	$n \in \{1, 2, \ldots, N \}$
that
$
	\mathcal{X}^{p, w}_0 = x	
$
and
\begin{equation}
\label{Lorentz:eq2}
\begin{split} 
	\mathcal{X}^{p, w}_n
&= 
	\mathcal{X}^{p, w}_{n-1} 
	+
	\sqrt{T/N}  (\beta \odot  w_n )
	+
	\tfrac{T}{2N} 
	\big(
		\mu_\alpha(\mathcal{X}^{p, w}_{n-1}) 
		+
		\mu_\alpha\big(
			\mathcal{X}^{p, w}_{n-1}
			+
			\tfrac{T}{N}
			\mu_\alpha(\mathcal{X}^{p, w}_{n-1}) 
			+
			\sqrt{T/N}  (\beta \odot w_n)
		\big) 
	\big),
\end{split}
\end{equation}
let 
$
	\Phi_k \colon
	\mathfrak{P} \times \mathbb{R}^{N d} \to
	\mathbb{R}
$,
	$k \in \{0, 1\}$,
satisfy for all 
	$ p \in \mathfrak{P} $,
	$ w \in \mathbb{R}^{N d}$
that
\begin{equation}
\begin{split} 
	\Phi_0( p, w)
=
	g(\mathcal{X}^{p, w}_N)
\qquad\text{and}\qquad
	\Phi_1( p, w)
=
	\tfrac{1}{2}\big[g(\mathcal{X}^{p, w}_N) + g(\mathcal{X}^{p, -w}_N)\big],
\end{split}
\end{equation}
assume for all 
	$ p  \in \mathfrak{P} $,
	$ w = (w_1, \ldots, w_{\mathscr{M} N d}) \in \mathbb{R}^{\mathscr{M} N d}$,
	$ \theta = (\theta_1,\ldots, \theta_{\mathfrak{M} N d}) \in \mathbb{R}^{\mathfrak{M} N d}$
that
\begin{equation}
\label{Lorentz:eq3.0}
\begin{split} 
	\Xi(p, w)
=
	\frac{ 1 }{ \mathscr{M} }
	\left[ 
		{\textstyle \sum\limits_{ \mathfrak{m} = 1 }^{ \mathscr{M} }}
			\Phi_a\big( 
				p, 
				( 
				w_{ ( \mathfrak{m} - 1 ) Nd + k }
				)_{k \in \{1, 2, \ldots, Nd \}}         
			\big)
	\right]
\end{split}
\end{equation}
\begin{equation}
\label{Lorentz:eq3}
\begin{split} 
\text{and}\qquad
	\Psi(p, \theta)
=
	\frac{ 1 }{ \mathfrak{M} }
	\left[ 
		{\textstyle \sum\limits_{ \mathfrak{m} = 1 }^{ \mathfrak{M} } }
			\Phi_a\big( 
				p, 
				( 
				\theta_{ ( \mathfrak{m} - 1 ) Nd + k }
				)_{k \in \{1, 2, \ldots, Nd \}}         
			\big)
	\right],
\end{split}
\end{equation}
assume for all
	$x, y \in \mathbb{R}$
that
	$\mathfrak{C}(x, y) = | x-y |^2$,
let
	$\mathbf{a} = {\frac{9}{10}}$,
	$\mathbf{b}  = {\frac{999}{1000}}$,
	$ \varepsilon \in (0,\infty)$,
let
	$(\gamma_m)_{m \in \mathbb{N}} \subseteq (0,\infty)$
satisfy for all 
	$m \in \{1, 2, \ldots, 10000 \}$
that
$
	\gamma_m 
=  
	{\mathbbm{1}_{(0,5000]}(m)}{10^{-3}}
	+
	{\mathbbm{1}_{(5000,8000]}(m)}{10^{-4}}
	+
	{\mathbbm{1}_{(8000,10000]}(m)}{10^{-5}}
$,
assume for all
	$m \in \mathbb{N}$,
	$i \in \{1, 2,\ldots,\mathfrak{d} \}$,
	$g_1 = (g_1^{(1)},\ldots,g_1^{(\mathfrak{d})})$, 	$g_2 = (g_2^{(1)},\ldots,g_2^{(\mathfrak{d})})$,
	$\dots$, 
	$g_m = (g_m^{(1)},\ldots,g_m^{(\mathfrak{d})})
	\in 
		\mathbb{R}^{\mathfrak{d}}
	$
that
\begin{equation}
\label{Lorentz:eq4}
\begin{split} 
	\psi_m^{(i)}(g_1, g_2, \ldots, g_m)
=
	\gamma_m
	\left[
		\frac{
			\sum_{k = 1}^m \mathbf{a}^{m-k}(1-\mathbf{a}) g_k^{(i)}
		}{
			1-\mathbf{a}^m
		}
	\right]
	\left[
		\varepsilon +\bigg[
			\frac{
				\sum_{k = 1}^m \mathbf{b}^{m-k}(1-\mathbf{b}) |g_k^{(i)}|^2
			}{
				1- \mathbf{b}^m
			}
		\bigg]^{\nicefrac{1}{2}} 
	\right]^{-1},
\end{split}
\end{equation}	
assume that  
$ \mathfrak{W}$ is a standard normal random vector,
assume that  
$  
  P_{ 
    1, 1 
  }
$ is $\mathcal{U}_{\mathfrak{P}}$-distributed,
and
assume that 
$ 
  W^{ 1, 1 }
$
is a standard normal random vector.

Let us add some comments regarding the setup introduced above.
Observe that \eqref{Lorentz:eq2} corresponds to a Heun discretization (cf., e.g., Kloeden \cite[(1.4) in Section 15]{KloedenPlaten13}) of the SDE associated to the Kolmogorov PDE in \eqref{Lorentz:eq1}.
In the case $a = 0$ the proposal algorithm in \eqref{Lorentz:eq3} thus corresponds to the MC-Heun method and 
in the case $a = 1$ the proposal algorithm in \eqref{Lorentz:eq3} thus corresponds to the antithetic MC-Heun method.
In particular, we remark that for all
	$p = (T, \alpha, \beta, x) \in \mathfrak{P}$
we have that
\begin{equation}
\begin{split} 
	  \frac{ 1 }{ \mathfrak{M} }
	  \left[
	    \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
	    \Phi_0\big( 
	    p, 
	    W^{0,\mathfrak{m}}
	    \big)
	  \right]  
\approx
	u_{\alpha, \beta}(T,x) 
\approx
	  \frac{ 1 }{ \mathfrak{M} }
	  \left[
	    \sum_{ \mathfrak{m} = 1 }^{ \mathfrak{M} }
	    \Phi_1\big( 
	    p, 
	    W^{0,\mathfrak{m}}
	    \big)
	  \right].
\end{split}
\end{equation}
Furthermore, observe that \cref{Lorentz:eq4} describes the Adam optimizer in the setup of \cref{algo:general} (cf.\ Kingma \& Ba \cite{Kingma2014} and \cref{sect:Adam}).

In \cref{table:Lorentz_MC} we approximately present for
	$\mathfrak{M} \in \{2^5, 2^6, \ldots, 2^{13}\}$,
	$a \in \{0,1\}$
one random realization of the $L^1(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{Lorentz_lrv:eq1}
\begin{split} 
\textstyle
	\int_\mathfrak{P} 
		|u_{\alpha, \beta}(T, x) - \Psi((T, \alpha, \beta, x), \Theta_{10000})|
	\d \, (T, \alpha, \beta, x)
\end{split}
\end{equation}
(\nth{3} column in \cref{table:Lorentz_MC}),
one random realization of the $L^2(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{Lorentz_lrv:eq2}
\begin{split} 
\textstyle
	\bigl[
		\int_\mathfrak{P} 
			|u_{\alpha, \beta}(T, x) - \Psi((T, \alpha, \beta, x), \Theta_{10000})|^2
		\d \, (T, \alpha, \beta, x)
	\bigr]^{\nicefrac{1}{2}}
\end{split}
\end{equation}
(\nth{4} column in \cref{table:Lorentz_MC}),
one random realization of the $L^\infty(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation error
\begin{equation}
\label{Lorentz_lrv:eq3}
\begin{split} 
	\sup_{(T, \alpha, \beta, x) \in \mathfrak{P}}
		|u_{\alpha, \beta}(T, x) - \Psi((T, \alpha, \beta, x), \Theta_{10000})|
\end{split}
\end{equation}
(\nth{5} column in \cref{table:Lorentz_MC}), and
the time to compute $\Theta_{10000}$
(\nth{6} column in \cref{table:Lorentz_MC}).
For every 	
	$\mathfrak{M} \in \{2^5, 2^6, \ldots, 2^{13}\}$,
	$a \in \{0,1\}$
we approximated the integrals in \eqref{Lorentz_lrv:eq1} and \eqref{Lorentz_lrv:eq2} with the MC method based on $6400} \newcommand{\numRefAlgSamplesLorentz}{{256}$ samples
and we approximated the supremum in \eqref{Lorentz_lrv:eq3} based on $6400} \newcommand{\numRefAlgSamplesLorentz}{{256}$ random samples
(cf., e.g., Beck et al.\ \cite[Lemma 3.5]{BeckJafaari21} and Beck et al.\ \cite[Section 3.3]{Beck2019}).
In our approximations of \eqref{Lorentz_lrv:eq1}, \eqref{Lorentz_lrv:eq2}, and \eqref{Lorentz_lrv:eq3} we have approximately computed for all required sample points 
	$(T, \alpha, \beta, x) \in \mathfrak{P}$
the value $u_{\alpha, \beta}(T, x)$ of the unknown exact solution by means of an antithetic MC approximation with ${8\,196\,000}$ MC samples.

To compare the LRV strategy with existing approximation techniques from the literature, we also employ several other methods to approximate the function $\mathfrak{P} \ni (T, \alpha, \beta, x) \mapsto u_{\alpha, \beta}(T, x) \in \mathbb{R}$ in \eqref{Lorentz:eq1}.
Specifically, 
	in \cref{table:Lorentz_ANNs} we present numerical simulations for the deep learning method induced by Becker et al.\ \cite{BeckJafaari21} 
	(with Adam $160000$ training steps, 
	batch size $256$, and 
	learning rate schedule 
	$
		\mathbb{N} \ni j 
		\mapsto 
			{\mathbbm{1}_{(0,50000]}(m)}{10^{-1}}
			+
			{\mathbbm{1}_{(50000,100000]}(m)}{10^{-2}}
			+
			{\mathbbm{1}_{(100000,120000]}(m)}{10^{-3}}
			+
			{\mathbbm{1}_{(120000,140000]}(m)}{10^{-4}}
			+
			{\mathbbm{1}_{(140000,160000]}(m)}{10^{-5}}
	$),
	and
	in \cref{table:Lorentz_MC} we present numerical simulations for the standard and the antithetic MC method.
In Tables~\ref{table:Lorentz_ANNs} and \ref{table:Lorentz_MC} 
	we have approximated the $L^1(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the respective approximation methods with the MC method based on $6400} \newcommand{\numRefAlgSamplesLorentz}{{256}$ samples,
	we have approximated the $L^2(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the respective approximation methods with the MC method based on $6400} \newcommand{\numRefAlgSamplesLorentz}{{256}$ samples, and 
	we have approximated the $L^\infty(\lambda_{\mathfrak{P}};\mathbb{R})$-approximation errors
		of the respective approximation methods based on $6400} \newcommand{\numRefAlgSamplesLorentz}{{256}$ random samples
		(cf., e.g., Beck et al.\ \cite[Lemma 3.5]{BeckJafaari21} and Beck et al.\ \cite[Section 3.3]{Beck2019}).
In our approximations of the above mentioned approximation errors we have approximately computed for all required sample points 
	$(T, \alpha, \beta, x) \in \mathfrak{P}$
the value $u_{\alpha, \beta}(T, x)$ of the unknown exact solution by means of an antithetic MC approximation with ${8\,196\,000}$ MC samples












\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|,
     table head=
     \hline $\mathfrak{M}$ &  $a$ &\thead{Number\\of\\trainable\\parameters}& \thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} & \thead{Training \\ time \\ in \\seconds} 
     \\\hline,
    late after line=\\\hline]{Table_13.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	num-params = \nparam,
	antithetic=\anti, 
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	rel-l1-error = \rel, 
	rel-l-2-error = \rell,
	rel-l-inf-error = \relinf,
	train-time = \train,
	eval-time=\eval, 
	dtype = \dtype
}
{\nsamp& \anti &\nparam & \lll & \llll & \linf &\train
}
}
\caption{\label{table:Lorentz_LRV}Numerical simulations for the LRV strategy in case of the parametric stochastic Lorentz equation described in \cref{simul:Lorentz} (10-dimensional approximation problem)}
\end{table}







\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|c|,
     table head=
	\hline \thead{Number \\ of \\ hidden \\ layers}& \thead{Number  \\of neurons \\ on each \\ hidden \\ layer}& \thead{Number \\ of \\ trainable \\ parameters} & \thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} & \thead{Training time \\ in seconds} 
	\\\hline,
	 late after line=\\\hline]{Table_14.csv}
{
	layers = \lyrs,
	inner-dim = \inner,
	num-weights = \weights,
	mc-samples = \mcsamp, 
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	rel-l1-error = \rel, 
	rel-l-2-error = \rell, 
	rel-l-inf-error = \relinf,
	train-time=\train, 
	time=\eval  ,
	dtype = \dtype
}
{\pgfmathparse{int(\lyrs )}\pgfmathresult & \inner &\weights & \lll & \llll & \linf &\train 
}}
\caption{\label{table:Lorentz_ANNs}Numerical simulations for the deep learning method induced by Becker et al.\ \cite{BeckJafaari21} in case of the parametric stochastic Lorentz equation described in \cref{simul:Lorentz} (10-dimensional approximation problem)}
\end{table}


\begin{table}[H] \tiny
\resizebox{\textwidth}{!}{
\csvreader[tabular=|c|c|c|c|c|c|c|c|c|,
     table head=
     \hline \thead{Number \\ of \\MC \\samples}&\thead{Number of \\ scalar random \\ variables \\per evaluation} &\thead{MC Method \\ 0: standard \\ 1: antithetic} & \thead{$L^1$-approx. \\ error}& \thead{$L^2$-approx. \\error}&\thead{ $L^\infty$-approx. \\ error} 
     \\\hline,
    late after line=\\\hline]{Table_15.csv}
{
	num-samples=\nsamp, 
	mc-samples = \mcsamp, 
	num-inputs = \ninputs,
	antithetic=\anti, 
	l1-error=\lll, 
	l-2-error=\llll, 
	l-inf-error=\linf, 
	rel-l1-error = \rel, 
	rel-l-2-error = \rell, 
	rel-l-inf-error = \relinf,
	time=\train, 
	dtype = \dtype
}
{\nsamp& \ninputs &\anti& \lll & \llll & \linf 
}
}
\caption{\label{table:Lorentz_MC}Numerical simulations for the standard and the antithetic MC method in case of the parametric stochastic Lorentz equation described in \cref{simul:Lorentz} (10-dimensional approximation problem)}
\end{table}











\section{{\sc Python} source code}
\label{sect:source}
In this section we present one representative {\sc Python} code for the LRV strategy. 
Specifically {\sc Python }~code~\ref{code:BS1lrv} was employed to produce the results in \cref{table:BS_LRV} in \cref{simul:BS1}.
Note that the source code can be downloaded as a separate file  {\it lrv\_bs.py} on arXiv by clicking, first, on \emph{Other formats} and, thereafter, on \emph{Download source}.


\lstinputlisting[language=Python,
caption={\it Source code associated to the LRV simulations in \cref{simul:BS1}},
label={code:BS1lrv}]
{lrv_bs.py}




\section*{Acknowledgments}
Nor Jaafari is gratefully acknowledged for his useful assistance regarding some numerical simulations.
The second author acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2044-390685587, Mathematics Muenster: Dynamics-Geometry-Structure.

\bibliographystyle{acm}

