New Preconditioning Techniques for Saddle Point Problems Arising from the Time-Harmonic Maxwell Equations

We study two parameterized preconditioners for iteratively solving the saddle point linear systems arising from finite element discretization of the mixed formulation of the time-harmonic Maxwell equations in electromagnetic problems. We establish some spectral properties of the preconditioned saddle pointmatrices, involving the choice of the parameter.Meanwhile, we provide some results of numerical experiments to show the effectiveness of the proposed parameterized preconditioners.


Introduction
The following time-harmonic Maxwell equations in lossless media and perfectly conduction boundaries [1][2][3][4][5] are partial differential equations of the mixed form with constant coefficients: find the vector field  and the multiplier  such that ∇ × ∇ ×  + ∇ =  in Ω, ∇ ⋅  = 0 in Ω,  ×  = 0 on Ω,  = 0 on Ω. (1) Here, Ω ⊂ R 3 is a simply connected polyhedron domain with a connected boundary Ω, and  denotes the outward unit normal on Ω.The datum  is a given generic source (not necessarily divergence free), and the wave number satisfies  2 =  2 ∈ , where  ≥ 0 is the frequency, and  is a positive permeability parameter.
If the lowest order Nédélec elements of the first kind [1] are used for the approximation of the electric field and the standard nodal elements for the multiplier, discretizing the problem (1), we derive the approximation solution of (1) through solving the saddle point linear system of the form where  ∈   and  ∈   are finite arrays representing the finite element approximations, and  ∈   is the load vector associated with .The matrix  ∈  × is symmetric positive semidefinite with nullity  and corresponds to the discrete curl-curl operator;  ∈  × is a discrete divergence operator with full row rank.For large, sparse, or structured matrices iterative methods are the interesting alternative.Stationary iterative methods solve a linear system (2) with an operator approximating the original one; based on a measurement of the error in the result (the residual), form a "correction equation" for which this process is repeated.While these methods are simple to derive, implement, and analyze, convergence is only guaranteed for a limited class of matrices.Examples of stationary iterative methods are the Jacobi method, Gauss-Seidel method and SOR method [1,2,[6][7][8][9], and HSS splitting iterative method [7,10].However, Krylov subspace methods work by forming an orthogonal basis of the sequence of successive matrix powers times the initial residual (the Krylov sequence).The approximations to the solution are then formed by minimizing the residual over the subspace formed.Krylov subspace methods apply techniques that involve orthogonal projections onto subspaces of the form where K  (,  0 ) is the th Krylov subspace associated with  and  0 (see, e.g., [9,11,12]).

ISRN Applied Mathematics
The basic algorithm within this class is the conjugate gradient (CG) method which has the nice properties that it uses only three vectors in memory and minimizes the error in the -norm.However, the algorithm mainly performs well if the matrix  is symmetric and positive definite.In cases where one of these two properties is violated, CG may break down.For nonsymmetric or indefinite linear systems, CG can be applied to the normal equations since the resulting linear system becomes (positive) definite.Upon application of CG to the normal equations, CGNE and CGNR are obtained [12].The CGNE and CGNR methods, transforming the system to a symmetric definite one and then applying the conjugate gradient method, is attractive for its coding simplicity.The iterations are guaranteed to converge.The drawback is that the condition number of the normal equations equals the square of the condition number of , slowing down the convergence drastically.
MINRES can also be used to solve indefinite symmetric linear systems, as well as it generalization to the nonsymmetric case, GMRES [9,[12][13][14].Both algorithms have the minimization property, but GMRES uses long recurrences.GMRES has the advantage that theoretically the algorithm does not break down unless convergence has been reached.The main problem in GMRES is that the amount of storage increases as the iteration number increases.Therefore, the application of GMRES may be limited by the computer storage.To remedy this problem, a restarted version, GMRES(), can be utilized.Since restarting removes the previous convergence history, GMRES() is not guaranteed to converge.There is no specific rule to determine the restart parameter .In cases characterized by superlinear convergence,  should often be chosen very large which makes restarting much less attractive.Another way to remedy the storage problem in GMRES is by including a so-called "inner iteration" as in GMRESR and FGMRES [9,12,14].
The approximating operator that appears in stationary iterative methods can also be incorporated in Krylov subspace methods such as GMRES (alternatively, preconditioned Krylov methods can be considered as accelerations of stationary iterative methods), where they become transformations of the original operator to a presumably better conditioned one.The construction of preconditioners is a large research area [3,4,8,9,13,[15][16][17][18][19][20][21][22][23].Generally, preconditioning attempts to improve the spectral properties of the system matrix.For symmetric problems, the rate of convergence of Krylov subspace methods like CG or MINRES depends on the distribution of the eigenvalues of A. A key for the rapid convergence of an iterative method for a linear system of the form A =  is the availability of an effective preconditioner P. Each step of an outer iteration for solving the preconditioned linear system P −1 A = P −1  requires the solution of an inner linear system whose coefficient matrix is A. Therefore, convergence of the outer iteration is fast if the eigenvalues of the preconditioned matrix P −1 A are clustered, but careful attention must be paid to the conditioning and eigenvalue distribution of the matrix A itself, which determine the speed of convergence of the inner iteration; see [8] for a comprehensive survey.
For solving the saddle point linear system (2), Greif and Schötzau [19] presented a block diagonal preconditioner where  ∈ R × is a symmetric positive definite.In 2007, Greif and Schötzau [3] applied this type block diagonal preconditioner to the saddle point linear systems arising from the discretized time-harmonic Maxwell equations in mixed form (1). Rees and Greif [20] studied a block triangular preconditioner which has the property that, the more the ill-conditioned (1, 1) block of the saddle point matrix is, the faster a minimum residual solver, such as MINRES, converges.In 2008, Cao [21] gave two augmentation block triangular preconditioners And he has shown that if the nullity of  is , then the preconditioned matrices T −1 + A and T −1 − A have only three distinct eigenvalues 1, (−1 ± √ 5)/2 and 1, (1 ±  √ 3)/2, respectively.Thus, the preconditioned GMRES with either augmentation block triangular preconditioner converges within three iterations.
Based on the preconditioners T + and T − , we study two generalized block triangular preconditioners for iteratively solving the saddle point linear systems arising from finite element discretization of the mixed formulation of the timeharmonic Maxwell equations in electromagnetic problems.Spectral properties and computational performance of the generalized block triangular preconditioners are discussed in detail, involving the choice of the parameters that are considered.Meanwhile, we give the optimal parameter in practice.Finally, numerical experiments show the effectiveness of the proposed preconditioners.
The remainder of this paper is organized as follows.In Section 2, the new block triangular preconditioners with one parameter are proposed, and the spectral distribmtion of the new preconditioned matrices is analyzed.In Section 3, numerical experiments are provided to validate our results in Section 2. Finally, we draw some conclusions.

Block Triangular Preconditioners
It is well known that characterizing the rate of convergence of nonsymmetric preconditioned iterations can be a difficult task.In particular, eigenvalue information alone may not be sufficient to give meaningful estimates of the convergence rate of a method like preconditioned GMRES [9,24].Nevertheless, experience shows that, for many linear systems arising in practice, a well-clustered spectrum (away from zero) usually results in rapid convergence of the preconditioned iteration.
Next we develop some estimates for the eigenvalues of the preconditioned matrices P −1  A and F −1  A, assuming exact solves for the (1, 1) block .We show that, for this "ideal" version of the preconditioner, the eigenvalues of the preconditioned matrix become tightly clustered around 1 when the nullity of  is .

Spectral Properties of the Preconditioned Matrix. We let
represent the saddle point matrices of (10).We assume that  is symmetric and positive semidefinite with nullity  and that  is of size  ×  ( ≤ ) and has full row rank.Note that the assumption that A is nonsingular implies that null() ∩ null() = 0 and null() ∪ null() =   , which we use in our analysis below.
We present the following block triangular preconditioner: We will study spectral properties and computational performance of the preconditioner P  .Also numerical experiments show that the preconditioner P  with a parameter is more efficient than both T + and T − if it is proper to choose the value of the parameter.Following the spirit of the proof of [11,25] we give the following result to describe the spectral distribution.
Theorem 1.Let A be nonsingular, and its (1, 1) block  is singular with nullity .Then  = 1 is an eigenvalue of P −1  A of geometric multiplicity −, and  = (−1± √ 5 − 8)/2(1−2) is an eigenvalue of geometric multiplicity .The remaining 2( − ) eigenvalues satisfy the relation where  are some  −  generalized eigenvalues of the following generalized eigenvalue problem: Proof.Suppose that  is an eigenvalue of P −1  A, whose eigenvalue is (  ,   )  .Then we have Furthermore, it satisfies the generalized eigenvalue problem or As A is nonsingular,  ̸ = 0.The second equality gives  = (1/(1 − 2)) −1 , and substituting it into the first equality gives If  ∈ null(), then (20) implies that from which it follows that  = 1 is an eigenvalue of P −1  A of geometric multiplicity  − .
If  = 1 and P  apparently reduce to T − , then we have the following corollary.

We next consider another block triangular preconditioner
We give the following result to describe the spectral distribution of the preconditioned matrix F −1  A.

Theorem 6.
Let A be nonsingular, and its (1, 1) block  is singular with nullity .Then  = 1 is an eigenvalue of F −1  A of geometric multiplicity , and  = 1/(2 − 1) is an eigenvalue of geometric multiplicity .The remaining − eigenvalues satisfy the relation where  are some  −  generalized eigenvalues of the following generalized eigenvalue problem: Proof.Suppose that  is an eigenvalue of F −1  A, whose eigenvalue is (  ,   )  .Then we have Furthermore, it satisfies the generalized eigenvalue problem or As A is nonsingular,  ̸ = 0.The second equality gives  = (1/(1 − 2)) −1 , and substituting it into the first equality gives It is straightforward to see that any vector  ∈   satisfies (29) with  = 1, and thus  = 1 is an eigenvalue of F −1  A. We claim that the eigenvalue  = 1 has geometric multiplicity .

Corollary 7.
Let A be nonsingular, and its (1, 1) block  is singular with nullity .Then  = 1 is an eigenvalue of M −1 A of geometric multiplicity , and  = −1 is an eigenvalue of geometric multiplicity .The remaining − eigenvalues satisfy the relation where  are some  −  generalized eigenvalues of the following generalized eigenvalue problem: Corollary 8. Let A be nonsingular, and its (1, 1) block  is singular with nullity .Then  = 1 is an eigenvalue of M −1 A of geometric multiplicity , and  = −1 is an eigenvalue of geometric multiplicity .
If  = 1, then we have the following corollary.

Corollary 9.
Let A be nonsingular, and its (1, 1) block  is singular with nullity .Then  = 1 is an eigenvalue of F −1 1 A of geometric multiplicity  + .The remaining  −  eigenvalues satisfy the relation where  are some  −  generalized eigenvalues of the following generalized eigenvalue problem: Corollary 10.Let A be nonsingular, and its (1, 1) block  is singular with nullity .Then the preconditioned matrix F −1 1 A has only precisely one eigenvalue:  = 1 of geometric multiplicity  + .
Remark 11.From Corollaries 8 and 10, if the nullity of  is , it is readily seen that the preconditioned matrix M −1 A has precisely two distinct eigenvalues and that the preconditioned matrix F −1 1 A has only precisely one eigenvalue.Thus, we know that any preconditioned Krylov subspace method such as GMRES terminates in at most two steps if round-off errors are ignored.
Remark 12. From the proof of Theorems 1 and 6, we see that we have not relied on specific properties; Theorems 1 and 6 hold for any positive definite matrix .In the next section, we will discuss the choices of .

Choices of the Weight Matrix 𝑊.
In this section we will discuss the choice of the weight matrix .Given a weight matrix, an efficient method of factoring or iteratively solving systems with the preconditioner must be sought.These considerations are motivated by the fact that each iteration of a preconditioned Krylov subspace method requires solutions to linear systems of the forms P   =  and F   = ; based on the block structure of P  and F  , this requires solving systems with  +    −1  and .
If  is diagonal or block diagonal with small blocks, the matrix  +    −1  is also going to be sparse.A simple, oneparameter choice is a scaled identity.Letting  = (1/) with  = ‖‖/‖‖ 2 ,  could be chosen so that the augmenting term    is of norm comparable to .See, for example, [18,20] for a general algebraic discussion.For solving P   =  and F   = , iterative method is possible.In an iterative scheme the inner iteration can be solved using the preconditioned GMRES method based on incomplete Cholesky factorization.
In fact, for  = (1/), it follows from two identities that the action of the preconditioner on a given vector requires one application of ( +   ) −1 and one sparse matrix-vector product with   .Clearly, the main issue is how to solve linear systems with coefficient matrix  +   .
For large problems these have to be solved by an inner iterative method.Even though the inner solvers need not be performed to high accuracy, developing a robust and efficient iterative method for such problems is a nontrivial task.An effective multigrid method has been developed in [12].We will further discuss the issue of inexact solvers in the section on numerical experiments.Finally, we mention that other choices of  are possible.For example, setting  = (  ) −1 has the advantage that the matrix (  ) −1   is an orthogonal projector onto the range of , which is orthogonal to the null space of   , and since the null space of  does not intersect with the null space of   either, such a choice of  is viable.See [11,18,20] for other choices of .

Numerical Experiments
In this section, we illustrate the performance of our preconditioning approach on Maxwell equation in which the (1, 1) blocks of the associated matrices are highly singular.Our results include iteration counts, condition counts, and some timings and eigenvalue plots.All the numerical experiments were performed with MATLAB 7.0.The machine we have used is a PC-AMD, CPU T7400 2.2 GHz process.We first consider a finite element discretization of the time-harmonic Maxwell equations (2) with wave numbers in a square domain (−1 ≤  ≤ 1, −1 ≤  ≤ 1) [17] and derive the saddle point linear system of the form (2).
In our first numerical examples the matrix  in the augmentation block preconditioners is taken as  = (1/)  , whereas the positive parameter  is taken as (cf.[18]) Figures 1, 2, 3, and 4 display the eigenvalues of the preconditioned matrices P −1  A and F −1  A for different values of .As predicted by the theories and Corollaries 2-5 and Corollaries 7-10, from these figures we can see that the higher the nullity of the (1, 1) block is, the stronger the eigenvalues of the preconditioned matrices are clustered.They are extremely close to 1.Note that the clustering of the spectrum away from the zero eigenvalue suggests that GMRES with the block triangular preconditioners will converge fast.3.1.Exact Solvers.In this section we first study the effectiveness of the preconditioner with "exact" solvers: that is, linear systems with coefficient matrix  +    are solved by a direct sparse LU factorization in combination with appropriate sparsity-preserving orderings.
In Tables 1-4, Dcond, PreDcond, and PreDcond1 denote the condition numbers of the system matrices A, P −1  A, and F −1  A, respectively.Iter1, Iter2, and Iter3 denote the iteration counts for GMRES of A, P −1  A, and F −1  A, respectively.As is evident, our solver performs extremely well.From Tables 1 and 2, we can see that two types of iteration numbers with the preconditioners P −1  and F −1  are slightly changed by the change of parameter .From Table 1, augmentation block triangular preconditioners P −  ( > 1/2) show more effective performance than those of augmentation block triangular preconditioners P +  (0 <  < 1/2).From Tables 3 and 4, we can see that the preconditioners P −1  and F −1  greatly reduce the condition counts of coefficient matrix A, and iteration numbers of the preconditioned GMRES with the preconditioners P −1  and F −1  are greatly reduced by the change of grids.Moreover, from Tables 3 and 4, we can know that the preconditioned GMRES method with  the preconditioner F −1  is more effective than the preconditioned GMRES method with the preconditioner P −1  .

Inexact Solvers.
In practice, using exact solvers in the application of the preconditioner may be too expensive.Here we consider replacing the exact solvers with inexact ones, obtained via an inner preconditioned GMRES iteration.In this paper we explore the use of drop tolerance-based incomplete LU as the preconditioner for the inner iteration.During implementation of our augmentation block preconditioners, we need the operation ( +   ) −1  for a given vector  or, equivalently, to solve the following equation: ( +   )V =  for which we use an incomplete LU factorization of  +    = + with drop tolerance , where  = ( , ), | , | ≤ .

Conclusions
In this paper, we have considered preconditioned iterative methods applied to discretizations of the time-harmonic Maxwell equations ( = 0) in electromagnetic problems.
We have analyzed the spectral properties as well as the computational performance of two types of block triangular augmentation preconditioners.Complete theoretical analysis shows that all eigenvalues of the preconditioned matrices are strongly clustered.A good parameter choice may substantially reduce the iteration numbers and condition counts.Especially, we have shown that, in cases where the (1, 1) block has high nullity, convergence for each of the two preconditioned GMRES iterative methods is guaranteed to be almost immediate.We have compared the performance of various types of preconditioners with regard to the grids, condition counts, the time step, and other problems parameters.
Finally, the approach that we have investigated is parameter dependent, and it would be desirable to explore choices that reduce the overall computational cost.

Table 1 :
Comparison of iteration counts for GMRES of A and P −1  A for different .

Table 2 :
Comparison of iteration counts for GMRES of A and F −1  A for different .