Successive Coordinate Search and Component-by-Component Construction of Rank-1 Lattice Rules

TL;DR

This paper introduces a successive coordinate search (SCS) method for constructing rank-1 lattice rules that improves upon the traditional CBC algorithm, achieving near-optimal worst-case errors with reduced computational cost.

Contribution

It generalizes the CBC algorithm to include a successive coordinate search, leading to better generating vectors and maintaining optimal convergence rates.

Findings

SCS finds better generating vectors than CBC, especially with slower decaying weights.

The error bounds are independent of dimension under certain weight conditions.

The fast SCS algorithm reduces computational complexity to O(d n log n).

Abstract

The (fast) component-by-component (CBC) algorithm is an efficient tool for the construction of generating vectors for quasi-Monte Carlo rank-1 lattice rules in weighted reproducing kernel Hilbert spaces. We consider product weights, which assigns a weight to each dimension. These weights encode the effect a certain variable (or a group of variables by the product of the individual weights) has. Smaller weights indicate less importance. Kuo (2003) proved that the CBC algorithm achieves the optimal rate of convergence in the respective function spaces, but this does not imply the algorithm will find the generating vector with the smallest worst-case error. In fact it does not. We investigate a generalization of the component-by-component construction that allows for a general successive coordinate search (SCS), based on an initial generating vector, and with the aim of getting closer to the smallest worst-case error. The proposed method admits the same type of worst-case error bounds as the CBC algorithm, independent of the choice of the initial vector. Under the same summability conditions on the weights as in [Kuo,2003] the error bound of the algorithm can be made independent of the dimension $d$ and we achieve the same optimal order of convergence for the function spaces from [Kuo,2003]. Moreover, a fast version of our method, based on the fast CBC algorithm by Nuyens and Cools, is available, reducing the computational cost of the algorithm to $O(d \, n \log(n))$ operations, where $n$ denotes the number of function evaluations. Numerical experiments seeded by a Korobov-type generating vector show that the new SCS algorithm will find better choices than the CBC algorithm and the effect is better when the weights decay slower.

Tables

Table 1. Table 1: Weighted unanchored Sobolev space with d = 5 , β j = 1 , γ j = ( 0.95 ) j , q = 100 formulae-sequence 𝑑 5 formulae-sequence subscript 𝛽 𝑗 1 formulae-sequence subscript 𝛾 𝑗 superscript 0.95 𝑗 𝑞 100 d=5,\beta_{j}=1,\gamma_{j}=(0.95)^{j},q=100

$n$	$e_{n,d}\left({𝒛}_{\text{kor}}^{\ast }\right)$	$e_{n,d}\left({𝒛}_{\text{rand}}^{\ast }\right)$	$e_{n,d}\left({𝒛}_{\text{cbc}}\right)$	$e_{n,d}\left({𝒛}_{\text{full}}\right)$
101	2.6003e-02	2.6000e-02	2.6022e-02	2.6000e-02
127	2.1794e-02	2.1834e-02	2.2180e-02	2.1751e-02
139	2.0016e-02	2.0010e-02	2.0493e-02	1.9999e-02
151	1.8886e-02	1.8893e-02	1.9175e-02	1.8843e-02
181	1.5963e-02	1.5937e-02	1.6453e-02	1.5928e-02
199	1.4813e-02	1.4808e-02	1.5368e-02	1.4802e-02

Table 2. Table 2: Weighted unanchored Sobolev space with d = 5 , β j = 1 , γ j = ( 0.7 ) j , q = 100 formulae-sequence 𝑑 5 formulae-sequence subscript 𝛽 𝑗 1 formulae-sequence subscript 𝛾 𝑗 superscript 0.7 𝑗 𝑞 100 d=5,\beta_{j}=1,\gamma_{j}=(0.7)^{j},q=100

$n$	$e_{n,d}\left({𝒛}_{\text{kor}}^{\ast }\right)$	$e_{n,d}\left({𝒛}_{\text{rand}}^{\ast }\right)$	$e_{n,d}\left({𝒛}_{\text{cbc}}\right)$	$e_{n,d}\left({𝒛}_{\text{full}}\right)$
101	1.0721e-02	1.0695e-02	1.0878e-02	1.0695e-02
127	8.7079e-03	8.6296e-03	8.6700e-03	8.6275e-03
139	8.0567e-03	8.0439e-03	8.0724e-03	8.0439e-03
151	7.4913e-03	7.4913e-03	7.5295e-03	7.4913e-03
181	6.26793e-03	6.2594e-03	6.3898e-03	6.2421e-03
199	5.7456e-03	5.7682e-03	5.8758e-03	5.7352e-03

Table 3. Table 3: Weighted Korobov space with d = 100 , α = 1 , β j = 2 3 , γ j = 2 3 ​ ( 0.95 ) j , q = 100 formulae-sequence 𝑑 100 formulae-sequence 𝛼 1 formulae-sequence subscript 𝛽 𝑗 2 3 formulae-sequence subscript 𝛾 𝑗 2 3 superscript 0.95 𝑗 𝑞 100 d=100,\alpha=1,\beta_{j}=\frac{2}{3},\gamma_{j}=\frac{2}{3}(0.95)^{j},q=100

$n$	$\overline{e_{n,d}\left({𝒛}_{\text{kor}}\right)}$	$e_{n,d}\left({𝒛}_{\text{kor}}^{\ast }\right)$	$e_{n,d}\left({𝒛}_{\text{cbc}}\right)$
1009	1.6554e-02	1.6221e-02	1.6566e-02
2003	1.1759e-02	1.1474e-02	1.1719e-02
4001	8.3025e-03	8.1204e-03	8.2869e-03
8009	5.8655e-03	5.7730e-03	5.8500e-03
32003	2.9320e-03	2.8874e-03	2.9301e-03

Table 4. Table 4: Weighted Korobov space with d = 100 , α = 1 , β j = 1 , γ j = ( 0.7 ) j , q = 100 formulae-sequence 𝑑 100 formulae-sequence 𝛼 1 formulae-sequence subscript 𝛽 𝑗 1 formulae-sequence subscript 𝛾 𝑗 superscript 0.7 𝑗 𝑞 100 d=100,\alpha=1,\beta_{j}=1,\gamma_{j}=(0.7)^{j},q=100

$n$	$\overline{e_{n,d}\left({𝒛}_{\text{kor}}\right)}$	$e_{n,d}\left({𝒛}_{\text{kor}}^{\ast }\right)$	$e_{n,d}\left({𝒛}_{\text{cbc}}\right)$
1009	3.1185e-01	3.0834e-01	3.0931e-01
2003	2.0902e-01	2.0661e-01	2.0708e-01
4001	1.3894e-01	1.3713e-01	1.3658e-01
8009	9.1757e-02	9.0445e-02	8.9611e-02
32003	3.9467e-02	3.8763e-02	3.8528e-02

Equations

$$I(f)={\int }_{[0,1{]}^d}f(\mathbf{x})\mathrm{d}\mathbf{x}$$

$$Q_{n,d}(f)=\frac{1}{n}\sum _{k=0}^{n-1}f({\mathbf{x}}_k),$$

$$e_{n,d}(P_n,H)=\underset{\Vert f{\Vert }_H\le 1}{\sup }\left|{\int }_{[0,1{]}^d}f(\mathbf{x})\mathrm{d}\mathbf{x}-\frac{1}{n}\sum _{k=0}^{n-1}f({\mathbf{x}}_k)\right|.$$

$$P_n$$

$${\mathbb{Z}}_n=\left\{0,1,2,\dots ,n-1\right\}.$$

$$e_{n,d}^2(P_n,H)=e_{n,d}^2(\mathbf{z})$$

$$e_{0,d}^2(0,H)=e_{0,d}^2=\prod _{j=1}^d{\beta }_j.$$

$$\frac{e_{n,d}^2(\mathbf{z})}{e_{0,d}^2}$$

$$e_{n,d}^2(\tilde{\mathbf{z}})$$

$$=-\prod _{j=1}^d{\beta }_j+\frac{1}{n}\sum _{k=0}^{n-1}\prod _{j=1}^s\left({\beta }_j+{\gamma }_j\omega \left(\left\{\frac{k{z}_j}{n}\right\}\right)\right)\prod _{j=s+1}^d\left({\beta }_j+{\gamma }_j\omega (0)\right)$$

$$=-\prod _{j=1}^d{\beta }_j+\frac{C_s}{n}\sum _{k=0}^{n-1}\prod _{j=1}^s\left({\beta }_j+{\gamma }_j\omega \left(\left\{\frac{k{z}_j}{n}\right\}\right)\right)$$

$$=-\prod _{j=1}^d{\beta }_j+C_s\left(e_{n,s}^2(\mathbf{z})+\prod _{j=1}^s{\beta }_j\right),$$

$$e_{n,d}^2(z_1,\dots ,z_{s-1},z_s,0^{d-s})$$

$$e_{n,d}({\mathbf{z}}^1)\le e_{n,d}({\mathbf{z}}^0),$$

$${\overline{\beta }}_s=\prod _{j=1}^s{\beta }_j,{\Omega }_n={\left[\omega \left(\left\{\frac{k\mathbf{z}}{n}\right\}\right)\right]}_{\begin{array}{c}z=0,\dots ,n-1\\ k=0,\dots ,n-1\end{array}}$$

$$K_{d,\mathbf{\gamma },\mathbf{\beta }}(\mathbf{x},\mathbf{y})$$

$$r_{\alpha }(\mathbf{h})=\prod _{j=1}^d{r}_{\alpha }(h_j),r_{\alpha }(h)=\{\begin{array}{ll}|h{|}^{2\alpha },& \text{if h=0},\\ 1,& \text{otherwise}.\end{array}$$

$$\Vert f{\Vert }_{K_{d,\mathbf{\gamma },\mathbf{\beta }}}^2=\sum _{\mathbf{h}\in {\mathbb{Z}}^d}|{\hat{f}}_{\mathbf{h}}{|}^2{r}_{\alpha }(\mathbf{h})\prod _{\begin{array}{c}j=1\\ h_j\ne 0\end{array}}^d{\gamma }_j\prod _{\begin{array}{c}j=1\\ h_j=0\end{array}}^d{\beta }_j,$$

$$0\le \frac{1}{n}\sum _{z_s=0}^{n-1}\sum _{\begin{array}{c}\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d\\ \mathbf{h}\cdot \mathbf{z}\equiv 0(\mathrm{mod}n)\end{array}}{r}^{-1}(\mathbf{h})\le \frac{2}{n}\sum _{\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d}{r}^{-1}(\mathbf{h}).$$

$$\frac{1}{n}\sum _{z_s=0}^{n-1}\sum _{\begin{array}{c}\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d\\ \mathbf{h}\cdot \mathbf{z}\equiv 0(\mathrm{mod}n)\end{array}}{r}^{-1}(\mathbf{h})=\frac{1}{n}\sum _{z_s=0}^{n-1}\sum _{\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d}{r}^{-1}(\mathbf{h})\left[\frac{1}{n}\sum _{k=0}^{n-1}\prod _{j=1}^d{\mathrm{e}}^{2\pi \mathrm{i}k{h}_j{z}_j/n}\right].$$

$$\frac{1}{n}\sum _{z=0}^{n-1}{\mathrm{e}}^{2\pi \mathrm{i}khz/n}=\{\begin{array}{ll}1,& \text{if k=0 or }h\equiv 0(\mathrm{mod}n),\\ 0,& \text{otherwise},\end{array}$$

$$\frac{1}{n}\sum _{k=0}^{n-1}\frac{1}{n}\sum _{z=0}^{n-1}{\mathrm{e}}^{2\pi \mathrm{i}khz/n}=\{\begin{array}{ll}1,& \text{if }h\equiv 0(\mathrm{mod}n),\\ 1/n,& \text{if }h\not\equiv 0(\mathrm{mod}n).\end{array}$$

$$\sum _{\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d}{r}^{-1}(\mathbf{h})\frac{1}{n}\sum _{k=0}^{n-1}\left(\prod _{s\ne j=1}^d{\mathrm{e}}^{2\pi \mathrm{i}k{h}_j{z}_j/n}\right)\frac{1}{n}\sum _{z_s=0}^{n-1}{\mathrm{e}}^{2\pi \mathrm{i}k{h}_s{z}_s/n}$$

$$\le \sum _{\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d}{r}^{-1}(\mathbf{h})\frac{1}{n}\sum _{k=0}^{n-1}\left(\prod _{s\ne j=1}^d\left|{\mathrm{e}}^{2\pi \mathrm{i}k{h}_j{z}_j/n}\right|\right)\left|\frac{1}{n}\sum _{z_s=0}^{n-1}{\mathrm{e}}^{2\pi \mathrm{i}k{h}_s{z}_s/n}\right|$$

$$\le \sum _{\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d}{r}^{-1}(\mathbf{h})\frac{1}{n}\sum _{k=0}^{n-1}\frac{1}{n}\sum _{z_s=0}^{n-1}{\mathrm{e}}^{2\pi \mathrm{i}k{h}_s{z}_s/n}$$

$$=\sum _{\begin{array}{c}\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d\\ h_s\equiv 0(\mathrm{mod}n)\end{array}}{r}^{-1}(\mathbf{h})+\frac{1}{n}\sum _{\begin{array}{c}\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d\\ h_s\not\equiv 0(\mathrm{mod}n)\end{array}}{r}^{-1}(\mathbf{h})$$

$$\le \frac{1}{n}\sum _{\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d}{r}^{-1}(\mathbf{h})+\frac{1}{n}\sum _{\begin{array}{c}\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d\\ h_s\not\equiv 0(\mathrm{mod}n)\end{array}}{r}^{-1}(\mathbf{h})\le \frac{2}{n}\sum _{\mathbf{h}\in (\mathbb{Z}\setminus \{0\}{)}^d}{r}^{-1}(\mathbf{h}).$$

$$e_{n,d}^2({\mathbf{z}}^{\ast })$$

$$C_{d,\lambda }$$

$$\sum _{j=1}^{\infty }\frac{{\gamma }_j^{1/\lambda }}{{\beta }_j^{1/\lambda }}<\infty$$

$$e_{n,d}^2(\mathbf{z})=\sum _{0\ne \mathbf{h}\in {\mathbb{Z}}^d}{\left|\frac{1}{n}\sum _{k=0}^{n-1}{\mathrm{e}}^{2\pi \mathrm{i}k\mathbf{h}\cdot \mathbf{z}/n}\right|}^2{r}_{\alpha }^{-1}(\mathbf{h})\prod _{\begin{array}{c}j=1\\ h_j\ne 0\end{array}}^d{\gamma }_j$$

$$g_{\mathfrak{u}}({\mathbf{z}}_{\mathfrak{u}})$$

$$e_{n,d}^2(\mathbf{z})$$

$${\theta }_s(\mathbf{z})$$

Full text

1 Introduction

In this article we study the numerical approximation of integrals of the form

[TABLE]

for $d$-variate functions $f$ via quasi-Monte Carlo quadrature rules. Quasi-Monte Carlo rules are equal-weight quadrature rules of the form

[TABLE]

where the quadrature points ${\bm{x}}_{0},\ldots,{\bm{x}}_{n-1}\in[0,1]^{d}$ are chosen deterministically. Here, we consider integrands $f:[0,1]^{d}\to{\mathbb{R}}$ which belong to some normed function space $(H,\|\cdot\|_{H})$. In order to assess the quality of a particular QMC rule $Q_{n,d}$ with underlying point set $P_{n}=\{{\bm{x}}_{0},\ldots,{\bm{x}}_{n-1}\}$, we introduce the notion of the so-called worst-case error, see, e.g., [1], defined by

[TABLE]

In other words, $e_{n,d}(P_{n},H)$ is the worst error that is attained over all functions in the unit ball of $H$ using the quasi-Monte Carlo rule with quadrature points in $P_{n}$. It is often possible to obtain explicit expressions to calculate $e_{n,d}(P_{n},H)$, see, e.g., [4]. In particular, we consider weighted Korobov and weighted shift-averaged Sobolev spaces, which are both reproducing kernel Hilbert spaces, for details see, e.g., [8, 9, 3, 1]. In this paper we will limit ourselves to the original choice of “product weights”. In essence, the idea is to quantify the varying importance of the coordinate directions $x_{j}$ with $j=1,\ldots,d$ w.r.t. the function values by a sequence ${\bm{\gamma}}=\{\gamma_{j}\}_{j=1}^{d}$ of positive weights.

There are many ways to choose the underlying point set $P_{n}$ of a QMC rule, ranging from lattice rules and sequences, digital nets and sequences and more recent constructions such as interlaced polynomial lattice rules. In this paper, however, we will restrict ourselves to rank-1 lattice rules. This type of QMC rules has an underlying point set $P_{n}\subseteq[0,1]^{d}$ of the form

[TABLE]

where ${\bm{z}}\in{\mathbb{Z}}^{d}$ is the generating vector of the rank-1 lattice rule and $\{\cdot\}$ denotes the fractional part, componentwise if applied to a vector. It is clear that any vector congruent modulo $n$ is equivalent and so we only consider values modulo $n$. The components of ${\bm{z}}$ are often restricted to the set of integers in $\{1,\ldots,n-1\}$ that are relatively prime to $n$, see, e.g., [3], such that one obtains $n$ distinct points for all one-dimensional projections, and as such for any projection. Here we consider generating vectors ${\bm{z}}\in{\mathbb{Z}}_{n}^{d}$ with n prime and

[TABLE]

For rank-1 lattice rules in a weighted shift-invariant tensor-product reproducing kernel Hilbert space $H(K)$ with reproducing kernel $K({\bm{x}},{\bm{y}})=\prod_{j=1}^{d}(\beta_{j}+\gamma_{j}\,\omega(x_{j}-y_{j}))$ the squared worst-case error can be written as

[TABLE]

with positive weights ${\bm{\gamma}}=\{\gamma_{j}\}_{j=1}^{d}$ and ${\bm{\beta}}=\{\beta_{j}\}_{j=1}^{d}$, and where $\int_{0}^{1}\omega(t)\,{\mathrm{d}}t=0$, see, e.g., [5]. We note that the weights $\beta_{j}$ are to easily accommodate for some types of shift-averaged Sobolev spaces. Moreover, the initial squared worst-case error in this function space, i.e., using $n=0$ samples and with the convention that $Q_{0,d}(f)=0$, is given by

[TABLE]

Remark 1**.**

It is always possible to consider the normalized worst-case error by dividing by the initial worst-case error for the zero-algorithm. The squared normalized worst-case error then takes the form

[TABLE]

with $e^{2}_{0,d}=\prod_{j=1}^{d}\beta_{j}$. This is equivalent to considering the worst-case error $e_{n,d}({\bm{z}})$ with modified weight sequences $\hat{\beta}_{j}=1$ and $\hat{\gamma}_{j}=\gamma_{j}/\beta_{j}$.

One of the most commonly considered methods to construct good rank-1 lattice rules is the component-by-component (CBC) construction, see, e.g., [4], which extends the generating vector one component at a time by selecting the next components $z_{s}$ which minimizes the worst-case error of the $s$-dimensional rule. The pseudo-code of the CBC algorithm is given below.

It was shown in [3] that the component-by-component construction generates lattice rules which achieve optimal rates of convergence in weighted Korobov and Sobolev function spaces. Additionally, a fast construction method is available, see [5, 6], that reduces the construction cost to $O(d\,n\,\log(n))$ operations.

Even though the CBC algorithm constructs generating vectors ${\bm{z}}$ which exhibit the optimal error asymptotics, the constructed vector is not necessarily the one minimizing the worst-case error $e_{n,d}({\bm{z}})$. We will therefore introduce and investigate a different construction method which can generate lattice rules with a smaller worst-case error than the CBC construction.

The article is structured as follows. In Section 2 we introduce the successive coordinate search (SCS) algorithm and analyse some properties. In Section 3 we prove that the SCS construction achieves optimal rates of convergence in the weighted Korobov and weighted shift-averaged Sobolev space. To get dimension-independent bounds, i.e., achieve tractability, we show that the summability condition on the weights is the same as for the normal CBC construction. Finally we report on various numerical experiments in Section 4.

2 Formulation of the successive coordinate search algorithm

In this section we introduce an algorithm of similar nature to the component-by-component construction. One advantage of the component-by-component construction is that the algorithm is extensible in the dimension $d$, i.e., to find the $(d+1)$-dimensional generating vector, the algorithm does not need to restart but just starts from the generating vector of dimension $d$. In our setting this also implies that a $d$-dimensional vector, with $d$ large enough, could be constructed and used for all problems with less than $d$ dimensions. This allows us to fix the maximum number of dimensions to some large enough $d$ and successively try to find the best $s$-th component of a $d$-dimensional generating vector, keeping all other $d-1$ choices fixed. The pseudocode of the successive coordinate search (SCS) algorithm is given below.

Instead of increasing the dimension in each step of the algorithm, we keep $d$ fixed during all calculations. Based on a starting vector ${\bm{z}}^{0}\in{\mathbb{Z}}_{n}^{d}=\{0,1,\ldots,n-1\}^{d}$, the algorithm successively selects the coordinate $z_{s}\in{\mathbb{Z}}_{n}$ which minimizes the squared worst-case error $e^{2}_{n,d}(z_{1},\ldots,z_{s-1},z_{s},z_{s+1}^{0},\ldots,z_{d}^{0})$ while keeping all other coordinates of ${\bm{z}}$ fixed. Thus, in the process of the SCS algorithm the coordinates of the starting vector ${\bm{z}}^{0}$ are altered in each step of the algorithm. Our construction is very similar to the component-by-component construction, with the only difference being that an initial vector ${\bm{z}}^{0}$ is required as input for the algorithm. In fact, we can prove that the successive coordinate search algorithm is a generalized version of the CBC algorithm as the following theorem shows by starting from an initial vector ${\bm{z}}^{0}=(0,\ldots,0)\in{\mathbb{Z}}_{n}^{d}$. We note that this is a degenerate vector as it generates only a $1$-point rule, and thus is in some sense the worst possible choice for any $n\geq 1$.

Theorem 1**.**

The component-by-component (CBC) algorithm and the successive coordinate search (SCS) algorithm with starting vector ${\bm{z}}^{0}=(0,\ldots,0)$ both yield the same generating vector as outcome (with equivalent choices selected in the same way in both algorithms).

Proof.

Denote by $0^{r}$ the $r$-dimensional zero vector, where $1\leq r\leq d$. For an arbitrary ${\bm{z}}\in{\mathbb{Z}}_{n}^{s}$ with $1\leq s\leq d$ and with $\tilde{{\bm{z}}}=({\bm{z}},0^{d-s})\in{\mathbb{Z}}_{n}^{d}$, the squared worst-case error equals

[TABLE]

where $C_{s}=\prod_{j=s+1}^{d}\,(\beta_{j}+\gamma_{j}\,\omega(0))$. Note that due to the non-negativity of the squared worst-case error $e_{n,d}^{2}$ the function $\omega$ is such that $\omega(0)\geq 0$ and so the constants $C_{s}$ are positive for all $s=1,\ldots,d$.

Now, in each step $1\leq s\leq d$ of the SCS algorithm with initial vector ${\bm{z}}^{0}=(0,\ldots,0)$, we search for the $z_{s}\in{\mathbb{Z}}_{n}$ that minimizes $e_{n,d}^{2}(z_{1},\ldots,z_{s-1},z_{s},0^{d-s})$, where $z_{1},\ldots,z_{s-1}$ have been determined in the previous steps of the algorithm. By the above identity we have that

[TABLE]

and so, since the remaining terms on the right-hand side are independent of $z_{s}$, this is equivalent to finding $z_{s}\in{\mathbb{Z}}_{n}$ such that $e_{n,s}^{2}(z_{1},\ldots,z_{s-1},z_{s})$ is minimized. As this is exactly the same quantity which is minimized in each step of the component-by-component construction algorithm, we see that the CBC algorithm and the SCS algorithm with starting vector ${\bm{z}}^{0}=(0,\ldots,0)$ yield exactly the same outcome under the assumption that both algorithms select the same minimizer whenever multiple choices occur in a minimization step. ∎

Furthermore, the formulation of the successive coordinate search construction guarantees that the generating vector ${\bm{z}}$ obtained by the SCS algorithm with initial vector ${\bm{z}}^{0}$ is never worse than the input vector ${\bm{z}}^{0}$.

Proposition 1**.**

Let ${\bm{z}}^{0}\in{\mathbb{Z}}_{n}^{d}$ be an arbitrary generating vector for a rank-1 lattice rule and denote by ${\bm{z}}^{1}\in{\mathbb{Z}}_{n}^{d}$ the generating vector constructed by the SCS algorithm with starting vector ${\bm{z}}^{0}$. Then we have that

[TABLE]

i.e., the SCS method constructs a generating vector with worst-case error smaller than or equal to the worst-case error of the initial vector.

Proof.

The statement follows directly from the formulation of the algorithm. ∎

Similar to the case of the component-by-component construction there is a fast version available that allows for the construction of generating vectors with time complexity $O(d\,n\,\log(n))$. In case $n$ is a prime number this results in the following algorithm.

Here we used the notations

[TABLE]

and $\Omega_{n}^{\langle g\rangle}$ denotes the reordering of $\Omega_{n}$ w.r.t. a generator $g$ for the cyclic group of ${\mathbb{Z}}_{n}$. For more details see [5, 4]. The symbols $\mathbin{.*}$ and $\mathbin{./}$ denote componentwise vector multiplication and division, respectively, and $\Omega_{n}^{\langle g\rangle}\!(j,:)$ stands for the $j$-th row of $\Omega_{n}^{\langle g\rangle}$. Note that the computation is only slightly more expensive than the fast CBC algorithm since $\mathbf{q}$ has to be initialized and updated using ${\bm{z}}^{0}$, the computational complexity is still $O(d\,n\log(n))$.

3 Error Bounds for the SCS Algorithm

In this section we derive worst-case error bounds and show that the previously introduced successive coordinate search construction achieves optimal convergence rates in the respective function space. Here, we consider two of the most common function spaces in QMC theory, the weighted Korobov space and the weighted shift-averaged (anchored) Sobolev space.

3.1 The weighted Korobov space

Let ${\bm{\gamma}}=\{\gamma_{j}\}$ and ${\bm{\beta}}=\{\beta_{j}\}$ be two weight sequences. The reproducing kernel of the corresponding $d$-dimensional weighted Korobov space is then given by

[TABLE]

where $\alpha>\tfrac{1}{2}$ is referred to as the smoothness parameter and we define

[TABLE]

For integer $\alpha$ the smoothness can be interpreted as the number of mixed partial derivatives $f^{(\tau_{1},\ldots,\tau_{d})}$ with $(\tau_{1},\ldots,\tau_{d})\leq(\alpha,\ldots,\alpha)$ that exist and are square-integrable. The space consists of functions which can be represented as absolutely summable Fourier series with norm

[TABLE]

where the $\hat{f}_{\bm{h}}$ denote the Fourier coefficients of $f$.

We prove that the successive coordinate search (SCS) algorithm achieves the optimal rate of convergence for multivariate integration in the weighted Korobov space for any initial vector. As is usual practice, we restrict ourselves to a prime number of points to simplify the needed proof techniques. We need the following lemma in the proof of the theorem.

Lemma 1**.**

For $s\in\{1,\ldots,d\}$, $n$ prime, arbitrary integers $z_{j}$, $j\in\{1\mathpunct{:}d\}\setminus\{s\}$, and $r({\bm{h}})=\prod_{j=1}^{d}r(h_{j})$ with $r(h)>0$ such that for $h\neq 0$ we have $r(nh)\geq n^{c}\,r(h)$ for $c\geq 1$, then

[TABLE]

Proof.

The condition ${\bm{h}}\cdot{\bm{z}}\equiv 0\pmod{n}$ can be written equivalently by

[TABLE]

Further, for $0\leq k<n$ and $n$ prime

[TABLE]

and hence for $h\in{\mathbb{Z}}$

[TABLE]

Thus

[TABLE]

Which completes the proof. ∎

Theorem 2**.**

Let $n$ be a prime number and ${\bm{z}}^{0}=(z_{1}^{0},\ldots,z_{d}^{0})\in{\mathbb{Z}}_{n}^{d}$ be an arbitrary initial vector. Furthermore, denote by ${\bm{z}}^{*}=(z_{1}^{*},\ldots,z_{d}^{*})\in{\mathbb{Z}}_{n}^{d}$ the generating vector constructed by the successive coordinate search method with initial vector ${\bm{z}}^{0}$. Then the squared worst-case error $e_{n,d}^{2}({\bm{z}}^{*})$ in the Korobov space with kernel $K_{d,{\bm{\gamma}},{\bm{\beta}}}$ satisfies

[TABLE]

where the constant $C_{d,\lambda}$ is given by

[TABLE]

with $\mu_{\alpha,\lambda}=2\zeta(2\alpha/\lambda)\,.$ Additionally, if the weights satisfy the summability condition

[TABLE]

then $C_{d,\lambda}\leq C_{\lambda}<\infty$, and the constant $C_{\lambda}$ is bounded independent of the dimension $d$. Hence, the worst-case error $e_{n,d}({\bm{z}}^{*})$ can be taken arbitrarily close to $O(n^{-\alpha})$, with the implied constant independent of $n$, and independent of $d$ if the summability condition holds.

Proof.

We use the notation from [4]: for a subset ${\mathfrak{u}}\subseteq\{1\mathpunct{:}d\}=\{1,\ldots,d\}$, we define the set ${\mathbb{Z}}_{\mathfrak{u}}=\{{\bm{h}}\in{\mathbb{Z}}^{d}:h_{j}\neq 0\text{ for all }j\in{\mathfrak{u}}\text{ and }h_{j}=0\text{ for }j\not\in{\mathfrak{u}}\}$, and define the dual lattice $L^{\perp}_{\mathfrak{u}}({\bm{z}},n)=L^{\perp}_{\mathfrak{u}}({\bm{z}}_{\mathfrak{u}},n)=\{{\bm{h}}\in{\mathbb{Z}}_{\mathfrak{u}}:{\bm{h}}_{\mathfrak{u}}\cdot{\bm{z}}_{\mathfrak{u}}\equiv 0\pmod{n}\}$ where we write ${\bm{z}}_{\mathfrak{u}}$ and ${\bm{h}}_{\mathfrak{u}}$ to refer only to those components in ${\bm{z}}$ and ${\bm{h}}$. For ${\bm{h}}\in{\mathbb{Z}}_{\mathfrak{u}}$ we will write ${\bm{h}}_{\mathfrak{u}}\in{\mathbb{Z}}_{\mathfrak{u}}$ and $r_{\alpha}({\bm{h}}_{\mathfrak{u}})$ to explicitly denote the dependence on the dimensions in ${\mathfrak{u}}$ only. We also write $\gamma_{u}=\prod_{j\in u}\gamma_{j}$ and set $\gamma_{\emptyset}=1$. Now, without loss of generality, we consider the case where $\beta_{j}=1$ for all $j$ and correct the final expression afterwards, see Remark 1. Then from (1), or, see, e.g., [4, p. 5, eq. (6)] with $q=2$ and $\varphi({\bm{x}}_{k})=\mathrm{e}^{2\pi\mathrm{i}\,k\,{\bm{h}}\cdot{\bm{z}}/n}$, we have

[TABLE]

Now define

[TABLE]

which gives

[TABLE]

Minimizing $e_{n,d}({\bm{z}})$ over $z_{s}\in{\mathbb{Z}}_{n}$ is equivalent to minimizing only those parts which depend on $z_{s}$, resulting in the auxiliary target function

[TABLE]

We note that in the standard CBC proofs this quantity only depends on the dimensions up to $s$ while here it depends on all $d$ dimensions. Obviously

[TABLE]

where the tilde on top of ${\bm{z}}$ means that in replacing the sum over $s\in{\mathfrak{u}}\subseteq\{1\mathpunct{:}s\}$ by the sum over $s\in{\mathfrak{u}}\subseteq\{1\mathpunct{:}d\}$ we choose arbitrary $z_{j}$ for $j>s$. We are free to do so since we are just adding positive quantities. Furthermore, for $1\leq\lambda<\infty$, using the so-called Jensen’s inequality, we obtain

[TABLE]

which holds for all choices of ${\bm{w}}$. Since in minimizing $e^{2}_{n,d}(z_{1}^{*},\ldots,z_{s-1}^{*},z_{s},z_{s+1}^{0},\ldots,z_{d}^{0})$ we are in fact minimizing $\theta_{s}(z_{1}^{*},\ldots,z_{s-1}^{*},z_{s},z_{s+1}^{0},\ldots,z_{d}^{0})$ we now use the standard reasoning that the best choice $z_{s}=z_{s}^{*}$ makes $\theta_{s}$ at least as small as the average over all choices, and the same reasoning holds if we raise $\theta_{s}$ to the power $1/\lambda$. Therefore we obtain

[TABLE]

where we used Jensen’s inequality to obtain the second line (and where $\bar{{\bm{z}}}_{\mathfrak{u}}$ means we take $\bar{\bm{z}}=(z_{1}^{*}$,…,$z_{s-1}^{*},z_{s},z_{s+1}^{0},\ldots,z_{s+1}^{d})$) and Lemma 1, relabeling the set ${\mathfrak{u}}$ to be $\{1,\ldots,|{\mathfrak{u}}|\}$, $d=|{\mathfrak{u}}|$, and with $r(h)=|h|^{2\alpha/\lambda}$ and $c=2\alpha/\lambda\geq 1$, to obtain the last line. For convenience we define

[TABLE]

from which it follows that $2\alpha/\lambda>1$ and we thus need $\lambda<2\alpha$. Since for ${\bm{h}}_{\mathfrak{u}}\in{\mathbb{Z}}_{\mathfrak{u}}$ we have $r_{\alpha}^{-1/\lambda}({\bm{h}}_{\mathfrak{u}})=\prod_{j\in{\mathfrak{u}}}r_{\alpha}^{-1/\lambda}(h_{j})$ we find $\sum_{{\bm{h}}_{\mathfrak{u}}\in{\mathbb{Z}}_{{\mathfrak{u}}}}r_{\alpha}^{-1/\lambda}({\bm{h}}_{\mathfrak{u}})=\mu_{\alpha,\lambda}^{|{\mathfrak{u}}|}$.

In each step of the SCS algorithm we now have a bound on $\theta_{s}^{1/\lambda}$ which we insert in our bound for the worst-case error, each time choosing the components of ${\bm{z}}^{0}$ for ${\bm{w}}$, to obtain

[TABLE]

To show that the summability condition $\sum_{j=1}^{\infty}\gamma_{j}^{1/\lambda}<\infty$ gives a bound independent of $d$ we note that

[TABLE]

Now using that $\log(1+x)\leq x$ for $x>-1$, we find that

[TABLE]

which implies the result. ∎

3.2 The weighted Sobolev space

Again, let ${\bm{\gamma}}=\{\gamma_{j}\}$ and ${\bm{\beta}}=\{\beta_{j}\}$ be two weight sequences. There is a close relationship between the weighted Korobov space with smoothness parameter $\alpha=1$ and the shift-averaged weighted Sobolev space. The shift-invariant kernel of the weighted Sobolev space with anchor ${\bm{a}}=(a_{1},\ldots,a_{d})$ of $d$-variate functions is given by

[TABLE]

where $\hat{\beta}_{j}=\beta_{j}+\gamma_{j}\left(a_{j}^{2}-a_{j}+\frac{1}{3}\right)$ and $\hat{\gamma}_{j}=\frac{\gamma_{j}}{2\pi^{2}}$ for anchor values $a_{j}\,$. Furthermore, for $c_{j}=a_{j}^{2}-a_{j}+\frac{1}{3}$ the shift-averaged squared worst-case error $\hat{e}_{n,d}^{2}({\bm{z}})$ with generating vector ${\bm{z}}$ takes the following form

[TABLE]

Additionally, the initial worst-case error in the weighted Sobolev space is given by

[TABLE]

Since these are precisely the worst-case error expressions as for the weighted Korobov space with $\alpha=1$ and weights $\hat{\beta_{j}}$ and $\hat{\gamma_{j}}$, we obtain similar error bounds as before.

Theorem 3**.**

Let $n$ be a prime number and ${\bm{z}}^{0}=(z_{1}^{0},\ldots,z_{d}^{0})\in{\mathbb{Z}}_{n}^{d}$ be an arbitrary initial vector. Furthermore, denote by ${\bm{z}}^{*}=(z_{1}^{*},\ldots,z_{d}^{*})\in{\mathbb{Z}}_{n}^{d}$ the generating vector constructed by the successive coordinate search method with initial vector ${\bm{z}}^{0}$. Then the squared worst-case error $\hat{e}^{2}_{n,d}({\bm{z}}^{*})$ in the shift-averaged (anchored) Sobolev space with kernel $K^{*}_{d,{\bm{\gamma}},{\bm{\beta}}}$ satisfies

[TABLE]

where the constant $\hat{C}_{d,\lambda}$ is given by the expression for $C_{d,\lambda}$ from Theorem 2 with $\alpha=1$ and weights $\hat{\beta_{j}}=\beta_{j}+\gamma_{j}\left(a_{j}^{2}-a_{j}+\frac{1}{3}\right)$ and $\hat{\gamma}_{j}=\frac{\gamma_{j}}{2\pi^{2}}$.

Additionally, if the weights satisfies the summability condition

[TABLE]

then $\hat{C}_{d,\lambda}\leq\hat{C}_{\lambda}<\infty$, and the constant $C_{\lambda}$ is bounded independent of the dimension $d$. Hence, the worst-case error $\hat{e}_{n,d}({\bm{z}}^{*})$ can be taken arbitrarily close to $O(n^{-1})$, with the implied constant independent of $n$, and independent of $d$ if the summability condition holds.

Proof.

The theorem follows directly from the previous result in Theorem 2. ∎

4 Numerical results and experiments

The idea regarding the SCS algorithm is to obtain generating vectors with smaller error values than obtained by the CBC algorithm, provided we choose a suitable initial vector ${\bm{z}}^{0}\in{\mathbb{Z}}_{n}^{d}$. The formulation of the algorithm suggests that the performance of the SCS construction strongly depends on the starting vector ${\bm{z}}^{0}$ which we select beforehand. In this section we conduct some numerical experiments in the same setting as for the CBC algorithm in order to assess the performance of the SCS algorithm. For the experiments we prefer to have $n$ distinct points in each dimension and so restrict our generating vector choices to exclude the choice $z_{s}=0$ for the components of ${\bm{z}}$ for prime $n$, i.e., ${\bm{z}}\in\{1,\ldots,n-1\}^{d}$. Allowing the choice $z_{s}=0$ has effect on the results which depend on the weights since the CBC algorithm can now pick a zero component when the weights $\gamma_{j}$ decay too slow.

4.1 Construction methods

As we do not know how to best choose the initial vectors for the SCS algorithm, we propose to start from randomly selected initial vectors. This is different from the randomized CBC construction, see, e.g., [7], where in each minimization step the number of possible candidates $z_{s}$ is restricted to $r$ random integers in $\{1,\ldots,n-1\}$. We consider the following two methods.

1. Uniform random vectors + SCS algorithm: Choose $q$ initial vectors ${\bm{z}}^{0}\in{\mathbb{Z}}_{n}^{d}$ at random, apply the fast SCS algorithm to them and then select the one with the smallest worst-case error $e_{n,d}({\bm{z}})$.

2. Korobov-type generating vector + SCS algorithm: Take $q$ randomly chosen Korobov-type generating vectors ${\bm{z}}^{0}={\bm{z}}(a)\equiv(a^{0},a^{1},\ldots,a^{d-1})\pmod{n}$, with $a\in\{1,\ldots,n-1\}$, as initial vectors, apply the fast SCS algorithm to them and then select the one with the smallest worst-case error $e_{n,d}({\bm{z}})$.

As the successive coordinate search algorithm has time complexity $O(d\,n\,\log(n))$, both proposed construction methods have time complexity $O(q\,d\,n\,\log(n))$.

Remark 2**.**

The obvious candidate for the initial vector ${\bm{z}}^{0}$ would of course be the generating vector constructed by the CBC method since by Proposition 1 one would construct ${\bm{z}}^{1}$ such that $e_{n,d}({\bm{z}}^{1})\leq e_{n,d}({\bm{z}}^{0})$. However, experiments show that in most cases the CBC vector is a fixed point with respect to the SCS method, i.e., applying the SCS algorithm to the CBC vector ${\bm{z}}^{0}$ leaves the coordinates of ${\bm{z}}^{0}$ unchanged. Thus, this approach yields usually no further improvement.

4.2 Exhaustive search in low dimensions

In order to test the effectivity of our method we perform some numerical experiments in low dimensions and for a low number of points. Here, we can compute the best generating vector for the respective function space via an exhaustive search over the full set ${\mathbb{Z}}_{n}^{d}$ and then compare its worst-case error to the error values of the generating vectors obtained by our method.

For the weighted unanchored Sobolev space the squared worst-case error is given by

[TABLE]

where $B_{2}(x)=x^{2}-x+\frac{1}{6}$ denotes the Bernoulli polynomial of degree two. Furthermore, ${\bm{z}}_{\text{full}}$ denotes the generating vector obtained by the full exhaustive search, ${\bm{z}}_{\text{cbc}}$ denotes the generating vector obtained via the component-by-component construction and ${\bm{z}}^{*}_{\text{rand}}$ and ${\bm{z}}^{*}_{\text{kor}}$ are the best generating vectors obtained out of $q=100$ initial random choices by the above construction methods 1 and 2, respectively. For two different weight sequences $\gamma_{j}$ and a selection of prime $n$ we obtain the results in Tables 1 and 2, where $\gamma_{j}=(0.95)^{j}$ and $\gamma_{j}=(0.7)^{j}$, respectively. To be able to find the global minimum ${\bm{z}}_{\text{full}}$ over the whole set we limited the dimensionality to $d=5$ and the number of points to $n\leq 199$. This leads to exhaustive searches over about $6$ to $96$ million possible choices for ${\bm{z}}$, where we used the symmetry of the kernel and the fact that we only need to consider generating vectors with $z_{1}=1$ since multiplication by the multiplicative inverse of the first component normalizes any generating vector to have $z_{1}=1$ and this is just a reordering of the cubature nodes.

The results in Table 1 and 2 show that, even for a moderate value of $q$, the randomized SCS method generates lattice rules which have a smaller worst-case error than the one obtained via the CBC construction. Additionally, we see that our method generates worst-case errors that lie in the region of the smallest worst-case error $e_{n,d}({\bm{z}}_{\text{full}})$ and sometimes even constructs the best possible lattice rule. Although we only show two small tables here, similar results were observed for other test cases as well. In particular, we considered weight sequences of the form $\gamma_{j}=q^{j}$ with $0<q<1$ and $\gamma_{j}=j^{-k}$ with $k\in\{2,3,4\}$, for additional results see [2]. The experiments showed that the SCS algorithm outperforms the CBC construction when the decay of the weight sequence $\gamma_{j}$ is slow.

4.3 Numerical experiments for higher dimensions

In higher dimensions and/or for higher number of points it is not possible to perform an exhaustive search in order to obtain a reference value to measure the quality of the constructed generating vectors. Thus, we compare the outcome of the SCS method with the generating vector constructed by the CBC algorithm. Additionally, the empirical numerical results suggested that the use of Korobov-type initial vectors is to be preferred over uniform random vectors and we will therefore only consider Korobov-type initial vectors in this section. We denote by $\overline{e_{n,d}({\bm{z}}_{\text{kor}})}$ the average over the $q$ random choices of the worst-case errors of the SCS constructed vectors ${\bm{z}}_{\text{kor}}$ and with $e_{n,d}({\bm{z}}_{\text{kor}}^{*})$ the best over the $q$ random choices.

The numerical results presented in Tables 3 and 4 are for a Korobov space with $d=100,\alpha=1$ and two different choices of weights, being $\beta_{j}=\frac{2}{3}$ and $\gamma_{j}=\frac{2}{3}(0.95)^{j}$, and $\beta_{j}=1$ and $\gamma_{j}=(0.7)^{j}$, respectively, both with $q=100$ random initial Korobov-type vectors. Our experiments show that the SCS method can construct good lattice rules for high dimensions and large $n$. For our choice of parameters, the SCS algorithm performs moderately better than the CBC construction when the weight sequence ${\bm{\gamma}}=\{\gamma_{j}\}_{j=1}^{d}$ is slowly decaying, as can be seen by comparing the relative difference between $e_{n,d}({\bm{z}}_{\text{kor}}^{*})$ and $e_{n,d}({\bm{z}}_{\text{cbc}})$ for the two different weight sequences in Table 3 and 4. For a more extensive analysis of this behaviour we refer again to [2] where a wider range of weight sequences is considered.

Fig. 1 illustrates the performance of the SCS method compared to the CBC method. The blue dots represent the worst-case error values of lattice rules with $n=4001$ points constructed by the SCS method with $q=300$ Korobov-type initial vectors. The minimal error amongst the constructed lattice rules and the average over the $q$ random seed choices is indicated by the red or black line, respectively. The error corresponding to the CBC method is indicated by the green line. From the figure it becomes evident that the CBC algorithm outperforms the average of the SCS algorithm applied to randomly selected Korobov-type rules, but the best SCS results clearly win over the generating vector constructed by the CBC method.

5 Conclusion

The results and experiments in the previous section, see [2] for additional results, showed that it is possible to use the successive coordinate search algorithm to construct good generating vectors for rank-1 lattice rules. They also confirmed that randomized methods based on the SCS construction can provide generating vectors with smaller worst-case errors than the CBC vector. However, the computational cost of the SCS method can be several times higher while the gained improvement depends on the weight sequence ${\bm{\gamma}}$. Future research could help to find a selection criterion for the starting vector ${\bm{z}}^{0}$ in order to reduce the construction cost of the SCS algorithm. The SCS algorithm should further be regarded as a generalization of the existing component-by-component construction rather than a completely new algorithm. Due to the formulation of the successive coordinate search method it can also be used to improve existing lattice rules. Numerical experiments show that the improvements of the SCS method are higher when the decay of the weights ${\bm{\gamma}}=\{\gamma_{j}\}_{j=1}^{d}$ is slow.

Acknowledgements We thank Peter Kritzer for some useful comments and discussions about the manuscript and we acknowledge financial support from the KU Leuven research fund (OT:3E130287 and C3:3E150478).