Projected Krylov Methods for Unsymmetric Augmented Systems - PowerPoint PPT Presentation

Projected Krylov Methods for Unsymmetric Augmented Systems Dominique Orban dominique.orban@gerad.ca GERAD and ´ Ecole Polytechnique de Montr´ eal, Canada CERFACS Sparse Days 2008

Outline Problem Statement and Assumptions Direct Application of Iterative Methods Nullspace Methods Krylov Menu Krylov Methods in the Nullspace Projected Krylov Methods Numerical Experience

Problem Statement and Context � � u � b B T � A � � = 0 B p d ◮ Optimization and Control, incl. least squares ( A symmetric), ◮ Multi immiscible fluid flow with free surface, ◮ Fictitious domains method for flow around obstacle, ◮ Electromagnetism, image restoration, . . . Denote the augmented matrix K ( A ). ◮ Do not want to assemble A , ◮ B “flat” and contributes marginally to density of K ( A ), ◮ Can assemble B , ◮ Can compute Ap but not A T p .

Typical Block Structure Density( K ( A )): 0 . 315%, A accounts for 0 . 279%, B for 0 . 036%

Assumptions and Basic Results Do not assume A symmetric or positive (semi-)definite, but: 1. B full row rank, 2. N ( A ) ∩ N ( B ) = { 0 } , 3. R ( A | N ( B )) ∩ R ( B T ) = { 0 } . where R ( A | N ( B )) = { v = Aw for some w ∈ N ( B ) } Theorem K ( A ) is nonsingular ⇐ ⇒ assumptions 1–3 Theorem (Preconditioners) G = G T positive definite over N ( B ) = ⇒ K ( G ) nonsingular

Direct Application of Iterative Methods System size: n + m = 7761, itmax = 6( n + m ), no preconditioner Solver Iterations Residual Workspace Time (s) BCG 46566 0 . 1329880e+02 54327 21 . 95 DBCG 46567 0 . 4489614e+01 85371 23 . 02 0 . 1537764e − 04 CGNR 46567 38805 21 . 25 BCGSTAB 46567 0 . 3111281e − 04 62088 21 . 95 0 . 8147863e − 04 TFQMR 46567 85371 23 . 11 GMRES 46566 0 . 6374488e+97 132108 28 . 33 FGMRES 46566 0 . 3499084e+96 248519 28 . 63 DQGMRES 46566 0 . 1387398e+95 256177 44 . 45 F77 implementations from SPARSKIT [Saad]

Nullspace Methods Au + B T p = b , and Bu = d . Let Z = orthonormal basis for N ( B ). Any solution u ∗ has the form u ∗ = Zu ∗ Z + B T u ∗ B .

Nullspace Methods Au + B T p = b , and Bu = d . Let Z = orthonormal basis for N ( B ). Any solution u ∗ has the form u ∗ = Zu ∗ Z + B T u ∗ B . Then d = Bu ∗ = BB T u ∗ ← unique solution , B

Nullspace Methods Au + B T p = b , and Bu = d . Let Z = orthonormal basis for N ( B ). Any solution u ∗ has the form u ∗ = Zu ∗ Z + B T u ∗ B . Then d = Bu ∗ = BB T u ∗ ← unique solution , B and Z T AZ u ∗ Z = Z T ( b − AB T u ∗ B ) ← concentrate on this system .

Nullspace Methods Au + B T p = b , and Bu = d . Let Z = orthonormal basis for N ( B ). Any solution u ∗ has the form u ∗ = Zu ∗ Z + B T u ∗ B . Then d = Bu ∗ = BB T u ∗ ← unique solution , B and Z T AZ u ∗ Z = Z T ( b − AB T u ∗ B ) ← concentrate on this system . Worry about p ∗ later. Problem: Computing Z is too costly. Worry about it in a minute.

Krylov Menu ◮ K ( A ) is always indefinite (unless B vanishes): eliminate CG, ◮ K ( A ) is not symmetric: eliminate MINRES and SYMMLQ, ◮ Do not want to compute A T p : eliminate Bi-CG and QMR. Leaves CGNE, CGNR, CGS, TFQMR, Bi-CGSTAB, GMRES( m ). ◮ CGNE and CGNR “square the conditioning”, ◮ CGS exhibits erratic convergence, ◮ GMRES is quite memory hungry. At this point, concentrate efforts on Bi-CGSTAB and TFQMR.

Preconditioned Bi-CGSTAB for R − T MR − 1 2 x = R − T y 1 1 1. x 0 ∈ R n , r 0 = b − Mx 0 , ¯ r T 0 r 0 � = 0, d 0 = r 0 , k = 0. r T ¯ 0 r k 2. Solve C ¯ d k = d k , set α k = , 0 M ¯ r T ¯ d k 3. s k = r k − α k M ¯ d k , solve C ¯ s k = s k , 4. ω k = ( R − T s k ) T ( R − T M ¯ s k ) 1 1 , � R − T s k � 2 M ¯ 1 2 5. x k +1 = x k + α k ¯ d k + ω k ¯ s k , 6. r k +1 = s k − ω k M ¯ s k , r T ¯ 7. β k = α k 0 r k +1 , r T ω k ¯ 0 r k 8. d k +1 = r k +1 + β k ( d k − ω k M ¯ d k ). C = R T 1 R 2 is the preconditioner.

Preconditioned Bi-CGSTAB for R − T MR − 1 2 x = R − T y 1 1 1. x 0 ∈ R n , r 0 = b − Mx 0 , ¯ r T 0 r 0 � = 0, d 0 = r 0 , k = 0. r T ¯ 0 r k 2. Solve C ¯ d k = d k , set α k = , 0 M ¯ r T ¯ d k 3. s k = r k − α k M ¯ d k , solve C ¯ s k = s k , 4. ω k = ( R − T s k ) T ( R − T ≈ s T k M ¯ M ¯ s k ) s k 1 1 , � R − T s k � 2 s k � 2 � M ¯ M ¯ 2 1 2 5. x k +1 = x k + α k ¯ d k + ω k ¯ s k , 6. r k +1 = s k − ω k M ¯ s k , r T ¯ 7. β k = α k 0 r k +1 , r T ω k ¯ 0 r k 8. d k +1 = r k +1 + β k ( d k − ω k M ¯ d k ). C = R T 1 R 2 is the preconditioner.

Apply Bi-CGSTAB to Nullspace System Z T AZ u = Z T ( b − AB T u ∗ C = Z T GZ B ) Choose an initial u Z 0 .

Apply Bi-CGSTAB to Nullspace System Z T AZ u = Z T ( b − AB T u ∗ C = Z T GZ B ) Choose an initial u Z 0 . The initial residual is Z T ( b − AB T u ∗ B ) − Z T AZ u Z r Z = 0 0 Z T � � 0 + B T u ∗ b − A ( Zu Z = B ) Z T ( b − Au 0 ) = Z T r 0 . =

Apply Bi-CGSTAB to Nullspace System Z T AZ u = Z T ( b − AB T u ∗ C = Z T GZ B ) Choose an initial u Z 0 . The initial residual is Z T ( b − AB T u ∗ B ) − Z T AZ u Z r Z = 0 0 Z T � � 0 + B T u ∗ b − A ( Zu Z = B ) Z T ( b − Au 0 ) = Z T r 0 . = k = Z T d k , Z ¯ k = ¯ k = Z T r k , d Z Similarly, r Z d Z d k , so Z T GZ ¯ d Z k = d Z k

Apply Bi-CGSTAB to Nullspace System Z T AZ u = Z T ( b − AB T u ∗ C = Z T GZ B ) Choose an initial u Z 0 . The initial residual is Z T ( b − AB T u ∗ B ) − Z T AZ u Z r Z = 0 0 Z T � � 0 + B T u ∗ b − A ( Zu Z = B ) Z T ( b − Au 0 ) = Z T r 0 . = k = Z T d k , Z ¯ k = ¯ k = Z T r k , d Z Similarly, r Z d Z d k , so ¯ ( Z T GZ ) − 1 d Z k = d Z k

Apply Bi-CGSTAB to Nullspace System Z T AZ u = Z T ( b − AB T u ∗ C = Z T GZ B ) Choose an initial u Z 0 . The initial residual is Z T ( b − AB T u ∗ B ) − Z T AZ u Z r Z = 0 0 Z T � � 0 + B T u ∗ b − A ( Zu Z = B ) Z T ( b − Au 0 ) = Z T r 0 . = k = Z T d k , Z ¯ k = ¯ k = Z T r k , d Z Similarly, r Z d Z d k , so ¯ ( Z T GZ ) − 1 Z T d k d Z k =

Apply Bi-CGSTAB to Nullspace System Z T AZ u = Z T ( b − AB T u ∗ C = Z T GZ B ) Choose an initial u Z 0 . The initial residual is Z T ( b − AB T u ∗ B ) − Z T AZ u Z r Z = 0 0 Z T � � 0 + B T u ∗ b − A ( Zu Z = B ) Z T ( b − Au 0 ) = Z T r 0 . = k = Z T d k , Z ¯ k = ¯ k = Z T r k , d Z Similarly, r Z d Z d k , so ¯ k = Z ( Z T GZ ) − 1 Z T d k d Z Z

Apply Bi-CGSTAB to Nullspace System Z T AZ u = Z T ( b − AB T u ∗ C = Z T GZ B ) Choose an initial u Z 0 . The initial residual is Z T ( b − AB T u ∗ B ) − Z T AZ u Z r Z = 0 0 Z T � � 0 + B T u ∗ b − A ( Zu Z = B ) Z T ( b − Au 0 ) = Z T r 0 . = k = Z T d k , Z ¯ k = ¯ k = Z T r k , d Z Similarly, r Z d Z d k , so ¯ d k = Z ( Z T GZ ) − 1 Z T d k

Apply Bi-CGSTAB to Nullspace System Z T AZ u = Z T ( b − AB T u ∗ C = Z T GZ B ) Choose an initial u Z 0 . The initial residual is Z T ( b − AB T u ∗ B ) − Z T AZ u Z r Z = 0 0 Z T � � 0 + B T u ∗ b − A ( Zu Z = B ) Z T ( b − Au 0 ) = Z T r 0 . = k = Z T d k , Z ¯ k = ¯ k = Z T r k , d Z Similarly, r Z d Z d k , so ¯ d k = Z ( Z T GZ ) − 1 Z T d k i.e. ¯ d k = P G ( d k ) oblique projection onto N ( B )

Projected Preconditioned Bi-CGSTAB 1. u 0 ∈ R n , r 0 = b − Au 0 , ( Z ¯ r 0 ) T r 0 � = 0, d 0 = r 0 , k = 0. r 0 ) T r k ( Z ¯ 2. Compute ¯ d k = P G ( d k ), set α k = , r 0 ) T A ¯ ( Z ¯ d k 3. s k = r k − α k A ¯ d k , compute ¯ s k = P G ( s k ), 4. ω k = ( ZZ T s k ) T A ¯ s k , � Z T A ¯ s k � 2 2 5. u k +1 = u k + α k ¯ d k + ω k ¯ s k , 6. r k +1 = s k − ω k A ¯ s k , r 0 ) T r k +1 ( Z ¯ 7. β k = α k , r 0 ) T r k ω k ( Z ¯ 8. d k +1 = r k +1 + β k ( d k − ω k A ¯ d k ).

Projected Preconditioned Bi-CGSTAB 1. u 0 ∈ R n , r 0 = b − Au 0 , ( Z ¯ r 0 ) T r 0 � = 0, d 0 = r 0 , k = 0. r 0 ) T r k ( Z ¯ 2. Compute ¯ d k = P G ( d k ), set α k = , r 0 ) T A ¯ ( Z ¯ d k 3. s k = r k − α k A ¯ d k , compute ¯ s k = P G ( s k ), 4. ω k = ( ZZ T s k ) T A ¯ = P I ( s k ) T A ¯ s k s k , � Z T A ¯ s k � 2 s k ) � 2 � P I ( A ¯ 2 2 5. u k +1 = u k + α k ¯ d k + ω k ¯ s k , 6. r k +1 = s k − ω k A ¯ s k , r 0 ) T r k +1 ( Z ¯ 7. β k = α k , r 0 ) T r k ω k ( Z ¯ 8. d k +1 = r k +1 + β k ( d k − ω k A ¯ d k ).

Choice of Fixed Vector ◮ ¯ r 0 = Z T ˜ r 0 ∈ R n . Then r 0 for some ˜ r 0 ) T v = ( ZZ T ˜ r 0 ) T v = P I (˜ r 0 ) T v . ( Z ¯

Choice of Fixed Vector ◮ ¯ r 0 = Z T ˜ r 0 ∈ R n . Then r 0 for some ˜ r 0 ) T v = ( ZZ T ˜ r 0 ) T v = P I (˜ r 0 ) T v . ( Z ¯ ◮ ¯ r 0 = ( Z T GZ ) − 1 Z T ˜ r 0 ∈ R n . Then r 0 for some ˜ r 0 ) T v = P G (˜ r 0 ) T v . ( Z ¯

Choice of Fixed Vector ◮ ¯ r 0 = Z T ˜ r 0 ∈ R n . Then r 0 for some ˜ r 0 ) T v = ( ZZ T ˜ r 0 ) T v = P I (˜ r 0 ) T v . ( Z ¯ ◮ ¯ r 0 = ( Z T GZ ) − 1 Z T ˜ r 0 ∈ R n . Then r 0 for some ˜ r 0 ) T v = P G (˜ r 0 ) T v . ( Z ¯ r 0 = r 0 , we ask that r 0 �⊥ N ( B ). Upon picking ˜

Computing Projections � I � � ¯ � v B T � � v v = P I ( v ) ¯ ⇐ ⇒ = K ( I ) 0 0 B w � � ¯ � v B T � G � � v ⇐ ⇒ ¯ v = P G ( v ) = K ( G ) 0 0 B w Rely on symmetric indefinite factorization of K ( I ) and/or K ( G ).