Stochastic methods for solving high-dimensional partial differential equations

TL;DR

This paper introduces novel stochastic algorithms that leverage probabilistic representations and sparse interpolation techniques to efficiently solve high-dimensional PDEs, addressing computational challenges in complex systems.

Contribution

It develops new algorithms combining control variates and adaptive sparse interpolation for high-dimensional PDEs, enhancing computational efficiency and accuracy.

Findings

Algorithms effectively solve high-dimensional PDEs

Numerical examples demonstrate improved accuracy

Methods outperform traditional approaches in complexity

Abstract

We propose algorithms for solving high-dimensional Partial Differential Equations (PDEs) that combine a probabilistic interpretation of PDEs, through Feynman-Kac representation, with sparse interpolation. Monte-Carlo methods and time-integration schemes are used to estimate pointwise evaluations of the solution of a PDE. We use a sequential control variates algorithm, where control variates are constructed based on successive approximations of the solution of the PDE. Two different algorithms are proposed, combining in different ways the sequential control variates algorithm and adaptive sparse interpolation. Numerical examples will illustrate the behavior of these algorithms.

Figures

Tables

Table 1. Table 1: Algorithm 3.2 computed on the exact solution of ( TC1 ): evolution of # ​ Λ n # subscript Λ 𝑛 \#\Lambda_{n} , error criterion ε n superscript 𝜀 𝑛 \varepsilon^{n} and interpolation errors in norms L 2 superscript 𝐿 2 L^{2} and L ∞ superscript 𝐿 L^{\infty} at each step n 𝑛 n .

${\Lambda }_n$	$\#{\Lambda }_n$	${\epsilon }^n$	${\Vert u-u_n\Vert }_2$	${\Vert u-u_n\Vert }_{\infty }$
	1	6.183372e-01	1.261601e+00	4.213566e+00
	10	2.792486e-02	1.204421e-01	3.602629e-01
	20	2.178450e-05	9.394419e-04	3.393999e-03
${\Lambda }^6$	26	9.632815e-07	4.270457e-06	1.585129e-05
	30	9.699704e-08	2.447475e-06	8.316435e-06
${\Lambda }^8$	33	4.114730e-09	2.189518e-08	9.880306e-08
	40	1.936050e-10	6.135776e-10	1.739848e-09
${\Lambda }^{10}$	41	1.008412e-11	9.535433e-11	4.781375e-10
	50	1.900248e-14	1.004230e-13	4.223288e-13
	55	7.453467e-15	2.905404e-14	1.254552e-13

Table 2. Table 2: (TC1) Comparison of the algorithmic complexity to reach the precision 3 ⋅ 10 − 5 ⋅ 3 superscript 10 5 3\cdot 10^{-5} , with ( Δ ​ t , M ) = ( 10 − 4 , 1000 ) Δ 𝑡 𝑀 superscript 10 4 1000 (\Delta t,M)=(10^{-4},1000) .

	Adaptive Algorithm 2.3	Algorithm 4.1	Full-grid Algorithm 2.3
Th. Complexity	$M{\left(\Delta t\right)}^{-1}\left({\sum }_k\#{\Lambda }_k\right)$	$M{\left(\Delta t\right)}^{-1}\left({\sum }_n\#{\Lambda }_n{N}_n\right)$	$M{\left(\Delta t\right)}^{-1}\#{\Lambda }_{max}N$
Est. Complexity	$4\cdot {10}^9$ operations	$16\cdot {10}^9$ operations	${10}^{12}N$ operations

Table 3. Table 3: (TC2) Comparison of Algorithm 3.2 (first four columns) and Algorithm 4.1 (last four columns).

$\#{\Lambda }_n$	${\epsilon }^n$	${\Vert u-u_n\Vert }_{\infty }$	${\Vert u-u_n\Vert }_2$	$\#{\Lambda }_n$	${\tilde{\epsilon }}^n$	${\Vert u-{\tilde{u}}_{{\Lambda }_n}\Vert }_{\infty }$	${\Vert u-{\tilde{u}}_{{\Lambda }_n}\Vert }_2$
1	4.0523e-01	3.0151e+00	1.2094e+00	1	3.9118e-01	8.3958e-01	6.9168e-01
17	1.6243e-01	1.8876e+00	5.9579e-01	17	1.6259e-01	5.2498e-01	3.4420e-01
36	5.4494e-02	7.0219e-01	2.0016e-01	36	5.4699e-02	1.9209e-01	1.2594e-01
46	1.2767e-02	1.6715e-01	4.9736e-02	46	1.2806e-02	4.6904e-02	2.8524e-02
53	9.6987e-04	2.9343e-02	4.8820e-03	53	1.0350e-03	7.8754e-03	2.8960e-03
60	7.6753e-04	1.5475e-02	4.1979e-03	61	7.0354e-04	3.0365e-03	1.7610e-03
71	3.2532e-04	8.4575e-03	2.1450e-03	71	3.1998e-04	2.3486e-03	1.2395e-03
77	1.7434e-16	3.9968e-15	1.5784e-15	77	7.3621e-16	6.2172e-15	1.2874e-15

Table 4. Table 4: (TC3) Comparison of Algorithm 3.2 (first four columns) and Algorithm 4.1 (last four columns).

$\#{\Lambda }_n$	${\epsilon }^n$	${\Vert u-u_n\Vert }_{\infty }$	${\Vert u-u_n\Vert }_2$	$\#{\Lambda }_n$	${\tilde{\epsilon }}^n$	${\Vert u-{\tilde{u}}_{{\Lambda }_n}\Vert }_{\infty }$	${\Vert u-{\tilde{u}}_{{\Lambda }_n}\Vert }_2$
1	7.0155e-01	3.9361e+00	1.2194e+00	1	5.5582e-01	7.2832e-01	7.0771e-01
6	1.4749e-01	2.2705e+00	5.4886e-01	6	7.4253e-02	2.7579e-01	5.1539e-01
11	2.1902e-02	2.8669e-01	1.0829e-01	11	1.4929e-02	4.4614e-02	4.1973e-02
15	7.6086e-03	1.6425e-01	4.7394e-02	15	1.2916e-02	1.5567e-02	2.5650e-02
20	2.2275e-04	2.7715e-03	7.2230e-04	20	3.4446e-04	5.6927e-04	5.3597e-04
24	1.4581e-05	1.5564e-04	7.5314e-05	24	1.6036e-05	2.5952e-05	3.0835e-05
30	1.8263e-06	8.0838e-06	2.1924e-06	30	9.0141e-07	2.8808e-06	1.9451e-06
35	3.9219e-09	8.9815e-08	2.4651e-08	35	8.1962e-09	2.1927e-08	1.5127e-08
40	1.7933e-10	2.0152e-09	6.9097e-10	40	1.6755e-10	2.8455e-10	2.6952e-10
45	5.0775e-12	2.4783e-10	4.1600e-11	45	1.4627e-11	3.3188e-11	1.7911e-11
49	1.7722e-14	4.6274e-13	8.5980e-14	49	1.7938e-14	8.6362e-14	5.0992e-14
54	3.9609e-15	2.2681e-13	3.1952e-14	54	3.2195e-15	4.8142e-14	2.6617e-14
56	4.5746e-16	8.4376e-15	3.0438e-15	56	8.2539e-16	8.4376e-15	6.3039e-15

Equations

$$\begin{array}{rl}\mathcal{A}(u)& =g\text{in}\mathcal{D},\\ u& =f\text{on}\partial \mathcal{D},\end{array}$$

$$\mathcal{L}(u)(x)=\frac{1}{2}\sum _{i,j=1}^d(\sigma (x)\sigma (x{)}^T{)}_{ij}{\partial }_{x_i{x}_j}^{2}u(x)+\sum _{i=1}^d{b}_i(x){\partial }_{x_i}u(x)$$

$$d{X}_t^x=b(X_t^x)dt+\sigma (X_t^x)d{W}_t,X_0^x=x\in \overline{\mathcal{D}},$$

$$\mathbb{P}({\tau }^x=0)=1,x\in \partial \mathcal{D},$$

$$\sum _{i,j=1}^d{\left(\sigma (x)\sigma (x{)}^T\right)}_{ij}{\xi }_i{\xi }_j\ge c\sum _{i=1}^d{\xi }_i^2,\xi \in {\mathbb{R}}^d,x\in \overline{\mathcal{D}}.$$

$$\begin{array}{rl}u(x)& =\mathbb{E}\left[F(u,X^x)\right]\end{array}$$

$$F(u,{{X}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{x}}})=u({{X}}^{{{x}}}_{\tau^{{{x}}}})\exp\left(-\int_{0}^{\tau^{{{x}}}}k({{X}}^{{{x}}}_{t})dt\right)+\int_{0}^{\tau^{{{x}}}}{\mathcal{A}}(u)({{X}}^{{{x}}}_{t})\exp\left(-\int_{0}^{t}k({{X}}^{{{x}}}_{s})ds\right)dt,$$

$$\begin{array}{rl}{X}_{n+1}^{x,\Delta t}& =X_n^{x,\Delta t}+\Delta tb(X_n^{x,\Delta t})+\sigma (X_n^{x,\Delta t})\Delta W_n,\\ X_0^{x,\Delta t}& =x.\end{array}$$

$$\begin{split}u_{\Delta t,M}({{x}})&:=\dfrac{1}{M}\sum_{m=1}^{M}F\left(u,{{X}}^{{{x}},\Delta t}(\omega_{m})\right)\\ &=\dfrac{1}{M}\sum_{m=1}^{M}\bigg{[}f({{X}}^{{{x}},\Delta t}_{\tau^{{{x}},\Delta t}}(\omega_{m}))\exp\left(-\int_{0}^{\tau^{{{x}},\Delta t}}k({{X}}^{{{x}},\Delta t}_{t}(\omega_{m}))dt\right)\\ &+\int_{0}^{\tau^{{{x}},\Delta t}}g({{X}}^{{{x}},\Delta t}_{t}(\omega_{m}))\exp\left(-\int_{0}^{t}k({{{X}}}^{{{x}},\Delta t}_{s}(\omega_{m}))ds\right)dt\bigg{]}\end{split}$$

$${\tau }^{x,\Delta t}=\inf \left\{t>0:X_t^{x,\Delta t}\notin \mathcal{D}\right\}=\min \left\{t_n>0:X_{t_n}^{x,\Delta t}\notin \mathcal{D}\right\}.$$

$$\begin{array}{rl}u(x)-u_{\Delta t,M}(x)=& \stackrel{{\epsilon }_{\Delta t}}{\overbrace{u(x)-\mathbb{E}\left[F\left(u,X^{x,\Delta t}\right)\right]}}\\ & +\underset{{\epsilon }_{MC}}{\underbrace{\mathbb{E}\left[F\left(u,X^{x,\Delta t}\right)\right]-\frac{1}{M}\sum _{m=1}^{M}F\left(u,X^{x,\Delta t}({\omega }_m)\right)}},\end{array}$$

$${\mathcal{I}}_{\Lambda }(w)(x_{\nu })=w(x_{\nu }),\nu \in \Lambda .$$

$$\displaystyle\begin{array}[]{rcll}&{\mathcal{A}}(e^{k})({{x}})=g({{x}})-{\mathcal{A}}(\tilde{u}^{k-1})({{x}}),&{{x}}\in{\mathcal{D}},\\ &e^{k}({{x}})=f({{x}})-\tilde{u}^{k-1}({{x}}),&{{x}}\in\partial{\mathcal{D}}.\end{array}$$

$$e^{\Delta t}(h,x)=\mathbb{E}[F(h,X^{\Delta t,x})]-\mathbb{E}[F\left(h,X^x\right)].$$

$$\underset{\nu \in \Lambda }{\sup }\left|\mathbb{E}\left[{\tilde{u}}^{n+1}(x_{\nu })-u(x_{\nu })\right]\right|⩽C(\Delta t,\Lambda )\underset{\nu \in \Lambda }{\sup }\left|\mathbb{E}\left[{\tilde{u}}^n(x_{\nu })-u(x_{\nu })\right]\right|+C_1(\Delta t,\Lambda )$$

$$\underset{n\to \infty }{\mathrm{limsup}}\underset{\nu \in \Lambda }{\sup }\left|\mathbb{E}\left[{\tilde{u}}^n(x_{\nu })-u(x_{\nu })\right]\right|⩽\frac{C_1(\Delta t,\Lambda )}{1-C(\Delta t,\Lambda )}.$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{L}}_{\Lambda}}}=\sup_{v\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{C}^{0}(\overline{\mathcal{D}})}}}\dfrac{\|{\mathcal{I}}_{\Lambda}(v)\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}}{\|v\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}},$$

$$\|v-{\mathcal{I}}_{\Lambda}(v)\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}\leq(1+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{L}}_{\Lambda}}})\inf_{w\in V_{\Lambda}}\|v-w\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}.$$

$$\limsup_{n\to\infty}~{}\|\mathbb{E}\left[\tilde{u}^{n}-u\right]\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}\leqslant\dfrac{C_{1}(\Delta t,\Lambda)}{1-C(\Delta t,\Lambda)}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{L}}_{\Lambda}}}+\|u-{\mathcal{I}}_{\Lambda}(u)\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}.$$

$$\|{\mathbb{E}}\left[\tilde{u}^{n}-u\right]\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}\leqslant\|{\mathbb{E}}\left[\tilde{u}^{n}-{\mathcal{I}}_{\Lambda}(u)\right]\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}+\|{\mathcal{I}}_{\Lambda}(u)-u\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}.$$

$$\Vert w{\Vert }_{\infty }=\underset{\nu \in \Lambda }{\sup }|w(x_{\nu })|=\underset{\nu \in \Lambda }{\sup }\left|\mathbb{E}\left[{\tilde{u}}^n(x_{\nu })-u(x_{\nu })\right]\right|.$$

$$\displaystyle\|{\mathbb{E}}\left[\tilde{u}^{n}-{\mathcal{I}}_{\Lambda}(u)\right]\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}=\|{\mathcal{I}}_{\Lambda}(w)\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}\leq{\mathcal{L}}_{\Lambda}\|w\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}.$$

$$\Vert v{\Vert }_2\le |\mathcal{D}{|}^{1/2}\Vert v{\Vert }_{\infty },$$

$${\phi }_{\nu }(x)=\prod _{i=1}^d{\phi }_{{\nu }_i}^{(i)}(x_i).$$

$$\forall \nu \in \Lambda ,\mu \le \nu \Rightarrow \mu \in \Lambda .$$

$${\Gamma }_{\Lambda }=\{z_{\nu }=(z_{{\nu }_1}^{(1)},\dots ,z_{{\nu }_d}^{(d)}):\nu \in \Lambda \}\subset \overline{\mathcal{D}}.$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{L}}_{\Lambda}}}\leq\left(\#\Lambda\right)^{s+1}.$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{L}}_{\Lambda_{n}}}}\inf_{w\in{\mathcal{P}}_{\Lambda_{n}}}\|u-w\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}\longrightarrow 0\text{ as }n\to\infty$$

$$\inf_{w\in{\mathcal{P}}_{\Lambda_{n}}}\|u-w\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}=O((\#\Lambda_{n})^{-r})$$

$${\mathcal{M}}_{{\Lambda }_n}=\{\nu \in {\mathbb{N}}^d\setminus {\Lambda }_n:\forall j\text{s.t.}{\nu }_j>0,\nu -e_j\in {\Lambda }_n\},$$

Full text

1 Introduction

We consider the solution of an elliptic partial differential equation

[TABLE]

where $u:\overline{{\mathcal{D}}}\rightarrow{\mathbb{R}}$ is a real-valued function, and ${\mathcal{D}}$ is an open bounded domain in ${\mathbb{R}}^{d}$. ${\mathcal{A}}$ is an elliptic linear differential operator and $f:\partial{\mathcal{D}}\to{\mathbb{R}}$, $g:\overline{{\mathcal{D}}}\rightarrow{\mathbb{R}}$ are respectively the boundary condition and the source term of the PDE.

We are interested in approximating the solution of (1) up to a given precision. For high dimensional PDEs ($d\gg 1$), this requires suitable approximation formats such as sparse tensors [5, 26] or low-rank tensors [24, 18, 19, 1, 23]. Also, this requires algorithms that provide approximations in a given approximation format. Approximations are typically provided by Galerkin projections using variational formulations of PDEs. Another path consists in using a probabilistic representation of the solution $u$ through Feynman-Kac formula, and Monte-Carlo methods to provide estimations of pointwise evaluations of $u$ (see e.g., [17]). This allows to compute approximations in a given approximation format through classical interpolation or regression [3, 27, 4]. In [14, 15], the authors consider interpolations on fixed polynomial spaces and propose a sequential control variates method for improving the performance of Monte-Carlo estimation. In this paper, we propose algorithms that combine this variance reduction method with adaptive sparse interpolation [6, 7].

The outline is as follows. In section 2, we recall the theoretical and numerical aspects associated to probabilistic tools for estimating the solution of (1). We also present the sequential control variates algorithm introduced in [14, 15]. In section 3 we introduce sparse polynomial interpolation methods and present a classical adaptive algorithm. In section 4, we present two algorithms combining the sequential control variates algorithm from section 2 and adaptive sparse polynomial interpolation. Finally, numerical results are presented in section 4.

2 Probabilistic tools for solving PDEs

We consider the problem (1) with a linear partial differential operator defined by ${\mathcal{A}}(u)=-{\mathcal{L}}(u)+ku$, where $k$ is a real valued function defined on $\overline{\mathcal{D}}$, and where

[TABLE]

is the infinitesimal generator associated to the $d$-dimensional diffusion process ${{X}}^{{{x}}}$ solution of the stochastic differential equation

[TABLE]

where $W$ is a $d$-dimensional Brownian motion and $b:=(b_{1},\ldots,b_{d})^{T}:{\mathbb{R}}^{d}\to{\mathbb{R}}^{d}$ and $\sigma:{\mathbb{R}}^{d}\to{\mathbb{R}}^{d\times d}$ stand for the drift and the diffusion respectively.

2.1 Pointwise evaluations of the solution

The following theorem recalls the Feynman-Kac formula (see [11, Theorem 2.4] or [12, Theorem 2.4] and the references therein) that provides a probabilistic representation of $u({{x}})$, the solution of (1) evaluated at ${{x}}\in\overline{{\mathcal{D}}}$.

Theorem 2.1** (Feynman-Kac formula).**

Assume that

1. $(H{1})$

${\mathcal{D}}$* is an open connected bounded domain of ${\mathbb{R}}^{d}$, regular in the sense that, if $\tau^{{{x}}}=\inf\left\{s>0~{}:~{}{{X}}^{{{x}}}_{s}\notin{\mathcal{D}}\right\}$ is the first exit time of ${\mathcal{D}}$ for the process ${{X}}^{{{x}}}$, we have*

[TABLE] 2. $(H{2})$

$b,\sigma$* are Lipschitz functions,* 3. $(H{3})$

$f$* is continuous on $\partial{\mathcal{D}}$, $g$ and $k\geq 0$ are Hölder-continuous functions on $\overline{{\mathcal{D}}}$,* 4. $(H{4})$

(uniform ellipticity assumption) there exists $c>0$ such that

[TABLE]

Then, there exists a unique solution of (1) in ${\mathcal{C}}\left(\overline{{\mathcal{D}}}\right)\cap{\mathcal{C}}^{2}\left({\mathcal{D}}\right)$, which satisfies for all ${{x}}\in\overline{{\mathcal{D}}}$

[TABLE]

where

[TABLE]

with $u({{X}}^{{{x}}}_{\tau^{{{x}}}})=f({{X}}^{{{x}}}_{\tau^{{{x}}}})$ and ${\mathcal{A}}(u)({{X}}^{{{x}}}_{t})=g({{X}}^{{{x}}}_{t})$.

Note that $F(u,X^{x})$ in (4) only depends on the values of $u$ on $\partial D$ and ${\mathcal{A}}(u)$ on $D$, which are the given data $f$ and $g$ respectively. A Monte-Carlo method can then be used to estimate $u(x)$ using (4), which relies on the simulation of independent samples of an approximation of the stochastic process ${{X}}^{{x}}$. This process is here approximated by an Euler-Maruyama scheme. More precisely, letting $t_{n}=n\Delta t$, $n\in\mathbb{N}$, ${{X}}^{{{x}}}$ is approximated by a piecewise constant process ${{X}}^{{{x}},\Delta t}$, where ${{X}}^{{{x}},\Delta t}_{t}={{X}}^{{{x}},\Delta t}_{n}$ for $t\in[t_{n},t_{n+1}[$ and

[TABLE]

Here $\Delta W_{n}=W_{n+1}-W_{n}$ is an increment of the standard Brownian motion. For details on time-integration schemes, the reader can refer to [20]. Letting $\{{{X}}^{{{x}},\Delta t}(\omega_{m})\}_{m=1}^{M}$ be independent samples of ${{X}}^{{{x}},\Delta t}$, we obtain an estimation $u_{\Delta t,M}({{x}})$ of $u({{x}})$ defined as

[TABLE]

where $\tau^{{{x}},\Delta t}$ is the first exit time of $D$ for the process ${{X}}^{{{x}},\Delta t}(\omega_{m})$, given by

[TABLE]

Remark 2.2**.**

In practice, $f$ has to be defined over $\mathbb{R}^{d}$ and not only on the boundary $\partial\cal D$. Indeed, although $X^{{{x}}}_{\tau^{{{x}}}}\in\partial\cal D$ with probability one, $X^{{{x}},\Delta t}_{\tau^{{{x}},\Delta t}}\in\mathbb{R}^{d}\setminus\overline{{\mathcal{D}}}$ with probability one.

The error can be decomposed in two terms

[TABLE]

where $\varepsilon_{\Delta t}$ is the time integration error and $\varepsilon_{MC}$ is the Monte-Carlo estimation error. Before discussing the contribution of each of both terms to the error, let us introduce the following additional assumption, which ensures that ${\mathcal{D}}$ does not have singular points***Note that together with $(H{4})$ , assumption $(H{5})$ implies $(H{1})$ (see [15, §4.1] for details), so that the set of hypotheses $(H{1})$-$(H{5})$ could be reduced to $(H{2})$-$(H{5})$..

$(H{5})$ Each point of $\partial{\mathcal{D}}$ satisfies the exterior cone condition which means that, for all ${{x}}\in\partial{\mathcal{D}}$, there exists a finite right circular cone $K$, with vertex ${{x}}$, such that $\overline{K}\cap\overline{\mathcal{D}}=\{{{x}}\}$.

Under assumptions $(H{1})$-$(H{5})$, it can be proven [15, §4.1] that the time integration error $\varepsilon_{\Delta t}$ converges to zero. It can be improved to $O(\Delta t^{1/2})$ by adding differentiability assumptions on the boundary [16]. The estimation error $\varepsilon_{MC}$ is a random variable with zero mean and standard deviation converging as $O(M^{-1/2})$. The computational complexity for computing a pointwise evaluation of $u_{\Delta t,M}({{x}})$ is in $O\left(M\Delta t^{-1}\right)$ in expectation for $\Delta t$ sufficiently small†††A realization of $X^{x,\Delta t}$ over the time interval $[0,\tau^{x,\Delta t}]$ can be computed in $O\left(\tau^{x,\Delta t}\Delta t^{-1}\right)$. Then, the complexity to evaluate $u_{\Delta t,M}(x)$ is in $O({\mathbb{E}}(\tau^{x,\Delta t})M\Delta t^{-1})$ in expectation. Under $(H{1})$-$(H{5})$, it is stated in the proof of [15, Theorem 4.2] that $\sup_{x}{\mathbb{E}}[\tau^{x,\Delta t}]\leq C$ with $C$ independent of $\Delta t$ for $\Delta t$ sufficiently small., so that the computational complexity for achieving a precision $\epsilon$ (root mean squared error) behaves as $O(\epsilon^{-4})$. This does not allow to obtain a very high accuracy in a reasonable computational time. The convergence with $\Delta t$ can be improved to $O(\Delta t)$ by suitable boundary corrections [16], therefore yielding a convergence in $O(\epsilon^{-3}).$ To further improve the convergence, high-order integration schemes could be considered (see [20] for a survey). Also, variance reduction methods can be used to further improve the convergence, such as antithetic variables, importance sampling, control variates (see [17]). Multilevel Monte-Carlo [13] can be considered as a variance reduction method using several control variates (associated with processes $X^{{{x}},\Delta t_{k}}$ using different time discretizations). Here, we rely on the sequential control variates algorithm proposed in [14] and analyzed in [15]. This algorithm constructs a sequence of approximations of $u$. At each iteration of the algorithm, the current approximation is used as a control variate for the estimation of $u$ through Feynman-Kac formula.

2.2 A sequential control variates algorithm

Here we recall the sequential control variates algorithm introduced in [14] in a general interpolation framework. We let $V_{\Lambda}\subset{\mathcal{C}}^{2}(\overline{\mathcal{D}})$ be an approximation space of finite dimension $\#\Lambda$ and let ${\mathcal{I}}_{\Lambda}:\mathbb{R}^{\mathcal{D}}\to V_{\Lambda}$ be the interpolation operator associated with a unisolvent grid $\Gamma_{\Lambda}=\{{{x}}_{{\nu}}:{{\nu}}\in\Lambda\}$. We let $({{l}}_{{\nu}})_{{{\nu}}\in\Lambda}$ denote the (unique) basis of $V_{\Lambda}$ that satisfies the interpolation property ${{l}}_{{\nu}}({{x}}_{{\mu}})=\delta_{{{\nu}}{{\mu}}}$ for all ${{\nu}},{{\mu}}\in\Lambda$. The interpolation ${\mathcal{I}}_{\Lambda}(w)=\sum_{\nu\in\Lambda}w(x_{\nu})l_{\nu}(x)$ of function $w$ is then the unique function in $V_{\Lambda}$ such that

[TABLE]

The following algorithm provides a sequence of approximations $(\tilde{u}^{k})_{k\geq 1}$ of $u$ in $V_{\Lambda}$, which are defined by $\tilde{u}^{k}=\tilde{u}^{k-1}+\tilde{e}^{k}$, where $\tilde{e}^{k}$ is an approximation of $e^{k}$, solution of

[TABLE]

Note that $e^{k}$ admits a Feyman-Kac representation $e^{k}(x)=\mathbb{E}(F(e^{k},X^{{{x}}}))$, where $F(e^{k},X^{{{x}}})$ depends on the residuals $g-{\mathcal{A}}(\tilde{u}^{k-1})$ on $\cal D$ and $f-\tilde{u}^{k-1}$ on $\partial{\mathcal{D}}$. The approximation $\tilde{e}^{k}$ is then defined as the interpolation ${\mathcal{I}}_{\Lambda}(e^{k}_{\Delta t,M})$ of the Monte-Carlo estimate $e^{k}_{\Delta t,M}(x)$ of $e^{k}_{\Delta t}(x)=\mathbb{E}(F(e^{k},X^{{{x}},\Delta t}))$ (using $M$ samples of $X^{x,\Delta t}$).

Algorithm 2.3**.**

(Sequential control variates algorithm)

1:  Set $\tilde{u}^{0}=0$, $k=1$ and $S=0$.

2:  while $k\leq K~{}\text{and}~{}S<n_{s}$ do

3:     Compute $e^{k}_{\Delta t,M}({{x}}_{{\nu}})$ for ${{x}}_{{\nu}}\in\Gamma_{\Lambda}$.

4:     Compute $\tilde{e}^{k}={\mathcal{I}}_{\Lambda}(e^{k}_{\Delta t,M})=\sum_{{{\nu}}\in\Lambda}e^{k}_{\Delta t,M}({{x}}_{{\nu}}){{l}}_{{\nu}}({{x}})$.

5:     Update $\tilde{u}^{k}=\tilde{u}^{k-1}+\tilde{e}^{k}$.

6:     If $\|\tilde{u}^{k}-\tilde{u}^{k-1}\|_{2}\leq\epsilon_{tol}\|\tilde{u}^{k-1}\|_{2}$ then $S=S+1$ else $S=0$.

7:     Set $k=k+1$.

8:  end while

For practical reasons, Algorithm 2.3 is stopped using an heuristic error criterion based on stagnation. This criterion is satisfied when the desired tolerance $\epsilon_{tol}$ is reached for $n_{s}$ successive iterations (in practice we chose $n_{s}=5$).

Now let us provide some convergence results for Algorithm 2.3. To that goal, we introduce the time integration error at point ${{x}}$ for a function $h$

[TABLE]

Then the following theorem [15, Theorem 3.1] gives a control of the error in expectation.

Theorem 2.4**.**

Assuming $(H{2})$-$(H{5})$, it holds

[TABLE]

*with $\displaystyle C(\Delta t,\Lambda)=\underset{{{\nu}}\in\Lambda}{\sup}\sum_{{{\mu}}\in\Lambda}|e^{\Delta t}({{l}}_{{\mu}},{{x}}_{{\nu}})|$ and $C_{1}(\Delta t,\Lambda)=\sup_{{{\nu}}\in\Lambda}\left|e^{\Delta t}(u-{\mathcal{I}}_{\Lambda}(u),{{x}}_{{\nu}})\right|$.

Moreover if $C(\Delta t,\Lambda)~{}<1$, it holds*

[TABLE]

The condition $C(\Delta t,\Lambda)<1$ implies that in practice $\Delta t$ should be chosen sufficiently small [15, Theorem 4.2]. Under this condition, the error at interpolation points uniformly converges geometrically up to a threshold term depending on time integration errors for interpolation functions $l_{\nu}$ and the interpolation error $u-{\mathcal{I}}_{\Lambda}(u)$.

Theorem 9 provides a convergence result at interpolation points. Below, we provide a corollary to this theorem that provides a convergence result in L^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}({{\mathcal{D}}}). This result involves the Lebesgue constants in L^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}-norm associated to ${\mathcal{I}}_{\Lambda}$, defined by

[TABLE]

and such that for any v\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{C}}^{0}(\overline{\mathcal{D}})}},

[TABLE]

Throughout this article, we adopt the convention that supremum exclude elements with norm [math]. We recall also that the $L^{\infty}$ Lebesgue constant can be expressed as ${\mathcal{L}}_{\Lambda}=\sup_{x\in\overline{\mathcal{D}}}\sum_{\nu\in\Lambda}|{{l}}_{{\nu}}({{x}})|$.

Corollary 2.5** (Convergence in L^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}).**

Assuming $(H{2})$-$(H{5})$, it holds

[TABLE]

Proof.

By triangular inequality, we have

[TABLE]

We can build a continuous function $w$ such that $w(x_{\nu})={\mathbb{E}}\left[\tilde{u}^{n}(x_{\nu})-u(x_{\nu})\right]$ for all $\nu\in\Lambda$, and such that

[TABLE]

We have then

[TABLE]

The result follows from the definition of the function $w$ and Theorem 9.

Remark 2.6**.**

Since for bounded domains ${\mathcal{D}}$, we have

[TABLE]

*for all $v$ in ${\mathcal{C}}^{0}(\overline{{\mathcal{D}}})$, where $|{\mathcal{D}}|$ denotes the Lebesgue measure of ${\mathcal{D}}$, we can deduce the convergence results in $L^{2}$ norm from those in $L^{\infty}$ norm. *

3 Adaptive sparse interpolation

We here present sparse interpolation methods following [6, 7].

3.1 Sparse interpolation

For $1\leq i\leq d$, we let \{\varphi_{k}^{(i)}\}_{k\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathbb{N}}_{0}}}} be a univariate polynomial basis, where $\varphi^{(i)}_{k}(x_{i})$ is a polynomial of degree $k$. For a multi-index \nu=(\nu_{1},\dots,\nu_{d})\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{N}_{0}^{d}}}, we introduce the multivariate polynomial

[TABLE]

For a subset $\Lambda\subset\mathbb{N}^{d}$, we let ${\mathcal{P}}_{\Lambda}=\mathrm{span}\{\varphi_{{\nu}}:\nu\in\Lambda\}$. A subset $\Lambda$ is said to be downward closed if

[TABLE]

If $\Lambda$ is downward closed, then the polynomial space ${\mathcal{P}}_{\Lambda}$ does not depend on the choice of univariate polynomial bases and is such that ${\mathcal{P}}_{\Lambda}=\mathrm{span}\{x^{\nu}:\nu\in\Lambda\}$, with $x^{\nu}=x_{1}^{\nu_{1}}\ldots x_{d}^{\nu_{d}}$.

In the case where ${\mathcal{D}}={\mathcal{D}}_{1}\times\ldots\times{\mathcal{D}}_{d}$, we can choose for \{\varphi^{(i)}_{k}\}_{k\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathbb{N}}_{0}}}} an orthonormal basis in $L^{2}({\mathcal{D}}_{i})$ (i.e. a rescaled and shifted Legendre basis). Then \{\varphi_{\nu}\}_{\nu\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{N}_{0}^{d}}}} is an orthonormal basis of $L^{2}(\cal D).$ To define a set of points $\Gamma_{\Lambda}$ unisolvent for ${\mathcal{P}}_{\Lambda}$, we can proceed as follows. For each dimension $1\leq i\leq d$, we introduce a sequence of points \{z^{(i)}_{k}\}_{k\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathbb{N}}_{0}}}} in $\overline{\mathcal{D}}_{i}$ such that for any $p\geq 0$, $\Gamma^{(i)}_{p}=\{z^{(i)}_{k}\}_{k=0}^{p}$ is unisolvent for ${\mathcal{P}}_{p}=\mathrm{span}\{\varphi^{(i)}_{k}:0\leq k\leq p\}$, therefore defining an interpolation operator ${\mathcal{I}}^{(i)}_{p}$. Then we let

[TABLE]

This construction is interesting for adaptive sparse algorithms since for an increasing sequence of subsets $\Lambda_{n}$, we obtain an increasing sequence of sets $\Gamma_{\Lambda_{n}}$, and the computation of the interpolation on ${\mathcal{P}}_{\Lambda_{n}}$ only requires the evaluation of the function on the new set of points $\Gamma_{\Lambda_{n}}\setminus\Gamma_{\Lambda_{n-1}}.$ Also, with such a construction, we have the following property of the Lebesgue constant of ${\mathcal{I}}_{\Lambda}$ in L^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}-norm. This result is directly taken from [7, Section 3].

Proposition 3.1**.**

If for each dimension $1\leq i\leq d$, the sequence of points \{z_{k}^{(i)}\}_{k\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathbb{N}}_{0}}}} is such that the interpolation operator ${\mathcal{I}}_{p}^{(i)}$ has a Lebesgue constant {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{L}}_{p}}}\leq(p+1)^{s} for some $s>0$, then for any downward closed set $\Lambda$, the Lebesgue constant {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathcal{L}}_{\Lambda}}} satisfies

[TABLE]

Leja points or magic points [21] are examples of sequences of points such that the interpolation operators ${\mathcal{I}}_{p}^{(i)}$ have Lebesgue constants not growing too fast with $p$. For a given $\Lambda$ with $\rho_{i}:=\max_{\nu\in\Lambda}\nu_{i}$, it is possible to construct univariate interpolation grids $\Gamma^{(i)}_{\rho_{i}}$ with better properties (e.g., Chebychev points), therefore resulting in better properties for the associated interpolation operator ${\mathcal{I}}_{\Lambda}$. However for Chebychev points, e.g., $\rho_{i}\leq\rho^{\prime}_{i}$ does not ensure $\Gamma^{(i)}_{\rho_{i}}\subset\Gamma^{(i)}_{\rho^{\prime}_{i}}$. Thus with such univariate grids, an increasing sequence of sets $\Lambda_{n}$ will not be associated with an increasing sequence of sets $\Gamma_{\Lambda_{n}}$, and the evaluations of the function will not be completely recycled in adaptive algorithms. However, for some of the algorithms described in Section 4, this is not an issue as evaluations can not be recycled anyway.

Note that for general domains ${\mathcal{D}}$ which are not the product of intervals, the above constructions of grids $\Gamma_{\Lambda}$ are not viable since it may yield to grids not contained in the domain ${\mathcal{D}}$. For such general domains, magic points obtained through greedy algorithms could be considered.

3.2 Adaptive algorithm for sparse interpolation

An adaptive sparse interpolation algorithm consists in constructing a sequence of approximations $(u_{n})_{n\geq 1}$ associated with an increasing sequence of downward closed subsets $(\Lambda_{n})_{n\geq 1}$. According to (11), we have to construct a sequence such that the best approximation error and the Lebesgue constant are such that

[TABLE]

for obtaining a convergent algorithm. For example, if

[TABLE]

holds‡‡‡see e.g. [9] for conditions on $u$ ensuring such a behavior of the approximation error. for some $r>1$ and if ${\mathcal{L}}_{\Lambda_{n}}=O((\#\Lambda_{n})^{k})$ for $k<r$, then the error \|u-u_{n}\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\infty}}}=O(n^{-r^{\prime}}) tends to zero with an algebraic rate of convergence $r^{\prime}=r-k>0$. Of course, the challenge is to propose a practical algorithm that constructs a good sequence of sets $\Lambda_{m}$.

We now present the adaptive sparse interpolation algorithm with bulk chasing procedure introduced in [6]. Let $\theta$ be a fixed bulk chasing parameter in $(0,1)$ and let ${\mathcal{E}}_{\Lambda}(v)=\|P_{\Lambda}(v)\|^{2}_{2}$, where $P_{\Lambda}$ is the orthogonal projector over ${\mathcal{P}}_{\Lambda}$ for any subset \Lambda\subset{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathbb{N}}_{0}^{d}}}.

Algorithm 3.2**.**

(Adaptive interpolation algorithm)

1:  Set $\Lambda_{1}=\{\boldsymbol{0}_{d}\}$ and $n=1$.

2:  while $n\leq N$ and $\varepsilon^{n-1}>\varepsilon$ do

3:     Compute ${\mathcal{M}}_{\Lambda_{n}}$.

4:     Set $\Lambda_{n}^{\star}=\Lambda_{n}\cup{\mathcal{M}}_{\Lambda_{n}}$ and compute ${\mathcal{I}}_{\Lambda_{n}^{\star}}(u)$.

5:     Select $N_{n}\subset{\mathcal{M}}_{\Lambda_{n}}$ the smallest such that ${\mathcal{E}}_{N_{n}}({\mathcal{I}}_{\Lambda_{n}^{\star}}(u))\geq\theta{\mathcal{E}}_{{\mathcal{M}}_{\Lambda_{n}}}({\mathcal{I}}_{\Lambda_{n}^{\star}}(u))$

6:     Update $\Lambda_{n+1}=\Lambda_{n}\cup N_{n}$.

7:     Compute $u_{n+1}={\mathcal{I}}_{\Lambda_{n+1}}(u)$ (this step is not necessary in practice).

8:     Compute $\varepsilon^{n}$.

9:     Update $n=n+1$.

10:  end while

At iteration $n$, Algorithm 3.2 selects a subset of multi-indices $N_{n}$ in the reduced margin of $\Lambda_{n}$ defined by

[TABLE]

where $({{e}}_{j})_{k}=\delta_{kj}.$ The reduced margin is such that for any subset $S\subset{\mathcal{M}}_{\Lambda_{n}}$, ${\Lambda_{n}}\cup S$ is downward closed. This ensures that the sequence $(\Lambda_{n})_{n\geq 1}$ generated by the algorithm is an increasing sequence of downward closed sets. Finally, Algorithm 3.2 is stopped using a criterion based on

[TABLE]

4 Combining sparse adaptive interpolation with sequential control variates algorithm

We present in this section two ways of combining Algorithm 2.3 and Algorithm 3.2. First we introduce a perturbed version of Algorithm 3.2 and then an adaptive version of Algorithm 2.3. At the end of the section, numerical results will illustrate the behavior of the proposed algorithms.

4.1 Perturbed version of Algorithm 3.2

As we do not have access to exact evaluations of the solution $u$ of (1), Algorithm 3.2 can not be used for interpolating $u$. So we introduce a perturbed version of this algorithm, where the computation of the exact interpolant ${\mathcal{I}}_{\Lambda}(u)$ is replaced by an approximation denoted $\tilde{u}_{\Lambda}$, which can be computed for example with Algorithm 2.3 stopped for a given tolerance $\epsilon_{tol}$ or at step $k$. This brings the following algorithm.

Algorithm 4.1**.**

(Perturbed adaptive sparse interpolation algorithm)

1:  Set $\Lambda_{1}=\{\boldsymbol{0}_{d}\}$ and $n=1$.

2:  while $n\leq N$ and $\tilde{\varepsilon}^{n-1}>\varepsilon$ do

3:     Compute ${\mathcal{M}}_{\Lambda_{n}}$.

4:     Set $\Lambda_{n}^{\star}=\Lambda_{n}\cup{\mathcal{M}}_{\Lambda_{n}}$ and compute $\tilde{u}_{\Lambda_{n}^{\star}}$.

5:     Select $N_{n}$ as the smallest subset of ${\mathcal{M}}_{\Lambda_{n}}$ such that ${\mathcal{E}}_{N_{n}}(\tilde{u}_{\Lambda_{n}^{\star}})\geq\theta{\mathcal{E}}_{{\mathcal{M}}_{\Lambda_{n}}}(\tilde{u}_{\Lambda_{n}^{\star}})$

6:     Update $\Lambda_{n+1}=\Lambda_{n}\cup N_{n}$.

7:     Compute $\tilde{u}_{\Lambda_{n+1}}$.

8:     Compute $\tilde{\varepsilon}^{n}$.

9:     Update $n=n+1$.

10:  end while

4.2 Adaptive version of Algorithm 2.3

As a second algorithm, we consider the sequential control variates algorithm (Algorithm 2.3) where at step 4, an approximation $\tilde{e}^{k}$ of $e^{k}$ is obtained by applying the adaptive interpolation algorithm (Algorithm 4.1) to the function $e^{k}_{\Delta t,M}$, which uses Monte-Carlo estimations $e^{k}_{\Delta t,M}(x_{\nu})$ of $e^{k}(x_{\nu})$ at interpolation points. At each iteration, $\tilde{e}^{k}$ therefore belongs to a different approximation space ${\mathcal{P}}_{\Lambda_{k}}.$ In the numerical section, we will call this algorithm adaptive Algorithm 2.3.

4.3 Numerical results

In this section, we illustrate the behavior of algorithms previously introduced on different test cases. We consider the simple diffusion equation

[TABLE]

were ${\mathcal{D}}=]-1,1[^{d}$. The source terms and boundary conditions will be specified later for each test case.

The stochastic differential equation associated to (15) is the following

[TABLE]

where $(W_{t})_{t\geq 0}$ is a $d$-dimensional Brownian motion.

We use tensorized grids of magic points for the selection of interpolation points evolved in adaptive algorithms.

Small dimensional test case. We consider a first test case (TC1) in dimension $d=5$. Here the source term and the boundary conditions in problem (15) are chosen such that the solution is given by

[TABLE]

We first test the influence of $\Delta t$ and $M$ on the convergence of Algorithm 2.3 when $\Lambda$ is fixed. In that case, $\Lambda$ is selected a priori with Algorithm 3.2 using samples of the exact solution $u$ for (TC1), stopped for $\varepsilon\in\{10^{-6},10^{-8},10^{-10}\}$. In what follows, the notation $\Lambda^{i}$ stands for the set obtained for $\varepsilon=10^{-i},i\in\{6,8,10\}$. We represent on Figure 1 the evolution of the absolute error in $L^{2}$-norm (similar results hold for the $L^{\infty}$-norm) between the approximation and the true solution with respect to step $k$ for $\Lambda=\Lambda^{6}$. As claimed in Corollary 1, we recover the geometric convergence up to a threshold value that depends on $\Delta t$. We also notice faster convergence as $M$ increases and when $\Delta t$ decreases. We fix $M=1000$ in the next simulations.

We study the impact of the choice of $\Lambda^{i}$ on the convergence of Algorithm 2.3. Again we observe on Figure 2 that the convergence rate gets better as $\Delta t$ decreases. Moreover as $\#\Lambda$ increases the threshold value decreases. This is justified by the fact that interpolation error decreases as $\#\Lambda^{i}$ increases (see Table 1). Nevertheless, we observe that it may also deteriorate the convergence rate if it is chosen too large together with $\Delta t$ not sufficiently small. Indeed for the same number of iterations $k=10$ and the same time-step $\Delta t=2.5\cdot 10^{-3}$, we have an approximate absolute error equal to 10^{-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{7}}} for $\Lambda^{8}$ against $10^{-4}$ for $\Lambda^{10}$.

We present now the behavior of Algorithm 4.1. Simulations are performed with a bulk-chasing parameter $\theta=0.5$. At each step $n$ of Algorithm 4.1, we use Algorithm 2.3 with $(\Delta t,M)=(10^{-4},1000)$, stopped when a stagnation is detected. As shown on the left plot of Figure 3, for \#\Lambda_{n}=5{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{5}} we reach approximately a precision of $10^{-14}$ as for Algorithm 3.2 performed on the exact solution (see Table 1). According to the right plot of Figure 3, we also observe that the enrichment procedure behaves similarly for both algorithms ($\tilde{\varepsilon}^{n}$ and $\varepsilon^{n}$ are almost the same). Here using the approximation provided by Algorithm 2.3 has a low impact on the behavior of Algorithm 3.2.

We present then results provided with the adaptive Algorithm 2.3. The parameters chosen for the adaptive interpolation are $\varepsilon=5\cdot 10^{-2}$, $\theta=0.5$. $K=30$ ensures the stopping of Algorithm 2.3. As illustrated by Figure 4, we recover globally the same behavior as for Algorithm 2.3 without adaptive interpolation. Indeed as $k$ increases, both absolute errors in $L^{2}$-norm and $L^{\infty}$-norm decrease and then stagnate. Again, we notice the influence of $\Delta t$ on the stagnation level. Nevertheless, the convergence rates are deteriorated and the algorithm provides less accurate approximations than Algorithm 4.1. This might be due to the sparse adaptive interpolation procedure, which uses here pointwise evaluations based on Monte-Carlo estimates, unlike Algorithm 4.1 which relies on pointwise evaluations resulting from Algorithm 2.3 stopping for a given tolerance.

Finally in Table 2, we compare the algorithmic complexity of these algorithms to reach a precision of $3\cdot 10^{-5}$ for $(\Delta t,M)=(10^{-4},1000)$. For adaptive Algorithm 2.3, $\Lambda_{k}$ refers to the set of multi-indices considered at step $k$ of Algorithm 2.3. For Algorithm 4.1, $N_{n}$ stands for the number of iteration required by Algorithm 2.3 to reach tolerance $\varepsilon_{tol}$ at step $n$. Finally, Algorithm 2.3 is run with full-grid $\Lambda=\Lambda_{max}$ where $\Lambda_{max}=\{\nu\in{\mathbb{N}}^{d}~{}:~{}\nu_{i}\leq 10\}$ is the set of multi-indices allowing to reach the machine precision. In this case, $N$ stands for the number of steps for this algorithm to converge.

We observe that both the adaptive version of Algorithm 2.3 and Algorithm 4.1 have a similar complexity, which is better than for the full-grid version of Algorithm 2.3. Moreover, we observed that while adaptive version of Algorithm 2.3 stagnates at a precision of $3\cdot 10^{-5}$, Algorithm 4.1, with the same parameters $\Delta t$ and $M$, converges almost up to the machine precision. This is why the high-dimensional test cases will be run only with Algorithm 4.1.

Higher-dimensional test cases. Now, we consider two other test cases noted respectively (TC2) and (TC3) in higher dimension.

1. (TC2)

As second test case in dimension $d=10$, we define (15) such that its solution is the Henon-Heiles potential

[TABLE]

We set (\Delta t,M)=({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}10^{-4}},1000) and $K=30$ for Algorithm 2.3. 2. (TC3)

We also consider the problem (15) whose exact solution is a sum of non-polynomial functions, like (TC1) but now in dimension $d=20$, given by

[TABLE]

Here, the Monte-Carlo simulations are performed for $(\Delta t,M)=(10^{-4},1000)$ and $K=30$.

Since for both test cases the exact solution is known, we propose to compare the behavior of Algorithm 4.1 and Algorithm 3.2. Again, the approximations $\tilde{u}_{n}$, at each step $n$ of Algorithm 4.1, are provided by Algorithm 2.3 stopped when a stagnation is detected. In both cases, the parameters for Algorithm 4.1 are set to $\theta=0.5$ and $\varepsilon=10^{-15}$.

In Table 3 and Table 4, we summarize the results associated to the exact and perturbed sparse adaptive algorithms for (TC2) and (TC3) respectively. We observe that Algorithm 4.1 performs well in comparison to Algorithm 3.2, for (TC2). Indeed, we get an approximation with a precision below the prescribed value $\varepsilon$ for both algorithms.

Similar observation holds for (TC3) in Table 4 and this despite the fact that the test case involves higher dimensional problem.

5 Conclusion

In this paper we have introduced a probabilistic approach to approximate the solution of high-dimensional elliptic PDEs. This approach relies on adaptive sparse polynomial interpolation using pointwise evaluations of the solution estimated using a Monte-Carlo method with control variates.

Especially, we have proposed and compared different algorithms. First we proposed Algorithm 2.3 which combines the sequential algorithm proposed in [14] and sparse interpolation. For the non-adaptive version of this algorithm we recover the convergence up to a threshold as the original sequential algorithm [15]. Nevertheless it remains limited to small-dimensional test cases, since its algorithmic complexity remains high. Hence, for practical use, the adaptive Algorithm 2.3 should be preferred. Adaptive Algorithm 2.3 converges but it does not allow to reach low precision with reasonable number of Monte-Carlo samples or time-step in the Euler-Maruyama scheme. Secondly, we proposed Algorithm 4.1. It is a perturbed sparse adaptive interpolation algorithm relying on inexact pointwise evaluations of the function to approximate. Numerical experiments have shown that the perturbed algorithm (Algorithm 4.1) performs well in comparison to the ideal one (Algorithm 3.2) and better than the adapted Algorithm 2.3 with a similar algorithmic complexity. Here, since only heuristic tools have been provided to justify the convergence of this algorithm, the proof of its convergence, under assumptions on the class of functions to be approximated, should be addressed in a future work.