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

Authors: Qingbing Liu
Journal: ISRN Applied Mathematics
 (2013)

---

## Abstract

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 point matrices, involving the choice of the parameter. Meanwhile, we provide some results of numerical experiments to show the effectiveness of the proposed parameterized preconditioners.


---

## Body

## 1. Introduction

The following time-harmonic Maxwell equations in lossless media and perfectly conduction boundaries [1–5] are partial differential equations of the mixed form with constant coefficients: find the vector field u and the multiplier p such that
(1)∇×∇×u+∇p=finΩ,∇·u=0inΩ,u×n=0on∂Ω,p=0on∂Ω.Here,Ω⊂ℝ3 is a simply connected polyhedron domain with a connected boundary ∂Ω, and n denotes the outward unit normal on ∂Ω. The datum f is a given generic source (not necessarily divergence free), and the wave number satisfies k2=ω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
(2)𝒜x=[ABTB0][up]=[g0]=b,
where u∈Rn and p∈Rm are finite arrays representing the finite element approximations, and g∈Rn is the load vector associated with f. The matrix A∈Rn×n is symmetric positive semidefinite with nullity s and corresponds to the discrete curl-curl operator; B∈Rm×n 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–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
(3)𝒦j(A,r0)=span{r0,Ar0,A2r0,…,Aj-1r0},
where 𝒦j(A,r0) is the jth Krylov subspace associated with A and r0 (see, e.g., [9, 11, 12]).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 theA-norm. However, the algorithm mainly performs well if the matrix A 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 A, 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–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(m), can be utilized. Since restarting removes the previous convergence history, GMRES(m) is not guaranteed to converge. There is no specific rule to determine the restart parameter m. In cases characterized by superlinear convergence, m 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–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 key for the rapid convergence of an iterative method for a linear system of the form 𝒜x=b is the availability of an effective preconditioner 𝒫. Each step of an outer iteration for solving the preconditioned linear system 𝒫-1𝒜x=𝒫-1b requires the solution of an inner linear system whose coefficient matrix is 𝒜. Therefore, convergence of the outer iteration is fast if the eigenvalues of the preconditioned matrix 𝒫-1𝒜 are clustered, but careful attention must be paid to the conditioning and eigenvalue distribution of the matrix 𝒜 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
(4)ℳ=[A+BTW-1B00W],
where W∈ℝm×m 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
(5)𝒩=[A+BTW-1BBT0W]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
(6)𝒯+=[A+BTW-1BBT0W],𝒯-=[A+BTW-1BBT0-W].
And he has shown that if the nullity of A is m, then the preconditioned matrices 𝒯+-1𝒜 and 𝒯--1𝒜 have only three distinct eigenvalues 1, (-1±5)/2 and 1, (1±i3)/2, respectively. Thus, the preconditioned GMRES with either augmentation block triangular preconditioner converges within three iterations.Based on the preconditioners𝒯+ and 𝒯-, 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 time-harmonic 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 Section2, 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.


## 2. 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𝒫k-1𝒜 and ℱk-1𝒜, assuming exact solves for the (1,1) block A. 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 A is m.


### 2.1. Spectral Properties of the Preconditioned Matrix

We let(7)𝒜=[ABTB0]
represent the saddle point matrices of (10). We assume that A is symmetric and positive semidefinite with nullity s and that B is of size m×n (m≤n) and has full row rank. Note that the assumption that 𝒜 is nonsingular implies that null(A)∩null(B)=0 and null(A)∪null(B)=Rn, which we use in our analysis below.We present the following block triangular preconditioner:(8)𝒫k=k𝒯-+(1-k)𝒯+=[A+BTW-1BBT0(1-2k)W],k≠12.
We will study spectral properties and computational performance of the preconditioner 𝒫k. Also numerical experiments show that the preconditioner 𝒫k with a parameter is more efficient than both 𝒯+ and 𝒯- 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𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of 𝒫k-1𝒜 of geometric multiplicity n-m, and λ=(-1±5-8k)/2(1-2k) is an eigenvalue of geometric multiplicity s. The remaining 2(m-s) eigenvalues satisfy the relation
(9)λ=(1-μ-2k)±4μ2-8kμ2-4k+2μ+12(μ+1)(1-2k),
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(10)BTW-1Bx=μAx.Proof.
Suppose thatλ is an eigenvalue of 𝒫k-1𝒜, whose eigenvalue is (xT,yT)T. Then we have
(11)𝒫k-1𝒜[xy]=λ[xy].
Furthermore, it satisfies the generalized eigenvalue problem
(12)[ABTB0][xy]=λ[A+BTW-1BBT0(1-2k)W][xy]
or
(13)Ax+BTy=λ(A+BTW-1B)x+λBTy,Bx=λ(1-2k)Wy.
As 𝒜 is nonsingular, λ≠0. The second equality gives y=(1/λ(1-2k))W-1Bx, and substituting it into the first equality gives
(14)(λ2-λ)(1-2k)Ax+[λ2(1-2k)+λ-1]BTW-1Bx=0.
If x∈null(B), then (20) implies that
(15)(λ2-λ)(1-2k)Ax=0
from which it follows that λ=1 is an eigenvalue of 𝒫k-1𝒜 of geometric multiplicity n-m.
Ifx∈null(A), then (14) implies that
(16)[λ2(1-2k)+λ-1]BTW-1Bx=0
from which we obtain that λ=(-1+5-8k)/2(1-2k) are two eigenvalues of 𝒫k-1𝒜 of geometric multiplicity s.
We have determinedn-m+2s eigenvalues. Now we consider the remaining 2(m-s) eigenvalues of 𝒫k-1𝒜.
Suppose thatλ≠(-1+5-8k)/2(1-2k); then from (14) we have
(17)BTW-1Bx=μAx,
where μ=-(λ2-λ)(1-2k)/(λ2(1-2k)+λ-1), which implies that
(18)λ=--(2k+μ-1)±4μ2-8kμ2-4k+2μ+12(μ+1)(1-2k).
Since BTW-1B+A is assumed to be nonsingular, the matrix pencil BTW-1B+εA is regular (cf. [11]). Thus, the generalized eigenvalue problem (17) is well posed. It is easy to see that μ=0 is a generalized eigenvalue of geometric multiplicity n-m and μ=∞ is a generalized eigenvalue of geometric multiplicity s. The remaining m-s generalized eigenvalues μ of (17) determine the remaining 2(m-s) eigenvalues λ of 𝒫k-1𝒜 by (18).Ifk=0 and 𝒫k apparently reduce to 𝒯+, then we have the following corollary.Corollary 2.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of 𝒯+-1𝒜 of geometric multiplicity n-m, and λ=(-1±5)/2 is an eigenvalue of geometric multiplicity s. The remaining 2(m-s) eigenvalues satisfy the relation
(19)λ=(1-μ)±4μ2+2μ+12(μ+1),
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(20)BTW-1Bx=μAx.Corollary 3.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity m. Then the preconditioned matrix 𝒯+-1𝒜 has precisely three distinct eigenvalues: λ=1 of geometric multiplicity n-m and λ=(-1+5)/2 and λ=(-1-5)/2 both of multiplicity m.Ifk=1 and 𝒫k apparently reduce to 𝒯-, then we have the following corollary.Corollary 4.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of 𝒯--1𝒜 of geometric multiplicity n-m, and λ=(1±i3)/2 is an eigenvalue of geometric multiplicity s. The remaining 2(m-s) eigenvalues satisfy the relation
(21)λ=(μ+1)±-4μ2+2μ-32(μ+1),
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(22)BTW-1Bx=μAx.Corollary 5.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity m. Then the preconditioned matrix 𝒯--1𝒜 has precisely three distinct eigenvalues: λ=1 of geometric multiplicity n-m and λ=(1+i3)/2 and λ=(1-i3)/2 both of multiplicity m.We next consider another block triangular preconditioner(23)ℱk=[A+BTW-1B2kBT0(1-2k)W],k≠12.
We give the following result to describe the spectral distribution of the preconditioned matrix ℱk-1𝒜.Theorem 6.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of ℱk-1𝒜 of geometric multiplicity n, and λ=1/(2k-1) is an eigenvalue of geometric multiplicity s. The remaining m-s eigenvalues satisfy the relation
(24)λ=-μ(μ+1)(1-2k),
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(25)BTW-1Bx=μAx.Proof.
Suppose thatλ is an eigenvalue of ℱk-1𝒜, whose eigenvalue is (xT,yT)T. Then we have
(26)ℱk-1𝒜[xy]=λ[xy].
Furthermore, it satisfies the generalized eigenvalue problem
(27)[ABTB0][xy]=λ[A+BTW-1B2kBT0(1-2k)W][xy]
or(28)Ax+BTy=λ(A+BTW-1B)x+2λkBTy,Bx=λ(1-2k)Wy.
As 𝒜 is nonsingular, λ≠0. The second equality gives y=(1/λ(1-2k))W-1Bx, and substituting it into the first equality gives
(29)(λ-1)[λ(1-2k)Ax+(λ+1-2kλ)BTW-1Bx]=0.
It is straightforward to see that any vector x∈Rn satisfies (29) with λ=1, and thus λ=1 is an eigenvalue of ℱk-1𝒜. We claim that the eigenvalue λ=1 has geometric multiplicity n.
Ifλ≠1, then from (29) we have
(30)λ(1-2k)Ax+(λ+1-2kλ)BTW-1Bx=0.
Ifx∈null(A), then from (30) we have
(31)(λ+1-2kλ)BTW-1Bx=0.
From (31), we must have λ+1-2kλ=0. Since A is singular with nullity s, hence, we know that λ=1/(2k-1) is an eigenvalue of ℱk-1𝒜 of geometric multiplicity s.
We have determinedn+s eigenvalues. Now we consider the remaining m-s eigenvalues of ℱk-1𝒜.
Suppose thatλ≠1,1/(2k-1); then from (30) we have
(32)BTW-1Bx=μAx,
where μ=-λ(1-2k)/(λ(1-2k)+1), which implies that λ=-μ/(μ+1)(1-2k). Since BTW-1B+A is assumed to be nonsingular, the matrix pencil BTW-1B+εA is regular (cf.  [11]). Thus, the generalized eigenvalue problem (32) is well posed. It is easy to see that μ=0 is a generalized eigenvalue of geometric multiplicity n-m, and μ=∞ is a generalized eigenvalue of geometric multiplicity s. The remaining m-s generalized eigenvalues μ of (32) determine the remaining m-s eigenvalues λ of ℱk-1𝒜 by -μ/(μ+1)(1-2k).Ifk=0 and 𝒫k apparently reduce to ℳ, then we have the following corollary.Corollary 7.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of ℳ-1𝒜 of geometric multiplicity n, and λ=-1 is an eigenvalue of geometric multiplicity s. The remaining m-s eigenvalues satisfy the relation
(33)λ=-μμ+1,
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(34)BTW-1Bx=μAx.Corollary 8.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity m. Then λ=1 is an eigenvalue of ℳ-1𝒜 of geometric multiplicity n, and λ=-1 is an eigenvalue of geometric multiplicity m.Ifk=1, then we have the following corollary.Corollary 9.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of ℱ1-1𝒜 of geometric multiplicity n+s. The remaining m-s eigenvalues satisfy the relation
(35)λ=μμ+1,
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(36)BTW-1Bx=μAx.Corollary 10.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity m. Then the preconditioned matrix ℱ1-1𝒜 has only precisely one eigenvalue: λ=1 of geometric multiplicity n+m.Remark 11.
From Corollaries8 and 10, if the nullity of A is m, it is readily seen that the preconditioned matrix ℳ-1𝒜 has precisely two distinct eigenvalues and that the preconditioned matrix ℱ1-1𝒜 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 Theorems1 and 6, we see that we have not relied on specific properties; Theorems 1 and 6 hold for any positive definite matrix W. In the next section, we will discuss the choices of W.


### 2.2. Choices of the Weight MatrixW

In this section we will discuss the choice of the weight matrixW. 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 𝒫kx=b and ℱkx=b; based on the block structure of 𝒫k and ℱk, this requires solving systems with A+BTW-1B and W.IfW is diagonal or block diagonal with small blocks, the matrix A+BTW-1B is also going to be sparse. A simple, one-parameter choice is a scaled identity. Letting W=(1/γ)I with γ=∥A∥/∥B∥2, γ could be chosen so that the augmenting term γBTB is of norm comparable to A. See, for example, [18, 20] for a general algebraic discussion. For solving 𝒫kx=b and ℱkx=b, 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, forW=(1/γ)I, it follows from two identities
(37)𝒫k-1=((A+γBTB)-100Im)(InBT0-Im)(In00γ2k-1Im),ℱk-1=((A+γBTB)-100Im)(In2kBT0-Im)(In00γ2k-1Im)
that the action of the preconditioner on a given vector requires one application of (A+γBTB)-1 and one sparse matrix-vector product with BT. Clearly, the main issue is how to solve linear systems with coefficient matrix A+γBTB. 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 ofW are possible. For example, setting W=(BTB)-1 has the advantage that the matrix B(BTB)-1BT is an orthogonal projector onto the range of B, which is orthogonal to the null space of BT, and since the null space of A does not intersect with the null space of BT either, such a choice of W is viable. See [11, 18, 20] for other choices of W.


## 2.1. Spectral Properties of the Preconditioned Matrix

We let(7)𝒜=[ABTB0]
represent the saddle point matrices of (10). We assume that A is symmetric and positive semidefinite with nullity s and that B is of size m×n (m≤n) and has full row rank. Note that the assumption that 𝒜 is nonsingular implies that null(A)∩null(B)=0 and null(A)∪null(B)=Rn, which we use in our analysis below.We present the following block triangular preconditioner:(8)𝒫k=k𝒯-+(1-k)𝒯+=[A+BTW-1BBT0(1-2k)W],k≠12.
We will study spectral properties and computational performance of the preconditioner 𝒫k. Also numerical experiments show that the preconditioner 𝒫k with a parameter is more efficient than both 𝒯+ and 𝒯- 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𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of 𝒫k-1𝒜 of geometric multiplicity n-m, and λ=(-1±5-8k)/2(1-2k) is an eigenvalue of geometric multiplicity s. The remaining 2(m-s) eigenvalues satisfy the relation
(9)λ=(1-μ-2k)±4μ2-8kμ2-4k+2μ+12(μ+1)(1-2k),
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(10)BTW-1Bx=μAx.Proof.
Suppose thatλ is an eigenvalue of 𝒫k-1𝒜, whose eigenvalue is (xT,yT)T. Then we have
(11)𝒫k-1𝒜[xy]=λ[xy].
Furthermore, it satisfies the generalized eigenvalue problem
(12)[ABTB0][xy]=λ[A+BTW-1BBT0(1-2k)W][xy]
or
(13)Ax+BTy=λ(A+BTW-1B)x+λBTy,Bx=λ(1-2k)Wy.
As 𝒜 is nonsingular, λ≠0. The second equality gives y=(1/λ(1-2k))W-1Bx, and substituting it into the first equality gives
(14)(λ2-λ)(1-2k)Ax+[λ2(1-2k)+λ-1]BTW-1Bx=0.
If x∈null(B), then (20) implies that
(15)(λ2-λ)(1-2k)Ax=0
from which it follows that λ=1 is an eigenvalue of 𝒫k-1𝒜 of geometric multiplicity n-m.
Ifx∈null(A), then (14) implies that
(16)[λ2(1-2k)+λ-1]BTW-1Bx=0
from which we obtain that λ=(-1+5-8k)/2(1-2k) are two eigenvalues of 𝒫k-1𝒜 of geometric multiplicity s.
We have determinedn-m+2s eigenvalues. Now we consider the remaining 2(m-s) eigenvalues of 𝒫k-1𝒜.
Suppose thatλ≠(-1+5-8k)/2(1-2k); then from (14) we have
(17)BTW-1Bx=μAx,
where μ=-(λ2-λ)(1-2k)/(λ2(1-2k)+λ-1), which implies that
(18)λ=--(2k+μ-1)±4μ2-8kμ2-4k+2μ+12(μ+1)(1-2k).
Since BTW-1B+A is assumed to be nonsingular, the matrix pencil BTW-1B+εA is regular (cf. [11]). Thus, the generalized eigenvalue problem (17) is well posed. It is easy to see that μ=0 is a generalized eigenvalue of geometric multiplicity n-m and μ=∞ is a generalized eigenvalue of geometric multiplicity s. The remaining m-s generalized eigenvalues μ of (17) determine the remaining 2(m-s) eigenvalues λ of 𝒫k-1𝒜 by (18).Ifk=0 and 𝒫k apparently reduce to 𝒯+, then we have the following corollary.Corollary 2.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of 𝒯+-1𝒜 of geometric multiplicity n-m, and λ=(-1±5)/2 is an eigenvalue of geometric multiplicity s. The remaining 2(m-s) eigenvalues satisfy the relation
(19)λ=(1-μ)±4μ2+2μ+12(μ+1),
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(20)BTW-1Bx=μAx.Corollary 3.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity m. Then the preconditioned matrix 𝒯+-1𝒜 has precisely three distinct eigenvalues: λ=1 of geometric multiplicity n-m and λ=(-1+5)/2 and λ=(-1-5)/2 both of multiplicity m.Ifk=1 and 𝒫k apparently reduce to 𝒯-, then we have the following corollary.Corollary 4.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of 𝒯--1𝒜 of geometric multiplicity n-m, and λ=(1±i3)/2 is an eigenvalue of geometric multiplicity s. The remaining 2(m-s) eigenvalues satisfy the relation
(21)λ=(μ+1)±-4μ2+2μ-32(μ+1),
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(22)BTW-1Bx=μAx.Corollary 5.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity m. Then the preconditioned matrix 𝒯--1𝒜 has precisely three distinct eigenvalues: λ=1 of geometric multiplicity n-m and λ=(1+i3)/2 and λ=(1-i3)/2 both of multiplicity m.We next consider another block triangular preconditioner(23)ℱk=[A+BTW-1B2kBT0(1-2k)W],k≠12.
We give the following result to describe the spectral distribution of the preconditioned matrix ℱk-1𝒜.Theorem 6.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of ℱk-1𝒜 of geometric multiplicity n, and λ=1/(2k-1) is an eigenvalue of geometric multiplicity s. The remaining m-s eigenvalues satisfy the relation
(24)λ=-μ(μ+1)(1-2k),
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(25)BTW-1Bx=μAx.Proof.
Suppose thatλ is an eigenvalue of ℱk-1𝒜, whose eigenvalue is (xT,yT)T. Then we have
(26)ℱk-1𝒜[xy]=λ[xy].
Furthermore, it satisfies the generalized eigenvalue problem
(27)[ABTB0][xy]=λ[A+BTW-1B2kBT0(1-2k)W][xy]
or(28)Ax+BTy=λ(A+BTW-1B)x+2λkBTy,Bx=λ(1-2k)Wy.
As 𝒜 is nonsingular, λ≠0. The second equality gives y=(1/λ(1-2k))W-1Bx, and substituting it into the first equality gives
(29)(λ-1)[λ(1-2k)Ax+(λ+1-2kλ)BTW-1Bx]=0.
It is straightforward to see that any vector x∈Rn satisfies (29) with λ=1, and thus λ=1 is an eigenvalue of ℱk-1𝒜. We claim that the eigenvalue λ=1 has geometric multiplicity n.
Ifλ≠1, then from (29) we have
(30)λ(1-2k)Ax+(λ+1-2kλ)BTW-1Bx=0.
Ifx∈null(A), then from (30) we have
(31)(λ+1-2kλ)BTW-1Bx=0.
From (31), we must have λ+1-2kλ=0. Since A is singular with nullity s, hence, we know that λ=1/(2k-1) is an eigenvalue of ℱk-1𝒜 of geometric multiplicity s.
We have determinedn+s eigenvalues. Now we consider the remaining m-s eigenvalues of ℱk-1𝒜.
Suppose thatλ≠1,1/(2k-1); then from (30) we have
(32)BTW-1Bx=μAx,
where μ=-λ(1-2k)/(λ(1-2k)+1), which implies that λ=-μ/(μ+1)(1-2k). Since BTW-1B+A is assumed to be nonsingular, the matrix pencil BTW-1B+εA is regular (cf.  [11]). Thus, the generalized eigenvalue problem (32) is well posed. It is easy to see that μ=0 is a generalized eigenvalue of geometric multiplicity n-m, and μ=∞ is a generalized eigenvalue of geometric multiplicity s. The remaining m-s generalized eigenvalues μ of (32) determine the remaining m-s eigenvalues λ of ℱk-1𝒜 by -μ/(μ+1)(1-2k).Ifk=0 and 𝒫k apparently reduce to ℳ, then we have the following corollary.Corollary 7.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of ℳ-1𝒜 of geometric multiplicity n, and λ=-1 is an eigenvalue of geometric multiplicity s. The remaining m-s eigenvalues satisfy the relation
(33)λ=-μμ+1,
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(34)BTW-1Bx=μAx.Corollary 8.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity m. Then λ=1 is an eigenvalue of ℳ-1𝒜 of geometric multiplicity n, and λ=-1 is an eigenvalue of geometric multiplicity m.Ifk=1, then we have the following corollary.Corollary 9.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity s. Then λ=1 is an eigenvalue of ℱ1-1𝒜 of geometric multiplicity n+s. The remaining m-s eigenvalues satisfy the relation
(35)λ=μμ+1,
where μ are some m-s generalized eigenvalues of the following generalized eigenvalue problem:
(36)BTW-1Bx=μAx.Corollary 10.
Let𝒜 be nonsingular, and its (1,1) block A is singular with nullity m. Then the preconditioned matrix ℱ1-1𝒜 has only precisely one eigenvalue: λ=1 of geometric multiplicity n+m.Remark 11.
From Corollaries8 and 10, if the nullity of A is m, it is readily seen that the preconditioned matrix ℳ-1𝒜 has precisely two distinct eigenvalues and that the preconditioned matrix ℱ1-1𝒜 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 Theorems1 and 6, we see that we have not relied on specific properties; Theorems 1 and 6 hold for any positive definite matrix W. In the next section, we will discuss the choices of W.


## 2.2. Choices of the Weight MatrixW

In this section we will discuss the choice of the weight matrixW. 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 𝒫kx=b and ℱkx=b; based on the block structure of 𝒫k and ℱk, this requires solving systems with A+BTW-1B and W.IfW is diagonal or block diagonal with small blocks, the matrix A+BTW-1B is also going to be sparse. A simple, one-parameter choice is a scaled identity. Letting W=(1/γ)I with γ=∥A∥/∥B∥2, γ could be chosen so that the augmenting term γBTB is of norm comparable to A. See, for example, [18, 20] for a general algebraic discussion. For solving 𝒫kx=b and ℱkx=b, 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, forW=(1/γ)I, it follows from two identities
(37)𝒫k-1=((A+γBTB)-100Im)(InBT0-Im)(In00γ2k-1Im),ℱk-1=((A+γBTB)-100Im)(In2kBT0-Im)(In00γ2k-1Im)
that the action of the preconditioner on a given vector requires one application of (A+γBTB)-1 and one sparse matrix-vector product with BT. Clearly, the main issue is how to solve linear systems with coefficient matrix A+γBTB. 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 ofW are possible. For example, setting W=(BTB)-1 has the advantage that the matrix B(BTB)-1BT is an orthogonal projector onto the range of B, which is orthogonal to the null space of BT, and since the null space of A does not intersect with the null space of BT either, such a choice of W is viable. See [11, 18, 20] for other choices of W.


## 3. 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≤x≤1,-1≤y≤1) [17] and derive the saddle point linear system of the form (2).In our first numerical examples the matrixW in the augmentation block preconditioners is taken as W=(1/γ)Im, whereas the positive parameter γ is taken as (cf. [18])
(38)γ=∥A∥1∥B∥12.Figures1, 2, 3, and 4 display the eigenvalues of the preconditioned matrices 𝒫k-1𝒜 and ℱk-1𝒜 for different values of k. 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.Figure 1
Eigenvalues of the preconditioned matrix𝒫k-1𝒜 with k=0.4.Figure 2
Eigenvalues of the preconditioned matrix𝒫k-1𝒜 with k=0.4.Figure 3
Eigenvalues of the preconditioned matrixℱk-1𝒜 with k=0.4.Figure 4
Eigenvalues of the preconditioned matrixℱk-1𝒜 with k=0.4.


### 3.1. Exact Solvers

In this section we first study the effectiveness of the preconditioner with “exact” solvers: that is, linear systems with coefficient matrixA+γBTB are solved by a direct sparse LU factorization in combination with appropriate sparsity-preserving orderings.In Tables1–4, Dcond, PreDcond, and PreDcond1 denote the condition numbers of the system matrices 𝒜, 𝒫k-1𝒜, and ℱk-1𝒜, respectively. Iter1, Iter2, and Iter3 denote the iteration counts for GMRES of 𝒜, 𝒫k-1𝒜, and ℱk-1𝒜, 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 𝒫k-1 and ℱk-1 are slightly changed by the change of parameter k. From Table 1, augmentation block triangular preconditioners 𝒫k- (k>1/2) show more effective performance than those of augmentation block triangular preconditioners 𝒫k+ (0<k<1/2). From Tables 3 and 4, we can see that the preconditioners 𝒫k-1 and ℱk-1 greatly reduce the condition counts of coefficient matrix 𝒜, and iteration numbers of the preconditioned GMRES with the preconditioners 𝒫k-1 and ℱk-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 ℱk-1 is more effective than the preconditioned GMRES method with the preconditioner 𝒫k-1.Table 1
Comparison of iteration counts for GMRES of𝒜 and 𝒫k-1𝒜 for different k.
k
0
0.3
0.9
1
Iter1
10
10
10
10
Iter2
4
3
2
2
Dcond
9.3418
e
+
016
9.3418
e
+
016
9.3418
e
+
016
9.3418
e
+
016
PreDcond
7.8317
e
+
05
1.0174
e
+
06
7.2132
e
+
05
3.3704
e
+
05Table 2
Comparison of iteration counts for GMRES of𝒜 and ℱk-1𝒜 for different k.
k
0
0.3
0.9
1
Iter1
10
10
10
10
Iter3
1
2
2
1
Dcond
9.3418
e
+
016
9.3418
e
+
016
9.3418
e
+
016
9.3418
e
+
016
PreDcond
3.1740
e
+
05
2.2504
e
+
04
8.1205
e
+
04
1.0528
e
+
05Table 3
Iteration counts for GMRES of𝒜, 𝒫k-1𝒜, and ℱk-1𝒜 for k=0.01 and tol=1.0000e-04.
Grid
8
×
8
16
×
16
32
×
32
64
×
64
Iter1
89
835
1000
>1000
Iter2
3
4
6
8
Iter3
1
2
2
2Table 4
Iteration counts for GMRES of𝒜, 𝒫k-1𝒜, and ℱk-1𝒜 for k=0.01 and tol=1.0000e-08.
Grid
8
×
8
16
×
16
32
×
32
64
×
64
Iter1
89
835
1000
>1000
Iter2
6
11
25
32
Iter3
2
4
7
8


### 3.2. 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(A+γBTB)-1u for a given vector u or, equivalently, to solve the following equation: (A+γBTB)v=u for which we use an incomplete LU factorization of A+γBTB=LU+R with drop tolerance τ, where R=(ri,j),|ri,j|≤τ.


## 3.1. Exact Solvers

In this section we first study the effectiveness of the preconditioner with “exact” solvers: that is, linear systems with coefficient matrixA+γBTB are solved by a direct sparse LU factorization in combination with appropriate sparsity-preserving orderings.In Tables1–4, Dcond, PreDcond, and PreDcond1 denote the condition numbers of the system matrices 𝒜, 𝒫k-1𝒜, and ℱk-1𝒜, respectively. Iter1, Iter2, and Iter3 denote the iteration counts for GMRES of 𝒜, 𝒫k-1𝒜, and ℱk-1𝒜, 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 𝒫k-1 and ℱk-1 are slightly changed by the change of parameter k. From Table 1, augmentation block triangular preconditioners 𝒫k- (k>1/2) show more effective performance than those of augmentation block triangular preconditioners 𝒫k+ (0<k<1/2). From Tables 3 and 4, we can see that the preconditioners 𝒫k-1 and ℱk-1 greatly reduce the condition counts of coefficient matrix 𝒜, and iteration numbers of the preconditioned GMRES with the preconditioners 𝒫k-1 and ℱk-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 ℱk-1 is more effective than the preconditioned GMRES method with the preconditioner 𝒫k-1.Table 1
Comparison of iteration counts for GMRES of𝒜 and 𝒫k-1𝒜 for different k.
k
0
0.3
0.9
1
Iter1
10
10
10
10
Iter2
4
3
2
2
Dcond
9.3418
e
+
016
9.3418
e
+
016
9.3418
e
+
016
9.3418
e
+
016
PreDcond
7.8317
e
+
05
1.0174
e
+
06
7.2132
e
+
05
3.3704
e
+
05Table 2
Comparison of iteration counts for GMRES of𝒜 and ℱk-1𝒜 for different k.
k
0
0.3
0.9
1
Iter1
10
10
10
10
Iter3
1
2
2
1
Dcond
9.3418
e
+
016
9.3418
e
+
016
9.3418
e
+
016
9.3418
e
+
016
PreDcond
3.1740
e
+
05
2.2504
e
+
04
8.1205
e
+
04
1.0528
e
+
05Table 3
Iteration counts for GMRES of𝒜, 𝒫k-1𝒜, and ℱk-1𝒜 for k=0.01 and tol=1.0000e-04.
Grid
8
×
8
16
×
16
32
×
32
64
×
64
Iter1
89
835
1000
>1000
Iter2
3
4
6
8
Iter3
1
2
2
2Table 4
Iteration counts for GMRES of𝒜, 𝒫k-1𝒜, and ℱk-1𝒜 for k=0.01 and tol=1.0000e-08.
Grid
8
×
8
16
×
16
32
×
32
64
×
64
Iter1
89
835
1000
>1000
Iter2
6
11
25
32
Iter3
2
4
7
8


## 3.2. 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(A+γBTB)-1u for a given vector u or, equivalently, to solve the following equation: (A+γBTB)v=u for which we use an incomplete LU factorization of A+γBTB=LU+R with drop tolerance τ, where R=(ri,j),|ri,j|≤τ.


## 4. Conclusions

In this paper, we have considered preconditioned iterative methods applied to discretizations of the time-harmonic Maxwell equations(k=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.



---

*Source: 905723-2013-09-15.xml*