Construction-free median quasi-Monte Carlo rules for function spaces with unspecified smoothness and general weights

We study quasi-Monte Carlo (QMC) integration of smooth functions defined over the multi-dimensional unit cube. Inspired by a recent work of Pan and Owen, we study a new construction-free median QMC rule which can exploit the smoothness and the weights of function spaces adaptively. For weighted Korobov spaces, we draw a sample of r independent generating vectors of rank-1 lattice rules, compute the integral estimate for each, and approximate the true integral by the median of these r estimates. For weighted Sobolev spaces, we use the same approach but with the rank-1 lattice rules replaced by high-order polynomial lattice rules. A major advantage over the existing approaches is that we do not need to construct good generating vectors by a computer search algorithm, while our median QMC rule achieves almost the optimal worst-case error rate for the respective function space with any smoothness and weights, with a probability that converges to 1 exponentially fast as r increases. Numerical experiments illustrate and support our theoretical findings.


page 1

page 2

page 3

page 4


A universal median quasi-Monte Carlo integration

We study quasi-Monte Carlo (QMC) integration over the multi-dimensional ...

Component-by-component construction of randomized rank-1 lattice rules achieving almost the optimal randomized error rate

We study a randomized quadrature algorithm to approximate the integral o...

Stability of lattice rules and polynomial lattice rules constructed by the component-by-component algorithm

We study quasi-Monte Carlo (QMC) methods for numerical integration of mu...

Exploration of a Cosine Expansion Lattice Scheme

In this article, we combine a lattice sequence from Quasi-Monte Carlo ru...

A note on the CBC-DBD construction of lattice rules with general positive weights

Lattice rules are among the most prominently studied quasi-Monte Carlo m...

Weighted integration over a cube based on digital nets and sequences

Quasi-Monte Carlo (QMC) methods are equal weight quadrature rules to app...

1 Introduction

We consider numerical integration of functions defined over the -dimensional unit cube . For an integrable function , we denote the integral of by

The quasi-Monte Carlo (QMC) method approximates by the equally-weighted average of function evaluations over a deterministic -element point set :

The worst-case error for a given normed function space and point set is

where is the norm of in this space. The key to success of the QMC method lies in a proper choice of the point set depending on a target class of functions. One wants to construct point sets for which this worst-case error is small and converges at the fastest possible rate as a function of , for the given space . In randomized QMC, the point set is randomized in a way that

becomes an unbiased estimator of

and one wishes to minimize its variance

[vLEC18a, vLEC02a, Lemieux09]. In this paper, we focus on deterministic QMC point sets only.

There are two main families of QMC point sets: digital nets and sequences [DickPilli10, Niederreiter92] and lattice point sets [Niederreiter92, SloanJoe94]. We refer the reader to [DKS13, Lemieux09, LeobPilli14] for further introductory details. In this paper, we consider rank-1 lattice point sets for Korobov spaces of periodic functions, and high-order polynomial lattice point sets [DP07, Niederreiter92b, Pillichshammer12] (which are a special type of digital nets) for Sobolev spaces of non-periodic functions. Each point set from these types is defined by an -dimensional generating vector, with integer coordinates in the ordinary lattice case and with polynomial coordinates in the polynomial lattice case.

In both cases, the weighted spaces of functions are defined by selecting a positive smoothness parameter and a positive weight for each subset of coordinates , with . The parameter determines how smooth the admissible functions are required to be. For the Korobov spaces, it tells the minimal rate at which the Fourier coefficients of are required to decay, and when it is an integer it corresponds to the minimal number of square-integrable mixed partial derivatives of with respect to each coordinate; see [NW_Tract1, Appendix A]. For the Sobolev spaces, is a positive integer which also imposes integrability conditions on the partial derivatives of . The weights act as constant importance factors given to the subsets of coordinates [DSWW06, SW98]. A larger means that the projection of over the subset of coordinates in can have a larger variation in some sense, so that more importance should be given to the uniformity of the points over this projection.

It is known that the best possible QMC point sets cannot provide a better convergence rate than for the worst-case error for these two function spaces. On the other hand, there are effective search algorithms which, for a given and a given selection of weights , can construct good rank-1 lattice or polynomial lattice point sets for the function spaces determined by these parameters, and for which the worst-case error converges as for any [DickPilli10, Kuo03]. Software that implement such algorithms is also freely available [LatNetbuilder, Latticebuilder, iNUY20m]. These algorithms typically use a greedy component-by-component (CBC) construction approach proposed originally by [Korobov59], then re-introduced and popularized by [SR02]. With the CBC approach, the generating vector is determined one coordinate at a time by optimizing a figure of merit that depends only on this new coordinate and the previous ones, and where the previous coordinates can no longer be changed.

In general, the number of weights to specify is . When is large, specifying all these weights individually becomes impractical, so it is common practice to “parameterize” the weights by a smaller number of parameters, usually linear instead of exponential in . The most popular forms of parameterizations are the product weights, the order-dependent weights, and their combination. For the product weights, one specifies a weight for each coordinate , and the ’s are defined as . For the order-dependent weights, depends only on the cardinality of : , where are selected positive constants. Their combination gives the product and order-dependent (POD) weights, for which [vNIC14a]. The main reason why the most popular choices of weights have this form is that the existing search algorithms are truly efficient for large only when the weights have this specific POD form [DSWW06, KSS12, LatNetbuilder, Latticebuilder, vNIC14a, NC06]

. Then, by using a fast-CBC approach that speeds up the search by exploiting a fast Fourier transform

[Latticebuilder, NC06], one can find a generating vector that gives a worst-case error of in operations for ordinary lattices and for polynomial lattices (with interlacing).

Although this form is convenient, the restriction to POD weights is limitative: for a given application, the appropriate weights may be quite far from the POD form. In this case, imposing POD implies that the point sets are constructed with the wrong weights. Moreover, even without constraints imposed on the form of weights, finding or approximating appropriate weights and the appropriate for a given application is generally very difficult [vLEC12a]. When the points are constructed with the wrong weights, the QMC method can be quite ineffective in general. These drawbacks have been addressed very partially in recent papers. In [EKNO21], the authors introduces a construction algorithm that does not require the knowledge of . It uses a CBC construction algorithm with a figure of merit that assumes , and for each coordinate it also constructs the generating vector one binary digit at a time. The method provides a convergence rate of . In [DG21], the authors study the stability of rank-1 lattice rules and polynomial lattice rules to a (limited) misspecification of and the weights, for product and POD weights. They obtain worst-case error bounds for function spaces determined by parameters and when the rules are constructed using parameters and instead, under certain conditions on those weights. These results are interesting but they do not completely eliminate the need to specify the weights.

The method studied in this paper requires no knowledge at all on and the weights . No value needs to be specified for any of these parameters. The algorithm is inspired by recent work from Pan and Owen [PanOwen2021]

, and works as follows. For a fixed odd integer

, we draw generating vectors independently and uniformly from the set of all admissible generating vectors. For each of them, we compute the corresponding QMC approximation , then we take the median of these approximations as our final estimate of . Since the method does not require the explicit construction of a good point set, we call it a construction-free median QMC rule.

Our main contribution is to prove that for representing either a weighted Korobov or weighted Sobolev space determined by parameters and , the error obeys the following type of probabilistic bound: For any and , there is a constant (which depends on , the ’s and ) such that

In other words, the worst-case error of the median estimator is bounded by a quantity that decreases almost at the best possible rate of , with a probability that converges to 1 exponentially fast as a function of . That is, we have a simple method that provides essentially the best possible convergence rate with very high probability, without requiring any knowledge of and the weights. The key reason why this is possible is that the vast majority of the choices of generating vectors turn out to be quite good and give a QMC approximation which is quite close to . Only a small minority give a large error. For the vector giving the median value to be in that small minority, there must be at least generating vectors in the sample of size that belong to this small minority, and the probability that this happens decreases towards zero exponentially in .

The remainder is organized as follows. In Section 2, we recall some basic facts on lattice rules for Korobov spaces, and we prove our main result for the median estimator in this setting. In Section 3, we do the same for high-order polynomial lattice rules in Sobolev spaces. In Section 4, we report numerical experiments to support our theoretical findings.

2 Lattice rules for Korobov spaces

2.1 Definitions

Lattice point sets are well suited for performing numerical integration of smooth periodic functions. A rank-1 lattice point set is defined as follows:

Definition 1 (rank-1 lattice point set)

Let be the number of points and . The rank-1 lattice point set defined by and the generating vector is

where denotes the fractional part of a real . The QMC algorithm using as a point set is called the rank-1 lattice rule with generating vector .

Let be periodic with an absolutely convergent Fourier series

where the dot product denotes the usual inner product of two vectors on the Euclidean space and denotes the -th Fourier coefficient of :

Note that coincides with the integral . As a class of periodic functions, we consider the following weighted Korobov space.

Definition 2 (weighted Korobov space)

Let and be a set of positive weights with . For a non-empty subset and a vector , we denote by the vector such that if and otherwise, and define

and set . The weighted Korobov space, denoted by , is a reproducing kernel Hilbert space with reproducing kernel

and inner product

We denote the induced norm by .

One wishes to have a good generating vector such that the worst-case error of the corresponding lattice rule for , defined by

is small. No good explicit construction scheme for such a is known for , so that we usually resort on a computer search algorithm as mentioned earlier. By restricting each to be in the set

we ensure that each projection of on a single coordinate contains the distinct values (no superposed points). The CBC construction algorithm for a good generating vector starts with , then for it searches for the best component from the set while keeping the earlier components unchanged.

For our median QMC rank-1 lattice rule for weighted Korobov spaces, we select an independent random sample from the set , and we approximate by the median

The worst-case error in this case is the random variable

In this random expression, we assume that are first picked randomly, then is taken as the worst-case function for the median, for these given . The index for which gives the median generally depends on .

2.2 Our main results on lattice rules for Korobov spaces

To prove our main result, we need a few more definitions.

Definition 3 (dual lattice)

For and , the set

is called the dual lattice of the rank-1 lattice point set .

The following character property of the rank-1 lattice rule is well-known, see for instance [DHP15, Lemmas 4.2 and 4.3].

Lemma 1 (character property)

For , and , we have

As our first main result, we prove a probabilistic upper bound on the worst-case error of our median rank-1 lattice rule for weighted Korobov spaces.

Theorem 1

Let be an integer, be an odd integer, and be chosen independently and randomly from the set (with replacement). Then, for any and , the worst-case error of the median rule obeys the following bound:

with a probability of at least

for any , where and denote the Euler totient function and the Riemann zeta function, respectively.


Since any has the absolutely convergent Fourier series, by applying Lemma 1 and the Cauchy–Schwarz inequality, it holds that


in which

For , by using the subadditivity


which holds for non-negative reals , see [DHP15, Theorem 2.2], and noting that the cardinality of is equal to , we have

where we write

Since it follows from [KJ02, Lemmas 2.1 & 2.2] that, for any positive integer

we obtain

This gives an upper bound on the average of over all of the admissible , which holds for any .

Then, Markov’s inequality ensures that for any , the probability of having

is at most for a random choice of . For the median estimator to be larger than this bound , we must have for at least vectors among . The probability that this happens is bounded above by

Combining this with the bound shown in (1) completes the proof.

Remark 1

One can easily prove by induction on that for . Indeed, this is true for , and for , one has

Then, for any odd , we have

Thus, the probability given in Theorem 1 must be larger than , which converges to 1 exponentially fast as a function of for .

By taking and and using the previous remark, we obtain the following corollary as a simplified version of Theorem 1.

Corollary 1

For any odd , , and , there is a constant (which depends on , the ’s and ) such that


For prime , we have and the corollary follows. For a general , we know from [RS62, Theorem 15] that

for any , where is the Euler’s constant. From this, the probabilistic bound follows.

Remark 2

Under some additional conditions on the weights , this bound depends only polynomially on the dimension , and can even be independent on the dimension . We refer to [DSWW06, Theorem 3] for the case of general weights and [DSWW06, Theorem 4] for the case of product weights.

Remark 3

As mentioned in Section 1, our median rank-1 lattice rule is motivated by the observation that most of the possible generating vectors are a good choice, but the remaining ones are bad. Although it has been already shown in the proof of Theorem 1 that most of the generating vectors yield a small value of of order , we have not shown that some of the remaining one are bad. To show this, let us first point out that coincides with the worst-case error of the rank-1 lattice rule with the given , see [DKS13, Theorem 5.12]. The averaging argument in the above proof with the case leads to an upper bound on the average of the squared worst-case error, which is only of order . Thus, a small fraction of the possible generating vectors must have a quite large of order at best.

Remark 4

It is known that rank-1 lattice rules also work for non-periodic functions by applying the tent transformation

component-wise to every point in the set [Hickernell02, DNP14, CKNS16, GSY19]. The same probabilistic upper bound, shown in Theorem 1, holds for the worst-case error of the median rule built up of the tent-transformed rank-1 lattice rules in the so-called weighted half-period cosine spaces with any parameter and weights . As shown in [DNP14, Lemma 1], the half-period cosine space coincides with an unanchored Sobolev space with smoothness 1 when .

3 High-order polynomial lattice rules for Sobolev spaces

We now consider high-order polynomial lattice point sets as defined in [DP07]. These point sets are well suited for performing numerical integration of smooth non-periodic functions. In what follows let be a prime, and be the finite field of order , which we identify with the set . Let be the set of positive integers and . For having the -adic finite expansion , we write . With these ingredients, we have the following definition from [DP07].

Definition 4 (high-order polynomial lattice point sets)

Let with , with and with . The high-order polynomial lattice point set defined by and consists of points and is given by

where is defined by

The QMC algorithm using as a point set is called the high-order polynomial lattice rule with modulus and generating vector . The order of this rule is defined as .

Typically, will be a multiple of . When , this gives the digital net construction introduced in [Niederreiter92b] and called polynomial lattice rule in [DickPilli10]. Note that [vLEC02a, vLEM03a] introduced the term “polynomial lattice rule” with a slightly different definition, in which the coordinates of the points have an infinite periodic expansion and the modulus has degree . The construction in Definition 4 essentially builds a polynomial lattice point set with points and uses only the first points.

Instead of the weighted Korobov space , we consider the following Sobolev-type Banach space as our target space for high-order polynomial lattice rules.

Definition 5 (weighted Sobolev space)

Let , , and let be a set of positive weights with . The weighted Sobolev space, denoted by , is a Banach space consisting of non-periodic (in the sense of not necessarily periodic) smooth functions with the norm

where denotes the vector such that if , if , and otherwise, and denotes the mixed derivative of order of . Moreover, we write and .

This Sobolev space was introduced by [DKLNS14]

in the context of partial differential equations with random coefficients. (The original function space in

[DKLNS14] contains the additional parameter , and the definition of the norm has been corrected in In this paper we choose , which makes the norm smallest over .) The parameter determines the differentiability of the non-periodic functions. As for the rank-1 lattice rules for the weighted Korobov spaces, it is desirable to have good modulus and generating vector such that the worst-case error of the corresponding high-order polynomial lattice rule for is small. Originally in [DKLNS14], interlaced polynomial lattice rules [G15, GD15] were used instead of high-order polynomial lattice rules, and it was shown that the worst-case error bound of order with arbitrarily small can be achieved by the CBC algorithm applied to interlaced polynomial lattice rules. The major advantage of interlaced polynomial lattice rules over high-order polynomial lattice rules lies in the construction cost for the CBC algorithm: for the product weights, constructing an interlaced rule requires operations with memory [G15], whereas constructing a high-order rule require operations with memory [BDLNP12].

However, interlaced polynomial lattice rules are not necessarily a better choice than high-order polynomial lattice rules. To construct an interlaced polynomial lattice rule, which relies on the digit interlacing method due to Dick [Dick08], we must select an integer interlacing factor , the construction cost increases linearly with , and the resulting rule cannot exploit the smoothness of functions beyond . This means that if , the worst-case error bound is only of order . High-order polynomial lattice rules do not explicitly require such an factor. We only need to specify the maximum precision of the points. This can be set as large as possible, for instance, with , we can take for the double-precision floating-point format. This way, high-order polynomial lattice rules can be possibly made adaptive to the smoothness of functions, addressing the drawback of interlaced polynomial lattice rules. Since we do not apply CBC in this paper, we prefer high-order rules over the interlaced ones.

In what follows, we assume that the polynomial is irreducible and we write

Analogously to the rank-1 lattice case in Section 2, we consider the following median high-order polynomial lattice rule for weighted Sobolev spaces. For an odd integer , we draw randomly and independently from the set , and we approximate by

The worst-case error is the random variable

3.1 Main results for polynomial lattice point sets in Sobolev spaces

We first need a few definitions and lemmas.

Definition 6 (dual polynomial lattice)

Let with , with and . The set

is called the dual net of the high-order polynomial lattice point set , where, for with the -adic finite expansion , we have defined

Definition 7 (Walsh functions)

Let us write . For , we denote the -adic expansion of by . The -th Walsh function is defined by

where the -adic expansion of is denoted by , which is understood to be unique in the sense that infinitely many of the are different from .

For and , the -dimensional -th Walsh function is defined by

It is well-known that the system of Walsh functions is a complete orthogonal system in , see [DickPilli10, Appendix A]. The following character property of the high-order polynomial lattice point set is analogous to what is stated in Lemma 1.

Lemma 2 (character property)

For with , with and , we have

For any , we have the following absolutely convergent Walsh series