Function integration, reconstruction and approximation using rank-1 lattices

TL;DR

This paper develops new methods for function integration, reconstruction, and approximation using rank-1 lattices, leveraging transformations between periodic and non-periodic function spaces to improve efficiency and accuracy.

Contribution

It introduces novel theoretical insights and algorithms for rank-1 lattice-based function reconstruction in non-periodic settings, including bi-orthonormal basis methods and CBC construction strategies.

Findings

Fast discrete cosine transform enables efficient reconstruction.

New bounds and algorithms for CBC construction reduce auxiliary set size.

Methods improve accuracy and efficiency in non-periodic function approximation.

Abstract

We consider rank-1 lattices for integration and reconstruction of functions with series expansion supported on a finite index set. We explore the connection between the periodic Fourier space and the non-periodic cosine space and Chebyshev space, via tent transform and then cosine transform, to transfer known results from the periodic setting into new insights for the non-periodic settings. Fast discrete cosine transform can be applied for the reconstruction phase. To reduce the size of the auxiliary index set in the associated component-by-component (CBC) construction for the lattice generating vectors, we work with a bi-orthonormal set of basis functions, leading to three methods for function reconstruction in the non-periodic settings. We provide new theory and efficient algorithmic strategies for the CBC construction. We also interpret our results in the context of general function approximation and discrete least-squares approximation.

Tables

Table 1. Table 1: Summary of CBC algorithms for function reconstruction

Fourier space Cosine space and Chebyshev space Plan A Plan B Plan C (a) $𝒉\cdot 𝒛{\not\equiv }_{n}0$  for all $𝒉\cdot 𝒛{\not\equiv }_{n}0$  for all $𝒉\cdot 𝒛{\not\equiv }_{n}0$  for all NA $𝒉\in {\left(\Lambda \ominus \Lambda \right)}_s\setminus \left\{\mathrm{𝟎}\right\}$ $𝒉\in {\left(ℳ\left(\Lambda \right)\oplus ℳ\left(\Lambda \right)\right)}_s\setminus \left\{\mathrm{𝟎}\right\}$ $𝒉\in {\left(\Lambda \oplus ℳ\left(\Lambda \right)\right)}_s\setminus \left\{\mathrm{𝟎}\right\}$ (b) $𝒉\cdot 𝒛{\not\equiv }_n{𝒉}^{\prime }\cdot 𝒛$ ${𝝈}^{\prime }\left({𝒌}^{\prime }\right)\cdot 𝒛{\not\equiv }_n𝝈\left(𝒌\right)\cdot 𝒛$ $𝝈\left({𝒌}^{\prime }\right)\cdot 𝒛{\not\equiv }_n𝒌\cdot 𝒛$ $𝝈\left({𝒌}^{\prime }\right)\cdot 𝒛{\not\equiv }_n𝒌\cdot 𝒛$ for all $𝒉,{𝒉}^{\prime }\in {\Lambda }_s,$ for all $𝒌,{𝒌}^{\prime }\in {\Lambda }_s,$ for all $𝒌,{𝒌}^{\prime }\in {\Lambda }_s,$ for all $𝒌,{𝒌}^{\prime }\in {\Lambda }_s,$ $𝒉\ne {𝒉}^{\prime }$ $𝝈\in {𝒮}_{𝒌}$, ${𝝈}^{\prime }\in {𝒮}_{{𝒌}^{\prime }},$ $𝝈\in {𝒮}_{{𝒌}^{\prime }},$ $𝝈\in {𝒮}_{{𝒌}^{\prime }},$ ${𝝈}^{\prime }\left({𝒌}^{\prime }\right)\ne 𝝈\left(𝒌\right)$ $𝝈\left({𝒌}^{\prime }\right)\ne 𝒌$ $𝒌\ne {𝒌}^{\prime }$ (c) $n\sim \#{\left(\Lambda \ominus \Lambda \right)}_s$ $n\sim \#{\left(ℳ\left(\Lambda \right)\oplus ℳ\left(\Lambda \right)\right)}_s$ $n\sim \#{\left(\Lambda \oplus ℳ\left(\Lambda \right)\right)}_s$ $n\sim \#{\Lambda }_s\#ℳ\left({\Lambda }_s\right)$ (d) $\#{\left(\Lambda \ominus \Lambda \right)}_s$ $\#{\left(ℳ\left(\Lambda \right)\oplus ℳ\left(\Lambda \right)\right)}_s$ $\#{\left(\Lambda \oplus ℳ\left(\Lambda \right)\right)}_s$ NA (e) $n_{\mathrm{fail}}\#{\left(\Lambda \ominus \Lambda \right)}_s$ $n_{\mathrm{fail}}\#{\left(ℳ\left(\Lambda \right)\oplus ℳ\left(\Lambda \right)\right)}_s$ $n_{\mathrm{fail}}\#{\left(\Lambda \oplus ℳ\left(\Lambda \right)\right)}_s$ NA (f) ${\left(\#{\Lambda }_s\right)}^2$ ${\left(\#ℳ\left({\Lambda }_s\right)\right)}^2$ $\#{\Lambda }_s\#ℳ\left({\Lambda }_s\right)$ $\#{\Lambda }_s\#ℳ\left({\Lambda }_s\right)$ (g) $n_{\mathrm{fail}}{\left(\#{\Lambda }_s\right)}^2$ $n_{\mathrm{fail}}{\left(\#ℳ\left({\Lambda }_s\right)\right)}^2$ $n_{\mathrm{fail}}\#{\Lambda }_s\#ℳ\left({\Lambda }_s\right)$ (h) $n_{\mathrm{fail}}\#{\Lambda }_s$ $n_{\mathrm{fail}}\#ℳ\left({\Lambda }_s\right)$ $n_{\mathrm{fail}}\#ℳ\left({\Lambda }_s\right)$ (i) $n_{\mathrm{fail}}\sim \#{\Lambda }_s$ $n_{\mathrm{fail}}\sim \#ℳ\left({\Lambda }_s\right)$ $n_{\mathrm{fail}}\sim \#ℳ\left({\Lambda }_s\right)$ (j) $d{\left(\#\Lambda \right)}^2$ $d{\left(\#ℳ\left(\Lambda \right)\right)}^2$ $d\#\Lambda \#ℳ\left(\Lambda \right)$

(a) Standard formulation of the reconstruction condition at step $s$ with full/zero projection (b) Equivalent formulation of the reconstruction condition at step $s$ with full projection (c) Required size of $n$ to guarantee success at step $s$ (also need $n$ to cover spread of index set) (d) Cost of elimination approach at step $s$ based on (a) with full/zero projection (e) Cost of brute force approach at step $s$ based on (a) with full/zero projection (f) Cost of elimination approach at step $s$ based on (b) with full projection (g) Cost of brute force approach at step $s$ based on (b) with full projection (h) Cost of brute force approach at step $s$ based on (b) with full projection and smart lookup (i) Switching point on $n_{\mathrm{fail}}$ from brute force (h) to elimination (d)/(f) (j) Total cost of mixed CBC: brute force until $n_{\mathrm{fail}}$ reaches switching point then elimination

Equations

$$f=\sum _{\mathbf{k}}{\hat{f}}_{\mathbf{k}}{\alpha }_{\mathbf{k}}.$$

$$f=\sum _{\mathbf{k}\in \Lambda }{\hat{f}}_{\mathbf{k}}{\alpha }_{\mathbf{k}}.$$

$${\mathbf{t}}_i=\frac{i\mathbf{z}\mathrm{mod}n}{n}\in [0,1{]}^d\text{mbox}fori=0,\dots ,n-1.$$

$$Q_n(f):=\frac{1}{n}\sum _{i=0}^{n-1}f({\mathbf{t}}_i)$$

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

$${\mathcal{S}}_{\boldsymbol{k}}\,:=\,\left\{\boldsymbol{\sigma}\in\{\pm 1\}^{d}:\sigma_{j}=+1\mbox{ for each $j$ for which $k_{j}=0$}\right\}$$

$$\max (\Lambda ):=\underset{\mathbf{k}\in \Lambda }{\max }\underset{1\le j\le d}{\max }|k_j|,$$

$$\Lambda \oplus \Lambda$$

$$\Lambda \ominus \Lambda$$

$$\mathcal{M}(\Lambda )$$

$$\#(\Lambda \oplus \Lambda )\le (\#\Lambda {)}^2,\#(\Lambda \ominus \Lambda )\le (\#\Lambda {)}^2,\text{mbox}and\#\mathcal{M}(\Lambda )\le \sum _{\mathbf{k}\in \Lambda }{2}^{|\mathbf{k}{|}_0}\le 2^d\#\Lambda .$$

$$\sum _{\mathbf{k}\in \Lambda }\prod _{j=1}^{d}p(|k_j|)\le (\#\Lambda {)}^{\eta +1}.$$

$$\sum _{\mathbf{k}\in \Lambda }{2}^{|\mathbf{k}{|}_0}\le (\#\Lambda {)}^{\ln 3/\ln 2}.$$

$$\underset{\mathbf{k}\in \Lambda }{\max }{2}^{|\mathbf{k}{|}_0}\le \#\Lambda ,\sum _{\mathbf{k}\in \Lambda }{2}^{|\mathbf{k}{|}_0}\le (\#\Lambda {)}^{\ln 3/\ln 2},\text{mbox}and\#\mathcal{M}(\Lambda )\le \min \left(2^d\#\Lambda ,(\#\Lambda {)}^{\ln 3/\ln 2}\right).$$

$$\underset{1\le j\le d}{\max }\frac{k_j}{{\beta }_j},\sum _{1\le j\le d}\frac{k_j}{{\beta }_j},\text{mbox}or\prod _{j=1}^d\max \left(1,\frac{k_j}{{\beta }_j}\right),$$

$$\frac{\#{{\mathcal{M}}(\Lambda)}}{\#{\Lambda}}\,=\,\prod_{j=1}^{d}\left(1+\frac{\lfloor\beta_{j}m\rfloor}{1+\lfloor\beta_{j}m\rfloor}\right)\,\leq\,\min\bigg{(}2^{J},\exp\bigg{(}m\sum_{j=1}^{\infty}\beta_{j}\bigg{)}\bigg{)},$$

$$m+1\le \#\Lambda$$

$$2m+1\le \#\mathcal{M}(\Lambda )$$

$${\mathcal{F}}^{\rm Four}:=\Bigl{\{}f\in L^{2}\mathbin{\Big{|}}f:[0,1]^{d}\to\mathbb{C},\;f({\boldsymbol{x}})=\sum_{{\boldsymbol{h}}\in{\mathbb{Z}}^{d}}\widehat{f}_{\boldsymbol{h}}\,e^{2\pi{\mathrm{i}}\,{\boldsymbol{h}}\cdot{\boldsymbol{x}}}\text{ and }\sum_{{\boldsymbol{h}}\in{\mathbb{Z}}^{d}}|\widehat{f}_{\boldsymbol{h}}|<\infty\Bigr{\}},$$

$$⟨f_1,f_2⟩:={\int }_{[0,1{]}^d}{f}_1(\mathbf{x})\overline{f_2(\mathbf{x})}\mathrm{d}\mathbf{x}.$$

$$e_{\mathbf{h}}(\mathbf{x}):=e^{2\pi \mathrm{i}\mathbf{h}\cdot \mathbf{x}}$$

$${\hat{f}}_{\mathbf{h}}:=⟨f,e_{\mathbf{h}}⟩={\int }_{[0,1{]}^d}f(\mathbf{x})e^{-2\pi \mathrm{i}\mathbf{h}\cdot \mathbf{x}}\mathrm{d}\mathbf{x},\mathbf{h}\in {\mathbb{Z}}^d.$$

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

$$⟨f_1,f_2{⟩}_n:=Q_n(f_1\overline{f_2})$$

$$\text{for }f\in {\mathcal{F}}_{\Lambda }^{\mathrm{Four}}:f(\mathbf{x})=\sum _{\mathbf{h}\in \Lambda }{\hat{f}}_{\mathbf{h}}{e}^{2\pi \mathrm{i}\mathbf{h}\cdot \mathbf{x}}.$$

$$Q_n(e_{\mathbf{h}})=I(e_{\mathbf{h}})={\delta }_{\mathbf{h},0}\text{mbox}forall\mathbf{h}\in \Lambda ,$$

$$f^a(\mathbf{x})=\sum _{\mathbf{h}\in \Lambda }{\hat{f}}_{\mathbf{h}}^a{e}^{2\pi \mathrm{i}\mathbf{h}\cdot \mathbf{x}},$$

$${\hat{f}}_{\mathbf{h}}^a={\hat{f}}_{\mathbf{h}}\text{mbox}forall\mathbf{h}\in \Lambda \text{mbox}andf\in {\mathcal{F}}_{\Lambda }^{\mathrm{Four}},$$

$${\hat{f}}_{\mathbf{h}}^a$$

$$\displaystyle\,=\,Q_{n}\bigg{(}\bigg{(}\sum_{{\boldsymbol{h}}^{\prime}\in\Lambda}\widehat{f}_{{\boldsymbol{h}}^{\prime}}\,e_{{\boldsymbol{h}}^{\prime}}\bigg{)}e_{-{\boldsymbol{h}}}\bigg{)}\,=\,\sum_{{\boldsymbol{h}}^{\prime}\in\Lambda}\widehat{f}_{{\boldsymbol{h}}^{\prime}}\,Q_{n}(e_{{\boldsymbol{h}}^{\prime}-{\boldsymbol{h}}})=\widehat{f}_{\boldsymbol{h}}\quad\mbox{for all }{\boldsymbol{h}}\in\Lambda\mbox{ and }f\in{\mathcal{F}}^{\rm Four}_{\Lambda}.$$

$$Q_n(e_{{\mathbf{h}}^{\prime }-\mathbf{h}})=⟨e_{{\mathbf{h}}^{\prime }},e_{\mathbf{h}}⟩={\delta }_{\mathbf{h},{\mathbf{h}}^{\prime }}\text{mbox}forall\mathbf{h},{\mathbf{h}}^{\prime }\in \Lambda ,$$

$$Q_n(e_{\mathbf{h}})=I(e_{\mathbf{h}})={\delta }_{\mathbf{h},0}\text{mbox}forall\mathbf{h}\in \Lambda \ominus \Lambda ,$$

Full text

Keywords: Exact integration and approximation on finite index sets, Quasi-Monte Carlo methods, Rank-$1$ lattice points, Fourier space, Cosine space, Chebyshev space, Component-by-component construction.

AMS Subject classifications: 41A10 (Approximation by polynomials), 42A10 (Trigonometric approximation), 41A63 (Multidimensional problems), 42B05 (Fourier series and coefficients), 65D30 (Numerical integration), 65D32 (Quadrature and cubature formulas), 65D15 (Algorithms for functional approximation).

1 Introduction

In this paper we consider function integration, reconstruction and approximation in the periodic and non-periodic settings using rank-$1$ lattices. We explore the connection between three function space settings to transfer known results on rank-$1$ lattices from the periodic setting to the non-periodic settings. We obtain necessary and sufficient conditions on rank-$1$ lattices to achieve the exactness properties we require in each setting, and we develop efficient algorithms to construct the generating vectors for rank-$1$ lattices that satisfy these conditions.

More precisely, we consider functions with absolutely convergent series expansions with respect to an orthonormal basis, written in the generic form

[TABLE]

A large part of this paper is devoted to functions which are fully supported on a finite index set $\Lambda$, i.e.,

[TABLE]

We develop methods based on rank-$1$ lattices to exactly integrate such functions (2), and to exactly reconstruct all series coefficients $\widehat{f}_{\boldsymbol{k}}$ in (2). We also consider the approximation problem for functions (1) which are not finitely supported on $\Lambda$.

The three function space settings we consider are as follows:

• •

The Fourier space contains all absolutely convergent Fourier series in the unit cube $[0,1]^{d}$, with exponential basis functions $e_{\boldsymbol{h}}({\boldsymbol{x}})=e^{2\pi{\mathrm{i}}\,{\boldsymbol{h}}\cdot{\boldsymbol{x}}}$ and indices ${\boldsymbol{h}}\in{\mathbb{Z}}^{d}$.

• •

The cosine space contains all absolutely convergent cosine series in $[0,1]^{d}$, with half-period cosine basis functions $\phi_{\boldsymbol{k}}$ (see (10) below) and nonnegative indices ${\boldsymbol{k}}\in{\mathbb{N}}_{0}^{d}$.

• •

The Chebyshev space consists of all absolutely convergent Chebyshev series in the larger domain $[-1,1]^{d}$, under the Chebyshev measure, with Chebyshev basis functions $\eta_{\boldsymbol{k}}$ (see (27) below) and also nonnegative indices ${\boldsymbol{k}}\in{\mathbb{N}}_{0}^{d}$.

To avoid excessive notation we keep to generic notation for the three spaces wherever possible, including the same ‘hat’ notation for series coefficients. However, to effectively describe the connection between spaces, we often distinguish the basis functions $e_{\boldsymbol{h}}$, $\phi_{\boldsymbol{k}}$, $\eta_{\boldsymbol{k}}$, and we often use ${\boldsymbol{h}}$ and ${\boldsymbol{k}}$ to contrast indices containing integers ${\mathbb{Z}}$ or only nonnegative integers ${\mathbb{N}}_{0}$.

The Fourier space contains periodic functions while the cosine and Chebyshev spaces contain nonperiodic functions. The Fourier space is often referred to as the Wiener algebra; it is the standard setting for analyzing periodic functions, see, e.g., [11, 12, 13, 27]. The cosine space is connected to the Fourier space by the tent transform which is defined by $\varphi_{\rm tent}(x):=1-|2x-1|$ for $x\in[0,1]$ and is applied componentwise in $d$ dimensions, see, e.g., [5, 7, 29]. We show that the composition $\phi_{\boldsymbol{k}}\circ\varphi_{\rm tent}$ is the average over all of those exponential basis functions $e_{\boldsymbol{h}}$ for which $(|h_{1}|,\ldots,|h_{d}|)={\boldsymbol{k}}$ (see (14) below). Consequently, the tent-transformed cosine space is a subspace of the Fourier space. Thus we can apply results from the Fourier space to the cosine space via tent transform.

The Chebyshev space is related to the cosine space by the cosine transform, given by ${\boldsymbol{x}}=\cos(\pi{\boldsymbol{x}}^{\prime})\in[-1,1]^{d}$ for ${\boldsymbol{x}}^{\prime}\in[0,1]^{d}$, where the cosine function is applied componentwise, and we have $\eta_{\boldsymbol{k}}({\boldsymbol{x}})=\eta_{\boldsymbol{k}}(\cos(\pi{\boldsymbol{x}}^{\prime}))=\phi_{\boldsymbol{k}}({\boldsymbol{x}}^{\prime})$. Thus the cosine transform provides an isomorphism between the Chebyshev space and the cosine space. Trivially all results from the cosine space can be carried over to the Chebyshev space.

Rank-$1$ lattices have been well studied for integration, reconstruction and approximation in the Fourier space; see, e.g., [4, 28] for integration, [9, 10] for reconstruction, and [1, 11, 12, 13, 14, 15, 18, 19, 35, 36] for approximation. Given the generating vector ${\boldsymbol{z}}\in{\mathbb{Z}}^{d}$, the $n$ points of a rank-$1$ lattice are specified by

[TABLE]

For a Fourier space function $f$, the average of function values at the lattice points

[TABLE]

is known as a rank-$1$ lattice rule which is an equal-weight cubature rule for approximating the integral

[TABLE]

Rank-$1$ lattices have an important property known as the “character property” (see (9) below) which states that the cubature sum of the exponential basis functions $Q_{n}(e_{\boldsymbol{h}})$ can only take the value of $1$ or [math], depending on whether or not the dot product ${\boldsymbol{h}}\cdot{\boldsymbol{z}}$ is a multiple of $n$. Since the integral of the basis function $I(e_{\boldsymbol{h}})$ is $1$ if ${\boldsymbol{h}}=\boldsymbol{0}$ and is [math] otherwise, we easily deduce that a rank-$1$ lattice rule can exactly integrate a function $f=\sum_{{\boldsymbol{h}}\in\Lambda}\widehat{f}_{\boldsymbol{h}}\,e_{\boldsymbol{h}}$ whose Fourier series is supported on a finite set $\Lambda\subset{\mathbb{Z}}^{d}$ if and only if ${\boldsymbol{h}}\cdot{\boldsymbol{z}}$ is not a multiple of $n$ for all nonzero vectors ${\boldsymbol{h}}\in\Lambda$. This condition in turn leads to an efficient algorithm to construct a generating vector ${\boldsymbol{z}}$ with the exactness property in a component-by-component fashion. This result is stated later in Lemma 4, see also [4], and it can be said to be the starting point of all results in this paper. Indeed, the result extends to function reconstruction on $\Lambda$ where we evaluate all the Fourier coefficients $\widehat{f}_{\boldsymbol{h}}$ for ${\boldsymbol{h}}\in\Lambda$ by a rank-$1$ lattice rule, and the evaluations can be done using the fast Fourier transform. Using the character property one can deduce a necessary and sufficient condition when these Fourier coefficients can be recovered exactly, thus leading to a constructive algorithm to find suitable generating vectors by working with the “difference set” $\Lambda\ominus\Lambda$ which is obtained by forming all differences of indices in $\Lambda$. We state this result later in Lemma 5, which was first proved with varying generality in [9, 10, 27]. The idea has been further extended to the construction of “multiple rank-$1$ lattices” in [13], where the benefits of multiple reconstruction lattices are combined strategically to achieve the same goal with a reduced overall number of sampling nodes; we do not go down this path.

The connection between the Fourier space and the cosine space allows us to apply the theory of rank-$1$ lattices to the cosine space by tent transform. We can obtain necessary and sufficient conditions for tent-transformed rank-$1$ lattices to achieve the integral exactness and function reconstruction properties in the cosine space, see Lemmas 10 and 11 below (see also [29] for part of Lemma 11). In the case of function reconstruction on a finite index set $\Lambda\subset{\mathbb{N}}_{0}^{d}$, we end up having to work with quite a large auxiliary index set ${\mathcal{M}}(\Lambda)\oplus{\mathcal{M}}(\Lambda)$ in our component-by-component construction of the lattice generating vector, where ${\mathcal{M}}(\Lambda)$ denotes the “mirrored set” obtained from $\Lambda$ by including all sign changes of indices in $\Lambda$, while the $\oplus$ then means that we form the sum of all indices from ${\mathcal{M}}(\Lambda)$; we call this plan A. To improve on the computational efficiency of plan A, we show by working with a bi-orthonormal set of basis functions that we can achieve function reconstruction on $\Lambda$ with a weaker condition which means working with a smaller auxiliary index set $\Lambda\oplus{\mathcal{M}}(\Lambda)$, see Lemma 15 below; we call this plan B. We also relax the algorithm to not necessarily recover the normalization of the basis functions to arrive at plan C, which achieves the same reconstruction property at a lower computational cost, see Lemma 16 below. All three plans for function reconstruction in the cosine space can be computed using the fast discrete cosine transform.

The isomorphism between the cosine space and the Chebyshev space allows us to take all results from the cosine space to the Chebyshev space, including plans A, B, C. We arrive at tent-transformed and then cosine-transformed rank-$1$ lattices, which in the case of $n$ being even is also known as “Chebyshev lattices”, see, e.g., [6, 26], although we do not adopt this terminology. Our plan C for the Chebyshev space with even $n$ is essentially the approach in [26]; in this paper we do not require $n$ to be even.

1.1 Layout of the paper and highlight of new results

In Section 2 we review results on rank-$1$ lattices for integration and function reconstruction on a finite index set in the Fourier space, referencing essential results from [4, 9, 10, 27].

In Section 3 we introduce the cosine space and consider integration and function reconstruction, with three plans for achieving exact function reconstruction using rank-$1$ lattices with varying costs. Except for the if-part of Lemma 11 and Lemma 12 which was proved in [29], all remaining results in this section are new, including plan B and plan C for function reconstruction and the applicability of fast discrete cosine transform.

In Section 4 we present the corresponding results for the Chebyshev space. Lemma 21 for plan C with even $n$ turns out to be equivalent to the approach in [26]. However, the precise connection to the cosine space via cosine transform and in turn the precise connection to the Fourier space via tent transform are both new interpretations here, and they lead to broader implications in the Chebyshev space. In particular, the multiplicity of the transformed points under these interpretations are known explicitly for $n$ both even and odd, and fast discrete cosine transform can be applied for all $n$.

Section 5 is devoted to the theory and algorithmic aspect of the component-by-component (CBC) construction for lattice generating vectors achieving various conditions needed for the exactness properties. As the theoretical justification for the CBC construction, Theorem 23 generalises previous results proved in [4, 9, 10] and provides a cheaper variant of the algorithm when building up the index set, while Theorem 26 is new and specific to plan C. The systematic way to combine two different approaches (namely, the “brute force” approach and the “elimination” approach, to be explained in Section 5) in a mixed CBC construction is new. Strategies for storage and a “smart lookup” to efficiently search through difference and/or sum involving mirrored sets are also new.

Finally in Section 6 we interpret our results in the context of approximation of general functions that are not necessarily supported on a finite index set, and compare them with discrete least-squares approximation as analysed in [2, 3, 21, 22, 25]. We mention other known results in function approximation based on rank-$1$ lattices (see, e.g., [16, 17, 33, 34] for general results and [1, 5, 11, 13, 14, 15, 18, 19, 35, 36] for rank-$1$ lattices).

We end the introduction with setting the notation on multiindices and introducing some special index sets.

1.2 Notation on multiindices and special index sets

Throughout this paper we use $\#{}$ to denote the cardinality of a set. For any multiindex ${\boldsymbol{k}}\in{\mathbb{Z}}^{d}$, we write $|{\boldsymbol{k}}|_{0}:=\#{\{1\leq j\leq\ d:k_{j}\neq 0\}}$ for the number of nonzero indices in ${\boldsymbol{k}}$. For ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in{\mathbb{Z}}^{d}$, ${\boldsymbol{k}}^{\prime}\leq{\boldsymbol{k}}$ is to be interpreted componentwise, i.e., $k_{j}^{\prime}\leq k_{j}$ for all $j$.

For $\boldsymbol{\sigma}\in\{\pm 1\}^{d}$ and ${\boldsymbol{k}}\in{\mathbb{Z}}^{d}$, we write $\boldsymbol{\sigma}({\boldsymbol{k}}):=(\sigma_{1}k_{1},\ldots,\sigma_{d}k_{d})$ to mean that we apply the sign changes in $\boldsymbol{\sigma}$ componentwise to ${\boldsymbol{k}}$. For any ${\boldsymbol{k}}\in{\mathbb{Z}}^{d}$ we use

[TABLE]

to denote a set of unique sign changes for ${\boldsymbol{k}}$. Then clearly we have $\#{{\mathcal{S}}_{\boldsymbol{k}}}=2^{|{\boldsymbol{k}}|_{0}}$.

We will consider index sets with some special properties:

• •

An index set $\Lambda\subset{\mathbb{N}}_{0}^{d}$ is downward closed if ${\boldsymbol{k}}^{\prime}\in\Lambda$ whenever ${\boldsymbol{k}}^{\prime}\leq{\boldsymbol{k}}$ and ${\boldsymbol{k}}\in\Lambda$. This means that from every ${\boldsymbol{k}}\in\Lambda$ we can move towards $\boldsymbol{0}$ along the coordinate axes without finding a ${\boldsymbol{k}}^{\prime}\not\in\Lambda$. Analogous definition holds with ${\mathbb{N}}_{0}^{d}$ replaced by ${\mathbb{Z}}^{d}$.

• •

An index set $\Lambda\subset{\mathbb{Z}}^{d}$ is centrally symmetric if $-{\boldsymbol{k}}\in\Lambda$ whenever ${\boldsymbol{k}}\in\Lambda$.

• •

An index set $\Lambda\subset{\mathbb{Z}}^{d}$ is fully sign symmetric if $\boldsymbol{\sigma}({\boldsymbol{k}})\in\Lambda$ whenever ${\boldsymbol{k}}\in\Lambda$ and $\boldsymbol{\sigma}\in\{\pm 1\}^{d}$.

• •

An index set $\Lambda\subset{\mathbb{N}}_{0}^{d}$ is an (anisotropic) tensor product set if there exist ${\boldsymbol{a}},{\boldsymbol{b}}\in{\mathbb{N}}_{0}^{d}$ such that $\Lambda=\{{\boldsymbol{k}}\in{\mathbb{N}}_{0}^{d}:{\boldsymbol{a}}\leq{\boldsymbol{k}}\leq{\boldsymbol{b}}\}$. Analogous definition holds with ${\mathbb{N}}_{0}^{d}$ replaced by ${\mathbb{Z}}^{d}$.

For any index set $\Lambda\subset{\mathbb{N}}_{0}^{d}$ or $\Lambda\subset{\mathbb{Z}}^{d}$, we denote its largest component in magnitude by

[TABLE]

and we define

[TABLE]

If $\Lambda\subset{\mathbb{Z}}^{d}$ is centrally symmetric, then $\Lambda\oplus\Lambda=\Lambda\ominus\Lambda$. If $\Lambda\in{\mathbb{Z}}^{d}$ is fully sign symmetric, then ${\mathcal{M}}(\Lambda)=\Lambda$. For $\Lambda\subset\mathbb{Z}^{d}$, $\Lambda\ominus\Lambda$ is always centrally symmetric (since both ${\boldsymbol{k}}-{\boldsymbol{k}}^{\prime}$ and ${\boldsymbol{k}}^{\prime}-{\boldsymbol{k}}$ belong to $\Lambda\ominus\Lambda$ when ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in\Lambda$).

Trivially we have

[TABLE]

The squaring effect in the upper bounds for the sum/difference set cannot be avoided in general, since even if all multiindices in $\Lambda$ fall on the first two axes (i.e., all components of ${\boldsymbol{k}}$ are zero except for one of $k_{1}$ and $k_{2}$), the sum/difference sets will contain a large rectangle (so there is a lower bound of the same order). On the other hand, the $2^{d}$ factor in the upper bound for ${\mathcal{M}}(\Lambda)$ can sometimes be improved, as shown in the forthcoming Lemma 2 and Example 3 below. We also need the next Lemma 1, whose proof uses induction arguments from [23, 2].

Lemma 1**.**

In any dimension $d$, given any $\Lambda\subset\mathbb{Z}^{d}$ downward closed and any polynomial $p(n)=\sum_{k=0}^{\eta}b_{k}n^{k}$ of degree $\eta\geq 0$ with nonnegative coefficients $b_{0}\leq 1$ and all $b_{k}\leq\binom{\eta+1}{k}$, it holds

[TABLE]

Moreover, it holds

[TABLE]

Proof.

When $\eta=0$ the result holds true. Consider then the case $\eta\geq 1$. Every downward closed set $\Lambda\subset\mathbb{Z}^{d}$ can be seen as a set $\widetilde{\Lambda}\subset\mathbb{N}_{0}^{2d}$ constructed in the following way. Start with $\widetilde{\Lambda}=\emptyset$. For any ${\boldsymbol{k}}\in\Lambda$, define the sets $C_{\boldsymbol{k}}:=\{1\leq j\leq d:k_{j}<0\}\subseteq\{1,\ldots,d\}$ and $U_{\boldsymbol{k}}:=\{1,\ldots,d\}\setminus C_{\boldsymbol{k}}$. Then define $\widetilde{{\boldsymbol{k}}}\in\mathbb{N}_{0}^{2d}$ by setting $\widetilde{k}_{j}=k_{j}$ and $\widetilde{k}_{j+d}=0$ for all $j\in U_{\boldsymbol{k}}$, and $\widetilde{k}_{j}=0$ and $\widetilde{k}_{j+d}=-k_{j}$ for all $j\in C_{\boldsymbol{k}}$. Finally add $\widetilde{{\boldsymbol{k}}}$ to $\widetilde{\Lambda}$. Notice that this algorithm establishes a one-to-one correspondence between the elements of $\Lambda$ and $\widetilde{\Lambda}$, and therefore $\#\Lambda=\#\widetilde{\Lambda}$. By construction $\widetilde{\Lambda}$ is also downward closed in $\mathbb{N}_{0}^{2d}$. Applying Theorem 1 from [23] to the set $\widetilde{\Lambda}\in\mathbb{N}_{0}^{2d}$ we obtain (3).

For the proof of (4), as above, starting from $\Lambda$ we construct the downward closed set $\widetilde{\Lambda}\subset\mathbb{N}_{0}^{2d}$ such that $\#\widetilde{\Lambda}=\#\Lambda$, and then apply Lemma 3.3 from [2] to the set $\widetilde{\Lambda}$. $\Box$

Lemma 2**.**

If $\Lambda\subset{\mathbb{Z}}^{d}$ is downward closed then

[TABLE]

Proof.

For the first bound, since $\Lambda$ is downward closed, for any ${\boldsymbol{k}}\in\Lambda$, the set $\Lambda$ will include the hyper-rectangle with ${\boldsymbol{k}}$ and the origin as corners. Thus $\#{\Lambda}\geq\prod_{1\leq j\leq d,\,k_{j}\neq 0}(1+|k_{j}|)\geq 2^{|{\boldsymbol{k}}|_{0}}$. The second bound is proved in Lemma 1. The third bound is an immediate consequence of the second bound. $\Box$

Since $\#{\Lambda}$ most likely grows with $d$, in general it is not obvious which of $2^{d}\,\#{\Lambda}$ or $(\#{\Lambda})^{\ln 3/\ln 2}$ is a better bound for $\#{{\mathcal{M}}(\Lambda)}$. If $\#{\Lambda}$ can be bounded independently of $d$, then most likely so can $\#{{\mathcal{M}}(\Lambda)}$.

Example 3**.**

Consider a “weighted” index set of “degree” $m\in{\mathbb{N}}$ defined by (see, e.g., [4]) $\Lambda\,=\,\{{\boldsymbol{k}}\in{\mathbb{N}}_{0}^{d}:r({\boldsymbol{k}})\leq m\}$, where $r({\boldsymbol{k}})$ is given by

[TABLE]

with $1=\beta_{1}\geq\beta_{2}\geq\cdots>0$ and $\sum_{j=1}^{\infty}\beta_{j}<\infty$. The first example is an anisotropic tensor product set and is the largest of the three examples. We have $\#{\Lambda}=\prod_{j=1}^{d}(1+\lfloor\beta_{j}m\rfloor)\leq\exp(m\sum_{j=1}^{\infty}\beta_{j})$ and $\#{{\mathcal{M}}(\Lambda)}=\prod_{j=1}^{d}(1+2\lfloor\beta_{j}m\rfloor)\leq\exp(2m\sum_{j=1}^{\infty}\beta_{j})$, so both are bounded independently of $d$. Their ratio satisfies

[TABLE]

where $J$ is the “truncation dimension” such that $\beta_{j}m<1$ for all $j\geq J$. If we have $\lambda\in(0,1]$ such that $\sum_{j=1}^{\infty}\beta_{j}^{\lambda}<\infty$, then $\beta_{J}<J^{-1/\lambda}(\sum_{j=1}^{\infty}\beta_{j}^{\lambda})^{1/\lambda}$ and it suffices to take $J=m^{\lambda}(\sum_{j=1}^{\infty}\beta_{j}^{\lambda})$. Both upper bounds on the ratio grow exponentially with $m$.

The third example in (5) is the smallest of the three. Its mirror set is commonly referred to as the “Zaremba cross” or “hyperbolic cross”, see, e.g., [4, 9]. For all $\tau>1$ we have

[TABLE]

where $\zeta(\tau):=\sum_{k=1}^{\infty}k^{-\tau}$ is the Riemann zeta function. Since $\tau$ can be arbitrarily close to $1$, $\#{{\mathcal{M}}(\Lambda)}$ is essentially of the same order as $\#{\Lambda}$, both are bounded independently of $d$. The upper bound on $\#{{\mathcal{M}}(\Lambda)}$ is proved in [14].

2 Periodic setting based on trigonometric polynomials

2.1 Fourier series

We start by considering periodic functions on $[0,1]^{d}$. Let ${\mathcal{F}}^{\rm Four}$ denote the space of complex-valued functions defined on $[0,1]^{d}$ with absolutely converging Fourier series:

[TABLE]

where ${\boldsymbol{h}}\cdot{\boldsymbol{x}}:=h_{1}x_{1}+\cdots+h_{d}x_{d}$ is the usual dot product and $\widehat{f}_{\boldsymbol{h}}$ are the Fourier coefficients. We equip ${\mathcal{F}}^{\rm Four}$ with the usual $L^{2}$ inner product

[TABLE]

The exponential functions form an orthonormal basis

[TABLE]

satisfying $\langle e_{\boldsymbol{h}},e_{{\boldsymbol{h}}^{\prime}}\rangle=\delta_{{\boldsymbol{h}},{\boldsymbol{h}}^{\prime}}$, where the Kronecker delta function yields $1$ if ${\boldsymbol{h}}={\boldsymbol{h}}^{\prime}$ and [math] if ${\boldsymbol{h}}\neq{\boldsymbol{h}}^{\prime}$. The Fourier coefficients are given by

[TABLE]

The norm of $f$ satisfies $\|f\|^{2}\,=\,\int_{[0,1]^{d}}|f({\boldsymbol{x}})|^{2}\,\mathrm{d}{\boldsymbol{x}}\,=\,\sum_{{\boldsymbol{h}}\in{\mathbb{Z}}^{d}}|\widehat{f}_{\boldsymbol{h}}|^{2}$.

2.2 Fourier coefficients by cubature

For $f\in{\mathcal{F}}^{\rm Four}$ we define the integral operator

[TABLE]

Later we will seek a cubature formula $Q_{n}(f)$ which uses linear combinations of $n$ evaluations of $f$ to approximate $I(f)$. We define a discrete inner product

[TABLE]

as an approximation to (6).

Given an arbitrary finite index set $\Lambda\subset{\mathbb{Z}}^{d}$, we consider the subspace ${\mathcal{F}}^{\rm Four}_{\Lambda}$ of all functions whose Fourier series is supported solely on $\Lambda$, i.e.,

[TABLE]

Implicitly, this means that all other Fourier coefficients of $f$ are zero, i.e., $\widehat{f}_{\boldsymbol{h}}$ = 0 for ${\boldsymbol{h}}\notin\Lambda$.

In this paper we will demand one or both of the following related properties on the cubature formula:

• •

Integral exactness. We want our cubature formula to be exact for all functions which are supported solely on $\Lambda$, i.e., we want $Q_{n}(f)=I(f)$ for all $f\in{\mathcal{F}}^{\rm Four}_{\Lambda}$. This holds if and only if

[TABLE]

i.e., our cubature formula integrates exactly all basis functions $e^{2\pi{\mathrm{i}}\,{\boldsymbol{h}}\cdot{\boldsymbol{x}}}$ with ${\boldsymbol{h}}\in\Lambda$.

• •

Function reconstruction. Instead of (7) we consider

[TABLE]

where each Fourier coefficient $\widehat{f}_{\boldsymbol{h}}=\langle f,e_{\boldsymbol{h}}\rangle=I(f\,e_{-{\boldsymbol{h}}})$ in (7) is replaced by the cubature formula $\widehat{f}_{\boldsymbol{h}}^{a}:=\langle f,e_{\boldsymbol{h}}\rangle_{n}=Q_{n}(f\,e_{-{\boldsymbol{h}}})$. We demand the “non-aliasing” condition that

[TABLE]

so that $f^{a}$ is a reconstruction of $f$. (If other coefficients $\widehat{f}_{{\boldsymbol{h}}^{\prime}}$ with ${\boldsymbol{h}}^{\prime}\neq{\boldsymbol{h}}$ contribute to $\widehat{f}_{\boldsymbol{h}}^{a}$ then this is called “aliasing”.) Using the linearity of $Q_{n}$, we then have

[TABLE]

This holds if and only if

[TABLE]

which is equivalent to

[TABLE]

i.e., our cubature formula integrates exactly all basis functions $e^{2\pi{\mathrm{i}}\,{\boldsymbol{h}}\cdot{\boldsymbol{x}}}$ with ${\boldsymbol{h}}\in\Lambda\ominus\Lambda$.

2.3 Rank-1 lattice rules

Consider now the cubature formula given by rank-$1$ lattices

[TABLE]

where ${\boldsymbol{z}}\in{\mathbb{Z}}^{d}$ is an integer vector known as the generating vector. It is easy to verify the “character property” that for any ${\boldsymbol{h}}\in{\mathbb{Z}}^{d}$,

[TABLE]

where the notation $a\equiv_{n}b$ means that $(a\bmod n)=(b\bmod n)$. This leads to the well-known lattice cubature error formula for $f\in{\mathcal{F}}^{\rm Four}$

[TABLE]

The set of integer vectors $\{{\boldsymbol{h}}\in{\mathbb{Z}}^{d}:{\boldsymbol{h}}\cdot{\boldsymbol{z}}\equiv_{n}0\}$ is known as the dual lattice. Clearly the cubature rule is exact for a function $f$ solely supported on $\Lambda$ if and only if the dual lattice does not contain any index from $\Lambda\setminus\{\boldsymbol{0}\}$. We know how to obtain such a lattice rule generating vector using a component-by-component construction.

Lemma 4** (Integral exactness).**

Let $\Lambda\subset{\mathbb{Z}}^{d}$ be an arbitrary index set. A lattice rule with $n$ points and generating vector ${\boldsymbol{z}}$ integrates exactly all functions $f\in{\mathcal{F}}^{\rm Four}_{\Lambda}$ solely supported on $\Lambda$ if and only if

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

with $\kappa=2$ if $\Lambda$ is centrally symmetric and $\kappa=1$ otherwise.

Proof.

The result for some standard anisotropic, downward closed and centrally symmetric sets $\Lambda$ can be found in Cools, Kuo & Nuyens [4]. A proof for general index sets is provided later in Section 5, see Theorem 23 and Remark 24. A similar proof can be found in Kämmerer [10]. $\Box$

Lemma 5** (Function reconstruction).**

Let $\Lambda\subset{\mathbb{Z}}^{d}$ be an arbitrary index set. A lattice rule $Q_{n}$ with $n$ points and generating vector ${\boldsymbol{z}}$ reconstructs exactly the Fourier coefficients of all functions $f\in{\mathcal{F}}^{\rm Four}_{\Lambda}$ solely supported on $\Lambda$, by

[TABLE]

if and only if

[TABLE]

which is equivalent to

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Proof.

The result follows directly from Theorem 23, noting that the difference set $\Lambda\ominus\Lambda$ is centrally symmetric and contains $\boldsymbol{0}$, and therefore $\frac{1}{2}\#((\Lambda\ominus\Lambda)\setminus\{\boldsymbol{0}\})+1=\frac{1}{2}(\#(\Lambda\ominus\Lambda)+1)$. Alternatively, the result for $\Lambda$ a hyperbolic cross index set can be found in Kämmerer [9], while the result for any arbitrary index set $\Lambda$ can be found in Kämmerer [10] and Potts & Volkmer [27, Theorem 2.1]. $\Box$

We end this section by the very interesting property that mapping from function values to Fourier coefficients and the other way around can be done using a one-dimensional fast Fourier transform.

Lemma 6**.**

Let ${\boldsymbol{z}}$ be a generating vector for an $n$-point rank-$1$ lattice satisfying the reconstruction property on an arbitrary index set $\Lambda\subset\mathbb{Z}^{d}$ according to Lemma 5. For a function $f\in{\mathcal{F}}^{\rm Four}_{\Lambda}$ solely supported on $\Lambda$ we can compute

[TABLE]

where ${\boldsymbol{f}}\in\mathbb{C}^{n}$ is a vector containing function values and ${\boldsymbol{F}}\in\mathbb{C}^{n}$ is a vector containing Fourier coefficients. Here $\mathrm{FFT}$ and $\mathrm{IFFT}$ are the one-dimensional fast Fourier transform and its inverse, respectively, with a normalization $1/n$ for $\mathrm{FFT}$ and $1$ for $\mathrm{IFFT}$; both mappings have cost ${\mathcal{O}}(n\log(n))$.

Proof.

This follows from expanding the formula for $\widehat{f}_{\boldsymbol{h}}^{a}$ in Lemma 5. Each ${\boldsymbol{h}}\in\Lambda$ will correspond to a unique value of ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\bmod{n}$ by the non-aliasing condition in Lemma 5. The other Fourier coefficients are zero by the assumption that $f\in{\mathcal{F}}^{\rm Four}_{\Lambda}$ is solely supported on $\Lambda$. $\Box$

Remark 7**.**

If the function $f$ has wider support in the Fourier space than just $\Lambda$, then the vector ${\boldsymbol{F}}$ resulting from the evaluation ${\boldsymbol{F}}=\mathrm{FFT}({\boldsymbol{f}})$ will not necessarily be zero at positions $F_{\kappa}$ when $\kappa$ does not correspond to a value of ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\bmod{n}$ for some ${\boldsymbol{h}}\in\Lambda$. This is due to the aliasing effect from ${\boldsymbol{h}}$ outside of $\Lambda$ and this will also contaminate all other components of ${\boldsymbol{F}}$. It is possible to extend the index set to full size $n$ while still keeping the reconstruction property on the extended index set such that all values in ${\boldsymbol{F}}$ can be interpreted as Fourier coefficients. This technique has been used, e.g., in [19, 24, 30, 31].

3 Nonperiodic setting based on half-period cosines

3.1 Cosine series

The cosine basis functions are a complete and orthonormal basis for $L^{2}([0,1]^{d})$:

[TABLE]

where $|{\boldsymbol{k}}|_{0}$ denotes the count of the nonzero entries in the vector ${\boldsymbol{k}}$, and we have $\langle\phi_{\boldsymbol{k}},\phi_{{\boldsymbol{k}}^{\prime}}\rangle=\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$. The “cosine space” ${\mathcal{F}}^{\rm cos}$ consists of nonperiodic real-valued functions on $[0,1]^{d}$ with absolutely converging cosine series:

[TABLE]

where the cosine coefficients are

[TABLE]

This space was studied for integration and approximation in [7, 29, 5]. Even though the cosine basis is a complete orthonormal system for $L^{2}([0,1]^{d})$, it does not allow the representation of arbitrary polynomials.

3.2 Cosine coefficients by cubature

As in Section 2, for a given finite index set $\Lambda\subset{\mathbb{N}}_{0}^{d}$ we consider the subspace ${\mathcal{F}}^{\rm cos}_{\Lambda}$ of all functions whose cosine series is supported solely on $\Lambda$, i.e.,

[TABLE]

and we are interested in two related properties on the cubature formula:

• •

Integral exactness. We want $Q_{n}(f)=I(f)$ for all $f\in{\mathcal{F}}^{\rm cos}_{\Lambda}$, which holds if and only if

[TABLE]

• •

Function reconstruction. We replace each cosine coefficient $\widehat{f}_{\boldsymbol{k}}=\langle f,\phi_{\boldsymbol{k}}\rangle=I(f\,\phi_{\boldsymbol{k}})$ in (11) by the cubature formula $\widehat{f}_{\boldsymbol{k}}^{a}:=\langle f,\phi_{\boldsymbol{k}}\rangle_{n}=Q_{n}(f\,\phi_{\boldsymbol{k}})$, and demand the non-aliasing condition

[TABLE]

which holds if and only if

[TABLE]

Unlike the Fourier case where a product of two basis functions is another basis function, here the condition (13) is not straightforward to simplify, except when the index set $\Lambda$ is downward closed. In the next section we will obtain necessary and sufficient conditions for function reconstruction by connecting with the Fourier space, without the restriction to downward closed index sets.

3.3 Connection with the Fourier case via tent transform

Below we will obtain sufficient conditions to achieve (12) and (13) in the cosine space by utilizing a known connection with the Fourier case via the so-called “tent transform”(see, e.g., [8])

[TABLE]

The tent transform is a Lebesgue preserving transformation and therefore a componentwise mapping of

[TABLE]

yields

[TABLE]

To get a sense of how this transformation works, it is informative to consider the univariate case:

[TABLE]

In the following, we recall the definition of the “mirrored” index set associated with the index set $\Lambda$, ${\mathcal{M}}(\Lambda):=\{\boldsymbol{\sigma}({\boldsymbol{k}}):{\boldsymbol{k}}\in\Lambda,\boldsymbol{\sigma}\in\{\pm 1\}^{d}\}=\bigcup_{{\boldsymbol{k}}\in\Lambda}\{\boldsymbol{\sigma}({\boldsymbol{k}}):\boldsymbol{\sigma}\in{\mathcal{S}}_{\boldsymbol{k}}\}$, where ${\mathcal{S}}_{\boldsymbol{k}}$ is the set of all unique sign changes of ${\boldsymbol{k}}$.

Lemma 8** (Integral exactness – sufficiency).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. If a cubature rule $Q_{n}^{*}(f)=\sum_{i=0}^{n-1}w_{i}^{*}\,f({\boldsymbol{t}}_{i}^{*})$ integrates exactly all Fourier basis functions $e_{\boldsymbol{h}}$ with ${\boldsymbol{h}}\in{\mathcal{M}}(\Lambda)$, then the cubature rule $Q_{n}(f)=\sum_{i=0}^{n-1}w_{i}\,f({\boldsymbol{t}}_{i})$ with $w_{i}=w_{i}^{*}$ and ${\boldsymbol{t}}_{i}=\varphi_{\rm tent}({\boldsymbol{t}}_{i}^{*})$ integrates exactly all cosine space functions $f\in{\mathcal{F}}^{\rm cos}_{\Lambda}$ solely supported on $\Lambda$.

Proof.

For any ${\boldsymbol{k}}\in{\mathbb{N}}_{0}^{d}$ we can write

[TABLE]

which follows from expanding the product of $\cos(\theta_{j})=(e^{{\mathrm{i}}\theta_{j}}+e^{-{\mathrm{i}}\theta_{j}})/2$ for those $\theta_{j}\neq 0$. Furthermore, since $\cos(\pi k\,\varphi_{\rm tent}(x))=\cos(2\pi kx)$ for all $k\in{\mathbb{N}}_{0}$, we also have

[TABLE]

Thus if we have a cubature rule $Q^{*}_{n}$ which integrates exactly all Fourier basis functions $e_{\boldsymbol{\sigma}({\boldsymbol{k}})}$ for all sign changes of ${\boldsymbol{k}}\in\Lambda$, then

[TABLE]

as required for integral exactness in (12). The cubature rule $Q_{n}$ is obtained from $Q_{n}^{*}$ by applying the tent-transform to the points. $\Box$

Lemma 9** (Function reconstruction – sufficiency).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. If a cubature rule $Q_{n}^{*}(f)=\sum_{i=0}^{n-1}w_{i}^{*}\,f({\boldsymbol{t}}_{i}^{*})$ integrates exactly all Fourier basis functions $e_{\boldsymbol{h}}$ with ${\boldsymbol{h}}\in{\mathcal{M}}(\Lambda)\oplus{\mathcal{M}}(\Lambda)$, then the cubature rule $Q_{n}(f)=\sum_{i=0}^{n-1}w_{i}\,f({\boldsymbol{t}}_{i})$ with $w_{i}=w_{i}^{*}$ and ${\boldsymbol{t}}_{i}=\varphi_{\rm tent}({\boldsymbol{t}}_{i}^{*})$ reconstructs exactly the cosine coefficients of all cosine space functions $f\in{\mathcal{F}}^{\rm cos}_{\Lambda}$ solely supported on $\Lambda$.

Proof.

For any ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in{\mathbb{N}}_{0}^{d}$ we have from (14) that

[TABLE]

Thus if we have a cubature rule $Q^{*}_{n}$ which integrates exactly all Fourier basis functions $e_{\boldsymbol{\sigma}({\boldsymbol{k}})+\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}^{\prime})}$ for all sign changes of ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in\Lambda$, then

[TABLE]

which is the reconstruction property (13). In the penultimate step we used the property that $\boldsymbol{\sigma}({\boldsymbol{k}})+\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}^{\prime})=\boldsymbol{0}$ if and only if ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime}$ and $\sigma_{j}^{\prime}=-\sigma_{j}$ whenever $k_{j}\neq 0$. $\Box$

Now we consider the situation where the cubature rule $Q_{n}^{*}$ in Lemma 8 and Lemma 9 is a rank-1 lattice rule (8). In this case, the corresponding cubature rule $Q_{n}(f)=Q_{n}^{*}(f\circ\varphi_{\rm tent})$ is often called a tent-transformed lattice rule, given explicitly by

[TABLE]

The character property (9) of lattice rules enables us to conclude that the implications in Lemma 8 and Lemma 9 also hold in the opposite direction, and we obtain necessary and sufficient conditions for tent-transformed lattice rules to achieve our desired properties. Lemma 10 and the “only if” part of Lemma 11 have not been explicitly stated in the literature.

Lemma 10** (Integral exactness).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. A tent-transformed lattice rule $Q_{n}$ of a lattice rule $Q_{n}^{*}$ with $n$ points and generating vector ${\boldsymbol{z}}$ integrates exactly all cosine space functions $f\in{\mathcal{F}}^{\rm cos}_{\Lambda}$ solely supported on $\Lambda$ if and only if

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Proof.

The “if” direction follows by combining Lemma 4 with Lemma 8. To prove the “only if” direction, we observe from the character property (9) that the terms $Q_{n}^{*}(e_{\boldsymbol{\sigma}({\boldsymbol{k}})})$ on the right-hand side of (15) can only take the values of $1$ or [math] so there can be no cancelation. In particular, when ${\boldsymbol{k}}\neq\boldsymbol{0}$, if $Q_{n}(\phi_{\boldsymbol{k}})$ is [math] on the left-hand side of (15) then necessarily all terms $Q_{n}^{*}(e_{\boldsymbol{\sigma}({\boldsymbol{k}})})$ are [math] on the right-hand side of (15), which implies $\boldsymbol{\sigma}({\boldsymbol{k}})\cdot{\boldsymbol{z}}\not\equiv_{n}0$. When ${\boldsymbol{k}}=\boldsymbol{0}$ both sides of (15) are equal to $1$, and trivially $Q_{n}(\phi_{\boldsymbol{0}})=1$ implies $Q_{n}^{*}(e_{\boldsymbol{0}})=1$. The CBC result follows from Theorem 23, noting that ${\mathcal{M}}(\Lambda)$ is centrally symmetric. $\Box$

Lemma 11** (Function reconstruction – plan A).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. A tent-transformed lattice rule $Q_{n}$ of a lattice rule $Q_{n}^{*}$ with $n$ points and generating vector ${\boldsymbol{z}}$ reconstructs exactly the cosine coefficients of all cosine space functions $f\in{\mathcal{F}}^{\rm cos}_{\Lambda}$ solely supported on $\Lambda$, by

[TABLE]

if and only if

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Proof.

As in the previous proof, the “if” direction follows by combining Lemma 5 with Lemma 9. When ${\boldsymbol{k}}\neq{\boldsymbol{k}}^{\prime}$, if $Q_{n}(\phi_{\boldsymbol{k}}\,\phi_{{\boldsymbol{k}}^{\prime}})=0$ on the left-hand side of (16), then necessarily all terms $Q_{n}^{*}(e_{\boldsymbol{\sigma}({\boldsymbol{k}})+\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}^{\prime})})$ are [math] on the right-hand side of (16), since the only permissible values are [math] or $1$ due to the character property (9). When ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime}$, if $Q_{n}(\phi_{\boldsymbol{k}}\,\phi_{\boldsymbol{k}})=1$ on the left-hand side of (16), then we have

[TABLE]

Necessarily, all terms $Q_{n}^{*}(e_{\boldsymbol{\sigma}({\boldsymbol{k}})+\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}})})$ must be zero for $\boldsymbol{\sigma}({\boldsymbol{k}})+\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}})\neq\boldsymbol{0}$. Hence we conclude that $(\boldsymbol{\sigma}({\boldsymbol{k}})+\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}))\cdot{\boldsymbol{z}}\not\equiv_{n}0$ for all ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime},\boldsymbol{\sigma},\boldsymbol{\sigma}^{\prime}$ satisfying $\boldsymbol{\sigma}({\boldsymbol{k}})+\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}^{\prime})\neq\boldsymbol{0}$. The CBC result follows from Theorem 23. $\Box$

The tent transformation “stretches and folds the domain” so that essentially one half of the lattice points will land on top of the other half. This is given precisely by the property that

[TABLE]

There is one point at the origin ${\boldsymbol{t}}_{0}$ which will not be duplicated. When $n$ is even, there is one other point $\varphi_{\rm tent}({\boldsymbol{t}}_{n/2})$ which will not be duplicated. All other points have a multiplicity of two provided that the generating vector includes at least one component $z_{j}$ such that $\gcd(n,z_{j})=1$. Typically in a CBC construction we set $z_{1}=1$. This is sufficient to ensure uniqueness.

Lemma 12**.**

A rank-$1$ lattice rule with $n$ points and generating vector ${\boldsymbol{z}}\in\mathbb{Z}_{n}^{d}$, where $\gcd(n,z_{j})=1$ for some $j$, has $\lfloor n/2+1\rfloor$ unique points after tent transform.

Proof.

See Suryanarayana, Nuyens & Cools [29]. $\Box$

3.4 Alternative approach for function reconstruction

In this subsection we use an alternative approach for function reconstruction. There are two essential ingredients, which we will separate into plan B and plan C below. Firstly, we make use of the tent transform and bi-orthonormality to switch to a simpler set of functions. Secondly, we allow “self-aliasing” of the cubature rule to relax bi-normality and correct for this normalization afterward. We shall see in the next section that this alternative approach has a connection with the method in Potts & Volkmer [26].

Remark 13** (Orthonormal and bi-orthonormal families).**

For function reconstruction we demand that the inner product of all basis functions in our support set are exactly represented by replacing the integral by a cubature rule. An alternative is to be able to exactly represent the inner product of all basis functions with another set of orthogonal functions which have the bi-orthonormal property. In general, if $\{{\mathrm{u}}_{\boldsymbol{k}}\}$ is an orthonormal basis with $\langle{\mathrm{u}}_{{\boldsymbol{k}}},{\mathrm{u}}_{{\boldsymbol{k}}^{\prime}}\rangle=\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$, and $\{{\mathrm{v}}_{\boldsymbol{k}}\}$ is an orthogonal set with the bi-orthonormal property $\langle{\mathrm{u}}_{\boldsymbol{k}},{\mathrm{v}}_{{\boldsymbol{k}}^{\prime}}\rangle=\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$, while $\langle{\mathrm{v}}_{{\boldsymbol{k}}},{\mathrm{v}}_{{\boldsymbol{k}}^{\prime}}\rangle=d_{\boldsymbol{k}}\,\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$ with $d_{\boldsymbol{k}}$ not necessarily equal to $1$, then the coefficients of a function $f\in\operatorname{span}\{{\mathrm{u}}_{\boldsymbol{k}}\}$ can be calculated by the inner product against either $\{{\mathrm{u}}_{\boldsymbol{k}}\}$ or $\{{\mathrm{v}}_{\boldsymbol{k}}\}$ since $\langle f,{\mathrm{u}}_{\boldsymbol{k}}\rangle=\langle f,{\mathrm{v}}_{\boldsymbol{k}}\rangle$. A cubature rule which can exactly calculate the inner products $\langle{\mathrm{u}}_{\boldsymbol{k}},{\mathrm{v}}_{{\boldsymbol{k}}^{\prime}}\rangle=\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$ for all ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}$ in our support set then also has the reconstruction property. We will make use of such a bi-orthonormal property below.

Lemma 14**.**

For any $f\in{\mathcal{F}}^{\rm cos}$ and ${\boldsymbol{k}}\in\mathbb{N}_{0}^{d}$, we can write the cosine coefficients in multiple ways

[TABLE]

Proof.

Using the Lebesgue preserving property of the tent transform and (14), we can write

[TABLE]

For each integral, we apply the change of variables $x_{j}^{\prime}=x_{j}$ if $\sigma_{j}=1$ and $x_{j}^{\prime}=1-x_{j}$ if $\sigma_{j}=-1$, and use the properties $\exp(2\pi{\mathrm{i}}\,\boldsymbol{\sigma}({\boldsymbol{k}})\cdot{\boldsymbol{x}})=\exp(2\pi{\mathrm{i}}\,{\boldsymbol{k}}\cdot\boldsymbol{\sigma}({\boldsymbol{x}}))=\exp(2\pi{\mathrm{i}}\,{\boldsymbol{k}}\cdot{\boldsymbol{x}}^{\prime})$ and $\varphi_{\rm tent}({\boldsymbol{x}})=\varphi_{\rm tent}({\boldsymbol{x}}^{\prime})$ to deduce that

[TABLE]

Thus all integrals are equal regardless of the sign changes on ${\boldsymbol{k}}$. Furthermore, since $f$ is a real-valued function, all its cosine coefficients will be real. Hence we may replace the exponential function $e_{\boldsymbol{\sigma}({\boldsymbol{k}})}$ by its real part \cos(2\pi\boldsymbol{\sigma}({\boldsymbol{k}})\cdot\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}). $\Box$

By considering the special case of $f=\phi_{{\boldsymbol{k}}^{\prime}}$ in Lemma 14 we obtain

[TABLE]

i.e., the functions $\{{\mathrm{u}}_{\boldsymbol{k}}=\phi_{\boldsymbol{k}}\circ\varphi_{\rm tent}\}$ and \{{\mathrm{v}}_{\boldsymbol{k}}=\sqrt{2}^{|{\boldsymbol{k}}|_{0}}\cos(2\pi{\boldsymbol{k}}\cdot\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}})\} are bi-orthonormal in $L^{2}$. Instead of demanding a cubature rule with exactness for the first inner product in (17) (as we did in (13)), below we seek a cubature rule with exactness for the second inner product in (17) (see (21) below), thus preserving bi-orthonormality.

Lemma 15** (Function reconstruction – plan B).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. A lattice rule $Q_{n}^{*}$ with $n$ points and generating vector ${\boldsymbol{z}}$ reconstructs exactly the cosine coefficients of all cosine space functions $f\in{\mathcal{F}}^{\rm cos}_{\Lambda}$ solely supported on $\Lambda$, by

[TABLE]

if and only if

[TABLE]

which is equivalent to

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Proof.

Substituting the cosine series of $f$ into (18), it follows that we have exact reconstruction of the cosine coefficients, i.e.,

[TABLE]

if and only if

[TABLE]

It remains to prove that (21) holds if and only if (19) holds.

Using (14) and \cos(2\pi{\boldsymbol{k}}\cdot\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}})=(e_{\boldsymbol{k}}+e_{-{\boldsymbol{k}}})/2, we can write

[TABLE]

where it is valid to replace one $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})$ by $-\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})$ since we sum over all unique sign changes.

Consider first the case ${\boldsymbol{k}}\neq{\boldsymbol{k}}^{\prime}$. Then $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})\neq{\boldsymbol{k}}$. By the character property (9) we know that $Q_{n}^{*}(e_{\pm(\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})-{\boldsymbol{k}})})$ can only take the values of $1$ or [math]. Thus (3.4) is equal to [math] if and only if all terms $Q_{n}^{*}(e_{\pm(\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})-{\boldsymbol{k}})})$ are [math], which holds following the character property if and only if $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})\cdot{\boldsymbol{z}}\not\equiv_{n}{\boldsymbol{k}}\cdot{\boldsymbol{z}}$.

Consider now the case ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime}$. Then we can rewrite (3.4) as

[TABLE]

Using the character property as before, we conclude that (23) is equal to $1$ if and only if all terms $Q_{n}^{*}(e_{\pm(\boldsymbol{\sigma}({\boldsymbol{k}})-{\boldsymbol{k}})})$ are [math] whenever $\boldsymbol{\sigma}({\boldsymbol{k}})\neq{\boldsymbol{k}}$, and in turn this means that $\boldsymbol{\sigma}({\boldsymbol{k}})\cdot{\boldsymbol{z}}\not\equiv_{n}{\boldsymbol{k}}\cdot{\boldsymbol{z}}$ except for when $\boldsymbol{\sigma}({\boldsymbol{k}})={\boldsymbol{k}}$. Combining all conditions, we conclude that (21) holds if and only if (19) holds.

Finally (19) is clearly equivalent to (20). The condition on $n$ then follows from Theorem 23, noting that the set $\Lambda\oplus{\mathcal{M}}(\Lambda)$ includes $\boldsymbol{0}$ but is not centrally symmetric. $\Box$

In the next lemma we propose another modification which allows “self-aliasing” in the lattice rule with respect to sign changes (see $\boldsymbol{\sigma}({\boldsymbol{k}})\cdot{\boldsymbol{z}}\equiv_{n}{\boldsymbol{k}}\cdot{\boldsymbol{z}}$ in (24) below). Consequently, the right-hand side of (21) for the case ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime}$ can be an integer $c_{\boldsymbol{k}}$, not necessarily $1$ (see (26) below). In other words, the cubature rule no longer preserves bi-normality, with normalization to be corrected by this factor $c_{\boldsymbol{k}}$.

Lemma 16** (Function reconstruction – plan C).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. A lattice rule $Q_{n}^{*}$ with $n$ points and generating vector ${\boldsymbol{z}}$ reconstructs exactly the cosine coefficients of all cosine space functions $f\in{\mathcal{F}}^{\rm cos}_{\Lambda}$ solely supported on $\Lambda$, by

[TABLE]

if and only if

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Proof.

Following the argument in the proof of Lemma 15, we now have exact reconstruction of the cosine coefficients if and only if (instead of (21))

[TABLE]

The case ${\boldsymbol{k}}\neq{\boldsymbol{k}}^{\prime}$ is the same as in Lemma 15. It suffices to reconsider the case ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime}$. Instead of separating out the term $\boldsymbol{\sigma}({\boldsymbol{k}})={\boldsymbol{k}}$ as in (23), we apply the character property (9) for $Q_{n}^{*}(e_{\pm(\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})-{\boldsymbol{k}})})$ in (3.4) with ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime}$ to arrive at

[TABLE]

which is equal to $c_{\boldsymbol{k}}$ as required. The CBC result is proved in Theorem 26 later. $\Box$

3.5 Fast calculation of cosine coefficients and function values

Here we can also make use of a one-dimensional fast Fourier transform to map cosine coefficients to function values on the tent-transformed lattice points and vice versa.

Lemma 17**.**

Let ${\boldsymbol{z}}$ be a generating vector for an $n$-point rank-$1$ lattice satisfying the reconstruction property on an arbitrary index set $\Lambda\subset\mathbb{N}_{0}^{d}$ according to Lemma 11 (plan A), Lemma 15 (plan B) or Lemma 16 (plan C). For a function $f\in{\mathcal{F}}^{\rm cos}_{\Lambda}$ solely supported on $\Lambda$ we can compute

[TABLE]

where ${\boldsymbol{f}}\in\mathbb{R}^{n}$ is a vector containing function values with the symmetry $f_{i}=f_{n-i}$ and ${\boldsymbol{F}}\in\mathbb{R}^{n}$ is a vector containing cosine coefficients with the symmetry $F_{\kappa}=F_{n-\kappa}$. Here $\mathrm{FFT}$ and $\mathrm{IFFT}$ are the one-dimensional fast Fourier transform and its inverse, respectively, with a normalization $1/n$ for $\mathrm{FFT}$ and $1$ for $\mathrm{IFFT}$; both mappings have cost ${\mathcal{O}}(n\log(n))$. For plan A and B we set $c_{\boldsymbol{k}}=1$ for all ${\boldsymbol{k}}$.

Alternatively, for even $n=2m$ a length $m+1$ DCT-I can be used, while for odd $n=2m-1$ a length $m$ DCT-V can be used. In this case the memory requirement and computational effort is halved (w.r.t. a real-to-real FFT).

Proof.

We show the result by using the calculation of Lemma 16 (plan C). Plan A and plan B are essentially the same with $c_{\boldsymbol{k}}=1$. For all plans we have the option to use the inner product with respect to $\cos(2\pi\,i\,{\boldsymbol{k}}\cdot{\boldsymbol{z}}/n)$ as given in Lemma 14 since the reconstruction property is guaranteed by the conditions of plan C which is a subset of the conditions of plan B which in turn is a subset of the conditions of plan A.

Let $f_{i}:=f(\varphi_{\rm tent}((i{\boldsymbol{z}}\bmod{n})/n))$ for $i=0,\ldots,n-1$. In all three plans we have for ${\boldsymbol{k}}\in\Lambda$,

[TABLE]

which follows because of the symmetry of $f_{i}=f_{n-i}$ due to the tent transform and the odd and even properties of the sine and cosine functions respectively. The last expression is a scaled discrete Fourier transform of length $n$ in terms of $i$ and $\kappa={\boldsymbol{k}}\cdot{\boldsymbol{z}}\bmod{n}$. This shows the calculation of coefficients from function values by one-dimensional FFT.

Next we consider the evaluation of function values from coefficients. For each $i=0,\ldots,n-1$, we have from (14) that

[TABLE]

where it is unfortunately not possible to avoid considering all sign changes of ${\boldsymbol{k}}$. For plan A and plan B there is only one $\kappa$ associated with each $\sigma({\boldsymbol{k}})$. For plan C it might occur that different sign changes on the same ${\boldsymbol{k}}$ map to the same value of $\kappa$, hence the summation in the algorithm to prepare the coefficient vector ${\boldsymbol{F}}$.

Now we explain how to make use of DCT using symmetry. We have $f_{i}=f_{n-i}$ due to the symmetry of $\varphi_{\rm tent}$. In the formula below, we will write $f_{n/2}$ which is to be interpreted in the way we just defined for even $n$, and to be considered equal to zero for odd $n$. Then, by making use of the symmetry, we have for general $n$ (odd or even)

[TABLE]

where $f_{n/2}\,\cos(\pi\,\kappa)$ is only present for even $n$. Now for $n=2m$ we have $\lfloor(n-1)/2\rfloor=n/2-1=m-1$ and we find

[TABLE]

which is the formula for the one-dimensional DCT-I of length $m+1$ to turn the sequence $f_{0},f_{1},\ldots,f_{m}$ into the sequence $F_{0},F_{1},\ldots,F_{m}$. For odd $n=2m-1$ we have $\lfloor(n-1)/2\rfloor=m-1$ and we find

[TABLE]

which is the formula for the one-dimensional DCT-V of length $m$ to turn the sequence $f_{0},f_{1}$, …, $f_{m-1}$ into the sequence $F_{0},F_{1},\ldots,F_{m-1}$. For DCT-I and DCT-V see Martucci [20] (definitions (A.1) and (A.5) therein).

Similarly, we have $F_{\kappa}=F_{n-\kappa}$ so we can write

[TABLE]

where $F_{n/2}:=0$ if $n$ is odd. Therefore DCT-I and DCT-V work in an analogous way. $\Box$

We note that all coefficients calculated by the FFT are real because of the symmetry in the input, but in a typical implementation one would take the real part in case of a complex FFT routine to remove possible numerical noise. Alternatively one can make use of a special real to real FFT implementation or use a specific implementation for the corresponding DCT.

Similarly to Remark 7, one can also extend the set $\Lambda$ such that $\sigma({\boldsymbol{k}})\cdot{\boldsymbol{z}}\bmod{n}$ covers as many values as possible in $\mathbb{Z}_{n}$ for functions which are not solely supported on $\Lambda$.

4 Nonperiodic setting based on Chebyshev polynomials

4.1 Chebyshev series

In the univariate case, the Chebyshev polynomials of the first kind for $|x|\leq 1$ can be written

[TABLE]

We have orthogonality with respect to the Chebyshev weight $(\sqrt{1-x^{2}})^{-1}$:

[TABLE]

To obtain an orthonormal basis on $[-1,1]$, we first normalize the measure to 1 by adjusting the Chebyshev weight to $(\pi\sqrt{1-x^{2}})^{-1}$. Then we define

[TABLE]

and for the multivariate case we define the tensor product basis functions

[TABLE]

where, as for cosine basis, $|{\boldsymbol{k}}|_{0}$ denotes the count of the nonzero entries in the vector ${\boldsymbol{k}}$.

Let ${\mathcal{F}}^{\rm Cheb}$ denote the space of real valued functions defined on $[-1,1]^{d}$ with absolutely converging Chebyshev series:

[TABLE]

where $\widehat{f}_{\boldsymbol{k}}$ are the Chebyshev coefficients of $f$. We equip ${\mathcal{F}}^{\rm Cheb}$ with the weighted $L^{2}$ inner product

[TABLE]

The orthonormal Chebyshev basis functions satisfy $\langle\eta_{\boldsymbol{k}},\eta_{{\boldsymbol{k}}^{\prime}}\rangle_{\mu}=\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$. The Chebyshev coefficients are

[TABLE]

Given an arbitrary finite index set $\Lambda\subset{\mathbb{N}}_{0}^{d}$, we consider the subspace ${\mathcal{F}}^{\rm Cheb}_{\Lambda}$ of all functions whose Chebyshev series is supported solely on $\Lambda$, i.e.,

[TABLE]

4.2 Isomorphism with the cosine space via cosine transform

The Chebyshev basis functions and the cosine basis functions are related by the mapping

[TABLE]

where the cosine function and its inverse are applied componentwise. This provides an isomorphism between the Chebyshev setting and the cosine space, with

[TABLE]

To get a sense of how this transformation works, it is informative to consider the univariate case:

[TABLE]

For the multivariate case we have the integral operator

[TABLE]

so that $I_{\mu}(\eta_{\boldsymbol{k}})=I(\phi_{\boldsymbol{k}})$, $I_{\mu}(\eta_{\boldsymbol{k}}\,\eta_{{\boldsymbol{k}}^{\prime}})=I(\phi_{\boldsymbol{k}}\,\phi_{{\boldsymbol{k}}^{\prime}})$, and

[TABLE]

i.e., the Chebyshev coefficients of $f$ are precisely the cosine coefficients of $f_{\rm cos}$.

4.3 Chebyshev coefficients by transformed rank-1 lattice rules

The isomorphism between the spaces means that we can bring all results from the cosine space over to the Chebyshev space. Noting the useful property

[TABLE]

we then arrive at the results for a tent-transformed and then cosine-transformed lattice rule, which is given explicitly by

[TABLE]

Since $\cos(2\pi\,i{\boldsymbol{z}}/n)=\cos(2\pi\,(n-i){\boldsymbol{z}}/n)$, the cubature points double up and we can write

[TABLE]

Thus the cubature rule can be computed with $\lfloor n/2+1\rfloor$ function evaluations, where the cubature weight for $i=0$ and $i=n/2$ (if $n$ is even) are $1/n$ and the others are $2/n$. In general there can be further duplication of points. However, if $\gcd(z_{j},n)=1$ for at least one $j=1,\ldots,d$, then all the $\lfloor n/2+1\rfloor$ points are distinct, see Lemma 12.

For even $n=2m$ this point set has previously been called a “Chebyshev lattice”, see, e.g., [6, 26], defined by \bigl{\{}\cos(\pi i{\boldsymbol{z}}/m):i=0,\ldots,m\bigr{\}}. However, we prefer the interpretation as a tent-transformed and then cosine-transformed lattice, since then we can also use odd $n$, and the cubature weights are automatically correct according to the multiplicity of the points.

We now state the analogous results to Lemmas 10, 11, 15, 16, 17 from the cosine space.

Lemma 18** (Integral exactness).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. A tent-transformed and then cosine transformed lattice rule with $n$ points and generating vector ${\boldsymbol{z}}$ integrates exactly (against the Chebyshev density) all Chebyshev space functions $f\in{\mathcal{F}}^{\rm Cheb}_{\Lambda}$ solely supported on $\Lambda$ if and only if

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Lemma 19** (Function reconstruction – plan A).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. A lattice rule $Q_{n}^{*}$ with $n$ points and generating vector ${\boldsymbol{z}}$ reconstructs exactly the Chebyshev coefficients of all Chebyshev space functions $f\in{\mathcal{F}}^{\rm Cheb}_{\Lambda}$ solely supported on $\Lambda$, by

[TABLE]

if and only if

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Lemma 20** (Function reconstruction – plan B).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. A lattice rule $Q_{n}^{*}$ with $n$ points and generating vector ${\boldsymbol{z}}$ reconstructs exactly the Chebyshev coefficients of all Chebyshev space functions $f\in{\mathcal{F}}^{\rm Cheb}_{\Lambda}$ solely supported on $\Lambda$, by

[TABLE]

if and only if

[TABLE]

which is equivalent to

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Lemma 21** (Function reconstruction – plan C).**

Let $\Lambda\subset\mathbb{N}_{0}^{d}$ be an arbitrary index set. A lattice rule $Q_{n}^{*}$ with $n$ points and generating vector ${\boldsymbol{z}}$ reconstructs exactly the Chebyshev coefficients of all Chebyshev space functions $f\in{\mathcal{F}}^{\rm Cheb}_{\Lambda}$ solely supported on $\Lambda$, by

[TABLE]

if and only if

[TABLE]

Such a generating vector ${\boldsymbol{z}}$ can be constructed component-by-component if $n$ is a prime satisfying

[TABLE]

Lemma 22**.**

We can use FFTs or DCTs to map Chebyshev coefficients to function values on tent-transformed and then cosine-transformed lattice points, and the other way round, for an $n$-point rank-$1$ lattice rule with generating vector ${\boldsymbol{z}}$ satisfying the non-aliasing conditions of Lemma 19 (plan A), Lemma 20 (plan B) or Lemma 21 (plan C) on an arbitrary index set $\Lambda\subset\mathbb{N}_{0}^{d}$ by replacing \varphi_{\rm tent}(\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}) by \cos(\pi\,\varphi_{\rm tent}(\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}})) in the statement of Lemma 17.

Lemma 21 with even $n=2m$ in combination with the DCT-I in Lemma 22 is essentially the approach in Potts & Volkmer [26].

5 The CBC construction

5.1 Induction proof for the component-by-component construction

Let ${\mathbb{Z}}_{n}^{*}:=\{1,\ldots,n-1\}$ for $n$ prime. The necessary and sufficient conditions on the lattice rule generating vector ${\boldsymbol{z}}\in({\mathbb{Z}}_{n}^{*})^{d}$ for integral exactness and function reconstruction in most cases boil down to the same generic form of verifying for a given index set ${\mathcal{A}}\subset{\mathbb{Z}}^{d}$ that

[TABLE]

The following theorem justifies a generic component-by-component (CBC) algorithm to find a ${\boldsymbol{z}}$ satisfying this condition. The inductive argument needs to work with projections of the index set ${\mathcal{A}}$ down to the lower coordinates. We consider two definitions for the projections since each has its advantages:

[TABLE]

Both definitions yield ${\mathcal{A}}_{d}={\mathcal{A}}$. The definition (30) includes only the indices whose higher components are zero; we shall refer to this as the ‘zero’ projection. The definition (31) includes all indices obtained by truncating the original indices; we shall refer to this as the ‘full’ projection. If the index set ${\mathcal{A}}$ is downward closed then they are the same; otherwise (30) is a subset of (31). The full projection (31) was used in [4, 9, 10]; the zero projection (30) is new in this paper.

The condition $n>\max({\mathcal{A}})$ in the theorem guarantees that the components of ${\boldsymbol{h}}\in{\mathcal{A}}$ all satisfy $|h_{j}|<n$, and thus $h_{j}\equiv_{n}0$ if and only if $h_{j}=0$. This condition can be replaced by the direct assumption that there is no ${\boldsymbol{h}}\in{\mathcal{A}}$ with a nonzero component $h_{j}$ that is a multiple of $n$.

Theorem 23**.**

Let ${\mathcal{A}}\subset{\mathbb{Z}}^{d}$ be an arbitrary index set, and let $n$ be a prime number satisfying

[TABLE]

with $\kappa=2$ if ${\mathcal{A}}$ is centrally symmetric and $\kappa=1$ otherwise. Define the projections ${\mathcal{A}}_{s}$ by (30) or (31). Then a generating vector ${\boldsymbol{z}}^{*}=(z_{1},\ldots,z_{d})\in(\mathbb{Z}^{*}_{n})^{d}$ can be constructed component-by-component such that for all $s=1,\ldots,d$ and ${\boldsymbol{z}}=(z_{1},\ldots,z_{s})$ we have

[TABLE]

Proof.

The result for centrally symmetric and downward closed index sets (e.g. hyperbolic cross) or more general index sets with the full projection (31) has been proved in [4, 9, 10]. So we focus on proving the general result with the zero projection (30).

The proof is by induction on $s$. We will attempt to derive the condition (32) rather than assuming it from the beginning.

For $s=1$, the condition $h_{1}z_{1}\not\equiv_{n}0$ holds for all $z_{1}\in\mathbb{Z}^{*}_{n}$ if $h_{1}\not\equiv_{n}0$, and fails for all $z_{1}$ if $h_{1}\equiv_{n}0$. To avoid the latter scenario we assume that $n>|h_{1}|$ always holds.

Suppose we already obtained the generating vector ${\boldsymbol{z}}\in(\mathbb{Z}^{*}_{n})^{s-1}$ satisfying (33) for ${\mathcal{A}}_{s-1}$. For each $({\boldsymbol{h}},h_{s})\in{\mathcal{A}}_{s}\!\setminus\!\{\boldsymbol{0}\}$, we will eliminate any ‘bad’ $z_{s}\in{\mathbb{Z}}_{n}^{*}$ that satisfies

[TABLE]

We stress that $({\boldsymbol{h}},h_{s})\in{\mathcal{A}}_{s}$ does not imply ${\boldsymbol{h}}\in{\mathcal{A}}_{s-1}$ under the zero projection (30). Depending on the value of $h_{s}$ we have the following scenarios:

If $h_{s}=0$ then ${\boldsymbol{h}}\in{\mathcal{A}}_{s-1}\!\setminus\!\{\boldsymbol{0}\}$; in turn the induction hypotheses (33) for ${\mathcal{A}}_{s-1}$ guarantees that ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\not\equiv_{n}0$ and so (34) has no solution for $z_{s}$. There are $\#({\mathcal{A}}_{s-1}\!\setminus\!\{\boldsymbol{0}\})$ such cases. 2. 2.

If $h_{s}\neq 0$ but $h_{s}\equiv_{n}0$, then (34) has no solution for $z_{s}$ if ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\not\equiv_{n}0$, or it holds for all $z_{s}$ if ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\equiv_{n}0$. To avoid the latter scenario we assume that $n>|h_{s}|$ always holds. 3. 3.

If $h_{s}\not\equiv_{n}0$ then, since $n$ is prime, (34) has a unique solution for $z_{s}\in\mathbb{Z}^{*}_{n}$ if ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\not\equiv_{n}0$, or has no solution for $z_{s}\in\mathbb{Z}^{*}_{n}$ if ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\equiv_{n}0$. The latter scenario includes ${\boldsymbol{h}}=\boldsymbol{0}$ so there are at least $\#{{\mathcal{A}}_{s}^{\dagger}}$ such cases, where ${\mathcal{A}}_{s}^{\dagger}:=\{h_{s}\in{\mathbb{Z}}:h_{s}\not\equiv_{n}0\mbox{ and }(\boldsymbol{0},h_{s})\in{\mathcal{A}}_{s}\}$.

Thus, provided that $n>|h_{s}|$, there is at most one bad $z_{s}$ to be eliminated for each $({\boldsymbol{h}},h_{s})\in{\mathcal{A}}_{s}\!\setminus\!\{\boldsymbol{0}\}$, and the total number of bad $z_{s}$ we eliminate is at most $\#({\mathcal{A}}_{s}\!\setminus\!\{\boldsymbol{0}\})-\#({\mathcal{A}}_{s-1}\!\setminus\!\{\boldsymbol{0}\})-\#{{\mathcal{A}}_{s}^{\dagger}}$.

Hence, provided additionally that $\#{\mathbb{Z}^{*}_{n}}=n-1>\#{{\mathcal{A}}_{s}}-\#{{\mathcal{A}}_{s-1}}-\#{{\mathcal{A}}_{s}^{\dagger}}$, there is always a ‘good’ $z_{s}$ remaining such that $({\boldsymbol{z}},z_{s})$ will satisfy (33) for ${\mathcal{A}}_{s}$. Moreover if ${\mathcal{A}}$ is centrally symmetric, then all ${\mathcal{A}}_{s}$ are centrally symmetric, and both $({\boldsymbol{h}}^{\prime},h_{s})$ and $(-{\boldsymbol{h}},-h_{s})$ will eliminate the same $z_{s}$ if a solution for (34) exists. It then suffices to demand that $n-1>(\#{{\mathcal{A}}_{s}}-\#{{\mathcal{A}}_{s-1}}-\#{{\mathcal{A}}_{s}^{\dagger}})/2$. By induction, to ensure that a good ${\boldsymbol{z}}^{*}\in({\mathbb{Z}}_{n}^{*})^{d}$ exists, it suffices to assume that

[TABLE]

with $\kappa=2$ if ${\mathcal{A}}$ is centrally symmetric and $\kappa=1$ otherwise. This leads to the simplified condition on $n$ in the theorem.

Now for completeness we discuss briefly the case for the full projection (31). The proof is almost identical to the case for the zero projection but is slightly simpler. The subtle difference is that for each $({\boldsymbol{h}},h_{s})\in{\mathcal{A}}_{s}$ we now have ${\boldsymbol{h}}\in{\mathcal{A}}_{s-1}$ regardless of the value of $h_{s}$, and the induction hypothesis (33) for ${\mathcal{A}}_{s-1}$ guarantees that ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\not\equiv_{n}0$ if $h_{s}\neq 0$, thus simplifying the second and third scenarios above. $\Box$

5.2 Algorithmic aspects of the CBC construction

Remark 24**.**

For the projections ${\mathcal{A}}_{s}$ defined by either (30) or (31), Theorem 23 and its proof justify two different approaches to carry out the component-by-component construction:

• •

Brute force approach: At step $s$, we search through $z_{s}\in\mathbb{Z}^{*}_{n}$ until we find one that satisfies (33) for all ${\boldsymbol{h}}\in{\mathcal{A}}_{s}\!\setminus\!\{\boldsymbol{0}\}$. The cost is ${\mathcal{O}}(n_{\rm fail}\,\#{{\mathcal{A}}_{s}})$, where $n_{\rm fail}$ is the number of different $z_{s}$ that was checked. So the cost is at worst ${\mathcal{O}}(n\,\#{{\mathcal{A}}_{s}})$, leading to a total cost of ${\mathcal{O}}(d\,n\,\#{{\mathcal{A}}})$.

• •

Elimination approach: At step $s$, we loop through every ${\boldsymbol{h}}\in{\mathcal{A}}_{s}\!\setminus\!\{\boldsymbol{0}\}$ and eliminate the corresponding $z_{s}\in\mathbb{Z}^{*}_{n}$ that fails (33), if any. Then we take any remaining $z_{s}$. The cost is only ${\mathcal{O}}(\#{{\mathcal{A}}_{s}})$, leading to a total cost of ${\mathcal{O}}(d\,\#{{\mathcal{A}}})$.

In both approaches we have left out the ${\mathcal{O}}(s)$ factor in step $s$ that arises from the evaluation of dot products; this is valid if we store and update the values of ${\boldsymbol{h}}\cdot{\boldsymbol{z}}$ for all ${\boldsymbol{h}}$ in each step. In the elimination approach we mark the bad choices of $z_{s}$ in an array of length $n-1$ with pointers linking the previous and next good choices of $z_{s}$, so that it is ${\mathcal{O}}(1)$ cost to obtain a good $z_{s}$ at the end.

The two approaches may be used for different steps in the algorithm if it is advantageous to mix them. Both approaches are guaranteed to succeed provided $n$ is sufficiently large, see (32). We can run the algorithm with smaller values of $n$ (or even composite values of $n$ in the brute force approach), and be prepared to increase $n$ when the algorithm fails. Once a ${\boldsymbol{z}}$ is found, we can systematically test and reduce the value of $n$ by verifying whether (33) still holds for ${\mathcal{A}}$.

In general the zero projections (30) are subsets of the full projections (31), and consequently the condition (33) is weaker and faster to check for (30) than for (31). There are also algorithmic advantages in the data structure for iterating the sets based on (30), namely, that the indices can be ordered according to the number of zeros at the end.

We now apply Theorem 23 to the situation where the input set ${\mathcal{A}}$ is a difference set, i.e., ${\mathcal{A}}=\Lambda\ominus\Lambda$. Then the condition (33) is now explicitly given by

[TABLE]

where $(\Lambda\ominus\Lambda)_{s}$ denote the projection of the difference set $\Lambda\ominus\Lambda$ defined according to (30) or (31). Since $\#(\Lambda\ominus\Lambda)\leq(\#{\Lambda})^{2}$, Remark 24 indicates that the cost of CBC construction for ${\mathcal{A}}=\Lambda\ominus\Lambda$ is ${\mathcal{O}}(d\,n\,(\#{\Lambda})^{2})$ or ${\mathcal{O}}(d\,(\#{\Lambda})^{2})$ for the two approaches, respectively.

Similarly, for all sign changes on an index set $\Lambda$ we have $\#{{\mathcal{M}}(\Lambda)}\leq 2^{d}\,\#{\Lambda}$ and therefore the cost of CBC construction for ${\mathcal{A}}={\mathcal{M}}(\Lambda)\oplus{\mathcal{M}}(\Lambda)$ is ${\mathcal{O}}(d\,n\,2^{2d}\,(\#{\Lambda})^{2})$ or ${\mathcal{O}}(d\,2^{2d}\,(\#{\Lambda})^{2})$; the cost for ${\mathcal{A}}=\Lambda\oplus{\mathcal{M}}(\Lambda)$ is ${\mathcal{O}}(d\,n\,2^{d}\,(\#{\Lambda})^{2})$ or ${\mathcal{O}}(d\,2^{d}\,(\#{\Lambda})^{2})$.

Hence, we can apply Theorem 23 and Remark 24 to the Fourier space in Lemmas 4 and 5, noting that the difference set $\Lambda\ominus\Lambda$ is always centrally symmetric and contains the zero vector. Analogously, we can apply Theorem 23 and Remark 24 to the cosine space in Lemmas 10 and 11, as well as Lemma 15 – plan B, and correspondingly, to the Chebyshev space in Lemmas 18, 19, and 20.

However, plan C for the cosine space and Chebyshev space, see Lemmas 16 and 21, respectively, cannot be formulated in the same generic form (29). So we will need to develop a separate justification for it. We will return to this later in Subsection 5.5.

5.3 Smart lookup for the brute force approach with full projection

When ${\mathcal{A}}_{s}$ are full projections (31), we have the important property that the projection of the difference set equals the difference set of the projections, i.e.,

[TABLE]

Thus the condition (35) becomes

[TABLE]

which is equivalent to

[TABLE]

In other words, every dot product needs to have a unique value.

The following code snippet shows that it is possible to verify condition (37) for a given ${\boldsymbol{z}}$ with cost ${\mathcal{O}}(\#{\Lambda_{s}})$ rather than ${\mathcal{O}}((\#{\Lambda_{s}})^{2})$, by marking a bit string of length $n$ for the values of dot product modulo $n$ that have occurred.

⬇

// Fourier space

// INPUT: ${\boldsymbol{z}}$, $n$ and $\Lambda_{s}$

// VERIFY: ${\boldsymbol{h}}\cdot{\boldsymbol{z}}\not\equiv_{n}{\boldsymbol{h}}^{\prime}\cdot{\boldsymbol{z}}$ for all ${\boldsymbol{h}},{\boldsymbol{h}}^{\prime}\in\Lambda_{s}$, ${\boldsymbol{h}}\neq{\boldsymbol{h}}^{\prime}$

// COST: ${\mathcal{O}}(\#{\Lambda_{s}})$

$S=0$

for ${\boldsymbol{h}}\in\Lambda_{s}$:

 $\alpha={\boldsymbol{h}}\cdot{\boldsymbol{z}}\bmod{n}$

 **if** $S[\alpha]=1$: **return** FALSE

 $S[\alpha]=1$

return TRUE

Consequently, we can reduce the cost of the brute force CBC construction for ${\mathcal{A}}=\Lambda\ominus\Lambda$ from ${\mathcal{O}}(d\,n\,(\#{\Lambda})^{2})$ to ${\mathcal{O}}(d\,n\,\#{\Lambda})$. We shall refer to this as the “smart lookup” trick.

We stress once again that (36) only holds when we have the full projection (31). Under the zero projection (30) we would have in general $(\Lambda\ominus\Lambda)_{s}\supseteq\Lambda_{s}\ominus\Lambda_{s}$; in this case the alternative formulation (37) would miss out on some indices.

Similar reduction in cost can be achieved for ${\mathcal{A}}={\mathcal{M}}(\Lambda)\oplus{\mathcal{M}}(\Lambda)$ as we show in the code snippet below.

⬇

// Cosine space and Chebyshev space -- plan A

// INPUT: ${\boldsymbol{z}}$, $n$ and $\Lambda_{s}$

// VERIFY: $\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}^{\prime})\cdot{\boldsymbol{z}}\not\equiv_{n}\boldsymbol{\sigma}({\boldsymbol{k}})\cdot{\boldsymbol{z}}$ for all ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in\Lambda_{s}$, $\boldsymbol{\sigma}\in{\mathcal{S}}_{\boldsymbol{k}}$, $\boldsymbol{\sigma}^{\prime}\in{\mathcal{S}}_{{\boldsymbol{k}}^{\prime}}$, $\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}^{\prime})\neq\boldsymbol{\sigma}({\boldsymbol{k}})$

// COST: ${\mathcal{O}}(\#{{\mathcal{M}}(\Lambda_{s}}))$

$S=0$

for ${\boldsymbol{k}}\in\Lambda_{s}$:

 **for** ${\boldsymbol{h}}\in\{\boldsymbol{\sigma}({\boldsymbol{k}}):\boldsymbol{\sigma}\in{\mathcal{S}}_{\boldsymbol{k}}\}$:

     $\alpha={\boldsymbol{h}}\cdot{\boldsymbol{z}}\bmod{n}$

     **if** $S[\alpha]=1$: **return** FALSE

     $S[\alpha]=1$

return TRUE

We can save on half of the calculations, since if ${\boldsymbol{k}}\neq\boldsymbol{0}$ we can fix one of the signs for a non-zero element of ${\boldsymbol{k}}$ to get half of the sign changes and multiply by $-1$ to get the other half as shown below.

⬇

// Cosine space and Chebyshev space -- plan A -- halved

// INPUT: ${\boldsymbol{z}}$, $n$ and $\Lambda_{s}$

// VERIFY: $\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}^{\prime})\cdot{\boldsymbol{z}}\not\equiv_{n}\boldsymbol{\sigma}({\boldsymbol{k}})\cdot{\boldsymbol{z}}$ for all ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in\Lambda_{s}$, $\boldsymbol{\sigma}\in{\mathcal{S}}_{\boldsymbol{k}}$, $\boldsymbol{\sigma}^{\prime}\in{\mathcal{S}}_{{\boldsymbol{k}}^{\prime}}$, $\boldsymbol{\sigma}^{\prime}({\boldsymbol{k}}^{\prime})\neq\boldsymbol{\sigma}({\boldsymbol{k}})$

// COST: ${\mathcal{O}}(\#{{\mathcal{M}}(\Lambda_{s}})/2)$

$S=0$

if $\boldsymbol{0}\in\Lambda_{s}$: $S[0]=1$

for ${\boldsymbol{k}}\in\Lambda_{s}\setminus\{\boldsymbol{0}\}$:

 $i=\min\{j:k_{j}\neq 0\}$

 **for** ${\boldsymbol{h}}\in\{\boldsymbol{\sigma}({\boldsymbol{k}}):\boldsymbol{\sigma}\in{\mathcal{S}}_{\boldsymbol{k}},\sigma_{i}=+1\}$:

     $\alpha={\boldsymbol{h}}\cdot{\boldsymbol{z}}\bmod{n}$        *// note: $n-\alpha=-{\boldsymbol{h}}\cdot{\boldsymbol{z}}\bmod{n}$*

     **if** $S[\alpha]=1$: **return** FALSE

     $S[\alpha]=1$

     **if** $S[n-\alpha]=1$: **return** FALSE   *// note: it could happen that $n-\alpha=\alpha$*

     $S[n-\alpha]=1$

return TRUE

We can achieve a similar reduction in cost for ${\mathcal{A}}=\Lambda\oplus{\mathcal{M}}(\Lambda)$, but this is more complicated because we need to distinguish between the dot products coming from the original indices and the dot products coming from sign changes of the indices. We do this by keeping two bit strings of length $n$ as shown in the code snippet below: $S_{1}$ marks the original dot products, while $S_{2}$ marks the dot products from all sign changes thus including $S_{1}$. (We can also half the cost as above but we do not include that here.)

⬇

// Cosine space and Chebyshev space -- plan B

// INPUT: ${\boldsymbol{z}}$, $n$ and $\Lambda_{s}$

// VERIFY: $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})\cdot{\boldsymbol{z}}\not\equiv_{n}{\boldsymbol{k}}\cdot{\boldsymbol{z}}$ for all ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in\Lambda_{s}$, $\boldsymbol{\sigma}\in{\mathcal{S}}_{{\boldsymbol{k}}^{\prime}}$, $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})\neq{\boldsymbol{k}}$

// COST: ${\mathcal{O}}(\#{{\mathcal{M}}(\Lambda_{s}}))$

$S_{1}=0$

$S_{2}=0$

for ${\boldsymbol{k}}\in\Lambda_{s}$:

 $\alpha={\boldsymbol{k}}\cdot{\boldsymbol{z}}\bmod{n}$

 **if** $S_{2}[\alpha]=1$: **return** FALSE     *// note: the value of $S_{1}[\alpha]$ is also checked since $S_{1}\subseteq S_{2}$*

 $S_{2}[\alpha]=1$

 $S_{1}[\alpha]=1$

 **for** ${\boldsymbol{h}}\in\{\boldsymbol{\sigma}({\boldsymbol{k}}):\boldsymbol{\sigma}\in{\mathcal{S}}_{\boldsymbol{k}}\}$ with ${\boldsymbol{h}}\neq{\boldsymbol{k}}$:

     $\alpha^{\prime}={\boldsymbol{h}}\cdot{\boldsymbol{z}}\bmod{n}$

     **if** $S_{1}[\alpha^{\prime}]=1$: **return** FALSE

     $S_{2}[\alpha^{\prime}]=1$        *// note: it does not matter if $S_{2}[\alpha^{\prime}]$ is already set*

return TRUE

The previous algorithm can now be modified to allow self-aliasing and to keep track of the constant $c_{\boldsymbol{k}}$, see Lemma 16 for cosine space and Lemma 21 for Chebyshev space.

⬇

// Cosine space and Chebyshev space -- plan C

// INPUT: ${\boldsymbol{z}}$, $n$ and $\Lambda_{s}$

// VERIFY: $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})\cdot{\boldsymbol{z}}\not\equiv_{n}{\boldsymbol{k}}\cdot{\boldsymbol{z}}$ for all ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in\Lambda_{s}$, $\boldsymbol{\sigma}\in{\mathcal{S}}_{{\boldsymbol{k}}^{\prime}}$, ${\boldsymbol{k}}^{\prime}\neq{\boldsymbol{k}}$

// COST: ${\mathcal{O}}(\#{{\mathcal{M}}(\Lambda_{s}}))$

$S_{1}=0$

$S_{2}=0$

for ${\boldsymbol{k}}\in\Lambda_{s}$:

 $\alpha={\boldsymbol{k}}\cdot{\boldsymbol{z}}\bmod{n}$

 **if** $S_{2}[\alpha]=1$: **return** FALSE

 $S_{2}[\alpha]=1$

 $c_{\boldsymbol{k}}=1$

 **for** ${\boldsymbol{h}}\in\{\boldsymbol{\sigma}({\boldsymbol{k}}):\boldsymbol{\sigma}\in{\mathcal{S}}_{\boldsymbol{k}}\}$ with ${\boldsymbol{h}}\neq{\boldsymbol{k}}$:

     $\alpha^{\prime}={\boldsymbol{h}}\cdot{\boldsymbol{z}}\bmod{n}$

     **if** $\alpha^{\prime}=\alpha$: $c_{\boldsymbol{k}}$ += $1$

     **if** $S_{1}[\alpha^{\prime}]=1$: **return** FALSE

     $S_{2}[\alpha^{\prime}]=1$

 $S_{1}[\alpha]=1$

return TRUE with $\{c_{\boldsymbol{k}}:{\boldsymbol{k}}\in\Lambda_{s}\}$

The crucial difference between the last two code snippets is that in plan B the bit $S_{1}[\alpha]$ is marked before the dot products $\alpha^{\prime}$ from the sign changes are checked against $S_{1}$, thus not allowing $\alpha^{\prime}=\alpha$ (no aliasing), while in plan C the bit $S_{1}[\alpha]$ is marked only after all $\alpha^{\prime}$ have been checked against $S_{1}$, thus allowing $\alpha^{\prime}=\alpha$ and indeed counts the number of times this occurs in $c_{\boldsymbol{k}}$ (self-aliasing).

We summarize this subsection in the following remark.

Remark 25**.**

The cost for brute force CBC in step $s$ with full projection (31) and smart lookup is

[TABLE]

for the index sets ${\mathcal{A}}=\Lambda\ominus\Lambda$, ${\mathcal{A}}={\mathcal{M}}(\Lambda)\oplus{\mathcal{M}}(\Lambda)$ and ${\mathcal{A}}=\Lambda\oplus{\mathcal{M}}(\Lambda)$, respectively.

Similar strategies have been implemented in the sparseFFTr1l software library of Toni Volkmer [32].

5.4 Mixed CBC construction

Combining Remark 25 with Remark 24, we see that there is advantage in mixing the two different approaches. As long as $n_{\rm fail}$ remains small it is advantageous to follow the brute force approach with full projection (31) and smart lookup. We anticipate this to be the case for the initial dimensions.

Starting from $z_{1}=1$, at step $s$ we begin our brute force search with the value $z_{s}=z_{s-1}+1$. If this $z_{s}$ fails then we increment again by $1$ and do this repeatedly (if $n-1$ is reached then we continue from $1$) until a valid $z_{s}$ is found, while keeping a count on $n_{\rm fail}$. Then gradually as the dimension increases and as we run out of choices, we expect the value of $n_{\rm fail}$ to increase until at some point the balance tips over the other way and it becomes cheaper to follow the elimination approach. From then on we switch over to the elimination approach in the generic formulation (29) with the zero projection (30) so that the sets are smaller (except for the case of plan C which we discuss in the next subsection).

We summarize our results for the different spaces in Table 1.

5.5 A new CBC proof for plan C

Recall that the condition (25) is weaker than the condition (19), which is in turn equivalent to (20). Thus when we have the full projection (31), the condition on $n$ in Theorem 23 guarantees the existence of $z_{s}$ with the required property in step $s$. However, to prove that the CBC construction can find this vector, we need a new CBC proof.

Theorem 26**.**

Let $\Lambda\subset{\mathbb{N}}^{d}_{0}$ be an arbitrary index set, and let $n$ be a prime number satisfying

[TABLE]

Define $\Lambda_{s}$ to be the full projection of $\Lambda$ as in (31). Then a generating vector ${\boldsymbol{z}}^{*}=(z_{1},\ldots,z_{d})\in(\mathbb{Z}^{*}_{n})^{d}$ can be constructed component-by-component such that for all $s=1,\ldots,d$ and ${\boldsymbol{z}}=(z_{1},\ldots,z_{s})$ we have

[TABLE]

Proof.

The proof is by induction on $s$. For $s=1$ and $k_{1}\neq k_{1}^{\prime}$, the condition $\sigma_{1}(k_{1}^{\prime})z_{1}\not\equiv_{n}k_{1}z_{1}$ holds for all $z_{1}\in\mathbb{Z}^{*}_{n}$ if $\sigma_{1}(k_{1}^{\prime})-k_{1}\not\equiv_{n}0$, and fails for all $z_{1}$ if $\sigma_{1}(k_{1}^{\prime})-k_{1}\equiv_{n}0$. To avoid the latter scenario we assume that $n>2\max_{k_{1}\in\Lambda_{1}}|k_{1}|$.

Suppose we already obtained the generating vector ${\boldsymbol{z}}\in(\mathbb{Z}^{*}_{n})^{s-1}$ satisfying (39) for $\Lambda_{s-1}$. For each distinct pair $({\boldsymbol{k}},k_{s}),({\boldsymbol{k}}^{\prime},k_{s}^{\prime})\in\Lambda_{s}$ and each $(\boldsymbol{\sigma},\sigma_{s})\in{\mathcal{S}}_{({\boldsymbol{k}}^{\prime},k_{s}^{\prime})}$, we will eliminate any ‘bad’ $z_{s}\in{\mathbb{Z}}_{n}^{*}$ that satisfies

[TABLE]

From the definition (31) we have ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}\in\Lambda_{s-1}$. We have the following scenarios:

If ${\boldsymbol{k}}\neq{\boldsymbol{k}}^{\prime}$ then the induction hypotheses (39) for $\Lambda_{s-1}$ guarantees that $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})\cdot{\boldsymbol{z}}\not\equiv_{n}{\boldsymbol{k}}\cdot{\boldsymbol{z}}$ for all $\boldsymbol{\sigma}$. Thus (40) has no solution for $z_{s}$ if $\sigma_{s}(k_{s}^{\prime})-k_{s}\equiv_{n}0$, and (40) has a unique solution for $z_{s}$ if $\sigma_{s}(k_{s}^{\prime})-k_{s}\not\equiv_{n}0$. 2. 2.

If ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime}$ (thus $k_{s}\neq k_{s}^{\prime}$) and $\boldsymbol{\sigma}$ satisfies $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})\cdot{\boldsymbol{z}}\not\equiv_{n}{\boldsymbol{k}}\cdot{\boldsymbol{z}}$ then, as in the previous scenario, (40) has no solution for $z_{s}$ if $\sigma_{s}(k_{s}^{\prime})-k_{s}\equiv_{n}0$, and (40) has a unique solution for $z_{s}$ if $\sigma_{s}(k_{s}^{\prime})-k_{s}\not\equiv_{n}0$. 3. 3.

If ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime}$ (thus $k_{s}\neq k_{s}^{\prime}$) and $\boldsymbol{\sigma}$ satisfies $\boldsymbol{\sigma}({\boldsymbol{k}}^{\prime})\cdot{\boldsymbol{z}}\equiv_{n}{\boldsymbol{k}}\cdot{\boldsymbol{z}}$, then (40) has no solution for $z_{s}$ if $\sigma_{s}(k_{s}^{\prime})-k_{s}\not\equiv_{n}0$, and (40) holds for all $z_{s}$ if $\sigma_{s}(k_{s}^{\prime})-k_{s}\equiv_{n}0$. To avoid the latter scenario we assume that $n>|\sigma_{s}(k_{s}^{\prime})-k_{s}|$. Note that it is not possible to have $\sigma_{s}(k_{s}^{\prime})=k_{s}$ when $k_{s}\neq k_{s}^{\prime}$ since both $k_{s}$ and $k_{s}^{\prime}$ are nonnegative integers.

Thus, provided that $n>2\max_{({\boldsymbol{k}},k_{s})\in\Lambda_{s}}|k_{s}|$, there is at most one bad $z_{s}$ to be eliminated for each distinct pair $({\boldsymbol{k}},k_{s}),({\boldsymbol{k}}^{\prime},k_{s}^{\prime})\in\Lambda_{s}$ and each $(\boldsymbol{\sigma},\sigma_{s})\in{\mathcal{S}}_{({\boldsymbol{k}}^{\prime},k_{s}^{\prime})}$, so the total number of bad $z_{s}$ we eliminate is at most $\#{\Lambda_{s}}\,(\#{{\mathcal{M}}(\Lambda_{s})}-1)$.

Hence, provided additionally that $\#{\mathbb{Z}^{*}_{n}}=n-1>\#{\Lambda_{s}}\,(\#{{\mathcal{M}}(\Lambda_{s})}-1)$, there is always a ‘good’ $z_{s}$ remaining such that $({\boldsymbol{z}},z_{s})$ will satisfy (39) for $\Lambda_{s}$. By induction, to ensure that a good ${\boldsymbol{z}}^{*}\in({\mathbb{Z}}_{n}^{*})^{d}$ exists, it suffices to assume that $n$ satisfies (38). This completes the proof. $\Box$

6 Approximation

We now discuss function approximation for all three settings under a unified framework. A major difference of this section compared to the previous sections is that the function $f$ under consideration is no longer supported only on a finite index set. We cannot achieve exact function reconstruction and therefore an error analysis is needed.

6.1 Function approximation under a unified framework

We have an orthonormal basis $\{\alpha_{\boldsymbol{k}}\}$ for $L_{\mu}^{2}(\Omega)$, where $\mu(\Omega)=1$, and consider functions with absolutely converging series expansions

[TABLE]

where the sum is over ${\mathbb{Z}}^{d}$ (for Fourier space) or ${\mathbb{N}}_{0}^{d}$ (for cosine and Chebyshev spaces). Now consider a subset $\Lambda$ of the indices and represent the exact $L_{\mu}^{2}$ projection of $f$, i.e., the best $L_{\mu}^{2}$ approximation on $\Lambda$, by

[TABLE]

We cannot calculate these coefficients $\widehat{f}_{\boldsymbol{k}}$ exactly and will have to approximate them, leading to

[TABLE]

with $\square\in\{a,b,c\}$ denoting the approximation by plan A (including the Fourier case), plan B, or plan C.

We write

[TABLE]

Then we have $\langle{\rm u}_{\boldsymbol{k}},{\mathrm{u}}_{{\boldsymbol{k}}^{\prime}}\rangle=\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$, $\langle{\mathrm{u}}_{\boldsymbol{k}},{\mathrm{v}}_{{\boldsymbol{k}}^{\prime}}\rangle=\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$, and $\langle{\mathrm{v}}_{\boldsymbol{k}},{\mathrm{v}}_{{\boldsymbol{k}}^{\prime}}\rangle=d_{\boldsymbol{k}}\,\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$, with $d_{\boldsymbol{k}}=2^{|{\boldsymbol{k}}|_{0}-1}$ for ${\boldsymbol{k}}\neq\boldsymbol{0}$ in the case of cosine or Chebyshev space plan B or C and $d_{\boldsymbol{k}}=1$ otherwise.

We demand from our lattice rule that

[TABLE]

where $c_{\boldsymbol{k}}$ in the case of cosine or Chebyshev space plan C (see (24) or (28)) can be a positive integer up to the number of unique sign changes of ${\boldsymbol{k}}$, i.e., $1\leq c_{\boldsymbol{k}}\leq 2^{|{\boldsymbol{k}}|_{0}}$, and $c_{\boldsymbol{k}}=1$ otherwise. Note that we do not necessarily have $Q_{n}({\mathrm{u}}_{\boldsymbol{k}}\,\overline{{\mathrm{u}}_{{\boldsymbol{k}}^{\prime}}})=\delta_{{\boldsymbol{k}},{\boldsymbol{k}}^{\prime}}$ (except for when ${\mathrm{v}}_{\boldsymbol{k}}={\mathrm{u}}_{\boldsymbol{k}}$).

With the above unifying notation, and with $w_{i}=1/n$ and ${\boldsymbol{t}}_{i}$ our lattice points, we can write our approximate coefficient as

[TABLE]

In comparison, the exact coefficient is given by $\widehat{f}_{\boldsymbol{k}}=\langle f,\alpha_{\boldsymbol{k}}\rangle_{\mu}=\langle f\circ\varphi,{\mathrm{u}}_{\boldsymbol{k}}\rangle=\langle f\circ\varphi,{\mathrm{v}}_{\boldsymbol{k}}\rangle$.

6.2 Connection to discrete least squares

With a prescribed ordering of the elements in $\Lambda$, the approximate coefficients (43) for ${\boldsymbol{k}}\in\Lambda$ can be written in matrix-vector notation as

[TABLE]

with column vectors ${\boldsymbol{f}}_{\varphi}=[f(\varphi({\boldsymbol{t}}_{i}))]_{i}$, $\widehat{{\boldsymbol{f}}}^{\square}=[\widehat{f}_{\boldsymbol{k}}^{\square}]_{{\boldsymbol{k}}}$, and matrices $V=[{\mathrm{v}}_{\boldsymbol{k}}({\boldsymbol{t}}_{i})]_{i,{\boldsymbol{k}}}$, $C={\rm diag}(c_{\boldsymbol{k}})$, and $W={\rm diag}(w_{i})$. The conditions (41) and (42) can be expressed as

[TABLE]

with matrices $U=[{\mathrm{u}}_{\boldsymbol{k}}({\boldsymbol{t}}_{i})]_{i,{\boldsymbol{k}}}$ and $D={\rm diag}(d_{\boldsymbol{k}})$.

For plan A (or for the Fourier case) we have $U=V$ and $C=D=I$, so

[TABLE]

which is precisely the solution to the normal equations

[TABLE]

which in turn solves the discrete least-squares problem

[TABLE]

Plan B and plan C do not have the least-squares interpretation.

6.3 Stability to perturbation

Suppose that there is perturbation error in the function evaluations in (43) so that instead of $f(\varphi({\boldsymbol{t}}_{i}))$ we have

[TABLE]

We denote the corresponding perturbed approximate coefficients by $\widehat{f}_{{\boldsymbol{k}},{\rm pert}}^{\square}$ and the corresponding approximate function over $\Lambda$ by $f_{\Lambda,{\rm pert}}^{\square}$ for $\square\in\{a,b,c\}$. Using (43), we can write

[TABLE]

with column vector ${\boldsymbol{\varepsilon}}=[{\varepsilon}_{i}]_{i}$. We have

[TABLE]

where

[TABLE]

so that $\|D^{-1/2}\,V^{*}\,W^{1/2}\|_{2}=1$. Thus

[TABLE]

Here $\rho_{\Lambda}^{\square}$ is the stability constant, and we have

[TABLE]

where $1_{\boldsymbol{0}\in\Lambda}$ is 1 if $\boldsymbol{0}\in\Lambda$ and is [math] otherwise.

For plan A the approximation is perfectly stable.

For plan B and plan C we have the general upper bound $\rho_{\Lambda}^{c}\leq\rho_{\Lambda}^{b}\leq 2^{d-1}$ which might be too pessimistic. If we have a weighted index set with decaying weights (see Example 3) $\rho_{\Lambda}^{b}$ could be much smaller. Alternatively, if $\Lambda$ is downward closed then by Lemma 2 we have $\rho_{\Lambda}^{c}\leq\rho_{\Lambda}^{b}\leq\#{\Lambda}$. In any case, the stability constant for plan B is likely to be much bigger than $1$ even if it is independent of $d$.

For plan C the values of $c_{\boldsymbol{k}}$ depend on the lattice rule and can potentially be as large as $2^{|{\boldsymbol{k}}|_{0}}$, giving hope that one may attempt to minimize the stability constant $\rho_{\Lambda}^{c}$ as part of the CBC construction of the lattice generating vector. Unfortunately, numerical experiments show that not much improvement can be obtained because “self-aliasing” does not happen often enough.

6.4 Error analysis

We have for our three plans A (including the Fourier space), B, and C, annotated by $\square\in\{a,b,c\}$,

[TABLE]

The first part is the truncation error for the finite index set $\Lambda$, and thus represents the best $L_{\mu}^{2}$ approximation error for the choice of $\Lambda$. In the Information Based Complexity (IBC) error analysis, this would be a complexity result using arbitrary linear information: if we know more about the smoothness class of our functions, then this bound is known in terms of a set $\Lambda$ which is constructed according to the decay of the singular values of the approximation operator and this error is exactly the next singular value, see, e.g., [33]. The second part is how well we approximate this best possible approximation by our numerical algorithm which only uses function values; in IBC this is known as standard information, see, e.g., [34].

We proceed to analyze the second error $\|f_{\Lambda}-f_{\Lambda}^{\square}\|^{2}_{L^{2}_{\mu}}$. Since $f_{\Lambda}$ is supported only on $\Lambda$, our reconstruction lattice can exactly compute its coefficients on $\Lambda$. Thus for ${\boldsymbol{k}}\in\Lambda$ we have

[TABLE]

Using (45) and (43) and following the same argument as for the stability analysis, we obtain

[TABLE]

with column vector ${\boldsymbol{g}}=[(f_{\Lambda}-f)(\varphi({\boldsymbol{t}}_{i}))]_{i}$, and we arrive at

[TABLE]

We summarize our combined result for function approximation in the following theorem.

Theorem 27**.**

Consider the problem of approximating a function $f\in{\mathcal{F}}$ from the Fourier, cosine, or Chebyshev space by $f^{\square}_{\Lambda}$, $\square\in\{a,b,c\}$, using a finite index set $\Lambda$ and an $n$-point rank-$1$ lattice under plans A, B, or C as described in this paper. For sufficiently large $n$ we have

[TABLE]

with stability constant $\rho_{\Lambda}^{\square}$ given in (44) and discrete seminorm $\|\cdot\|_{n}$ defined in (46). We have $\rho_{\Lambda}^{a}=1$, and if $\Lambda$ is downward closed then $\rho_{\Lambda}^{c}\leq\rho_{\Lambda}^{b}\leq\min(2^{d-1},\#{\Lambda})$. A loose upper bound is

[TABLE]

The requirement on $n$ to achieve (47) is proportional to $\#(\Lambda\ominus\Lambda)$ for the Fourier case, while for the cosine/Chebyshev case it is $\#({\mathcal{M}}(\Lambda)\oplus{\mathcal{M}}(\Lambda))$ with plan A, $\#(\Lambda\oplus{\mathcal{M}}(\Lambda))$ with plan B, and $\#{\Lambda}\,\#{{\mathcal{M}}(\Lambda)}$ with plan C.

6.5 Comparison with previous results from the literature and conclusions

Approximation by discrete least-squares has been analysed for different measures $\mu$ and bases $\{\alpha_{{\boldsymbol{k}}}\}$ in several works. We mention for instance [2, 3, 21] when using evaluations at random points, and in [22, 25] when random points are replaced by deterministic point sets. A common denominator in all the aforementioned analyses is the equivalence of the norm $\|\cdot\|_{L^{2}_{\mu}}$ and a suitably defined discrete seminorm $\|\cdot\|_{n}$, on the finite-dimensional space ${\mathcal{F}}_{\Lambda}$. More precisely, there exists $\delta\in[0,1)$ such that, under appropriate conditions on $n$, $\#\Lambda$ and $\delta$, it holds that

[TABLE]

Under the same conditions between $n$, $\#\Lambda$ and $\delta$ that ensure (49), the discrete least-squares approximation $\Pi_{n}f$ of any $f\in{\mathcal{F}}$ satisfies

[TABLE]

see [21, Proposition 1] for a proof.

Our results for the Fourier space and for plan A of the cosine and Chebyshev spaces achieve exactly $\delta=0$; see also (48) with $\rho_{\Lambda}^{a}=1$. For the Fourier case we obtain essentially a scaling of $n\geq(\#{\Lambda})^{2}$. For the cosine and Chebyshev spaces we obtain essentially $n\geq 4^{d}\,(\#{\Lambda})^{2}$ in general, and $n\geq\min(4^{d}\,(\#{\Lambda})^{2},(\#{\Lambda})^{2\ln 3/\ln 2})$ for downward closed index sets. However, if the mirrored index set is a weighted hyperbolic cross with sufficiently fast decaying weights (see Example 3), then we obtain essentially $n\geq c_{\tau}\,(\#{\Lambda})^{2\tau}$ for $\tau>1$ arbitrarily close to $1$.

With the Chebyshev space and for any downward closed set $\Lambda$, these results improve on [25] where it is proven that (49) holds true if $n\geq 2^{2d+1}d^{2}(\#\Lambda)^{2}$.

Moreover, we mention that in the case of uniform measure $\mu$ and expansion on the Legendre basis, the results in [22] show a scaling of $n/(\ln n)^{d}$ as $(\#\Lambda)^{4}$ for general downward closed sets, and a scaling of $n/(\ln n)^{d}$ as $(\#\Lambda)^{2}$ when ${\mathcal{F}}_{\Lambda}$ is an anisotropic tensor product space.

Our results for plans B and C in the cosine and Chebyshev spaces do not have the discrete least-squares interpretation. All three plans give exact function reconstruction in ${\mathcal{F}}_{\Lambda}$, but for a general $f\in{\mathcal{F}}$ not finitely supported on $\Lambda$, there is a trade-off between the approximation error and the requirement on $n$ (e.g., plan $A$ requires $n$ to be larger but also has the smallest constant $\rho_{\Lambda}^{a}$). Therefore it is not easy to directly compare them without further analysis.

To proceed further from the very general result in Theorem 27, one would need to make further assumptions on, for example, the smoothness properties of the function space, and the knowledge of a corresponding index set that has been chosen to take advantage of such properties. Starting from the loose upper bound (48), if we know that the best $L^{\infty}$ approximation error satisfies $\|f-f_{\Lambda}\|_{L^{\infty}}\leq c_{q}\,(\#{\Lambda})^{-q}$ for some $q>0$, see, e.g., [16], then we arrive at $\|f-f_{\Lambda}^{\square}\|_{L^{2}_{\mu}}\,\leq\,\sqrt{1+\rho_{\Lambda}^{\square}}\;c_{q}\,(\#{\Lambda})^{-q}$. For the Fourier space we have $\rho_{\Lambda}^{a}=1$ and $n$ needs to be proportional to $\#(\Lambda\ominus\Lambda)\leq(\#{\Lambda})^{2}$, leading to $\|f-f_{\Lambda}^{a}\|_{L^{2}_{\mu}}={\mathcal{O}}(n^{-q/2})$, where the implied constant is independent of $d$ if $c_{q}$ is independent of $d$. For the cosine or Chebyshev space, the result is more complicated because it depends on the size of the mirrored index set. If we have a weighted hyperbolic cross with sufficiently decaying weights (see Example 3) then the mirrored set itself is not of concern. However, for plan B or C we need to further take into account the value of $\rho_{\Lambda}^{b}$ or $\rho_{\Lambda}^{c}$. In general $\rho_{\Lambda}^{b}$ and $\rho_{\Lambda}^{c}$ can be much worse than $\rho_{\Lambda}^{a}=1$, but depending on the actual index set they might also be manageable.

Finally we stress that the $L^{2}$ approximation result based on the estimate (48) is not sharp, and neither is (50), because the best $L^{2}$ approximation error $\|f-f_{\Lambda}\|_{L^{2}_{\mu}}$, i.e., the first term on the right-hand side of (47), has been estimated by the best $L^{\infty}$ approximation error $\|f-f_{\Lambda}\|_{L^{\infty}}$, which is generally half an order worse in the convergence rate (e.g., rate $p$ for $L^{2}$ versus rate $p-1/2=q$ for $L^{\infty}$), see, e.g., [17]. Moreover, a direct analysis on the discrete norm $\|f-f_{\Lambda}\|_{n}$, i.e., the second term on the right-hand side of (47), based on properties of the lattice points, has a chance to improve upon the best $L^{\infty}$ approximation error too. Indeed, function approximation based on rank-$1$ lattices has been analyzed in [5, 14, 15, 18] where the lattice generating vectors were constructed to minimize the approximation error directly, without the reconstruction property. It is known that if $p$ is the rate of convergence for the best $L^{2}$ approximation error (rather than $L^{\infty}$) then lattice generating vectors can be constructed to achieve $\|f-f_{\Lambda}^{a}\|_{L^{2}_{\mu}}={\mathcal{O}}(n^{-p/2})$. There are also other approximation results using rank-$1$ or multiple rank-$1$ lattices, see, e.g., [1, 11, 12, 13, 19, 35, 36].

Rank-$1$ lattices are very attractive due to their simplicity and stability, and the availability of fast computation methods compared to other approximation algorithms.