A compact Arnoldi algorithm for polynomial eigenvalue problems Yangfeng Su, Junyi Zhang, Zhaojun Bai School of Mathematical Sciences Fudan University January, 2008, Taiwan RANMEP2008
Goals Polynomial Eigenvalue Problems (PEPs): � � A 0 + λ A 1 + ... + λ d A d x = 0 ( λ, x ) is called an eigenpair. d=1, SEP,GEP d=2, QEP ( Quadratic Eigenvalue Problem) Goals Solving large scale polynomial eigenvalue problems with Implicitly Restarted Arnoldi (IRA) algorithm with less memory.
Goals Polynomial Eigenvalue Problems (PEPs): � � A 0 + λ A 1 + ... + λ d A d x = 0 ( λ, x ) is called an eigenpair. d=1, SEP,GEP d=2, QEP ( Quadratic Eigenvalue Problem) Goals Solving large scale polynomial eigenvalue problems with Implicitly Restarted Arnoldi (IRA) algorithm with less memory.
Outline 1 Background 2 Algorithm 3 Numerical Comparison
Outline 1 Background 2 Algorithm 3 Numerical Comparison
Linearization Linearization ⇒ A larger GEP: − A 1 − A 2 . . . − A d A 0 x I I λ x − λ = 0 . ... ... . . λ d − 1 x 0 I I Sleijpen, van der Vorst, and van Gijzen. Quadratic eigenproblems are no problem, SIAM News , 1996 Can be solved with ARPACK Not true: structure is neither used nor preserved. size and therefore memory are d times! This work: explore the structure to save memory in ARPACK. It is a byproduct of SOAR [Bai & S. SIMAX 2005].
Linearization Linearization ⇒ A larger GEP: − A 1 − A 2 . . . − A d A 0 x I I λ x − λ = 0 . ... ... . . λ d − 1 x 0 I I Sleijpen, van der Vorst, and van Gijzen. Quadratic eigenproblems are no problem, SIAM News , 1996 Can be solved with ARPACK Not true: structure is neither used nor preserved. size and therefore memory are d times! This work: explore the structure to save memory in ARPACK. It is a byproduct of SOAR [Bai & S. SIMAX 2005].
Linearization Linearization ⇒ A larger GEP: − A 1 − A 2 . . . − A d A 0 x I I λ x − λ = 0 . ... ... . . λ d − 1 x 0 I I Sleijpen, van der Vorst, and van Gijzen. Quadratic eigenproblems are no problem, SIAM News , 1996 Can be solved with ARPACK Not true: structure is neither used nor preserved. size and therefore memory are d times! This work: explore the structure to save memory in ARPACK. It is a byproduct of SOAR [Bai & S. SIMAX 2005].
Linearization Linearization ⇒ A larger GEP: − A 1 − A 2 . . . − A d A 0 x I I λ x − λ = 0 . ... ... . . λ d − 1 x 0 I I Sleijpen, van der Vorst, and van Gijzen. Quadratic eigenproblems are no problem, SIAM News , 1996 Can be solved with ARPACK Not true: structure is neither used nor preserved. size and therefore memory are d times! This work: explore the structure to save memory in ARPACK. It is a byproduct of SOAR [Bai & S. SIMAX 2005].
Arnoldi decomposition An Arnoldi decomposition of order- j : OP Q j = Q j H j + h j +1 , j q j +1 e T j where OP : a matrix Q j = [ q 1 , q 2 , . . . , q j ]: orthonormal H j : upper Hessenberg Arnoldi process is used to compute the Arnoldi decomposition with numerical stability
Implicitly Restarted Arnoldi (IRA) for SEP Sorensen [SIMAX92] Given an Arnoldi decomposition of order- p , OP Q p = Q p H p + β p q p +1 e T p . 1 Extend Arnoldi decomposition from order p to order k : AQ k = Q k H k + β k q k +1 e T k . 2 Divide the eigenvalues of H k as “good” ones µ 1 , . . . , µ p and “bad” ones µ p +1 , . . . , µ k . 3 For H k do implicitly QR steps with shifts µ k +1 , . . . , µ k , get H k = U ˜ H k U H 4 Take first p columns of AQ k U = Q k U ˜ H k + β k q k +1 e T k U as a restarted Arnoldi decomposition of order p .
ARPACK ARPACK is an implementation of IRA algorithm a well-coded, well-documented package produced by Lehoucq, Sorensen and Yang during 1992-1997 used in MATLAB as eigs and arpackc
IRA for QEP For simplicity, we only discuss QEPs. For QEP: ( λ 2 M + λ D + K ) x = 0 1 shift-and-invert: for shift σ , let λ = σ + 1 /µ ( µ 2 I − µ A − B ) x = 0 where − ( σ 2 M + σ D + K ) − 1 (2 σ M + D ) A = − ( σ 2 M + σ D + K ) − 1 M = B 2 linearize � A � � µ x � � µ x � B = µ I 0 x x 3 apply IRA
IRA for QEP Easy use of ARPACK How to utilize the Frobenius structure to save memory?
Outline 1 Background 2 Algorithm 3 Numerical Comparison
Arnoldi decomposition for QEP An Arnoldi decomposition with order- j OP Q j = Q j +1 � H j For QEP: � A � � Q j , 1 � � Q j +1 , 1 � B � = H j I 0 Q j , 2 Q j +1 , 2 Since Q j , 1 = Q j +1 , 2 � H j we have Theorem rank ([ Q j , 1 , Q j , 2 ]) ≡ r j ≤ j + 1 Observed by many people, e.g. Meerbergen [SIMAX06], Bai & S. [SIMAX05]. The key is how to use it with numerical stability.
Arnoldi Decomposition for QEP Theorem rank ([ Q j , 1 , Q j , 2 ]) ≡ r j ≤ j + 1 Let V j ∈ C n × r j = orth [ Q j , 1 , Q j , 2 ] then � Q j , 1 � � V j R j , 1 � � V j � � R j , 1 � Q j = = = Q j , 2 V j R j , 2 V j R j , 2 Two levels of orthonormality: V j is orthonormal � R j , 1 � is orthonormal R j , 2
Compact ARnoldi Decomposition (CARD) Compact ARnoldi Decomposition (CARD) � V j � � R j , 1 � � V j +1 � � R j +1 , 1 � � = OP H j V j R j , 2 V j +1 R j +1 , 2 V j ∈ C n × r j R j ∈ C r j × j j ≤ r j +1 ≤ j + 1 Memory cost: Arnoldi: 2 n ( j + 1) ( for PEPs: dn ( j + 1)) CARD: nr j +1 ≤ n ( j + 2) (for PEPs: ≤ n ( j + d + 1))
CARD process CARD process is to compute the CARD with numerical stability! CARD of order j : � VR 1 � � VR 1 � � ˆ � V ˆ R 1 ˆ = H + β qe 1 = OP H V ˆ ˆ VR 2 VR 2 R 2 Expand it to a CARD of order j + 1 ⇒ next two pages
Expand CARD process 1 compute q 1 = Aq 1 + Bq 2 ; 2 decompose q 1 = ˆ V x + v α with MGS ˆ V T q 1 , x = q 1 − ˆ = V x , v α = � v � , = v /α v 3 update � � ˆ ˆ V = V , v , r j +1 = r j + 1 , � ˆ � ˆ � � ˆ R 2 (: , j + 1) R 1 x R 2 ˆ , ˆ R 1 = R 2 = 0 α 0 0
Expand CARD process ˆ R 1 x � ˆ � R 1 0 α := ˆ ˆ ˆ R 1 (: , j + 1) R 2 R 2 0 0 4 decompose with MGS: � ˆ � � ˆ � R 1 (: , j + 2) R 1 (: , 1 : j + 1) = H 1: j +1 , j +1 ˆ ˆ R 2 (: , j + 2) R 2 (: , 1 : j + 1) old � ˆ � R 1 (: , j + 2) + H j +2 , j +1 ˆ R 2 (: , j + 2) new
5 update the current Arnoldi vector q : � q 1 � q = q 2 q 1 = ˆ V ˆ R j +1 , 1 [: , j + 2] q 2 = ˆ V ˆ R j +1 , 2 [: , j + 2] Only GMS (with re-orthogonalization), no inversion � CARD is numerically stable!
IRA with CARD Given a CARD of k -order: � ˆ � VR 1 � � � � V ˆ R 1 H OP = V ˆ ˆ β e T VR 2 R 2 k IRA does ( m − p ) QR steps on H with shifts µ p +1 , ..., µ m , i.e. H = U � HU H Then � VR 1 � � ˆ � � � V ˆ U ˜ R 1 H U = OP V ˆ ˆ β e T VR 2 k U R 2 Denote � U � ˆ U = 1
IRA with CARD Then � ˆ � VR 1 U � � � � V ˆ R 1 ˆ ˜ U H OP = V ˆ ˆ R 2 ˆ β e T VR 2 U k U U Its first p columns , denoted by � V k R 1 , p � � V k +1 R 1 , p +1 � � H p � = OP ˜ β e T V k R 2 , p V k +1 R 2 , p +1 p still form an Arnoldi decomposition of order p However, the V k has r k (instead of r p ) columns, it is not a CARD!
IRA with CARD Since � V k +1 R 1 , p +1 � V k +1 R 2 , p +1 is the orthonormal factor of an Arnoldi decomposition, from previous theorem, rank [ V k +1 R 1 , p +1 , V k +1 R 2 , p +1 ] = rank [ R 1 , p +1 , R 2 , p +1 ] = r p +1 ≤ p + 2 , we have a compact SVD: P r k +1 × r p +1 Σ r p +1 × r p +1 [ G r p +1 × ( p +1) , G r p +1 × ( p +1) [ R 1 , p +1 , R 2 , p +1 ] = ] 1 2 ≡ P [ R 1 , R 2 ]
IRA with CARD Therefore, � � � ( V k +1 P ) n × r p +1 R 1 � V n × r k +1 R 1 , p +1 k +1 = ( V k +1 P ) R 2 V k +1 R 2 , p +1 The Arnoldi decomposition is expressed in CARD again! This process can also be implemented by a compact “QR” decomposition, which is similar with the compact Arnoldi decomposition (CARD). Details omitted.
POLYAR POLYAR: modified ARPARK for polynomial eigenvalue problems (not only QEPs) 1 znaitr p : compute CARD 2 znapps p : IRA with CARD; use LAPACK routine zgesdd to compute SVD decomposition 3 znaupd p , znaupd2 p , zgetv0 p : slightly revised (arguments, storage) 4 zgemip(added) : compute inner product in compact form
Outline 1 Background 2 Algorithm 3 Numerical Comparison
Recommend
More recommend
