Optimal Monte Carlo Methods for $L^2$-Approximation

TL;DR

This paper develops Monte Carlo algorithms for $L^2$-approximation in Hilbert spaces that achieve optimal convergence rates using limited function samples, matching the performance of any linear information-based method.

Contribution

It introduces Monte Carlo methods that attain optimal error rates for $L^2$-approximation with minimal function evaluations, matching the best linear information algorithms.

Findings

Monte Carlo methods achieve optimal convergence rates.

Error behavior matches that of any linear information sampling algorithm.

Methods are effective for multivariate functions in Hilbert spaces.

Abstract

We construct Monte Carlo methods for the $L^2$-approximation in Hilbert spaces of multivariate functions sampling no more than $n$ function values of the target function. Their errors catch up with the rate of convergence and the preasymptotic behavior of the error of any algorithm sampling $n$ pieces of arbitrary linear information, including function values.

Equations

$$e^{\mathrm{ran}}(A_n)=\underset{f\in F_{\circ }}{\sup }{\left(\mathbb{E}{\int }_D{\left|f-A_n(f)\right|}^2\text{leavevmode}\mathrm{d}\mu \right)}^{1/2}.$$

$$e(n)=\underset{A_n}{\inf }{e}^{\mathrm{ran}}(A_n).$$

$$e(n)\ge \frac{\sigma (2n)}{\sqrt{2}}$$

$$\sigma (n)\preccurlyeq n^{-p}(\ln n{)}^q\Rightarrow e(n)\preccurlyeq n^{-p}(\ln n{)}^q(\ln \ln n{)}^{p+1/2}$$

$$\sigma (n)\preccurlyeq n^{-p}\Rightarrow e(n)\preccurlyeq n^{-p}$$

$$L(m)\ge 2^{-r}L(n)\text{whenever}n\le m\le 2n.$$

$$L(n)=n^{-p}{\left(1+{\log }_{2}n\right)}^q$$

$$L(n)={\left(1+{\log }_{2}n\right)}^{-q}$$

$$\sigma (n)\preccurlyeq L(n)\Rightarrow e(n)\preccurlyeq L(n).$$

$$\sigma (n)$$

$$e^{\mathrm{ran}}(A_n^r)$$

$${\int }_{D}f\text{leavevmode}\mathrm{d}\mu ={\int }_D{A}_n^{r}f\text{leavevmode}\mathrm{d}\mu +{\int }_D(f-A_n^{r}f)\text{leavevmode}\mathrm{d}\mu .$$

$${⟨f,g⟩}_2={\int }_{D}f\cdot \overline{g}\text{leavevmode}\mathrm{d}\mu .$$

$$T:F\to L^2,Tf=f$$

$$e^{\det }(A)=\underset{f\in F_{\circ }}{\sup }{\Vert f-A(f)\Vert }_2.$$

$$A:F\times \Omega \to L^2,(f,\omega )\to A^{\omega }(f),$$

$$e^{\mathrm{ran}}(A)=\underset{f\in F_{\circ }}{\sup }{\left({\mathbb{E}}_{\omega }{\Vert f-A^{\omega }(f)\Vert }_2^2\right)}^{1/2}.$$

$$e^{\det }(n,T,{\Lambda }^{\mathrm{all}})=\underset{A\in {\mathcal{A}}_n^{\det ,\mathrm{all}}}{\inf }{e}^{\det }(A)$$

$$e^{\mathrm{ran}}(n,T,{\Lambda }^{\mathrm{all}})=\underset{A\in {\mathcal{A}}_n^{\mathrm{ran},\mathrm{all}}}{\inf }{e}^{\mathrm{ran}}(A)$$

$$e^{\det }(n,T,{\Lambda }^{\mathrm{std}})=\underset{A\in {\mathcal{A}}_n^{\det ,\mathrm{std}}}{\inf }{e}^{\det }(A)$$

$$e(n)=e^{\mathrm{ran}}(n,T,{\Lambda }^{\mathrm{std}})=\underset{A\in {\mathcal{A}}_n^{\mathrm{ran},\mathrm{std}}}{\inf }{e}^{\mathrm{ran}}(A)$$

$$e^{\mathrm{ran}}(n,T,{\Lambda }^{\mathrm{std}})\ge e^{\mathrm{ran}}(n,T,{\Lambda }^{\mathrm{all}})\ge \frac{1}{\sqrt{2}}{e}^{\det }(2n-1,T,{\Lambda }^{\mathrm{all}}).$$

$$P_n:F\to L^2,P_n(f)=\sum _{k=1}^{n\wedge N}{⟨f,b_k⟩}_2{b}_k.$$

$$u_m=\frac{1}{m}\sum _{j=1}^m{\left|b_j\right|}^2.$$

$${\mu }_m:\mathcal{A}\to [0,1],{\mu }_m(E)={\int }_E{u}_m\text{leavevmode}\mathrm{d}\mu$$

$$M_{\mathbf{n},\mathbf{m}}^{(k)}f=M_{\mathbf{n},\mathbf{m}}^{(k-1)}f+\sum _{j=1}^{m_k}\left[\frac{1}{n_k}\sum _{i=1}^{n_k}\frac{\left(f-M_{\mathbf{n},\mathbf{m}}^{(k-1)}f\right)\overline{b_j}}{u_{m_k}}\left(X_i^{(k)}\right)\right]b_j.$$

$$\sigma (m_k+1{)}^2\le e^{\mathrm{ran}}{\left(M_{\mathbf{n},\mathbf{m}}^{(k)}\right)}^2\le \frac{m_k}{n_k}{e}^{\mathrm{ran}}{\left(M_{\mathbf{n},\mathbf{m}}^{(k-1)}\right)}^2+\sigma (m_k+1{)}^2.$$

$${\mathbb{E}}_{\left\{1\dots k\right\}}{\Vert f-M_{\mathbf{n},\mathbf{m}}^{(k)}f\Vert }_2^2=\sum _{j=1}^N{\mathbb{E}}_{\left\{1\dots k\right\}}{\left|{⟨f-M_{\mathbf{n},\mathbf{m}}^{(k)}f,b_j⟩}_2\right|}^2.$$

$$\sum _{j=m_k+1}^N{\mathbb{E}}_{\left\{1\dots k\right\}}{\left|{⟨f-M_{\mathbf{n},\mathbf{m}}^{(k)}f,b_j⟩}_2\right|}^2=\sum _{j=m_k+1}^N{\left|{⟨f,b_j⟩}_2\right|}^2=\sum _{j=m_k+1}^N{\left|{⟨f,W{b}_j⟩}_F\right|}^2$$

$$=\sum _{j=m_k+1}^N{\left|{⟨f,\sigma (j)b_j⟩}_F\right|}^2\sigma (j{)}^2\le \sigma (m_k+1{)}^2{\Vert f\Vert }_F^2\le \sigma (m_k+1{)}^2.$$

$$g_j=\frac{\left(f-M_{\mathbf{n},\mathbf{m}}^{(k-1)}f\right)\overline{b_j}}{u_{m_k}}$$

Full text

AMS classification: 41A25, 41A63, 65C05, 65D15, 65D30, 68Q25, 65Y20.

Key words: Approximation of multivariate functions, Monte Carlo methods, optimal order of convergence, preasymptotic estimates, multivariate integration.

1 Introduction

Assume we want to approximate an unknown real or complex valued function on a set $D$ based on a finite number $n$ of function values which may be evaluated at randomly and adaptively chosen points. In general, these function values do not determine the function uniquely and so we cannot expect our approximation to be correct. We make an approximation error which we measure in the space $L^{2}(D,\mathcal{A},\mu)$ of quadratically integrable functions on $D$ with respect to an arbitrary measure $\mu$. In order to make any meaningful statement regarding this error, we need to have additional a priori knowledge of the unknown function. Here, we assume structural knowledge of the form that it is contained in the unit ball $F_{\circ}$ of a Hilbert space $F$ which is compactly embedded in $L^{2}(D,\mathcal{A},\mu)$. For instance, it may be bounded with respect to some Sobolev norm on a compact manifold $D$. The error of the randomized algorithm or Monte Carlo method $A_{n}$ is the quantity

[TABLE]

The error of an optimal randomized algorithm that ask for at most $n$ function values is denoted by

[TABLE]

While it seems impossible to provide such algorithms, the optimal deterministic algorithm evaluating $n$ arbitrary linear functionals is well known. It is given by the orthogonal projection $P_{n}$ onto the span of the first $n$ functions in the singular value decomposition of the embedding $T:F\hookrightarrow L^{2}$. Its worst case error is the $(n+1)$-st largest singular value or approximation number $\sigma(n+1)$ of that embedding, the square root of the $(n+1)$-st largest eigenvalue of the operator $W=T^{*}T$.

The algorithm $P_{n}$ asks for the first $n$ coefficients of $f$ with respect to the singular value decomposition of the embedding $T$. In most applications, however, it is not possible to sample these coefficients and we may only make use of function values. This leads to the following questions:

• •

How does the error $e(n)$ of optimal randomized algorithms using $n$ function values compare to the the error $\sigma(n+1)$ of the orthogonal projection $P_{n}$?

• •

If possible, find a randomized algorithm $A_{n}$ whose error is close to $\sigma(n+1)$.

These are not new questions in the fields of Monte Carlo methods and information-based complexity. There are several results for particular spaces $F$ where $e(n)$ behaves similarly to the error of $P_{n}$. See, for instance, Traub, Wasilkowski and Woźniakowski [20], Mathé [13] and Heinrich [5]. Results by Cohen, Davenport and Leviatan [2] and Cohen and Migliorati [3] contain a similar message, see Remark 3. In 1992, Novak [16] proved that

[TABLE]

holds for arbitrary spaces $F$. This means that optimal randomized algorithms using $n$ function values are never much better than the orthogonal projection $P_{n}$. On the other hand, Wasilkowski and Woźniakowski [23] proved in 2006 that

[TABLE]

for all $p>0$ and $q\geq 0$. Here, we write $x_{n}\preccurlyeq y_{n}$ if there is some $C>0$ and $n_{0}\in\mathbb{N}$ such that $x_{n}\leq Cy_{n}$ for all $n\geq n_{0}$. If $x_{n}\preccurlyeq y_{n}$ and $y_{n}\preccurlyeq x_{n}$, we write $x_{n}\asymp y_{n}$. This means that optimal randomized algorithms using function values are always almost as good as the orthogonal projection $P_{n}$. The proof of this result is constructive. It raises the question whether the additional power of the double logarithm is necessary or not. In fact, Novak and Woźniakowski showed in 2012 that this is not the case for $q=0$, that is

[TABLE]

for all $p>0$. The proof of this result, however, is not constructive. Both proofs can be found in their monograph [18, Chapter 22]. In the present paper, we prove the corresponding statement for $q>0$. More generally, we consider upper bounds with the following property. We say that the sequence $L:\mathbb{N}\to(0,\infty)$ is regularly decreasing if there is some $r\geq 0$ such that

[TABLE]

If there is some $n_{0}\in\mathbb{N}$ such that $L(n)$ is nonincreasing for $n\geq n_{0}$, this is equivalent to $L(2n)\asymp L(n)$. Property (1) is satisfied if $L(n)n^{r}$ is nondecreasing. The sequence

[TABLE]

is regularly decreasing for any $p>0$ and $q\geq 0$. It satisfies (1) for $r=p$. Another example is

[TABLE]

for any $q>0$, which satisfies (1) for $r=q$. The sequence is not regularly decreasing if it decays exponentially or has huge jumps. We obtain the following result.

Theorem 1**.**

If $L:\mathbb{N}\to(0,\infty)$ is regularly decreasing, then

[TABLE]

This solves Open Problem 99 as posed by Novak and Woźniakowski in [18]. One problem with this result is that it does not provide any algorithm, it only states the existence of good algorithms. Another problem is that the error bound is only asymptotic. The preasymptotic behavior of $e(n)$ may, however, be very different from its asymptotic behavior. This is typically the case if the set $D$ is a domain in high dimensional euclidean space.

These problems are tackled by Theorem 2. In Section 3, we provide a randomized algorithm $A_{n}^{r}$ for any $n\in\mathbb{N}$ and $r\geq 0$. This algorithm is a refinement of the algorithm proposed by Wasilkowski and Woźniakowski [23]. It asks for at most $n$ function values and satisfies the following error bound.

Theorem 2**.**

Assume that $L:\mathbb{N}\to(0,\infty)$ satisfies (1) and let $c_{r}=2^{r\lceil 2r+3\rceil+1}$.

[TABLE]

The constant $c_{r}$ only depends on the order $r$. If $D$ is a domain in $d$-dimensional euclidean space, this order is often independent of $d$ or even strictly decreasing with $d$. See Section 3 for the definition of this algorithm and several examples.

We find that the error of randomized algorithms using $n$ function values of the target function can get very close to the error of the orthogonal projection $P_{n}$ and that this is achieved by the algorithm $A_{n}^{r}$.

In Section 4, we use these algorithms for the integration of functions $f$ in $F$ with respect to probability measures $\mu$. We simply exploit the relation

[TABLE]

We compute the integral of $A_{n}^{r}f$ and use a direct simulation to approximate the integral of $(f\,$–$\,A_{n}^{r}f)$, which has a small variance. This technique is called variance reduction and widely used for Monte Carlo integration. See Heinrich [5, Theorem 5.3] for another example. Even if $D$ is a high dimensional domain, the resulting method can significantly improve on the error of a sole direct simulation for a relatively small number of samples.

These results are based on the a priori knowledge that our target function is contained in the unit ball of the space $F$. In Section 5, we discuss how this assumption can be weakened.

2 The Setting

Let $(D,\mathcal{A},\mu)$ be a measure space and $\mathbb{K}\in\left\{\mathbb{R},\mathbb{C}\right\}$. The space $L^{2}=L^{2}(D,\mathcal{A},\mu)$ is the space of quadratically integrable $\mathbb{K}$-valued functions on $(D,\mathcal{A},\mu)$, equipped with the scalar product

[TABLE]

Let $F$ be a second Hilbert space and $F_{\circ}$ be its unit ball. We assume that $F$ is a subset of $L^{2}$ and that

[TABLE]

is compact. With the embedding $T$ we associate a positive semi-definite and compact operator $W=T^{*}T$ on the space $F$. By the spectral theorem, there is a (possibly finite) orthogonal basis $\mathcal{B}=\left\{b_{1},b_{2},\dots\right\}$ of $F$, consisting of eigenvectors corresponding to a nonincreasing zero sequence $(\lambda_{n})_{n\in\mathbb{N}}$ of eigenvalues of $W$. Let $N$ be the cardinality of $\mathcal{B}$. One can easily check that $\mathcal{B}$ is orthogonal in $L^{2}$, as well. We take the eigenvectors $b_{n}$ to be normalized in $L^{2}$. We call this basis the singular value decomposition of $T$.111 This term is more commonly used to refer to the representation $Tf=\sum_{b\in\mathcal{B}}\left\langle f,b\right\rangle Tb$ of the compact operator. Here, the altered terminology shall ease the notation. The number $\sigma(n)=\sqrt{\lambda_{n}}$ is called its $n$-th singular value or approximation number.

The worst case error of a deterministic algorithm $A:F\to L^{2}$ is the quantity

[TABLE]

The worst case error of a measurable randomized algorithm

[TABLE]

where $\Omega$ is the sample space of some probability space $(\Omega,\mathcal{F},\mathbb{P})$, is the quantity

[TABLE]

We usually skip the $\omega$ in the notation. See Novak and Woźniakowski [17, Chapter 4] for a precise definition of such algorithms. We furthermore define the following minimal worst case errors within certain classes of algorithms.

The quantity

[TABLE]

is the minimal worst case error within the class $\mathcal{A}_{n}^{\rm det,all}$ of all deterministic algorithms evaluating at most $n$ linear functionals of the input function.

The quantity

[TABLE]

is the minimal worst case error within the class $\mathcal{A}_{n}^{\rm ran,all}$ of all measurable randomized algorithms evaluating at most $n$ linear functionals.

The quantity

[TABLE]

is the minimal worst case error within the class $\mathcal{A}_{n}^{\rm det,std}$ of all deterministic algorithms evaluating at most $n$ function values of the input function.

The quantity

[TABLE]

finally is the minimal worst case error within the class $\mathcal{A}_{n}^{\rm ran,std}$ of all measurable randomized algorithms evaluating at most $n$ function values. This is the error to be analyzed. It was proven by Novak [16] that

[TABLE]

The error $e^{\rm det}(n,T,\Lambda^{\rm all})$ is known to coincide with $\sigma(n+1)$. We refer to Novak and Woźniakowski [17, Section 4.2.3]. The infimum is attained for the nonadaptive linear algorithm

[TABLE]

Here, $\log_{2}x$ denotes the logarithm of $x>0$ in base 2, whereas $\ln x$ denotes its natural logarithm. The minimum of $a\in\mathbb{R}$ and $b\in\mathbb{R}$ is denoted by $a\land b$. Recall that we write $x_{n}\preccurlyeq y_{n}$, if there is a positive constant $C$ and some $n_{0}\in\mathbb{N}$ such that $x_{n}\leq Cy_{n}$ for all $n\geq n_{0}$. We write $x_{n}\asymp y_{n}$ if $x_{n}\preccurlyeq y_{n}$ and $y_{n}\preccurlyeq x_{n}$.

3 A Method for Multivariate Approximation

Let us keep the notation of the previous section. For any $m\in\mathbb{N}$ with $m\leq N$, we define

[TABLE]

This is a probability density with respect to $\mu$. We consider the probability measure

[TABLE]

on $(D,\mathcal{A})$. In view of optimal algorithms in $\mathcal{A}_{n}^{\rm det,all}$, we introduce the following family of algorithms in $\mathcal{A}_{n}^{\rm ran,std}$.

Algorithm**.**

Let $\boldsymbol{n}=\left(n_{1},n_{2},\dots\right)$ and $\boldsymbol{m}=\left(m_{1},m_{2},\dots\right)$ be sequences of nonnegative integers such that $\boldsymbol{m}$ is nondecreasing and bounded above by $N=\left|\mathcal{B}\right|$. We define the algorithms $M^{(k)}_{\boldsymbol{n},\boldsymbol{m}}:L^{2}\to L^{2}$ for $k\in\mathbb{N}_{0}$ as follows.

• •

Set $M^{(0)}_{\boldsymbol{n},\boldsymbol{m}}=0$.

• •

For $k\geq 1$ and $f\in L^{2}$, let $X_{1}^{(k)},\ldots,X_{n_{k}}^{(k)}$ be random variables with distribution $\mu_{m_{k}}$ that are each independent of all the other random variables and set

[TABLE]

Note that the expectation of each term in the inner sum is $\langle f-M^{(k-1)}_{\boldsymbol{n},\boldsymbol{m}}f,b_{j}\rangle_{2}$. The algorithm $M^{(k)}_{\boldsymbol{n},\boldsymbol{m}}$ hence approximates $f$ in $k$ steps. In the first step, $n_{1}$ function values of $f$ are used for standard Monte Carlo type approximations of its $m_{1}$ leading coefficients with respect to the orthonormal system $\mathcal{B}$. In the second step, $n_{2}$ values of the residue are used for standard Monte Carlo type approximations of its $m_{2}$ leading coefficients and so on. In total, $M^{(k)}_{\boldsymbol{n},\boldsymbol{m}}$ uses $\sum_{j=1}^{k}n_{j}$ function values of $f$. The total number of approximated coefficients is $m_{k}$.

Algorithms of this type have already been studied by Wasilkowski and Woźniakowski in [23]. The simple but crucial difference with the above algorithms is the variable number $n_{j}$ of nodes in each approximation step. Note that this stepwise approximation is similar to several multilevel Monte Carlo methods as introduced by Heinrich in 1998, see [4].

The benefit from the $k$-th step is controlled by $m_{k}$ and $n_{k}$ as follows.

Lemma 1**.**

For all nondecreasing sequences $\boldsymbol{n}$ and $\boldsymbol{m}$ of nonnegative integers and all $k\in\mathbb{N}$, we have

[TABLE]

Lemma 1 corresponds to Theorem 22.14 by Novak and Woźniakowski [18]. The setting of the present paper is slightly more general, but the proof is the same. Since Lemma 1 is essential for the following investigation, I present the proof.

Proof.

The lower bound holds true, since $M^{(k)}_{\boldsymbol{n},\boldsymbol{m}}(b_{m_{k}+1})$ is perpendicular to $b_{m_{k}+1}$. To prove the upper bound, let $f\in F_{\circ}$. By $\mathbb{E}_{I}$ we denote the expectation with respect to the random variables $X_{i}^{(j)}$ for $j\in I$ and $i=1\ldots n_{j}$. We need to estimate

[TABLE]

On the one hand, we have

[TABLE]

We use the abbreviation

[TABLE]

for each $j\leq m_{k}$. Note that $u_{m_{k}}=0$ implies $b_{j}=0$ and we set $g_{j}=0$ in this case. We then obtain on the other hand for each $j\leq m_{k}$ that

[TABLE]

and hence

[TABLE]

With Fubini’s theorem this yields that

[TABLE]

and the upper bound is proven. ∎

We now define the algorithm of Theorem 2. We consider such algorithms $M^{(k)}_{\boldsymbol{n},\boldsymbol{m}}$, where the number of nodes $n_{j}$ is doubled in each step and the ratio $\frac{m_{j}}{n_{j}}$ of approximated coefficients and computed function values is constant, say $2^{-\ell}$. This way, the total number $m_{k}$ of approximated coefficients is linear in the total number $n$ of computed function values. This is necessary to achieve an error of the same order as with optimal algorithms using arbitrary linear information, which precisely compute the first $n$ coefficients. The algorithms by Wasilkowski and Woźniakowski [23] do not have this property. If the ratio is small enough, Lemma 1 ensures that $M^{(k)}_{\boldsymbol{n},\boldsymbol{m}}$ inherits optimal error bounds from $M^{(k-1)}_{\boldsymbol{n},\boldsymbol{m}}$.

Algorithm**.**

Given $r\geq 0$, we set $\ell_{r}=\lceil 2r+1\rceil$ and define the sequences $\boldsymbol{n}$ and $\boldsymbol{m}$ by

[TABLE]

For $n\in\mathbb{N}$, we choose $k\in\mathbb{N}_{0}$ such that $2^{k}\leq n<2^{k+1}$ and set

[TABLE]

The algorithm $A_{n}^{r}$ obviously performs less than $n$ function evaluations.

Proof of Theorem 2.

Let $\boldsymbol{n}$ and $\boldsymbol{m}$ be defined as above and $k\in\mathbb{N}_{0}$. We first show that

[TABLE]

where $\bar{c}_{r}=2^{r(\ell_{r}+1)+1}$. We use induction on $k$. If $k\leq\ell_{r}$, we have $M^{(k)}_{\boldsymbol{n},\boldsymbol{m}}=0$ and

[TABLE]

For $k>\ell_{r}$, we inductively obtain with Lemma 1 that

[TABLE]

where the term in brackets is smaller than 1. This shows (3). For $n\in\mathbb{N}$, we choose $k\in\mathbb{N}_{0}$ with $2^{k}\leq n<2^{k+1}$ and obtain

[TABLE]

as it was to be proven. ∎

Note that Theorem 1 is a direct consequence of Theorem 2. Of course, the best possible upper bound for $\sigma(n)$ is $\sigma(n)$ itself. If we combine Theorem 1 for $L(n)=\sigma(n)$ with Novak’s lower bound (2), we obtain the following statement on the order of convergence.

Corollary 1**.**

Assume that $\sigma(2n)\asymp\sigma(n)$. Then

[TABLE]

Note that the error $e^{\rm det}\left(n,F\hookrightarrow L^{2},\Lambda^{\rm std}\right)$ of optimal deterministic algorithms based on function values may perform much worse, as shown by Hinrichs, Novak and Vybíral [7], see also Novak and Woźniakowski [18, Section 26.6.1]. It is a very interesting question whether the condition on the decay of the singular values can be relaxed. Note that we use this condition both to prove the upper and the lower bound of Corollary 1. On the other hand, if we combine Theorem 2 for $L(n)=\sigma(n)$ and the lower bound (2), we obtain the following optimality result.

Corollary 2**.**

Assume that there is some $r\geq 0$ such that $\sigma(2n)\geq 2^{-r}\sigma(n)$ holds for all $n\in\mathbb{N}$. We set $\tilde{c}_{r}=2^{r\lceil 2r+4\rceil+3/2}$. Then we have

[TABLE]

Let us now consider some examples. In each example, we first discuss the order of convergence of $e^{\rm ran}\left(n,F\hookrightarrow L^{2},\Lambda^{\rm std}\right)$. We then talk about explicit upper bounds.

Example 1** (Approximation of mixed order Sobolev functions on the torus).**

Let $D$ be the $d$-dimensional torus $\mathbb{T}^{d}$, represented by the unit cube $[0,1]^{d}$, where opposite faces are identified. Let $\mathcal{A}$ be the Borel $\sigma$-algebra on $\mathbb{T}^{d}$ and $\mu$ the Lebesgue measure. Let $F$ be the Sobolev space of complex valued functions on $D$ with dominating mixed smoothness $r\in\mathbb{N}$, equipped with the scalar product

[TABLE]

We know that

[TABLE]

This classical result goes back to Babenko [1] and Mityagin [14]. Corollary 1 yields

[TABLE]

This is a new result. The optimal order is achieved by the algorithm $A_{n}^{r}$ and the author does not know of any other algorithm with this property. It is still an open problem whether the same rate can be achieved with deterministic algorithms based on function values. So far, it is only known that

[TABLE]

The upper bound is achieved by Smolyak’s algorithm, see Sickel and Ullrich [19].

We now turn to explicit estimates. We know that there is some $C_{r,d}>0$ such that

[TABLE]

This upper bound is optimal as $n$ tends to infinity. However, it is not useful to describe the error numbers for small values of $n$. Simple calculus shows that the right hand side in (5) is increasing for $n\leq e^{d-1}$. The error numbers, on the other hand, are decreasing. Moreover, the right hand side attains its minimum for $n=2$ if restricted to $n\leq(d-1)^{d-1}$ and is hence larger than $e^{\rm ran}\left(2,F\hookrightarrow L^{2},\Lambda^{\rm std}\right)$. This means that the trivial upper bound

[TABLE]

is better than (5) for all $n\leq(d-1)^{d-1}$ and regardless of the value of $C_{r,d}$. For these reasons, it is important to consider different error bounds, if the dimension $d$ is large. See also the paper of Kühn, Sickel and Ullrich [12]. Based on this paper, it is shown by the author [9] that

[TABLE]

We obtain with Theorem 2 that

[TABLE]

Example 2** (Approximation of mixed order Sobolev functions on the cube).**

Now, let $D$ be the $d$-dimensional unit cube $[0,1]^{d}$ with the induced topology and let $\mathcal{A}$ be the Borel $\sigma$-algebra and $\mu$ the Lebesgue measure. Let $F$ be the Sobolev space of complex valued functions on $[0,1]^{d}$ with dominating mixed smoothness $r\in\mathbb{N}$, equipped with the scalar product (4). Just like on the torus, we have

[TABLE]

where the optimal rate is achieved by $A_{n}^{r}$. Like in Example 1, the corresponding upper bounds are bad for $n\leq(d-1)^{d-1}$. In this range, we need different estimates for the approximation numbers. It is known that

[TABLE]

This estimate cannot be improved significantly for $n\leq 2^{d}$, even if $r=\infty$. See the author’s paper [9] for more details. With Theorem 2, we obtain the upper bound

[TABLE]

Example 3** (Approximation in tensor product spaces).**

This example is more general than the previous ones. By $H_{1}\otimes H_{2}$ we denote the tensor product of two Hilbert spaces $H_{1}$ and $H_{2}$. For $j=1\ldots d$ let $(D_{j},\mathcal{A}_{j},\nu_{j})$ be a $\sigma$-finite measure space and $F_{j}$ be a Hilbert space of $\mathbb{K}$-valued functions which is compactly embedded in $L^{2}(D_{j},\mathcal{A}_{j},\nu_{j})$. The $\sigma$-finity of the measure spaces ensures that

[TABLE]

where $D$ is the Cartesian product of the sets $D_{j}$ and $\mu$ is the unique product measure of the measures $\nu_{j}$ on the tensor product $\mathcal{A}$ of the $\sigma$-algebras $\mathcal{A}_{j}$. The tensor product space

[TABLE]

is compactly embedded in $L^{2}(D,\mathcal{A},\mu)$. Assuming that the approximation numbers of the univariate embeddings $F_{j}\hookrightarrow L^{2}(D_{j},\mathcal{A}_{j},\nu_{j})$ are of polynomial decay, that is

[TABLE]

for some $r_{j}>0$, it can be derived from Mityagin [14] and Nikol’skaya [15] that

[TABLE]

where $r$ is the minimum among all numbers $r_{j}$ and $d_{0}$ is its multiplicity. Corollary 1 implies

[TABLE]

where the optimal order is achieved by $A_{n}^{r}$. We do not discuss explicit estimates in this abstract setting.

Example 4** (Approximation of isotropic Sobolev functions on the torus).**

Let $D$ again be the $d$-torus, this time represented by $[0,2\pi]^{d}$. Let $F$ be the Sobolev space of complex valued functions on $D$ with isotropic smoothness $r\in\mathbb{N}$, equipped with the scalar product

[TABLE]

This example is not a tensor product problem. For this classical problem, it is known that

[TABLE]

for $r>d/2$. In the case $r\leq d/2$, where function values are only defined almost everywhere, the last three relations stay valid. See Jerome [8], Triebel [21], Mathé [13] and Heinrich [6]. For $n\leq 2^{d}$, however, the function $n^{-r/d}$ is not suited to describe the behavior of $\sigma(n)$. It has been proven by Kühn, Mayer and Ullrich [11] that there are positive constants $b_{r}$ and $B_{r}$ that do not depend on $d$ such that

[TABLE]

for all $d>1$ and $n\in\mathbb{N}$ with $d\leq n\leq 2^{d}$. If we apply Relation (2) and Theorem 2222We take $L(n)$ as the right hand side in (7) for $d\leq n\leq 2^{d}$, $L(n)=L(2^{d})$ for $n>2^{d}$ and $L(n)=\max\left\{1,L(d)\right\}$ for $n<d$. Then $\sigma(n)\leq L(n)$ for $n\in\mathbb{N}$ and $L(n)n^{r}$ is nondecreasing., we obtain the existence of $d$-independent positive constants $\tilde{b}_{r}$ and $\widetilde{B}_{r}$ such that

[TABLE]

for all $d>1$ and $n\in\mathbb{N}$ with $d\leq n\leq 2^{d-1}$. This optimal behavior is achieved by the algorithm $A_{n}^{r}$.

Remark 1** (Implementation of these algorithms).**

The construction of the algorithms $A_{n}^{r}$ is completely explicit. We are able to implement these algorithms, if we know the singular value decomposition $\mathcal{B}$ of the embedding $F\hookrightarrow L^{2}$ and if we are able to sample from the probability distributions $\mu_{m}$. This task may be very hard. In Example 1 and 4, however, it is not. Here, $\mathcal{B}$ is the Fourier basis of $L^{2}$ and all the random variables are independent and uniformly distributed on the unit cube. Also the case of general tensor product spaces $F$ and $L^{2}$ can be handled, if the singular value decompositions $\mathcal{B}_{j}$ of the univariate embeddings $F_{j}\hookrightarrow L^{2}(D_{j},\mathcal{A}_{j},\nu_{j})$ are known. Then, the singular value decomposition of the embedding $F\hookrightarrow L^{2}$ is given by

[TABLE]

and the probability measure $\mu_{m}$ is the average of $m$ product densities, that is

[TABLE]

where ${\rm d}\eta_{i,j}=|b_{i,j}|^{2}{\rm d}\nu_{j}$ with some $b_{i,j}\in\mathcal{B}_{j}$. A random sample $x$ from this distribution can be obtained as follows:

• (1)

Get $i$ from the uniform distribution on $\left\{1,\dots,m\right\}$.

• (2)

Get $x_{1},\dots,x_{d}$ independently from the probability distributions $\eta_{i,1},\dots,\eta_{i,d}$.

The second step can for example be done by rejection sampling, if the measures $\eta_{i,j}$ have a bounded Lebesgue density. This way, the total sampling costs are linear in $d$. Another method of sampling from $\mu_{m}$ is proposed by Cohen and Migliorati in [3, Section 5].

4 A Method for Multivariate Integration

In this section, we require the measure $\mu$ to be finite. This ensures that the integral operator

[TABLE]

is well defined and continuous on $F$. Let us assume that $\mu$ is a probability measure. We want to approximate $I(f)$ for an unknown function $f\in F_{\circ}$ by a randomized algorithm $Q_{n}$ which evaluates at most $n$ function values of $f$. The worst case error of $Q_{n}$ is the quantity

[TABLE]

The minimal worst case error among such algorithms is denoted by

[TABLE]

Like any method for $L^{2}$-approximation, the algorithm $A_{n}^{r}$ from Section 3 can also be used for numerical integration.

Algorithm**.**

For all $r>0$, any $n\in\mathbb{N}$ and $f\in L^{2}$, let

[TABLE]

where $X_{1},\dots,X_{n}$ are random variables with distribution $\mu$ which are independent of each other and the random variables in $A_{n}^{r}$.

It is easy to verify that $Q_{2n}^{r}$ is unbiased, evaluates at most $2n$ function values of $f$ and satisfies

[TABLE]

for each $f$ in $L^{2}$. We thus obtain the following corollary.

Corollary 3**.**

Assume that $L:\mathbb{N}\to(0,\infty)$ satisfies (1) and let $c_{r}=2^{r\lceil 2r+3\rceil+1}$.

[TABLE]

In particular:

[TABLE]

The result on the order of convergence is quite general but not always optimal. An example is given by integration with respect to the Lebesgue measure $\mu$ on the Sobolev space $F$ with dominating mixed smoothness $r$ on the $d$-dimensional unit cube, as treated by Novak and the author [10] and Ullrich [22]. In this case, we have

[TABLE]

The main strength of Corollary 3 is that it provides us with unbiased methods for high dimensional integration achieving a small error with a modest number of function values.

Example 5** (Integration of mixed order Sobolev functions on the torus).**

Like in Example 1, let $F$ be the Sobolev space of dominating mixed smoothness $r$ on the $d$-torus and let $\mu$ be the Lebesgue measure. Among all randomized algorithms for multivariate integration in $F$ the randomized Frolov algorithm $Q_{n}^{*}$ is known to have the optimal error rate. It is shown by Ullrich [22] that there is some constant $c>2^{d}$ such that

[TABLE]

However, this estimate is trivial, if $n$ is not exponentially large in $d$. For smaller values of $n$, an error less than one is guaranteed by the direct simulation

[TABLE]

with independent and uniformly distributed random variables $X_{j}$. It satisfies

[TABLE]

However, this error bound converges only slowly, as $n$ tends to infinity. It does not reflect the smoothness of the integrands at all. The above method also guarantees nontrivial error bounds for smaller values of $n$, but converges faster than $S_{n}$. Relation (6) immediately yields that

[TABLE]

with $p=\frac{r}{2+\ln d}$ and $C=2^{p\lceil 2p+4\rceil+1}$. For example, let $d=500$ and $r=8$. For one million function values, the estimate (8) for the Frolov algorithm is larger than one, the estimate (9) for the direct simulation gives the error $10^{-3}$ and the estimate (10) for our new algorithm gives an error smaller than $5\cdot 10^{-7}$.

Remark 2** (Implementation of these algorithms).**

We are able to implement the algorithms $Q_{2n}^{r}$ under the following assumptions:

• •

We are able to implement $A_{n}^{r}$. This issue is discussed in Remark 1.

• •

We know the integrals $I(b_{j})$ of the eigenfunctions $b_{j}\in\mathcal{B}$ for all $j\leq 2^{-\ell_{r}}n$.

• •

We can sample from the probability distribution $\mu$.

In the above example, the implementation is particularly easy, since $\mathcal{B}$ is the Fourier basis and all the random variables are independent and uniformly distributed on the unit cube.

5 A weaker type of a priori knowledge

In the previous sections, we assumed that the target function $f$ is contained in the unit ball of a Hilbert space $F$ which is compactly embedded into $L^{2}$, that is

[TABLE]

As we have seen in Section 2, the space $F$ induces a nonincreasing sequence $\sigma$, the singular numbers

[TABLE]

of the embedding $F\hookrightarrow L^{2}$. This sequence is either finite or tends to zero. It also induces a nested sequence $V$ of subspaces

[TABLE]

where $V_{m}$ is spanned by the first $m$ elements of the singular value decomposition.

In turn, any such pair $\left(\sigma,V\right)$ induces a Hilbert space $F$ which is compactly embedded into $L^{2}$. We choose $b_{m}$ as an element of the orthogonal complement of $V_{m-1}$ in $V_{m}$ with $\left\|b_{m}\right\|_{2}=1$ and define $F$ by its orthonormal basis $\{\sigma(1)b_{1},\sigma(2)b_{2},\ldots\}$. It has the scalar product

[TABLE]

where we take the sum over the whole sequence $\sigma$. It is not hard to see that the correspondence between $F$ and the pair $\left(\sigma,V\right)$ is bijective up to the choice of the spaces $V_{m}$ for which we have $\sigma(m+1)=\sigma(m)$.

Let $P_{m}$ denote the orthogonal projection onto $V_{m}$ in $L^{2}$. It is readily verified that our assumption (11) on the target function $f$ implies that

[TABLE]

In general, however, (12) is strictly weaker than (11). For example, if $\sigma(k)=1/k$ for $k\in\mathbb{N}$, the function

[TABLE]

satisfies (12) but is not even contained in the space $F$. In Section 3, we constructed a randomized algorithm $A_{n}^{r}:L^{2}\to V_{m}$ and proved upper bounds on the mean square error $\mathbb{E}\left\|f-A_{n}^{r}(f)\right\|_{2}^{2}$ for any $f$ from (11). In fact, the same error bounds hold for any $f$ from (12). We state this as Theorem 3.

Theorem 3**.**

Let $(D,\mathcal{A},\mu)$ be a measure space and $L^{2}=L^{2}(D,\mathcal{A},\mu)$. For any $m\in\mathbb{N}_{0}$ let $V_{m}$ be an $m$-dimensional subspace of $L^{2}$ such that $V_{m}\subset V_{m+1}$ and let $P_{m}:L^{2}\to V_{m}$ be the orthogonal projection onto $V_{m}$. Assume that $f\in L^{2}$ satisfies

[TABLE]

with some $\varepsilon:\mathbb{N}_{0}\to(0,\infty)$. Then the randomized algorithm $Q_{m}:L^{2}\to V_{m}$ as defined below satisfies

[TABLE]

for any $m=2^{k}$ and $k\in\mathbb{N}_{0}$. The number of requested function values is at most

[TABLE]

To define the algorithm $Q_{m}$ we choose $b_{n}$ in the orthogonal complement of $V_{n-1}$ in $V_{n}$ with $\left\|b_{n}\right\|_{2}=1$ for all $n\in\mathbb{N}$. For $j\in\mathbb{N}$, we set

[TABLE]

Then the method $M^{(k)}_{\boldsymbol{n},\boldsymbol{m}}:L^{2}\to V_{m_{k}}$ for $k\in\mathbb{N}_{0}$ can be defined as in Section 3. Given $m=2^{k}$ for some $k\in\mathbb{N}_{0}$, we define $Q_{m}=M^{(k+1)}_{\boldsymbol{n},\boldsymbol{m}}:L^{2}\to V_{m}$.

Proof.

We only sketch the proof since it is very similar to the proof of Theorem 2. Just like in Lemma 1, we can show for any $k\in\mathbb{N}_{0}$ that

[TABLE]

The statement follows by induction on $k\in\mathbb{N}_{0}$. ∎

Note that we did not impose any condition on the upper bound $\varepsilon:\mathbb{N}_{0}\to(0,\infty)$. If $\varepsilon$ is regularly decreasing, the maximum in (14) is bounded by a constant which does not depend on $m$. Roughly speaking, the algorithm $Q_{m}$ admits a mean square error of order $\varepsilon(m)$ with a sample size of order $m$ for any $f$ from (13).

Remark 3** (Optimal approximation within a subspace).**

Let $D$ be a Borel subset of $\mathbb{R}^{d}$ with positive Lebesgue measure, $\mathcal{A}$ be the Borel sigma algebra on $D$ and $\mu$ be a probability measure on $(D,\mathcal{A})$. The best approximation of $f\in L^{2}(D,\mathcal{A},\mu)$ within the subspace $V_{m}$ is given by $P_{m}f$. Its error is given by the number

[TABLE]

In general, we cannot find $P_{m}f$ by sampling only a finite number of function values of $f$. What we can provide, is a random approximation $v_{m}$ within $V_{m}$ whose root mean square error

[TABLE]

is close to $e_{m}(f)$. If we know the numbers $e_{m}(f)$ for all $m\in\mathbb{N}$ (or some good upper bound) and if they are regularly decreasing, we can choose $v_{m}$ as the output of the method $Q_{m}$ from Theorem 3, which uses a sample size of order $m$. But even if we do not know anything about $f\in L^{2}$, we can still find an approximation $v_{m}$ like above. We only need the mild assumption that $V_{m}$ consists of functions defined everywhere on $D$ and that for each $x\in D$, there is some $v\in V_{m}$ with $v(x)\neq 0$. We can then choose $v_{m}$ as the output of a weighted least squares method, see Cohen and Migliorati [3, Theorem 2.1 (iv)]. The sample size of this method, however, is at least of order $m\ln m$. In both cases, the involved proportionality constants are independent of the dimension of the domain $D$.

Acknowledgements**.**

I wish to thank Erich Novak, Robert Kunsch, Winfried Sickel and two anonymous referees, whose comments and questions led to the present generality of the theorems.