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
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.
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
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
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.
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.
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.
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.
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.
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 https://arxiv.org/abs/1309.4624. 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