Construction of interlaced polynomial lattice rules for infinitely differentiable functions

# Construction of interlaced polynomial lattice rules for infinitely differentiable functions

Josef Dick, Takashi Goda, Kosuke Suzuki, Takehito Yoshiki School of Mathematics and Statistics, The University of New South Wales, Sydney, NSW 2052, Australia (josef.dick@unsw.edu.au)Graduate School of Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan (goda@frcer.t.u-tokyo.ac.jp)School of Mathematics and Statistics, The University of New South Wales, Sydney, NSW 2052, Australia (kosuke.suzuki1@unsw.edu.au)School of Mathematics and Statistics, The University of New South Wales, Sydney, NSW 2052, Australia (takehito.yoshiki1@unsw.edu.au)
July 8, 2019
###### Abstract

We study multivariate integration over the -dimensional unit cube in a weighted space of infinitely differentiable functions. It is known from a recent result by Suzuki that there exists a good quasi-Monte Carlo (QMC) rule which achieves a super-polynomial convergence of the worst-case error in this function space, and moreover, that this convergence behavior is independent of the dimension under a certain condition on the weights.

In this paper we provide a constructive approach to finding a good QMC rule achieving such a dimension-independent super-polynomial convergence of the worst-case error. Specifically, we prove that interlaced polynomial lattice rules, with an interlacing factor chosen properly depending on the number of points and the weights, can be constructed using a fast component-by-component algorithm in at most arithmetic operations to achieve a dimension-independent super-polynomial convergence. The key idea for the proof of the worst-case error bound is to use a variant of Jensen’s inequality with a purposely-designed concave function.

Keywords:  Quasi-Monte Carlo integration, super-polynomial convergence, interlaced polynomial lattice rules, infinitely differentiable functions
MSC classifications:  65C05, 65D30, 65D32

## 1 Introduction

We study the approximation of multivariate integrals of real-valued functions defined over the -dimensional unit cube ,

 I(f)=∫[0,1]sf(x)dx.

Quasi-Monte Carlo (QMC) integration approximates by using a deterministically chosen finite point set as

 I(f;P)=1|P|∑x∈Pf(x),

where denotes the cardinality of . Note that we interpret here as a set in which the multiplicity of elements matters. In order to make the integration error small for a class of functions , needs to be carefully designed depending on the class to which the function belongs. Digital nets and sequences are a well-known choice for constructing good quadrature points for several classes of functions [10, 20].

A classical criterion for measuring the distribution properties of point sets is the so-called star-discrepancy. The Koksma–Hlawka inequality bounds the integration error using a point set by the star-discrepancy of this point set times the total variation in the sense of Hardy and Krause, see for instance [17, Chapter 2, Section 5]. Thus a low-discrepancy point set of points yields a small integration error bound, typically of order with arbitrarily small , assuming that the function has bounded total variation in the sense of Hardy and Krause. Regarding explicit constructions of low-discrepancy digital nets and sequences, we refer to [10, Chapter 8] and [20, Chapter 4]. Polynomial lattice point sets, first introduced in [21], are a special construction method for digital nets and have been extensively studied in the literature, see for instance [10, Chapter 10] and [24]. Polynomial lattice rules are QMC rules using a polynomial lattice point set as quadrature points. While we usually resort to some computer search algorithm to find good polynomial lattice rules for , the major advantage of polynomial lattice rules lies in their flexibility, that is, we can design a suitable QMC rule for the problem at hand.

In order to achieve a faster convergence of the integration error, explicit constructions of point sets, referred to as higher order digital nets, have been established by Dick [2, 3] which can fully exploit the smoothness of an integrand. Specifically QMC rules using higher order digital nets achieve the optimal convergence rate of the integration error of order with arbitrarily small , when the function has square integrable partial mixed derivatives up to order in each variable. We remark that recent applications in the area of uncertainty quantification, in particular partial differential equations with random coefficients, are in need of using these types of quadrature rules, see for instance [7]. The above result by Dick is based chiefly on analyzing the decay of the Walsh coefficients of smooth functions [3, 4].

Numerical integration of infinitely many times differentiable functions in certain function spaces has recently been considered in [8, 11, 15, 16]. However, the results on higher order digital nets in [2, 3] do not improve if one assumes that the integrand is infinitely many times differentiable. More precisely, if one sets in [2, 3] one obtains constants which are infinite and the error bounds become trivial. To improve the error bounds in these papers for function spaces consisting of infinitely many times differentiable functions using higher order digital nets requires new bounds on the Walsh coefficients. Such an analysis of the Walsh coefficients was recently done in [28, 30], where they obtained a space of infinitely differentiable functions whose Walsh coefficients decay with a certain order. The worst-case error in by a digital net is closely related to the Walsh figure of merit (WAFOM) introduced in [18, 26], which is one of the computable quality criteria of digital nets, although WAFOM was originally derived in a different way from [28, 30]. Moreover, Suzuki [27] considered a weighted space of infinitely differentiable functions and studied tractability of multivariate integration in , where the positive real numbers are the weights. His result can be summarized as follows: There exists a good QMC rule using a digital net which achieves an super-polynomial convergence of the worst-case error in as , and moreover, the convergence can be independent of the dimension as for some under a certain condition on the weights .

In this paper, beyond the existence result of [27], we provide a constructive approach to finding good QMC rules achieving a dimension-independent super-polynomial convergence of the worst-case error. Specifically we prove that interlaced polynomial lattice rules can be constructed using a fast component-by-component (CBC) algorithm, in at most arithmetic operations, to achieve a dimension-independent super-polynomial convergence. As first studied in [12, 13, 14], interlaced polynomial lattice rules belong to the family of higher order digital nets and therefore achieve a higher order polynomial convergence of the integration error. We use them as QMC rules achieving a super-polynomial convergence in this paper. For this purpose, we are required to choose an interlacing factor depending on the number of points and the weights, instead of keeping it fixed (as for instance in [2, 3]). Furthermore, in order to show the worst-case error bound with a super-polynomial convergence, we purposely design a concave function to modify Jensen’s inequality which has been often used in the literature to obtain error bounds with an improved rate of convergence.

Our approach requires to set the weights for constructing a tailored QMC rule, as often encountered in this type of construction algorithms. In practical applications, however, it is not always the case where one can know in advance to which function class the functions of interest belong. To work around this drawback, it must be interesting to study whether a good convergence property which such a tailored QMC rule holds for a specific function class can be also established for other function classes, as discussed for instance in [14, Remark 1]. We observe in Section 5 that our constructed rules empirically work even for some functions not belonging to the target space. In another direction for constructing a robust QMC rule working for many different function classes, one can implement a more elaborate construction algorithm as given in [5]. However, theoretical analysis of these issues is beyond the scope of this paper and we leave them open for further research.

The remainder of this paper is organized as follows. In the next section, we introduce the necessary background and notation, namely Walsh functions, a weighted space of infinitely differentiable functions, our considering super-polynomial convergence and interlaced polynomial lattice rules. We also describe the main results of this paper. Namely, we introduce a component-by-component algorithm, state a result on the convergence behavior of interlaced polynomial lattice rules and discuss the dependence of the worst-case error bound on the dimension. In Section 3, we study the worst-case error in for QMC rules using a digital net and derive a computable upper bound. We prove in Section 4 that the CBC algorithm can be used to obtain good interlaced polynomial lattice rules which achieve a dimension-independent super-polynomial convergence of the worst-case error. Thereafter we describe the fast CBC algorithm using the fast Fourier transform as in [22, 23], and show that interlaced polynomial lattice rules achieving a dimension-independent super-polynomial convergence can be constructed in at most arithmetic operations using memory. Finally, we conclude this paper with numerical experiments in Section 5.

## 2 Background, notation and results

Throughout this paper, we use the following notation. Let be the set of positive integers and let . For a positive integer , let be a finite ring with elements, which we identify with the set equipped with addition and multiplication modulo . For , we denote its -adic expansion by with for all , which is unique in the sense that infinitely many are different from . The operators and denote digitwise addition and subtraction modulo , respectively. That is, for whose unique -adic expansions are and , and are defined as

 x⊕x′=∞∑i=1ηib−i% andx⊖x′=∞∑i=1η′ib−i,

where and , respectively. Similarly, we define digitwise addition and subtraction for non-negative integers based on their -adic expansions. In case of vectors in or , the operators and are applied componentwise.

### 2.1 Walsh functions

Walsh functions were first introduced in [29] for the case and were later generalized to arbitrary base , see for instance [1]. We refer to [10, Appendix A] for more information on Walsh functions in the context of numerical integration. We first give the definition for the one-dimensional case.

###### Definition 1.

Let be a positive integer and let be a -th root of unity. We denote the -adic expansion of by with . The -th -adic Walsh function is defined as

 bwalk(x):=ωκ0ξ1+⋯+κa−1ξab,

for with its unique -adic expansion .

This definition can be generalized to the higher-dimensional case.

###### Definition 2.

Let be a positive integer. For a dimension , let and . The -th -adic Walsh function is defined as

 bwalk(x):=s∏j=1bwalkj(xj).

Since we always use Walsh functions in a fixed base , we omit the subscript and simply write or in the remainder of this paper. From the fact that the system is a complete orthonormal system in for any [10, Theorem A.11], we have a Walsh series expansion for any

 ∑k∈Ns0^f(k)walk,

where denotes the -th Walsh coefficient of , which is defined as

 ^f(k):=∫[0,1]sf(x)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯walk(x)dx.

For continuous functions for which , the Walsh series of converges to pointwise absolutely. In fact, for any function in a weighted space which we consider in this paper, its Walsh series converges to pointwise absolutely.

### 2.2 Weighted function space Fs,u

We first define the function for a real number .

###### Definition 3.

Let be a real number. For , we denote its -adic expansion by such that and . The function is defined as

 μ(a;k):=v∑i=1(ci+a), (1)

and .

###### Remark 1.

Let us consider the case . If the sum on the right-hand side of (1) which runs over is replaced by the sum which runs over for a fixed , we recover the definitions by Niederreiter, Rosenbloom and Tsfasman in [19, 25] for and by Dick in [3] for . Our function with has been used in [18, 26]. The parameter was included in the definition originally by Yoshiki [30] for and later by Suzuki [27] for an arbitrary real number .

For the higher-dimensional case, we consider a vector of real numbers and define the function as follows.

###### Definition 4.

Let be a vector of real numbers, and let . The function is defined as

 μ(a;k):=s∑j=1μ(aj;kj).

We are now ready to introduce a weighted space of infinitely differentiable functions. Let be a sequence of positive real numbers which we call weights, and we assume that throughout this paper.

###### Definition 5.

Let be a sequence of weights. We define a weighted space as

 Fs,u:=⎧⎨⎩f∈C∞([0,1]s):∥f∥Fs,u:=sup(α1,…,αs)∈Ns0∥f(α1,…,αs)∥L1∏sj=1uαjj<∞⎫⎬⎭,

where denotes the -th mixed partial derivative of , i.e., .

In the function space , small means that higher order partial mixed derivatives associated with the -th coordinate must be relatively small. Thus, the weights play a role in moderating the importance of different variables. Owing to the refined analyses of the Walsh coefficients in [28, 30], it was shown that the Walsh coefficients of any function in decay with a certain order, as we describe in the following. Let

 mb:=minc=1,2,…,b−1|1−¯¯¯¯¯ωbc|=2sin(π/b),

and

 Mb:=maxc=1,2,…,b−1|1−¯¯¯¯¯ωbc|={2if b is even,2sin((b+1)π/2b)if b is odd.

Moreover, let

 Cb={2if b=2,Mb+bmbb−Mbif b≠2.

Then we have the following.

###### Proposition 1 ([28, 30]).

Let be a sequence of weights, and let and be constants depending only on as above. For any function in and , the -th Walsh coefficient of is bounded by

 |^f(k)|≤∥f∥Fs,ub−μ(a;k),

where is a sequence given by for all .

### 2.3 Super-polynomial convergence

From [27], it is known that there exists a good QMC rule which achieves a dimension-independent super-polynomial convergence of the worst-case error in under a certain condition on the weights . Here we briefly recall the result of [27].

The initial error in is given by the error of the zero algorithm, i.e.,

 ewor(Fs,u;∅)=supf∈Fs,u∥f∥Fs,u≤1|I(f)|,

which indeed equals 1 for any and . Hence the integration problem in is well normalized. The worst-case error in for a QMC rule using a point set is defined as

 ewor(Fs,u;P)=supf∈Fs,u∥f∥Fs,u≤1|I(f;P)−I(f)|.

We are interested in a dimension-independent super-polynomial convergence of the worst-case error of the form

 ewor(Fs,u;P)≤Ce−c(logn)pfor all n,s∈N, (2)

where and are positive constants independent of and . The following existence result is from [27].

###### Theorem 1 ([27]).

Consider the integration problem in the weighted function space for a sequence of weights . If satisfies for , then there exists a QMC rule which achieves a dimension-independent super-polynomial convergence of the worst-case error in as (2) with .

### 2.4 Interlaced polynomial lattice rules

Here we give the definition of interlaced polynomial lattice rules, which are based on polynomial lattice rules, introduced by Niederreiter [21], and a digit interlacing composition, introduced by Dick [2, 3].

We first introduce polynomial lattice rules. In this subsection, let be a prime number, and let be the finite field with elements. We denote by the set of all polynomials over , and denote by the field of formal Laurent series over . Every element of can be represented as

 L=∞∑l=wtlx−l,

for some integer and for all . For a given integer , we define the mapping from to the interval by

 vm(∞∑l=wtlx−l)=m∑l=max(1,w)tlb−l.

A non-negative integer whose -adic expansion is given by will be identified with the polynomial . For and , we define the inner product as

 k⋅q:=s∑j=1kjqj∈Zb[x], (3)

and we write if divides in . Using this notation, polynomial lattice rules are constructed as follows.

###### Definition 6.

Let . Let such that and let . A polynomial lattice point set is a set consisting of points that are defined as

 xn:=(vm(n(x)q1(x)p(x)),…,vm(n(x)qs(x)p(x)))∈[0,1)s,

for . A QMC rule using this point set is called a polynomial lattice rule with generating vector and modulus .

We add one more notation and introduce the concept of the so-called dual polynomial lattice of a polynomial lattice point set. For with its -adic expansion , let be the polynomial of degree at most obtained by truncating the associated polynomial as

 trm(k)=κ0+κ1x+⋯+κm−1xm−1,

where we set if . For a vector , we define . With this notation, we introduce the following definition of the dual polynomial lattice .

###### Definition 7.

The dual polynomial lattice of a polynomial lattice point set with modulus , , and generating vector is given by

 P⊥(q,p)={k∈Ns0: trm(k)⋅q≡0(modp)},

where the inner product is in the sense of (3).

The following important lemma relates the dual polynomial lattice to numerical integration of Walsh functions, see [10, Lemmas 4.75 and 10.6] for the proof.

###### Lemma 1.

Let be a polynomial lattice point set with modulus , , and generating vector , and let be its dual polynomial lattice. Then we have

 1bmbm−1∑n=0walk(xn)={1if k∈P⊥(q,p),0otherwise.

We introduce the digit interlacing composition next. Let be a positive integer called interlacing factor, and let be a generic point in whose unique -adic expansions are given by . Then the digit interlacing function is defined as

 Dd(x):=∞∑i=1d∑j=1ξi,jb−d(i−1)−j.

We also define such a function for -dimensional vectors by applying to every consecutive components, that is,

 Dd(x):=(Dd(x1,…,xd),Dd(xd+1,…,x2d),…,Dd(xd(s−1)+1,…,xds)).

Now we are ready to introduce the definition of interlaced polynomial lattice rules [12, 13, 14].

###### Definition 8.

Let . Let such that and let . An interlaced polynomial lattice point set of order is a set consisting of points defined as

 Dd(P(q,p)):={Dd(x):x∈P(q,p)}.

A QMC rule using this point set is called an interlaced polynomial lattice rule of order with generating vector and modulus .

### 2.5 The results

We now describe the main results of this paper. In the following, let be a prime number and let . For with and , we denote the polynomial lattice point set by with , and denote the -adic expansion of by for and . Moreover, we denote the interlaced polynomial lattice point set by , where for .

Let be a sequence of weights, and as in Proposition 1, let be the sequence given by

 aj=−logb(Cbm−1buj),j∈N.

In Section 3, we show that the worst-case error in by a QMC rule using as quadrature points is bounded by

 ewor(Fs,u;Dd(P(q,p)))≤Cu−1+CuBu(q,p),

where and are given by

 Cu=s∏j=1d∏h=1∞∏i=m+1{1+b−1bd(i−1)+h+aj},

and

 Bu(q,p)=−1+1bmbm−1∑n=0s∏j=1d∏h=1m∏i=1{1+η(ξi,n,d(j−1)+h)bd(i−1)+h+aj},

respectively, where is defined as

 η(ξ):={b−1if ξ=0,−1otherwise.

Since is independent of the modulus and generating vector , can be used as a quality criterion for searching for good and . In the following we introduce the CBC algorithm.

We restrict , , to non-zero polynomials over with its degree less than , where . Provided that is irreducible, we can set without loss of generality. We denote by the set of all non-zero polynomials over with degree less than , i.e.,

 Rb,m={q∈Zb[x]:deg(q)

We note that . Further, we write for . The idea is now to search for the polynomials component-by-component. To do so, we need to define for arbitrary . This is done in the following way. Let and . Then

 Bu(qτ,p) =−1+1bmbm−1∑n=0β−1∏j=1d∏h=1m∏i=1{1+η(ξi,n,d(j−1)+h)bd(i−1)+h+aj} ×τ−d(β−1)∏h=1m∏i=1{1+η(ξi,n,d(β−1)+h)bd(i−1)+h+aβ}. (4)

The CBC construction proceeds as follows.

###### Algorithm 1.

Let be as above.

1. Choose an irreducible polynomial with .

2. Set .

3. For , find by minimizing as a function of .

In Subsection 4.3, we show that one can also use the fast CBC algorithm of [22, 23] to find good generating vectors.

Next we show that the generating vector found by Algorithm 1 satisfies the following bound.

###### Theorem 2.

Let be a prime and be irreducible with . Let be a concave and unbounded monotonic increasing function. Suppose that is constructed using Algorithm 1. Then we have

 Bu(q,p)≤ϕ−1⎡⎢ ⎢ ⎢ ⎢⎣1bm−1∑k∈Ns0∖{0}kj

The proof of this result is presented in Subsection 4.1.

The function , , has been often used to obtain these types of error bounds in the literature. In this case, one may apply so-called Jensen’s inequality

 ϕ(∑ncn)≤∑nϕ(cn), (5)

for any sequence of non-negative real numbers . The inequality (5), however, also holds for any concave function [6, Section 2.3]. In our case, the function is not a good choice because it does not give us the worst-case error bound with a super-polynomial convergence. Instead we use which maps to for . Such a map can be designed as follows. For , let for . Then

 ϕ(x)=⎧⎨⎩2−(log2(1/x))1/λif 0

For ,

 ϕ(x)=⎧⎨⎩b−(logb(1/x))1/λif 0

Note that we set for any and , and that the function is concave and unbounded monotonic increasing on . As above we need a slight modification for the case since the function is concave over the interval but not over the interval . Using this function and under the same condition on the weights with Theorem 1, we have the following corollary of Theorem 2.

###### Corollary 1.

Assume that satisfies for . Let be a prime and be irreducible with . Suppose that is constructed using Algorithm 1. Then there exist constants both independent of such that we have

 Bu(q,p)≤Er,λb−(logb(bm−1Dr,λ))1/λ,

for any . Moreover, by setting , the worst-case error satisfies the bound

 ewor(Fs,u;Dd(P(q,p)))≤E′r,λb−(logb(bm−1Dr,λ+1))1/λ,

where is a constant independent of .

The proof of this result is presented in Subsection 4.2.

This result means that we can construct a QMC rule which achieves a dimension-independent super-polynomial convergence of the worst-case error in as (2) with . This is a bit weaker than Theorem 1 (shown by Suzuki in [27]), since we do not have an error bound for the endpoint . Under an additional assumption, however, it is even possible to include the case in Corollary 1, see Remark 2. The most important advantage of our approach is that a good QMC rule can be explicitly constructed by using a CBC algorithm.

## 3 The worst-case error in Fs,u

To analyze the worst-case error of interlaced polynomial lattice rules, we introduce a digit interlacing composition for non-negative integers. Let be an interlacing factor, and let whose -adic expansions are given by , which is actually a finite expansion. Then the digit interlacing function is defined as

 Ed(k):=∞∑i=0d∑j=1κi,jbdi+j−1.

It is obvious to show that is bijective. We also define such a function for -dimensional vectors by applying to every consecutive components, that is,

 Ed(k):=(Ed(k1,…,kd),Ed(kd+1,…,k2d),…,Ed(kd(s−1)+1,…,kds)).

The following lemma relates an interlaced polynomial lattice point set to numerical integration of Walsh functions, see [12, Lemma 1] for the proof.

###### Lemma 2.

Let be an interlaced polynomial lattice point set of order with modulus , , and generating vector . For , we have

 1bmbm−1∑n=0walEd(k)(yn)={1if k∈P⊥(q,p),0otherwise.

We introduce another function for a real number and an integer . For , we denote its -adic expansion by such that and . The function is defined as

 ~μ(a,h;k):=v∑i=1[d(ci−1)+h+a],

and . For vectors of real numbers and , we define

 ~μ(a;k):=s∑j=1d∑h=1~μ(aj,h;kd(j−1)+h).

With a slight abuse of notation, for and , we write , where the vector denotes the -dimensional vector whose -th component is for and otherwise. From Definition 4 and the definition of , we have

 μ(a;Ed(k))=~μ(a;k). (8)

Now the worst-case error for numerical integration in using an interlaced polynomial lattice rule is given as follows.

###### Proposition 2.

Let be an interlaced polynomial lattice point set of order with modulus , , and generating vector . For a sequence of the weights , we have

 ewor(Fs,u;Dd(P(q,p)))≤∑k∈P⊥(q,p)∖{0}b−~μ(a;k),

where is the dual polynomial lattice of , and is a sequence of real numbers given as in Proposition 1.

###### Proof.

We write . Let us consider a function . Given the Walsh series expansion of and the fact that is bijective, the signed integration error becomes

 I(f;Dd(P(q,p)))−I(f) =1bmbm−1∑n=0f(yn)−^f(0) =1