On Multiresolution Methods

in Numerical Analysis

Gregory Beylkin

Abstract. As a way to emphasize several distinct features of the multiresolution methods based on wavelets, we describe connections between the multiresolution LU decomposition, multigrid and multiresolution reduction/homogenization for self-adjoint, strictly elliptic operators. We point out that the multiresolution LU decomposition resembles a direct multigrid method (without W-cycles) and that the algorithm scales properly in higher dimensions.

Also, the exponential of these operators is sparse where sparsity is defined as that for a finite but arbitrary precision. We describe time evolution schemes for advection-diffusion equations, in particular the Navier-Stokes equation, based on using sparse operator-valued coefficients. We point out a significant improvement in the stability of such schemes.

1991 Mathematics Subject Classification: 65M55, 65M99, 65F05, 65F50, 65R20, 35J, 76D05 Keywords and Phrases: multigrid methods, fast multipole method, wavelet bases, multiresolution analysis, multiresolution LU decomposition, time evolution schemes, exponential of operators, advectiondiffusion equations

1 Introduction

Multiresolution methods have a fairly long history in numerical analysis, going back to the introduction of multigrid methods [10], [18] and even earlier [22]. A renewed interest in multiresolution methods was generated recently by the development of wavelet bases and other bases with controlled time-frequency localization [23], [20], [13], [19], [12], [2], [1], etc.. The introduction of these new tools allows us to relate numerical analysis with harmonic analysis and signal processing by the fundamental need of an efficient representation of operators and functions.

It is useful to compare the wavelet approach with the multigrid method (MG) and the Fast Multipole Method (FMM). For most problems the wavelet approach, FMM, and MG provide the same asymptotic complexity. The differences are typically in the “constants” of the complexity estimates. These differences will, most likely, diminish in the future.

A typical MG is a fast iterative solver based on a hierarchical subdivision. Hierarchical subdivision is also used in FMM which was initially proposed for computing potential interactions [21], [17]. This algorithm requires order N operations to compute all the sums

p_{j} = ∑_{i≠j} \frac{q_{i}q_{j}}{|x_{i} − x_{j}|}, where x_{i} ∈ R^{3} i, j = 1, …, N,  (1)

and the number of operations is independent of the configuration of charges. In the FMM, the reduction of the complexity of computing the sums in (1) from order N^{2} to −N log ϵ, where ϵ is the desired accuracy, is achieved by approximating the far field effect of a cloud of charges located in a box by the effect of a single multipole at the center of the box.

Although both MG and FMM have been extended well beyond their original applications, neither of these methods use the notion of bases in their development and, specifically, orthonormal bases^{1}. On the conceptual level using bases makes it easier to consider efficient representations of functions and operators that handle smooth, oscillatory, and scaling behavior.

In particular, to emphasize several distinct features of the wavelet approach, we consider two topics. First, we describe connections between the multiresolution LU decomposition, MG, and multiresolution reduction/homogenization for selfadjoint, strictly elliptic operators. Second, we describe the effects of computing the exponential of such operators on numerical properties of time evolution schemes for advection-diffusion equations.

The essence of the first topic is that multiresolution LU decomposition (the usual LU decomposition interlaced with projections) is equivalent to the direct MG, i.e., a MG without W-cycles. The reason for the absence of W-cycles is that on every scale we construct equations for the orthogonal projection of the true solution. Once these equations are solved, there is no need to return to a coarser scale to correct the solution (which is the role of W-cycles in MG). Moreover, equations obtained in this manner on coarser scales are of interest by themselves, since they can be interpreted as “homogenized” or reduced equations, leading to (numerical) multiresolution reduction and homogenization.

The essence of the second topic is that we can drastically improve properties of time evolution schemes for advection-diffusion equations by using the exponential of operators. As it turns out, for self-adjoint, strictly elliptic operators L the exponential exp (−tL) is sparse in wavelet bases (for a finite but arbitrary precision) for all t ≥ 0. This observation makes the construction of exp (−tL) feasible in two and three spatial dimensions. Given a proper choice of basis and several additional algorithms, we are led to adaptive numerical schemes for the solution of advection-diffusion equations [8].

1We note that the representation of functions via their values and via coefficients in an expansion are closely related. In fact if one uses interpolating bases functions then there is a way to simplify this relation (see [3]).

Documenta Mathematica · Extra Volume ICM 1998 · III · 481–490

2 Multiresolution Direct Solvers

Direct solvers are not used for problems in multiple dimensions since the standard LU decomposition will fill most of the matrix and, thus, render the method inefficient. This is even without considering additional difficulties due to the high condition numbers typical in these problems. It turns out that both difficulties can be overcome for self-adjoint, strictly elliptic operators by using wavelet bases and multiresolution LU decomposition [7], [16].

As usual, we consider multiresolution analysis (MRA), a chain of subspaces

… ⊂ V_{2} ⊂ V_{1} ⊂ V_{0} ⊂ V_{−1} ⊂ V_{−2} ⊂ …

such that

⋂_{j} V_{j} = \{0\} and \overline{⋃_{j} V_{j}} = L^{2}(R^{d}).

Let the subspace V_{j} be spanned by an orthonormal basis formed by the tensor product of scaling functions {ϕ_{k}^{j}(·) = 2^{−j/2}ϕ(2^{−j} · −k)}_{k∈Z}, where ϕ satisfies the two-scale difference equation (see e.g. [13] for details). Let us denote by W_{j} the orthogonal complement of V_{j} in V_{j−1}, V_{j−1} = Vj ⊕ Wj. We use P_{j} and Q_{j} to denote the projection operators onto V_{j} and W_{j}. If x ∈ V_{j}, we write s_{x} = P_{j+1}x and d_{x} = Q_{j+1}x, where s_{x} ∈ V_{j+1} and d_{x} ∈ W_{j+1}.

Given a bounded linear operator S on L^{2}(R^{d}), let us consider its projection S_{j} on V_{j}, S_{j} = P_{j}SP_{j} and represent the operator S_{j} as a (possibly infinite) matrix in that basis. With a slight abuse of notation, we will use the same symbol S_{j} to represent both the operator and its matrix. Since V_{j} = V_{j+1} ⊕ W_{j+1}, we may also write S_{j} : V_{j} → V_{j} in a block form

S_{j} = ( \begin{matrix}A_{S_{j}} & B_{S_{j}} \\ C_{S_{j}} & T_{S_{j}} \end{matrix} ) : V_{j+1} ⊕ W_{j+1} → V_{j+1} ⊕ W_{j+1},  (2)

where A_{S_{j}} = Q_{j+1}S_{j}Q_{j+1}, B_{S_{j}} = Q_{j+1}S_{j}P_{j+1}, C_{S_{j}} = P_{j+1}S_{j}Q_{j+1}, and T_{S_{j}} = S_{j+1} = P_{j+1}S_{j}P_{j+1}. Each of the operators may be considered as a matrix

and in the matrix form the transition from S_{j} in (2) to ( \begin{matrix}A_{S_{j}} & B_{S_{j}} \\ C_{S_{j}} & T_{S_{j}} \end{matrix} ) requires

application of the wavelet transform. We refer to A_{S_{j}}, B_{S_{j}}, C_{S_{j}} and T_{S_{j}} as the A, B, C, and T blocks of S_{j}.

Consider a bounded linear operator S_{j} : V_{j} → V_{j} and a linear equation

S_{j}x = f,  (3)

which we may write as

( \begin{matrix}A_{S_{j}} & B_{S_{j}} \\ C_{S_{j}} & T_{S_{j}} \end{matrix} ) ( \begin{matrix}d_{x} \\ s_{x} \end{matrix} ) = ( \begin{matrix}d_{f} \\ s_{f} \end{matrix} ) .  (4)

Formally eliminating d_{x} from (4) by substituting d_{x} = A_{S_{j}}^{−1}(d_{f} − B_{S_{j}}s_{x}) (Gaussian elimination) yields

(T_{S_{j}} − C_{S_{j}}A_{S_{j}}^{−1}B_{S_{j}})s_{x} = s_{f} − C_{S_{j}}A_{S_{j}}^{−1}d_{f}.  (5)

We call (5) the reduced equation, and the operator

R_{S_{j}} = T_{S_{j}} − C_{S_{j}}A_{S_{j}}^{−1}B_{S_{j}}  (6)

the one-step reduction of the operator S_{j}. The right-hand side of (6) is also known

the Schur complement of the block-matrix ( \begin{matrix}A_{S_{j}} & B_{S_{j}} \\ C_{S_{j}} & T_{S_{j}} \end{matrix} ).

Note that the solution s_{x} of the reduced equation is exactly P_{j+1}x, the projection of the solution of the original equation in V_{j+1}. The solution of the reduced equation is the same on the subspace V_{j+1} as the solution of the original equation (3). Once we have obtained the reduced equation, it may be reduced again to produce an equation on V_{j+2}. Likewise, we may reduce n times to produce an equation on V_{j+n} the solution of which is the projection of the solution of (3) on V_{j+n}. We note that in the finite-dimensional case, the reduced equation (5) has 1/2^{d} as many unknowns as the original equation (3). Reduction, therefore, preserves the coarse-scale behavior of solutions while reducing the number of unknowns.

The critical questions are: (i) can we control the sparsity (for any finite but arbitrary precision) of the matrix C_{S_{j}}A_{S_{j}}^{−1}B_{S_{j}}? and, (ii) can we repeat the reduction step for R_{S_{j}}? In MG literature the Schur complement appears in a number of papers but these questions were not answered. In [7] and [16] these questions were answered affirmatively for a finite number of reduction steps. The key property that makes this affirmative answer possible is the vanishing moments property of the basis functions.

The sparsity (for any finite but arbitrary precision) of the multiresolution LU factorization does not depend on dimension. This is in a sharp contrast with the usual practice, where LU factorization is not recommended as an efficient approach in problems of dimension two or higher. For example, if we consider the Poisson equation, then LU decomposition is not considered as a practical option since the fill-ins will yield dense LU factors.

A close examination of the algorithm in [16] reveals a striking resemblance of the multiresolution LU decomposition coupled with the multiresolution forward and backward substitution to a MG technique. The important difference, however, is that there are no W-cycles.

As described above, reduction is an algebraic procedure carried out on matrices over a finite number of scales. It relies on the explicit hierarchy of scales provided by the MRA to algebraically eliminate the fine-scale variables, leaving only the coarse-scale variables and can be cast as a multiresolution reduction procedure for the corresponding ODEs and PDEs [11]. The classical homogenization of partial differential equations is the process of finding “effective” coefficients. In classical homogenization, the fine scale is associated with a small parameter, and the limit is considered as this small parameter goes to zero. In dimension one a connection has been established [15],[14] between multiresolution reduction and classical homogenization (see e.g. [4]). It is important to point out that reduction approximately preserves small eigenvalues of elliptic operators, and the accuracy of this approximation depends on the order of the wavelets [7].

3 Sparsity of Exponential Operators

If L is a self-adjoint, strictly elliptic operator then the operator e^{Lt} is sparse in wavelet bases (for a finite but arbitrary precision) for all t ≥ 0. This observation has a significant effect on the methods for solving PDEs.

Let us consider a class of advection-diffusion equations of the form

u_{t} = Lu + N (u),  x ∈ Ω ⊂ R^{d},  (7)

where u = u(x, t), x ∈ R^{d}, d = 1, 2, 3 and t ∈ [0, T] with the initial conditions,

u(x, 0) = u_{0}(x),  x ∈ Ω,  (8)

and the linear boundary conditions

Bu(x, t) = 0,  x ∈ ∂Ω, t ∈ [0, T].  (9)

In (7) L represents the linear and N (·) the nonlinear terms of the equation, respectively.

Using the semigroup approach we rewrite the partial differential equation (7) as a nonlinear integral equation in time,

u(x, t) = e^{(t−t_{0})L}u_{0}(x) + ∫_{t_{0}}^{t} e^{(t−τ)L}N (u(x, τ)) dτ,  (10)

and describe a new class of time-evolution schemes based on its discretization. A distinctive feature of these new schemes is exact evaluation of the contribution of the linear part. Namely, if the nonlinear part is zero, then the scheme reduces to the evaluation of the exponential function of the operator (or matrix) L representing the linear part.

We note that the incompressible Navier-Stokes equations can be written in the form (7). Let us start with the usual form of the Navier-Stokes equations for x ∈ Ω ⊂ R^{3},

u_{t} = νΔu − (u_{1}∂_{1} + u_{2}∂_{2} + u_{3}∂_{3}) u − ∇p,  (11)
∂_{1}u_{1} + ∂_{2}u_{2} + ∂_{3}u_{3} = 0,  (12)
u(x, 0) = u_{0},  (13)

where p denotes the pressure and u = ( \begin{matrix}u_{1} \\ u_{2} \\ u_{3} \end{matrix} ), x = ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ) and ∂_{k} = \frac{∂}{∂x_{k}}. In

addition, we impose the boundary condition

u(x, t) = 0 x ∈ ∂Ω, t ∈ [0, T],  (14)

Let us introduce the Riesz transforms which are defined in the Fourier domain as

ˆ (R_{j}f)(ξ) = \frac{ξ_{j}}{|ξ|}f̂(ξ), j = 1, 2, 3,  (15)

where f̂ denotes the Fourier transform of the function f. It is not difficult to show that the projection operator on the divergence free functions (the Leray projection) may be written with the help of the Riesz transforms,

P = ( \begin{matrix}I & 0 & 0 \\ 0 & I & 0 \\ 0 & 0 & I \end{matrix} ) − ( \begin{matrix}R_{1}^{2} & R_{1}R_{2} & R_{1}R_{3} \\ R_{2}R_{1} & R_{2}^{2} & R_{2}R_{3} \\ R_{3}R_{1} & R_{3}R_{2} & R_{3}^{2} \end{matrix} ) .  (16)

Applying the divergence operator to (11), we obtain −Δp = ∑_{k,l=1}^{3} ∂_{k}∂_{l}u_{k}u_{l} and an expression for pressure in terms of the Riesz transforms, p = − ∑_{k,l=1}^{3} R_{k}R_{l}(u_{k}u_{l}). Substituting the expression for the pressure into (11) and taking into consideration that the Riesz transforms commute with derivatives and, moreover, R_{k}∂_{l} = R_{l}∂_{k}, we obtain

u_{t} = νΔu − P( ∑_{m=1}^{3} u_{m}∂_{m}u),  (17)

instead of (11) and (12). Equations (17) are now in the form (7), where L = νΔ and N (u) = −P(∑_{m=1}^{3} u_{m}∂_{m}u). The transformation from (11) and (12) to (17) is well known and appears in a variety of forms in the literature. Here we followed a derivation presented by Yves Meyer at Summer School at Luminy in 1997.

The apparent problem with (17) for use in numerical computations is that the Riesz transforms are integral operators (which makes (17) into an integrodifferential equation). Let us point out that the presence of the Riesz transforms does not create serious difficulties if we represent operators R_{j}, j = 1, 2, 3 in a wavelet basis with a sufficient number of vanishing moments (for a given accuracy). The reason is that these operators are nearly local on wavelets, and thus, have a sparse representation. This approximate locality follows directly from the vanishing moments property. Vanishing moments imply that the Fourier transform of the wavelet and its several first derivatives vanish at zero, and therefore, the discontinuity of the symbol of the Riesz transform at zero has almost no effect. The precise statements about such operators can be found in [6] and [5].

Finally, in rewriting (17) as u_{t} = Lu + N (u), we incorporate the boundary conditions into the operator L. For example, u = L^{−1}v means that u solves Lu = v with the boundary conditions Bu = 0. Similarly, u(x, t) = e^{Lt}u_{0}(x) means that u solves u_{t} = Lu, u(x, 0) = u_{0}(x) and Bu(x, t) = 0.

Computing and applying the exponential or other functions of operators in the usual manner typically requires evaluating dense matrices and is highly inefficient unless there is a fast transform that diagonalizes the operator. For example, if L is a circulant matrix, then computing functions of operators can be accomplished using the FFT. It is clear that in this case the need of the FFT for diagonalization prevents one from extending this approach to the case of variable coefficients.

In the wavelet system of coordinates computing the exponential of self-adjoint, strictly elliptic operators always results in sparse matrices, and therefore, using the exponential of operators for numerical purposes is an efficient option [8].

Further development of the approach of [8] can be found in [9], where issues of stability of time-discretization schemes with exact treatment of the linear part (ELP) schemes are considered. The ELP schemes are shown to have distinctly different stability properties as compared with the usual implicit-explicit schemes. The stability properties of traditional time-discretization schemes for advectiondiffusion equations are controlled by the linear term and, typically, these equations require implicit treatment to avoid choosing an unreasonably small time step. As it is shown in [9], using an explicit ELP scheme, it is possible to achieve stability usually associated with implicit predictor-corrector schemes.

If an implicit ELP scheme is used, as it is done in [8], an approximation is used only for the nonlinear term. This changes the behavior of the corrector step of implicit schemes. The corrector step iterations of usual implicit schemes for advection-diffusion equations involve either both linear and nonlinear terms or only the linear term. Due to the high condition number of the matrix representing the linear (diffusion) term, convergence of the fixed point iteration requires a very small time step, making the fixed point iteration impractical. Implicit ELP schemes do not involve the linear term and, typically, the fixed point iteration is sufficient as in [8].

We would like to note, that (10) in effect reduces the problem to an ODE–type setup, and for that reason, a variety of methods can be used for its solution. We present operator valued coefficients of multistep ELP schemes and our main point is that these coefficients can be represented by sparse matrices and applied in an efficient manner.

Let us consider the function u(x, t) at the discrete moments of time t_{n} = t_{0} + nΔt, where Δt is the time step so that u_{n} ≡ u(x, t_{n}) and N_{n} ≡ N (u(x, t_{n})). Discretizing (10) yields

u_{n+1} = e^{lLΔt}u_{n+1−l} + Δt (γN_{n+1} + ∑_{m=0}^{M−1} β_{m}N_{n−m}) ,  (18)

where M +1 is the number of time levels involved in the discretization, and l ≤ M. The expression in parenthesis in (18) may be viewed as the numerical quadrature for the integral in (10). The coefficients γ and β_{m} are functions of LΔt. In what follows we restrict our considerations to the case l = 1. We observe that the algorithm is explicit if γ = 0 and it is implicit otherwise. Typically, for a given M, the order of accuracy is M for an explicit scheme and M + 1 for an implicit scheme due to one more degree of freedom, γ.

For l = 1 we provide Tables 1 and 2 for M = 1, 2, 3 with expressions for the coefficients of the implicit (γ ≠ 0) and the explicit (γ = 0) schemes in terms of Q_{k} = Q_{k}(LΔt), where

Q_{k}(LΔt) = \frac{e^{LΔt} − E_{k}(LΔt)}{(LΔt)^{k}},  (19)

and

E_{k}(LΔt)) = ∑_{l=0}^{k−1} \frac{(LΔt)^{l}}{l!}  (20)

Documenta Mathematica · Extra Volume ICM 1998 · III · 481–490

M γ  β_{0}  β_{1}  β_{2}
1 Q_{2}  Q_{1} − Q_{2}  0  0
2 \frac{1}{2}Q_{2} + Q_{3}  Q_{1} − 2Q_{3}  Q_{3} − \frac{1}{2}Q_{2}  0
3 \frac{1}{3}Q_{2} + Q_{3} + Q_{4} Q_{1} + \frac{1}{2}Q_{2} − 2Q_{3} − 3Q_{4} Q_{3} − Q_{2} + 3Q_{4} \frac{1}{6}Q_{2} − Q_{4}

Table 1: Coefficients of implicit ELP schemes for l = 1, where Q_{k} = Q_{k}(LΔt)).

M β_{0}  β_{1}  β_{2}
1 Q_{1}  0  0
2 Q_{1} + Q_{2}  −Q_{2}  0
3 Q_{1} + 3Q_{2}/2 + Q_{3} −2(Q_{2} + Q_{3}) Q_{2}/2 + Q_{3}

Table 2: Coefficients of explicit ELP schemes for l = 1, where Q_{k} = Q_{k}(LΔt).

In Tables 1 and 2 we have presented examples of the so-called “bare” coefficients. Modified coefficients [8] differ in high order terms: these terms do not affect the order of accuracy but do affect the stability properties. Modified coefficients depend on a particular form of the nonlinear term.

Let us describe a method to compute operators Q_{0}, Q_{1}, Q_{2}, …. without computing (LΔt)^{−1}. In computing the exponential, Q_{0}, we use the scaling and squaring method which is based on the identity

Q_{0}(2x) = (Q_{0}(x))^{2} .  (21)

First we compute Q_{0}(LΔt2^{−l}) for some l chosen so that the largest singular value of LΔt2^{−l} is less than one. This computation is performed using the Taylor expansion. Using (21), the resulting matrix is then squared l times to obtain the final answer. In all of these computations it is necessary (and possible) to remove small matrix elements to maintain sparsity, and at the same time, maintain a predetermined accuracy.

A similar algorithm may be used for computing Q_{j}(LΔt), j = 1, 2, … for any finite j. Let us illustrate this approach by considering j = 1, 2. It is easy to verify that

Q_{1}(2x) = \frac{1}{2} (Q_{0}(x)Q_{1}(x) + Q_{1}(x)) ,  (22)
Q_{2}(2x) = \frac{1}{4} (Q_{1}(x)Q_{1}(x) + 2Q_{2}(x)) .  (23)

Thus, a modified scaling and squaring method for computing operator-valued quadrature coefficients for ELP schemes starts by the computation of Q_{0}(LΔt2^{−l}), Q_{1}(LΔt2^{−l}) and Q_{2}(LΔt2^{−l}) for some l selected so that the largest singular value of all three operators is less than one. For these evaluations we use the Taylor expansion. We then proceed by using identities in (21), (22) and (23) l times to compute the operators for the required value of the argument.

As an example consider Burgers’ equation

u_{t} + uu_{x} = νu_{xx},  0 ≤ x ≤ 1, t ≥ 0,  (24)

for ν > 0, together with an initial condition,

u(x, 0) = u_{0}(x),  0 ≤ x ≤ 1,  (25)

and periodic boundary conditions u(0, t) = u(1, t). Burgers’ equation is the simplest example of a nonlinear partial differential equation incorporating both linear diffusion and nonlinear advection. In [8] a spatially adaptive approach is used to compute solutions of Burgers’ equation via

u_{n+1} = Q_{0}(LΔt)u_{n} − \frac{Δt}{2}Q_{1}(LΔt) [u_{n}∂_{x}u_{n+1} + u_{n+1}∂_{x}u_{n}] .  (26)

We refer to [9] for the analysis of stability of ELP schemes.

4 Conclusions

The wavelet based algorithms described above are quite efficient in dimension one. Although algorithms described above scale properly with size in all dimensions, establishing ways of reducing operation counts remains an important task in dimensions two and three. This is an area of the ongoing research and the progress will be reported elsewhere.

References

[1] B. Alpert. A Class of Bases in l^{2} for the Sparse Representation of Integral Operators. SIAM J. Math. Anal, 24(1):246–262, 1993.

[2] B. Alpert, G. Beylkin, R. R. Coifman, and V. Rokhlin. Wavelet-like bases for the fast solution of second-kind integral equations. SIAM Journal of Scientific and Statistical Computing, 14(1):159–174, 1993.

[3] B. Alpert, G. Beylkin, D. Gines, and L. Vozovoi. Toward adaptive solution of partial differential equations in multiwavelet bases. 1998. in progress.

[4] A. Bensoussan, J.L. Lions, and G. Papanicolaou. Asymptotic Analysis for Periodic Structures. North-Holland Pub. Co., New York, 1978.

[5] G. Beylkin. On the representation of operators in bases of compactly supported wavelets. SIAM J. Numer. Anal., 29(6):1716–1740, 1992.

[6] G. Beylkin, R. R. Coifman, and V. Rokhlin. Fast wavelet transforms and numerical algorithms I. Comm. Pure and Appl. Math., 44:141–183, 1991.

[7] G. Beylkin and N. Coult. A multiresolution strategy for reduction of elliptic PDE’s and eigenvalue problems. Applied and Computational Harmonic Analysis, 5:129–155, 1998.

[8] G. Beylkin and J.M. Keiser. On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases. J. Comp. Phys., 132:233–259, 1997.

[9] G. Beylkin, J.M. Keiser, and L. Vozovoi. A new class of stable time discretization schemes for the solution of nonlinear PDEs. PAM Report 347, 1998. submitted to JCP.

[10] A. Brandt. Multi-level adaptive solutions to boundary value problems. Math. Comp., 31:333–390, 1977.

[11] M. E. Brewster and G. Beylkin. A Multiresolution Strategy for Numerical Homogenization. ACHA, 2:327–349, 1995.

[12] R. R. Coifman and Y. Meyer. Nouvelles bases orthogonales. C.R. Acad. Sci., Paris, 1990.

[13] I. Daubechies. Orthonormal bases of compactly supported wavelets. Comm. Pure and Appl. Math., 41:909–996, 1988.

[14] M. Dorobantu and B. Engquist. Wavelet-based numerical homogenization. SIAM J. Numer. Anal., 35(2):540–559, 1998.

[15] A. C. Gilbert. A comparison of multiresolution and classical one-dimensional homogenization schemes. to appear in Appl. and Comp. Harmonic Analysis.

[16] D.L. Gines, G. Beylkin, and J. Dunn. LU factorization of non-standard forms and direct multiresolution solvers. Applied and Computational Harmonic Analysis, 5:156–201, 1998.

[17] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J. Comp. Phys., 73(1):325–348, 1987.

[18] W. Hackbusch. On multi-grid method applied to difference equations. Computing, 20:291–306, 1978.

[19] H. S. Malvar. Lapped Transforms for Efficient Transform/Subband Coding. IEEE Trans. Acoust., Speech, Signal Processing, 38(6):969–978, 1990.

[20] Y. Meyer. Principe d’incertitude, bases hilbertiennes et algébres d’opérateurs. In Séminaire Bourbaki,, page 662. Société Mathématique de France, 1985-86. Astérisque.

[21] V. Rokhlin. Rapid solution of integral equations of classical potential theory. J. Comp. Phys., 60(2), 1985.

[22] R.P.Fedorenko. The speed of convergence of one iterative process. USSR Comp. Math. and Math. Physics, 4:227–235, 1964.

[23] J. O. Stromberg. A Modified Franklin System and Higher-Order Spline Systems on R^{n} as Unconditional Bases for Hardy Spaces. In Conference in harmonic analysis in honor of Antoni Zygmund, Wadworth math. series, pages 475–493, 1983.

Gregory Beylkin

Dept. of Applied Mathematics University of Colorado at Boulder Boulder, CO 80309-0526, USA