A new multilevel smoothing method for wavelet-based algebraic multigrid poisson problem solver

A new multilevel smoothing method for wavelet-based algebraic multigrid poisson problem solver

Fabio Henrique Pereira^I; Kleber Rogério Moreira Prado^II; Silvio Ikuyo Nabeta^III

^IIndustrial Engineering Post Graduation Program – UNINOVE, São Paulo, SP, Brazil, [e-mail]

^IIIndustrial Engineering Post Graduation Program – UNINOVE, São Paulo, SP, Brazil, [e-mail]

^IIIElectrical Machine and Drives Lab – GMAcq – USP, Escola Politécnica, São Paulo, SP, Brazil, [e-mail]

ABSTRACT

In contrast to the standard algebraic multigrid, the Wavelet-based Algebraic Multigrid method relies more strongly on the smoothing method because the coarse spaces are chosen a priori. So, it is very important to develop new smoother methods, especially for those cases where the classical Gauss-Seidel smoothing method does not give good results. This paper proposes a new multilevel smoothing approach based on projection technique. The proposed smoothing method was applied to smoothing the error in a linear systems issued from finite element solutions of the elliptic equation and the results compared with those obtained from the Gauss-Seidel method.

Index Terms: Algebraic Multigrid, Gauss-Seidel Method, Multilevel Smoother, Poisson Problem, Projection Technique, Smoothing Method.

I. INTRODUCTION

The algebraic multigrid method is a well-grounded numerical technique for solving large sparse linear systems of algebraic equations. The method can be used as an iterative solver or as a preconditioner for one appropriate iterative method. This technique explores the use of some levels to eliminate the high-frequency or oscillatory components of the error, which are not efficiently removed by basic iterative methods, such as Jacobi or Gauss-Seidel [1], [2].

However, the use of the standard multigrid procedure presents some difficulties. The first one is related to the choice of the strength threshold in the coarsening process. The coarsening scheme, for traditional multigrid, can lead to computational complexity growth as the problem size increases, resulting in an elevated memory use and execution time, and in a reduced scalability [3]. Moreover, the coarsening process is inherently sequential in nature that makes difficult its implementation on distributed memory machines [4]. On the other hand, the classical algebraic procedure cannot directly be applied to problems in which the coefficient matrix violates the M-matrix property like, for example, ungauged edge-based finite-element (FE) analysis. In this context, simple approaches such as the use of a shifted coefficient matrix have a limited efficiency [5]. Additionally, the method performs poorly or may break down in elliptic problems with highly discontinuous, oscillatory or anisotropic coefficients [6]-[8].

The investigation of the wavelet decomposition to obtain coarse grid, restriction and prolongation operators, has been the main subject of many multilevel approaches that are aimed to overcome those difficulties [9]-[16]. The similarities between the multigrid methods and wavelets were first brought out by Briggs and Henson [9]. They observed that the space of fine grid vectors in multigrid scheme and the space of highest resolution in the wavelet multiresolution analysis are correlated. Thenceforward, several works have exploited these similarities [10]-[17]. In [15] it is made an exploration of the Briggs's idea to develop wavelet base interpolation and restriction operators in the geometric multigrid context. It compares the results with the geometric multigrid schemes, with conventional intergrid operators, and obtains results similar to the ones obtained with the standard procedure. More recently, other wavelet-based multigrid approaches have shown the efficiency and potentiality of this combination in different areas. Especially in problems related to the computation of electromagnetic fields, the combination of multigrid and wavelets has revealed to be very efficient and promising, presenting in many cases a better performance than some of the most advanced and current traditional multigrid algorithms [19].

As in any multilevel approach, it is very important to ensure an efficient interaction between the smoother and the coarse-grid correction in the wavelet-based algebraic multigrid method (WAMG). In contrast to the standard algebraic multigrid, which uses some simple relaxation scheme and enforces this interaction by choosing the coarser levels and the transfer operators appropriately [1], [21], the WAMG relies more strongly on the smoothing method because the coarse spaces are chosen a priori. So, the WAMG should select a suitable smoothing method to ensure an accurate error representation in the coarsest level. Concerning to this point of view, the WAMG are similar to the traditional geometric multigrid method [1], [21].

In fact, the difficult in the choice of the suitable smoother for geometric multigrid method are not new and many works have already deal with it [22], [23]. In [22], four smoothing methods were compared in the multigrid context to solve the Navier-Stokes equation. Analyses of robustness and efficiency were carried out and the authors concluded that the most robust smoother is also less efficient and has less parallelization potential. Moreover, the paper states that the choice of the smoother method has a strong influence on the robustness and efficiency of a multigrid method and, consequently, further development is necessary. Other works have presented similar conclusions, especially in parallel computation sense [24].

Specifically in the WAMG context and in the computation of electromagnetic fields problems, investigations about the choice of the most appropriate smoother are still very incipient [25]-[27]. In particular, common smoother types do not work properly in the case of the single layer potential [28], [29].

In face of this, it is very important to develop new smoothing methods with suitable parallel properties as an alternative for the classical Gauss-Seidel smoothing method. This paper proposes a new multilevel smoothing approach based on projection techniques which search an approximation in the high frequency subspace created by the high-pass filters in the discrete wavelet transform. The proposed method was applied for smoothing the error in a linear system issued from finite element solutions of the elliptic equation and the results compared with those obtained from the Gauss-Seidel (GS) method. An analysis of the high frequency error after some steps of the smoother is also presented showing the efficiency of the multilevel smoother, mainly for large problems.

II. THE SMOOTHING PROPERTY OF THE STATIONARY ITERATIVE METHODS

The stationary iterative methods like Damped Jacobi and Gauss-Seidel have a smoothing property which is very well exploited by the multigrid technique. Those iterative methods are very efficient for removing the high frequency components of the error while keeping unchanged the smooth components [21].

The convergence of those iterative approaches is defined by the eigenvalues of its iteration matrix. For the Damped Jacobi, in a one-dimensional elliptic problem discretized by Finite Difference Methods, for example, the iteration matrix P_ω is given by [1]

in which D is the main diagonal part of A, I is the identity matrix and M = A - D.

The eigenvalues of the matrix P_ω for this example are

where n is the dimension of A.

Then, if h=1/n denotes the mesh spacing:

It means that the first and largest eigenvalue of P_ω will be very close to one when h→0,

As the eigenvalue is related to λ₁ an eigenvector in the form

the convergence for Damped Jacobi is slow for the lowest frequency. This property is shared with the others stationary iterative methods [21].

As the smoothing components look like more oscillatory in a coarse grid, the transfer of the problem to a grid with larger mesh size (or larger value of h in equation (3) either a lower value of n in equation (5)) will improve the convergence of the stationary iterative method. Then, if we are sure that the error of our approximation has become smooth after some iteration steps, we can approximate this error by a suitable procedure on a coarser grid [1].

The recursion of this two grid process produces the multigrid method.

III. THE REQUIREMENTS OVER THE SMOOTHER IN WAMG

In order to analyze the requirements over the smoothing methods in WAMG, we have to understand the DWT behavior and the effects of the use of this technique in the algebraic multigrid context. Usually in the literature, the discrete wavelet transform is defined using the concept of filter bank [30]. In such cases low-pass and high-pass filters are chosen to obtain the low frequency and the high frequency terms of the input signal, respectively. In this procedure, if the discrete input signal x = {x[n]}_n∈ is sufficiently smooth then the high frequency terms are close to zero and the output low frequency component will be a good approximation of the signal.

To better understand the effects of this technique in the WAMG context, we will suppose that there are only two grids: the original grid and a coarse grid.

The WAMG method uses only the low-pass filters to produce an approximation of the original coefficient matrix. This approximation will represent the system of equation in the coarse grid and it will be used to solve the residual equation on the coarse grid. Then, the residual equation solution is used in the correction of the error in the original grid.

As the WAMG coarse grid is created by the low-pass filters, the coarse grid error correction will only be accurate if the error is smooth. The fact that the error is smooth means that the high frequency components become small after a few stationary iterative iterations whereas the low frequency components hardly change [21]. Therefore, the smoothing method should to be able to remove or, at least, to lessen the high frequency components of the error.

IV. THE PROPOSED MULTILEVEL SMOOTHER METHOD

In a projection process an approximation for the solution of a linear system is extracted from a subspace which is the subspace of candidate approximations, or search subspace. If m is the dimension of , then the residual vector b – Ax is constrained to be orthogonal to m linearly independent vectors in order to obtain such an approximation. This orthogonal condition defines another subspace , of dimension m, called the subspace of constraints [32], [33].

Mathematically, if we know an initial guess, x₀∈ , the approximate solutions can be defined as

in which r₀ = b – Ax₀.

As an ideal smoother should remove the high frequency error components it is interesting that this smoother searches an approximation to the linear system solution in a high frequency subspace.

In other hand, the multiresolution analysis in discrete wavelet transform decompose the original grid (space) V in two subspaces V₁ and W₁ such that [30]

such as

where Φ and ψ are the scaling and wavelets function, respectively. Therefore, we can develop a smoothing method based on a projection approach with,

which will search the approximations to the solution in the high frequency subspace W₁. As result the obtained error should be smooth in the sense of the Discrete Wavelet Transform.

In a matrix representation, denoting by W = [w₁ , ..., w_m] an n x m matrix whose column-vectors form a basis of W₁, a two-level smoothing method based on projection technique can be represented by the following algorithm:

Algorithm: Two-level smoothing method

1. Get the current approximation x_i

2. r^h = b - Ax_i

3. r^2h = W^T r^h

4. e^2h = (W^T AW)^-1r^2h

5. x_i+1 = x_i + We^2h

6. Return x_i+1

The recursion of this two-level process produces a new multilevel smoothing method (MS) which has shown to be efficient in smoothing the high frequency error.

V. DESCRIPTION OF THE TEST PROBLEMS

In order to analyze the efficiency of the proposed smoother the method was applied to solve linear systems issued from finite element solutions of the elliptic equation (12), defined in a n × n – point uniform grid on the square domain, with Dirichlet boundary conditions:

Three different discretizations were applied and the resulting matrices are presented in Table I.

For analyzing the WAMG convergence behavior we assume that the exact solution is known (13), and we use it to supply the right hand side vector b.

VI. NUMERICAL RESULTS AND COMMENTS

First of all, the error is analyzed after one, two and three smoother iterations, in order to identify the efficiency of the proposed method in smoothing the high frequency error. So, the smoother was used alone, without the WAMG method. As in the WAMG context is supposed to there is no predefined geometric information about the problem, the error is analyzed in an algebraic sense, what means the error is slowly reduced by the smoothing method [21]. The exact solution presented in (13) also was defined in the algebraic point of view. The corresponding results are presented in Fig. 1. The classical Gauss-Seidel was used for comparison.

It is possible to see from Fig. 1 that the proposed multilevel smoother is more affective in the reduction of the high frequency error for this case.

As a zero initial guess was adopted the initial error is equal to the exact solution .

The numerical reduction of the high frequency components of the error was also evaluated. In this case, the projection of the error in the high frequency subspace W, denoted as We, was analyzed and the results are reported in Fig. 2. Again here, only the application of the smoother has been considered.

As the WAMG method accomplishes the coarse grid correction on the low frequency subspace V, an ideal smoother method should remove the high frequency projection We. The Euclidian norm of the We vector presented in Fig. 2 shows the efficiency of the proposed approach. In this case, the smaller the norm of the We vector the better is the smoothing result.

Finally, the multilevel smoother and the Gauss-Seidel method were compared within the original WAMG algorithm. The WAMG method using the two smoothers was applied to solve the linear systems with the matrices presented in Table I. The variations forward, backward and symmetric of the smoothing method Gauss-Seidel were tested for the WAMG. In the forward (backward) approach the smoothing process is realized in the increasing (decreasing) order of the matrix rows. For the symmetric version, the smoother is applied in an alternate form in even and odd rows. The number of WAMG iterations needed for convergence, and the setup and solver times (in seconds) for each configuration are presented in Table II. The objective of this test was to observe the effect of the reduction of the high frequency error components in the WAMG convergence behavior.

The number of iterations corresponds to the number of V-cycles necessary to reduce the Euclidean norm of the residual vector to the order of 10^-6.

From the results in Table II, two main aspects must be highlighted: the multilevel smoother (MS) is computationally more expensive than the Gauss-Seidel (GS) method. It is an expected result because in MS all the levels are visited even though the method was applied only on the first level. Moreover, the WAMG setup time is about two times higher, because it involves the construction of a coarse matrix for the smoother, as defined in step 4 of the algorithm for the high frequency subspace W. However, the reduction number of WAMG iterations means we can get an overall reduction in computational time for large systems. The second aspect is related to the WAMG convergence which preserves the typical mesh-independent convergence behavior of the multigrid cycles.

VII. CONCLUSIONS

The proposed method seems to be very efficient in smoothing the high frequency error. The results in Fig. 1 show that the errors produced by the multilevel approach appear less oscillatory and smaller than those created by the classical Gauss-Seidel method in the three first iterations.

Although being computationally more expensive than the Gauss-Seidel method, the multilevel smoother produces a reduction in the number of iterations that means we could get an overall reduction in computational time. In the tested cases, the WAMG convergence was better when the magnitude of the high frequency components of the error was reduced by the smoother, as presented in Fig. 2. In this aspect the multilevel approach produced the best results, mainly for the larger problem.

Additionally, the multilevel approach has suitable characteristics for parallelization in distributed memory computers, since it can be applied independently in the rows and columns of the matrix exactly in the same way presented in [19] for the WAMG method. So, it can be a good alternative to inherently sequential approaches such as the smoothers based on incomplete LU decompositions, which are robust but has less parallelization potential [22]. The application of the proposed method in the parallel context should be addressed in further work.

Received 9 June 2011; for review 9 June 2011; accepted 30 Sept. 2011

• [1] W. L. Briggs, V. E. Henson, and S. F. McCormick, A Multigrid Tutorial, second ed., SIAM, California, 2000.
• [2] Q. Chang, Y. S. Wong, and H. Fu, "On the algebraic multigrid method", J. Comput. Phys. 125:279292, 1996.
• [3] G. Haase, M. Kuhn, and S. Reitzinger, "Parallel Algebraic Multigrid Methods on Distributed Memory Computers", SIAM J. Sci. Comput., 24:410-427, 2002.
• [4] V. E. Henson, and U. M. Yang, "BoomerAMG: a Parallel Algebraic Multigrid Solver and Preconditioner", Appl. Num. Math., 41(1):155-177, 2002.
• [5] T. Mifune, T. Iwashita, and M. Shimasaki, "New algebraic multigrid preconditioning for iterative solvers in electromagnetic finite edge-element analyses," IEEE Trans. Magn., 39(3):6771680, 2003.
• [6] B. Engquist and E. Luo, "Convergence of a multigrid method for elliptic equations with highly oscillatory coefficients", SIAM J. on Numerical Analysis, 34:22542273, 1997.
• [7] W. L. Wan, "Interface preserving coarsening multigrid for elliptic problems with highly discontinuous coefficients", Numerical Linear Algebra with Applications, 7(8):727742, 2000.
• [8] Y. Xiao, P. Zhang and S. Shu, "Algebraic multigrid methods for elastic structures with highly discontinuous coefficients", Mathematics and Computers in Simulation, 76(4):249-262, 2007.
• [9] W. L. Briggs and V. E. Henson, "Wavelets and multigrid", SIAM J. Sci. Comput. 14:506510, 1993.
• [10] B. Engquist and E. Luo, "The multigrid method based on a wavelet transformation and Schur complement", Unpublished.
• [11] A. Rieder, "Multilevel methods based on wavelet decomposition", East-West J. Numer. Math., 2(4):313-330, 1994.
• [12] D. De Leon, Wavelet Operators Applied to Multigrid Methods (Ph.D. Thesis), UCLA Mathematics Department CAM Report 00-22, June 2000.
• [13] D. De Leon, "A new wavelet multigrid method", J. Comput. Appl. Math., 220:674-685, 2008.
• [14] D. De Leon, "A Wavelet Multigrid Method Applied to the Stokes and Incompressible Navier-Stokes Problems", Journal of Mathematical Sciences: Advances and Applications, 1(3):601-622, 2008.
• [15] A. Avudainayagam, and C. Vani, "Wavelet based multigrid methods for linear and nonlinear elliptic partial differential equations", Appl. Math. Comput 148:307320, 2004.
• [16] F. H. Pereira, S. L. L. Verardi, S. I. Nabeta, "A Wavelet-based Algebraic Multigrid preconditioner for sparse linear systems", Appl. Math. Comput., 182:1098-1107, 2006.
• [17] V. M. Garcia, L. Acevedo, and A. M. Vidal, "Variants of algebraic wavelet-based multigrid methods: Application to shifted linear systems," Appl. Math. Comput., 202(1): 287299, 2008.
• [18] F. H. Pereira, et. al., "A Wavelet-based Algebraic Multigrid Preconditioning for Iterative Solvers in Finite Element Analysis". IEEE Trans. on Magn., 43(4): 1553-1556, 2007.
• [19] F. H. Pereira, M. F. Palin, S. L. L. Verardi, S. I. Nabeta, "A Parallel Wavelet-based Algebraic Multigrid black-box Solver and Preconditioner", 16th Compumag - Conference on the Computation of Electromagnetic Fields, 2007, Aachen.
• [20] F. H. Pereira, S. I. Nabeta, "Wavelet-Based Algebraic Multigrid Method Using the Lifting Technique", J. of Microwaves, Optoelectronics and Electromagnetic Applications, 9(1): 1-9, 2010.
• [21] U. Trottenberg, C. W. Oosterlee, and A. Schuller, Multigrid, Academic Press, New York, 2001, pp. 45-56.
• [22] S. Zeng and P. Wesseling, "Numerical Stydy of a Multigrid Method with Four Smoothing Methods for the Incompressible Navier-Stokes Equations in General Coordinates". In N. D. Melson,Th. A. Manteuffel, and S. F. McCormick, editors, Sixth Copper Mountain Conference on Multigrid Methods, 691-708. NASA Langley Research Center, 1993.
• [23] J. Blazek, C.-C. Rossow, N. Kroll and R.C. Swanson, "A comparison of several implicit residual smoothing methods in combination with multigrid". Lecture Notes in Physics, 414:386-390, 1993.
• [24] M. F. Adams, M. Brezina, J. J. Hu, and R. S. Tuminaro, "Parallel multigrid smoothing: Polynomial versus Gauss-Seidel". J. Comput. Phys, 188:593610, 2003.
• [25] K. Watanabe, S. Fujino, H. Igarashi, "Multigrid Method With Adaptive IDR-Based Jacobi Smoother", IEEE Trans. on Magn., 47(5):1210-1213, 2011.
• [26] F. H. Pereira, A. C. M. Junqueira, R. Perrussel, S. I. Nabeta, "An investigation of the requirements over smoother methods in Wavelet-based Algebraic Multigrid". In: 13th Biennial IEEE Conference on Electromagnetic Field Computation, Athens, 13:1-4, 2008.
• [27] F. H. Pereira, S. I. Nabeta, "A New Multilevel Smoothing Method for the Wavelet-Based Algebraic Multigrid." In: 17th Conference on the Computation of Electromagnetic Fields, Florianópolis, 2009.
• [28] D. Pusch, Efficient Algebraic Multigrid Preconditioners for Boundary Element Matrices, Phd thesis. Johannes Kepler University Linz, Austria, 1997.
• [29] F. H. Pereira, M. M. Afonso, J. A. Vasconcelos, and S. I. Nabeta, "An Efficient Two-Level Preconditioner based on Lifting for FEM-BEM Equations", J. of Microwaves, Optoelectronics and Electromagnetic Applications, 9(2):78-88, 2010.
• [30] A. Jensen, A. la Cour-Harbo, The Discrete Wavelet Transform, Ripples in Mathematics, Springer, Berlin, 2001.
• [31] T. K. Sarkar, M. Salazar-Palma, C. W. Michael, Wavelet Applications in Engineering Electromagnetics, Artech House, Boston, 2002.
• [32] Y. Saad, Iterative Methods for Sparse Linear Systems, 2^nd ed, SIAM, Philadelphia, 2003.
• [33] O. Axelsson, Iterative solution methods, Cambridge University Press, New York, NY, 1995.

Publication Dates

• Publication in this collection
  09 Mar 2012
• Date of issue
  Dec 2011

History

• Reviewed
  09 June 2011
• Received
  09 June 2011
• Accepted
  30 Sept 2011