CI. OMPUT
Vol. 37, No. 5, pp. A2403–A2423
A ROBUST NUMERICAL ALGORITHM FOR COMPUTING
MAXWELL’S TRANSMISSION EIGENVALUE PROBLEMS∗
TSUNG-MING HUANG†, WEI-QIANG HUANG‡,AND WEN-WEI LIN‡
Abstract. We study a robust and efficient eigensolver for computing a few smallest positive
eigenvalues of the three-dimensional Maxwell’s transmission eigenvalue problem. The discretized governing equations by the N´ed´elec edge element result in a large-scale quadratic eigenvalue problem (QEP) for which the spectrum contains many zero eigenvalues and the coefficient matrices consist of patterns in the matrix form XY−1Z, both of which prevent existing eigenvalue solvers from being
ef-ficient. To remedy these difficulties, we rewrite the QEP as a particular nonlinear eigenvalue problem and develop a secant-type iteration, together with an indefinite locally optimal block preconditioned conjugate gradient (LOBPCG) method, to sequentially compute the desired positive eigenvalues. Furthermore, we propose a novel method to solve the linear systems in each iteration of LOBPCG. Intensive numerical experiments show that our proposed method is robust, although the desired real eigenvalues are surrounded by complex eigenvalues.
Key words. transmission eigenvalues, Maxwell’s equations, quadratic eigenvalue problems, secant-type iteration, LOBPCG
AMS subject classifications. 78A46, 65N30, 65N25, 65F15 DOI. 10.1137/15M1018927
1. Introduction. The transmission eigenvalue problem has recently attracted much attention in the area of inverse scattering theory, as it is important for the study of the direct/inverse scattering problem for nonabsorbing inhomogeneous me-dia [6, 8, 9, 10, 11, 12, 13, 20, 30]. As shown in [3, 4, 5, 6, 7, 8, 31], transmission eigenvalues can be determined from the far-field pattern of the scattered wave or from the near-field data, and used to estimate the material properties of the scat-tering object. In addition, transmission eigenvalues are also related to the validity of some recently developed reconstruction methods for scattering problems such as the linear sampling method and factorization method [11]. For recent progress in the theories and applications of transmission eigenvalue problems, we refer to [10] and the references therein.
Efficient numerical methods to determine transmission eigenvalues are required in estimating the index of refraction [6, 31], and numerical evidence from the discrete system may contribute to the progress of further theoretical developments such as the distribution of real eigenvalues for the original infinite dimensional system. Nonethe-less, numerical techniques for solving the transmission eigenvalues are limited and only a few papers have addressed the issues of numerical computation on this topic in the past few years, partly because the transmission eigenvalue problem is neither elliptic nor self-adjoint and as a consequence, it cannot be addressed by the standard theory of elliptic partial differential equations.
∗Submitted to the journal’s Methods and Algorithms for Scientific Computing section April 28,
2015; accepted for publication (in revised form) July 16, 2015; published electronically September 23, 2015. This work was supported by the Ministry of Science and Technology, the National Center for Theoretical Sciences, and the ST Yau Center at the National Chiao Tung University in Taiwan.
http://www.siam.org/journals/sisc/37-5/M101892.html
†Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan (min@
ntnu.edu.tw).
‡Department of Applied Mathematics and ST Yau Center, National Chiao Tung University,
Hsinchu 300, Taiwan (wqhuang@math.nctu.edu.tw, wwlin@math.nctu.edu.tw). A2403
Recently, there have been some papers [12, 15, 18, 19, 21, 25, 28, 32, 33] ad-dressing numerical computations in transmission eigenvalue problems. In [12], three finite element methods (FEMs) were proposed for solving the two-dimensional (2D) transmission eigenvalue problem. A coupled boundary element method and FEM was introduced for the interior transmission problem in [15]. Then, Sun [32] proposed two iterative methods together with convergence analysis based on the existence theory of the fourth-order reformulation for the transmission eigenvalues [9, 30]. A mixed FEM for 2D transmission eigenvalue problems was proposed in [18] and the corresponding non-Hermitian quadratic eigenvalue problem (QEP) was solved by the classical secant iteration with an adaptive Arnoldi method. In [19], Ji, Sun, and Xie used the multi-level correction method to transform the solution of the transmission problem into a series of solutions corresponding to linear boundary value problems and solved them by the multigrid method. Li et al. in [25] rewrote the QEP as a particular param-eterized generalized eigenvalue problem (GEP) for which the eigenvalue curves are arranged in a monotonic order so that the desired curves can be sequentially solved with a new secant-type iteration.
For a three-dimensional (3D) Maxwell’s transmission eigenvalue problem, two FEMs with an adaptive Arnoldi method were proposed in [28]. The resulting GEPs are large, sparse, and non-Hermitian. The numerical challenges for solving the cor-responding GEPs are (i) a few of the smallest positive eigenvalues, which may be surrounded by complex eigenvalues, are of interest; (ii) the number of zero eigenval-ues of the GEP is huge because the nullity of the discrete double curl operator equals the number of edges in the spanning tree of a finite element mesh [2]; (iii) efficient solu-tion of the associated large sparse linear system in each iterasolu-tion of the eigensolver. To tackle drawbacks (i) and (ii), in [33], a mixed FEM was applied to an equivalent quad-curl eigenvalue problem, and the resulting QEP can be solved by a classical secant iterative method by introducing a sequence of the parameterized GEPs with symmet-ric positive definite and semidefinite coefficient matsymmet-rices. However, in [33], there is no theoretical guarantee for why the desired positive transmission values would not be lost. Moreover, due to the complexity of the matrix structures, the mesh is rather coarse, and thus more efficient eigensolvers for solving the QEP and the associated parameterized GEPs are desirable for larger problems [33]. Note that, for the vector case, Kleefeld [21] presented an accurate numerical method, based on a surface integral formulation of the interior transmission problem, for solving corresponding nonlinear eigenvalue problems for many different obstacles in three dimensions. However, only constant index of refraction and smooth domains can be treated.
In this paper, we focus on the 3D Maxwell’s transmission eigenvalue problem and make the following contributions.
• We show that the QEP in [33] can be deduced from the GEP in [28] via a
suitable equivalence transformation. In fact, the QEP and GEP have the same spectrum except for nonphysical zero eigenvalues.
• Rewriting the QEP as a particular parameterized GEP with symmetric and
symmetric positive semidefinite coefficient matrices, we then use the secant-type iteration (SecTypIt) method in [25] to sequentially compute the desired positive eigenvalues.
• To efficiently solve the parameterized GEP, we introduce the locally optimal
block preconditioned conjugate gradient (LOBPCG) method [1, 22, 23] with some modification schemes to accelerate the convergence rate. Numerical results show that the convergence of LOBPCG is not affected by the huge nullity.
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2405
• To solve the linear system appearing in LOBPCG, due to the complicated
matrix formulations of the parameterized GEP, we propose a new augmented linear system so that it can be solved by the direct/iterative method for a large-size problem.
• In practice, we propose some adaptive strategies for determining initial data
and stopping tolerance. Intensive numerical experiments show that our method is robust although the desired eigenvalues are surrounded by complex eigen-values.
Throughout this paper, the notations· and ·∗ are used to represent the trans-pose and conjugate transtrans-pose of vectors or matrices, respectively. Given a real square matrix A, we write A 0 (A 0) if A is symmetric and positive definite (semidefi-nite).
We organize this paper as follows. In section 2, we review the 3D Maxwell’s transmission eigenvalue problem and two discretization schemes proposed in [28, 33]. In section 3, we introduce the SecTypIt in [25] to address a parameterized GEP of the QEP for computing a few desired positive eigenvalues of the QEP. Sections 4 and 5 focus on the LOBPCG method and its detailed implementation for the purpose of providing an efficient and robust eigensolver to address the parameterized GEP. Numerical experiments with different indices of refraction on the unit ball and the unit square are presented in section 6. Finally, we give concluding remarks in section 7.
2. The 3D Maxwell’s transmission eigenvalue problem and its dis-cretization. Let D ⊂ R3 be a bounded simply connected domain with a piecewise smooth boundary ∂D and ν denote the unit outer normal vector to ∂D. Following [14], we introduce the Hilbert spaces
H(curl, D) :=
u∈L2(D)3:∇ × u ∈L2(D)3
, H(curl2, D) := {u ∈ H(curl, D) : ∇ × u ∈ H(curl, D)},
equipped with the scalar products
(u, v)curl:= (u, v) + (∇ × u, ∇ × v),
(u, v)curl2:= (u, v) + (∇ × u, ∇ × v)curl,
respectively. Here, (·, ·) is the L2 scalar product on D. Furthermore, H0(curl, D) and
H0(curl2, D) are, respectively, two subspaces of H(curl, D) and H(curl2, D) defined by
H0(curl, D) := {u ∈ H(curl, D) : u × ν = 0 on ∂D} ,
H0(curl2, D) := {u ∈ H0(curl, D) : ∇ × u ∈ H0(curl, D)} .
Assuming that N, N−1, and either (N − I)−1or (I − N )−1 are bounded positive definite real matrix fields on D, then, in terms of the electric field, the so-called transmission eigenvalue problem for the Maxwell’s equations is to find 0= λ ∈ C and nontrivial fields E, E0∈L2(D)3with E− E0∈ H0(curl2, D) satisfying
∇ × ∇ × E − λNE = 0 in D, (2.1a) ∇ × ∇ × E0− λE0= 0 in D, (2.1b) E× ν = E0× ν on ∂D, (2.1c) (∇ × E) × ν = (∇ × E0)× ν on ∂D. (2.1d)
The nonzero (complex) values λ such that (2.1) has nontrivial solutions E and E0are
called Maxwell’s transmission eigenvalues.
Multiplying (2.1a) and (2.1b) by suitable test functions and applying the inte-gration by parts, a variational formulation for (2.1) can be stated as follows: Find 0= λ ∈ C and E0, E ∈ H(curl, D) such that
(∇ × E, ∇ × φ) − λ(NE, φ) = 0, (2.2a) (∇ × E0, ∇ × φ) − λ(E0, φ) = 0, (2.2b) (∇ × (E − E0), ∇ × ψ) − λ(N E − E0, ψ) = 0 (2.2c)
for all φ ∈ H0(curl, D) and ψ ∈ H(curl, D) with the essential boundary condition
E× ν = E0× ν on ∂D [28]. Note that in (2.2c), the boundary condition (2.1d) has
been enforced weakly.
On the other hand, as shown in [9], (2.1) is equivalent to a quad-curl problem for E− E0∈ H0(curl2, D) satisfying
(2.3) (∇ × ∇ × −λN) (N − I)−1(∇ × ∇ × −λ) (E − E0) = 0.
A variational form of (2.3) is to find a 0 = λ ∈ C and a nontrivial field u ∈
H0(curl2, D) satisfying
(2.4) ((N −I)−1(∇×∇×u−λu), (∇×∇×φ−λφ))+λ2(u, φ)−λ(∇×u, ∇×φ) = 0 for all φ ∈ H0(curl2, D). Following the approach of a mixed formulation proposed
in [33], (2.4) can be further transformed into another weak formulation for finding 0= λ ∈ C, p ∈ H0(curl, D), and v∈ H(curl, D) such that
(∇ × v, ∇ × φ) − λ(v, φ) + λ2(p, φ) = λ(∇ × p, ∇ × φ), (2.5a)
(∇ × p, ∇ × ξ) − λ(p, ξ) = ((N − I)v, ξ) (2.5b)
for allφ ∈ H0(curl, D) and ξ ∈ H(curl, D).
Now, we use the lowest order curl-conforming N´ed´elec edge elements [27, 29] to discretize (2.1) and (2.3). Given a regular tetrahedral mesh of D, we define the space
Sh and the subspace Sh0 of Sh as
Sh={the lowest order edge elements on D} ⊂ H(curl, D), Sh0= Sh∩ H0(curl, D) ⊂ H0(curl, D)
={the functions in Shthat have vanishing DoFs on ∂D} ,
where DoFs are the degrees of freedom. Let {φ1, . . . , φn} be a basis of Sh0 and {φ1, . . . , φn, ψ1, . . . , ψm} a basis for Sh. In addition, we define ShB = span{ψj}mj=1.
Then, the mass and stiffness matrices based on linear edge elements are given by
(2.6) K = K E E H , M1= M1 F1 F1 G1 , MN = MN FN FN GN ,
where the block matrix entries are given in Table 1. Moreover, we let
S := K E, T1:= M1 F1, (2.7) M = M F F G := MN − M1 FN − F1 FN− F1 GN − G1 =MN − M1. (2.8)
Note that dim(Null(S)) > 0 as the matrices K and E are assembled from the discretization of the degenerate curl operators. Here, Null(S) denotes the null space of the matrixS. Moreover,M 0, M 0, and G 0 because of the positivity of
N and (N − I)−1.
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2407
Table 1
Stiffness and mass matrices.
Matrix Dimension Definition
interior space stiffness matrix.
K n × n
Kij= (∇ × φi, ∇ × φj)
E n × m interior/boundary stiffness matrix. Eij= (∇ × φi, ∇ × ψj)
boundary space stiffness matrix.
H m × m
Hij= (∇ × ψi, ∇ × ψj)
interior space mass matrices.
M1, MN n × n
(M1)ij= (φi, φj), (MN)ij= (N φi, φj)
interior/boundary mass matrices.
F1, FN n × m (F
1)ij= (φi, ψj), (FN)ij= (N φi, ψj)
boundary space mass matrices.
G1, GN m × m (G
1)ij= (ψi, ψj), (GN)ij= (N ψi, ψj)
2.1. The resulting generalized eigenvalue problem from (2.2). Based on the N´ed´elec edge elements, we let u0,h= ni=1uiφi ∈ Sh0, v0,h= ni=1viφi∈ S0h, and uB,h= mi=1wiψi∈ SB
h so that uh= u0,h+uB,hand vh= v0,h+uB,hare the discrete
approximations for E and E0, respectively. In addition, we set u = [u1, . . . , un],
v = [v1, . . . , vn], and w = [w1, . . . , wm] and, then, the discretization of (2.2) gives
rise to a GEP (2.9) L(λ)z := ⎛ ⎝ ⎡ ⎣K0 K0 EE E −E 0 ⎤ ⎦ − λ ⎡ ⎣M0N M01 FFN1 FN −F1 GN − G1 ⎤ ⎦ ⎞ ⎠ ⎡ ⎣uv w ⎤ ⎦ = 0. 2.2. The resulting quadratic eigenvalue problem from (2.5). Let ph = n
i=1piφi and vh =
n
j=1vjφj +
m
j=1wjψj. Moreover, we set the vectors p =
[p1, . . . , pn] and v = [v1, . . . , vn, w1, . . . , wm]. Then, with the notation in (2.6),
(2.8), and Table 1, the matrix problem corresponding to (2.5) is given by
Sv − λT1v + λ2M1p = λKp, (2.10a)
Sp− λT
1 p =Mv,
(2.10b)
whereS and T1are the matrices given in (2.7). Expressingv in terms of p by (2.10b)
and plugging it into (2.10a), we end up with the QEP
(2.11) λ2M1+ (S − λT1)M−1(S − λT1)p = λKp.
2.3. Relation between GEP (2.9) and QEP (2.11). In this subsection, we first present explicit representations for the coefficient matrices of (2.11). Then, we show that (2.11) can be deduced from the GEP (2.9) via a suitable equivalence transformation.
To make the following discussion more concise, we first introduce some convenient notation. Let
M1:= M1− F1G−1F, M := M − F G −1F, K := K − EG −1F,
where M , F , and G are defined as in (2.8). Note that M is symmetric positive definite
becauseM 0 and G 0.
Lemma 2.1. The QEP (2.11) can be expressed as (2.12a) Q(λ)p :=λ2A2+ λA1+ A0p = 0,
where A2, A1, and A0 are all n × n symmetric matrices given by A2= M1+T1M−1T1 (2.12b) = M1+ M1M−1M1+ F1G−1F1, A1=−K − SM−1T1− T1M−1S (2.12c) =−K − K M−1M1− M1M−1K− EG−1F1− F1G−1E, A0=SM−1S (2.12d) = K M−1K+ EG−1E.
In particular, A2 is positive definite and A0 is positive semidefinite. Proof. Rewriting (2.11) as
(2.13) [λ2M1+T1M−1T1
+λ−K − SM−1T1− T1M−1S+SM−1S]p = 0
and using the fact that
M−1= M F F G −1 = M−1 0 −G−1FM−1 G−1 I −F G−1 0 I ,
we can show, by routine calculation, that the coefficient matrices in (2.13) are equal to those of (2.12) (see also section 2.2 in [25] for the related results). Moreover, it is obvious to see that the coefficient matrices of (2.12) are symmetric. In addition,
A2 and A0 are positive definite and positive semidefinite, respectively, which follows
from the fact that M1 0, M 0, and dim(Null(S)) > 0.
Theorem 2.2. Let L(λ) and Q(λ) be defined in (2.9) and (2.12), respectively.
Then
σ(L(λ)) = {0, . . . , 0} m
∪ σ(Q(λ)). Here, σ(·) denotes the spectrum of the associated matrix pencil.
Proof. We first note from (2.8) that MN = M + M1, FN = F + F1, and G = GN− G1. The λ-matrix L(λ) in (2.9) can then be rewritten as
(2.14) L(λ) = ⎡ ⎣ K − λ(M + M0 1) K − λM0 1 E − λ(F + FE − λF1 1) E− λ(F+ F1) −E+ λF1 −λG ⎤ ⎦ . Letting J := ⎡ ⎣I0n −I0n 00 0 0 Im ⎤ ⎦ , P := ⎡ ⎣−I0n IInn 00 0 0 Im ⎤ ⎦ , we can further transformL(λ) in (2.14) to a symmetric λ-matrix: (2.15) J PL(λ)P = ⎡ ⎣ −K + λMK − λM11 K − λM−λM 1 E − λF−λF 1 ET− λF 1 −λF −λG ⎤ ⎦ =−K + λM1 S − λT1 (S − λT1) −λM ,
whereS, T1, and M are the matrices defined in (2.7) and (2.8).
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2409 Next, we will show that (2.15) can be reduced to a block diagonal form using Gaussian eliminations. In fact, by considering the λ-matrix
C(λ) := In 0 1 λM−1(S − λT1) In+m ,
and settingE(λ) := (C(λ))J P and F(λ) := PC(λ), we can compute that
E(λ)L(λ)F(λ) = (C(λ))(J PL(λ)P)C(λ) (2.16) = −K + λM1+λ1(S − λT1)M−1(S − λT1) 0 0 −λM = 1 λQ(λ) 0 0 −λM ,
where the last equality follows from (2.11) and (2.12).
Thanks to det(E(λ)) = 1 = det(F(λ)) and the nonsingularity of M, we have det(L(λ)) = det(E(λ)L(λ)F(λ)) = det(λ1Q(λ)) det(−λM) = 0
⇔ 1
λndet (Q(λ)) λ
n+mdet(−M) = 0 ⇔ λmdet (Q(λ)) = 0.
This implies that Q(λ) preserves 2n eigenvalues of L(λ) and throws away m nonphys-ical zero eigenvalues.
Remark 2.3. A similar result as in Theorem 2.2 for the 2D transmission
eigen-value problems has been discussed in [25]. Due to the singularity of S, we know that the matrix K in (2.6) is singular. However, for the 2D transmission eigenvalue problems, K obtained from the discretization of the Laplacian operator is nonsingu-lar. Therefore, the proof technique in [25] based on the nonsingularity of K cannot be directly applied to Theorem 2.2. In Theorem 2.2, we provide a more general proof.
The result in (2.16) indicates that the QEP (2.12) obtained by applying the mixed FEM for the quad-curl problem (2.3) can be directly deduced from the GEP (2.9) discretized by a curl-conforming FEM of (2.1). It is worth considering the QEP (2.12) compared with the GEP (2.9) as the former eliminates m nonphysical zero eigenvalues and maintains the other ones of the latter equation. However, the QEP (2.12) still contains a huge number of zero eigenvalues due to the large null space of S in (2.12d) associated with the curl operator [2]. Because the smallest positive eigenvalues are interesting, these zero eigenvalues lead to numerical difficulties in computing the desired eigenpairs. Additionally, to find the desired positive eigenvalues surrounded by complex eigenvalues is another challenge.
To remedy these difficulties, in what follows, we will introduce a secant-type iterative method [25] in section 3 so that we can sequentially compute the wanted positive eigenvalues without computing any complex ones. In addition, the LOBPCG method [23] will be introduced in section 4 to prevent the disturbance from the huge presence of zero eigenvalues.
3. A secant-type method for computing positive transmission eigen-values. In this section, we focus on the numerical method for finding a few smallest positive transmission eigenvalues of (2.1), which are of great interest for estimating the index of refraction in inverse scattering theory.
To avoid the influence of complex and zero eigenvalues, we first consider a par-ticular symmetric definite GEP with a parameter μ (μ-SDGEP) for the QEP (2.12)
(3.1) A(μ)p(μ) = β(μ)A0p(μ), A(μ) := −A1− μA2,
where A(μ) is symmetric and A0 is symmetric positive semidefinite.
Theorem 3.1. Consider the μ-SDGEP (3.1). Let βi(μ) be the eigenvalue curve
of the matrix pair (A(μ), A0), i = 1, . . . , n. Then
(i) βi(μ) is either real or infinity for any μ ∈ R, i = 1, . . . , n;
(ii) each real eigenvalue curve βi(μ) is strictly decreasing in μ;
(iii) (λ, p) is a real eigenpair of the QEP (2.12) with pA0p = 1 if and only if
(β(λ), p) is a real eigenpair of the μ-SDGEP (3.1) and
β(λ) = 1 λ.
Proof. The proof is similar to Lemma 1 in [25] but with the positive definiteness
of A0 replaced by A0 0.
Remark 3.2. There is another μ-SDGEP of the form in [33]
A(μ)p(μ) = α(μ)Kp(μ), A(μ) := μ 2M1+ (S − μT1)M−1(S − μT1) 0, to be considered for solving the QEP (2.12). From this viewpoint, μ is an eigen-value of (2.12) if and only if it is a fixed point of the eigeneigen-value curve α(μ), i.e.,
α(μ) = μ. Although the eigenvalue curves α(μ) are still real, so that solving the
cor-responding fixed-point problem can avoid capturing complex eigenvalues, it cannot guarantee, in this case, that α(μ) is monotonically increasing because the differentia-tion of p(μ)A(μ)p(μ) with respect to μ is, in general, indefinite. This indicates that eigenvalue curves α(μ) could cross each other and the fixed points may not appear in order. Such uncertainty makes the associated fixed-point problem much more compli-cated, and the traditional secant iteration or Newton’s method may lose some desired real eigenvalues.
Based on Theorem 3.1, we see that any real eigenvalue λ of the QEP (2.12) is a fixed point of the eigenvalue curve 1/β(μ), which means β(λ) = 1/λ. In addition, the monotonicity of β(μ) motivates us to exploit the SecTypIt in [25] for sequentially computing desired positive transmission eigenvalues.
We simply explain the idea of the SecTypIt and summarize this update process in Algorithm 1. For details on the SecTypIt algorithm and its implementation, we refer to [25].
Suppose that 0 < μl< μrare two approximate values for a positive eigenvalue λ
of (2.12). Let βl:= β(μl) and βr:= β(μr) be the corresponding points on the strictly decreasing eigenvalue curve β(μ) passing through the point (λ, 1/λ). Here, μlβl< 1
is required to ensure that (μl, βl) can always converge to (λ, 1/λ).
• Update of (μ+l , βl+). At each iteration, SecTypIt first updates (μl, βl)
accord-ing to the location of (μr, βr). If μrβr< 1, the new (μ+l , β +
l ) is set to be (μr, βr) (see
Figure 1(a)); otherwise, for μrβr > 1, μ+l is updated by a fixed-point iteration from
(μl, βl) to the hyperbola curve, that is μ+l = 1/βl, while βl+ is left to be determined
by solving (3.1) with μ = μ+l (see Figure 1(c)). • Update of (μ+
r, βr+). The correction of μ+r depends on the secant line through
the points (μl, βl) to (μr, βr). When this secant line intersects with the hyperbola
curve, μ+r is shifted to the μ-coordinate of the intersection point closer to the vertical
axis (see also Figure 1(a)). For the case in which the secant line and the hyperbola do not intersect each other, we solve the intersection point μ× > μr from the point
(μr, βr) tangent to the hyperbola curve and modify μ+r by μ×. Finally, we compute
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2411 Algorithm1. [μ+ l , μ + r, β + l , flag] = SecTypIt(μl, μr, βl, βr) [25].
Input: Two approximate solutions (μl, βl) and (μr, βr) to the fixed point (λ, 1/λ).
Output: The updated values (μ+l , βl+) and μ+r. 1: if μrβr< 1. then 2: Set flag = 0. 3: μ+l = μr and βl+= βr. 4: else 5: Set flag = 1. 6: μ+l = 1/βland β+l = []. 7: end if 8: Compute α2= βr−βl μr−μl and α1= βl− α2μl. 9: Set Δ = α21+ 4α2. 10: if Δ > 0 then 11: Set μ+r =−α1+sign(α1) √Δ 2α2 . 12: else 13: Set μ+r = 1+√1−βrμr βr . 14: end if
the associated βr+ by solving (3.1) with μ = μ+r so that we end up with a one-step
iteration for capturing the fixed point (λ, 1/λ) (see Figure 1(b)).
Note that, no matter what the case may be, we have to solve a corresponding
μ-SDGEP (3.1) with an updated μ parameter, and the cost as well as the technique
for solving (3.1) dominate the efficiency and accuracy of this iterative method for capturing the desired positive transmission eigenvalues. In fact, for any fixed μ > 0, one can see that the desired positive eigenvalues of (3.1) suffer from the disturbance of a cluster of infinite eigenvalues. This is because (3.1) consists of an indefinite matrix A(μ) and a positive semidefinite matrix A0, and the nullity of A0 is quite
large. To study this issue, in the following two sections, we will introduce an efficient and robust eigensolver, called LOBPCG [22, 23], that can exclude the disturbance of infinite eigenvalues when solving the μ-SDGEP (3.1).
Remark 3.3. The μ-SDGEP (3.1) has been studied in [25]. As stated in
Re-mark 2.3, the matrices A0 and K in [25] are nonsingular. So, one can solve the
μ-SDGEP (3.1) by the invert Lanczos method and the associated linear system by
the direct method with the Sherman-Morrison-Woodbury formula. However, these techniques fail when A0 and K are singular matrices. That is why we need to
intro-duce the LOBPCG method for solving (3.1).
4. LOBPCG method. Solving (3.1) is a very crucial point for using the Sec-TypIt to update the approximate eigenvalue μ. An appropriate choice of the solver will help to improve efficiency and effectiveness for capturing the desired eigen-values.
At first glance, the shift-and-invert Lanczos method (SILM) seems a feasible ap-proach as we are interested in finding a few desired eigenvalues of (3.1). However, we can immediately note that applying the SILM to solve (3.1) has some drawbacks. (i) The nullity of A0 in (3.1) is huge, and the large dimension of the null space leads to several numerical difficulties [16, 17]. (ii) When the desired eigenpairs of (3.1) are convergent, it is natural to use the associated eigenvectors as the initial vectors for the next μ-SDGEP to accelerate the convergence. However, only one vector in the convergent eigensubspace can be used as an initial vector when the SILM is applied to solve (3.1).
(a) classical secant update
(b) pseudosecant update
(c) mixed secant update
Fig. 1. SecTypIt.
To settle these drawbacks, we apply the LOBPCG method to solve (3.1). LOBPCG was proposed by Knyazev [22] to compute the smallest eigenvalues of matrix pencil
A − λB, where A is Hermitian and B is Hermitian positive definite. For the case in
which B is an indefinite matrix, two variants of LOBPCG were recently studied in [23]. In what follows, we will show that the LOBPCG method can dramatically ex-clude the influence of the infinite eigenvalues and efficiently find some largest positive eigenvalue of (3.1). To begin with, we briefly recall some fundamental properties.
Definition 4.1. Let A and B be n × n Hermitian matrices.
(i) In(A) = (s+, s−, s0) is defined to be the inertia of A, i.e., s+, s−, and s0 are the numbers of positive, negative, and zero eigenvalues of A, respectively.
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2413 (ii) The matrix pencil A − λB is called a positive definite matrix pencil if there
is a shift λ0∈ R such that A − λ0B is positive definite.
Theorem 4.2 (see [23, 24, 26]). Let A − λB be a positive definite matrix pencil.
Then, there is an invertible matrix W such that
W∗AW − λW∗BW = diag (Λ+, −Λ−, Is0)− λ diag Is+, −Is−, 0s0 , (4.1)
where Λ+= diag(λ+1, . . . , λ+s+) and Λ−= diag(λ−1, . . . , λ−s−) with
λ−s− ≤ · · · ≤ λ
−
1 < λ+1 ≤ · · · ≤ λ+s+.
(4.2)
From the factorization of (4.1), it is clear that A − λ0B is positive definite if and
only if λ−1 < λ0 < λ+1. The next theorem is an extension of the classical Cauchy
interlacing theorem for definite pencils.
Theorem 4.3 (see [23, Theorem 2.3]). Let A − λB be a positive definite matrix
pencil and U ∈ Cn×p have full column rank. Then, the eigenvalues of the matrix pencil (U∗AU, U∗BU ) are real and can be ordered as
θ−p−≤ · · · ≤ θ−1 < θ
+
1 ≤ · · · ≤ θ+p+ with In(U∗BU ) = (p+, p−, p0). Moreover,
λ+i ≤ θi+≤ λ+i+n−p for 1≤ i ≤ p+,
λ−j ≥ θ−j ≥ λ−j+n−p for 1≤ j ≤ p−.
Based on the results in Theorem 4.3, Kressner, Pandur, and Shao [23] obtained a corresponding Ky Fan-type theorem (trace minimization principle) and used it to develop two indefinite variants of the LOBPCG method. Algorithm 2 is the indefi-nite LOBPCG method with one preconditioner for computing the smallest positive eigenvalues of a positive definite pencil A − λB.
Note that some largest, although noninfinite, positive eigenvalues of (3.1) are of interest for the modification of μ. So, to solve it with LOBPCG, which is of benefit for computing some smallest eigenvalues, we need to rewrite (3.1) as follows:
(4.3) A0p(μ) = λ(μ)A(μ)p(μ), λ(μ) := 1
β(μ).
Suppose that we are interested in finding smallest positive eigenvalues of (4.3), which has a large number of zero eigenvalues due to the singularity of A0. To satisfy
the requirement for using the LOBPCG method, we assume, for a given μi > 0,
that there exists a sufficiently small λi,0 > 0 such that A0− λi,0A(μi) is a positive
definite matrix pencil. In general, this assumption is reasonable because the norm of
A0 dominates those of A1 and A2. Let λ−i,s− ≤ · · · ≤ λ−i,1 < λ+i,1 ≤ · · · ≤ λ+i,s+ be
the eigenvalues of A0− λ(μi)A(μi). As shown in Lemma 2.1 and its proof, we know
that A2 0, A1=−K − SM−1T1− T1M−1S, and A0=SM−1S 0. Consider
the matrix U for which the columns form a basis of Null(S) with the orthogonality condition UA2U = I. Because Null(S)⊆ Null(K) (by the definition of S in (2.7))
and Null(S)⊆ Null(A0), we obtain
UA0U = 0 and UA(μi)U = U(−A1− μiA2)U = −μiUA2U = −μiI.
By Theorem 4.2, it holds that λ−i,j0 = 0 for some j0.
Algorithm2. [λ1, . . . , λ, x1, . . . , x] = LOBPCG(A, B, , X0, ε) [23].
Input: Coefficient matrices A and B, Hermitian positive definite preconditioned T , the number of desired smallest positive eigenvalues , an initial matrix X0∈ Rn×,
and stopping tolerance ε.
Output: The desired eigenpairs (λi, xi) for i = 1, . . . , with λ1≤ · · · ≤ λ. 1: B-orthonormalize X0 such that X0∗BX0= diag(±1).
2: Compute the eigendecomposition (X0∗AX0)V0= (X0∗BX0)V0Λ0. 3: Update X0= X0V0 and set k = 0, P0= [].
4: repeat
5: Compute Rk = AXk− BXkΛk.
6: if ∃j such that Rk(:, j)/((A + |Λk(j, j)|B)Xk(:, j)) ≥ ε then
7: Compute Wk = T Rk.
8: Set Uk = [Xk, Wk, Pk].
9: B-orthonormalize Uk such that Uk∗BUk= diag(±1).
10: Compute the desired eigenpairs (Λk+1, Vk+1) of the matrix pencil
(Uk∗AUk, Uk∗BUk). Here Λk+1∈ R×, Vk+1∈ Rn×.
11: Compute Pk+1 = Uk,2Vk+1,2 and Xk+1 = Uk,1Vk+1,1+ Pk+1, where Vk+1 =
[Vk+1,1∗ , Vk+1,2∗ ]∗and Uk= [Uk,1, Uk,2].
12: Set k = k + 1. 13: end if
14: until all desired eigenpairs are convergent
15: Set Λk= diag(λ1, . . . , λ) and Xk = [x1, . . . , x].
From the assumption that a sufficiently small λi,0 > 0 can always be found,
we observe, by Theorem 4.2 again, that the zero and positive eigenvalues of (4.3) are separated by λi,0, i.e., λ−i,1 = λ−i,j0 = 0 < λi,0 < λ
+
i,1. This also shows that λ+i,1, . . . , λ
+
i, are the desired smallest positive eigenvalues. The above discussion
leads to the following theorem.
Theorem 4.4. Suppose, for a given μi> 0, A0−λ(μi)A(μi) is a positive definite
matrix pencil with eigenvalues ordered as in (4.2). Then, there is a sufficiently small λi,0> 0 such that A0− λi,0A(μi) 0 and
λ−i,s−≤ · · · ≤ λ
−
i,1= 0 < λi,0< λ+i,1≤ · · · ≤ λ +
i,s+.
Together with the results in Theorem 4.3, we conclude that the Ritz values
θ+1, . . . , θ++ for each iteration of the LOBPCG method (Algorithm 2) satisfy
0 < λi,0< λ+i,j≤ θ +
j ≤ λ
+
i,j+n−, 1≤ j ≤ +.
This indicates that the zero eigenvalues will not degrade the computational efficiency. 5. Practical implementation. Suppose that we want to find smallest posi-tive eigenvalues λ1, . . . , λ of the QEP (2.12). Applying the SecTypIt approach
(Al-gorithm 1) to compute λd, we need to solve a sequence of μ-SDGEPs
(5.1) A0p = λA(μ(d)i )p, A(μi(d)) :=−A1− μ(d)i A2
for i = 0, 1, 2, . . . . The sequence of μ-SDGEPs (5.1) is then solved by LOBPCG. In this section, we will propose heuristic strategies for (i) the choice of the preconditioner, (ii) the setting of the initial vectors, and (iii) the criterion of the stopping tolerance to accelerate the convergence of LOBPCG and SecTypIt.
Table 2 collects the notation employed in the next two sections. Note that if the SecTypIt converges to λd at the idth step, we have λ(d)id,d= μ(d)d = λd.
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2415
Table 2
Notation for the μ-SDGEP (5.1).
λd the dth desired positive eigenvalue of the QEP (2.12), d = 1, . . . . μ(d)i the approximate value for λdat the ith SecTypIt, i = 0, 1 . . . .
(λ(d)i,j, p(d)i,j) the eigenpairs of (5.1) with a given μ(d)i , j = 1, . . . , n.
5.1. Solving linear systems. As presented in line 7 of Algorithm 2, we have to solve a linear system Wk = T Rkwith an appropriate preconditioner T , which is an
essential factor dominating the convergence of the LOBPCG method. As mentioned in [1], we take T as T = A0− τA(μ(d)i ) −1 = A0− τ(−A1− μ(d)i A2) −1 , (5.2)
where τ is a shift value. That is, the modified directions from computing Wk are
parallel to those obtained from one step of the inverse power iteration on the residual
Rk. In SecTypIt, we need to compute d smallest positive eigenvalues λ(d)i,1 ≤ · · · ≤ λ(d)i,d
for (5.1) and use λ(d)i,d to produce the new μ (d)
i+1. This indicates that we can focus on
the improvement of the convergence for computing λ(d)i,d. Therefore, the shift value τ in (5.2) can be chosen as closer to the desired eigenvalue λ(d)i,d.
Here, we take τ = 0.85θk,d, where θk,1 ≤ · · · ≤ θk,d are the smallest positive
Ritz values in the kth iteration of LOBPCG for solving (5.1). In computing the first eigenvalue λ1 of (2.12), the initial vectors of (5.1) with the initial guess μ(1)0 are
randomly constructed. We can see that, in the first few iterations of LOBPCG, the Ritz values are far away from λ(1)1,1. In practice, τ is kept fixed as a given target value
for the first few iterations of LOBPCG.
From (5.2), computing Wk = T Rk is equivalent to solving linear systems
A0+ τ (A1+ μ(d)i A2)
y = (Rk)j,
(5.3)
where (Rk)j is the jth column of the residual matrix Rk. However, due to the
com-plexity of A2, A1, and A0in (2.12b)–(2.12d), it is still challenging to solve (5.3). The
difficulties are (i) the direct methods can hardly be directly applied to solve (5.3) because the matrices A0, A1, and A2 are fully dense; (ii) the iterative methods for
solving (5.3) are not efficient because a suitable preconditioner is, in general, not avail-able. To remedy these drawbacks, we enlarge the linear system (5.3) to augmented linear systems (see (5.6) and (5.7), respectively, below) according to the cases of the refractive indices so that the augmented systems can then be solved by the direct methods.
We first consider N (x) = n0I3 for some positive constant n0 > 1. In this case,
the QEP (2.12) can be further simplified as follows: [λ2(n 0M1) A2 + λ(−(n 0+ 1)K) A1 +SM−11 S A0 ]p = 0 and we have A0+ τ (A1+ μ(d)i A2) =SM1−1S− τ(n0+ 1)K + τ μ(d)i n0M1.
Let
(5.4) u := M−11 Sy ⇒ M1u − Sy = 0.
Then, (5.3) implies that
Su + −τ(n0+ 1)K + τ μ(d)i n0M1
y = (Rk)j.
(5.5)
Combining (5.4) and (5.5), we obtain the augmented linear system as −τ(n0+ 1)K + τ μ(d)i n0M1 S −S M1 y u = (Rk)j 0 . (5.6)
Proceeding similarly, for general nonconstant index of refraction, we enlarge (5.3) into the augmented linear system
⎡ ⎣ M 0 −S 0 M −T1 S − τT1 τ (μ(d)i T1− S) τ(μ (d) i M1− K) ⎤ ⎦ ⎡ ⎣uv y ⎤ ⎦ = ⎡ ⎣ 00 (Rk)j ⎤ ⎦ . (5.7)
5.2. Initializations of the LOBPCG. When the LOBPCG is applied to solve (5.1), we compute the first d = min{d + 2, } smallest positive eigenvalues λ(d)i,1 ≤ · · · ≤ λ(d)i,d and the associated eigenvectors p(d)i,1, . . . , p(d)i,
d of (5.1). Because the new
μ(d)i+1 is produced by SecTypIt with μ(d)i , we naturally use p(d)i,1, . . . , p(d)i,d as the initial vectors of LOBPCG for solving the μ-SDGEP (5.1) with the new μ(d)i+1. Moreover,
when the sequence{λ(d)i,d} converges to λd at i = id, the eigenvalue λ(d)i
d,d+1 and the eigenvector matrix [p(d)i
d,1· · · p
(d)
id,d] can also be chosen as good initial approximations for μlin SecTypIt and for X0in LOBPCG, respectively, for finding the next λd+1.
5.3. Stopping tolerance for theμ-SDGEP (5.1). The stopping criteria can be divided into outer (SecTypIt) and inner (LOBPCG) criteria.
For each SecTypIt procedure, we inspect if the sequence{μ(d)i } converges to the desired eigenvalue λd by checking
(5.8) |μ (d) i − μ(d)i−1| μ(d)i = |λ (d) i,d − λ (d) i−1,d| λ(d)i,d < tol := 10 −8.
On the other hand, the LOBPCG method, as an inner iteration, aims to compute the dth desired eigenpair (λ(d)i,d, p(d)i,d) of (5.1) for the update of (μ(d)i , βi(d)) with β(d)i =
1/λ(d)i,d. So, at each step of LOBPCG, we only need to measure the magnitude of
the relative residual of (λ(d)i,d, p (d)
i,d) in Line 6 of Algorithm 2. From the definition of
the matrices A2, A1, and A0 in (2.12b)–(2.12d), the quantity A0F +|λ|(A1F + |μ|A2F) is roughly approximated by [K E]2F +|λ|[K E]F, where · F is
the Frobenius norm. Therefore, to verify the convergence of the dth Ritz eigenpair (λ(d)i,d, p(d)i,d), we use the normalized residual norms (NRes) defined by
NRes(d)i,d = A0p (d) i,d+ λ (d) i,d(A1+ μ (d) i A2)p (d) i,dF
([K E]F+ λ(d)i,d)[K E]Fp(d)i,dF.
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2417 We then introduce an adaptive stopping criterion for the LOBPCG method ac-cording to the step number i of the SecTypIt approach. Given a suitable initial guess
μ(d)0 , we will choose a corresponding tolerance ε(d)0 for computing β0(d)= 1/λ(d)0,1 from
(5.1) with i = 0. For the subsequent iterations, we tighten the tolerance of LOBPCG according to the outer iteration number i given by
ε(d)i = max{10−13, ε(d)i−1/10} for i = 1, 2, . . . .
In other words, when the sequence of approximate eigenvalues{μ(d)i } is getting closer to the exact solution λd, the stopping criterion becomes increasingly tight to ensure
that the SecTypIt is applied on the exact eigencurve passing through (λd, 1/λd). So,
how to determine the initial tolerance ε(d)0 for each d ≥ 1?
At the very beginning of the SecTypIt, μ10can be selected by any positive number
sufficiently small that it may be far away from the exact eigenvalue λ1, and to save the
computational cost of LOBPCG, we only need a rough approximation of β0(1)= 1/λ(1)0,1
from (5.1) with d = 1 and i = 0. For d ≥ 2, to correct the accuracy of μ(d)0 ,
the corresponding tolerance is dependent on the NRes(d−1)i
d−1,d, where id−1denotes the iteration number i of SecTypIt satisfying (5.8). Based on the above description, we set ε(d)0 = 10−7 if d = 1, min(10−7, NRes(d−1)i d−1,d) if d ≥ 2.
Now, we have Algorithm 3, which summarizes the practical procedure for solving a few positive eigenvalues of the QEP (2.12) by SecTypIt [25] combined with the indefinite LOBPCG with one preconditioner [23].
6. Numerical results. In this section, we demonstrate some numerical results for computing the 6 smallest positive eigenvalues on two domains [33]: (i) the unit ball D1centered at the origin and (ii) the unit cube D2defined as [0, 1] × [0, 1] × [0, 1]. The tetrahedra mesh is used to construct the meshes for D1 and D2.
All computations in this section are carried out in MATLAB 2014b. The systems in (5.6) and (5.7) are solved by the direct method. For the hardware configuration, we use an HP workstation that is equipped with two Intel Quad-Core Xeon E5-2643 3.33 GHz CPUs, 96 GB of main memory, and the RedHat Linux operating system.
The benchmark problems contain three different types, N1, N2, and N3, for the
index of refraction N (x). N1 and (N2, N3) correspond to isotropic and anisotropic
mediums, respectively, with constant index of refraction given by
N1(n0) = n0I3, N2(n0) = ⎡ ⎣n10 n10 00 0 0 n0− 2 ⎤ ⎦ , N3(n0) = ⎡ ⎣n10 n01− 1 0.50.8 0.5 0.8 n0− 1.5 ⎤ ⎦ . The full and zoom-in spectra of the QEP in (2.12) with different N (x) are shown in Figure 2.
6.1. Numerical correctness validation. We validate the correctness of the proposed algorithm by solving the benchmark problem for domain D1with mesh size h ≈ 0.05 and N (x) = N1(n0) = 16I3. The matrix sizes of the associated matrices K and G are 216468 and 12705, respectively. The number of nonzeros of each
ma-trix can be found in Table 3. The values √λ of the 6 smallest positive eigenvalues
Algorithm3. The SecTypIt with LOBPCG for solving the QEP (2.12).
Input: Matrices (A2, A1, A0) in (2.12), the number of desired smallest positive
eigen-values , an initial matrix P0 ∈ Rn×3, tolerance tol, initial values μ0 > 0 and τ > 0.
Output: The desired eigenpairs (λd, pd) for d = 1, . . . , with 0 < λ1≤ · · · ≤ λ. 1: for d = 1, . . . , do
2: Set d = min(d + 2, ) and ε = min(10−7, 50 × NResd), where NRes1= 1.
3: % Fixed point iteration to generate μ1 and μ2 4: for i = 0, 1 do
5: Set B = −A1− μiA2. 6: Call [λ+1, . . . , λ+
d, p
+
1, . . . , p+d] = LOBPCG(A0, B, d, Pi, ε)with solving an augmented linear system (5.6) or (5.7) to compute Wk in line 7 of
Algo-rithm 2.
7: Set μi+1= λ+d, Pi+1 = [p+1, . . . , p+
d] and ε = max{10 −13, ε/10}. 8: if i = 0 then 9: Set μl= λ+d and βl= 1/λ+d. 10: else 11: Set μr= λ+d, βr= 1/λ+d, and i = 2. 12: end if 13: end for
14: while (|μi−1− μi−2|/μi−1 ≥ tol) do
15: Call [μl, μr, βl, flag] = SecTypIt(μl, μr, βl, βr).
16: if flag = 1 then
17: Set μi= μl and B = −A1− μiA2. 18: Call [λ+1, . . . , λ+
d, p
+
1, . . . , p+d] = LOBPCG(A0, B, d, Pi, ε)with solving an augmented linear system (5.6) or (5.7) to compute Wk in line 7 of
Algo-rithm 2.
19: Set βl= 1/λ+d, Pi+1 = [p+1, . . . , p+ d].
20: end if
21: Set μi= μr and B = −A1− μiA2. 22: Call [λ+1, . . . , λ+
d, p
+
1, . . . , p+d] = LOBPCG(A0, B, d, Pi, ε)with solving an augmented linear system (5.6) or (5.7) to compute Wk in line 7 of
Algo-rithm 2.
23: Set βr= 1/λ+d, Pi+1 = [p+1, . . . , p+
d], and i = i + 1.
24: Set ε = max{10−13, ε/10}. 25: end while
26: Compute the normalized residual norm NResd+1 for (λ+d+1, p+d+1).
27: Set μ0= λ+d+1 and P0= [p +
1, . . . , p+d, p] with a given random vector p.
28: end for
produced by Algorithm 3 are 1.1669, 1.1670, 1.1670, 1.4623, 1.4623 and 1.4624. Monk and Sun in [28] show the 6 smallest positive eigenvalues by locating the zeros of the determinants in TM and TE modes as 1.1654 with multiplicity 3 and 1.4608 with multiplicity 3, respectively. This shows that our results coincide rather well with these exact transmission eigenvalues.
6.2. Convergence of LOBPCG. We apply indefinite LOBPCG to compute some smallest positive eigenvalues of (5.1) and use one step of the inverse power method to accelerate the convergence. The convergence of LOBPCG will affect the efficiency of Algorithm 3. Now, we demonstrate the convergence from the views of the
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2419
(a) N2(8) for D1 (b) N3(8) for D2
(c) N3(4) for D1 (d) N2(4) for D2
Fig. 2. Full and zoom-in spectra of the QEP in (2.12) with various indices of refraction. The
mesh sizes are h ≈ 0.2 and h = 1/8 for D1 and D2, respectively. The matrix sizes of K are 2,777
and 1,909, respectively, and the matrix sizes of G are 753 and 954, respectively.
Table 3
The matrix dimension m, n (K ∈ Rn×n, E ∈ Rn×m) and the numbers of the nonzero elements of the matrices for the benchmark problems.
N1(n0)
(m, n) Number of the nonzero elements
K E M1(MN) F1(FN) G1(GN)
D1 (12705, 216468) 3476268 75405 3476268 75405 63525
N2(n0)
(m, n) Number of the nonzero elements
K E M1(MN) F1(FN) (G1, GN) D1 (9312, 136833) 2189097 55410 2189097 55410 (46560, 46560)
D2 (13434, 130989) 1442037 37896 2074161 74394 (66738, 72166)
N3(n0)
(m, n) Number of the nonzero elements
K E M1(MN) F1(FN) (G1, GN) D1 (9312, 136833) 2189097 55410 2189097 55410 (46560, 46560)
D2 (13434, 130989) 1442037 37896 2074161 74394 (66738, 73482)
Ritz values and NRes. The matrix sizes of K and G with N3(8) in this benchmark
problem are 130989 and 13434, respectively, for the domain D2. The number of
nonzeros of each matrix can be found in Table 3. The 6 smallest positive eigenvalues are computed.
(a) Ritz values for λ1 (b) Normalized residual norms NRes(1)i,j
(c) Ritz values for λ2 (d) Normalized residual norms NRes(2)i,j
Fig. 3. The smallest and second positive Ritz values and the associated NRes produced by
LOBPCG in solving (5.1) with μ(d)0 , . . . , μ(d)i . The matrix size of K with N3(8) is 130989 for
domain D2.
The Ritz values and the associated NRes for computing the first and second smallest positive eigenvalues are shown in Figure 3. In this figure, we demonstrate the iteration number of SecTypIt, i.e., the number i of μ(d)0 , μ(d)1 , . . . , μ(d)i , and stack
the iteration number of LOBPCG for solving (5.1) with μ(d)0 , . . . , μ(d)i on the horizontal
axis. Because the initial data for (5.1) with μ(1)0 are randomly constructed, as shown
in Figure 3(a), it needs 25 iterations of LOBPCG to compute the first approximate eigenvalue λ(1)1,1. After μ(1)0 has been computed, according to the initialization scheme
in subsection 5.2, the good initial vectors obviously reduce the iteration number of LOBPCG as shown on the horizontal axis of Figure 3. The associated NRes of the Ritz pairs in Figures 3(b) and 3(d) are monotonically convergent to the stopping tolerance in a few iterations for other μ(d)i .
Note that when λ1 is computed, we add an extra random initial vector to the initial subspace. This random vector leads to corresponding NRes larger than those of others, as shown in Figure 3(d). On the other hand, the locking technique is also applied to deflate the convergent eigenpair if needed. When the Ritz vector is deflated, a random vector is added to the searching subspace. (See the first smallest Ritz value in Figure 3(c).)
For different indices of refraction N (x), we also obtain the same behavior about the iteration numbers of LOBPCG. The results are presented in Figure 4. From these numerical results, we see that it is efficient to solve the μ-SDGEP (5.1) by using LOBPCG with our adaptive strategies proposed in section 5.
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2421
(a) N2(16) for D1 (b) N3(16) for D2
(c) N2(8) for D1 (d) N3(8) for D2
(e) N3(4) for D1 (f) N2(4) for D2
Fig. 4. Iteration numbers of LOBPCG for computing each λ1, . . . , λ6 with various N (x). The
matrix sizes of K are 136833 and 130989 for D1and D2, respectively, and the matrix sizes of G are
9312 and 13434, respectively. The number of nonzeros of each matrix can be found in Table 3.
6.3. Convergence of the SecTypIt method. In this subsection, we will dis-cuss the convergence of the SecTypIt with stopping tolerance 10−8 as introduced in subsection 5.3. First, the iteration numbers of SecTypIt in computing the 6 smallest positive eigenvalues for D1 with N2(16), N2(8), N3(4) and D2 with N3(16), N3(8),
N2(4) are shown in Figure 5. For each index of refraction, the iteration number is less than or equal to 16 to compute one desired eigenvalue. This shows that re-gardless of whether the desired eigenvalues are obviously far away from the complex eigenvalues (see Figures 2(a)–2(b)) or are surrounded by the complex eigenvalues (see Figures 2(c)–2(d)), Algorithm 3 can be used to compute the desired eigenpairs efficiently and robustly.
Fig. 5. Iteration numbers of the secant method for computing each λ1, . . . , λ6 with various
N (x). The matrix sizes of K are 136833 and 130989 for D1 and D2, respectively, and the matrix
sizes of G are 9312 and 13434, respectively.
7. Conclusions. This paper focuses on computing a few smallest positive eigen-values of the three-dimensional Maxwell’s transmission eigenvalue problem, which plays an important role in inverse scattering theory. Its discretized matrix eigenvalue problems are related to a non-Hermitian GEP in (2.9) and a symmetric QEP in (2.12), which are deduced from two finite element methods in [28] and [33], respectively, us-ing the lowest N´ed´elec edge elements. We first show that these two problems have the same spectrum, except for the nonphysical zero eigenvalues. However, owing to the degenerate double-curl operator, the QEP still has a large number of zeros. To compute the desired smallest positive eigenvalues, we propose a SecTypIt method by rewriting the QEP as a sequence of μ-SDGEPs.
To avoid the effect of the large nullity of the μ-SDGEP inherited from the QEP, we apply LOBPCG with one preconditioner [23] to solve the μ-SDGEP. Due to the complexity of the coefficient matrices of the QEP, solving the preconditioning linear system becomes a challenging problem. To this end, we propose a novel method to en-large the preconditioning linear system so that one can solve it by the direct/iterative method. Furthermore, some important heuristic strategies for the determination of initial data and stopping tolerances for the SecTypIt and LOBPCG are introduced to accelerate the convergence. The numerical results demonstrate that Algorithm 3 is robust, although the desired eigenvalues are surrounded by complex eigenvalues.
Acknowledgment. The authors thank the anonymous referees for their valuable comments and suggestions.
REFERENCES
[1] Z. Bai and R.-C. Li, Minimization principles for the linear response eigenvalue problem II:
Computation, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 392–416.
[2] A. Bossavit, Computational Electromagnetism: Variational Formulations, Complementarity,
Edge Elements, Academic Press, San Diego, CA, 1998.
[3] F. Cakoni, M. C¸ ay¨oren, and D. Colton, Transmission eigenvalues and the nondestructive
testing of dielectrics, Inverse Problems, 24 (2008), 065016.
[4] F. Cakoni, D. Colton, and H. Haddar, On the determination of Dirichlet or transmission
eigenvalues from far field data, C. R. Math. Acad. Sci. Paris, 348 (2010), pp. 379–383.
AN EFFICIENT EIGENSOLVER FOR MAXWELL’S TEPs A2423
[5] F. Cakoni, D. Colton, and P. Monk, On the use of transmission eigenvalues to estimate the
index of refraction from far field data, Inverse Problems, 23 (2007), pp. 507–522.
[6] F. Cakoni, D. Colton, P. Monk, and J. Sun, The inverse electromagnetic scattering problem
for anisotropic media, Inverse Problems, 26 (2010), 074004.
[7] F. Cakoni, D. Colton, and J. Sun, Estimation of Dirichlet and transmission eigenvalues by
near field linear sampling method, in Proceedings of the 10th International Conference on
the Mathematical and Numerical Aspects of Waves, Vancouver, Canada, 2011, pp. 431–434. [8] F. Cakoni, D. Gintides, and H. Haddar, The existence of an infinite discrete set of
trans-mission eigenvalues, SIAM J. Math. Anal., 42 (2010), pp. 237–255.
[9] F. Cakoni and H. Haddar, On the existence of transmission eigenvalues in an inhomogeneous
medium, Appl. Anal., 88 (2009), pp. 475–493.
[10] F. Cakoni and H. Haddar, Transmission eigenvalues in inverse scattering theory, in Inverse Problems and Applications: Inside Out II, G. Uhlmann, ed., Math. Sci. Res. Inst. Publ. 60, Cambridge University Press, Cambridge, 2012, pp. 527–578.
[11] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, 3rd ed., Appl. Math. Sci. 93, Springer, New York, 2013.
[12] D. Colton, P. Monk, and J. Sun, Analytical and computational methods for transmission
eigenvalues, Inverse Problems, 26 (2010), 045011.
[13] D. Colton, L. P¨aiv¨arinta, and J. Sylvester, The interior transmission problem, Inverse Probl. Imaging, 1 (2007), pp. 13–28.
[14] H. Haddar, The interior transmission problem for anisotropic Maxwell’s equations and its
applications to the inverse problem, Math. Methods Appl. Sci., 27 (2004), pp. 2111–2129.
[15] G. C. Hsiao, F. Liu, J. Sun, and L. Xu, A coupled BEM and FEM for the interior
transmis-sion problem in acoustics, J. Comput. Appl. Math., 235 (2011), pp. 5213–5221.
[16] T.-M. Huang, H.-E. Hsieh, W.-W. Lin, and W. Wang, Eigenvalue solvers for three
di-mensional photonic crystals with face-centered cubic lattice, J. Comput. Appl. Math., 272
(2014), pp. 350–361.
[17] Y.-L. Huang, T.-M. Huang, W.-W. Lin, and W.-C. Wang, A null space free Jacobi–Davidson
iteration for Maxwell’s operator, SIAM J. Sci. Comput., 37 (2015), pp. A1–A29.
[18] X. Ji, J. Sun, and T. Turner, Algorithm 922: A mixed finite element method for Helmholtz
transmission eigenvalues, ACM Trans. Math. Software, 38 (2012), 29.
[19] X. Ji, J. Sun, and H. Xie, A multigrid method for Helmholtz transmission eigenvalue problems, J. Sci. Comput., 60 (2014), pp. 276–294.
[20] A. Kirsch, On the existence of transmission eigenvalues, Inverse Probl. Imaging, 3 (2009), pp. 155–172.
[21] A. Kleefeld, A numerical method to compute interior transmission eigenvalues, Inverse Prob-lems, 29 (2013), 104012.
[22] A. V. Knyazev, Toward the optimal preconditioned eigensolver: Locally optimal block
precon-ditoned conjugate gradient method, SIAM J. Sci. Comput., 23 (2001), pp. 517–541.
[23] D. Kressner, M. M. Pandur, and M. Shao, An indefinite variant of LOBPCG for definite
matrix pencils, Numer. Algorithms, 66 (2014), pp. 681–703.
[24] P. Lancaster and L. Rodman, Canonical forms for Hermitian matrix pairs under strict
equivalence and congurence, SIAM Rev., 47 (2005), pp. 407–443.
[25] T. Li, W.-Q. Huang, W.-W. Lin, and J. Liu, On spectral analysis and a novel algorithm for
transmission eigenvalue problems, J. Sci. Comput., 64 (2015), pp. 83–108.
[26] X. Liang, R.-C. Li, and Z. Bai, Trace minimization principles for positive semi-definite
pencils, Linear Algebra Appl., 438 (2013), pp. 3085–3106.
[27] P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford University Press, Oxford, UK, 2003.
[28] P. Monk and J. Sun, Finite element methods for Maxwell’s transmission eigenvalues, SIAM J. Sci. Comput., 34 (2012), pp. B247–B264.
[29] J. C. N´ed´elec, Mixed finite elements inR3, Numer. Math., 35 (1980), pp. 315–341.
[30] L. P¨aiv¨arinta and J. Sylvester, Transmission eigenvalues, SIAM J. Math. Anal., 40 (2008), pp. 738–753.
[31] J. Sun, Estimation of transmission eigenvalues and the index of refraction from Cauchy data, Inverse Problems, 27 (2011), 015009.
[32] J. Sun, Iterative methods for transmission eigenvalues, SIAM J. Numer. Anal., 49 (2011), pp. 1860–1874.
[33] J. Sun and L. Xu, Computation of Maxwell’s transmission eigenvalues and its applications in
inverse medium problems, Inverse Problems, 29 (2013), 104013.
