Trigonometric interpolation on lattice grids

In this paper we construct non-aliasing interpolation spaces and Lagrange functions for lattice grids. We argue that lattice grids are good for trigonometric interpolation and support this claim by numerical experiments. A greedy algorithm allows us to embed hyperbolic crosses in our interpolation spaces, and numerical experiments indicate that lattice grids are at least as good as sparse grids for trigonometric interpolation. A straightforward FFT-algorithm for functions sampled on lattice grids allows for fast computation and good approximation.


Introduction
If N basis functions are needed for a sufficiently accurate approximation in one dimension, then the naive extension to s dimensions is a tensor product of N s basis functions. This exponential growth in terms is a show-stopper for high dimensional computing, but in many applications the expansion terms of a large number of these basis functions Communicated by Lothar Reichel. are insignificant, and can be omitted without serious degradation of the accuracy. For example, the hyperbolic cross approximation space (whose name refers to the shape of its Fourier index set) has been suggested as an alternative [1,21,23] for functions of bounded mixed derivatives.
Sampling grids and approximation spaces come in pairs, in the sense that one desires to be able to go quickly and accurately between function values and expansion coefficients. For hyperbolic cross approximation spaces, a specially tailored version of the FFT-algorithm is available [5][6][7] for sparse grids [1], but Kämmerer and Kunis prove in [9] that the condition number for computing expansion coefficients grows disturbingly with the number of grid points. Numerical experiments also show that the computational overhead is much higher for sparse grid FFT than for regular multidimensional FFT [13].
For trigonometric interpolation, lattice grids have been suggested as an alternative to sparse grids [8,10,11,18]. By a simple variable transformation, sampling on lattice grids permits a standard multidimensional FFT-algorithm for computing the expansion coefficients [18], and the condition number of the matrix of the discrete Fourier transform grows only modestly [10] with the number of grid points.
For any combination of sampling grid and approximation space, the Lebesgue constant is a measure of the quality of this combination. One famous example is polynomial interpolation, which is terrible on the cartesian grid, but excellent on many different grids with grid points clustering near the end points of the interval. In this paper we study trigonometric interpolation on lattice grids in general, propose a conjecture regarding the Lebesgue constant of lattice grids, and conduct experiments indicating the behaviour of the Lebesgue constant for lattice grids and sparse grids. We believe that the modest growth of the Lebesgue constants of lattice grids explains why these grids seem excellent for trigonometric approximation.
In the current section we give the necessary background on multidimensional trigonometric interpolation, and introduce the Lebesgue constant. In Sect. 2 we compute the known Lebesgue constant for the cartesian grid in several dimensions. In Sects. 3 and 4 we discuss interpolation and Fast Fourier Transform on lattice grids. In Sect. 5 we illustrate our findings by numerical experiments.
Let Ω be a set of N grid points x j ∈ [0, 1) s . We are interested in approximating a periodic function f on [0, 1) s by an s−dimensional trigonometric polynomial Here k = (k 1 , k 2 , . . . , k s ) is a multi-index, commonly called wave numbers or Fourier indices, and S ⊂ ZZ s is a finite set with |S| = N . Associated with S is the approximation space H S = {e 2πik·x } k∈S . We write where I denotes the interpolation operator. The expansion coefficients c k are usually found by Fast Fourier Transform. If the system (1.1) is non-singular, the grid is said to be unisolvent with respect to H S . Approximation spaces of interest are: -The tensor product space: Setting all c k = 1 for all k defines an important function.
Another useful tool for interpolation are the Lagrange functions.
If Ω is unisolvent with respect to H S , a set of N Lagrange functions can always be found, and the interpolation can be written . A basic question is how well the interpolant approximates the interpolated function. Let f * be the best approximation of f onto H S . Since f * ∈ H S , we have and thus the norm of the interpolation operator tells us how far I f is from the best approximation of f in H S . From the definition of operator norm it follows that This motivates the definition of the Lebesgue constant.

Dirichlet kernels and Lebesgue constants on equidistant grids
In this section we compute Lebesgue constants for some equidistant grids. These are known results, but we include them for completeness, and for the intuition these proofs provide regarding our conjecture in Sect. 3. It was shown by Ehlich and Zeller [4] that the Lebesgue constant for trigonometric interpolation at N equidistant nodes equals the Lebesgue constant for standard polynomial interpolation at the Chebyshev points (the zeroes of the Chebyshev polynomials), which are known to grow as 2 π log N + C where C is a small constant [14].

Theorem 2.1 For trigonometric interpolation on the N
and a bit of elementary calculus it follows that Since D n (x) is 1-periodic, and D n (x − x j ) is defined by a regular shift of x j = j/(2n + 1), the sum n j=−n |D(x − x j )| is 1 2n+1 -periodic, and we can write For j = 0, D n (0) = 2n + 1 is the maximum value. For j = 0, bounds can be established by writing and using (2.1) evaluated at the interval endpoint closest to the origin, We can now write In general it is difficult to find closed form expressions for the Dirichlet kernel in several dimensions, but it is easy to show that D T s d factors into s univariate Dirichlet kernels. Lagrange functions can thus be constructed for the Cartesian grid by scaling and translation, and an application of the above theorem gives the following corollary.

Interpolation on lattice grids
An s-dimensional lattice Λ is simply the integer linear combinations of s linearly independent basis vectors. Arranging these vectors as rows in a matrix A, we have The matrix A is called a generator matrix for Λ. If Λ is periodic on [0, 1) s , the lattice will be denoted an integration lattice [19]. Associated with every lattice is the dual lattice, Λ ⊥ , defined by The dual lattice is itself a lattice, and has B = (A −1 ) T as a generator matrix. If Λ is an integration lattice, then Λ ⊥ is an integer lattice, B is an integer matrix [15], and the number of lattice points in [0, 1) s is N = | det(B)|. An alternative form for describing an integration lattice is by its tZ -form [16]. Let By the notation {x} we understand the fractional part of x. With this form the lattice points in Λ ∩ [0, 1) s can be written t coincide with the rank of the lattice rule [16]. Consequently if N is prime there is a tZ-form with t = 1.
Two Fourier modes e 2πik·x and e 2πih·x will be indistinguishable on the lattice grid points if k − h ∈ Λ ⊥ . This phenomenon is called aliasing.

Lemma 3.1 Let Λ be an integration lattice with N points in
By the definition of dual lattices, We have to prove that this implies (k − h) ∈ Λ ⊥ . Since N is prime, z/N may be used as generator for Λ (this is not true for any lattice point z/N unless N is prime), and the lattice points in [0, 1) s may be written is an integer for all j, and hence k − h ∈ Λ ⊥ .
For a function f with convergent Fourier series it follows that for x j ∈ Λ and S a non-aliasing index set, we have: This generalizes the standard result for trigonometric interpolation to functions sampled on lattice grids.  The following theorem proves that if N is prime, we can use shifted Dirichlet kernels as Lagrange functions for interpolation on the corresponding lattice grid.

Theorem 3.1 The trigonometric interpolating polynomial on a non-aliasing index set S is
Proof Since any linear combination of lattice points is also a lattice point, we have x j − x l = z/N , where z ∈ ZZ s . Thus k · (x j − x l ) = k · z/N = m N for all k ∈ S, and we can write The lower entry on the right is a finite geometric sum, since the complex exponentials take at most N different values (by construction N = |S|), and none of the m's are equal (mod N) (due to Lemma 3.1).

Remark 3.1
This theorem holds for general lattices, but cannot be proved with Lemma 3.1, since this does not hold unless N is prime. The general proof requires an approach based on group theory. We will describe this approach and present the general version of the theorem in a forthcoming paper.

Corollary 3.1 For a given integration lattice Λ and a non-aliasing index set, a complete set of Lagrange functions for trigonometric interpolation on H S is
The expression for the Lebesgue constant can now be written Since every Lagrange function is a shifted Dirichlet kernel, computing H is equivalent to summing the local maxima of one Dirichlet kernel over the N unit cells of the lattice.
where U is the parallellotope defined by 2 s adjacent lattice points. The last sum can be interpreted as an upper Riemann sum of the integral Similarly, we conjecture that the total volume of N 'hyperpyramids' with base U is a lower bound for ||D S (x)|| 1 , leading to the following bounds on H In Travaglini [22] proved the following bounds on the Dirichlet kernel for convex polyhedral interpolation spaces Here C 1 , C 2 are constants depending on s, but not on N . Combining (3.3) and (3.4) establishes the following conjecture.

Conjecture 3.1 The Lebesgue constant for trigonometric interpolation on a lattice grid Λ onto an associated non-aliasing interpolation space H S , is bounded bỹ
We now proceed to the practical contruction of non-aliasing index sets for integration lattices. Note that the interpolation space, as defined by S O , is not unique, since the choice of generator matrix is not unique. If U is a unimodular matrix, then U B generates Note that the transposed lattice Λ T of Theorem 3.3 changes with B.

Theorem 3.3 Let Λ be the integration lattice generated by A, Λ T generated by A T and B = (A
As illustrated in the figure, S-domains obtained using Theorem 3.3 are always parallellotopes, being affine mappings of the unit cube. If one wishes instead to mimic the spaces listed in Sect. 1, a simple greedy algorithm can be utilized to produce an interpolation space H Λ containing any specific space as a sub-space. Given a lattice, choose N non-aliasing points k ∈ ZZ s of increasing order. If we want S to mimic a hyperbolic cross we use the "Zaremba"-index ρ(k) = s i=1 max(1, |k i |) as the ordering criterion. Choosing the 1-norm ||k|| 1 = s i=1 |k i | as ordering criterion, gives an index-space resembling P s d . These two ordering criterions are illustrated in Fig. 2, where we display two different versions of S for the lattice of Fig. 1. None of these sets are unique, as there are many integer points having equal Zaremba index or 1-norm.

FFT on lattice grids
On equidistant s-dimensional Cartesian grid, the standard multi-dimensional FFTalgorithm provides fast and stable computation of the Fourier coefficients. On lattice grids, the FFT-algorithm is readily available, see [18]. We include it for completeness. Utilizing the t Z-form, the Fourier coefficients can be written  The expression for f (m) is a standard t-dimensional DFT, straightforwardly computed by your favorite FFT-algorithm. The components are identified by the m-index which are related to the k-index by the relation (4.1).

Numerical experiments
In this section we present some computational results in 2D and 3D. In Fig. 3 we show an estimate of the Lebesgue constants for lattice grids compared to sparse grids as a function of interpolation points. An estimate for the Lebesgue constant is obtained by evaluating the Lagrange functions on a fine grid. In 2D we have used the rank-1 lattices generated by z = [1, 2d − 1] and N = 2d 2 ; d = 2 : 2 : 24. These are known to be optimal lattices for the space P s d [2]. In 3D we have used lattices optimized for the space P s d , based on skew-circulant generator matrices [17]. For the computation of sparse grids and hyperbolic spaces we have used the matlab package nhcfft-0.91 freely available from TU-Chemnitz [3]. We note that in 2D as well in 3D the Lebesgue constant is significantly lower for lattice grids than for sparse grids. The Lebesgue constant seems for both grids to grow as C s (log n) s , with lattice grids having significantly lower constant C s .
To test the quality of the interpolation we have used the following test functions: sin 2π x j p ; p = 1, 2, 3.
The first two are found in [18], while the family g p is a scaled version of the family used for testing sparse grid interpolation in [6]. Note that in each direction the periodic extension of Since the asymptotic behavior of the Fourier coefficients of a one dimen- This indicates that the hyperbolic cross is a good approximation space for product functions in C p (R), and we expect the sparse grid to work very well in this case.
When f ∈ C ∞ (R) the decay rate for the 1-dimensional Fourier coefficients is exponential: |c k | ∼ a |k| ; 0 < a < 1. The corresponding s−dimensional coefficient is computed by evaluating the exact value of f i as well as the approximation on a fine, regular Cartesian grid. When comparing these results we have to keep in mind that the interpolation spaces are slightly different for the two grid types; for sparse grids we have used a dyadic hyperbolic cross, while for lattice grids we have used spaces which embed the largest possible symmetric hyperbolic cross. Nevertheless we find these results encouraging; lattice grids seem to perform similarly to, or better than, sparse grids for these examples, reflecting the fact that it has significantly lower Lebesgue constant. Where lattice grids really outperforms sparse grid is for f 2 . As explained above we contribute this to the fact that this function is in C ∞ (R). This is the same effect as seen for the test function g p . Increasing p means increasing the smoothness which should favor lattice grids. Our brute force approach for estimating the error and the Lebesgue constant becomes very expensive for large numbers of lattice points. This limits our possibilities for extensive experiments in higher dimensions. We have run some limited experiments where we have kept the number of lattice points constant and compared it with sparse grid interpolation of approximately the same size. The lattice rules we have used are rank-1 lattices constructed by the component-by-component technique introduced by Sloan and Reztsov [20]. Our lattices are those reported by Kuo et al. [12]. We have also restricted us to test function f 1 only. The results are reported in the table above. The most remarkable difference is the minor increase in Lebesgue constants for lattice grids, as opposed to the disturbing growth for sparse grids. We believe this to be the main reason why the relative error is smaller for lattice grids than sparse grids.

Conclusion
In this paper we have constructed non-aliasing interpolation spaces and Lagrange functions for lattice grids. Our family of approximation spaces does not include the hyperbolic cross, which is the natural Fourier index set for approximating functions with bounded mixed derivatives, but a simple greedy algorithm allows us to embed the hyperbolic cross as a sub-space in our interpolation spaces.
Both lattice grids and sparse grids seem to have quasi optimal Lebesgue constants.The quality of the lattice interpolation appears to be better than that of sparse grid interpolation, especially in several dimensions, as measured by estimation of the Lebesgue constant as well as in approximation of test functions.