DeepAI

# Sub-optimality of Gauss–Hermite quadrature and optimality of trapezoidal rule for functions with finite smoothness

A sub-optimality of Gauss–Hermite quadrature and an optimality of the trapezoidal rule are proved in the weighted Sobolev spaces of square integrable functions of order α, where the optimality is in the sense of worst-case error. For Gauss–Hermite quadrature, we obtain the matching lower and upper bounds, which turns out to be merely of the order n^-α/2 with n function evaluations, although the optimal rate for the best possible linear quadrature is known to be n^-α. Our proof on the lower bound exploits the structure of the Gauss–Hermite nodes; the bound is independent of the quadrature weights, and changing the Gauss–Hermite weights cannot improve the rate n^-α/2. In contrast, we show that a suitably truncated trapezoidal rule achieves the optimal rate up to a logarithmic factor.

• 6 publications
• 4 publications
• 18 publications
12/22/2022

### Randomizing the trapezoidal rule gives the optimal RMSE rate in Gaussian Sobolev spaces

Randomized quadratures for integrating functions in Sobolev spaces of or...
01/10/2018

### Optimal functional supervised classification with separation condition

We consider the binary supervised classification problem with the Gaussi...
11/19/2019

### Optimal Complexity and Certification of Bregman First-Order Methods

We provide a lower bound showing that the O(1/k) convergence rate of the...
09/01/2019

### Central limit theorems for discretized occupation time functionals

The approximation of integral type functionals is studied for discrete o...
01/30/2023

### Infinite-Variate L^2-Approximation with Nested Subspace Sampling

We consider L^2-approximation on weighted reproducing kernel Hilbert spa...
05/30/2021

### Optimality of the recursive Neyman allocation

We derive a formula for the optimal sample allocation in a general strat...
06/29/2021

### Exponential Weights Algorithms for Selective Learning

We study the selective learning problem introduced by Qiao and Valiant (...

## 1 Introduction

This paper is concerned with a sub-optimality of Gauss–Hermite quadrature and an optimality of the trapezoidal rule.

Given a function , Gauss–Hermite quadrature is arguably one of the standard numerical integration methods to compute the integral

 ∫Rf(x)1√2πe−x2/2dx.

It is a Gauss-type quadrature, i.e., the quadrature points are the zeros of the degree- orthogonal polynomial associated with the weight function, and the corresponding quadrature weights are readily defined. With the weight function , the orthogonal polynomial we have is the (so-called probabilist’s) Hermite polynomial.

Gauss–Hermite quadrature is widely used; here, we just mention spectral methods [Guo.BY_1999_Hermite, mao2017hermite, Canuto.C_etal_2006_book], and applications in aerospace engineering [braun2021satellite, BBM2020book], finance [Brandimarte2006book, FR2008book, GH2010], and physics [SS2016book, Gezerlis2020book]. Nevertheless, the limitation of this method seems to be less known.

We start with numerical results that illustrate this deficiency. Figure 1 shows a comparison of Gauss–Hermite rule and a suitably truncated trapezoidal rule on . Here the target function is , and the integral is with respect to the standard Gaussian measure. For the trapezoidal rule, we integrate , upon a suitable cut-off of the domain. We discuss the setting of this experiment in more detail at the end of Section 4.

What we observe in Figure 1 is that while the trapezoidal rule achieves , Gauss–Hermite quadrature achieves only a slower convergence rate , where is the number of quadrature points. A similar empirical inefficiency of Gauss–Hermite quadrature is reported in a paper by one of the present authors and Nuyens [NS2021], also in a recent paper by Trefethen [T2021], who pointed out the problems of the quadrature weights decaying too fast relative to the widely-spread quadrature points.

In this paper, we prove a sub-optimality of Gauss–Hermite quadrature. More precisely, we establish a sharp lower bound for the error decay of Gauss–Hermite quadrature in the sense of worst case error. Moreover, we show that a suitably truncated trapezoidal rule achieves the optimal rate of convergence, up to a logarithmic factor.

Integrands of our interest are functions with finite smoothness. In this regard, we will work under the assumption that the integrand function lives in the weighted -Sobolev space with smoothness . This function space is widely used; see for example books [B1998, Canuto.C_etal_2006_book, Shen.J_2011_book] and references therein. For , this space is equivalent to the so-called Hermite space of finite smoothness, which has been attracting attention in high-dimensional computations; see [IG2015_Hermite, DILP2018, gnewuch2021countable] for its use in high-dimensional computations, and see [B1998, 1.5.4. Proposition] together with Section 2 below for the equivalence.

In this setting, we prove that the worst case integration error of Gauss–Hermite quadrature is bounded from below by up to a constant. This rate matches the upper bound shown by Mastroianni and Monegato [MM1994], and thus cannot be improved. Moreover, this rate provides a rigorous verification of the numerical findings made in [DILP2018, Section 4], where the authors computed approximate values of the worst-case error for Gauss–Hermite quadrature in the Hermite space of finite smoothness, and observed the rate for and .

It turns out that this rate is merely half of the best possible: if we allow quadrature-points and weights to be arbitrary, then the best achievable using (linear) quadratures is of the rate ; see [DILP2018, Theorem 1] for a precise statement. Dick et al. [DILP2018] also show that a class of numerical integration methods based on so-called (scaled) higher order digital nets achieve the optimal rate up to a logarithmic factor in the multi-dimensional setting, including the one-dimensional setting as a special case.

Our results on the trapezoidal rule show that the trapezoidal rule, a method arguably significantly simpler than one-dimensional higher order digital nets, also achieves this optimal rate, up to a logarithmic factor, and thus nearly twice as fast as the error-decay rate of Gauss–Hermite rule. It is also worth mentioning that Gauss–Hermite quadrature requires a nontrivial algorithm to generate quadrature points, whereas for the trapezoidal rule we simply have equispaced points.

For analytic functions, efficiency of the trapezoidal rule is well known; related studies date back at least to the paper by Goodwin in 1949 [Goodwin1949], and this accuracy is not only widely known to, but also still actively studied by contemporary numerical analysts [S1997, MS2001, GilEtAl2007, Waldvogel2011, TW2014, T2021]. Our results show that this efficiency extends to Sobolev class functions, where we do not have tools from complex analysis such as contour integrals. Our proof uses the strategy recently developed by one of the present authors and Nuyens [NS2021] for a class of quasi-Monte Carlo methods.

We now mention other error estimates of Gauss–Hermite quadrature in the literature. Based on results by Freud

[Freud1972], Smith, Sloan, and Opie [SSO1983] showed an upper bound of the rate for -times continuously differentiable functions whose -th derivative satisfies a suitable growth condition for . Since the weighted Sobolev space seems to be more frequently used, our focus is on this class. Della Vecchia and Mastroianni [DM2003] showed an upper bound of the rate for -integrable a.e. differentiable functions whose derivative is also -integrable. Moreover, their result implies a matching lower bound in the sense of worst-case error, and thus in this sense their upper bound is sharp. It does not seem to be trivial if their bounds generalise, for example to the order or to , with the -times differentiability of the integrand for . In contrast, our results show that, by assuming the square -integrability, the rate improves to , for general .

Before moving on, we mention that our results motivate further studies of Gauss–Hermite based algorithms in high-dimensional problems. Integration problems in high-dimension arise, for example in computing statistics of solutions of partial differential equations parametrised by random variables. In particular, integration with respect to the Gaussian measure has been attracting increasing attention; see for example

[Graham.I_eta_2015_Numerische, C2018, YK2019, HS2019ML, D2021, dung2022analyticity]. The measure being Gaussian, algorithms that use Gauss–Hermite quadrature as a building block have gained popularity [C2018, D2021, dung2022analyticity]. The key to proving error estimates in this context is the regularity of the quantity of interest with respect to the parameter, and for elliptic and parabolic problems such smoothness, even an analytic regularity, has been shown [BNT2010SIREV, NT2009parabolic, dung2022analyticity]. In contrast, the solutions of parametric hyperbolic systems suffer from limited regularity under mild assumptions [MNT2013_StochasticCollocationMethod, MNT2015_AnalysisComputationElasticWave]. Hence, our results, which show a suboptimality of Gauss–Hermite rule for functions with finite smoothness, caution us and encourage further studies of the use of algorithms based on Gauss–Hermite rule for this class of problems.

Finally, we note that if we have the weight function instead of , the corresponding orthogonal polynomials are called physicist’s Hermite polynomials. Our results for Gauss–Hermite quadrature can be obtained for these polynomials by simply rescaling our results by . Likewise, results for physicist’s Hermite polynomials in the literature, e.g. in [Szego1975], are used throughout this paper.

The rest of this paper is organized as follows. In Section 2 we introduce necessary definitions such as Hermite polynomials, the weighted Sobolev spaces, and the Hermite spaces. We also discuss the norm equivalence between the weighted Sobolev space and the Hermite space. In Section 3 the sub-optimality of Gauss–Hermite quadrature is shown. In particular, we obtain the matching lower and upper bounds for the worst-case error. In Section 4 the optimality of the trapezoidal rule is shown. Finally, Section 5 concludes this paper.

## 2 Function spaces with finite smoothness

Throughout this paper, we use the weighted space , the normed space consisting of the equivalence classes of Lebesgue measurable functions satisfying , where the equivalence relation is given by if and only if .

### 2.1 Hermite polynomials

For , the -th degree probabilist’s Hermite polynomial is given by

 Hk(x)=(−1)k√k!ex2/2dkdxke−x2/2,x∈R, (1)

where they are normalised so that for all . The polynomials form a complete orthonormal system for .

The following properties are used throughout the paper.

 H′k(x)=√kHk−1(x),k≥1; (2)
 dτdxτ(Hk(x)ρ(x)) =(−1)τ√(k+τ)!k!Hk+τ(x)ρ(x),k≥0,τ≥0. (3)

### 2.2 Weighted Sobolev space

The function space we consider is the Sobolev space of square integrable functions, the integrability condition of which is imposed by the Gaussian measure.

###### Definition (Weighted Sobolev space)

For , the weighted Sobolev space (with the weight function ) is the class of all functions such that has weak derivatives satisfying for :

 Hα:={f∈L2ρ∣∣∣∥f∥α:=(α∑τ=0∥f(τ)∥2L2ρ)1/2<∞}.

Elements in for are in the standard local Sobolev space , and thus admit a continuous representative. In what follows, we always take the continuous representative of .

We recall another important class of functions, the so-called Hermite space. For this space we follow the definition of [DILP2018].

###### Definition (Hermite space with finite smoothness)

For , the Hermite space with finite smoothness is given by

where , and

 rα(k):={1,ifk=0,(∑ατ=0βτ(k))−1,ifk≥1,andβτ(k):={k!(k−τ)!,ifk≥τ,0,otherwise.

It turns out , where the equality here means the norm equivalence. Hence, results established for the Hermite space can be readily translated to up to a constant, which allows us to compare the results on higher order digital nets in [DILP2018] with ours.

A proof of this equivalence of the norm is outlined in [B1998, 1.5.4. Proposition]; see also [DILP2018, Lemma 6]. Here, for completeness we give a proof.

###### Lemma 1

Let with , then with the same smoothness parameter .

###### Proof

We first prove the claim for . Assume . Let and for . We have

 ∫RF′(x)ϕ(x)dx=−∫RF(x)ϕ′(x)dx,

for any function in the space of compactly supported infinitely differentiable functions . Since is in the standard Sobolev space , there exists a sequence that satisfies

 ∥G−ϕN∥2L2(R)+∥G′−ϕ′N∥2L2(R)=∥G−ϕN∥2W1,2(R)≤1/N.

Then, letting , the Cauchy–Schwarz inequality implies , and for we also have

 ∫R|F′(x)g(x)|dx ≤(∫R|f′(x)|2e−x2/2dx)1/2(∫R|g(x)|dx)1/2 +supt∈R|(ε+1/2)2t2e−εt2|(∫R|f(x)|2e−x2/2dx)1/2(∫R|g(x)|2dx)1/2<∞.

Therefore both and define a continuous functional on . Hence, we have and thus

 ∫Rf′(x)Hk(x)ρ(x)−(ε+1/2)xHk(x)f(x)ρ(x)dx=∫RF′(x)G(x)dx = −∫RF(x)G′(x)dx=−(∫Rf(x)√kHk−1(x)ρ(x)−(−ε+1/2)xHk(x)f(x)ρ(x)dx),

which is equivalent to

 ∫Rf′(x) Hk(x)ρ(x)=−∫Rf(x)√kHk−1(x)ρ(x)−xHk(x)f(x)ρ(x)dx (4) =−∫Rf(x)(Hk(x)ρ(x))′dx=∫Rf(x)√k+1Hk+1(x)ρ(x)dx,

where we used . Hence we obtain

 (f′,Hk)2L2ρ=(k+1)(f,Hk+1)2L2ρ,

and

 ∞∑k=0(k+1)(f,Hk+1)2L2ρ=∞∑k=0(f′,Hk)2L2ρ=∥f′∥2L2ρ<∞.

This implies since . For general , assuming and by repeating the same argument as above, we have for and thus

 ∥f(τ)∥2L2ρ =∞∑k=0(f,Hk+τ)2L2ρτ∏j=1(k+j)≥1ττ∞∑k=0(f,Hk+τ)2L2ρ(k+τ)τ =1ττ∞∑k=τkτ(f,Hk)2L2ρ.

## 3 Matching bounds for Gauss–Hermite quadrature

In this section, we prove the sub-optimality of Gauss–Hermite quadrature. We first introduce the following linear quadrature of general form

 Qn(f)=n∑j=1wjf(ξj) (5)

with arbitrary distinct quadrature points on the real line

 −∞<ξ1<ξ2<⋯<ξn<∞

and quadrature weights . Gauss–Hermite quadrature is given by the points being the roots of and the weights being , ; see for example [Shen.J_2011_book, Theorem 3.5].

Given a quadrature rule , it is convenient to introduce the notation

 ewor(Qn,Hα):=sup0≠f∈Hα|I(f)−Qn(f)|∥f∥α.

The quantity is commonly referred to as the worst-case error of in ; see for example [Dick.J_Kuo_Sloan_2013_ActaNumerica]. Now we can state our aim of this section more precisely: we prove the matching lower and upper bounds on for Gauss–Hermite quadrature.

### 3.1 Lower bound

We first derive the following lower bound on for the general quadrature (5).

###### Lemma 2

Let . For , let

 σ:={αif  minj=1,…,n−1(ξj+1−ξj)≤1.0otherwise. (6)

Then, there exists a constant , which depends only on , such that the worst-case error of a general function-value based linear quadrature (5) in the weighted Sobolev space is bounded below by

 ewor(Qn,Hα)≥cαmini=1,…,n−1(ξi+1−ξi)σ+1/2 ×n−1∑j=1e−max(ξ2j,ξ2j+1)/2(n−1∑k=1e−1≥0(ξkξk+1)min(ξ2k,ξ2k+1)/2)−1/2,

where is equal to 1 if and 0 otherwise.

###### Proof

The heart of the matter is to construct a function such that for all , resulting in , and that is small but is large. Define a function by

 h(x):=hn(x):=⎧⎪⎨⎪⎩(x−ξjξj+1−ξj)α(1−x−ξjξj+1−ξj)αif there % exists j∈{1,…,n−1}such that x∈[ξj,ξj+1],0otherwise.

Then, turns out to fulfill our purpose. This type of fooling function used to prove lower bounds for the worst-case error is called bump function (of finite smoothness), in quasi-Monte Carlo theory; see, for instance, [DHP2015, Section 2.7].

First we show . It follows from

 (x−ξjξj+1−ξj)α(1−x−ξjξj+1−ξj)α=α∑ℓ=0(−1)ℓ(α)ℓ(x−ξjξj+1−ξj)α+ℓ

that we have

 h(τ)(x)=1(ξj+1−ξj)τα∑ℓ=0(−1)ℓ(αℓ)(α+ℓ)!(α+ℓ−τ)!(x−ξjξj+1−ξj)α+ℓ−τ

for and any . As we have for , the function is -times continuously differentiable. Moreover, is continuous piecewise polynomial and thus weakly differentiable. Also, noting that only when and otherwise , we have

 ∫ξj+1ξj|h(τ)(x)|2ρ(x)dx ≤e−1≥0(ξjξj+1)min(ξ2j,ξ2j+1)/2√2π∫ξj+1ξj|h(τ)(x)|2dx =e−1≥0(ξjξj+1)min(ξ2j,ξ2j+1)/2√2π(ξj+1−ξj)2τ =e−1≥0(ξjξj+1)min(ξ2j,ξ2j+1)/2√2π(ξj+1−ξj)2τ−1 ×α∑ℓ1,ℓ2=0(−1)ℓ1+ℓ22(α−τ)+ℓ1+ℓ2+1(αℓ1)(αℓ2)(α+ℓ1)!(α+ℓ1−τ)!(α+ℓ2)!(α+ℓ2−τ)!,

for . The last sum over and does not depend on . Denoting this sum by , we obtain

 ∥h∥2α =α∑τ=0∫R|h(τ)(x)|2ρ(x)dx=α∑τ=0n−1∑j=1∫ξj+1ξj|h(τ)(x)|2ρ(x)dx ≤1√2πα∑τ=0Sα,τn−1∑j=1e−1≥0(ξjξj+1)min(ξ2j,ξ2j+1)/2(ξj+1−ξj)2τ−1 ≤1√2πα∑τ=0Sα,τmin1≤i≤n−1(ξi+1−ξi)2τ−1n−1∑j=1e−1≥0(ξjξj+1)min(ξ2j,ξ2j+1)/2 ≤1√2πmin1≤i≤n−1(ξi+1−ξi)2σ−1α∑τ=0Sα,τn−1∑j=1e−1≥0(ξjξj+1)min(ξ2j,ξ2j+1)/2<∞.

This proves .

By definition of , we have for all , and thus

 Qn(h)=0.

Moreover, we have

 I(h) =∫Rh(x)ρ(x)dx=n−1∑j=1∫ξj+1ξjh(x)ρ(x)dx ≥1√2πn−1∑j=1e−max(ξ2j,ξ2j+1)/2∫ξj+1ξj(x−ξjξj+1−ξj)α(1−x−ξjξj+1−ξj)αdx =1√2πn−1∑j=1e−max(ξ2j,ξ2j+1)/2(ξj+1−ξj)∫10xα(1−x)αdx =(α!)2(2α+1)!√2πn−1∑j=1e−max(ξ2j,ξ2j+1)/2(ξj+1−ξj) ≥(α!)2(2α+1)!√2πmin1≤i≤n−1(ξi+1−ξi)n−1∑j=1e−max(ξ2j,ξ2j+1)/2.

Using the above results, we obtain

 ewor(Qn,Hα) ≥|I(h)−Qn(h)|∥h∥α ≥(α!)2(2α+1)!(2π)1/4(α∑τ=0Sα,τ)−1/2min1≤i≤n−1(ξi+1−ξi)σ+1/2 ×n−1∑j=1e−max(ξ2j,ξ2j+1)/2(n−1∑k=1e−1≥0(ξkξk+1)min(ξ2k,ξ2k+1)/2)−1/2.

Now the proof is complete.

Using the general lower bound in Lemma 2, we obtain the following lower bound on the worst-case error for Gauss–Hermite quadrature.

###### Theorem 1

Let . For any , the worst-case error of the Gauss–Hermite quadrature in the weighted Sobolev space is bounded from below as

 ewor(QGHn,Hα)≥Cαn−α/2

with a constant that depends on but independent of .

###### Proof

Let , , be the roots of . For any , it holds that

 π√n+1/2

see, for instance, [Szego1975, Eq. (6.31.22)]. Thus, to invoke Lemma 2 we let

 σ:={αfor  n≥10,0for  2≤n<10,

so that the conditions in (6) are satisfied. Also, each node is bounded below and above as follows; see, for instance, [Szego1975, Eq. (6.31.19)]: for odd, we have with the positive zeros satisfying

 jπ√n+1/2<ξGH(n+1)/2+j<4j+3√n+1/2for  j=1,…,(n−1)/2, (8)

and for even,

 (j−1/2)π√n+1/2<ξGHn/2+j<4j+1√n+1/2for  j=1,…,n/2, (9)

with symmetricity , for odd and even.

Let be odd. Using the result in Lemma 2, equations (7), and (8), together with the symmetricity of the Hermite zeros, we obtain

 ewor(QGHn,Hα) ≥2cα(π2n+1/2)σ/2+1/4(n−1)/2∑j=1e−(4j+3)2/(2n+1)⎛⎝2+2(n−1)/2∑k=2e−π2(k−1)2/(2n+1)⎞⎠−1/2.

The sum over is further bounded below by

 (n−1)/2∑j=1e−(4j+3)2/(2n+1) ≥∫(n−1)/2+11e−(4x+3)2/(2n+1)dx =√n+1/24∫(2n+5)/√n+1/27/√n+1/2e−x2/2dx ≥√n+1/24∫11√2/√7√14e−x2/2dx =√π(n+1/2)4√2(erf(11/√7)−erf(√7)),

where denotes the error function, and the last inequality holds for any odd . The sum over is further bounded above by

 (n−1)/2∑k=2e−π2(k−1)2/(2n+1) ≤∫(n−1)/21e−π2(x−1)2/(2n+1)dx ≤√n+1/2∫∞0e−π2x2/2dx=√n+1/22π.

Using these bounds, we have

 ewor(QGHn,Hα) ≥2cα(π2n+1/2)σ/2+1/4√π(n+1/2)4√2(erf(11/√7)−erf(√7)) ×⎛⎝2+2√n+1/22π⎞⎠−1/2 ≥cαπσ+1/4erf(11/√7)−erf(√7)2√2√2+√21(n+1/2)σ/2 ≥cαπ1/4erf(11/√7)−erf(√7)2(α+5)/21nα/2.

Let be even. As in the odd case, but now using (9) instead of (8), we obtain

 ewor(QGHn,Hα) ≥cα(π2n+1/2)σ/2+1/4 ×(e−(ξGHn/2+1)2/2+2n/2∑j=2e−(ξGHn/2+j)2/2)(1+2n/2∑k=2e−(ξGHn/2+k−1)2/2)−1/2 ≥cα(π2n+1/