Fast Algorithms for Periodic Spline Wavelets on Sparse Grids.
We consider Boolean sums of univariate interpolation operators which define multivariate jth order blending interpolation operators on sparse grids. Sample spaces are defined as range of the blending operators. Sample and wavelet spaces have significantly lower dimension and good approximation order for certain function spaces. Fast decomposition and reconstruction algorithms for bivariate spline wavelets, based on algorithms for univariate functions, are described. Operation counts for the algorithms are given and it is shown that the complexity depends linearly on the dimension of sample spaces.
Introduction
. Wavelet analysis has evolved in the last few years as a powerful
method in numerical analysis and signal processing. An important question is the
development of efficient algorithms for multivariate wavelet transforms. The classical
approach for the construction of bivariate wavelets from given univariate functions by
tensor products leads to a fast increasing dimension of space of order O(2 2j ). There-
fore, one has to take into account a high complexity of the corresponding algorithms.
To avoid this problem the concept of Boolean sums and sparse grids is often used
in multivariate approximation and interpolation (see e.g. Delvos and Schempp [5]) as
well as in finite element methods (see e.g. Zenger [16]). Applying the Boolean sum
method sample spaces of dimension O(j2 j ) are described in [1] and [5]. It is shown
there that for functions from bivariate Korobov spaces (see [7]) one obtains only an
extra logarithmic factor in the approximation order, while the dimension of the sample
spaces is significantly lower. A unified approach to error estimates for Boolean sums
and tensor products in Korobov and Sobolev spaces can be found in [10].
The Boolean sum method is applied to wavelet theory by Sprengel in [14], where
corresponding sample and wavelet spaces of bivariate periodic functions are consid-
ered. A basis for sample spaces is given by fundamental interpolants on sparse grids.
It is crucial for this approach, that the underlying univariate wavelets and scaling
functions are fundamental interpolants. In [15], a first algorithm for complete decomposition
and reconstruction with complexity O(j described by using matrix-vector
notations and the fast Fourier transform (FFT).
The aim of this paper is to describe faster algorithms by using periodic spline
wavelets with minimal support. Therefore, we establish in x2 an approach, where
the underlying univariate wavelet and scaling functions need not be fundamental
interpolants. We demand only, that for the sample spaces a certain interpolation
problem has a unique solution. Basis functions for the bivariate sample spaces are
still the fundamental interpolants from [14], for which we establish the two-scale
relation in Theorem 2.1. However, for our applications, we mainly use a generating
system which consists of the translates of certain tensor products of univariate scaling
functions. This approach leads to algorithms, which are different from the algorithms
of Sprengel [15] and have some particular advantages.
Fachbereich Mathematik, Universit?t Rostock, Universit?tsplatz 1, D-18051 Rostock, Germany,
K. BITTNER
For our algorithms, we use periodic spline wavelets with minimal support as
described by Plonka and Tasche [8]. In x3, we give a short introduction to these
wavelets. In x4, we study fast algorithms for univariate periodic spline wavelets,
which are based on a method described in [2]. The idea of recursive filters used there
yields complexity O(2 j ) in the univariate setting.
In x5, we present new algorithms for wavelets on sparse grids, which do not require
the complete decomposition of a function until the lowest level as in [15] and
therefore improve the algorithms presented there. The main idea is to use the representation
of functions from the sample spaces by the generating system mentioned
above. In this way, we can apply ideas from the tensor product setting, such that
our algorithms make use of the algorithms for the univariate case from x4. Beyond
the decomposition and reconstruction algorithm, we give also an interpolation algorithm
and an algorithm for the computation of function values on the sparse grid.
For all algorithms, we count the floating point operations and compare it with the
corresponding numbers of the tensor product algorithms.
A numerical example is presented in x6. For this we compare the results of
Boolean sum and tensor product algorithms. We describe some special effects of the
Boolean sum method and show that for suited functions the Boolean sum method is
more efficient than tensor product methods.
2. Sample and Wavelet spaces. Bivariate interpolating scaling functions and
wavelets for sparse grids have been investigated for the first time by Sprengel in [14, 15]
using Boolean sums of univariate interpolation operators (see [1, 5]). In the following
we shortly recall this definitions and establish a more general approach for wavelets
on sparse grids.
By C(T n ) we denote the space of all continuous, n-variate 2-periodic functions
and by L 2 (T n ) the space of square-integrable, n-variate 2-periodic functions with
the inner product hf; gi T n := (2-) \Gamman
R
dx.
For a fixed d 2 N, we put d j := d . For each j 2 N 0 , we introduce
univariate scaling functions OE j;k as translates of a single function, i.e.,
In the following, we demand that the scaling functions OE j;k are linear indepentent.
Then we define corresponding univariate sample spaces by
Further, we demand the existence of a linear interpolation operator L j with the
precision set
prec
ae
oe
ae [0; 2-)
and the image ImL there have to exist
coefficients a j;k 2 R such that
a j;k OE j;k (t j;l
Obviously, the interpolation projectors L j satisfy the embedding relations
PERIODIC SPLINE WAVELETS ON SPARSE GRIDS 3
For the bivariate setting, we introduce the projector of j-th order blending interpolation
defined by
where A \Phi B := A+B \Gamma AB denotes the Boolean sum and
A\Omega B the tensor product
of projectors (cf. [5]). Then, the interpolation knots form a sparse grid (Figure 2.1)
given by
prec L r \Theta prec L
Fig. 2.1. Sparse grids T Bfor
Now we define the bivariate sample spaces as the images of B j , i.e.,
ImL
A generating system of V B
j is given by
This system consists of the translates of different tensor product functions and
therefore it is very useful in applications. But for the unique representation of a
a basis is necessary. To construct such a basis of V B
j we employ
the interpolation properties of B j . For this we consider the univariate fundamental
interpolants OE ffi
are uniquely determined (cf. (2.2)) by
and which form a basis of V j . Analogously, a basis of V B
j is given by the bivariate
fundamental interpolants for the sparse grid T B
4 K. BITTNER
and OE j;s j 0 for
2 Z. These
functions are not translates of each other and have a more difficult structure. But
it turns out, that for our algorithms, we need not use the functions OE
j;k and \Phi j;(k;l)
explicitly. However, it should be useful to investigate the functions \Phi j;(k;l) in more
detail. Therefore, we will establish and prove a refinement relation for these functions.
For this, we consider the following finer splitting of J j
K -;r
into the disjoint sets K -;r with K
For
s=r
s=r
and in particular for (k; l) 2 J j simplifies to
Now we establish the refinement equation of the bivariate fundamental interpolants

Theorem 2.1. Let (k; l)
we have
Proof. From the interpolation properties of OE ffi
j;k , we conclude
(2.
PERIODIC SPLINE WAVELETS ON SPARSE GRIDS 5
Let now j, k, l be fixed. We introduce the notations
oe 4
oe 5
For we have with (2.4)
s=r
s=r
From (2.7), it follows that
OE s;2 s\Gammaj
ds
ff s;m\Gamma2 s\Gammaj k OE s+1;2m+1
s
and
OE s;2 s\Gammaj
ds
ff s;m\Gamma2 s\Gammaj k OE s+1;2m+1
With (2.4), we obtain
s=r
oe 5
Taking into account (2.8) and (2.5), the assertion (2.6) follows.
Now we introduce the univariate wavelet space W j as the orthogonal complement
of denotes the orthogonal sum in L 2 (T n ).
We choose a wavelet function / j 2 W j , such that
6 K. BITTNER
is a basis of W j .
The bivariate wavelet spaces are given by
In [14], it has been proved that V B
. A basis of W B
is given in a natural way by
f/
By counting the basis functions we obtain
In [14], Theorem 2, it is shown that the bivariate wavelet basis (2.9) is stable, if the
univariate scaling functions and wavelets form stable bases.
3. Periodic spline wavelets. Now we will consider an example for suitable
univariate sample and wavelet spaces. By Sprengel [14, 15], interpolating trigonometric
polynomials as described in [11, 12] are used as univariate wavelet and scaling
functions. In this case one has OE
j;k and the wavelets satisfy an interpolation
property, too. This has some advantages in the notation. On the other hand, this
approach leads to algorithms with complexity O(j2 j ) such that the algorithms for
the bivariate case have complexity O(j To improve the algorithms we propose
another approach. Because of the wide applicability and the simple structure of spline
interpolation, we will restrict ourselves to spline wavelets in the sequel.
Here, we consider periodic spline wavelets with minimal support as described in
[8, 9]. We briefly recall some definitions. The k-th Fourier coefficient of f 2 L 2 (T) is
given by c k (f) := hf; e ik T Z. Further, we define the discrete Fourier transform
of length N by - xn :=
\Gamma2-i=N is a primitive N -th
root of unity.
The cardinal B-spline of order m - 1 is defined by
with
By the convolution property, the Fourier transform of
Nm reads as follows
Nm (u) :=Z
The scaling function OE j (j 2 N 0 ) is defined as 2-periodization of Nm (2 j \Delta), i.e.
Nm
By the Poisson summation formula, one obtains the Fourier coefficients
Nm
and c n (OE j;l
The two-scale relation is given by OE
or in the Fourier transformed
d j+1
with the two-scale symbol
Using the results of Locher [6] for construction of periodic fundamental interpolants
in shift invariant spaces, we obtain the Fourier coefficients of the fundamental
j;0 of even spline order
with the m-th Euler-Frobenius function \Phi m (z) :=
z k . In

Figure

3.1, we show bivariate basis functions constructed from these splines. For odd spline
order, there does not exist a unique spline interpolant on the grid T j . Consequently,
we will deal in the following only with splines of even order
Fig. 3.1. Bilinear
The wavelet / j is given as the 2-periodization of the the cardinal spline wavelet
introduced by C. K. Chui and J. Z. Wang [4, 3]. Therefore, we define /
by its two-scale relation
or in the Fourier transformed form
d j+1
with the two-scale symbol

Figure

3.2 shows bivariate wavelet functions constructed from the periodized Chui-
Wang wavelets.
8 K. BITTNER
-0.20.2
Fig. 3.2. Bilinear
4. Fast univariate spline algorithms. Fast algorithms for spline interpola-
tion, decomposition, reconstruction and computation of function values based on the
FFT are described in [8]. Unfortunately this algorithms have complexity O(j2 j ). Our
goal is to apply algorithms of complexity O(2 j ) to achieve the complexity O(j2 j ) in
the bivariate case.
Because we consider periodic functions we can assume that the coefficients are
extended periodically, e.g. a . For the reconstruction it is well known,
that in
a j+1;k OE
a j;k OE j;k
the coefficients a j+1;k can be computed by
a
a
with 4m d j real additions and multiplications. In the following we will count one real
addition and one real multiplication together as one real operation.
To obtain a decomposition algorithm of complexity O(2 j ) we refer to Berger and
[2], where an algorithm for fast decomposition based on recursive filters
is described. Let - m;- ! 0, the zeros of the Euler-Frobenius
polynomial inside the unit circle, which are all real and negative (see [13]). Further,
we consider the coefficients r k and s k of the numerators of the decomposition symbols
PERIODIC SPLINE WAVELETS ON SPARSE GRIDS 9
With these notations the decomposition algorithm of Berger and Str-omberg,
which determines the coefficients a j;k and b j;k in (4.1) from given a j+1;k reads as
follows.
Algorithm 4.1.
spline dec
1. For
l=1\Gamma2m
r l a j;2k\Gammal ; v (1)
l=1\Gamma2m
s l a j;2k\Gammal
2. For
(a) u (2-)
2m;- u (2- \Gamma1)
(b) For
(c) u (2-+1)
2m;- u (2-)
\Gammas
(d) For
3. For
Return
The number of real operations yields as follows. Let us assume that the values of
2m;- are precomputed. For the first step, we need 4m d j operations. In the second
step, for each - we need 2(d operations for (a) as well as for (c) and 2(d
operations for (b) and (d) each. Therefore, the second step needs (8m\Gamma8)d j operations
at all. Together these are (12m \Gamma 8)d j operations for the decomposition algorithm.
Furthermore, we have to consider an interpolation algorithm. For given function
values y we want to determine the coefficients a j;k such that
a j;k OE j;k (t j;l
In [8], the Fourier transformed solution for even spline order
- a
Because the denominator of the symbol w \Gammanr
is similar to the denominator
of the decomposition symbols, we can apply again the idea of Berger and Str-omberg,
which leads to the following interpolation algorithm.
Algorithm 4.2.
spline int
1. For
2. For
(a) u (2-)
(b) For
(c) u (2-+1)
\Gammas
(d) For
3. For
Return (a j;k )
k=0 .
For this algorithm, we count d j operations in the first step and (4r \Gamma 4)d j operations
in the second step. Hence, the absolute number of real operations for this
algorithm is (2m \Gamma 3)d j .
5. Algorithms for bivariate spline wavelets on sparse grids. In this section
we will describe algorithms for the decomposition of a bivariate function
j+1 into a scaling part f
j and a wavelet part
j as well as
for the reconstruction of f j+1 from f j and g j . The main difficulty of the Boolean sum
approach is the structure of the basis functions \Phi j;(k;l) . That is why we propose the
generating system (2.3) for our algorithms. In this way, we can use the ideas from
the tensor product case. Thus, we can establish fast algorithms, which have a simple
structure and are based on the univariate algorithms. The algorithms are independent
from the underlying univariate wavelets. For simplicity and because we want to give
exact operation counts we will restrict ourselves to spline wavelets.
5.1. Interpolation algorithm for sparse grids. To obtain a representation
of a function f
j as a linear combination of elements from the generating system
(2.3) we have to compute coefficients a r;(k;l) with
r;k;l a r;(k;l) OE
From that the coefficients a r;(k;l) are a solution of the following
interpolation problem. For all (-) 2 J j
dr
a r;(k;l) OE r;k (t j;- )OE j \Gammar;l (t j;-
Because we have not a basis the interpolation problem has several solutions a r;(k;l) ,
which represent all the same function. To obtain one of these solutions we will first
describe an interpolation algorithm for sparse grids. With
f j;r := L
dr
a
Now we write
gr (\Delta;y)
PERIODIC SPLINE WAVELETS ON SPARSE GRIDS 11
This means for each l we have to interpolate the values
Applying the univariate interpolation algorithm for y, we obtain for
and
l=0 u j;(k;l) OE j;l (y) with
By the univariate reconstruction formula (4.2), we compute coefficients v r;(k;l) , such
that
Finally, we have
a
and therefore
(a r;(k;l)
(v r;(-;l)
Hence, the spline interpolation algorithm for sparse grids reads as follows.
Algorithm 5.1.
bool int
1. For
For
2. For
For
For
3. For
For
(a r;(k;l)
(v r;(-;l)
Return (a r;(k;l) ) r;k;l .
In this algorithm, we need
real operations in the first as well
as in the third step and
operations in the second step. This
are altogether
operations.
5.2. Decomposition. Now we want to decompose a given function f
into functions f
we have to determine
coefficients a r;(k;l) and b r;(k;l) such that
dr
c r;(k;l) OE
r;k\Omega OE j+1\Gammar;l
dr
a r;(k;l) OE
for given c r;(k;l) . Applying the univariate decomposition algorithm first with respect
to y and then with respect to x we obtain
dr
~
j+1;k\Omega OE 0;l
dr
~ u r;(k;l) OE
\Gamma1;k\Omega / j;l
By the univariate reconstruction formula (4.2), we can then compute the values a r;(k;l)
such that
dr
a r;(k;l) OE
From this equalities, we derive the bivariate decomposition algorithm, which determines
the coefficients a r;(k;l) and b r;(k;l) .
Algorithm 5.2.
bool dec
1. For
For
spline dec
For
2. For
For
For
spline dec
3. For
For
a
a
For
Return
This algorithm needs
real operations for the first,
operations for the second step and finally
operations for the third step. Therefore, the decomposition
algorithm needs
operations.
5.3. Reconstruction. Analogously, we can derive a reconstruction algorithm,
which determines the coefficients c r;(k;l) from given coefficients a r;(k;l) and b r;(k;l) in
(5.1).
Algorithm 5.3.
bool rec
1. For
For
For
2. For
For
For
Return
14 K. BITTNER
This reconstruction algorithm needs
operations
in the first step and
operations in the second step. Summing
up, these are
real operations.
5.4. Computation of function values on the sparse grid. Finally, we consider
the question how to compute function values of f
j . In particular, we have
to determine the values for all points of the sparse grid T B
j to obtain the coefficients
Because of supp OE function values of
dr
a r;(k;l) OE
can be computed for (x; operations by
l=-
a r;(k;l) OE r;k (x)OE j \Gammar;l (y):
To determine the sum f(x;
need then In
this way, the computation of all coefficients fl k;l , (k; l) 2 J j would require O(j
operations. To reduce the complexity further to O(j2 j ) we will use the refinability
of OE j .
First we introduce index sets I j;r with I j;0 := f2
I j;r := f2 j \Gammar (2k j. By this definition we
have
I j;r and
Now we consider for I j;r the values
s=r
s=r
where
ds
k=0 a s;(k;l) OE s;k (t s;- ). By the reconstruction formula (4.2) we compute
coefficients v -;l , such that b
Analogously, we determine
the values c - :=
I j;r .
I j;r and - 2 I j;-\Gammar and we can
determine the values of f by
To obtain the function values OE j (t j;k ) we use the definition (3.2), such that we
only need the values Nm (k) of the cardinal B-spline for
the function values y on the sparse grid can be computed with the
following algorithm.
Algorithm 5.4.
sparse val
1. For
For
k=1 a r;(-\Gammak;-) Nm (k)
l=1 a r;(- \Gammal) Nm (l)
r;(- \Gammal) Nm (l)
2. For
For
For
3. For
s=1 u (1)
4. For
Return (y -)2J j .
K. BITTNER
Now we count the number of real operations for this algorithm. Obviously, the
first step needs d 2 (j operations. The number of operations for the
second step is given by/
d
d
Analogously, we obtain d 2 (j operations in the third step. In the fourth
step we have only real additions. Because each value w r;(-) is used once we obtain
the number of this additions by
dr
Therefore, we can estimate the number of real operations for this algorithm by
In the following table we sum up our operation counts for the Boolean sum algorithms
and compare it with the operation numbers for the corresponding tensor
product methods.
Algorithm Boolean sum Tensor product
Interpolation
Decomposition (39mj
Reconstruction
Function values (5m
Hence, the reconstruction and decomposition algorithms for spline wavelets on sparse
grids have the arithmetical complexity O(j2 j ), while the corresponding algorithms
for tensor products have complexity O(2 2j ).
Remark 5.5. From the structure of the sparse grids one would expect some
problems in the implementation of the algorithm. However, because of the simple
structure of the generating system (2.3) this problem does not occur in the Algorithms
5.2 and 5.3 for decomposition and reconstruction. In Algorithm 5.4, the second, third
and fourth step show how one can deal with the sparse grids using the index sets I j;r
and K -;r without a considerable loss of performance. Analogously, we can handle this
problem in step 1 of Algorithm 5.1. Therefore, the number of floating point operations
is suited for a comparison of tensor product methods and our Boolean sum algorithms.
6. Numerical examples. Finally, we present the results of numerical tests. The
detection of higher order discontinuities is an important property of wavelets. The
reason for this is, that the wavelet part of a function is the error of the orthogonal
projection, which depends on the smoothness of the function such that the large
wavelet coefficients are localized near the discontinuities.
Therefore, we consider first the approximation properties of the interpolation
projectors. Let f 2 C r;r (T 2 ) with
PERIODIC SPLINE WAVELETS ON SPARSE GRIDS 17
In [5], it is shown that for r - m the error of spline blending interpolation on sparse
grids is
while the error of tensor product interpolation on a regular grid is given by
That means, we obtain for Boolean sums only an extra logarithmic factor in approximation
error, while the cost for the algorithms is remarkably lower. On the other
hand, if
it can be shown, that for B j it holds that
In this case, the approximation order of the tensor product interpolation operator is
For our tests, we choose bivariate spline functions, because this functions have
discontinuities in its derivatives at well known lines. Therefore, we interpolate and
decompose the tensor product spline N
4\Omega N 4 as well as the box spline B (222) both
centered in the point (-) and periodically extended. Both functions are contained
in C 2 (T 2 ). On the other hand, N
4\Omega N 4 is contained in C 2;2 (T 2 ), too, whereas B (222) is
only contained in C 1;1 (T 2 ) but not in C 2;2 (T 2 ). With these examples, we demonstrate
for which kind of functions the Boolean sum method works well. The error estimates
given above let expect good results of the Boolean sum algorithm only for the tensor
product spline. We choose the underlying spline wavelets of order to rule out
effects from the spline structure of our test functions. Furthermore, we set
because with this setting the special properties of the Boolean sum approach should
appear most visible.
We compare the results of our algorithms with the tensor product methods.
For tensor products, we decompose the functions from level 9 to level 8. In the
Boolean sum case we choose once the decomposition from level 15 to 14. Because of
dim(V
15 we have around the same amount of memory in this way. On
the other hand, we perform the decomposition from level 11 to 10. The results show
that for this level the approximation error for N
4\Omega N 4 is roughly the same as in the
tensor product approximation to level 8.

Figure

6.1 shows density plots of the absolute values of the wavelet part of the
functions. For higher values the shade is darker. As expected the result for the tensor
product spline looks better than for the box spline. The pictures for the box spline
show, that only discontinuities on lines which are diagonal to the axes are detected.
The reason for this is that the function is not smooth here in the sense of C 2;2 (T 2 ).
Furthermore, for the Boolean sum approach we have other localization properties,
namely for the detection of singularities of derivatives in one direction we have good
localization in this variable and bad localization in the other variable. On the other
hand, for the detection of singularities in mixed derivatives the localization is equal
in both directions.
K. BITTNER
Fig. 6.1. Density plot of absolute values of the wavelet part for different decompositions. The
numbers in the brackets show the plot range (minimum=white; maximum=black).
PERIODIC SPLINE WAVELETS ON SPARSE GRIDS 19

Table
Number of wavelet coefficients which are greater than " for different decompositions.
tens. prod.
level 8
Bool. sum
level 14
Bool. sum
level
tens. prod.
level 8
Bool. sum
level 14
Bool. sum
level
" total % total % total % total % total % total %
In

Table

6.1, we see the number of wavelet coefficients with an absolute value
greater than some threshold ". Obviously, for the tensor product spline N
4\Omega N 4 the
Boolean sum method is more efficient than the tensor product algorithms. For the
box spline, the results of the Boolean sum method are similar as in the tensor product
case.
100.
CPU-Time in s
Fig. 6.2. CPU-time (on a Sun SPARCstation 20) of tensor product algorithm (dashed line)
and Boolean sum algorithm (continuous line) for
These examples show that the concept of Boolean sum seems to be more suitable
for function spaces of the type C r;r (T 2 ) whereas the tensor product methods appear
K. BITTNER
to be better for the spaces of the type C r (T 2 ). Therefore, a possible application for
wavelets on sparse grids could be seen in the numerical treatment of PDE's, where
the solution is contained in C r;r (T 2 ) but not in C r+1 (T 2 ). For such applications, it
should also be useful to investigate a Boolean sum approach for non-periodic wavelets
on the square.
Finally, we draw in Figure 6.2 a comparison of the time needed for tensor product
and Boolean sum algorithms. It shows that the Boolean sum method is remarkably
faster than tensor product algorithms with the same quality of results for suited
functions.

Acknowledgments

. The author thanks the referees and the editor for their
helpful comments and valuable suggestions to improve the manuscript. The research
of K. Bittner was supported by the Deutsche Forschungsgemeinschaft. The author
thanks J-urgen Prestin for suggesting the problem and the valuable comments.



--R

A discrete Fourier transform scheme for Boolean sums of trigonometric operators.

An Introduction to Wavelets
On compactly supported spline wavelets and a duality principle
Boolean methods in interpolation and approximation
Interpolation on uniform meshes by translates of one function and related attenuation factors
The approximate computation of multiple integrals
On the computation of periodic spline wavelets
A unified approach to periodic wavelets


Interpolatory and orthonormal trigonometric wavelets

Multivariate periodic interpolating wavelets
Periodic interpolation and wavelets on sparse grids
in Notes on Numerical Fluid Mechanics
--TR
sparse grids;multivariate periodic interpolation;wavelets;splines;boolean sums