NUMERICAL SOLUTION OF QUADRATIC EIGENVALUE PROBLEMS WITH STRUCTURE-PRESERVING METHODS∗
TSUNG-MIN HWANG†, WEN-WEI LIN‡, AND VOLKER MEHRMANN§
Abstract. Numerical methods for the solution of large scale structured quadratic eigenvalue problems are discussed. We describe a new extraction procedure for the computation of eigenvectors and invariant subspaces ofskew-Hamiltonian/Hamiltonian pencils using the recently proposed skew- Hamiltonian isotropic implicitly restarted Arnoldi method (SHIRA).
As an application we discuss damped gyroscopic systems. For this problem we first solve the eigenvalue problem for the undamped system using the structure-preserving method and then use the quadratic Jacobi–Davidson method as correction procedure. We also illustrate the properties of the new approach for several other application problems.
Key words. quadratic eigenvalue problems, skew-Hamiltonian/Hamiltonian pencils, invariant subspace, gyroscopic system, quadratic Jacobi–Davidson method, nonequivalence deflation technique
AMS subject classifications. 15A18, 47A15, 47J10 PII. S106482750139220X
1. Introduction. In this paper we study the numerical methods for computing eigenpairs or invariant subspaces of the quadratic eigenvalue problem
(λ2M + λ(G + εD) + K)x = 0, (1)
where M, C, K, D are square n × n real matrices with M = MT, G = −GT, D = DT, and K = KT. Typically, the matrices M and K or −K represent the mass matrix and the stiffness matrix, respectively; G and εD represent gyroscopic forces and the damping, respectively, of the system. We not only concentrate on the case in which ε is a small parameter, but also discuss in detail the undamped case that ε = 0.
Quadratic eigenvalue problems (1) arise in the solution of initial or boundary value problems for second order systems of the form
M ¨x + (G + εD) ˙x + Kx = f (2)
and in numerous other applications. These include finite element discretization in structural analysis [34], in acoustic simulation of poro-elastic materials [22, 32, 35], or in the elastic deformation of anisotropic materials [2, 19, 33]. See also [42] for a recent survey.
The classical approach in solving the quadratic eigenvalue problems is to turn them into linear eigenvalue problems by introducing a new vector y = λx. In the case
∗Received by the editors July 10, 2001; accepted for publication (in revised form) August 19, 2002;
published electronically February 6, 2003. This research was partially supported by the National Center for Theoretical Science, National Tsing Hua University, Hsinchu, Taiwan.
http://www.siam.org/journals/sisc/24-4/39220.html
†Department ofMathematics, National Taiwan Normal University, Taipei 116, Taiwan, Republic ofChina (min@math.ntnu.edu.tw).
‡Department ofMathematics, National Tsing Hua University, Hsinchu 300, Taiwan, Republic of China (wwlin@am.nthu.edu.tw).
§Institut f¨ur Mathematik, MA 4-5, Technische Universit¨at Berlin, Str. des 17. Juni 136, D-10623 Germany (mehrmann@math.tu-berlin.de). This author was also partially supported by Deutsche Forschungsgemeinschaft research grant Me 790/11-3.
1283
of (1) this leads to the linearized generalized eigenvalue problem
λ
M G + εD
0 I
−
0 −K
I 0
y x
= 0.
(3)
If (λ, [yx]) is an eigenpair of (3), then x is an eigenvector of (1) associated with the eigenvalue λ. The approach (3) allows us to determine eigenpairs numerically, since for the generalized eigenvalue problem like (3) the mathematical theory, numerical methods as well as the perturbation theory are well established [1, 11, 39]. How- ever, a difficulty of the linearization approach is that due to the embedding into the problem of double size, the condition number, i.e., the sensitivity of the eigenvalues and eigenvectors with respect to perturbations in the data matrices M, G, D, K, may increase. This is because the set of admissible perturbations for (3) is larger than for (1). If the perturbations would, however, respect the specific structure of the blocks in (3), i.e., the zero, the identity blocks, and the symmetries, then the perturbation results would be the same. Furthermore, there are many different ways to do the linearization with different conditioning, as has been demonstrated in [40]. In view of these remarks it would be ideal to have a numerical method that works directly with the original data of the quadratic eigenvalue problem and that avoids the problem of the increased condition numbers. It is difficult to develop such a method, and there are only few methods that partially fulfill these requirements, such as the quadratic or polynomial Jacobi–Davidson method [14, 36, 37].
In many applications the original quadratic eigenvalue problem has some extra structures that should be reflected in its linearization. A typical structure is the Hamiltonian eigenvalue symmetry, i.e., that the spectrum is symmetric with respect to the real and imaginary axes. This symmetry occurs if ε = 0 in (1), which is the case, for example, in the gyroscopic system [18] and in the elasticity problems [2, 26]; see [42] for a recent survey. A similar eigenvalue symmetry arises in the optimal control problems for which the variational principle leads to a skew-Hamiltonian/Hamiltonian (SHH) pencil that has similar properties as (3) in the undamped case; see [5, 25, 26].
In these applications the generalized eigenvalue problem arises directly and is not from any linearization of a quadratic eigenvalue problem.
Another structure that should be reflected in a proper method is preservation of the sparsity structure inherited from finite element or finite difference discretizations [3, 26, 34].
In this paper we discuss first the solution of large sparse generalized eigenvalue problems that have a Hamiltonian eigenvalue symmetry and then allow perturbations of this structure. The use of structure-preserving linearizations for problem (1) in the undamped case has first been suggested and successfully employed in [10, 26] to design structure-preserving methods for large sparse linearized quadratic eigenvalues problems with a Hamiltonian eigenvalue symmetry. We make use of these techniques, in particular, the skew-Hamiltonian isotropic implicitly-restarted Arnoldi algorithm (SHIRA) proposed in [26]. We will briefly recall this algorithm in section 2.
The key feature of the SHIRA algorithm is that it inherits the convergence prop- erties, the implicit restart, and the reorthogonalization techniques of the standard shift-and-invert Arnoldi algorithm [20, 31, 38], while being, in general, more efficient.
Furthermore, it generates isotropic Krylov subspaces to guarantee that the computed spectrum has the correct eigenvalue symmetry. In [41] the stability properties of the methods that preserve the Hamiltonian structure (although not for this algorithm) are analyzed.
A disadvantage of the SHIRA algorithm is that it does not directly generate the eigenvectors or invariant subspaces of the linearization. Instead, they are determined in [26] using a few steps of inverse or subspace iteration [11, 44], respectively. For these iterations one has to solve in every step a linear system, and hence extra sparse matrix factorizations are needed. This results in a bottleneck in the computation and also limits the possible system size. Even though, in practice, only a few iterations are needed, these extra factorizations add to the method substantial computational overhead.
In this paper we develop a new extraction method for invariant subspaces in the SHIRA algorithm that avoids extra sparse matrix factorizations. This method is a modification of an idea suggested in [45, 46] for the computation of the stable invariant subspace of a Hamiltonian matrix. We compare the performance of this extraction procedure with that of classical inverse or subspace iteration in section 5.
As an application we then discuss numerical methods for problems like (1) that have a damping term with ε = 0. In this problem the Hamiltonian eigenvalue symmetry is destroyed and therefore structure-preserving methods cannot be used. On the other hand, if the damping is small, then we can expect that the damped system has a spectrum and also invariant subspaces close to those of the structured system.
We discuss this topic in section 4 and introduce a method that first computes the eigenvalues, eigenvectors, or invariant subspaces of the structured system and then use them as initial vectors in the quadratic Jacobi–Davidson procedure [3, 36] to compute the eigenvalues, eigenvectors, or invariant subspaces of the damped system.
A similar approach for problems with small damping in structural mechanics was suggested in [27]. In order to treat eigenvalue problems with clusters of eigenvalues, a nonequivalence low-rank transformation technique [9, 12] is used to deflate the computed eigenpairs. We compare this approach with the classical shift-and-invert subspace iteration (SISI) for several application problems in section 5.
2. Structure-preserving algorithms. In this section we briefly recall the SHIRA algorithm of [26], which is designed for the computation of a small number of specified eigenvalues of real large scale generalized eigenvalue problems for SHH pencils of the form
αN − βH = α
F1 G1
H1 F1T
− β
F2 G2
H2 −F2T
, (4)
where G1= −GT1, H1= −H1T, G2= GT2, and H2= H2T.
Matrices N and H in (4) are called skew-Hamiltonian and Hamiltonian, respec- tively. If J = [−I0nI0n] with In being the n×n identity matrix, then skew-Hamiltonian matrices satisfy (N J)T = −N J and Hamiltonian matrices satisfy (HJ)T = HJ. T he special structure of SHH pencils ensures [23, 24, 25] that the eigenvalues occur in quadruples λ, λ, −λ, −λ, if λ is complex, or pairs λ, −λ if λ is real. Moreover, if x is a right eigenvector associated with an eigenvalue λ, then Jx is a left eigenvector associated with −λ.
If the skew-Hamiltonian matrix N in the SHH pencil αN −βH in (4) is invertible and given in the factored form N = Z1Z2 with Z2TJ = ±JZ1, then the pencil is equivalent to a pencil of the form αI − βW, where
W = ±Z1−1HZ2−1
is Hamiltonian. The factorization of N can be easily determined via a Cholesky-like decomposition for skew-symmetric matrices developed in [7]; see also [4].
Suppose that one wishes to compute the eigenvalues of W nearest to some target value λ0. Then typically one would use shift-and-invert of the form (W − λ0I)−1 to get the eigenvalues near λ0 via a Krylov subspace method. But this matrix is not Hamiltonian any longer. In order to preserve the structure, [26] suggests a rational transformation with four shifts (λ0, λ0, −λ0, −λ0) given by
R1(λ0, W) = (W − λ0I)−1(W + λ0I)−1(W − λ0I)−1(W + λ0I)−1. (5)
If the target λ0 is either real or purely imaginary, one may use the simpler trans- formation
R2(λ0, W) = (W − λ0I)−1(W + λ0I)−1. (6)
If the matrix W is real and Hamiltonian, then both matrices R1and R2are real and skew-Hamiltonian, and hence the eigenvalues have algebraic multiplicity of at least two, if they are not purely imaginary. For a real skew-Hamiltonian matrix K and a given vector q1, an Arnoldi iteration applied to K would generate the Krylov space
V(K, q1, k) = span{[q1, Kq1, . . . , Kk−1q1]}.
Using an appropriate orthogonal basis of this space given by the columns of an or- thogonal matrix Qk, one produces a “Ritz”-projection
Kk = QTkKQk.
The “Ritz”-values, i.e., the eigenvalues of Kk are then used as eigenvalue approxima- tions.
In order to obtain a structure-preserving method, one needs a “Ritz”-projection that is again skew-Hamiltonian. For this one would need an isotropic subspace V, i.e., a subspace for which xTJy = 0 for all x, y ∈ V; see [26]. Let Vk∈ R2n×k and let
QTVk =
Rk
T0k 0
(7)
be the symplectic QR-factorization [8] of Vk, where Q ∈ R2n×2n is orthogonal and symplectic and where Rk and Tk ∈ Rk×k are upper and strictly upper triangular, respectively. If Vk is isotropic and of full column rank, then it is shown in [8] that R−1k exists and that Tk= 0. This means, in particular, that if Vnis of full column rank, then there exists an orthogonal symplectic matrix Q ≡ [Qn| JQn] with Qn ∈ R2n×n, such that
QTn QTnJT
K [Qn JQn] =
Hn Kn 0 HnT
. (8)
Here Hn = QTnKQn is an upper Hessenberg matrix and Kn = QTnKJQn is skew- symmetric. This idea is the basis for the structure-preserving Arnoldi algorithm SHIRA introduced in [26], which is given by the recursion
KQk= QkHk+ qk+1hk+1,keTk, (9)
where Hk is the leading k × k principal submatrix of Hn and ek = [0, . . . , 0, 1]T ∈ Rk×1. Since Qk is orthogonal and QTkJQk = 0 for k < n, it is easily seen that [Qk|JQk] ∈ R2n×k is orthogonal. This implies that
QTk QTkJT
K [Qk | JQk] =
Hk Gk
0 HkT
∈ R2k×2k (10)
is skew-Hamiltonian and Gk= QTkKJQk is skew-symmetric.
Let (θi, vi) be an eigenpair of Hk, i.e., Hkvi = θivi, and let xi = Qkvi be a
“Ritz”-vector of the eigenvalue problem Kx = µx corresponding to the “Ritz” value θi. Then we have the following residual bound, see [31], for the “Ritz”-pair (θi, xi):
Kxi− θixi2= KQkvi− θiQkvi2
=
QkHk+ qk+1hk+1,keTk
vi− θiQkvi2
= qk+1hk+1,k(eTkvi) 2= |hk+1,k||eTkvi|.
Remark 1.
(a) The transformations R1(λ0, W) and R2(λ0, W) can be viewed as a structure- preserving shift-and-invert strategy to accelerate the convergence of the de- sired eigenvalues.
(b) In exact arithmetic, the values ti,k= (Jqi)TKqkin the kth step of the SHIRA algorithm are zero, but in practice roundoff errors cause them to be nonzero.
It has been shown in [26] how to design an isotropic orthogonalization scheme to ensure that the spaces span{q1, . . . , qk} are isotropic to working precision.
(c) To avoid the reorthogonalization, after a number of steps the SHIRA algo- rithm also uses implicit restarts as in [38].
(d) A minor disadvantage of the SHIRA method is that due to the rational trans- formations (5) or (6) the eigenvectors associated with pairs of eigenvalues λ and −λ in the real case or λ and −¯λ in the complex case are both mapped to the same invariant subspace associated with |λ|2. This has the consequence that in many applications the desired invariant subspace associated with the stable eigenvalues is not obtained directly from the method. In section 3 we discuss an extraction method that allows us to compute this subspace directly.
3. An extraction method for stable eigenspaces. In this section we discuss a method for computing specific invariant subspaces of a Hamiltonian matrix W directly from an isotropic invariant subspace of the skew-Hamiltonian matrix W2 without using inverse or subspace iteration. In many applications, we are interested in a subspace associated with eigenvalues in the left half plane. The same method can also be applied to the skew-Hamiltonian functions R1, R2introduced in section 2.
The method that we describe is a modification of a technique first suggested in [45, 46]. Suppose that W has no eigenvalues on the imaginary axis and let Qn be as in (8) for K = W2. Then
W2Qn= QnHn. (11)
Let Qk∈ R2n×k (k ≤ n) be a basis of an isotropic invariant subspace of W2such that W2Qk= QkΩk,
(12)
where for the spectrum σ we have σ(Ωk) ⊆ σ(Hn) and σ(Ωk) ∩ {σ(Hn)\σ(Ωk)} = ∅,
i.e., the complete multiplicity of a multiple eigenvalue is included. Since we have assumed that W has no purely imaginary eigenvalues, it follows that there exist two bases Vk− and Vk+, both of dimension k, associated with the stable and unstable isotropic, invariant subspaces of W, respectively, such that
WVk−= Vk−Λ−k, WVk+= Vk+Λ+k, (13)
where
σ(Λ−k) = {λ|Re(λ) < 0, λ2∈ σ(Ωk)}, σ(Λ+k) = {λ|Re(λ) > 0, λ2∈ σ(Ωk)}.
(14)
Since by assumption Ωk has no purely imaginary eigenvalue, it follows that Ωk
has a unique positive square root Xk[16] satisfying Xk2= Ωk such that all eigenvalues have positive real part. We call this root the positive square root. It can, for example, be computed via the MATLAB function sqrtm [13].
The following theorem shows how to extract the stable eigenspace span(Vk−) from the isotropic invariant subspace span(Qk).
Theorem 1. Let W ∈ R2n×2n be Hamiltonian and let Qk∈ R2n×k be as in (12).
Suppose that there exist nonsingular matrices L−k and L+k ∈ Rk×k such that Qk= Vk+L+k + Vk−L−k,
(15)
where Vk+ and Vk− are as in (13). Then for Xk, the positive square root of Ωk, we have
span(WQk− QkXk) = span(Vk−).
(16)
Proof. From (13) and (15) we have
WQk− QkXk = W(Vk+L+k + Vk−L−k) − (Vk+L+k + Vk−L−k)Xk
= Vk+
Λ+kL+k − L+kXk + Vk−
Λ−kL−k − L−kXk .
If we can prove that Λ+kL+k − L+kXk = 0 and that Λ−kL−k − L−kXkis nonsingular, then assertion (16) follows. To do this, we show first that
Xk = (L+k)−1Λ+kL+k. (17)
From (15) and (13) we have
W2Qk = W2
Vk+L+k + Vk−L−k
= Vk+ Λ+k 2
L+k + Vk− Λ−k 2
L−k. (18)
On the other hand, from (12) and (15) we have W2Qk = QkΩk =
Vk+L+k + Vk−L−k Ωk. (19)
Subtracting (19) from (18) we obtain Vk+ Vk−
Λ+k 2
L+k − L+kΩk
Λ−k 2
L−k − L−kΩk
= 0.
Since
Vk+ Vk−
has full rank it follows immediately that Λ+k 2
L+k − L+kΩk = 0 and (17) follows due to the uniqueness of Xk.
We also have Λ−k 2
L−k − L−kΩk = 0, and the uniqueness of the positive square root of Ωk implies that
Xk = L−k −1
ZkL−k, (20)
where Zk denotes the unique positive square root of Λ−k 2
. Hence, Λ−kL−k − L−kXk = Λ−k − Zk
L−k is nonsingular. This completes the proof.
Theorem 1 indicates a way to determine the desired stable invariant subspace.
We apply SHIRA to the skew-Hamiltonian operators R1(λ0, W) and R2(λ0, W) as in (5) and (6), respectively, and determine the associated isotropic invariant subspace Qk.
If the target λ0 /∈ σ(W) is real or purely imaginary, and if the entry hk+1,kof Hk is negligible, then the space
span(Qk) = span{q1, . . . , qk}
generated by SHIRA is a good approximation of an isotropic invariant subspace of R2(λ0, W). Then Qk satisfies (approximately)
R2(λ0, W)Qk= (W2− λ20I)−1Qk= QkHk, (21)
where Hk is an upper Hessenberg matrix. This implies that Ωk = QTkW2Qk = Hk−1+ λ20I.
(22)
Applying Theorem 1, we can extract the stable isotropic eigenspace Vk− of W from Ωk by computing its positive square root.
If the target λ0 /∈ σ(W) is neither real nor purely imaginary and if hk+1,k is negligible, then the space span(Qk) = span{q1, . . . , qk} generated by SHIRA is a good approximation to an isotropic invariant subspace of R1(λ0, W). However, as has been shown in [26], it may happen that this subspace fails to be invariant under W2.
It is therefore necessary to check whether or not span(Qk) is invariant under W2 by computing the residual
W2Qk− QkΩk F,
with Ωk = QTkW2Qk. If this residual is small, then the procedure to compute the invariant subspace is as before.
We summarize the extraction method in the following algorithm.
Algorithm 1 (Extraction method (EM) for the stable invariant subspace or eigenvectors of a Hamiltonian matrix).
Input Hamiltonian matrix W and a target value λ0 with negative real part.
Output Approximate invariant subspace V−of W associated with eigenvalues of negative real part nearest to λ0.
(i) IF λ0= α or λ0= iα for α ∈ R, THEN apply SHIRA to R2(λ0, W), i.e., (a) generate Arnoldi vectors as columns of Qk= [q1, . . . , qk] and an upper
Hessenberg matrix Hk∈ Rk×k such that
(W2− λ20I)−1Qk= QkHk. (b) Compute Ωk= Hk−1+ λ20I.
(ii) ELSE IF λ0= α + iβ, where α, β ∈ R, then apply SHIRA to R1(λ0, W), i.e., (a) generate Arnoldi vectors as columns of Qk= [q1, . . . , qk] and an upper
Hessenberg matrix Hk∈ Rk×k such that R1(λ0, W)Qk= QkHk. (b) Compute Ωk= QTkW2Qk.
(iii) EndCompute the real Schur decomposition of Ωk as Ωk= UkTkUkT.
(iv) Reorder the stable eigenvalues of Tk to the top of Tkusing the reordering method for the Schur form of [1], i.e.,
Tk= Vk
T11 T12
0 T22
VkT, where T11∈ R×has the stable eigenvalues.
(v) Set ˜Q:= QkUkVk
I
0
.
(vi) Compute the unique positive square root√
T11of T11using the method sqrtm of [21].
(vii) Compute the stable invariant subspace V−= W ˜Q− ˜Q√
T11.
The computational effort for the computation of the square root of Ωk is of order k3. If k is small compared to n, then this cost is negligible and it does not add much to the cost of the SHIRA iteration.
Remark 2. In a worst case analysis this procedure may suffer from a loss of accuracy due to the fact that we are computing the square root of the Schur form of a function of W, which at least contains quadratic terms in W. Thus, effects similar to the square reduced method for the Hamiltonian eigenvalue problem may occur, see [43], where the squaring of W leads to a loss of accuracy for small eigenvalues. See also the analysis in [45]. We will demonstrate this in Example 2 and also show that one step of inverse iteration fixes this accuracy loss.
Algorithm 1 can be applied directly to large sparse problems of the form (4).
An important application here is the linear quadratic optimal control problem for descriptor systems, where the pencil has the form [5, 6, 25]
α
E 0
0 ET
− β
A BBT
CTC −AT
. (23)
Such problems arise in the optimal control of semidiscretized parabolic partial differ- ential equations [15, 28, 29, 30].
Assume that the skew-Hamiltonian matrix N in (4) is given in factored form N = Z1Z2, where Z2TJ = ±JZ1.
(24)
Such a factorization (called J-Cholesky factorization) exists for all real skew-Hamiltonian matrices [5] and it can be obtained either trivially or numerically via a Cholesky-like factorization of skew-symmetric matrices [7].
If N is invertible, then using this factorization we can (at least formally) transform the SHH pencil αN − βH to a standard eigenvalue problem αI − βW, where W =
±Z1−1HZ2−1 is Hamiltonian. If in (23) the matrix E is invertible, then the resulting Hamiltonian eigenvalue problem is
α
I 0 0 I
− β
E−1A E−1BBTE−T CTC −ATE−T
. (25)
Numerical examples of this type are presented in section 5.
4. Damped gyroscopic systems. In this section we show how the extraction procedure given by Algorithm 1 can be used to solve quadratic eigenvalue problems for systems of the form (1). To do this we first study the undamped case (ε = 0), i.e.,
λ2Mx + λGx + Kx = 0, (26)
with M = MT and K = KT positive definite and G = −GT. Using the linearization (3) and the factorization
N = Z1Z2=
I 12G 0 M
M 12G
0 I
yields the Hamiltonian eigenvalue problem
(λI − W)x = (λI − Z1−1HZ2−1)x = 0.
(27)
The matrix (W − λI)−1 can be factored as (W − λI)−1
(28)
=
M 12G
0 I
I λI
0 I
0 M−1
−Q(λ)−1 0
I 12G + λM
0 M
,
where
Q(λ) = λ2M + λG + K.
(29)
To this problem we can directly apply the extraction procedure given by Algorithm 1.
If we include the damping term in the gyroscopic system as in (1), then we cannot use the structured SHIRA method directly. But we can still linearize the system (1) and obtain the perturbed SHH system (3) or, in the different form,
λ
M G
0 M
−
−εD −K
M 0
y x
= 0, (30)
which is now a pencil with one skew-Hamiltonian and one perturbed Hamiltonian matrix. We can still use the factorization (27) of the skew-Hamiltonian matrix and obtain the perturbed Hamiltonian eigenvalue problem
(λI − ˆW)x = (λI − Z1−1HZˆ 2−1)x = 0, (31)
with
W =ˆ
−(εD + 12G)M−1 −K +14GM−1G + εDM−1G
M−1 −12M−1G
= W + ε
−DM−1 DM−1G
0 0
.
If the perturbation ε is small, then eigenvalue/eigenvector pairs (or invariant subspaces) of the problem (1) can be regarded as small perturbations of eigenvalue/
eigenvector pairs (or invariant subspaces) of problem (26). In this situation it is natural to use an eigenvalue/eigenvector pair of the undamped problem (26) computed via Algorithm 1 as start for a correction method. This method could be any method that can be used for eigenvalue, eigenvector, or subspace correction, such as subspace iteration, inverse iteration, Newton’s method, or the Jacobi–Davidson method. Here we present results that are obtained with the quadratic Jacobi–Davidson method [3, 36] as correction. This method is typically very efficient and well suited for the given problem. However, if the desired eigenvalues of (1) form a cluster of nearby eigenvalues, then the quadratic Jacobi–Davidson method sometimes has difficulties in detecting and resolving such a cluster. The undesired effect is that, in this case for different starting values (eigenvalue/eigenvector pairs) obtained from the unperturbed problem, it converges to the same eigenvalue/eigenvector pair of (1). It is known that implicit deflation techniques based on Schur forms (see, e.g., section 4.7 and section 8.4 of [3]) combined with the Jacobi–Davidson method perform well for linear eigenvalue problems. However, in the quadratic eigenvalue problem, it is not clear how to incorporate an implicit deflation technique because a quadratic Schur form, in general, does not exist for quadratic matrix polynomials.
For this reason we have also analyzed the use of an explicit nonequivalence low- rank transformation deflation technique that was suggested in [9, 12] for quadratic eigenvalue problems. Let us briefly recall this technique for the polynomial
L(λ) = λ2M + λC + K with C = G + εD.
(32)
We study the two cases of real eigenvalues or complex conjugate pairs separately.
Suppose that we have computed a real eigenvalue λ1 as well as the associated right and left eigenvectors x1 and z1, respectively, with z1TKx1= 1 such that L(λ1)x1= 0 and zT1L(λ1) = 0. Let θ1= (z1TMx1)−1. We then introduce a new deflated quadratic eigenproblem as in [12] via
L(λ)x ≡
λ2M + λ C + K x = 0, (33)
where
M = M − θ 1Mx1z1TM, C = C + θ1
λ1(Mx1z1TK + Kx1z1TM), (34)
K = K − θ1
λ21Kx1zT1K.
Suppose that we have computed an complex eigenvalue λ1 = α1+ iβ1 as well as the associated right and left eigenvectors x1 = x1R+ ix1I and z1 = z1R+ iz1I, respectively, such that Z1TKX1= I2, where X1 = [x1R, x1I] and Z1= [z1R, z1I]. Let
Θ1= (Z1TMX1)−1. We then introduce a new deflated quadratic eigenproblem as in [9] via (33), where
M = M − MX 1Θ1Z1TM,
C = C + MX 1Θ1Λ−T1 Z1TK + KX1Λ−11 ΘT1Z1TM, (35)
K = K − KX 1Λ−11 Θ1Λ−T1 Z1TK in which Λ1= [−βα11βα11].
Note that if the matrices M and K are symmetric, then the matrices M and K in (34) or (35) are symmetric as well. The results in [9, 12] then imply the following proposition.
Proposition 1. (i) Let λ1 be a simple real eigenvalue of L(λ) as in (32) and let x1 and z1 with z1TKx1 = 1 be the right and left eigenvectors, respectively. Then the spectrum of L(λ) in (33) with coefficients as in (34) is given by (σ(L(λ)){λ1}) ∪ {∞} = σ(L(λ)) provided that λ21= θ1.
(ii) Let λ1 be a simple complex eigenvalue of L(λ) as in (32) and let X1 = [x1R, x1I] and Z1 = [z1R, z1I] with Z1TKX1 = I2, where x1 = x1R+ ix1I and z1 = z1R+iz1I are the associated right and left eigenvectors, respectively. Then the spectrum of L(λ) in (33) with coefficients as in (35) is given by
σ(L(λ)){λ1, ¯λ1}
∪{∞, ∞} = σ(L(λ)) provided that Λ1ΛT1 = Θ1.
Furthermore, in both cases (i) and (ii), if λ2= λ1 and (λ2, x2) is an eigenpair of L(λ), then the pair (λ2, x2) is also an eigenpair of L(λ).
Using this deflation procedure we have now presented all the ingredients for our algorithm. We solve the damped gyroscopic system (1) by the quadratic Jacobi–
Davidson method combined with the explicit nonequivalence deflation technique (33), (34) or (33), (35) by using eigenpairs computed via Algorithm 1 as starting values.
We summarize this approach in the following algorithm.
Algorithm 2 (Quadratic Jacobi–Davidson method with deflation).
Input Matrices M, G, D, and K and parameter ε as in (1). Target shift λ0 and number of desired eigenvalue/eigenvector pairs nearest to λ0. Tolerance Tol for stopping criterion.
Output The eigenpairs {(λj, xj)}j=1 of
L(λ)x = (λ2M + λ(G + εD) + K)x = 0, associated with eigenvalues {λj}j=1that are closest to the target λ0. (i) Compute the eigenpairs {(λ(0)j , x(0)j )}j=1of
(λ2M + λG + K)x = 0
by Algorithm 1, where {λ(0)j }j=1 are the closest eigenvalues to the target λ0
and (x(0)j )Hx(0)j = 1.
(ii) For j = 1, . . . ,
Compute the eigenpair (λj, (xj, zj)) of L(λ)by the quadratic Jacobi–Davidson method with target λ(0)j and initial vector x(0)j , where xjand zjare the associated right and left eigenvectors satisfying the relations as in (i)or (ii), respectively, of Proposition 1.
If xj− xi ≤ T ol for some i < j, then
Compute the eigenpair (λj, (xj, zj)) ofL(λ)as in (33), (34)if λ i is real, and as in (33), (35) if λi is complex, by the quadratic Jacobi–Davidson method with target λ(0)j and initial vector x(0)j .
End if End for
Note that a recently proposed locking technique for the solution of quadratic eigenvalue problem as suggested in [22] might also be adapted instead of our non- equivalence low-rank deflation technique in Algorithm 2.
5. Numerical results. In this section we present some numerical tests for the algorithms proposed in this paper. All computations were done in MATLAB 6.0 [21]
or in Fortran 90 on a Compaq DS20E workstation.
We first discuss the new extraction method and compare it with several other iterative methods. These methods are as follows:
EM: the SHIRA algorithm combined with the extraction method given by Algorithm 1;
SISI(q): the SHIRA algorithm followed by q steps of shift-and-invert subspace iteration;
EM SISI: the SHIRA algorithm combined with the extraction method EM adding one step of shift-and-invert subspace iteration;
IPI: one step of the inverse power iteration.
Here, in the SISI the target value λ0 with negative real part is taken as a fixed shift and the iteration starts with the subspace span(Qk) generated by SHIRA, which includes the eigenvectors associated with eigenvalues near λ0and −λ0. Then SISI con- verges to a subspace span(Q−k), associated with stable eigenvalues and the eigenpairs are computed from the “Ritz”-pairs of (Q−k)TWQ−k.
In the EM SISI variant, the target value with negative real part is again taken as a fixed shift and one extra step of SISI is performed starting with span(Vk−) obtained by the new extraction method EM. With the resulting subspace span(Vk), then the eigenpairs are computed from the “Ritz”-pairs of VkTWVk.
In the inverse power iteration (IPI) an approximate eigenvalue computed by SHIRA is taken as shift and the associated Arnoldi vector qi is used as starting vector.
In the following we present results for three problem classes: quadratic eigen- value problems from elasticity theory, linear quadratic optimal control problems, and damped gyroscopic systems.
5.1. Quadratic eigenproblems from elasticity theory. Consider the quadratic eigenvalue problem
λ2Mx + λGx + Kx = 0, (36)
where M, G, and K are defined as follows. As in [26], let
B =
0 0
1 ...
... ...
0 1 0
∈ Rm×m,
M =˜ 16(4Im+ B + BT), ˜G = B − BT, and ˜K = −(2Im− B − BT). Define M = c11Im⊗ ˜M + c12M ⊗ I˜ m,
G = c21Im⊗ ˜G + c22G ⊗ I˜ m, K = c31Im⊗ ˜K + c32K ⊗ I˜ m, (37)
−0.14 −0.12 −0.1 −0.08 −0.06 −0.04 (1a) eigenvalues
target value
−0.14 −0.12 −0.1 −0.08 −0.06 −0.04
10−15 10−10 10−5
(1b) eigenvalues
residuals
residuals by EM residuals by SISI(10) residuals by IPI
Fig. 1. Eigenvalues and the corresponding residual.
where cij are positive constants. Then we have M = MT > 0, G = −GT, and K = KT.
Example 1 (see [26, Example 6.2]). We take m = 90, so that W is a 16200×16200 matrix, and
c11= 1.00, c12= 1.30, c21= 0.10, c22= 1.10, (38)
c31= 1.00, c32= 1.20.
For this example, the twelve eigenvalues computed by the SHIRA algorithm clos- est to the target λ0= −0.1 are depicted in Figure 1(a). The residuals (λ2iM +λiG+
K)xi 2 for the eigenpairs (λi, xi), where xi is computed via EM, SISI(10), and IPI, respectively, are given in Figure 1(b).
We see from Figure 1(b) that the methods EM and IPI yield residuals of the same magnitude, ≈ 10−10, while SISI(10) yields smaller residual for some eigenvalues but the residuals of the eigenvalues further away from the target are much larger.
In order to evaluate the computational complexity for this problem class, we compare the major computational tasks per iteration step. One iteration of the inverse power iteration for (36) requires one forward/backward substitution for evaluating Q(λ)−1z (assuming that an LU-factorization for Q(λ) is available). On the other hand, from (29), we see that the matrix W can be factored as
W =
I −12GM−1
0 M−1
0 −K
I 0
I −12G
0 I
.
The dominant cost of the extraction method is to determine the eigenvector in step (vii) of Algorithm 1. It consists of matrix-vector products for evaluating Kz and Gz, respectively, forward/backward substitution for evaluating M−1z (assuming that an
Table 5.1
Computational costs for EM and IPI; here k = = 12.
EM IPI
0 evaluations Q(λi) evaluations Q(λi)
1 sparse LU factorization of M sparse LU factorizations of Q(λi)
f/b substitutions of M−1z f/b substitutions of Q(λi)−1z
matrix-vector products Kz 0 matrix-vector products Kz 2 matrix-vector products Gz 0 matrix-vector products Gz
(2 + k) SAXPY operations 0 SAXPY operations
LU or Cholesky-factorization for M is available) as well as 2 + k SAXPY operations.
The cost for the evaluation of Ωk, Uk and Vk in Algorithm 1 are negligible if k n.
In Table 5.1 we compare the computational cost for the computation of eigenvec- tors using the EM and one step of the IPI, respectively. This comparison shows that the computational cost for the new extraction method compares favorably with that of one step of inverse power iteration. The major savings in computational time arise from the fact that no further factorizations of Q(λ) are needed.
5.2. Optimal control problems. The second class of problems arises from continuous-time linear quadratic control problems of the form
minx,u
∞
0 (uTRu + xTCTCx)dt subject to
˙x = Ax + Bu, y = Cx, (39)
where A ∈ Rn×n, B ∈ Rn×m, C ∈ Rp×n, R ∈ Rm×mwith p, m n and where R is symmetric positive definite.
As test cases we consider the spatial (central difference) discretization in polar coordinates of a reaction-diffusion equation with Dirichlet boundary conditions on the two-dimensional unit disk; see [17] for details.
Performing a semidiscretization in space, one obtains a continuous-time system as in (39),
A = I + T, (40)
and there exists an orthogonal matrix W , see [17], such that the matrix A in (40) can be transformed to a tridiagonal matrix A ≡ WTAW . Thus, the LU-factorization of A is easily computed in O(n) operations [11]. Let
WT(A − λ0I)W = L1U1, −WT(AT+ λ0I)W = L2U2 (41)
be LU-factorizations of WT(A − λ0I)W and −WT(AT + λ0I)W , respectively. For this problem, (W − λ0I)−1 can be factored as
(W − λ0I)−1 =
W 0
0 W
F−1 0
0 I
I − ˜B ˜BT
0 I
F 0 0 D−1
×
I 0
− ˜CTC I˜
F−1 0
0 I
WT 0
0 WT
, (42)
−0.049 −0.048 −0.047
−4
−2 0 2 4x 10−4
target value
(2a) real parts
imaginary parts
−0.049 −0.048 −0.047
10−20 10−15 10−10 10−5
(2b) negative values of norms of eigenvalues
residuals
residuals by EM residuals by SISI(6) residuals by EM_SISI
Fig. 2. Eigenvalues and corresponding residuals.
where F = WT(A − λ0I)W , ˜B = WTB, ˜C = CW , and D = −WT(AT + λ0I)W − C˜TCF˜ −1B ˜˜BT. To evaluate D−1, we use the LU-factorizations in (41) combined with the Sherman–Morrison–Woodbury formula [11] to obtain
D−1= U2−1
I + L−12 C˜TS−1CU˜ 1−1L−11 B ˜˜BTU2−1 L−12 , (43)
where
S = I − ˜CU1−1L−11 B ˜˜BTU2−1L−12 C˜T. (44)
Example 2. With the data in [17] we get a system of size n = 8100. We chose B, C ∈ Rn×2 randomly with uniform distribution in the interval (0, 1) and used SHIRA to compute the 16 eigenvalues closest to λ0 = −0.048. The eigenvalues are depicted in Figure 2(a). The associated eigenvectors are then computed by EM, EM SISI, SISI(6), and IPI, respectively. The residuals for the eigenpairs computed by EM, EM SISI, and SISI(6), respectively, are shown in Figure 2(b).
In this example the IPI had convergence problems for the vectors xj associated with those eigenvalues with index j = 2, 3, 5, 6, 7, 8, 10, 11, 15, 16, because these com- plex conjugate pairs have relatively small imaginary parts and their use as shifts is not adequate. Here we see the loss of accuracy in the extraction method and we also observe that, as expected, this is easily compensated by adding a step of SISI.
Let us briefly compare the costs for computing an 5-dimensional invariant sub- space associated with stable eigenvalues for optimal control problems via EM SISI and SISI.
The matrix S in (44) is computed in SHIRA; it can be reused for SISI. It follows from the factorizations (42) and (43) that one iteration of SISI with re- orthogonalization requires 55 forward/backward substitutions by assuming that the
