Absolute value preconditioning for symmetric linear systems and - PowerPoint PPT Presentation

Absolute value preconditioning for symmetric linear systems and eigenvalue problems Eugene Vecharynski (joint work with Andrew Knyazev) Georgia State University Department of Mathematics and Statistics Scientific Computing Group Meeting Emory University, Atlanta, GA October 17, 2012 1 / 46

Two problems Symmetric indefinite linear system 1 Ax = b , A = A ∗ Mixed FE discretizations of PDEs in fluid and solid mechanics Discretizations of PDEs in acoustics Inner steps of interior point methods in optimization Inner steps of symmetric eigensolvers (e.g., Jacobi-Davidson) Symmetric interior eigenvalue problem 2 K = K ∗ , M = M ∗ > 0 Kv = λ Mv , Vibration modes analysis Electronic structure calculations 2 / 46

Ax = b and Kv = λ Mv General setting Large size, ill-conditioned Sparse matrices or matrix-free environment Direct methods are inapplicable. Iterate! Use preconditioners to improve convergence 3 / 46

Ax = b and Kv = λ Mv General setting Large size, ill-conditioned Sparse matrices or matrix-free environment Direct methods are inapplicable. Iterate! Use preconditioners to improve convergence ⇒ This talk concerns the construction of symmetric positive definite (SPD) preconditioners for both problems 4 / 46

Preconditioning Traditionally, For Ax = b , define T ≈ A − 1 For Kv = λ Mv , define T ≈ ( K − c 2 M ) − 1 ⇒ T is not SPD! 5 / 46

Outline In this talk ... Advocate SPD preconditioning Motivate and introduce absolute value (AV) preconditioning Construct an AV preconditioner for a model linear system and eigenvalue problem 6 / 46

AV preconditioning Symmetric indefinite linear systems 7 / 46

How to solve Ax = b , A = A ∗ ? Preconditioned system TAx = Tb 8 / 46

How to solve Ax = b , A = A ∗ ? Preconditioned system TAx = Tb T is symmetric indefinite or nonsymmetric (e.g., T ≈ A − 1 ) ⇒ TA is nonsymmetric in any inner product. The symmetry is lost. Methods: GMRES, BiCG, BiCGstab, QMR, etc. 9 / 46

How to solve Ax = b , A = A ∗ ? Preconditioned system TAx = Tb T is symmetric indefinite or nonsymmetric (e.g., T ≈ A − 1 ) ⇒ TA is nonsymmetric in any inner product. The symmetry is lost. Methods: GMRES, BiCG, BiCGstab, QMR, etc. T is SPD . ⇒ TA is symmetric in the T − 1 -based inner product ( · , · ) T − 1 = ( · , T − 1 · ) The symmetry is preserved. Method: PMINRES ( optimal, short-term recurrent ). 10 / 46

How to solve Ax = b , A = A ∗ ? Preconditioned system TAx = Tb T is symmetric indefinite or nonsymmetric (e.g., T ≈ A − 1 ) ⇒ TA is nonsymmetric in any inner product. The symmetry is lost. Methods: GMRES, BiCG, BiCGstab, QMR, etc. T is SPD . ⇒ TA is symmetric in the T − 1 -based inner product ( · , · ) T − 1 = ( · , T − 1 · ) The symmetry is preserved. Method: PMINRES ( optimal, short-term recurrent ). ⇒ How to correctly define an SPD preconditioner T ? 11 / 46

Matrix absolute value The eigenvalue decomposition A = V Λ V ∗ The matrix absolute value of A is defined as � � | A | = V | Λ | V ∗ , | Λ | = diag { � λ j � } . The matrix sign of A is defined as � � sign ( A ) = V sign (Λ) V ∗ , sign (Λ) = diag sign ( λ j ) . The polar decomposition A = | A | sign ( A ) = sign ( A ) | A | . 12 / 46

An ideal SPD preconditioner Let T = | A | − 1 . Then TAx = Tb is sign ( A ) x = | A | − 1 b . The matrix TA = sign ( A ) has two distinct eigenvalues: − 1 and 1. ⇒ PMINRES converges in at most two steps 13 / 46

An ideal SPD preconditioner Let T = | A | − 1 . Then TAx = Tb is sign ( A ) x = | A | − 1 b . The matrix TA = sign ( A ) has two distinct eigenvalues: − 1 and 1. ⇒ PMINRES converges in at most two steps T = | A | − 1 is an ideal SPD preconditioner Construction of the exact | A | − 1 is prohibitively expensive AV preconditioning: Construct an SPD preconditioner T as an approximation to | A | − 1 , i.e., T ≈ | A | − 1 . 14 / 46

Bounds on eigenvalues of TA Theorem Given a symmetric indefinite A, an SPD T, and constants δ 1 ≥ δ 0 > 0 , such that δ 0 ( v , T − 1 v ) ≤ ( v , | A | v ) ≤ δ 1 ( v , T − 1 v ) , ∀ v . Then � Λ( TA ) ⊂ [ − δ 1 , − δ 0 ] [ δ 0 , δ 1 ] . 15 / 46

Challenges For practical applications : Can we construct T ≈ | A | − 1 using A (matvecs with A )? Can we construct T ≈ | A | − 1 at an optimal cost O (nnz( A ))? 16 / 46

AV preconditioning Eigenvalue problems 17 / 46

Kv = λ Mv 18 / 46

Kv = λ Mv 19 / 46

Kv = λ Mv LOBPCG: ( V ( k ) ) ∗ KV ( k ) , ( V ( k ) ) ∗ MV ( k ) = I Λ ( k ) = � KV ( k ) − MV ( k ) Λ ( k ) � V ( k ) C ( k ) C ( k ) + P ( k ) C ( k ) V ( k + 1 ) = + T V R P V ( k + 1 ) − V ( k ) C ( k ) P ( k + 1 ) = V 20 / 46

Kv = λ Mv LOBPCG: ( V ( k ) ) ∗ KV ( k ) , ( V ( k ) ) ∗ MV ( k ) = I Λ ( k ) = � KV ( k ) − MV ( k ) Λ ( k ) � V ( k ) C ( k ) C ( k ) + P ( k ) C ( k ) V ( k + 1 ) = + T V R P V ( k + 1 ) − V ( k ) C ( k ) P ( k + 1 ) = V Choice of the preconditioner T ≈ K − 1 ⇒ emphasizes the extreme (smallest) eigenpairs � − 1 ⇒ indefinite preconditioning � K − c 2 M T ≈ � � � − 1 � K − c 2 M Our approach: T ≈ 21 / 46

The same challenges! For practical applications : Can we construct T ≈ | K − c 2 M | − 1 using K and M ? Can we construct T ≈ | K − c 2 M | − 1 at an optimal cost O (nnz( K )+nnz( M ))? 22 / 46

Model problems The “shifted Laplacian” equation − ∆ u ( x , y ) − c 2 u ( x , y ) = f ( x , y ) , ( x , y ) ∈ Ω = ( 0 , 1 ) × ( 0 , 1 ) u | Γ = 0 . The eigenvalue problem (find eigenpairs near c 2 ) − ∆ u ( x , y ) = λ u ( x , y ) , ( x , y ) ∈ Ω = ( 0 , 1 ) × ( 0 , 1 ) u | Γ = 0 . 23 / 46

Model problems After discretization (FD): Symmetric indefinite linear system ( L − c 2 I ) x = b ⇒ T ≈ | L − c 2 I | − 1 Symmetric eigenvalue problem (find eigenpairs near c 2 ). Lv = λ v ⇒ T ≈ | L − c 2 I | − 1 ⇒ Apply T by approximately solving | L − c 2 I | w = r . Use MG! 24 / 46

MG 25 / 46

Two-grid scheme for | L − c 2 I | w = r Algorithm (Two-grid scheme for | L − c 2 I | w = r ) Input r. Output w = Tr. Presmoothing . Apply ν smoothing steps, 1 w ( i + 1 ) = w ( i ) + τ ( r − | L − c 2 I | w ( i ) ) , i = 0 , . . . , ν − 1 , w ( 0 ) = 0 . Set w pre = w ( ν ) , ν ≥ 1 . Coarse grid correction . Restrict residual to the coarse grid (R), apply 2 � − 1 , prolongate to the fine grid (P), and add to w pre , � � L H − c 2 I H � w cgc = w pre + P | L H − c 2 I H | − 1 R � r − | L − c 2 I | w pre � . Postsmoothing . Apply ν smoothing steps, 3 w ( i + 1 ) = w ( i ) + τ ( r − | L − c 2 I | w ( i ) ) , i = 0 , . . . , ν − 1 , w ( 0 ) = w cgc . Return w = w ( ν ) . 26 / 46

Two-grid scheme for | L − c 2 I | w = r Algorithm (Two-grid scheme for | L − c 2 I | w = r ) Input r. Output w = Tr. Presmoothing . Apply ν smoothing steps, 1 w ( i + 1 ) = w ( i ) + τ ( r − | L − c 2 I | w ( i ) ) , i = 0 , . . . , ν − 1 , w ( 0 ) = 0 . Set w pre = w ( ν ) , ν ≥ 1 . Coarse grid correction . Restrict residual to the coarse grid (R), apply 2 � − 1 , prolongate to the fine grid (P), and add to w pre , � � L H − c 2 I H � w cgc = w pre + P | L H − c 2 I H | − 1 R � r − | L − c 2 I | w pre � . Postsmoothing . Apply ν smoothing steps, 3 w ( i + 1 ) = w ( i ) + τ ( r − | L − c 2 I | w ( i ) ) , i = 0 , . . . , ν − 1 , w ( 0 ) = w cgc . Return w = w ( ν ) . Idea: Replace | L − c 2 I | ≈ B , where y = Bu is easy to compute 27 / 46

Two-grid AV preconditioner Algorithm (Two-grid AV preconditioner) Input r, B ≈ | L − c 2 I | . Output w = Tr. Presmoothing . Apply ν smoothing steps, 1 w ( i + 1 ) = w ( i ) + τ ( r − Bw ( i ) ) , i = 0 , . . . , ν − 1 , w ( 0 ) = 0 . Set w pre = w ( ν ) , ν ≥ 1 . Coarse grid correction . Restrict residual to the coarse grid (R), apply 2 � − 1 , prolongate to the fine grid (P), and add to w pre , � � L H − c 2 I H � w cgc = w pre + P | L H − c 2 I H | − 1 R r − Bw pre � � . Postsmoothing . Apply ν smoothing steps, 3 w ( i + 1 ) = w ( i ) + τ ( r − Bw ( i ) ) , i = 0 , . . . , ν − 1 , w ( 0 ) = w cgc . Return w = w ( ν ) . 28 / 46

MG AV preconditioner (V-cycle) Algorithm (AV-MG( r l ): the MG AV preconditioner) Input r l , output w l . Set B l ≈ | L l − c 2 I l | . 1 Presmoothing . Apply ν smoothing steps, 2 w ( i + 1 ) = w ( i ) + τ l ( r l − B l w ( i ) l ) , i = 0 , . . . , ν − 1 , w ( 0 ) = 0 . l l l = w ( ν ) Set w pre , ν ≥ 1 . l l Coarse grid correction . Restrict residual to coarser grid (R l − 1 ), apply 3 AV-MG on this level, prolongate to the finer grid (P l ), add to w pre , l w cgc w pre r l − B l w pre � � = + P l AV-MG ( R l − 1 ) , if l > 1 ; l l l w cgc w pre + P 1 | L 0 − c 2 I 0 | − 1 R 0 r 1 − B 1 w pre � � = , if l = 1 . 1 1 1 Postsmoothing . Apply ν smoothing steps, 4 w ( i + 1 ) = w ( i ) + τ l ( r l − B l w ( i ) l ) , i = 0 , . . . , ν − 1 , w ( 0 ) = w cgc . l l l l Return w l = w ( ν ) . l 29 / 46