Critical Review of Preconditioners for Helmholtz Equation and their Spectral Analysis

Objectives: The study reviews the iterative solvers for solving Helmholtz equation. It also surveys the preconditioners for Krylov subspace methods for discrete Helmholtz equation. Numerical experiments based conclusions are drawn, after a thorough comparison of various preconditioners. Methods/Statistical Analysis: Five point finite difference scheme for two-dimensional Helmholtz equation is employed to obtain discrete analogue of Helmholtz equation. Krylov subspace method for is considered for iterations. Matrix based preconditioner ILU with variety of fill-in and operator based preconditioner complex shifted Laplace preconditioner (CSLP) are tested. Different shifts are proposed and experimented in CSLP. Optimization of shifts is recommended on basis of results. Spectral analysis is performed to validate results. Findings: For small wave number problem, ILU type preconditioned performs better and is fairly comparable with CSLP. However, increasing wave number makes unaffordable in terms of computational cost. The pure imaginary shift ( ) 0,1 and complex shift ( ) 1, / 4 π costsless time and number of iterations compared to other choices of shift. These choices are recommended when preconditioner is inverted directly. Approximate inversion may affect iterations, but will save solve time for large problems. Extensive spectral analysis, validating results, highlights this work. Application/Improvements: Dependency upon wave number is decreased with this recommended choice of shifts in CSLP. The number of iterations increases at very low rate when wave number is increased at rate of 2.


Introduction
Helmholtz equation often arises in the study of physical problems involving partial differential equation (PDE) in both space and time; in this paper the time-independent case considered.
Many problems related to steady state oscillations (mechanical, acoustical, and thermal) are modeled by the Helmholtz equation 1 . Wave scattering have many applications in physics, engineering and science. Examples include seismic imaging 2 , radars, electromagnetism 3 , bio medical imaging 4 (ultrasound), road-speed sensors etc. Wave scattering phenomena is mostly modeled by Mathematicians in the form of the Helmholtz equation 5 .
The linear system arising from a discretization of the Helmholtz equation in two dimensions (2D) or three dimensional (3D), domain, converges typically characterized by indefiniteness of the Eigen values of the corresponding coefficient matrix. With such a property, an iterative method either basic or advanced, encounters convergence problems. The method usually converges very slowly or divergences, for high frequency problems, as the linear system becomes extremely large and an indefinite. This makes the problem even harder to solve.
, c x y f λ = the wavelength and ( ) , c x y the speed of sound.
Two problems are considered here. The first one has Dirichlet's conditions and leads to symmetric positive definite system; the second one has Neumann (Sommerfeld) conditions 6,7 and leads to non-symmetric indefinite systems.

Problem Description
The Helmholtz problem considered this paper is nonhomogeneous, undammed (i.e.
This means that the waves propagate from the center of domain outwards. The Somerfield radiation conditions 8-10 of first order are imposed that is

Discretization
The problem under discussion can be solved by many methods including finite element method (FEM) and finite difference method (FDM). Higher order finite difference schemes have also been used. As, considered problem is defined on simple square domain, hence conventional FDM is chosen for discretization. Discretization is performed on a square grid, with uniform mesh size in both x and y directions problems.
In FDM, the differential terms are approximated by finite differences using Taylor's polynomials. The differential terms in Equation: (1) are approximated by finite differences as .
Using above differences for differential terms on the square grid, the discrete analogue of Equation: (1) will be ( ) The Sommerfeld radiation condition of first order given in Equation: (1) are descretized using by central difference at respective directional derivative, for example at right boundary, it is given as The stencil for nodes/points next to boundary will take form as The Equation: (A) represents stencil for node/point next to left side boundary, and Equation: (B) represents stencil for node next to left and lower boundary, respectively. The resultant linear system is written as

Properties of Linear System
The coefficient matrix in Equation: (9) is large, sparse, complex-valued due inclusion of Sommerfeld radiation condition and for bigger value of wave number, it is indefinite, as discussed in spectral analysis of the system. Furthermore the matrix is non-symmetric and non-Hermitian.

Iterative Methods for the Helmholtz Equation
The coefficient matrix A in Equation: (9) is mostly indefinite. By the "indefinite", we mean that the real part of the eigenvalues of A lie in both positive and negative half plane in the complex plane. The dimension of linear system increases when wavenumber increases, therefore more nodes are required to resolve higher wavenumber. The direct solvers are of impractical for many obvious reasons. Iterative methods are mostly preferred for large and particularly sparse linear systems. Few are discussed.

Basic Iterative Methods
The basic iterative methods are also called fixedpoint iterations which are based on the splitting , , After substitution into Equation: (9) and few steps of simplification will give iterative process for ( ) To improve the convergence of the fixed point iterations, a relaxation factor ω is often introduced in Equation: (10).
A standard approach to investigate the convergence of this type of methods is by analyzing the amplitude of the Fourier nodes between two successive iterations. For indefinite Helmholtz problems, it can be shown that there exist no r ω such that the error is reduced during the iterations. This means that Jacobi iteration always divergences if applied to the indefinite Helmholtz equation. Similarly result holds for the Gauss-Seidel iteration method.

Krylov Subspace Methods
For large, sparse matrix the Krylov subspace are very popular. The methods are developed on construction of iterants in the subspace.
The space is called the Krylov subspace of dimension j , associated with A and 0 r , and initial residual 0 Among methods which are based on construction the Krylov subspaces, a Conjugate Gradient (CG) 11,12 , COCG, GMRES, CGS, Bi-CGSTAB, and QMR are popular choices. For non-symmetric linear systems, Krylov subspace can be built from Arnoldi's process, which leads to GMRES. GMRES is optimal method; it reduces the 2-norm of the residual at the every iteration. GMRES, however, require long recurrences, which is usually limited by the available memory. A remedy is by restarting, which some-times lead to slow convergence or stagnation.
Since CG is standard formula for Krylov method, however CG is limited in application to only symmetric and definite problems. On the contrast, GMRES has the property it can be used for indefinite and non-symmetric problems as well. It is necessary to use preconditioner to get better results of Krylov subspace.

Eigen Decomposition
We can decompose a matrix into parts that give us some information about the matrix, which is not directly accessible from the one 2D array representation.
One such decomposition is Eigen-Decomposition, where we decompose a matrix, into a set of eigenvectors and Eigen values.
A matrix with all positive eigenvalues being positive is called positive definite. A matrix whose eigenvalues are all positive or zero is called positive semi-definite. If all eigenvalues are negative, matrix is called negative definite, and semi definite is similarly defined.

Numerical Experiments
For all the experiments, 0 u (zero vector) is used as initial guess, also from here and onwards, the mesh size h is such that for a wave number k , it satisfies 0.625 Kh = (equivalent to 10 grid points per wave length, and 2 k K = ) in all numerical experiments, unless mentioned.
The first result is presented in, the Table (1) shows the performance of the CGS and GMRES algorithm with respect to a variation of K . Increasing the wave number K severely affected the number of iterations and the actual residual.

Preconditioning
In an iterative method preconditioning is often vital component in enabling the solution of such system when the dimension is large. There are many iterative techniques for solving linear system. For large spare linear system 13 , the convergence rate is always a concern for researchers; one can improve significantly the convergent rate by applying appropriate preconditioner. The convergence of iterative methods depends on eigenvalues of the coefficient matrix, it is often advantageous to use preconditioner that transfer the system to one with a better distribution of eigenvalues. The preconditioner is the key to successful iterative solver. In simple words, the preconditioning is technique of transforming a linear system into favorable one for any iterative solver, by multiplying the preconditioner M (on left, right and split) preserving the solution. After choosing a preconditioner M for any system Equation: (9), who's invertible is cheap to compute, and then the transformed system is This is called left preconditioning, which is used in this paper.

ILU Preconditioner
One simple class of preconditioners is obtained by an Incomplete LU factorization of A [14,17], and there are many ways to obtain approximate ILU factorization of A, e.g fill-in ILU and ILU with some tolerance. ILU obtained by restricting the structure of L and U to equal that of A leads to preconditioner, known as ILU(0) which is easy to compute but not effective for Krylov subspace methods. A more accurate ILU factorization of A is obtained by allowing more fill-in the ILU-factorization. Dropping the elements less than some given value in ILU gives rise to ILU (tolerance), known as ILU(0). Some experiments are done in this paper using GMRES with precondition-ers ILU(0) and ILU(tolerance) with tolerance (0.01). Table 2 and 3 shows the numbers of iterations with actual residual as it seems to be decreased with decreasing mesh size and also increasing wave number K. Results in Table  2 and 3 indicate that ILU (0.01) works better than ILU (0) but at the cost of big storage, that is due to more fill-in is allowed, which is indicated by the number of non-zero elements in factorizations L and U.   Table 2     The increase in number of iterations can be seen, where number of iterations is plotted against number of grid points (size of problem). The convergence history for GMRES with ILU (0) and ILU (0.01) as a preconditioner is given in Figure 1, and compared with that of GMRES without preconditioners. Results with shifted Laplace preconditioner are presented in subsequent section. Further, Figure 2 clearly indicates that the number of iterations of GMRES with preconditioner ILU(0.01) are less than both GMRES as well as GMRES preconditioner ILU(0).

Shifted Laplace Preconditioner
Another approach is obtaining preconditioner based on operator, not on matrix. This preconditioner, in context of Helmholtz equation, is called shifted Laplace preconditioner [15][16][17] . This is developed by discretization of the operator respectively. This class starts with a simple Laplace opera-tor M = ∆ which was used as a preconditioner, it seems to be good, and works well, until mesh size is small, for large size of mesh this shift creates lot of problems, lot of unwanted scattered eigenvalues appear, as shown in Figure 3, which shows that for large mesh size this shift is not good choice. Later an additional real term shift was added in the Laplace operator, making this preconditioner resembling more the Helmholtz operator but with an opposite sign as investigated. Later Laplace operator with imaginary shift was introduced and found to be effective for the Helmholtz equation, as shown in Figure 3. The optimal choice of real and imaginary shifts with restriction to be SPD (with positive real parts of eigenvalues of preconditioner) is investigated in detail in 9,10 in the context of the condition number for CG, first for real shift and then generalized to complex shifts. β = -1, for (only) real shifted Laplace preconditioner β 1 = 0, and β 2 = 0 for shifted Laplace preconditioner are concluded as an optimal choice (for details, see 1,9,10 ). Table 4 shows the number of GMRES iterations using four different shifted Laplace preconditioners. An observation is that for very small wave number K , ( ) M − − is similar to that of ( )

1,1 M
, only turns counter clockwise as showed in Figure 5. The comparison of different shifts in CSLP is given in Figure   6, which also shows that ( ) 1,1 M has better spectrum. It is worth mentioning that ILU(0.01) perform comparable to CSLP for small wave numbers at the cost of huge memory, which is not affordable in case of large wave number.

Conclusion
In this study, Krylov subspace iterative methods has been apllied to the Helmholtz problem. Need of preconditioner is highlighted and a survey of preconditioner is presented.
ILU type preconditioners and CSLP are tried. The CSLP preconditioner is found to be very effective to enhance the convergence of Krylov subspace methods considered problem. CSLP with real and imaginary shifts as ( ) 1,1 works better for all ranges of wavenumber. Note that preconditioner is inverted exactly. A spectral analysis confims findings for said shift.