A reduced fast construction of polynomial lattice point sets with low weighted star discrepancy

TL;DR

This paper introduces a faster method for constructing polynomial lattice point sets with low weighted star discrepancy by limiting the polynomial set size, leveraging rapidly decaying weights to reduce computational costs.

Contribution

It proposes a reduced construction algorithm for polynomial lattice point sets that maintains low discrepancy while significantly lowering computational effort.

Findings

Construction cost is substantially reduced.

Maintains low weighted star discrepancy.

Applicable to high-dimensional integration problems.

Abstract

The weighted star discrepancy is a quantitative measure for the performance of point sets in quasi-Monte Carlo algorithms for numerical integration. We consider polynomial lattice point sets, whose generating vectors can be obtained by a component-by-component construction to ensure a small weighted star discre-pancy. Our aim is to significantly reduce the construction cost of such generating vectors by restricting the size of the set of polynomials from which we select the components of the vectors. To gain this reduction we exploit the fact that the weights of the spaces we consider decay very fast.

Equations

$$I_s(F):={\int }_{[0,1{)}^s}F(\mathbf{x})\mathrm{d}\mathbf{x}$$

$$Q_{N,s}(F):=\frac{1}{N}\sum _{n=0}^{N-1}F({\mathbf{x}}_n).$$

$$\Delta (\mathbf{t},\mathcal{P}):=\frac{1}{N}\sum _{n=0}^{N-1}{1}_{[0,\mathbf{t})}({\mathbf{x}}_n)-\prod _{j=1}^s{t}_j,$$

$$D_{N,\mathbf{\gamma }}^{\ast }(\mathcal{P}):=\underset{\mathbf{t}\in (0,1{]}^s}{\sup }\underset{\varnothing \ne \mathfrak{u}\subseteq [s]}{\max }{\gamma }_{\mathfrak{u}}|\Delta (({\mathbf{t}}_{\mathfrak{u}},1),\mathcal{P})|,$$

$$|Q_{N,s}(F)-I_s(F)|\le D_{N,\mathbf{\gamma }}^{\ast }(\mathcal{P})\Vert F{\Vert }_{\mathbf{\gamma }},$$

$$L=\sum _{l=\omega }^{\infty }{t}_l{x}^{-l},$$

$${\varphi }_m\left(\sum _{l=\omega }^{\infty }{t}_l{x}^{-l}\right)=\sum _{l=\max \{1,\omega \}}^m{t}_l{p}^{-l}.$$

$$n(x)=\sum _{r=0}^{m-1}{n}_r{x}^r\in {\mathbb{F}}_p[x],$$

$${\mathbf{x}}_n=\left({\varphi }_m\left(\frac{n(x)g_1(x)}{f(x)}\right),\dots ,{\varphi }_m\left(\frac{n(x)g_s(x)}{f(x)}\right)\right)\in [0,1{)}^s$$

$$G_{p,m}(f):=\{g\in G_{p,m}\mid \gcd (g,f)=1\}.$$

$${\mathcal{G}}_{p,m-w}(f):=\{\begin{array}{ll}{G}_{p,m-w}(f)& \text{mbox}ifw<m,\\ \{1\in {\mathbb{F}}_p[x]\}& \text{mbox}ifw\ge m\end{array}$$

$${\mathrm{disc}}_{\infty }(N,s):=\underset{\begin{array}{c}\mathcal{P}\subseteq [0,1{)}^s\\ \#\mathcal{P}=N\end{array}}{\inf }{D}_{N,\mathbf{\gamma }}^{\ast }(\mathcal{P}),$$

$$N^{\ast }(s,\epsilon ):=\min \{N\in \mathbb{N}\mid {\mathrm{disc}}_{\infty }(N,s)\le \epsilon \}.$$

$$N^{\ast }(s,\epsilon )\le C{\epsilon }^{-\tau }$$

$$\sum _{j=1}^{\infty }{\gamma }_j{p}^{w_j}<\infty$$

$$\mathcal{O}\left(N+\min \{s,t\}N+N\sum _{d=1}^{\min \{s,t\}}(m-w_d)p^{-w_d}\right)$$

$$\sum _{j=1}^{\infty }{\gamma }_j{p}^{w_j}\le \sum _{j=1}^{\infty }{j}^{-k}{j}^{k-\alpha }=\sum _{j=1}^{\infty }{j}^{-\alpha }=\zeta (\alpha )<\infty ,$$

$$D_{N,\mathbf{\gamma }}^{\ast }(\mathbf{g},f)\le \sum _{\begin{array}{c}\mathfrak{u}\subseteq [s]\\ \mathfrak{u}\ne \varnothing \end{array}}{\gamma }_{\mathfrak{u}}\left(1-{\left(1-\frac{1}{N}\right)}^{|\mathfrak{u}|}\right)+R_{\mathbf{\gamma }}^s(\mathbf{g},f),$$

$$R_{\mathbf{\gamma }}^s(\mathbf{g},f)=\sum _{\begin{array}{c}\mathbf{h}\in G_{p,m}^s\setminus \{0\}\\ \mathbf{h}\cdot \mathbf{g}\equiv 0\mathrm{mod}f\end{array}}\prod _{i=1}^s{r}_p(h_i,{\gamma }_i).$$

$$r_p(h,\gamma )=\{\begin{array}{ll}1+\gamma & \text{mbox}ifh=0,\\ \gamma r_p(h)& \text{mbox}otherwise,\end{array}$$

$$r_p(h)=\frac{1}{p^{a+1}{\sin }^2\left(\frac{\pi }{p}{h}_a\right)}.$$

$$\sum _{\begin{array}{c}\mathfrak{u}\subseteq [s]\\ \mathfrak{u}\ne \varnothing \end{array}}{\gamma }_{\mathfrak{u}}\left(1-{\left(1-\frac{1}{N}\right)}^{|\mathfrak{u}|}\right)\le \frac{\max (1,\Gamma ){\mathrm{e}}^{{\sum }_{i=1}^{\infty }{\gamma }_i}}{N},$$

$$R_{\mathbf{\gamma }}^{d+1}(x^{w_1}{g}_1,\dots ,x^{w_d}{g}_d,x^{w_{d+1}}{g}_{d+1})$$

$$R_{\mathbf{\gamma }}^d((x^{w_1}{g}_1,\dots ,x^{w_d}{g}_d),x^m)\le \frac{1}{p^m}\prod _{i=1}^d\left(1+{\gamma }_i+{\gamma }_{i}2p^{\min \{w_i,m\}}m\frac{p^2-1}{3p}\right).$$

$$D_{N,\mathbf{\gamma }}^{\ast }$$

$$\le \sum _{\begin{array}{c}\mathfrak{u}\subseteq [s]\\ \mathfrak{u}\ne \varnothing \end{array}}{\gamma }_{\mathfrak{u}}\left(1-{\left(1-\frac{1}{N}\right)}^{|\mathfrak{u}|}\right)+\frac{1}{N}\prod _{i=1}^s\left(1+{\gamma }_i+{\gamma }_{i}2p^{\min \{w_i,m\}}m\frac{p^2-1}{3p}\right).$$

$$\sum _{j=1}^{\infty }{\gamma }_j{p}^{w_j}<\infty$$

$$\sum _{\begin{array}{c}h\in G_{p,m}\setminus \{0\}\\ a\mid h\end{array}}{r}_p(h)=\left(m-\mathrm{deg}(a)\right)\frac{p^2-1}{3p}{p}^{-\mathrm{deg}(a)}.$$

$$\sum _{h\in G_{p,m}\setminus \{0\}}{r}_p(h)=m\frac{p^2-1}{3p}.$$

$$\sum _{v\in G_{p,m}}{X}_p\left(\frac{v}{f}g\right)=\{\begin{array}{ll}{p}^m& \text{mbox}iff\mid g,\\ 0& \text{mbox}otherwise.\end{array}$$

$$R_{\mathbf{\gamma }}^s(\mathbf{g},f)=-\prod _{i=1}^s(1+{\gamma }_i)+\frac{1}{p^m}\sum _{v\in G_{p,m}}\prod _{i=1}^s\left(1+{\gamma }_i+{\gamma }_i\sum _{h\in G_{p,m}\setminus \{0\}}{r}_p(h)X_p\left(\frac{v}{f}h{g}_i\right)\right).$$

Full text

Keywords: weighted star discrepancy, polynomial lattice point sets, quasi-Monte Carlo integration, component-by-component algorithm

MSC 2000: 11K06, 11K38, 65D30, 65D32

1 Introduction

A convenient way to approximate the value of an integral

[TABLE]

over the $s$-dimensional unit cube is to use a quasi-Monte Carlo rule of the form

[TABLE]

The integrand $F$ usually stems from some suitable (weighted) function space and the multiset $\mathcal{P}$ of integration nodes $\boldsymbol{x}_{0},\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N-1}$ in the algorithm $Q_{N,s}(F)$ is chosen deterministically from $[0,1)^{s}$. For comprehensive information on quasi-Monte Carlo algorithms consult, e.g., [5, 3, 9, 12]. The quality of a quasi-Monte Carlo rule is for instance measured by some notion of discrepancy. In this paper we consider the weighted star discrepancy, which has been introduced by Sloan and Woźniakowski in [21], exploiting the insight that the weights reflect the influence of different coordinates on the integration error. Let $[s]:=\{1,2,\dots,s\}$ and consider a weight sequence $\boldsymbol{\gamma}=(\gamma_{\mathfrak{u}})_{\mathfrak{u}\subseteq[s]}$ of nonnegative real numbers, i.e., every group of variables $(x_{i})_{i\in\mathfrak{u}}$ is equipped with a weight $\gamma_{\mathfrak{u}}$. Roughly speaking, a small weight indicates that the corresponding variables contribute little to the integration problem. For simplicity, throughout this paper we only consider product weights, defined as follows. Given a non-increasing sequence of positive real numbers $(\gamma_{j})_{j\geq 1}$ with $\gamma_{j}\leq 1$ we set $\gamma_{\mathfrak{u}}:=\prod_{j\in\mathfrak{u}}\gamma_{j}$ and $\gamma_{\emptyset}:=1$.

Definition 1

Let $\boldsymbol{\gamma}=(\gamma_{\mathfrak{u}})_{\mathfrak{u}\subseteq[s]}$ be a weight sequence and $\mathcal{P}=\left\{\boldsymbol{x}_{0},\dotsc,\boldsymbol{x}_{N-1}\right\}\subseteq[0,1)^{s}$ be an $N$-element point set. The local discrepancy of the point set $\mathcal{P}$ at $\boldsymbol{t}=(t_{1},\dotsc,t_{s})\in(0,1]^{s}$ is defined as

[TABLE]

where $\boldsymbol{1}_{[\boldsymbol{0},\boldsymbol{t})}$ denotes the characteristic function of $[\boldsymbol{0},\boldsymbol{t}):=[0,t_{1})\times\cdots\times[0,t_{s})$. The weighted star discrepancy of $\mathcal{P}$ is then defined as

[TABLE]

where $(\boldsymbol{t}_{\mathfrak{u}},\boldsymbol{1})$ denotes the vector $(\tilde{t}_{1},\dotsc,\tilde{t}_{s})$ with $\tilde{t}_{j}=t_{j}$ if $j\in\mathfrak{u}$ and $\tilde{t}_{j}=1$ if $j\notin\mathfrak{u}$.

A relation between the integration error of quasi-Monte Carlo rules and the weighted star discrepancy is given by the Koksma-Hlawka type inequality (see [21])

[TABLE]

where $\|\cdot\|_{\boldsymbol{\gamma}}$ is some norm which depends only on the weight sequence $\boldsymbol{\gamma}$ but not on the point set $\mathcal{P}$.

It turns out that lattice point sets (see, e.g., [12, Chapter 5], [9]) and polynomial lattice point sets (see, e.g., [12, Chapter 4], [13], [5, Chapter 10]) are often a good choice as sample points in (1). These two kind of point sets are strongly connected and have a lot of parallel tracks in their analysis. However, there are some situations were one type of point set is superior to the other in terms of error bounds or the size of the function classes where they yield good results for numerical integration. Thus it is beneficial to have constructions for lattice point sets as well as for polynomial lattice point sets at hand. For a detailed comparison of lattice point sets and polynomial lattice point sets see, e.g., [19]. Also, Ch. Schwab, in response to the first author’s talk about constructing lattice point sets at the MCQMC 2016 conference in Stanford, pointed out that it would be an interesting problem to extend the result in [7] to polynomial lattice point sets. Thus, in this paper we study polynomial lattice point sets, a special class of point sets with low weighted star discrepancy, introduced by Niederreiter in [12, Chapter 4], [13]. For a prime number $p$, let $\mathbb{F}_{p}$ be the finite field of order $p$. We identify $\mathbb{F}_{p}$ with the set $\{0,1,\ldots,p-1\}$ equipped with the modulo $p$ arithmetic. We denote by $\mathbb{F}_{p}[x]$ the set of polynomials over $\mathbb{F}_{p}$ and by $\mathbb{F}_{p}((x^{-1}))$ the field of formal Laurent series over $\mathbb{F}_{p}$ with elements of the form

[TABLE]

where $\omega\in\mathbb{Z}$ and $t_{l}\in\mathbb{F}_{p}$ for all $l\geq\omega$. For a given dimension $s\geq 2$ and some integer $m\geq 1$ we choose a so-called modulus $f\in\mathbb{F}_{p}[x]$ with $\deg(f)=m$ as well as polynomials $g_{1},\dots,g_{s}\in\mathbb{F}_{p}[x]$. The vector $\boldsymbol{g}=(g_{1},\dots,g_{s})$ is called the generating vector of the polynomial lattice point set. Further, we introduce the map $\phi_{m}:\mathbb{F}_{p}((x^{-1}))\to[0,1)$ such that

[TABLE]

With $n\in\{0,1,\dots,p^{m}-1\}$ we associate the polynomial

[TABLE]

as each such $n$ can uniquely be written as $n=n_{0}+n_{1}p+\dots+n_{m-1}p^{m-1}$ with digits $n_{r}\in\{0,1,\dots,p-1\}$ for all $r\in\{0,1,\dots,m-1\}$. With this notation, the polynomial lattice point set $\mathcal{P}(\boldsymbol{g},f)$ is defined as the set of $N:=p^{m}$ points

[TABLE]

for $0\leq n\leq p^{m}-1$. See also [5, Chapter 10].

In the following, by $G_{p,m}$ we denote the set of all polynomials $g$ over $\mathbb{F}_{p}$ with $\deg(g)<m$. Further we define

[TABLE]

For the weighted star discrepancy of a polynomial lattice point set we simply write $D_{N,\boldsymbol{\gamma}}^{\ast}(\boldsymbol{g},f)$.

Niederreiter [12] proved the existence of polynomial lattice point sets with low unweighted star discrepancy by averaging arguments. Generating vectors of good polynomial lattice point sets can be constructed by a component-by-component (CBC) construction. The standard structure of CBC constructions is as follows. We start by setting the first coordinate of the generating vector equal to 1. After this first step we proceed by increasing the dimension of the generating vector by one in each step until we have a generating vector $(g_{1},\dots,g_{s})$ of full size $s$. That is, all previously chosen components stay the same and one new component is added. This new coordinate is chosen from a given search set, most commonly from $G_{p,m}(f)$ given by (2). Usually it is determined such that the weighted star discrepancy of the lattice point set, corresponding to the generating vector, consisting of all previously chosen components plus one additional component, is minimized as a function of this last component.

Such constructions were provided in [4] for an irreducible modulus $f$ and in [1] for a reducible $f$. In these papers, the authors considered the unweighted star discrepancy as well as its weighted version, which we study here. It is the aim of the present paper to speed up these constructions by reducing the search sets for the components of the generating vector $\boldsymbol{g}$ according to each component’s importance. It is the nature of product weighted spaces that the components $g_{j}$ of the generating vector have less and less influence on the quality of the corresponding polynomial lattice point as $j$ increases. Roughly speaking this is due to the weights $(\gamma_{j})$ that are becoming ever smaller with increasing index $j$. We want to exploit this property in the following way. As the components’ influence is decreasing with their indices we want to use less and less time and computational cost to choose these components. To achieve this we choose them from smaller and smaller search sets, which are defined as follows. Let $w_{1}\leq w_{2}\leq\cdots$ be a non-decreasing sequence of nonnegative integers. This sequence of $w_{j}$’s is determined in accordance with the weight sequence $\boldsymbol{\gamma}$. Loosely speaking, the smaller $\gamma_{j}$, the bigger $w_{j}$ is chosen. For $w\in\mathbb{N}_{0}$ with $w<m$ we define $G_{p,m-w}$ and $G_{p,m-w}(f)$ analogously to $G_{p,m}$ and $G_{p,m}(f)$, respectively. Further we set

[TABLE]

for any $w\in\mathbb{N}_{0}$. For $w<m$ these sets have cardinality $p^{m-w}-1$ in the case of an irreducible modulus $f$ and $p^{m-w-1}(p-1)$ for the special case $f:\mathbb{F}_{p}\to\mathbb{F}_{p},x\mapsto x^{m}$. We will consider these two cases in what follows. Finally, for $d\in[s]$, we define $\mathcal{G}_{p,m-\boldsymbol{w}}^{d}(f):=\mathcal{G}_{p,m-w_{1}}(f)\times\dots\times\mathcal{G}_{p,m-w_{d}}(f)$. The idea is to choose the $i$th component of $\boldsymbol{g}$ of the form $x^{w_{i}}g_{i}$, where $g_{i}\in\mathcal{G}_{p,m-w_{i}}(f)$, i.e., the search set for the $i$th component is reduced by a factor $p^{-\min\{w_{i},m\}}$ in comparison to the standard CBC construction. We will show that under certain conditions on the weights $\boldsymbol{\gamma}$ and the parameters $w_{i}$ a polynomial lattice point set constructed according to our reduced CBC construction has a low weighted star discrepancy of order $N^{-1+\delta}$ for all $\delta>0$. The standard CBC construction (cf. [20]) can be done in $\mathcal{O}(sN^{2})$ operations. To speed up the construction, in a first step, making use of ideas from Nuyens and Cools [17, 18] on fast Fourier transformation, the construction cost can be reduced to $\mathcal{O}(sN\log{N})$, as for example done in [4]. Combining this with our reduced search sets we obtain a computational cost that is independent of the dimension eventually. Reduced CBC constructions have been introduced first by Dick et al. in [2] for lattice and polynomial lattice point sets with a small worst case integration error in Korobov and Walsh spaces, respectively, and have also been investigated in [7] for lattice point sets with small weighted star discrepancy.

An interesting aspect of the discrepancy of high dimensional point sets is the so-called tractability of discrepancy (see, e.g., [14, 15, 16] for detailed information). For $N,s\in\mathbb{N}$ let

[TABLE]

the $N$th minimal star discrepancy. To introduce the concept of tractability of discrepancy we define the information complexity (also called the inverse of the weighted star discrepancy) as

[TABLE]

Thus $N^{\ast}(s,\varepsilon)$ is the minimal number of points required to achieve a weighted star discrepancy of at most $\varepsilon$. To keep the construction cost of our generating vector low, it is, of course, beneficial to have a small information complexity and thus to stand a chance to have a polynomial lattice point set of small size. This is why we are interested in how fast the information complexity grows when $s$ and $\varepsilon^{-1}$ tend to infinity. Tractability describes this dependence of the information complexity on the dimension $s$ and the error demand $\varepsilon$. The best we can hope for is the case where $N^{\ast}(s,\varepsilon)$ is independent of $s$ and depends at most polynomially on $\varepsilon^{-1}$. To be more precise, we say that we achieve strong polynomial tractability if there exist constants $C,\tau>0$ such that

[TABLE]

for all $s\in\mathbb{N}$ and all $\varepsilon\in(0,1)$. Roughly speaking, a problem is considered tractable if its information complexity’s dependence on $s$ and $\varepsilon^{-1}$ is not exponential. Taking weights into account in the definition of discrepancy can sometimes overcome the so-called curse of dimensionality, i.e., an exponential dependence of $N^{\ast}(s,\varepsilon)$ on $s$. We will show that our reduced fast CBC algorithm finds a generating vector $\boldsymbol{g}$ of a polynomial lattice point set that achieves strong polynomial tractability provided that

[TABLE]

with a construction cost of

[TABLE]

operations, where $t=\max\{j\in\mathbb{N}\mid w_{j}<m\}$.

Before stating our main results we would like to discuss a motivating example. Consider first the standard CBC construction as treated in [1, 4], where $w_{j}=0$ for all $j\geq 0$. In this case, a sufficient condition for strong polynomial tractability is $\sum_{j=1}^{\infty}\gamma_{j}<\infty$, which for instance is satisfied for the special choices $\gamma_{j}=j^{-2}$ and $\gamma_{j}=j^{-1000}$. However, in the second example the weights decay much faster than in the first. We can make use of this fact by introducing the sequence $\boldsymbol{w}=(w_{j})_{j\geq 0}$ such that the condition $\sum_{j=1}^{\infty}\gamma_{j}p^{w_{j}}<\infty$ holds, while still achieving strong polynomial tractability (see Corollary 2). This way, we can reduce the size of the search sets for the components of the generating vector if the weights $\gamma_{j}$ decay very fast. Consider for example the weight sequence $\gamma_{j}=j^{-k}$ for some $k>1$. For $w_{j}=\lfloor(k-\alpha)\log_{p}{j}\rfloor$ with arbitrary $1<\alpha<k$ we find

[TABLE]

where $\zeta$ denotes the Riemann Zeta function. Observe that for large $k$, i.e., fast decaying weights, we may choose smaller search sets and thereby speed up the CBC algorithm. This paper is organized as follows. In the next section we give an algorithm for constructing polynomial lattice point sets and we derive an upper bound on the weighted star discrepancy of the point set constructed with this algorithm. We also give tractability results and analyze the computational cost of our algorithm. At first, we consider the case where $f:\mathbb{F}_{p}\to\mathbb{F}_{p},x\mapsto x^{m}$. Then we consider the case where the modulus of the polynomial lattice point set is irreducible.

2 A reduced CBC construction

In this section we present a CBC construction for the vector $\left(x^{w_{1}}g_{1},\dotsc,x^{w_{s}}g_{s}\right)$ and an upper bound for the weighted star discrepancy of the corresponding polynomial lattice point set.

First note that if $\boldsymbol{g}\in G_{p,m}^{s}$, then it is known (see [4]) that

[TABLE]

where in the case of product weights we have

[TABLE]

Here, for elements $\boldsymbol{h}=(h_{1},\dots,h_{s})$ and $\boldsymbol{g}=(g_{1},\dots,g_{s})$ in $G_{p,m}^{s}$ we define the scalar product by $\boldsymbol{h}\cdot\boldsymbol{g}:=h_{1}g_{1}+\dots+h_{s}g_{s}$. The numbers $r_{p}(h,\gamma)$ for $h\in G_{p,m}$ and $\gamma\in\mathbb{R}$ are defined as

[TABLE]

where for $h=h_{0}+h_{1}x+\dots+h_{a}x^{a}$ with $h_{a}\neq 0$ we set

[TABLE]

Thus, in order to analyze the weighted star discrepancy of a polynomial lattice point set it suffices to investigate the quantity $R^{s}_{\boldsymbol{\gamma}}(\boldsymbol{g},f)$. This is due to the result of Joe [10], who proved that for any summable weight sequence $(\gamma_{j})_{j\geq 1}$ we have

[TABLE]

with $\Gamma:=\sum_{i=1}^{\infty}\frac{\gamma_{i}}{1+\gamma_{i}}$.

Algorithm 1

Let $p\in\mathbb{P}$, $m\in\mathbb{N}$, $f\in\mathbb{F}_{p}[x]$ and let $(w_{j})_{j\geq 1}$ be a non-decreasing sequence of nonnegative integers and consider product weights $(\gamma_{j})_{j\geq 1}$. Construct $\left(g_{1},\dotsc,g_{s}\right)\in\mathcal{G}_{p,m-\boldsymbol{w}}^{s}(f)$ as follows:

Set $g_{1}=1$. 2. 2.

For $d\in[s-1]$ assume $(g_{1},\dotsc,g_{d})\in\mathcal{G}_{p,m-\boldsymbol{w}}^{d}(f)$ to be already found. Choose $g_{d+1}\in\mathcal{G}_{p,m-w_{d+1}}(f)$ such that

[TABLE]

is minimized as a function of $g_{d+1}$. 3. 3.

Increase $d$ by 1 and repeat the second step until $\left(g_{1},\dotsc,g_{s}\right)$ is found.

Remark 1

Of course we have $\mathcal{G}_{p,m-\boldsymbol{w}}^{s}(f)\subseteq G_{p,m}^{s}$, and thus in Algorithm 1 it indeed suffices to consider $R_{\boldsymbol{\gamma}}^{d+1}$ rather than the weighted star discrepancy.

In the algorithm above, the search set is reduced for each coordinate of $\left(g_{1},\dotsc,g_{s}\right)$ according to its importance, as with increasing $w_{j}$ the search set becomes smaller, as the weight $\gamma_{j}$ and thus the corresponding component’s influence on the quality of the generating vector decreases. For this reason we call Algorithm 1 a reduced CBC algorithm. We will now study Algorithm 1 for different choices of $f$.

2.1 Polynomial lattice point sets for $f(x)=x^{m}$

We will now study the interesting case where $f\colon\mathbb{F}_{p}\to\mathbb{F}_{p},x\mapsto x^{m}$. Throughout the rest of this section we write $x^{m}$ instead of $f$ to emphasize our special choice of $f$. Note that for $g\in\mathbb{F}_{p}((x^{-1}))$ the Laurent series $g/f$ can be easily computed in this case by shifting the coefficients of $g$ $m$ times to the left. It is the aim of this section to prove the following theorem:

Theorem 1

Let $\boldsymbol{\gamma}=(\gamma_{j})_{j\geq 1}$ and $\boldsymbol{w}$ with $0=w_{1}\leq w_{2}\leq\cdots$. Let further $(g_{1},\dots,g_{s})\in\mathcal{G}_{p,m-\boldsymbol{w}}^{s}(x^{m})$ be constructed using Algorithm 1. Then we have for every $d\in[s]$

[TABLE]

As a direct consequence we obtain the following discrepancy estimate.

Corollary 1

Let $N=p^{m}$ and $\boldsymbol{\gamma}$, $\boldsymbol{w}$ and $(g_{1},\dots,g_{s})$ as in Theorem 1. Then the polynomial lattice point set $\mathcal{P}\left((x^{w_{1}}g_{1},\dotsc,x^{w_{s}}g_{s}),x^{m}\right)$ has a weighted star discrepancy

[TABLE]

Knowing the above discrepancy bound, we are now ready to ask about the size of the polynomial lattice point set required to achieve a weighted star discrepancy not exceeding some $\varepsilon$ threshold. In particular, we would like to know how this size depends on the dimension $s$ and on $\varepsilon$.

Corollary 2

Let $N=p^{m}$, $\boldsymbol{\gamma}$ and $\boldsymbol{w}$ as in Theorem 1 and consider the problem of constructing generating vectors for polynomial lattice point sets with small weighted star discrepancy. Then

[TABLE]

is a sufficient condition for strong polynomial tractability. This condition further implies $D_{N,\boldsymbol{\gamma}}^{*}\left((x^{w_{1}}g_{1},\dotsc,x^{w_{s}}g_{s}),x^{m}\right)=\mathcal{O}(N^{-1+\delta})$, with the implied constant independent of $s$, for any $\delta>0$, where $(g_{1},\dots,g_{s})\in\mathcal{G}_{p,m-\boldsymbol{w}}^{s}(x^{m})$ is constructed using Algorithm 1.

• Proof.

Construct a generating vector $(g_{1},\dots,g_{s})\in\mathcal{G}_{p,m-w}^{s}(x^{m})$ applying Algorithm 1 and consider its weighted star discrepancy, bounded by (1). Following closely the lines of the argumentation in [7, Section 5] and noticing that $2m\frac{p^{2}-1}{3p}=\mathcal{O}(\log N)$ we obtain the result. More precisely, provided that the $\gamma_{j}p^{w_{j}}$’s are summable, we have a means to construct polynomial lattice point sets $\mathcal{P}(\boldsymbol{g},f)$ with $D^{*}_{N,\boldsymbol{\gamma}}(\boldsymbol{g},f)\leq\varepsilon$, whose sizes grow polynomially in $\varepsilon^{-1}$ and are independent of the dimension. As a result the problem is strongly polynomially tractable. The discrepancy result $D_{N,\boldsymbol{\gamma}}^{*}\left((x^{w_{1}}g_{1},\dotsc,x^{w_{s}}g_{s}),x^{m}\right)=\mathcal{O}(N^{-1+\delta})$ also follows directly from [7]. $\hfill\Box$

Remark 2

Recall that $t=\max\{j\in\mathbb{N}\colon w_{j}<m\}$ and note that setting $w_{j}=m$ for all $j>t$ does neither change the bound on the weighted star discrepancy nor the computational cost of Algorithm 1. It might change the generating vector though. If so, however, only components with very little influence on the quality of the point set are altered. Defining $w_{j}=m$ for all $j>t$, it suffices to have a summable weight sequence $\boldsymbol{\gamma}$ in order to achieve strong polynomial tractability, as long as $t$ is finite.

In order to show Theorem 1 we need several auxiliary results.

Lemma 1

Let $a\in\mathbb{F}_{p}[x]$ be monic. Then we have

[TABLE]

In particular, for $a=1$ this formula yields

[TABLE]

• Proof.

This fact follows from [1, p. 1055] (by setting $\gamma_{d+1}=1$). The special case $a=1$ also follows from [4, Lemma 2.2] by setting $s=1$. $\hfill\Box$

For our purposes, it is convenient to write $R^{s}_{\boldsymbol{\gamma}}(\boldsymbol{g},f)$ from (4) in an alternative way. To this end, we introduce some notation. For a Laurent series $L\in\mathbb{F}_{p}((x^{-1}))$ we denote by $c_{-1}(L)$ its coefficient of $x^{-1}$, i.e., its residuum. Further, we set $X_{p}(L):=\chi_{p}(c_{-1}(L))$, where $\chi_{p}$ is a non-trivial additive character of $\mathbb{F}_{p}$. One could for instance choose $\chi_{p}(n)=\mathrm{e}^{\frac{2\pi\mathrm{i}}{p}n}$ for $n\in\mathbb{F}_{p}$ (see, e.g., [11]). It is clear that $X_{p}(L)=1$ if $L$ is a polynomial and that $X_{p}(L_{1}+L_{2})=X_{p}(L_{1})X_{p}(L_{2})$ for $L_{1},L_{2}\in\mathbb{F}_{p}((x^{-1}))$. From [12, p. 78] we know that

[TABLE]

Lemma 2

We have

[TABLE]

• Proof.

We employ the properties of $X_{p}$ as stated above to obtain from (4)

[TABLE]

and the claimed formula is verified. $\hfill\Box$

Now we study a sum which will appear later in the proof of Theorem 1 and show an upper bound for it.

Lemma 3

Let $w\in\mathbb{N}_{0}$ and $v\in G_{p,m}$. Let

[TABLE]

where $x^{w}$ denotes the polynomial $f(x)=x^{w}$. Then we have

[TABLE]

• Proof.

Let us first assume that $w\geq m$. Then we have $\mathcal{G}_{p,m-w}(x^{m})=\{1\}$ and therefore

[TABLE]

with Lemma 1. This leads to

[TABLE]

in this case. For the rest of the proof let $w<m$ and additionally we abbreviate $\#\mathcal{G}_{p,m-w}(x^{m})$ by $\#\mathcal{G}$. We write

[TABLE]

In what follows, we refer to the latter sums as

[TABLE]

We may uniquely write any $v\in G_{p,m}\setminus\{0\}$ in the form $v=qx^{m-w}+\ell,$ where $q,\ell\in\mathbb{F}_{q}[x]$ with $\deg(q)<w$ and $\mathrm{deg}(\ell)<m-w$. Using the properties of $X_{p}$ it is clear that $Y_{p^{m},w}(v,x^{m})=Y_{p^{m},w}(\ell,x^{m})$ and hence

[TABLE]

We move on to $S_{2}$. Let $e(\ell):=\max\{k\in\{0,1,\dots,m-w-1\}:x^{k}\mid\ell\}$. With this definition we may display $S_{2}$ as

[TABLE]

In the following, we compute $Y_{p^{m},w}(\ell,x^{m})$ for $\ell\in G_{p,m-w}\setminus\{0\}$ with $e(\ell)=k$. Let $\mu_{p}$ be the Möbius function on the set of monic polynomials over $\mathbb{F}_{p}$, i.e., $\mu_{p}:\mathbb{F}_{p}[x]\rightarrow\{-1,0,1\}$ and

[TABLE]

The fact that $\mu_{p}(1)=1$, $\mu_{p}(x)=-1$ and $\mu_{p}(x^{i})=0$ for $i\in\mathbb{N}$, $i\geq 2$, yields the equivalence of $\sum_{t\mid\gcd(x^{m-w},g)}\mu_{p}(t)=1$ and $\gcd(x^{m-w},g)=1$. Therefore we can write

[TABLE]

The equivalence of the conditions $t\mid h\ell$ and $\frac{t}{\gcd(t,\ell)}\mid h$ yields

[TABLE]

We investigate the inner sum and use Lemma 1 with $a=\frac{t}{\gcd(t,\ell)}$ to find

[TABLE]

Now we have

[TABLE]

From the fact that $e(\ell)=k\leq m-w-1$ we obtain $\gcd(x^{m-w},\ell)=\gcd(x^{m-w-1},\ell)=x^{k}$. This observation leads to

[TABLE]

and

[TABLE]

Altogether we have

[TABLE]

Inserting this result into (7) yields

[TABLE]

Since

[TABLE]

we have

[TABLE]

Summarizing, we have shown

[TABLE]

which completes the proof.

$\hfill\Box$

Now we are ready to prove Theorem 1 using induction on $d$.

• Proof.

We show the result for $d=1$. From Lemma 2 we have

[TABLE]

If $w_{1}\geq m$, then

[TABLE]

If $w_{1}<m$, then we can write

[TABLE]

where we used (6) in the latter step. We regard Lemma 1 with $a=x^{m-w_{1}}$ to compute

[TABLE]

which leads to the desired result also in this case.

Now let $d\in[s-1]$. Assume that we have some $(g_{1},\dots,g_{d})\in\mathcal{G}_{p,m-\boldsymbol{w}}^{d}(x^{m})$ such that

[TABLE]

Let $g^{\ast}\in\mathcal{G}_{p,m-w_{d+1}}(x^{m})$ be such that $R_{\boldsymbol{\gamma}}^{d+1}((x^{w_{1}}g_{1},\dots,x^{w_{d}}g_{d},x^{w_{d+1}}g_{d+1}),x^{m})$ is minimized as a function of $g_{d+1}$ for $g_{d+1}=g^{\ast}$. Then we have

[TABLE]

where

[TABLE]

A minimizer $g^{\ast}$ of $R_{\boldsymbol{\gamma}}^{d+1}((x^{w_{1}}g_{1},\dots,x^{w_{d}}g_{d},x^{w_{d+1}}g_{d+1}),x^{m})$ is also a minimizer of $L(g_{d+1})$. Combining (4) and (Proof.) we obtain that $R_{\boldsymbol{\gamma}}^{d}(\boldsymbol{g},f)\in\mathbb{R}$ for all $d\in[s]$. Moreover with equation (10), established later on in Section 2.2, and the fact that $r_{p}(h,\gamma)>0$ for all $h\in G_{p,m}$ and $\gamma\in(0,1]$, we get that $L(g)\in\mathbb{R}^{+}$ for all $g\in\mathcal{G}_{p,m-w_{d+1}}(x^{m})$. Thus we may bound $L(g^{\ast})$ by the mean over all $g\in\mathcal{G}_{p,m-w_{d+1}}(x^{m})$, hence

[TABLE]

where we used the estimate $\left|X_{p}\left(\frac{v}{x^{m}}hx^{w_{i}}g_{i}\right)\right|\leq 1$ in the last step. With the induction hypothesis and Lemma 3 this leads to

[TABLE]

$\hfill\Box$

**The reduced fast CBC construction

**So far we have seen how to construct a generating vector $\boldsymbol{g}$ of the point set $\mathcal{P}(\boldsymbol{g},x^{m})$. In fact Algorithm 1 can be made much faster using results of [2, 18, 17]. In this section we are investigating and improving Algorithm 1 and additionally analyzing the computational cost of the improved algorithm.

As explained in the following lines Walsh functions are a suitable tool for analyzing the computational cost of CBC algorithms for constructing polynomial lattice point sets. Let $\omega=\mathrm{e}^{2\pi i/p}$, $x\in[0,1)$ and $h$ a nonnegative integer with base $p$ representation $x=x_{1}/p+x_{2}/p^{2}+\ldots$ and $h=h_{0}+h_{1}p+\ldots+h_{r}p^{r}$, respectively. Then we define

[TABLE]

The Walsh function system $\{\mathrm{wal}_{h}\leavevmode\nobreak\ |\leavevmode\nobreak\ h=0,1,\ldots\}$ is a complete orthonormal basis in $L_{2}([0,1))$ which has been used in the analysis of the discrepancy of digital nets (an important class of low-discrepancy point sets which contains polynomial lattice point sets) several times before, see for example [4, 6, 8]. For further information on Walsh functions see [5, Appendix A].

Let $d\geq 1$, $N=p^{m}$. For $P(\boldsymbol{g},f)=\{\boldsymbol{x}_{0},\ldots,\boldsymbol{x}_{p^{m}-1}\}$ with $\boldsymbol{x}_{n}=(x_{n}^{(1)},\ldots,x_{n}^{(s)})$ we have the formula (see [4, Section 4])

[TABLE]

which allows us to rewrite $R^{d}_{\boldsymbol{\gamma}}(\boldsymbol{g},x^{m})$ in the following way

[TABLE]

Note that $r_{p}(h,\gamma)$ is defined as in (4) and we identify the integer in base $p$ representation $h=h_{0}+h_{1}p+\ldots+h_{r}p^{r}$ with the polynomial $h(x)=h_{0}+h_{1}x+\ldots+h_{r}x^{r}$. If we set $\psi(\frac{nx^{w_{i}}g_{i}}{x^{m}}):=\sum_{h=1}^{p^{m}-1}r_{p}(h)\mathrm{wal}_{h}(\phi_{m}(\frac{nx^{w_{i}}g_{i}}{x^{m}}))$ we get that

[TABLE]

where $\eta_{d}(n)=\prod_{i=1}^{d}\left(1+\gamma_{i}+\gamma_{i}\psi(\frac{nx^{w_{i}}g_{i}}{x^{m}})\right)$.

In [4, Section 4] it is proved that we can compute the at most $N$ different values of $\psi(\frac{r}{x^{m}})$ for $r\in G_{p,m}$ in $\mathcal{O}(N)$ operations.

Let us now analyze one step of the reduced CBC Algorithm 1. Assuming we already have found $(g_{1},\ldots,g_{d})\in\mathcal{G}^{d}_{p,m-\boldsymbol{w}}(x^{m})$ we have to minimize

[TABLE]

as a function of $g_{d+1}\in\mathcal{G}_{p,m-w_{d+1}}(x^{m})$. If $w_{d+1}\geq m$ then $g_{d+1}=1$ and we are done. Let now $w_{d+1}<m$. From (9) we have that

[TABLE]

In order to minimize $R^{d+1}_{\boldsymbol{\gamma}}((x^{w_{1}}g_{1},\ldots,x^{w_{d+1}}g_{d+1}),x^{m})$ it is enough to minimize $T_{d}(g):=\sum_{n=0}^{p^{m}-1}\psi(\frac{nx^{w_{d+1}}g}{x^{m}})\eta_{d}(n)$. As in [2, Section 4] we can represent this quantity using some specific $(p^{m-w_{d+1}-1}(p-1)\times N)$-matrix $A$ and exploiting its additional structure. Let therefore

[TABLE]

First of all observe that we get $(T(g))_{g\in G_{p,m-w_{d+1}}(x^{m})}=A\boldsymbol{\eta}_{d}$. Secondly the matrix A is a block matrix and can be written in the following form

[TABLE]

If $\boldsymbol{x}$ is any vector of size $p^{m}$ then we compute

[TABLE]

With this representation we can apply the machinery of [17, 18] and get that multiplication with $\Omega^{(m-w_{d+1})}$ can be done in $\mathcal{O}((m-w_{d+1})p^{m-w_{d+1}})$ operations. Summarizing we have:

Algorithm 2

Compute $\psi(\frac{r}{x^{m}})$ for $r\in G_{p,m}$. 2. 2.

Set $\eta_{1}(n)=\psi(\frac{nx^{w_{1}}g_{1}}{x^{m}})$ for $n=0,\ldots,p^{m}-1$. 3. 3.

Set $g_{1}=1,\leavevmode\nobreak\ d=2$ and $t=\max\{j\in[s]\leavevmode\nobreak\ |\leavevmode\nobreak\ w_{j}<m\}$.

While $d\leq\min\{s,t\}$,

1. (a)

Partition $\eta_{d-1}$ into $p^{w_{d}}$ vectors $\eta_{d-1}^{(1)},\ldots,\eta_{d-1}^{(p^{w_{d}})}$ of length $p^{m-w_{d}}$ and let $\eta^{\prime}=\sum_{i=1}^{p^{w_{d}}}\eta_{d-1}^{(i)}$. 2. (b)

Let $T_{d}(g)=\Omega^{(m-w_{d})}\eta^{\prime}$. 3. (c)

Let $g_{d}=\mathrm{argmin}_{g}T_{d}(g)$. 4. (d)

Let $\eta_{d}(n)=\eta_{d-1}(n)\psi(\frac{nx^{w_{d}}z_{d}}{x^{m}})$ 5. (e)

Increase $d$ by 1. 4. 4.

If $s\geq t$ then set $g_{t}=g_{t+1}\ldots=g_{s}=1$.

Similar to [2] we obtain from the observations in this section the following theorem:

Theorem 2

The cost of Algorithm 2 is

[TABLE]

2.2 Polynomial lattice point sets for irreducible $f$

Finally we want to consider the special case where $f$ is an irreducible polynomial. So, for this section let $f$ be an irreducible polynomial over $\mathbb{F}_{p}$ with $\deg(f)=m$.

Theorem 3

Let $\boldsymbol{\gamma}$ and $\boldsymbol{w}$ as in Theorem 1 and let $f\in\mathbb{F}_{p}[x]$ be an irreducible polynomial with $\deg(f)=m$. Let further $(g_{1},\ldots,g_{s})\in\mathcal{G}_{p,m-\boldsymbol{w}}^{s}(f)$ be constructed according to Algorithm 1. Then we have for every $d\in[s]$

[TABLE]

• Proof.

We will prove this result by induction on $d$. According to Algorithm 1 we know that $g_{1}=1$ for $d=1$. Therefore $R_{\boldsymbol{\gamma}}^{1}((x^{w_{1}}g_{1}),f)=0$ since for all $h\in G_{p,m}$ we have $\deg(h)<m$ and hence the congruence $hx^{w_{1}}\equiv 0\pmod{f}$ has no solutions.

Let $d\in[s-1]$ and assume that we have already found $(g_{1},\ldots,g_{d})\in\mathcal{G}_{p,m-\boldsymbol{w}}^{d}(f)$. For $\boldsymbol{g}=(x^{w_{1}}g_{1},\ldots,x^{w_{d}}g_{d})$ we have from (4) that

[TABLE]

where we proceeded similarly as in the proof of Theorem 1. Here we have

[TABLE]

Let $g^{*}\in\mathcal{G}_{p,m-w_{d+1}}(f)$ be a minimizer of $R^{d+1}_{\boldsymbol{\gamma}}((\boldsymbol{g},x^{w_{d+1}}g_{d+1}),f)$ as a function of $g_{d+1}$. Therefore $g^{*}$ also minimizes $\theta(g_{d+1})$. Bounding $\theta(g^{*})$ by its mean we obtain

[TABLE]

Observe that $\gcd(f,h_{d+1}x^{w_{d+1}})=1$. Therefore the congruence $h_{d+1}x^{w_{d+1}}g_{d+1}\equiv-\boldsymbol{h}\cdot\boldsymbol{g}\pmod{f}$ has a unique solution in $G_{p,m}$ but not necessarily in $\mathcal{G}_{p,m-w_{d+1}}(f)$. In the case that $-\boldsymbol{h}\cdot\boldsymbol{g}\not\equiv 0\pmod{f}$ we conclude that the congruence has at most one solution in $\mathcal{G}_{p,m-w_{d+1}}(f)$. If $-\boldsymbol{h}\cdot\boldsymbol{g}\equiv 0\pmod{f}$ the congruence has no solution in $\mathcal{G}_{p,m-w_{d+1}}(f)$ since $0\not\in\mathcal{G}_{p,m-w_{d+1}}(f)$. Hence we find by an application of [4, Lemma 3.3]

[TABLE]

By (10) and the induction hypothesis we have that

[TABLE]

where we used in the latter step that $\frac{p^{m}}{\#\mathcal{G}_{p,m-w_{d+1}}(f)}\leq\frac{p}{p-1}p^{\min\{w_{d+1},m\}}$. This follows from the fact that $\#\mathcal{G}_{p,m-w_{d+1}}(f)=p^{m-w_{d+1}}-1$ if $w_{d+1}<m$ and $\#\mathcal{G}_{p,m-w_{d+1}}(f)=1$ if $w_{d+1}\geq m$. This finishes the proof of Theorem 3. $\hfill\Box$

As an immediate consequence of (3) and Theorem 3 we obtain the following result.

Corollary 3

Let $N=p^{m}$, $(w_{j})_{j\geq 1}$ be a non-decreasing sequence of nonnegative integers and let $(g_{1},\dots,g_{s})\in\mathcal{G}_{p,m-\boldsymbol{w}}^{s}(f)$ for irreducible $f\in G_{p,m}$ be constructed using Algorithm 1. Then the polynomial lattice point set $\mathcal{P}\left((x^{w_{1}}g_{1},\dotsc,x^{w_{s}}g_{s}),f\right)$ has a weighted star discrepancy

[TABLE]

Remark 3

Using the same argumentation as in Corollary 2 we again obtain the sufficient condition $\sum\limits_{j=1}^{\infty}\gamma_{j}p^{w_{j}}<\infty$ for strong polynomial tractability and for the discrepancy bound $D_{N,\boldsymbol{\gamma}}^{*}\left((x^{w_{1}}g_{1},\dotsc,x^{w_{s}}g_{s}),f\right)=\mathcal{O}(N^{-1+\delta})$, with the implied constant independent of $s$, for any $\delta>0$.

Acknowledgements. We would like to thank Peter Kritzer and Friedrich Pillichs-hammer for their valuable comments and suggestions which helped to improve our paper.