A Wavelet Multigrid Method Using Symmetric Biorthogonal Wavelets

In [1], the author introduced a wavelet multigrid method that used the wavelet transform to define the coarse grid, interpolation, and restriction operators for the multigrid method. In this paper, we modify the method by using symmetric biorthogonal wavelet transforms to define the requisite operators. Numerical examples are presented to demonstrate the effectiveness of the modified wavelet multigrid method for diffusion problems with highly oscillatory coefficients, as well as for advection-diffusion equations in which the advection is moderately dominant.


Introduction
It is well known that the multigrid method is very useful in increasing the efficiency of iterative methods used to solve systems of algebraic equations approximating partial differential equations.However, when confronted by certain problems, for example diffusion problems with discontinuous or highly oscillatory coefficients, the standard multigrid procedure converges slowly, with a rate dependent on the initial mesh size, or may even break down.
In [2], the wavelet transform, which uses both highand low-pass filter operators, is used to derive a new approach for the two-level multigrid method, under the assumption that the matrix on the fine grid is symmetric.Some one-dimensional examples are also examined in that paper.In [1,3], the author extended the results of this approach to two dimensions and to multiple-level multigrid, dropping the assumption of a symmetric fine grid operator.This approach was considered for several reasons.First, in [4], the authors investigate similarities between multiresolution analyses and multigrid methods, which similarities motivate investigation into waveletbased multigrid methods.In addition, in [5], for example, it is shown that a wavelet coarse grid operator defined by a Schur complement provides a good approximation to the homogenized coarse grid operator, and, as stated in [1,3], homogenization has been used to improve convergence of multigrid methods for diffusion problems with periodic coefficients (e.g., [6][7][8][9]) because the homogenized operator provides a very good approximation of the important properties (e.g., eigenvalues and eigenfunctions) of the original fine grid operator.Also a wavelet coarse grid operator defined by a Schur complement has a natural connection to the interpolation and restriction operators.Furthermore, wavelets can be applied to problems with periodic as well as non-periodic coefficients.Finally, the application of wavelet operators to vectors and matrices maintains the properties of the original problem.
Recently, symmetric biorthogonal wavelets have received attention for use in image compression (see, e.g., [10,11]), based on the lifting scheme developed by Sweldens in [12].In fact, symmetric biorthogonal wavelets are the basis for lossless compression in the JPEG 2000 standard (see, e.g., [13]).In addition, using biorthogonal wavelets to define multigrid methods has been explored in other papers.For example, in [14], the authors use biorthogonal wavelets adapted to a differential operator with constant coefficients, from which they develop a system of linear equations, and then use a biorthogonal two-grid method to solve.In [15], the authors consider solving ill-conditioned symmetric Toeplitz systems by using a two-grid method in which the interpolation and restriction operators are defined using the Cohen, Daubechies, and Feauveau (CDF) 9/7 symmetric biorthogonal wavelets.In [16], the authors use biorthogonal wavelets as preconditioners for an algebraic multigrid method which they use to solve various partial differential equations.
In this paper, the wavelet multigrid method introduced in [1] is modified by using symmetric biorthogonal wavelet transforms to define the requisite operators.The organization of this paper is as follows.In Section 2, some background on multigrid methods and wavelets is given.Section 3 discusses the wavelet multigrid method developed from the application of the biorthogonal wavelet transform to a general second order partial differential equation in two dimensions.In Section 4, numerical considerations are discussed.Section 5 presents some numerical results of applying the modified wavelet multigrid method using symmetric biorthogonal wavelets to two diffusion problems with highly oscillatory coefficients and to three advection-diffusion problems with moderately dominant advection.The rapid convergence, relatively independent of the initial mesh size, of the modified wavelet multigrid method using biorthogonal wavelets is demonstrated for these problems.For all numerical results in this paper, the V-cycle multigrid method is used with one iteration of the smoother for the coarsening and the correcting phases.

Multigrid Methods
The problem we are concerned with solving is the system of linear equations where A and arise from discretization of a (partial) differential equation on some grid , where is the step size.
We briefly describe the V-cycle method used in this paper.Given some interpolation operator, 2 , h h I where the superscript refers to the fine grid and the subscript refers to the coarse grid, and a restriction operator, 2 , h h I we can define a multigrid method recursively.For the two-level V-cycle method, we do the following.
Step 1: Relax a few (usually one or two) steps on the fine grid to get an initial guess .
Step 2: Compute the residual and restrict the residual to the coarse grid .
on the coarse grid.
Step 4: Set and again relax a few (usually one or two) steps on the fine grid.Based on this two-level method, the V-cycle multigrid scheme is defined recursively.Some good references for multigrid methods are [17][18][19].
One type of multigrid scheme is algebraic multigrid, which only uses the structure of the matrix in the problem to determine the coarsening process (choice of coarse grid and definition of interpolation/restriction operators).This process is performed in order to ensure that the range of interpolation approximates the errors not sufficiently reduced via relaxation.For a more detailed description of algebraic multigrid methods, see, e.g., [19][20][21][22].Note that in [21] the relationship between algebraic multigrid and Schur complements is discussed.Algebraic multigrid methods are of particular interest, in that they are the nearest methods to the approach used for the wavelet multigrid method.
It is expected that the multigrid method should converge at a rate independent of the fine mesh size.However, for certain problems, including elliptic problems with highly oscillatory coefficients, such convergence does not occur for the standard multigrid method.One difficulty is that the small eigenvalues of A are not necessarily associated with smooth eigenfunctions, a key assumption for the standard multigrid method.For such problems, it is not as simple to approximate the smooth eigenfunctions on the coarse grids.New methods for restriction and interpolation, or for treating the entire problem, must be found.One such approach is the wavelet multigrid method discussed in [1,3] and the modified wavelet multigrid method discussed in this paper.

Wavelets and Biorthogonal Wavelets
For background, a brief description of wavelets, and then of biorthogonal wavelets, follows.For more details, the reader is referred to [10,[23][24][25].

Wavelets
Wavelets basically separate data (or functions or operators) into different frequency components and analyze them by scaling.The wavelets can be chosen to form a complete orthonormal basis of   2 L  .Due to the scaling of the wavelet functions, they have time-or spacewidths that are related to their frequency: at high frequencies, they are narrow, and at low frequencies, they are broader.Therefore, they provide good localization of functions in both the frequency domain and physical space, and representation by wavelets seems natural to apply to the analysis of fine and coarse scales.
A multiresolution analysis (MRA) consists of a sequence of closed subspaces L  , the scaling spaces, that satisfy certain conditions.For every , j   j W is defined as the orthogonal complement of  G to be the operators that transform the basis of the space j V to the bases of the spaces : .So, the application of the wavelet transform requires only operations.In general, application of the wavelet transform requires operations, assuming a finite number of coefficients for the low-and high-frequency operators.
In two dimensions, the tensor product of one-dimensional multiresolution analyses is used.
Then, analogous to the one dimensional case, define the operators j H and j G so that H and j G have the same properties as their onedimensional analogues.As in the one-dimensional case, the wavelet transform, 1 :

Biorthogonal Wavelets
A biorthogonal multiresolution analysis consists of two dual multiresolution analyses.In other words, there are two dual sequences of closed subspaces of and , both of which satisfy the conditions of a multiresolution analysis.For every j   , 1 j V  can be written as the direct sum of j V and j W , and  can be written as the direct sum of j V and  The wavelet transform, :


and the dual transform : where we use  to denote the direct sum for ease of notation.Note that j  is orthogonal to j   due to the properties of j H , j G , j H  and j G  .Also, note that in terms of filter theory, j H and j H  are low-pass filters, and j G and j G  are high-pass filters.The discrete biorthogonal wavelet transforms are also computationally efficient, and application of these operators requires operations, assuming a finite number of coefficients for the low-and high-frequency operators.
In two dimensions, we use the tensor product in a similar way as for the standard multiresolution analysis.So, j V is defined by and j V is defined by L  .j W and j W are defined so that their basis functions are the tensor products of the one-dimensional wavelets and scaling functions.For each

V
can be written as the direct sum of j V and j W , and can be written as the direct sum of j V  and j W  .So, we have V , where we again use   to denote the direct sum for ease of notation.
Then, analogous to the one-dimensional case, define the operators j H , j G , j H  and j G  so that: and , and and .
where , V  and , , and have the same properties as their one-dimensional analogues.As in the one-dimensional case, the wavelet transform and the dual transform The main idea for the remainder of the method follows in a similar manner as the wavelet multigrid method using orthogonal wavelet transforms.First, ˆ is computed, and the resulting matrix, denoted by j L , is partitioned to obtain  , are defined by and .
j j j j j j j j j j j j j j j j j j

The Modified Wavelet Multigrid Method Using Biorthogonal Wavelets
Then, the block UDL decomposition of ˆj where U is block upper triangular with unit diagonal, D is block diagonal, and L is block lower triangular with unit diagonal, is computed and is used to find j L  as follows.

The block UDL decomposition of ˆj
L is determined to be The two-dimensional wavelet multigrid method was introduced in [3] and [1].Both of these works assumed orthogonal wavelet operators.Here, we describe the modified wavelet multigrid method using biorthogonal wavelet transforms.For the work in this paper, the biorthogonal wavelet transform is formed using symmetric biorthogonal wavelets.
The inverse of this factorization of ˆj L is then computed, which after multiplication of the factors gives where j L represents the operator obtained by discretizing a two-dimensional partial differential equation on the where we observe that Motivated by the work in [1], we denote and as the interpolation and restriction operators, respectively.
Plugging the interpolation and restriction operators defined in ( 9) and ( 10) into (8) gives, In multigrid, the error correction on the coarse grid is sought, i.e., the equation being solved is (11) with u replaced by the error, e, and f replaced by the residual, r: is small, i.e., r is almost in the error can be approximated by The above assumption is good for most of the classical iterative methods, like Jacobi and Gauss-Seidel.Therefore, the coarse grid operator is defined by which we again note is the Schur complement of j D in ˆ. j L Note: This coarse grid operator is the same as the operator obtained by solving for L u : The main numerical issue is fill-in in the matrices as the The above procedure may be repeatedly applied until the desired coarseness is reached.Notice that this method is applicable to problems which are not symmetric.However, regardless of the discretization scheme used, the operator matrix must be a square matrix whose row size is a multiple of four.

Numerical Considerations
multigrid method involves increasing levels in the Vcycle.Thresholding is used to reduce fill-in in the computation of 1 , , , and D is is dense due to fillin.How r, a significant amount of decay is observed, indicating that it is possible to increase the efficiency of the method in this area.
One step towards ach not dense, its inverse eve ieving this goal is to compute an approximate inverse using a factorized sparse approximate inverse.In this work, we use the approach suggested by Salkuyeh in [26].Salkuyeh's approach is based on Kolotilina and Yeremin's factorized sparse approximate inverse algorithm, FSAI [27,28].The goal is to determine a factorization of the form where i A denotes the ith principal submatrix of A ,

 
L ij G g  , and . In Salkuyeh's paper, these self-preconditioned minimum residual (MR) algorithm with dropping in the search direction (similar to that proposed by Chow and Saad in [29]), starting with an initial sparse guess for the solution and iterating while the solution has fewer than the specified number of nonzeros, denoted by the value lfil.In our work, lfil iterations of the MR algorithm are performed in order to avoid potential infinite loops.Then, structural requirements are enforced; namely, the inverse is required to have the same structure as the matrix systems are solved using a j D .In addition, thresholding is used to further reduce the l-in in the inverse.fil

Cost of Computing the Approximate Inverse ompute sparse appr
The bulk of the cost of the self-preconditioned MR algorithm with dropping in the search direction can be measured in terms of the sparse matrix-sparse vector products and the sparse vector-sparse vector products.One sparse matrix-sparse vector product is required to compute the initial residual, and each iteration of the algorithm requires one sparse matrix-sparse vector product and two sparse vector-sparse vector products.
The MR algorithm is used to c oximate solutions to the system in ( 13)-( 14).For each i, Copyright © 2013 SciRes.AJCM these systems have coefficient matrices that are i i  , where i ranges from 1 to n , with n the size of the - trix whose inverse is bein compu .Thus, the MR algorithm is run 2n times.
The final m r contribu ma g ted ajo tions to the cost of the approximate inverse are the computation of ( 1

Approximate Inverse
following briefly examine complexity of the approximate inverse for the wavelet multigrid method.The value of lfil used to construct the approximate inverse of j D depends on the type of problem, the type of wavelet used, the fine grid size, and the level in the multigrid method.For all problems and wavelets, post-processing is done to ensure that 1 j D  has the same structure as j D , and thresholding is use on all matrices (except for 1 d j D  for diffusion problems using Haar wavelets).Note this is an improvement on the work in [1].
Thresholding values that for the two diffusion examples in Section 5.1 are given in Tables 1 and 2. Note that the thresholding values for In term ize a s of the lfil chosen for FSAI, for the diffusion problem with oscillations in the x-direction, for the modified wavelet multigrid method using 9/7 symmetric biorthogonal wavelets lfil ranges from 3 to 6, and using the 10/6 symmetric biorthogonal wavelets lfil ranges from 4 to 10.For that same problem, for the wavelet multigrid method using Haar wavelets lfil ranges from 9 to 11, and using Daubechies 4 wavelets lfil ranges from 3 to 10.For the diffusion problem with oscillation along diagonals, for the 9/7 symmetric biorthogonal wavelets lfil ranges from 1 to 6, and for the 10/6 symmetric biorthogonal wavelets lfil ranges from 1 to 9. For that same problem, using Haar wavelets lfil ranges from 4 to 11, and using Daubechies 4 wavelets lfil ranges from 3 to 10.
Thresholding values for the three advection-diffusion examples in Section 5.2 are given in Tables 3-5.Note that the thresholding values for 1 j D  (and 1 j L  in the case of parabolic flow (9/7) and recirculant flow (10/6)) for the 9/7 and 10/6 symmetric biorthogonal wavelets depend on the fine grid size and the level of the multigrid method.Similarly, the thresholding values for L  in the case of recirculant flow) for Daubechies 4 wavelets depend on the fine grid size and the level of the multigrid method.
For the advection-diffusion problems with moderately dominant advection, for parabolic flow, for the modified wavelet multigrid method using 9/7 symmetric biorthogonal wavelets lfil ranges from 2 to 6, depending on the fine grid size and on the level of the multigrid, and using 10/6 symmetric biorthogonal wavelets, lfil ranges from 2 to 7. Using the wavelet multigrid method with Haar wavelets, lfil ranges from 2 to 12 and using it with Daubechies 4 wavelets, lfil ranges from 3 to 6.For recirculant flow, lfil ranges from 2 to 5 for 9/7 symmetric biorthogonal wavelets and from 2 to 4 for 10/6 symmetric biorthogonal wavelets, again depending on the fine grid size and on the level.For Haar wavelets, lfil ranges from 3 to 8, and for Daubechies 4 wavelets, it ranges from 2 to 10.For the problem with skewed vorticity, lfil ranges from 4 to 7 for 9/7 symmetric biorthogonal wavelets and from 2 to 6 for 10/6 symmetric biorthogonal wavelets.Using Haar wavelets, lfil ranges from 2 to 10, and using Daubechies 4 wavelets, it ranges from 4 to 9.

Storage and Other Computational Issues
The bulk of the remaining computational work occurs in the construction of the intergrid transfer and coarse grid operators.The construction of the intergrid transfer and coarse grid operators each requires two sparse matrixsparse matrix products and one sparse matrix-sparse matrix difference.
The storage requirements of the coarse grid and intergrid transfer operators are minimized by using sparse matrix storage techniques, resulting in storage requirements of the order of the number of nonzero elements in each matrix.

5
This section describes the numeri the modified wavelet multigrid method using symmetric biorthogonal wavelets to two diffusion problems with highly oscillatory coefficients and to three advectiondiffusion problems with moderately dominant advection.We compare the convergence of the modified wavelet multigrid method using both 9/7 and 10/6 symmetric biorthogonal wavelets with the wavelet multigrid method using both Haar wavelets and Daubechies 4 wavelets.For all problems, numerical results are analyzed using, for the fine grid in the interior, a 32 × 32 grid, leading to a 1024 × 1024 matrix; a 64 × 64 grid, leading to a 4096 × 16384 matrix.In all problems, the stopping criterion for all multigrid methods is where k r is the residual obtained from the kth iteration of the method.

Diffusion Problems
First, we look at the diffusion problem   , giving oscillation along diagonals.Results are examined using both 9/7 and 1


x y  0/6 symmetr s-Seidel ite ethod is used as the (16) where is the unit square centere ic biorthogonal wavelets.The standard Gaus rative m smoother.
Tables 6 and 7 compare the average convergence factor per cycle of the modified wavelet multigrid method using symmetric biorthogonal wavelets with that of the wavelet multigrid method.These tables demonstrate results using a fixed coarsest grid size of 8 × 8 for each fine mesh size.
For the diffusion problem with oscillation in x, as well as the diffusion problem with oscillation along diagonals, the convergence of the modified wavelet multigrid method using symmetric biorthogonal wavelets is relatively independent of the fine grid size, as it is for the wavelet multigrid method using orthogonal wavelets.Convergence results for the method using both the 9/7 symmetric biorthogonal wavelets and the 10/6 symmetric biorthogonal wavelets are comparable to the results using the orthogonal wavelets for both problems.

The Advection-Diffusion Problem
Here, we are investigating the problem .In this problem, difficulties with m me e encountered due to the of the nts haracteristics .W ultigrid fact that some compone of the solution oscillate along c [30,31] e apply the wavelet multigrid method to these problems to overcome this difficulty, since application of the wavelet operator preserves the characteristics of the original problem.
To discretize, we use the standard five-point centered discretization for the diffusion term and a first order upwind scheme for the advection part of the equation.Although using first order upwind introduces artificial diffusion into the solution of the order of the mesh size squared, it provides a convenient test of the effectiveness of the modified wavelet multigrid method.Symmetric Gauss-Seidel is used as the smoother in order to ensure that sweeps are performed in the direction of the characteristics over the entire flow field.Results shown use 2 10 .
First, we have a comparison of the methods for (16), where Note that the discontinuous boundary condition will give rise to a boundary layer near the left-hand boundary.Also, the characteristics are parabolic, resulting in flow entering and exiting through the left-hand boundary.For both the modified wavelet multigrid method with symmetric biorthogonal wavelets and the wavelet multigrid method, convergence appears to be relatively independent of the fine mesh size, as can be seen in Table 8, which compares the average convergence factor per cycle using a coarsest grid size of 8 × 8.
In a second example, giving recirculant flow (i.e., closed characteristics), and   f x is defined in (17).The modified wavelet multigrid method using symmetric biorthogonal wavelets and the wavelet multigrid method both have a convergence rate that is relatively independent of the fine grid size.For this problem, the modified wavelet multigrid method  ts.The are shown in able 9, which compares the average convergence factor the modified wavelet multigrid method using symmetric biorthogonal wavelets and the wavelet multigrid method converge rapidly, with a conv rgence rate that is relatively independent of the fine grid size, as sh e modificatio elet multig od to use mmetric bior al wavelets ha shown to be effective and utile.The results for the modified wavelet elet multigrid ently applied in many cases through et multigrid me g both s bechies 4 wavele results Th sy T per cycle using a coarsest grid size of 8 × 8.
The final example uses the boundary conditions given by (17), but the advection component has both closed characteristics and a skewed flow, so that it does not line up with the grid.Here,  10, which displays the average convergence factor per cycle for the modified wavelet multigrid method using both 9/7 and 10/6 symmetric biorthogonal wavelets and for the wavelet multigrid method using both Haar wavelets and Daubechies 4 wavelets.Again, the coarsest grid size for all trials is 8 × 8. multigrid method using the 9/7 and 10/6 symmetric biorthogonal wavelets are comparable to those obtained by using the wavelet multigrid method with Haar and Daubechies 4 wavelets.The properties of biorthogonal wavelets should permit the modified wav

Conclusion
n of the wav thogon rid meth s been method to be effici the use of compression.The results shown in this paper have demonstrated that it was worthwhile to further explore applying biorthogonal wavelets to the wavelet multigrid method.

j
is orthogonal due to the properties of j H and j G .In the discrete context, the wavelet operators are computationally efficient.With respect to the Haar multiresolution analysis, application of the low-frequency operator   j H to a vector in involves only operations.The same holds for the high-frequency operator n

2 LHH
satisfy certain conditions including forming a Riesz basis of    .For symmetric biorthogonal multiresolution analyses, these wavelet bases (and the bases for the scaling spaces) are symmetric.Define j H and j G to be the operators that transform the basis of the space j  and j G  to be the operators that transform the basis of the sp e ac j V  to the bases of the spaces  and j G  are as follows: fine grid, apply the wavelet transform j   to both sides of the equation and use the orthogonality of j

) 1 D
 , which requires n scalar inverses, and (2) U G D , which requires n z scalar products, where nnz number of

D
 forDaubechies 4  wavelets depend on the fine grid s nd the level of the multigrid method.
and a 128 × 128 grid, leading to a 16384 × We is the unit square centered at will examine two cases where