# Inner product computations using periodized daubechies wavelets

Juan Mario Restrepo, Gary Leaf

► To cite this version:

Juan Mario Restrepo, Gary Leaf. Inner product computations using periodized daubechies wavelets. International Journal for Numerical Methods in Engineering, 1998, 40 (19), 10.1002/(SICI)1097-0207(19971015)40:193.0.CO;2-A. hal-01323007

HAL Id: hal-01323007

https://hal.science/hal-01323007

Submitted on 29 May 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
INNER PRODUCT COMPUTATIONS USING
PERIODIZED DAUBECHIES WAVELETS

JUAN MARIO RESTREPO<fnref n="1" />* AND GARY K. LEAF<fnref n="2" />

Inner products of wavelets and their derivatives are presently known as connection coefficients. The numerical calculation of inner products of periodized Daubechies wavelets and their derivatives is reviewed, with the aim at providing potential users of the publicly-available numerical scheme, details of its operation. The numerical scheme for the calculation of connection coefficients is evaluated in the context of approximating differential operators, information which is useful in the solution of partial differential equations using wavelet-Galerkin techniques. Specific details of the periodization of inner products in the solution differential equations are included in this presentation. [wavelets; Galerkin; differential operators; connection coefficients]

# 1. INTRODUCTION

Wavelets have found a well-deserved niche in such areas of applied mathematics and engineering as approximation theory, signal analysis, and projection techniques for the solution of differential equations. While wavelets are not conceptually new,<fnref n="1" />⁻³ the past fifteen years have produced much of the theoretical underpinnings for the concept, as well as the generation of new wavelet families and the exploration of their potential in various areas of applied science (see References 4 and 5 for wavelets in signal processing, References 6 and 7 for wavelet applications in approximation theory, References 8 and 9 for wavelet applications in the solution of differential equations, and Reference 10 for other types of applications).

Wavelets have several advantages: (1) they have compact support or exponentially decaying support; (2) their continuity properties may easily be increased, albeit at the expense of a larger

<footnote n="1">¹Mathematics Department, University of California – Los Angeles, Los Angeles, CA 90095, U.S.A.</footnote>

<footnote n="2">² Mathematics and Computer Science Division, Argonne, National Laboratory, Argonne, IL 60439, U.S.A.</footnote>
domain of support; (3) for a given spline order, an complete basis may easily be generated by simple recurrence relations; (4) in the context of projection techniques, their convergence properties are as good as or better than Fourier methods, and they permit the analysis of extremely local functional behaviour without the need for windowing and with little or no bias from global behaviour; and (5) the manner in which the space is broken down into a family of multiply-enclosed subspaces enables spatial or temporal function multiresolution analysis.

A multiresolution analysis is a nested sequence

$$V_0 \subset V_1 \subset \cdots \subset L^2(\mathbb{R})$$

satisfying the following properties:

$$\bigcap_{j \in \mathbb{Z}} V_j = 0 \qquad (1)$$

$$\mathrm{clos}_{L^2} \left( \bigcup_{j \in \mathbb{Z}} V_j \right) = L^2(\mathbb{R}) \qquad (2)$$

$$f(x) \in V_j \Leftrightarrow f(2x) \in V_{j+1} \qquad (3)$$

and there is a function $\varphi \in V_0$ such that $\{\varphi_{0,k}(x) = \varphi(x-k)\}_{k \in \mathbb{Z}}$ forms a Riesz basis for $V_0$ (4)

The term $\varphi$ is called the *mother scaling function* since, from (3), there exist $\{h_k\} \in L^2$ such that

$$\varphi(x) = \sqrt{2} \sum_{k \in \mathbb{Z}} h_k \varphi(2x - k)$$

This relation, called the *scaling relation*, will also hold for $\varphi(2x)$ and, by induction, for $\varphi(2^j x)$. In accordance with the notation in (4), we denote the translates and dilations of $\varphi$ by

$$\varphi_{j,k}(x) = 2^{j/2} \varphi(2^j x - k)$$

The set $\{\varphi_{j,k}\}$ forms a Riesz basis for $V_j$. We define $W_j$ to be the orthogonal complement of $V_j$ with respect to $V_{j+1}$. Just as $V_j$ is spanned by dilations and translations of the mother scaling function, so are the $W_j$'s spanned by translations and dilations of the *mother wavelet*. The mother wavelet is defined by

$$\psi(x) = \sqrt{2} \sum_k (-1)^{k-1} h_{-k+1+2M} \varphi_{1,k}(x)$$

with $M$ a particular integer. Daubechies² constructed compactly supported wavelets and scaling functions using a finite set of non-zero $\{h_k\}_{k=0}^{N-1}$ scaling parameters with $N=2M$ and $\sum_{k=0}^{N-1} h_k = \sqrt{2}$. $M$ is the order of the wavelet function. With these scaling parameters, the recursion formulas generate the desired orthogonal wavelets and scaling functions with $\mathrm{supp}(\varphi) = [0, N-1]$. Henceforth, these will be the wavelets we shall use, which we refer to as Daubechies wavelets of genus $N$.

Values for $\varphi$ are calculated using the scaling relation as indicated in the following procedure. First, the values it takes are determined at integer points. Then at the dyadic rationals at level 1 (the dyadic rationals at level $j$ are $\mathbb{D}_j = \{k/2^j\}_{k \in \mathbb{Z}}$); using that information, we calculate the values at $\mathbb{D}_{j+1}$ and so on, until $\varphi$ is defined at the dyadic rationals at all levels. Since the dyadics are dense in the reals, we simply extend $\varphi$ continuously to $\mathbb{R}$. This procedure creates a function that is continuous but not differentiable for $N=4$, differentiable but not twice differentiable for $N=6$, and
with increasing regularity for increasing $N$.² The nature of the scaling relation guarantees that $\varphi$ will be discontinuous in some derivative.¹¹ This procedure is easily accomplished computationally.

Another pleasing feature of these wavelets is their compact support. Whereas Fourier methods return global results, with compactly supported wavelets one can easily analyse short-lived events or pulses. Wavelet projection methods avoid distortion that might result from a local analysis with a windowed Fourier transform. As will be shown in the next section, compact support also makes the periodization of these wavelets an elegant process.

In addition to items (1)–(4) mentioned above, wavelets have a number of other interesting properties. These will be given without proof. For further details, see Daubechies² or Chui:¹

$$ \{\varphi_{j,k}\}_{j \ge 0, k \in \mathbb{Z}} \text{ is an orthonormal basis for } L^2(\mathbb{R}) \quad (5) $$

$$ V_{j+1} = V_j \oplus W_j \quad (6) $$

$$ L^2(\mathbb{R}) = \operatorname{ clos}_{L^2} \left( V_0 \bigoplus_{j=0}^{\infty} W_j \right) \quad (7) $$

$$ \{\varphi_{0,k}, \psi_{j,k}\}_{j \ge 0, k \in \mathbb{Z}} \text{ is an orthonormal basis for } L^2(\mathbb{R}) \quad (8) $$

$$ \int_{-\infty}^{\infty} \varphi(x) \, dx = 1 \quad (9) $$

$$ \sum_{k \in \mathbb{Z}} \varphi_{0,k} = 1 \quad (10) $$

$$ \int_{-\infty}^{\infty} \psi(x) x^k \, dx = 0: k = 0, \dots, M-1 \quad (11) $$

$$ \{x^k\}_{k=0}^{M-1} \in V_N \quad (12) $$

Item (6) is really at the heart of multiresolution analysis and provides wavelet-based analysis with a distinctly different resolving quality in contrast to spectral methods: to go to a higher resolution of spatial scale, one simply adds on the next wavelet level (the next $W_j$) as implied by item (6). At some given level (say, with a representation in $V_j$), the multiresolution property guarantees all spatial scale information at all coarser levels. In contrast, with Fourier methods, information about one frequency gives no information about other frequencies.

This study is devoted to the numerical calculation of connection coefficients involving periodized Daubechies wavelets. Connection coefficients are matrix structures that result from the evaluation of inner products

$$ \Omega_{j,k} \equiv \int \varphi_{j,k_0}^{(d_0)} \varphi_{j,k_1}^{(d_1)} \cdots \varphi_{j,k_n}^{(d_n)} dx $$

where $d_i$ is the number of differentiations with respect to $x$ of the scaling function $\varphi = \varphi(x)$. Inner products arise naturally in the context of the Galerkin solution of differential equations.⁹,¹² The name for these inner products was coined by Latto et al.,¹³ who developed the computational method presented herein. It replaces a quadrature problem by a linear algebra problem. Since the method is simple, fast, and general, it was considered worthy of implementing numerically. The technique developed by Latto et al. is by no means the only alternative available for the calculation of inner products involving wavelets. Of note is Dahmen and Micchelli's¹⁴ algorithm
which is similar to the approach favoured in this study in that the quadrature problem is recast as a linear algebra problem. The algorithm is more simply explained in Reference 4 (p. 396–401). This latter algorithm was numerically implemented very recently by Kunoth and is also publicly available software.¹⁵ The code is capable of performing the inner products of up to four wavelets and their derivatives.

Sections 2 and 3 of this study describe how periodized Daubechies wavelets on the unit interval are assembled and present estimates of their approximating characteristics. Furthermore, they will be shown to inherit most of the properties of their more standard Daubechies wavelets cousins, and capable of serving as basis set for $L^2$ functions on the unit interval which are periodic. Readers interested in the construction of a basis sets defined over the semi-infinite interval, or with different edge conditions over a finite interval, may find Cohen et al.¹⁶ a useful point of reference. Section 4 describes the connection coefficient algorithm. This section also covers several useful strategies for the periodization of general inner product matrices which arise in Galerkin solutions of differential equations, as well as some properties of the connection coefficients, which are useful in testing numerical codes and results, and in the manipulation of these coefficients.

In Section 5 the connection coefficient algorithm is evaluated by examining how wavelet-Galerkin methods approximate several commonly used differential operators. Finally, several appedices are provided with information on how to obtain the connection coefficients numerical code as well as further details on the properties and approximating characteristics of the periodized Daubechies wavelets.

## 2. PERIODIZED WAVELETS

Daubechies wavelets are defined on the whole real line. For many applications, however, wavelets defined on a periodic domain are needed. Interestingly, the wavelets defined above can be periodized with a Poisson summation technique to give periodic wavelets² that possess many of the same properties of their non-periodic kin. Moreover, the periodized $\varphi_{j,k}$ and $\psi_{j,k}$ are identical to their non-periodic forms except for wrapping around the edges of the domain; for large enough $j$, this too can be reduced to the non-periodic case for most calculations. Hence, many of the above-mentioned properties are preserved in the periodic case as a result of the construction by the “scaling” property of the non-periodic functions and their compact support.

The wavelets are periodized as follows:

$$ \tilde{\varphi}_{j,k}(x) \equiv \sum_{l \in \mathbb{Z}} \varphi_{j,k}(x-l) \quad \text{and} \quad \tilde{\psi}_{j,k}(x) = \sum_{l \in \mathbb{Z}} \psi_{j,k}(x-l) $$

By construction $\tilde{\varphi}$ and $\tilde{\psi}$ are periodic and are well defined on $[0, 1]$, since $\varphi$ and $\psi$ have compact support. Note that $\varphi_{j,k} = \varphi_{j,k'}$ if $k \equiv k' \mod(2^j)$. Thus we shall restrict our attention to $0 \le k < 2^j$. The same holds for the $\psi$'s. In what follows, the properties of the periodic wavelets will be investigated in some detail.

Periodized wavelet bases are not generated in quite the same way as the non-periodic versions. In the non-periodic case, bases are generated by repeated translation and dilation of the mother functions; but this approach is not possible in the periodic case because periodization does not commute with dilation. Therefore, the wavelet must be first dilated, then periodized. Although proof can be shown for the general case, Appendix II includes a proof that shows instead that the elements in $V_j$ for $j \le 0$ are all constant functions. If dilation commuted with periodization, this would not be true.
Periodized Daubechies wavelets form an orthonormal basis that still has a scaling relation, as shown in Appendix II. This means that in comparison with the non-periodic case, we have a chain of spaces $\tilde{V}_0 \subset \tilde{V}_1 \subset \cdots \subset L^2[0, 1]$ with the following properties:

$$ \bigcup_{j \ge 0} \tilde{V}_j = L^2[0, 1] \quad \text{with} \quad \tilde{V}_j = \operatorname{span}\{\tilde{\phi}_{j,k}\}_{k=0}^{2^j-1} \qquad (13) $$

$$ \bigcap_{j \in \mathbb{Z}} \tilde{V}_j = \{\text{constant functions}\} \qquad (14) $$

$$ f(x) \in \tilde{V}_j \Leftrightarrow f(2x) \in \tilde{V}_{j+1} \qquad (15) $$

$$ \begin{gathered} \text{By defining } \tilde{W}_j = \operatorname{span}\{\tilde{\psi}_{j,k}\}_{k=0}^{2^j-1} \text{ we see that } \tilde{W}_j \text{ is the orthogonal complement of } \tilde{V}_j \text{ in } \tilde{V}_{j+1}. \text{ So then } \operatorname{clos}\left(\tilde{V}_0 \bigoplus_{j=0}^{\infty} \tilde{W}_j\right) = L^2[0, 1] \quad (16) \\ \text{Clearly, differences exist between the properties of the periodic case and the non-periodic case. While they are both multiresolution spaces, the basis functions in the non-periodic case are all formed by translations and dilations of the mother scaling function, } \phi, \text{ while in the periodic case it is often impossible to derive } \tilde{\phi}_{j+1} \text{ from } \tilde{\phi}_j \text{ (for example, consider } \tilde{\phi}_1 \text{ and } \tilde{\phi}_0; \text{ the latter is a constant function and thus unable to represent the former). It turns out, however, that there is no relation between } \tilde{\phi}_{j+1} \text{ and } \tilde{\phi}_j \text{ for very small } j \text{ only. For } j \text{ suitably large, the periodic case actually looks exactly the same as the non-periodic case. Thus, for large enough } j, \text{ the periodization will affect the functions only by 'wrapping them around' the edges of the domain. This is also a strong argument for using the scaling functions as trial-and-test functions when using periodized wavelets. By choosing } j \text{ well, calculations can be performed as in the non-periodic case; if desired, a multiresolution can then be performed easily. Calculations will not be so simple if } \tilde{V}_0 \bigoplus_{k=1}^j \tilde{W}_k \text{ is used as a test space, because, for low } k, \text{ the basis for } \tilde{W}_k \text{ is not equivalent to the basis for } W_k.} \end{gathered} $$

Clearly, differences exist between the properties of the periodic case and the non-periodic case. While they are both multiresolution spaces, the basis functions in the non-periodic case are all formed by translations and dilations of the mother scaling function, $\phi$, while in the periodic case it is often impossible to derive $\tilde{\phi}_{j+1}$ from $\tilde{\phi}_j$ (for example, consider $\tilde{\phi}_1$ and $\tilde{\phi}_0$; the latter is a constant function and thus unable to represent the former). It turns out, however, that there is no relation between $\tilde{\phi}_{j+1}$ and $\tilde{\phi}_j$ for very small $j$ only. For $j$ suitably large, the periodic case actually looks exactly the same as the non-periodic case. Thus, for large enough $j$, the periodization will affect the functions only by 'wrapping them around' the edges of the domain. This is also a strong argument for using the scaling functions as trial-and-test functions when using periodized wavelets. By choosing $j$ well, calculations can be performed as in the non-periodic case; if desired, a multiresolution can then be performed easily. Calculations will not be so simple if $\tilde{V}_0 \bigoplus_{k=1}^j \tilde{W}_k$ is used as a test space, because, for low $k$, the basis for $\tilde{W}_k$ is not equivalent to the basis for $W_k$.

### 3. APPROXIMATION RESULTS

A function $f \in L^2[0, 1]$ may be projected into the wavelet basis and expressed as $f = \sum_k a_{kj} \tilde{\phi}_{j,k}$, where $a_{kj} = \langle f, \tilde{\phi}_{j,k} \rangle$. The calculation of $a_{kj}$ is usually performed numerically. In Reference 7, a method was developed which employs Taylor series expansions to approximate $f$. The method requires use of the moment equations to make $\mathcal{O}(h^n)$ approximations for $f \in C^n$. Unfortunately, in applications such as signal processing or any area where only a finite number of samples of $f$ are provided, this method may produce unacceptable results. As an alternative, a function $f$ given as samples on $\mathbb{D}_j \cap [0, 1]$ may be approximated by a function $\tilde{f} \in \tilde{V}_j$.

Define the samples of $f$ as $\tilde{f} \in \mathbb{R}$, with the $k$th component of $\tilde{f} = f((k-1)/2^j)$. Construct $\tilde{\phi}_{j,k}$ from $\tilde{\phi}_{j,k}$ as shown previously. This yields, for $j > \log_2(N-2)$, a linearly independent spanning set for $\mathbb{R}^{2^j}$. Further, since $\text{supp}(\tilde{\phi}_{j,0}) = [0, 2^{-j}(N-1)]$, $\tilde{\phi}_{j,k}$ takes only $N-2$ values on $\mathbb{D}_j$, and $\tilde{\phi}_{j,k}$ is just the $k$th forward cyclic permutation of the elements of $\tilde{\phi}_{j,0}$. The problem is thus reduced to finding the unique representation of $\tilde{f}$ in terms of $\{\tilde{\phi}_{j,k}\}_{0 \le k < 2^j}$, which is simply the solution of

$$ A \tilde{v} = \tilde{f} $$
where $A$ is the very sparse transformation matrix from the standard basis to the $\tilde{\phi}$ basis. Define
$f(x) = \sum_k v_k (k) \tilde{\phi}_{j,k}(x)$. By construction, $\tilde{f}|_{\mathbb{D}_j} = f$.

In summary, this method involves solving the inverse problem with a sparse matrix, and results
in a function $\tilde{f}$ with $\tilde{f}|_{\mathbb{D}_j} = f$. In essence, $\tilde{f}$ will be equal to $f$ at all the sampled values.

We next show the manner in which the periodized wavelets may be used in the context of
functional approximation. Since the wavelets form an orthonormal basis, the orthogonal projection
operators onto $V_j$ and $W_j$ are defined respectively as

$$
P_j(f) = \sum_{k=0}^{2^j - 1} \langle f, \tilde{\phi}_{j,k} \rangle \tilde{\phi}_{j,k}, \quad Q_j(f) = \sum_{k=0}^{2^j - 1} \langle f, \tilde{\psi}_{j,k} \rangle \tilde{\psi}_{j,k}
$$

As we have already seen, periodized wavelets provide a basis for $L^2[0,1]$ so we have $\|f - P_j f\|_2 \to 0$, as $j \to \infty$. This is a property of any orthonormal basis of $L^2$, but this particular periodized basis has some additional properties. If $f$ is a continuous function on the torus, then $\|f - P_j f\|_\infty \to 0$ as $j \to \infty$. If $f \in L^1[0,1]$, then $\|f - P_j f\|_1 \to 0$ as $j \to \infty$. The proofs for these two results appear in Appendix III of this study. These two last statements are strong, in contrast to the convergence properties of Fourier functional approximations. In fact it has been shown that the continuous functions whose Fourier series do not uniformly converge are dense in $C(\mathbb{T})$.¹⁷ In this sense, wavelets provide a much more general basis than Fourier bases and hence have potentially broader applications. These results suggest that wavelets should do a better job at pointwise approximation, especially for continuous functions.

To measure the error to which a truncated projection will approximate a desired function, we
shall estimate its convergence. The natural choice of norms with which to measure convergence
is the Sobolev norms. The sth Sobolev norm of a function $f$ is defined as

$$
\|f\|_{H^s} = \left( \int (1 + \xi^2)^s \hat{f}^2(\xi) \, d\xi \right)^{1/2}
$$

where $H^s$ consists of those functions whose $s$ Sobolev norm exists and is finite. Daubechies$^2$ states
that the norm for $H^s$ is equivalent to

$$
\|f\|_{H^s[0,1]} = \left( \sum_{\substack{j \ge 0 \\ 0 \le k < 2^j}} (1 + 2^{2js}) \langle f, \hat{\psi}_{j,k} \rangle^2 \right)^{1/2}
$$

Using this result, we easily find a bound for $\|f - P_p f\|_2$. For $f \in H^s[0,1]$,

$$
\begin{align*}
\|f - P_p f\|_2 &= \left\| \sum_{j \ge p0 \le k < 2^j} \langle f, \hat{\psi}_{j,k} \rangle \right\|_2^{1/2} \\
&= \left[ \sum_{j,k} \langle f, \hat{\psi}_{j,k} \rangle^2 \right]^{1/2} = \left[ \sum_{j,k} \frac{2^{2js}}{2^{2js}} \langle f, \hat{\psi}_{j,k} \rangle^2 \right]^{1/2} \\
&\le \left[ \sum_{j,k} \frac{2^{2js}}{2^{2ps}} \langle f, \hat{\psi}_{j,k} \rangle^2 \right]^{1/2} \\
&= 2^{-ps} \left[ \sum_{j,k} 2^{2js} \langle f, \hat{\psi}_{j,k} \rangle^2 \right]^{1/2} \le 2^{-ps} \|f\|_{H^s[0,1]}^{1/2}
\end{align*}
$$
The same technique may be used to find error bounds in the $H^l$ norm for $l \le p$:

$$ \begin{aligned} \|f - P_p f\|_{H^l[0,1]} &= \left( \sum_{\substack{j \ge p \\ 0 \le k < 2^j}} (1 + 2^{2jl}) \langle f, \hat{\psi}_{j,k} \rangle^2 \right)^{1/2} \\ &= \left( \sum_{j,k} (1 + 2^{2jl}) \frac{(1 + 2^{2j(s-l)})}{(1 + 2^{2j(s-l)})} \langle f, \hat{\psi}_{j,k} \rangle^2 \right)^{1/2} \\ &\le (1 + 2^{2p(s-l)})^{-1/2} \left( \sum_{j,k} (1 + 2^{2sj}) \langle f, \hat{\psi}_{j,k} \rangle^2 \right)^{1/2} \\ &\le 2^{-p(s-l)} \|f\|_{H^s[0,1]} \end{aligned} $$

These bounds are similar to those on Fourier series.¹⁸

### 4. METHOD FOR COMPUTING CONNECTION COEFFICIENTS

In wavelet applications, one often must represent operators in terms of wavelets.¹⁹ An example of such an application is the Galerkin solution of differential equations. The formulation of solutions will require integrations of the form

$$ \Omega_{k_1, k_2, \dots, k_n}^{d_1, d_2, \dots, d_n} = \langle \tilde{\phi}_{j, k_0}^{(d_0)} \tilde{\phi}_{j, k_1}^{(d_1)} \dots \tilde{\phi}_{j, k_n}^{(d_n)} \rangle = \int_0^1 \tilde{\phi}_{j, k_0}^{(d_0)} \tilde{\phi}_{j, k_1}^{(d_1)} \dots \tilde{\phi}_{j, k_n}^{(d_n)} dx $$

where $\tilde{\phi}^{(d)} = d^d \tilde{\phi}/dx^d$. This expression is an n-term connection coefficient. Since $\tilde{\phi}$ cannot be represented in closed form for $N > 2$ and, by construction, has limited regularity, analytic calculation of the integral is impossible, and numerical quadrature is often inaccurate as a result of the wildly oscillating nature of the resulting kernels. An alternative approach developed by Latto et al.¹³ circumvents some of the difficulty by exploiting the scaling relation and the moment condition to reduce the calculation to an eigenvector problem. Their method is designed for non-periodic compactly supported wavelets. However, by invoking an extension of the earlier result regarding the equivalence of periodized and non-periodized wavelets in Section 2, one may infer that for $j \ge \log_2((N-1)n)$, the periodized case yields the same result as the non-periodized case. For illustration and for completeness in what follows we adopt closely the general procedure given in detail in Reference 13 to the 2-tuple case. Several tabulated connection coefficients using periodized wavelets are included in an appendix of this study.

First, integration by parts is performed repeatedly on the above integral to obtain

$$ \Omega_{k_1, k_2}^{d_1, d_2} = (-1)^{d_1} \Omega_{k_1, k_2}^{0, d_2+d_1} $$

where the periodicity of the wavelets has been invoked. By changing variables, we further reduce the equation to

$$ \Omega_{k_1, k_2}^{0, d} = \Omega_{0, k_2-k_1}^{0, d} \equiv \Lambda_{k_2-k_1}^{d} \quad \text{where } d = d_1 + d_2 $$

From these relations it is clear that any 2-tuple can be represented by a $\Lambda_k^d$.
To construct the eigenvector problem, fix $d$, then solve for $\{\Lambda_k^d\}_{0 \le k < 2^j}$, creating a system of $2^j$ homogeneous relations in $\Lambda_k^d$ and enough inhomogeneous equations to reduce the dimension of the associated eigenspace to 1. Although we are using the connection-coefficient method for the non-periodized case, we are computing them for the periodic case (by equivalence), which is where the bounds on $k$ come into play.

## 4.1. Connection coefficients algorithm

To generate the homogeneous equations fix $d$, $j \in \mathbb{N}$, such that $\tilde{\phi}_{j,k}^{(d)}$ is well defined. To simplify notation, denote $\tilde{\phi}_{j,k}^{(d)} \equiv \Phi_k^d$. In Reference 13 it is suggested, without proof, that this method also holds for the first $d$ for which $\Phi^d$ is discontinuous. For low-order differential equations, however, $N=6$ or $N=8$ wavelets should provide sufficient regularity. Since for every $0 \le k < 2^j$,

$$ 
\begin{aligned}
\Lambda_k^d &= \int \Phi_0(x) \Phi_k^d(x) \, dx = \int \left( \sum_{m=0}^{N-1} h_m \Phi_m(2x) \right) \left( \sum_{l=0}^{N-1} h_l \Phi_{l+2k}^d(2x) \right) 2^d d(2x) \\
&= 2^d \sum_m \sum_l h_m h_l \int \Phi_m(2x) \Phi_{l+2k}^d(2x) \, d(2x) \\
&= 2^d \sum_m \sum_l h_m h_l \int \Phi_0(\zeta) \Phi_{l+2k-m}^d(\zeta) \, d\zeta, \\
\Lambda_k^d &= 2^d \sum_{m=0}^{N-1} \sum_{l=0}^{N-1} h_m h_l \Lambda_{l+2k-m}^d
\end{aligned}
$$

In the above discussion, the integration is over the real line.

This linear homogeneous system can be represented as

$$ A \lambda \underset{\rightarrow}{=} 2^{-d} \lambda \underset{\rightarrow}{=} $$

where $\lambda \underset{\rightarrow}{=} \{\Lambda_k^d\}_{0 \le k < 2^j}$. It is worth noting here that if one needs to compute an $n$-tuple connection coefficient for $j < \log_2((N-1)n)$, then the periodic scaling relation can be used to resolve each $\phi_{j,k}$ into the sum of $\phi_{j',k}$'s with $j' \ge \log_2((N-1)n)$. Thus, the reduction to the non-periodic case is a universally applicable method.

To generate the inhomogeneous equations, we must first assume $d \le M-1$. The moment condition then guarantees that

$$ x^d = \sum_{l \in \mathbb{Z}} \tilde{M}_l^d $$

where $\tilde{M}_l^d = \langle x^d, \phi_{0,l} \rangle$. Setting $x = 2^{j\zeta}$ and defining $M_l^d = \langle x^d, \Phi_l \rangle$, we have

$$ \tilde{M}_l^d = 2^{dj} 2^{j/2} \langle \zeta^d, \Phi_l \rangle = 2^{dj} 2^{j/2} M_l^k $$

This gives the relation

$$ \zeta^d = \sum_{l \in \mathbb{Z}} M_l^d \Phi_l(\zeta) $$

which, when differentiated $d$ times, yields

$$ d! = \sum_{l \in \mathbb{Z}} M_l^d \Phi_l^d(\zeta) $$
Multiplying by $\Phi_0^0$ and integrating, we obtain

$$ \sum_{l \in \mathbb{Z}} M_l^d \int \Phi_0^0(\zeta) \Phi_l^d(\zeta) \, d\zeta = d! \int \varphi_{j,0}(\zeta) \, d\zeta = d! 2^{-j/2} $$

Thus $\sum_l M_l^d \Lambda_l^d = d! 2^{-j/2}$. The sum over $l$ is actually over $|l| \le N-2$ since the $\varphi$'s are compactly supported. Thus, by changing the indices of summation by $m = l+1+(N-2)$, the inhomogeneous equations are

$$ \sum_{m=1}^{2N-3} \Lambda_m^d M_{m-1-(N-2)}^d $$

with

$$ M_l^d = 2^{-j(2d+1)/2} \tilde{M}_l^d $$

The linear system formed by the $2^j$ homogeneous equations and the above inhomogeneous equations has eigenspace dimension equal to 1. Thus, all that remains to complete the system is to calculate $\tilde{M}_l^d$:

$$ \begin{aligned} \tilde{M}_l^k &= \int x^d \varphi(x-l) \, dx = \int (y+n)^k \varphi(y) \, dy \\ &= \int \sum_{j=0}^k \binom{k}{j} y^j n^{k-j} \varphi(y) \, dy \\ &= \sum_{j=0}^k \binom{k}{j} n^{k-j} \tilde{M}_0^j \end{aligned} $$

Since $\tilde{M}_l^0 = 1$ is a property of the $\varphi$, the above relation can be used to evaluate recursively $\tilde{M}_l^d$ for all $l$.

The linear system is now complete and fixes the values of $\Lambda_k^d$. Note that while the scaling relation, which is used to generate the homogeneous relations, exists for periodized wavelets, currently nothing is analogous to the moment condition that may be used to generate the necessary inhomogeneous equations. One possible approach is to use the scaling equation linked with the fact that $\tilde{V}_j = \{\text{constant functions}\}$ for $j \le 0$ to find additional inhomogeneous relations. Problems arise with relating $\int \tilde{\phi}_{j+1,0} \, dx$ to $\int \tilde{\phi}_{j,0} \, dx$, however, since dilation does not commute with periodization. While this method could probably be worked out, the periodic case can always be reduced to an equivalent non-periodic case for which the method is already well defined. Thus, to compute an $n$-term connection coefficient for periodized wavelets, one need only resolve the terms into $\tilde{V}_j$, for some $j \ge \log_2(N-1)n$ and apply the above method. This approach takes full advantage of the equivalence of periodic and non-periodic scaling functions and circumvents the need for a connection coefficient method particular to periodized wavelets.

## 4.2. Periodization of inner product matrices

The two most common situations that arise in the wavelet Galerkin solution of differential equations involve 2-tuples and 3-tuples. Periodic conditions require wraparound of the $\Lambda$ entries in $\Omega$ in the upper right-hand corner and the lower left-hand corner, assuming that the matrix row index increases from top to bottom corresponding to the dyadic points in $[0, 1)$. The wraparound
<figure><img src="image_2.png" /><figcaption>Figure 1. Structure of $\Lambda_{j'',k''}$ index pad for 3-tuples in the $(j,l)$ plane: (a) the index pad corresponds to $N=4$ in the region where no wraparound arises; (b) corresponding to $N=4$ for $k=k'+2$; (c) corresponding to $N=4$ for $k=N-3$, $\sigma=2$</figcaption></figure>

will produce at most $N(N+1)$ additional entries in the matrix, where $[0,N]$ is the support of the scaling functions. In this section we present explicit schemes that enable the construction of these matrices. The 2-tuple case is easily generated, since $\Omega$ is a circulant matrix of the form

$$ \Omega = \text{circ}(\Lambda_0, \Lambda_1, \dots, \Lambda_N, \dots, 0, 0, \dots, \Lambda_{-N}, \Lambda_{-N+1}, \dots, \Lambda_{-1}) $$

and $\Omega = \pm \Omega^T$, depending on whether the operator is symmetric or skew-symmetric (see Section 4.3).

The 3-tuple is somewhat more involved. The $\Lambda_{j'',k''}$ have the structure illustrated schematically in Figure 1(a). We call this structure the ‘index pad’. Let $K=2^p-1$, where the superscript $p$ corresponds to the resolution of the scaling functions. A typical situation in a Galerkin discretization when 3-tuples are involved would be the calculation of

$$ s_k = \sum_{l=0}^{K} a_l \omega_{k,l} $$
where

$$ \omega_{k,l} = \sum_{j=0}^{K} h_j \Omega_{k,j,l} $$

where $a_l$ and $h_j$ are the wavelet coefficients of the projection of real quantities $A(x, \cdot)$ and $H(x, \cdot)$ into $V^p$, and $\Omega$ is the inner product $\langle \tilde{\phi}_k^{d_1} \phi_j^{d_2} \phi_l^{d_3} \rangle$. For a given $k$ there is an index pad in $(j, l)$ over which these indices range when forming $\omega_{k,l}$ and $s_k$. Moreover, associated with each index $(j, l)$ in the pad there is a $\Lambda_{j'',k''}$, where $j'' = j-k$ and $l'' = l-k$, which gives the value of $\Omega_{k,j,l}$. It is easier to think of the index pad as having a rectangular structure in which the entries lying outside of the hexagon of the true index pad are zero.

First consider the case in which there is no wraparound, that is, for $N \le k < k'$, where $k' = K - (N - 1)$, the largest index $k$ for which the support of the scaling function lies entirely in $[0, 1]$. Associated with each $(j'', l'')$ in the index pad, there is a connection coefficient $\Lambda_{j'',l''}$. As $k$ is varied, the pad centred at $k$, moves along the line $l=j$ in the $(j, l)$ plane. This is the 'regular case' and is shown in Figure 1(a).

Next, consider the case $k \ge k'$. Let $k = k' + \sigma$, $0 \le \sigma \le N - 1$. The situation is illustrated in Figure 1(b). The index pad is beyond $(j, l) = (K, K)$. For convenience identify the following subregions.

$$ A: k - N \le j \le K $$

$$ k - N \le l \le K $$

$$ B: j = K + 1 + q, \quad 0 \le q \le \sigma $$

$$ l = K - t, \quad 0 \le t \le N - 1 $$

$$ C: j = K - t, \quad 0 \le t \le N - 1 $$

$$ l = K + 1 + q, \quad 0 \le q \le \sigma $$

$$ D: j = K + 1 + q, \quad 0 \le q \le \sigma $$

$$ l = K + 1 + q, \quad 0 \le q \le \sigma $$

Consider a given $k$, the corresponding index pad is subdivided into the four subsets $A, B, C, D$ defined above. When index pairs $(j, l)$ lie in subregions $B, C$ or $D$, periodization comes into play since the corresponding basis functions have been periodically extended. We now consider the consequences of this periodization on the $s_k$ in detail. Note that $s_k$ will be the sum of three partial sums, summed over subsets of the index pad. Next we observe that the range of $l$ is $k - N \le l \le K + 1 + \sigma$. We divide this range into two parts: $k - N \le l \le K$ and $K + 1 \le l \le K + 1 + \sigma$. For the range $k - N \le l \le K$, the index $l$ is not affected by the periodization since $l \le K$. The index $j$, on the other hand, produces a pair $(j, l)$ that ranges over the set $A$ and $B$. Let $\omega_{k,l} = \omega_{k,l}^r + \omega_{k,l}^p$. Then

$$ \omega_{k,l}^r = \sum_{j=k-N}^{K} h_j \Lambda_{j'',l''} \quad (\text{the regular part}) $$

$$ \omega_{k,l}^p = \sum_{j=K+1}^{K+1+\sigma} h_{j-K-1} \Lambda_{j'',l''} \quad (\text{the periodic adjustment}) $$
and thus

$$s_k^{(1)} = \sum_{l=0}^{K} a_l \omega_{k,l}$$

where, throughout this section, we set $j'' = j - k$ and $l'' = l - k$.

Next we deal with the range $K + 1 \le l \le K + 1 + \sigma$. In this case the pair $(j, l) \in C \cup D$. When $(j, l) \in C$,

$$\omega_{k,l} = \sum_{j=k-N}^{K} h_j \Lambda_{j'',l''}$$

Now the index $l$ is affected by periodization, so that

$$s_k^{(2)} = \sum_{l=K+1}^{K+1+\sigma} a_{l-K-1} \omega_{k,l}$$

When $(j, l) \in D$, both indices are affected by periodization, thus

$$\omega_{k,l} = \sum_{j=K+1}^{K+1+\sigma} h_{j-K-1} \Lambda_{j'',l''}$$

Then,

$$s_k^{(3)} = \sum_{l=K+1}^{K+1+\sigma} a_{l-K-1} \omega_{k,l}$$

So, for $k = k' + \sigma$, $0 \le \sigma \le N - 1$, combining, we have

$$s_k = s_k^{(1)} + s_k^{(2)} + s_k^{(3)}$$

We now consider the case when the index pad encounters the left boundary, $j=0$, $l=0$, which occurs for $0 \le k \le N - 1$. Schematically, the situation is portrayed in Figure 1(c). Denote the following subregions:

$$A: 0 \le j, l \le k + N$$

$$B: 0 \le l \le N - 1$$

$$-\sigma - 1 \le j \le -1$$

$$C: 0 \le j \le N - 1$$

$$-\sigma - 1 \le l \le -1$$

$$D: -\sigma - 1 \le j \le -1$$

$$-\sigma - 1 \le l \le -1$$

where $0 \le \sigma \le N - 1$.

Again we divide the entire $l$ range $-\sigma - 1 \le l \le k + N$ into two parts: $0 \le l \le k + N$ and $-\sigma - 1 \le l \le -1$. Consider the range $0 \le l \le k + N$ for which the index pair $(j, l)$ ranges over
B and A. When $(j,l) \in B$, the index $j$ is affected by periodization so that $\omega_{k,l}$ has a regular part and a periodic part. Let $\omega_{k,l} = \omega_{k,l}^r + \omega_{k,l}^p$. Then

$$ \omega_{k,l}^r = \sum_{j=0}^{k+N} h_j \Lambda_{j'',l''} \quad (\text{the regular part}) $$

$$ \omega_{k,l}^p = \sum_{j=-\sigma-1}^{-1} h_{K+1+j} \Lambda_{j'',l''} $$

thus

$$ s_k^{(1)} = \sum_{l=0}^{k+N} a_l \omega_{k,l} $$

Next we deal with the range $-\sigma - 1 \le l \le -1$. In this case $(j,l) \in C \cup D$. For $(j,l) \in C$, we have

$$ \omega_{k,l} = \sum_{j=0}^{N-1} h_j \Lambda_{j'',l''} $$

Then

$$ s_k^{(2)} = \sum_{l=-\sigma-1}^{-1} a_{K+1+l} \omega_{k,l} $$

When $(j,l) \in D$ we have

$$ \omega_{k,l} = \sum_{j=-\sigma-1}^{-1} h_{K+1+j} \Lambda_{j'',l''} $$

and

$$ s_k^3 = \sum_{l=-\sigma-1}^{-1} a_{K+1+l} \omega_{k,l} $$

Combining, for $k = N - \sigma - 1$, $0 \le \sigma \le N - 1$, we have

$$ s_k = s_k^{(1)} + s_k^{(2)} + s_k^{(3)} $$

The periodization is then complete. The result is the vector $\{s_k\}$. This procedure generalizes to the $n$-tuple case in a straightforward manner.

### 4.3. Useful connection coefficient relations

We list a number of identities and relations that are useful in the manipulation of connection coefficients. Many of these relations appear in References 13 and 20. Except for the operator inversion relationship, most of the inner product relations can easily be derived by invoking integration by parts, translation, and change of variables.

$$ \Omega_{l,m}^{d_1,d_2} = \Omega_{m,l}^{d_2,d_1} \qquad (17) $$

$$ \Lambda_l^{d_1, d_2} = -\Lambda_l^{d_1-1, d_2+1} \qquad (18) $$

$$ \Lambda^{1, d_2, d_3} = -\Lambda^{0, d_2+1, d_3} - \Lambda^{0, d_2, d_3+1} \qquad (19) $$
$$ \Lambda_{l,m}^{d_1,d_2,d_3} = \Lambda_{-l,-m-i}^{d_2,d_1,d_3} \qquad (20) $$

$$ \Lambda_{l,m}^{d_1,d_2,d_3} = \Lambda_{l-m,-m}^{d_3,d_2,d_1} \qquad (21) $$

$$ \Lambda_{l,m}^{d_1,d_2,d_3} = (-1)^{d_1} \sum_{i=0}^{d_1} \begin{pmatrix} d_1 \\ i \end{pmatrix} \Lambda_{l,m}^{0,d_2+i, d_3+d_2+d_3+i} \qquad (22) $$

Another relation that is very useful in checking the construction and accuracy of 3-tuple matrices is the check-sum procedure: the column sum of a 3-tuple matrix must equal a corresponding 2-tuple vector component for component. For example, $\sum_{l=0}^{K} \Omega_{l,m}^{1,0,0} = \Omega_{m}^{1,0}$.

Lastly, we mention an efficient procedure for the inversion of a first-order operator, which exploits the symmetric or skew-symmetric nature of the matrix, $\Omega = \pm \Omega^T$, so that the matrix may be diagonalized, $\Omega = \Phi D_{\Omega} \Phi^T$. A concrete example of operator inversion appears in Reference 20.

## 5. APPROXIMATION OF DIFFERENTIAL OPERATORS

The spectrum $\{\sigma_j\}$ of the continuous differential operator $\mathcal{L}$ with periodic boundary conditions on 0 and 1 is discrete. A discrete approximation $L$ of the operator may be found by projecting the operator onto the space spanned by the periodized scaling functions. The discrete and continuous operators can be compared by looking at their spectra. By the Galerkin procedure outlined in Section 4.1 one can form the matrix $\Omega$, an approximate representation of $L$, in the subspace of periodized scaling functions of resolution $p$ and of genus $N$. In this section we examine the usefulness and properties of wavelet projection techniques for differential operators. A comparison of the spectrum of the differential operator $L$ and its approximation using the connection coefficients algorithm reveals information important within the context of solving ordinary and partial differential equations using wavelet-Galerkin techniques.

The projection of the differential operator with periodic boundary conditions leads to a 2-tuple $\Omega$. Eigenvalues of this matrix may be found by noting²¹ that the circulant matrix $\Omega$, has eigenvalues $\lambda_j$ given by the expression

$$ \lambda_j = \sum_{k=0}^{K} a_k e^{i2\pi kj/K} $$

where $a_k$ is an entry in the first row in the circulant matrix, namely, an element of the set $\Lambda_0, \Lambda_1, \dots, \Lambda_N, 0, \dots, 0, \Lambda_{-N}, \Lambda_{-N+1}, \dots, \Lambda_{-1}$ and $j = 1, 2, \dots, K+1$.

The $\Lambda$'s were obtained by the procedure outlined in the preceding section. The overdetermined system was solved using a QR algorithm from LAPACK²² called dgess. For $\Omega^{0,1}$, the residual and the condition number as a function of the resolution $p$ and the genus $N$ appears in Table I.

Some of the eigenvalues of $\Omega^{0,1}$ are given in Tables II and III. Regarding these tables, one can make three general observations. First, the imaginary part of the eigenvalues grows only as $O(2^p)$, just as the Fourier-basis computation would, and this is in contrast to spectral approximations based on Chebychev or Legendre bases. Second, we observe that the real parts of the eigenvalues are comparable to the size of machine precision zero, since the real part of each eigenvalue is proportional to the $\Lambda_0$, again in sharp contrast to Chebychev and Legendre bases investigated by Trefethen et al.²³ Third, little evidence of contamination due to round-off error was seen even for high resolution (in our experiments, we tried values of $p$ as high as 11), in sharp contrast to the experience of Trefethen et al. with the above-mentioned bases functions. The spectrum of
<table><caption>Table I. Residual and condition number for the least-squares calculation of $\Lambda$'s used in $\Omega^{0,1}$, as a function of $p$ and $N$</caption><thead><tr><th>N</th><th>p</th><th>Residual</th><th>Condition number</th></tr></thead><tbody><tr><td>4</td><td>4</td><td>-2.1926904736347D - 15</td><td>137</td></tr><tr><td>4</td><td>5</td><td>-8.4099394115356D - 15</td><td>387</td></tr><tr><td>4</td><td>6</td><td>-8.3266726846887D - 15</td><td>1096</td></tr><tr><td>4</td><td>7</td><td>-2.6922908347160D - 14</td><td>3100</td></tr><tr><td>4</td><td>8</td><td>-2.1510571102112D - 16</td><td>8770</td></tr><tr><td>6</td><td>4</td><td>4.3905851176973D - 15</td><td>189</td></tr><tr><td>6</td><td>5</td><td>4.3905851176973D - 15</td><td>531</td></tr><tr><td>6</td><td>6</td><td>2.2332830029725D - 14</td><td>1503</td></tr><tr><td>6</td><td>7</td><td>3.2980562725271D - 14</td><td>4252</td></tr><tr><td>6</td><td>8</td><td>-1.6563139748627D - 14</td><td>12 029</td></tr><tr><td>8</td><td>4</td><td>-1.6766260279716D - 11</td><td>229</td></tr><tr><td>8</td><td>5</td><td>-3.3533817109376D - 11</td><td>642</td></tr><tr><td>8</td><td>6</td><td>6.7080697513231D - 11</td><td>1816</td></tr><tr><td>8</td><td>7</td><td>1.3417001243474D - 10</td><td>5138</td></tr><tr><td>8</td><td>8</td><td>2.6830899981444D - 10</td><td>14 532</td></tr></tbody></table>

<table><caption>Table II. Imaginary part of eigenvalues divided by $2\pi$ for increasing resolution $2^p$ calculated using $N=6$ periodized Daubechies wavelets</caption><thead><tr><th>p = 4<br>0.0D0</th><th>p = 6<br>0.0D0</th><th>p = 8<br>0.0D0</th></tr></thead><tbody><tr><td>0.99997222310553</td><td>0.99999999300555</td><td>0.99999999999689</td></tr><tr><td>1.9967784015701</td><td>1.9999991102504</td><td>1.999999997795</td></tr><tr><td>2.0423307943545</td><td>2.9999849535965</td><td>2.9999999962590</td></tr><tr><td>2.9532542889943</td><td>3.9998888924219</td><td>3.9999999720194</td></tr><tr><td>3.4758292699748</td><td>4.9994799193486</td><td>4.9999998667457</td></tr><tr><td>3.7208826878105</td><td>5.9981781285634</td><td>5.9999995232107</td></tr><tr><td>4.0025895152829</td><td>6.9947814463569</td><td>6.9999985996877</td></tr></tbody></table>

the skew-symmetric operator $d/dx$ is purely imaginary and equal to $2\pi k$, where $k \in \mathbb{Z}$. Periodicity will lead to the eigenvalue 0 having multiplicity 2. Table II shows the magnitude of the imaginary part of the first few eigenvalues, divided by $2\pi$, as a function of $p$. It is clear from the table that approximations to the eigenvalue improve as $p$ is increased and the number of reasonably correct eigenvalues grows as well with $p$. Parenthetically, we remark that with $p=6$ we can go as far $k=11$ for which $\lambda_{11} \doteq 10.8$ while retaining at least one digit of accuracy. In the case $p=8$ we can go as far as $k=36$ for which $\lambda_{36} \doteq 35.8$.

We examined the dependence of the eigenvalue convergence on both the resolution $p$ and the genus $N$. We did this by examining the imaginary part of several eigenvalues. We found that the rate of convergence was almost exactly $2^{-p}$, while the dependence of the rate of convergence on $N$ was essentially quadratic. Tables II and III illustrate the above-mentioned rates of convergence.
<table><caption>Table III. Imaginary part of eigenvalues of d/dx as a function of N for p = 6</caption><thead><tr><th>N = 4</th><th>N = 6</th><th>N = 8</th></tr></thead><tbody><tr><td>0.99999690699928</td><td>0.99999999300555</td><td>0.99999999998407</td></tr><tr><td>1.9999013641150</td><td>1.9999991102504</td><td>1.9999999916456</td></tr><tr><td>2.9992552692402</td><td>2.9999849535965</td><td>2.9999996834005</td></tr></tbody></table>

<table><caption>Table IV. Spectrum approximation of d³/dx³ with N = 6. Imaginary part of eigenvalues</caption><thead><tr><th>p = 4</th><th>p = 6</th><th>p = 8</th></tr><tr><th>0.0D0</th><th>0.0D0</th><th>0.0D0</th></tr></thead><tbody><tr><td>39.411980740054</td><td>39.478149261094</td><td>39.478416554374</td></tr><tr><td>308.17967443293</td><td>315.7932221828</td><td>315.82720643905</td></tr><tr><td>956.10277711077</td><td>1065.3407907411</td><td>1065.9149806290</td></tr><tr><td>1670.0988161963</td><td>2522.3667673633</td><td>2526.6015527133</td></tr><tr><td>1877.4681030848</td><td>4914.9287548852</td><td>4934.7204114901</td></tr><tr><td>2625.9035313813</td><td>8457.8432634425</td><td>8527.0455871351</td></tr><tr><td>2641.9768303508</td><td>13 342.439258378</td><td>13 540.237953821</td></tr></tbody></table>

<table><caption>Table V. Real part of eigenvalues for d²/dx² divided by 4π² for increasing resolution p. N = 6</caption><thead><tr><th>p = 4</th><th>p = 6</th><th>p = 8</th></tr></thead><tbody><tr><td>1.1518884962509D-14</td><td>2.7645323910020D-13</td><td>0.0D0</td></tr><tr><td>1.0033162940157</td><td>1.0000135268157</td><td>1.0000000529814</td></tr><tr><td>4.1847329915128</td><td>4.0008582644240</td><td>4.0000033890207</td></tr><tr><td>10.666776606377</td><td>9.0096362703930</td><td>9.0000385683104</td></tr><tr><td>22.726823972803</td><td>16.053060704248</td><td>16.000216429062</td></tr><tr><td>41.582358732169</td><td>25.197220008678</td><td>25.000824274322</td></tr></tbody></table>

Qualitatively, we find the same type of behaviour in higher odd-ordered differential operators. Table IV shows the approximation to the spectrum of the operator d³/dx³ with periodic boundary conditions. We found the rates of convergence very similar to those of the first derivative, with the proviso that the basis functions had to be chosen with sufficient smoothness.

Finally we consider approximation to even-ordered differential operators. For example, Table V shows the approximation of the spectrum of d²/dx² for periodic boundary conditions for the same values of p. This case corresponds to the eigenvalues of Ω⁰,². Each eigenvalue has a complex conjugate, for which the imaginary part is no greater than machine precision zero. Again, we found the same rates of convergence as in the odd-ordered cases.

The condition number (cn) of the vector calculation was found to have the following correspondence: for p = 5 it was cn = 52 778, p = 6 it was cn = 298 557, p = 7 it was cn = 1 688 892. In general it was found that the condition number grew as the size of p, N, and on the order of the derivatives was increased. The source of the high condition number is the ill-posed nature of the over-determined matrix. The high condition number is especially troublesome in the
calculation of the 3-tuples. There are few options available to improve the ill-conditioned nature of the matrix. One option could have been to use Chebychev polynomials in connection with the moment conditions. However, the ill-posedness is in the linear homogeneous equations, its exact nature being that several of the rows will have nearly the same entries for some values of $p, N$, and the order of the derivatives. Hence, little or no improvement would result in the use of the Chebychev polynomials as a way to mix the moment conditions.

Using the eigenvalue data corresponding to the three differential operators, we have also been able to assess numerically the size of the largest eigenvalues. This is a useful estimate, for example, in the determination of time stability in numerical schemes for the solution of differential equations. The estimate is that the largest eigenvalue $v/(2\pi)^n = O(2^{pn})$, for $d^n/dx^n$, where $n=1, 2, 3$. It was also possible to confirm that the eigenvalues of the discretized matrix of odd-ordered differential operators have real parts of magnitude no greater than machine precision zero. The same can be said of the imaginary part of discretizations of even-ordered operators.

## 6. SUMMARY

The primary aim of this study has been to complement the report of Latto et al.¹³ in which they propose a simple and useful technique for the calculation of inner products involving Daubechies wavelets. It is hoped that this study will be a helpful guide to a publicly available numerical code which implements their scheme. Details on how to obtain the code are given in Appendix I.

The periodization of matrices involving inner products of two and three wavelets and their derivatives was presented in detail. The three wavelet inner product matrix is not straightforward to implement, especially in periodized problems. In this study we have provided several useful algorithms to construct such matrices and have listed several properties of the inner products that are useful in checking the actual implementation of the wavelet-Galerkin solution of differential equations.

Lastly, we have shown numerically that common differential operators may be well approximated using wavelets. In light of other studies that compare several spectral approximations to the above operators,²³ the wavelet approximations compare favourably. Information on the approximation of differential operators is also useful in the solution of partial differential equations using wavelet-Galerkin techniques, for example, in making choices of time integrators in the solution of evolution equations, and in estimating the size of the time stepping in order to achieve time-stability in the calculation. It was found that the eigenvalues, which approximate commonly used differential operators, where fairly close in value and magnitude to those calculated by means of discrete Fourier projection techniques.

# APPENDIX I

## Tables of connection coefficients

Connection coefficients: 2-tuples for $p=0, N=6$. Computed in double precision on a Sun Sparc1 workstation using LAPACK solver routines.

$$ \Lambda_{-4}^{0,1} = -3.4246575342471D - 04 $$

$$ \Lambda_{-3}^{0,1} = -1.4611872146119D - 02 $$
$$ \Lambda_{-2}^{0,1} = +0.14520547945205 $$

$$ \Lambda_{-1}^{0,1} = -0.74520547945205 $$

$$ \Lambda_{0}^{0,1} = -3.2049276679778D - 15 $$

$$ \Lambda_{1}^{0,1} = +0.74520547945206 $$

$$ \Lambda_{2}^{0,1} = -0.14520547945205 $$

$$ \Lambda_{3}^{0,1} = +1.4611872146119D - 02 $$

$$ \Lambda_{4}^{0,1} = +3.4246575342476D - 04 $$

$$ \text{residual} = 1.9680979936043D - 16 $$

for the least-squares solution of the overdetermined system.

$$ \Lambda_{-4}^{1,1} = +5.3571428571412D - 03 $$

$$ \Lambda_{-3}^{1,1} = +0.11428571428572 $$

$$ \Lambda_{-2}^{1,1} = -0.87619047619048 $$

$$ \Lambda_{-1}^{1,1} = +3.3904761904762 $$

$$ \Lambda_{0}^{1,1} = -5.2678571428572 $$

$$ \Lambda_{1}^{1,1} = +3.3904761904762 $$

$$ \Lambda_{2}^{1,1} = -0.87619047619048 $$

$$ \Lambda_{3}^{1,1} = +0.11428571428571 $$

$$ \Lambda_{4}^{1,1} = +5.3571428571430D - 03 $$

$$ \text{residual} = 1.1362438767648D - 16 $$

$$ \Lambda_{-4}^{2,0} = +5.3571428571412D - 03 $$

$$ \Lambda_{-3}^{2,0} = +0.11428571428572 $$

$$ \Lambda_{-2}^{2,0} = -0.87619047619048 $$

$$ \Lambda_{-1}^{2,0} = +3.3904761904762 $$

$$ \Lambda_{0}^{2,0} = -5.2678571428572 $$

$$ \Lambda_{1}^{2,0} = +3.3904761904762 $$

$$ \Lambda_{2}^{2,0} = -0.87619047619048 $$

$$ \Lambda_{3}^{2,0} = +0.11428571428571 $$

$$ \Lambda_{4}^{2,0} = +5.3571428571430D - 03 $$

$$ \text{residual} = 1.1362438767648D - 16 $$
A code for the computation of connection coefficients is available. The URL is http://www.math.ucla.edu/~restrepo.

APPENDIX II

Proofs related to periodizing Daubechies wavelets are presented in this appendix.

*Proposition.* To prove that dilation and periodization of Daubechies wavelets do not commute it suffices to show that for $j \le 0$, $\tilde{\phi}_{j,k} = 2^{-j/2}$.

*Proof.* Since

$$ 
\begin{aligned}
\tilde{\phi}_{j,k} &= 2^{j/2} \tilde{\phi}(2^j x - k) \\
\tilde{\phi}_{j,k} &= \sum_{l \in \mathbb{Z}} 2^{j/2} \varphi(2^j(x-l) - k) \\
&= 2^{j/2} \sum_{b=0}^{2^{-j}-1} \sum_{l \in \mathbb{Z}} \varphi \left( 2^j x - \left( l + \frac{b}{2^{-j}} \right) - k \right)
\end{aligned}
$$

Letting $y = 2^j x$ and summing over $l$, we obtain

$$ \tilde{\phi}_{j,k} = 2^{j/2} \sum_{b=0}^{2^{-j}-1} 1 = 2^{-j/2} \quad \square $$

Periodized Daubechies wavelets form an orthonormal basis for $L^2[0, 1]$:

*Proposition.* {$\tilde{\phi}_{j,k}: j \ge 0, 0 \le k < 2^j$} forms an orthonormal basis for $L^2[0, 1]$.

*Definition.* $\langle f, g \rangle = \int f(x)g(x)dx$, the standard $L^2$ inner product.

*Proof.* We begin by showing that $\langle \tilde{\phi}_{j,k}, \tilde{\phi}^{j'}, k' \rangle = 0$:

$$ \langle \tilde{\phi}_{j,k}, \tilde{\phi}^{j'}, k' \rangle = 2^{(j+j')/2} \int_0^1 \sum_{l,l' \in \mathbb{Z}} \varphi(2^j(x-l)-k)\varphi(2^{j'}(x-l')-k') dx $$

Let $y = x - l'$, so that

$$ 
\begin{aligned}
\langle \tilde{\phi}_{j,k}, \tilde{\phi}^{j'}, k' \rangle &= 2^{(j+j')/2} \int_0^1 \sum_{l,l' \in \mathbb{Z}} \varphi(2^j(x-l)-k)\varphi(2^{j'}(x-l')-k') dx \\
&= 2^{(j+j')/2} \sum_{r \in \mathbb{Z}} \int_{-\infty}^{\infty} \varphi(2^j y + 2^j r - k)\varphi(2^{j'} y - k') dx \\
&= \sum_{r \in \mathbb{Z}} \langle \varphi_{j,k+2^j r}, \varphi_{j',k'} \rangle = \delta_{jj'} \delta_{kk'}
\end{aligned}
$$
with $l - l' = r$. Thus $\tilde{\phi}_{j,k}$ and $\tilde{\phi}_{j',k'}$ are orthonormal. Next we show that they form a basis for $L^2[0, 1]$.

Choose an arbitrary $f \in L^2[0, 1]$. Now consider

$$\tilde{f}(x) = f(x)x \in [0, 1] = 0x \notin [0, 1]$$

$\tilde{f} \in L^2(\mathbb{R})$ and $\{\phi_{j,k}\}$ form an orthonormal basis for $L^2(\mathbb{R})$, so we have

$$\tilde{f} = \sum_{\substack{0 \le j \\ k \in \mathbb{Z}}} \langle \tilde{f}, \phi_{j,k} \rangle$$

which, when periodized, becomes

$$f(x) = \sum_{l \in \mathbb{Z}} \tilde{f}(x-l) = \sum_{l \in \mathbb{Z}} \sum_{\substack{0 \le j \\ k \in \mathbb{Z}}} \langle \tilde{f}, \phi_{j,k} \rangle = \sum_{\substack{0 \le j \\ k \in \mathbb{Z}}} \langle \tilde{f}, \tilde{\phi}_{j,k} \rangle$$

This final result is actually a finite sum since, for $k \ge 2^j$ and $k \le 1-N$, $\text{supp}(\phi_{j,k}) \cap \text{supp}(\tilde{f}) = \emptyset$. Thus $f$ has a representation in the periodized wavelets. $\square$

The proof that $\{\tilde{\phi}_{j,k}, \tilde{\psi}_{j',k'} : j' \ge j \ge 0, 0 \le k < 2^j\}$ also form an orthonormal basis is nearly identical. Since $\phi(x) = \sum_{k=0}^{N-1} h_k \phi(2x - k)$, we can periodize both sides to get

$$\begin{aligned}
\tilde{\phi}(x) &= \sum_{l \in \mathbb{Z}} \phi(x-l) = \sum_{l} \sum_{k} h_k \phi(2(x-l) - k) \\
&= \sum_{k} \sum_{l} h_k \phi(2(x-l) - k) \\
&= \sum_{k} h_k \tilde{\phi}_{1,k}(x)
\end{aligned}$$

The connection coefficients code may be used for both periodized as well as the non-periodized Daubechies wavelets inner products. This result is formalized as follows.

**Proposition.** For $j \ge \log_2(N-1)$, $\tilde{\phi}_{j,0} = \phi_{j,0}$, where $\tilde{\phi}$ is extended to $\mathbb{R}$ by setting it to 0 away from the unit interval.

*Proof.* $\text{supp}(\phi_{j,0}) = [0, 2^{-j}(N-1)]$, so for $j \ge \log_2(N-1)$, $\text{supp}(\phi_{j,0}) = [0, \beta]$, $\beta \le 1$. Thus $\tilde{\phi}_{j,0}(x) = \sum_{l \in \mathbb{Z}} \phi_{j,0}(x-l) = \phi_{j,0}$. $\square$

# APPENDIX III

Several approximating characteristics in Section 3 are proved below. As we have already seen, periodized wavelets provide a basis for $L^2[0, 1]$ so we have $\|f - P_j f\|_2 \to 0$, as $j \to \infty$. This is a property of any orthonormal basis of $L^2$, but this particular periodized basis has some additional properties.
**Theorem.** If $f$ is a continuous function on the torus, then $\|f - P_j f\|_{\infty} \to 0$ as $j \to \infty$.

*Proof.* We begin the proof by showing that our projection operator is bounded. $P_j$ is an integral operator of the form $P_j f(x) = \int_0^1 \sum_{k=0}^{2^j-1} \tilde{\phi}_{j,k}(y) \tilde{\phi}_{j,k}(x) f(y) dy$. Thus,

$$ \begin{aligned} \|P_j\|_{\infty} &\le \sup_{x \in [0, 1]} \int_0^1 \left| \sum_{k=0}^{2^j-1} \tilde{\phi}_{j,k}(y) \tilde{\phi}_{j,k}(x) \right| dy \\ &\le \sup_{x \in [0, 1]} 2^{j/2} \left| \sum_{k=0}^{2^j-1} \tilde{\phi}(2^j x - k) \right| 2^{j/2} \|\tilde{\phi}\|_{\infty} 2^{j/2} [2^{-j}(N-1)] \\ &\le \sup_{x \in [0, 1]} \left| \sum_{k=0}^{2^j-1} \tilde{\phi}(2^j x - k) \right| \|\tilde{\phi}\|_{\infty} (N-1) \end{aligned} $$

Now, $|\sum_{k=0}^{2^j-1} \tilde{\phi}(2^j x - k)| < (N-1)\|\tilde{\phi}\|_{\infty}$ since for $j \ge \log_2(N-1)$ there are at most $(N-1)$ $k$'s such that for a given $x$, $\{x\} \cap \text{supp}(\tilde{\phi}_{j,k}) \ne \emptyset$. Hence, we have

$$ \|P_j\|_{\infty} \le \|\tilde{\phi}\|_{\infty}^2 (N-1)^2 $$

If we take $f \in \bigcup_{j \in \mathbb{N}} V_j$, then $\exists J$ such that $\forall j \ge J$, $Q_j f = 0$. Thus, $P_j f = f$ for $j \ge J$. $\bigcup_{j \in \mathbb{N}} V_j$ is dense in $L^2[0, 1]$ which is dense in $C(\mathbb{T})$, continuous functions of period 1 on the unit interval. Finally, by the boundedness of $P_j$, the theorem follows. $\square$

**Theorem (Daubechies).** If $f \in L^1[0, 1]$, then $\|f - P_j f\|_1 \to 0$ as $j \to \infty$ [10].

*Proof.* Since $L^1[0, 1] \subset (C[0, 1])^*$, we have

$$ \begin{aligned} \|P_j f\|_1 &= \sup\{|\langle P_j f, g \rangle|; g \text{ continuous, } \|g\|_{\infty} \le 1\} \\ &= \sup\{|\langle f, P_j g \rangle|; g \text{ continuous, } \|g\|_{\infty} \le 1\} \\ &\le \|f\|_1 \|P_j g\|_{\infty} \end{aligned} $$

$\|P_j\|_{\infty}$ is bounded by the previous theorem. Since $\bigcup_{j \in \mathbb{N}} V_j$ is dense in $L^2[0, 1]$, which is also dense in $L^1[0, 1]$, the uniform bound on $P_j$ is sufficient to prove our result. $\square$

ACKNOWLEDGEMENTS

Special thanks to George Schlossnagle for his help at various stages of this work. The authors thank John Weiss for making his preprints available to us. Gilbert Strang provided extremely helpful comments on a draft of this report. This research was supported in part by an appointment to the Distinguished Postdoctoral Research Program sponsored by the U.S. Department of Energy, Office of University and Science Education Programs, and administered by the Oak Ridge Institute for Science and Education. Further support was provided by the Atmospheric and Climate Research Division and the Office of Scientific Computing, U.S. Department of Energy, under Contract W-31-109-Eng-38.
REFERENCES

1. C. Chui, *An Introduction to Wavelets*, Academic Press, New York, 1992.

2. I. Daubechies, *Ten Lectures on Wavelets*, Regional Conference Series in Applied Mathematics, Vol. 61, SIAM, Philadelphia, 1992.

3. Y. Meyer, *Ondelettes et Operateurs I*, Editeurs des Sciences et des Arts, Hermann, 1990.

4. G. Strang and T. Nguyen, *Wavelets and Filter Banks*, Addison Wesley, Reading, MA, U.S.A., 1995.

5. A. N. Akansu and R. A. Haddad, *Multiresolution Signal Decomposition: Transforms, Subbands, and Wavelets*, Academic Press, New York, 1992.

6. R. DeVore and B. J. Lucier, 'Wavelets', *Acta Numerica* **1**, 1–56 (1991).

7. W. Sweldens and R. Piessens, 'Quadrature formulae for the calculation of the wavelet decomposition', preprint, 1991.

8. E. Bacry, S. Mallat and G. Papanicolau, 'A wavelet based space-time adaptive numerical method for partial differential equations', *Math. Modell. Numer. Anal.*, **26**, 793–834 (1992).

9. J. M. Restrepo and G. K. Leaf, 'Wavelet-Galerkin discretization of hyperbolic equations', *J. Comput. Phys.*, **122** (1995) 118–128.

10. M. Farge, J. C. R. Hunt and J. C. Vassilicos (eds.), 'Wavelets, fractals, and fourier transforms', in: *Proc. Conf. on Wavelets, Fractals and Fourier Transforms: New Developments and New Applications*, Institute of Mathematics and Applications and Société de Mathématiques Appliquées et Industrielles and held at Newnham College, Cambridge, December 1990, Clarendon Press, Oxford, 1993.

11. W. Press, S. Teukolsky, W. Vetterling and B. Flannery, 'Numerical recipes in Fortran', in *The Art of Scientific Computing*, Cambridge University Press, Cambridge, 1992.

12. K. Amaratunga, J. R. Williams, S. Qian and J. Weiss, 'Wavelet-Galerkin solution of one-dimensional pde's', *Int. J. Numer. Meth. Engng.*, **37**, 2703–2716 (1994).

13. A. Latto, H. L. Resnikoff and E. Tenenbaum, 'The evaluation of connection coefficients of compactly supported wavelets', *Proc. French-USA Workshop on Wavelets and Turbulence*, Princeton, 1991, Springer, New York.

14. W. Dahmen and C. A. Micchelli, 'Using the refinement equation for evaluating integrals of wavelets', *SIAM J. Numer. Anal.*, **30**, 507–537 (1993).

15. A. Kunoth, 'Computing refinable integrals: documentation of the program, version 1.1', Institute for Scientific Computation, Texas A & M University, *Technical Report ISC-95-02 MATH*, 1995.

16. A. Cohen, I. Daubechies and P. Vial, 'Wavelets on the interval and fast wavelet transforms', *Appl. Comput. Harmonic Anal.*, **1**, 54–81 (1993).

17. W. Rudin, *Real and Complex Analysis*, McGraw-Hill, New York, 1987.

18. C. Canuto, M. Y. Hussaini, A. Quarteroni and T. A. Zang, *Spectral Methods in Fluid Dynamics*, Springer Series in Computational Physics Springer, New York, 1988.

19. G. Beylkin, 'On the Representation of operators in bases of compactly supported wavelets', *SIAM J. Numer. Anal.* **6**, 1716–1740 (1992).

20. S. Qian and J. Weiss, 'Wavelets and the numerical solution of partial differential equations', *J. Comput. Phys.*, **106** (1993) 155–175.

21. P. J. Davis, *Circulant Matrices*, Wiley, New York, 1979.

22. E. Anderson et al., *LAPACK User's Guide*, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1992.

23. L. N. Trefethen and M. R. Trummer, 'An instability phenomenon in spectral methods', *SIAM J. Numer. Anal.*, **24**, 1008–1023 (1987).