Adaptive wavelet methods for elliptic operator equations: Convergence rates

Abstract

This paper is concerned with the construction and analysis of wavelet-based adaptive algorithms for the numerical solution of elliptic equations. These algorithms approximate the solution $u$ of the equation by a linear combination of $N$ wavelets. Therefore, a benchmark for their performance is provided by the rate of best approximation to $u$ by an arbitrary linear combination of $N$ wavelets (so called $N$-term approximation), which would be obtained by keeping the $N$ largest wavelet coefficients of the real solution (which of course is unknown). The main result of the paper is the construction of an adaptive scheme which produces an approximation to $u$ with error $O(N^{-s})$ in the energy norm, whenever such a rate is possible by $N$-term approximation. The range of $s>0$ for which this holds is only limited by the approximation properties of the wavelets together with their ability to compress the elliptic operator. Moreover, it is shown that the number of arithmetic operations needed to compute the approximate solution stays proportional to $N$. The adaptive algorithm applies to a wide class of elliptic problems and wavelet bases. The analysis in this paper puts forward new techniques for treating elliptic problems as well as the linear systems of equations that arise from the wavelet discretization.

1. Introduction

1.1. Background

Adaptive methods, such as adaptive finite elements methods (FEM), are frequently used to numerically solve elliptic equations when the solution is known to have singularities. A typical algorithm uses information gained during a given stage of the computation to produce a new mesh for the next iteration. Thus, the adaptive procedure depends on the current numerical resolution of $u$. Accordingly, these methods produce a form of nonlinear approximation of the solution, in contrast with linear methods in which the numerical procedure is set in advance and does not depend on the solution to be resolved.

The motivation for adaptive methods is that they provide flexibility to use finer resolution near singularities of the solution and thereby improve on the approximation efficiency. Since the startling papers 23 the understanding and practical realization of adaptive refinement schemes in a finite element context has been documented in numerous publications 3451336. A key ingredient in most adaptive algorithms are a-posteriori error estimators or indicators derived from the current residual or the solution of local problems. They consist of local quantities such as jumps of derivatives across the interface between adjacent triangles or simplices. One often succeeds in bounding the (global) error of the current solution with respect to the energy norm, say, by sums of these quantities from below and above. Thus refining the mesh where these local quantities are large is hoped to reduce the bounds and hence the error in the next computation. Computational experience frequently confirms the success of such techniques for elliptic boundary value problems in the sense that adaptively generated highly nonuniform meshes indeed give rise to an accuracy that would require the solution of much larger systems of equations based on uniform refinements. However, on a rigorous level the quantitative gain of adaptive techniques is usually not clear. The central question is whether the mesh refinements actually result, at each step, in some fixed error reduction. To our knowledge, only in 35 has convergence of an adaptive scheme been established for a rather special case, namely a piecewise linear finite element discretization of the classical Dirichlet problem for Laplace’s equation. There is usually no rigorous proof of the overall convergence of such schemes unless one assumes some quantitative information such as the saturation property about the unknown solution 13. Saturation properties are assumed but not proven to hold.

Moreover, the derivation of error indicators in conventional discretizations hinges on the locality of differential operators. Additional difficulties are therefore encountered when considering elliptic operators with nonlocal Schwartz kernel arising, for instance, in connection with boundary integral equations.

In summary, there seem to be at least two reasons for this state of affairs:

(i) There is an inherent difficulty, even for local operators, in utilizing the information available at a given stage in the adaptive computation to guarantee that a suitable reduction will occur in the residual error during the next adaptive step.

(ii) Finite element analysis is traditionally based on Sobolev regularity (see e.g., 14 or 15), which is known to govern the performance of linear methods. Only recent developments in the understanding of nonlinear methods have revealed that Besov regularity is a decidedly different and more appropriate smoothness scale for the analysis of adaptive schemes, see e.g., 31.

In view of the significant computational overhead and the severe complications caused by handling appropriate data structures for adaptive schemes, not only guaranteeing convergence but above all knowing its speed is of paramount importance for deciding whether or under which circumstances adaptive techniques actually pay off. To our knowledge nothing is known so far about the actual rate of convergence of adaptive FEM solvers, by which we mean the relation between the accuracy of the approximate solution and the involved degrees of freedom, or better yet the number of arithmetic operations.

1.2. Wavelet methods

An alternative to FEM are wavelet based methods. Similarily to mesh refinement in FEM, these methods offer the possibility to compress smooth functions with isolated singularities into high-order adaptive approximations involving a small number of basis functions. In addition, it has been recognized for some time 11 that for a large class of operators (including integral operators) wavelet bases give rise to matrix representations that are quasi-sparse (see Sections 2 and 3 for a definition of quasi-sparse) and admit simple diagonal preconditioners in the case of elliptic operators. Therefore, it is natural to develop adaptive strategies based on wavelet discretizations in order to solve elliptic operator equations numerically.

The state of wavelet-based solvers is still in its infancy, and certain inherent impediments to their numerical use remain. These are mainly due to the difficulty of dealing with realistic domain geometries. Nevertheless, these solvers show great promise, especially for adaptive approximation (see e.g., 112161825). Most adaptive strategies exploit the fact that wavelet coefficients convey detailed information on the local regularity of a function and thereby allow the detection of its singularities. The rule of thumb is that wherever wavelet coefficients of the currently computed solution are large in modulus, additional refinements are necessary. In some sense, this amounts to using the size of the computed coefficients as local a-posteriori error indicators. Note that here refinement has a somewhat different meaning than in the finite element setting. There the adapted spaces result from refining a mesh. The mesh is the primary controlling device and may create its own problems (of geometric nature) that have nothing to do with the underlying analytic task. In the wavelet context refinement means to add suitably selected further basis functions to those that are used to approximate the current solution. We refer to this as space refinement.

In spite of promising numerical performances, the problem remains (as in the finite element context) to quantify these strategies, that is, to decide which and how many additional wavelets need to be added in a refinement step in order to guarantee a fixed error reduction rate at the next resolution step. An adaptive wavelet scheme based on a-posteriori error estimators has been recently developed in 20, which ensures this fixed error reduction for a wide class of elliptic operators, including those of negative order. This shows that making use of the characteristic features of wavelet expansions, such as the sparsification and preconditioning of elliptic operators, allows one to go beyond what is typically known in the conventional framework of adaptive FEM. However, similar to FEM, there are so far no results about the rate of convergence of adaptive wavelet based solvers, i.e., the dependence of the error on the number of degrees of freedom.

1.3. The objectives

The purpose of the present paper is twofold. First, we provide analytical tools that can be utilized in studying the theoretical performance of adaptive algorithms. Second, we show how these tools can be used to construct and analyze wavelet based adaptive algorithms which display optimal approximation and complexity properties in the sense that we describe below.

The adaptive methods we analyze in this paper take the following form. We assume that we have in hand a wavelet basis $\{\psi _\lambda \}_{\lambda \in \nabla }$ to be used for numerically resolving the elliptic equation. Our adaptive scheme will iteratively produce finite sets $\Lambda _j\subset \nabla$, $j=1,2,\dots$, and the Galerkin approximation $u_{\Lambda _j}$ to $u$ from the space $S_{\Lambda _j}:=\mathrm{span}(\{\psi _\lambda \}_{\lambda \in \Lambda _j})$. The function $u_{\Lambda _j}$ is a linear combination of $N_j:=\#\Lambda _j$ wavelets. Thus the adaptive method can be viewed as a particular form of nonlinear $N$-term wavelet approximation, and a benchmark for the performance of such an adaptive method is provided by comparison with best $N$-term approximation (in the energy norm) when full knowledge of $u$ is available.

Much is known about $N$-term approximation. In particular, there is a characterization of the functions $v$ that can be approximated in the energy norm with accuracy $O(N^{-s})$ by using linear combinations of $N$ wavelets. As we already mentioned, this class $B^s$ is typically a Besov space, which is substantially larger than the corresponding Sobolev space $W^s$ which ensures $O(N^{-s})$ accuracy for uniform discretization with $N$ parameters. In several instances of the elliptic problems, e.g., when the right hand side $f$ has singularities, or when the boundary of $\Omega$ has corners, the Besov regularity of the solution will exceed its Sobolev regularity (see 19 and 21). So these solutions can be approximated better by best $N$-term approximation than by uniformly refined spaces, and the use of adaptive methods is suggested. Another important feature of $N$-term approximation is that a near best approximation is produced by thresholding, i.e., simply keeping the $N$ largest contributions (measured in the same metric as the approximation error) of the wavelet expansion of $v$.

Of course, since best $N$-term approximation requires complete information on the approximated function, it cannot be applied directly to the unknown solution. It is certainly not clear beforehand whether a concrete numerical scheme can produce at least asymptotically the same convergence rate. Thus ideally an optimal adaptive wavelet algorithm should produce a result similar to thresholding the exact solution. In more quantitative terms this means that whenever the solution $u$ is in $B^s$, the approximations $u_{\Lambda _j}$ should satisfy

$$\begin{equation} \|u-u_{\Lambda _j}\|\le C\|u\|_{B^s}N_j^{-s},\quad N_j:=\#\Lambda _j, {\tag{1.1}} \end{equation}$$

where $\|\cdot \|$ is the energy norm and $\|\cdot \|_{B^s}$ is the norm for $B^s$. Since in practice one is mostly interested in controlling a prescribed accuracy with a minimal number of parameters, we shall rather say that the adaptive algorithm is of optimal order $s>0$ if whenever the solution $u$ is in $B^s$, then for all $\epsilon >0$, there exists $j(\epsilon )$ such that

$$\begin{equation} \|u-u_{\Lambda _j}\|\le \epsilon , \;\; j\geq j(\epsilon ), {\tag{1.2}} \end{equation}$$

and such that

$$\begin{equation} \#(\Lambda _{j(\epsilon )})\le C\|u\|_{B^s}^{1/s}\epsilon ^{-1/s}. {\tag{1.3}} \end{equation}$$

Such a property ensures an optimal memory size for the description of the approximate solution.

Another crucial aspect of the adaptive algorithms is their computational complexity: we shall say that the adaptive algorithm is computationally optimal if, in addition to 1.2–1.3, the number of arithmetic operation needed to derive $u_{\Lambda _j}$ is proportional to $\#\Lambda _{j}$. Note that an instance of computational optimality in the context of linear methods is provided by the full multigrid algorithm when $N$ represents the number of unknowns necessary to achieve a given accuracy on a uniform grid. We are thus interested in algorithms that exhibit the same type of computational optimality with respect to an optimal adaptive grid which is not known in advance and should itself be generated by the algorithm.

The main accomplishment of this paper is the development and analysis of an adaptive numerical scheme which for a wide class of operator equations (including those of negative order) is optimal with regard to best $N$-term approximation and is also computationally optimal in the above sense. Let us mention that a simplified version of this algorithm has been developed and tested in 18, as well as a more elaborate version in 6. In this last version, appropriate object-oriented data structures play a crucial role for handling properly the sparse representation of the solution. In both cases, the numerical experiments confirm a similar behavior between the numerical error generated by such adaptive algorithms and by thresholding the exact solution. We note however, that depending on the concrete application at hand, these implementations still suffer from currently suboptimal schemes for a central task, namely evaluating the entries of the stiffness matrices. Since, in the case of piecewise polynomial wavelets and (piecewise) constant coefficients PDE’s, these elements are given by explicit formulae, the cost of this task is negligible. However, a much more sophisticated treatment appears to be necessary in the general case where, for instance, the techniques recently proposed in 10 should lead to significant improvements.

Although the present investigations are confined to symmetric elliptic problems, the results provide, in our opinion, a core ingredient for the treatment of more complex tasks. For instance, the first steps in this direction are the development of a-posteriori wavelet strategies for saddle point problems in 22, the use of stabilized variational formulations in 9, or least squares formulations in 24.

1.4. Organization of the paper

In Section 2, we introduce the general setting of elliptic operator equations where our results apply. In this context, after applying a diagonal preconditioner, wavelet discretizations allow us to view the equation as a discrete well conditioned $\ell _2$ linear system.

In Section 3, we review certain aspects of nonlinear approximation, quasi-sparse matrices and fast multiplication using such matrices. The main result of this section is an algorithm for the fast computation of the application of a quasi-sparse matrix to a vector.

In Section 4, we analyze the rate of convergence of the refinement procedure introduced earlier in 20. We will refer to this scheme here as Algorithm I. We show that this algorithm is optimal for a small range of $s>0$. However, the full range of optimality should be limited only by the properties of the wavelet basis (smoothness and vanishing moments) and the operator; this is not the case for Algorithm I. The analysis in Section 4, however, identifies the barrier that keeps Algorithm I from being optimal in the full range of $s$.

In Section 5, we introduce a second strategy—Algorithm II—for adaptively generating the sets $\Lambda _j$ that is shown to provide optimal approximation of order $s>0$ for the full range of $s$. The new ingredient that distinguishes Algorithm II from Algorithm I is the addition of thresholding steps which delete some indices from $\Lambda _j$. This would be the analogue of coarsening the mesh in FEM.

Although we have qualified so far both procedures in Sections 4 and 5 as “algorithms”, we have actually ignored any issue concerning practical realization. They are idealized in the sense that the exact assessment of residuals and Galerkin solutions is assumed. This was done in order to clearly identify the essential analytical tasks. Practical realizations require truncations and approximations of these quantities. Section 6 is devoted to developing the ingredients of a realistic numerical scheme. This includes quantitative thresholding procedures, approximate matrix/vector multiplication, approximate Galerkin solvers and the approximate evaluation of residuals.

In Section 7 we employ these ingredients to formulate a computable version of Algorithm II, which is shown to be computationally optimal for the full range of $s$. Recall that this means that it realizes for this range the order of best $N$-term approximation at the expense of a number of arithmetic operations that stays proportional to the number $N$ of significant coefficients. Computational optimality hinges to a great extent on the fast approximate matrix/vector multiplication from Section 3.

It should be noted however that an additional cost in our wavelet adaptive algorithm is incurred by sorting the coefficients in the currently computed solution. This cost at stage $j$ is of order $N\log N$, where $N=\#\Lambda _j$, thus slightly larger than the cost in arithmetic operations. It should be stressed that the complexity of the algorithm is analysed under the assumption that the solution exhibits a certain rate of best $N$-term approximation which is, for instance, implied by a certain Besov regularity. The algorithm itself does not require any a-priori assumption of that sort.

We have decided to carry out the (admittedly more technical) analysis of the numerical ingredients in some detail in order to substantiate our claim that the optimality analysis is not based on any hidden assumptions (beyond those hypotheses that are explicitly stated) such as accessing infinitely many data. Nevertheless the main message of this paper can be read in Sections 4 and 5: optimal adaptive approximations of elliptic equations can be computed by iterative wavelet refinements using a-posteriori error estimators, provided that the computed solution is regularly updated by appropriate thresholding procedures. This fact was already suggested by numerical experiments in 18 that show similar behavior between the numerical error generated by such adaptive algorithms and by thresholding the exact solution.

2. The setting

In this section, we shall introduce the setting in which our results apply. In essence, our analysis applies whenever the elliptic operator equation takes place on a manifold or domain which admits a biorthogonal wavelet basis.

2.1. Ellipticity assumptions

This subsection gives the assumptions we make on the operator equation to be numerically solved. These assumptions are quite mild and apply in great generality.

Let $\Omega$ denote a bounded open domain in the Euclidean space $\mathbb{R}^{d}$ with Lipschitz boundary or, more generally, a Lipschitz manifold of dimension $d$. In particular, $\Omega$ could be a closed surface which arises as a domain for a boundary integral equation. The space $L_{2}(\Omega )$ consists of all square integrable functions on $\Omega$ with respect to the (canonically induced) Lebesgue measure. The corresponding inner product is denoted by

$$\begin{equation} \langle \cdot ,\cdot \rangle _{L_2(\Omega )}. {\tag{2.1}} \end{equation}$$

Let $A$ be a linear operator mapping a Hilbert space $H$ into $H^{*}$ (its dual relative to the pairing $\langle \cdot ,\cdot \rangle _{L_2(\Omega )}$), where $H$ is a space with the property that either $H$ or its dual $H^{*}$ is embedded in $L_{2}(\Omega )$. The operator $A$ induces the bilinear form $a$ defined on $H\times H$ by

$$\begin{equation} a(u ,v):=\langle Au,v\rangle , {\tag{2.2}} \end{equation}$$

where $\langle \cdot ,\cdot \rangle$ denotes the $(H^*,H)$ duality product.

(A1): We assume that the bilinear form $a$ is symmetric positive definite and elliptic in the sense that

$$\begin{equation} a(v,v)\sim \|v\|_{H}^2 ,\quad v\in H. {\tag{2.3}} \end{equation}$$

Here, and throughout this paper, $\sim$ means that both quantities can be uniformly bounded by constant multiples of each other. Likewise $\lesssim$ indicates inequalities up to constant factors.

It follows that $H$ is also a Hilbert space with respect to the inner product $a$ and that this inner product induces an equivalent norm (called the energy norm) on $H$ by

$$\begin{equation} \|\cdot \|_a^{2} :=a(\cdot ,\cdot ). {\tag{2.4}} \end{equation}$$

By duality, $A$ thus defines an isomorphism from $H$ onto $H^{*}$. We shall study the equation

$$\begin{equation} Au=f {\tag{2.5}} \end{equation}$$

with $f\in H^{*}$. From our assumptions, it follows that for any $f\in H^{*}$, this equation has a unique solution in $H$, which will always be denoted by $u$. This is also the unique solution of the variational equation

$$\begin{equation} a(u,v)=\langle f,v\rangle ,\;\; \text{for}\;\;\text{all}\;\; v\in H. {\tag{2.6}} \end{equation}$$

The typical examples included in the above assumptions are the Poisson or the biharmonic equations on bounded domains in $\mathbb{R}^{d}$; single or double layer potentials, and hypersingular operators on closed surfaces arising in the context of boundary integral equations. In these examples $H$ is a Sobolev space, e.g., $H=H^{1}_{0}(\Omega )$, $H^{2}_{0}(\Omega )$, or $H^{-1/2}(\Omega )$; see 232041 for examples.

2.2. Wavelet assumptions

By now wavelet bases are available for various types of domains that are relevant for the formulation of operator equations. This covers, for instance, polyhedral surfaces of dimension two and three 29 as well as manifolds or domains that can be represented as a disjoint union of smooth regular parametric images of a simple parameter domain such as the unit cube 27.

There are many excellent accounts of wavelets on $\mathbb{R}^d$ (see e.g., 38 or 30). For the construction and description of wavelet bases on domains and manifolds, we refer the reader to the survey paper 23 and the references therein. This survey also sets forth the notation we shall employ below for indexing the elements in a wavelet basis. To understand this notation, it may be useful for the reader to keep in mind the case of wavelet bases on $\mathbb{R}^d$. In this setting, a typical biorthogonal wavelet basis of compactly supported functions is given by the shifted dilates of a set $\Gamma$ of $2^d-1$ functions. Namely, the collection of functions

$$\begin{equation} 2^{jd/2}\gamma (2^j\cdot -k),\quad j\in \mathbb{Z},\ k\in \mathbb{Z}^d,\ \gamma \in \Gamma , {\tag{2.7}} \end{equation}$$

forms a Riesz basis for $L_2(\mathbb{R}^d)$. The dual basis is given by

$$\begin{equation} 2^{jd/2}\tilde{\gamma }(2^j\cdot -k),\quad j\in \mathbb{Z},\ k\in \mathbb{Z}^d,\ \tilde{\gamma }\in \tilde{\Gamma }, {\tag{2.8}} \end{equation}$$

with $\tilde{\Gamma }$ again a set of $2^{d}-1$ functions. The integer $j$ gives the dyadic level ($2^j$ the frequency) of the wavelet. The multi-integer $k$ gives the position ($2^{-j}k$) of the wavelet. Namely, the wavelet has support contained in a cube of diameter $\lesssim 2^{-j}$ centered at the point $2^{-j}k$. Note that there are $2^d-1$ functions with the same dyadic level $j$ and position $2^{-j}k$.

Another way to construct a wavelet basis for $\mathbb{R}^d$ is to start the multiscale decomposition at a finite dyadic level $j_0$. In this case, the basis consists of the functions of 2.7 with $j\ge j_0$, together with a family of functions

$$\begin{equation} 2^{j_0d/2}\phi (2^{j_0}\cdot -k),\quad k\in \mathbb{Z}^d, {\tag{2.9}} \end{equation}$$

with $\phi$ a fixed (scaling) function. Wavelet bases for domains take a similar form except that some alterations are necessary near the boundary.

We shall denote wavelet bases by $\{\psi _\lambda \}_{\lambda \in \nabla }$. In the particular case above, this notation incorporates the three parameters $j,k,\gamma$ into the one $\lambda$. We use $|\lambda |:=j$ to denote the dyadic level of the wavelet. We let $\Psi _{j}=\{\psi _{\lambda }:\lambda \in \nabla _{j}\}$, $\nabla _{j}:=\{\lambda \in \nabla : |\lambda |=j\}$, consist of the wavelets at level $j$.

In all classical constructions of compactly supported wavelets, there exist fixed constants $C$ and $M$ such that $\mathrm{diam}(\operatorname {supp}(\psi _\lambda ))\leq C2^{-|\lambda |}$ and such that for all $\lambda \in \nabla _j$ there are at most $M$ indices $\mu \in \nabla _j$ such that $\mathrm{meas}(\operatorname {supp}(\psi _\lambda )\cap \operatorname {supp}(\psi _\mu ))\neq 0$.

Since we shall consider only bounded domains in this paper, the wavelet decomposition will begin at some fixed level $j_0$. For notational convenience only, we assume $j_0=1$. We define $\Psi _{0}$ to be the set of scaling functions in the wavelet basis. We shall assume that $\Omega$ is a domain or manifold which admits two sets of functions:

$$\begin{equation} \Psi =\{\psi _{\lambda }: \lambda \in \nabla \} \subset L_{2}(\Omega ), \quad \tilde{\Psi }=\{\tilde{\psi }_{\lambda }: \lambda \in \nabla \}\subset L_{2}(\Omega ) {\tag{2.10}} \end{equation}$$

that form a biorthogonal wavelet bases on $\Omega$: writing $\langle \Theta ,\Phi \rangle :=(\langle \theta ,\phi \rangle _{L_2(\Omega )})_{\theta \in \Theta , \phi \in \Phi }$ for any two collections $\Theta ,\Phi$ of functions in $L_{2}(\Omega )$, one has

$$\begin{equation} \langle \Psi ,\tilde{\Psi }\rangle =\mathbf{I}, {\tag{2.11}} \end{equation}$$

where $\mathbf{I}$ is the identity matrix.

A typical feature in the theory of biorthogonal bases is that the sequences $\Psi ,\tilde{\Psi }$ are Riesz bases. That is, using the shorthand notation $\mathbf{d}^T\Psi :=\sum _{\lambda \in \nabla }d_\lambda \psi _{\lambda }$, one has

$$\begin{equation} \|\mathbf{d}\|_{\ell _{2}(\nabla )}\sim \|\mathbf{d}^{T}\Psi \|_{L_{2}(\Omega )} \sim \|\mathbf{d}^{T}\tilde{\Psi }\|_{L_{2}(\Omega )} . {\tag{2.12}} \end{equation}$$

This property means that the wavelet bases characterize $L_2(\Omega )$. In the present context of elliptic equations, we shall need not 2.12 but rather the fact that these bases provide a characterization of $H$ and $H^{*}$ in terms of wavelet coefficients. This is expressed by the following specific assumption.

(A2): Let the energy space $H$ be equipped with the norm $\|\cdot \|_H$ and its dual space $H^{*}$ with the norm $\|v \|_{H^*}:=\sup _{\|w\|_H=1}|\langle v,w\rangle |$. We assume that the wavelets in $\Psi$ are in $H$, whereas those in $\tilde{\Psi }$ are in $H^*$ (in this context, we can assume that 2.11 simply holds in the sense of the duality $(H,H^*)$). We assume that each $v\in H$ has a wavelet expansion $v= \mathbf{d}^{T}\Psi$ (with coordinates $d_\lambda =\langle v,\tilde{\psi }_\lambda \rangle$) and that

$$\begin{equation} \|\mathbf{D}^{-1} \mathbf{d}\|_{\ell _{2}(\nabla )}\sim \|\mathbf{d}^{T}\Psi \|_H. {\tag{2.13}} \end{equation}$$

with $\mathbf{D}$ a fixed positive diagonal matrix.

Observe that 2.13 implies that $\mathbf{D}_{\lambda ,\lambda }\sim \|\psi _\lambda \|_H^{-1}$, and that $\Psi$ (resp. $D^{-1}\Psi$) is an unconditional (resp. Riesz) basis for $H$. By duality, one easily obtains that each $v\in H^*$ has a wavelet expansion $v=\mathbf{d}^{T}\tilde{\Psi }$ (with coordinates $d_\lambda =\langle v,\psi _\lambda \rangle$) that satisfies

$$\begin{equation} \|\mathbf{D}\mathbf{d}\|_{\ell _{2}(\nabla )}\sim \|\mathbf{d}^{T}\tilde{\Psi }\|_{H^*}. {\tag{2.14}} \end{equation}$$

One should keep in mind, though, that $\tilde{\Psi }$ is only needed for analysis purposes. The Galerkin schemes to be considered below only involve $\Psi$, while $\tilde{\Psi }$ never enters any computation and need not even be known explicitly.

It is well known (see e.g., 27) that wavelet bases provide such characterizations for a large variety of spaces (in particular the Sobolev and Besov spaces for a certain parameter range which depends on the smoothness of the wavelets). In the context of elliptic equations, $H$ is typically some Sobolev space $H^{t}$. In this case (A2) is satisfied whenever the wavelets are sufficiently smooth, with $\mathbf{D}_{\lambda ,\lambda }=2^{-|\lambda |t}$. For instance, when $A=-\Delta$, one has $t=1$.

2.3. Discretization and preconditioning of the elliptic equation

Using wavelets, we can rewrite 2.5 as an infinite system of linear equations. We take wavelet bases $\Psi$ and $\tilde{\Psi }$ satisfying (A2) and write the unknown solution $u=\mathbf{d}^{T}\Psi$ and the given right hand side $f$ in terms of the basis $\tilde{\Psi }$. This gives the system of equations

$$\begin{equation} \langle A\Psi ,\Psi \rangle ^{T}\mathbf{d}= \langle f,\Psi \rangle ^{T}. {\tag{2.15}} \end{equation}$$

The solution $\mathbf{d}$ to 2.15 gives the wavelet coefficients of the solution $u$ to 2.5.

An advantage of wavelet bases is that they allow for trivial preconditioning of the linear system 2.15. This preconditioning is given by the matrix $\mathbf{D}$ of (A2) and results in the system of equations

$$\begin{equation} \mathbf{D}\langle A\Psi ,\Psi \rangle ^{T}\mathbf{D}\mathbf{D}^{-1}\mathbf{d}= \mathbf{D}\langle f,\Psi \rangle ^{T}, {\tag{2.16}} \end{equation}$$

or, more compactly,

$$\begin{equation} \mathbf{A}\mathbf{u}={\mathbf{f}}, {\tag{2.17}} \end{equation}$$

where

$$\begin{equation} \mathbf{A}:=\mathbf{D}\langle A\Psi ,\Psi \rangle ^{T}\mathbf{D},\quad \mathbf{u}:= \mathbf{D}^{-1}\mathbf{d},\quad \mathbf{ f}:= \mathbf{D}\langle f,\Psi \rangle ^{T} \in \ell _2(\nabla ). {\tag{2.18}} \end{equation}$$

Let us briefly explain the effect of the above diagonal scaling with regard to preconditioning. To this end, note that by (A1), the matrix $\mathbf{A}$ is symmetric positive definite. We define its associated bilinear form $\mathbf{a}$ by

$$\begin{equation} \mathbf{a}(\mathbf{v},\mathbf{w}):=\langle \mathbf{A}\mathbf{v},\mathbf{w}\rangle _{\ell _2(\nabla )}, {\tag{2.19}} \end{equation}$$

where $\langle \cdot ,\cdot \rangle _{\ell _{2}(\nabla )}$ is the standard inner product in $\ell _{2}(\nabla )$, and denote the norm associated with this bilinear form by $\|\cdot \|$. In other words,

$$\begin{equation} \|\mathbf{v}\|^2:=\mathbf{a}(\mathbf{v},\mathbf{v}),\quad \mathbf{v}\in \ell _2(\nabla ). {\tag{2.20}} \end{equation}$$

Combining the ellipticity assumption (A1) together with the wavelet characterization of $H$ (A2), we obtain that $\|\cdot \|$ and $\|\cdot \|_{\ell _2(\nabla )}$ are equivalent norms, i.e., there exist constants $c_1,c_2>0$ such that

$$\begin{equation} c_1\|\mathbf{v}\|^2_{\ell _2(\nabla )}\le \|\mathbf{v}\|^2 \le c_2 \|\mathbf{v}\|^2_{\ell _2(\nabla )}. {\tag{2.21}} \end{equation}$$

It is immediate that these constants are also such that

$$\begin{equation} c_1\|\mathbf{v}\|_{\ell _2(\nabla )}\le \|\mathbf{A}\mathbf{v}\|_{\ell _2(\nabla )}\le c_2 \|\mathbf{v}\|_{\ell _2(\nabla )} {\tag{2.22}} \end{equation}$$

and

$$\begin{equation} c_2^{-1}\|\mathbf{v}\|_{\ell _2(\nabla )}\le \|\mathbf{A}^{-1}\mathbf{v}\|_{\ell _2(\nabla )}\le c_1^{-1} \|\mathbf{v}\|_{\ell _2(\nabla )}. {\tag{2.23}} \end{equation}$$

In particular, the condition number $\kappa :=\|\mathbf{A}\|\|\mathbf{A}^{-1}\|$ of $\mathbf{A}$ satisfies

$$\begin{equation} \kappa \le c_2c_1^{-1}. {\tag{2.24}} \end{equation}$$

The fact that the diagonal scaling turns the original operator into an isomorphism on $\ell _{2}(\nabla )$ will be a cornerstone of the subsequent development. Denoting by $a_{\lambda ,\lambda '}$ the entries of $\mathbf{A}$ and by $\mathbf{A}_{\Lambda }=(a_{\lambda ,\lambda '})_{\lambda ,\lambda '\in \Lambda }$ the section of $\mathbf{A}$ restricted to the set $\Lambda$, we see from the positive definiteness of $\mathbf{A}$ that

$$\begin{equation} \|\mathbf{A}_\Lambda \|\le \|\mathbf{A}\|, \quad \|\mathbf{A}_\Lambda ^{-1}\| \le \|\mathbf{A}^{-1}\|, {\tag{2.25}} \end{equation}$$

and that the condition numbers of the submatrices remain uniformly bounded for any subset $\Lambda \subset \nabla$, i.e.,

$$\begin{equation} \kappa (\mathbf{A}_{\Lambda }):= \|\mathbf{A}_{\Lambda }\|\|\mathbf{A}_{\Lambda }^{-1}\| \leq \kappa . {\tag{2.26}} \end{equation}$$

Finally, it is easy to check that the constants $c_1$ and $c_2$ also provide the equivalence

$$\begin{equation} c_1^{1/2} \|\mathbf{v}\| \le \|\mathbf{A}\mathbf{v}\|_{\ell _2(\nabla )} \le c_2^{1/2} \|\mathbf{v}\| . {\tag{2.27}} \end{equation}$$

Here and later, we adopt the following rule about denoting constants. We shall denote constants which appear later in our analysis by $c_1, c_2,\cdots$. Other constants, whose value is not so important for us, will be denoted by $C$ or incorporated into the $\lesssim , \sim$ notation.

A typical instance of the above setting involves Sobolev spaces $H=H^t$, in which case the entries of the diagonal matrix $\mathbf{D}$ can be chosen as $2^{-t|\lambda |}\delta _{\lambda ,\lambda '}$. Of course, the constants in 2.24 will then depend on the relation between the energy norm 2.20 and the Sobolev norm. In some cases such a detour through a Sobolev space is not necessary, and 2.13 can be arranged to hold for a suitable $\mathbf{D}$ when $\|\cdot \|_H$ already coincides with the energy norm. A simple example is $A u =-\epsilon \Delta u +u$, where $(\mathbf{D})_{\lambda ,\lambda '}:= \max \,\{1,\sqrt {\epsilon }2^{|\lambda |}\} \delta _{\lambda ,\lambda '}$ is an appropriate choice. In fact, 2.13 will then hold independently of $\epsilon$.

2.4. Quasi-sparsity assumptions on the stiffness matrix

Another advantage of the wavelet basis is that for a large class of elliptic operators, the resulting preconditioned matrix $\mathbf{A}$ exhibits fast decay away from the diagonal. This will later be crucial with regard to storage economy and efficiency of (approximate) matrix/vector multiplication.

Consider for example the (typical) case when $H$ is the Sobolev space $H^t$ of order $t$ or its subspace $H^t_0$. Then, for a large class of elliptic operators, we have

$$\begin{equation} 2^{-(|\lambda '|+|\lambda |)t}|\langle A\psi _{\lambda '},\psi _\lambda \rangle | \lesssim 2^{-||\lambda |-|\lambda '|| \sigma }(1+d(\lambda ,\lambda '))^{-\beta }, {\tag{2.28}} \end{equation}$$

with $\sigma >d/2$ and $\beta >d$ and

$$\begin{equation} d(\lambda ,\lambda '):=2^{\min (|\lambda |, |\lambda '|)}\mathop{dist}(\operatorname {supp}(\psi _\lambda ),\operatorname {supp}(\psi _{\lambda '})). {\tag{2.29}} \end{equation}$$

The validity of 2.28 has been established in numerous settings (see e.g., 23114042). Decay estimates of the form 2.28 were initially introduced in 37 in the context of Littlewood-Paley analysis. The constant $\sigma$ depends on the smoothness of the wavelets, whereas $\beta$ is related to the approximation order of the dual multiresolution (resp. the vanishing moments of the wavelets) and the order of the operator $A$. Estimates of the type 2.28 are known to hold for a wide range of cases, including classical pseudo-differential operators and Calderón-Zygmund operators (see e.g., 2641). In particular, the single and double layer potential operators fall into this category. We refer the reader to 23 for a full discussion of settings where 2.28 is valid.

We introduce the class $\mathcal{A}_{\sigma ,\beta }$ of all matrices $\mathbf{B}=(b_{\lambda ,\lambda ^{\prime }})_{\lambda ,\lambda ^{\prime }\in \nabla }$ which satisfy

$$\begin{equation} |b_{\lambda ,\lambda ^{\prime }}|\le c_\mathbf{B}2^{-||\lambda |-|\lambda ^{\prime }||\sigma } (1+d(\lambda ,\lambda ^{\prime }))^{-\beta }, {\tag{2.30}} \end{equation}$$

with $d(\lambda ,\lambda ^{\prime })$ defined by 2.29. We say that a matrix $\mathbf{B}$ is quasi-sparse if it is in the class $\mathcal{A}_{\sigma ,\beta }$ for some $\sigma >d/2$ and $\beta >d$. Properties of quasi-sparse matrices will be discussed in Section 3.

(A3): We assume that, for some $\sigma >d/2$, $\beta >d$, the matrix $\mathbf{A}$ of 2.17 is in the class $\mathcal{A}_{\sigma ,\beta }$.

Let us note that in the case $H=H^t$ discussed earlier, we obtain 2.30 from 2.28 because $\mathbf{D}=(2^{-t|\lambda |}\delta _{\lambda ,\lambda '})_{\lambda ,\lambda ' \in \nabla }$.

2.5. Wavelet Galerkin methods

A wavelet based Galerkin method for solving 2.5 takes the following form. We choose a finite set $\Lambda$ of wavelet indices and use the space $S_\Lambda :=\mathrm{span}\{\psi _\lambda : \lambda \in \Lambda \}$ as our trial and analysis space. The approximate Galerkin solution $u_\Lambda$ from $S_\Lambda$ is defined by the conditions