IDR( s ) A family of simple and fast algorithms for solving large - PowerPoint PPT Presentation

IDR( s ) A family of simple and fast algorithms for solving large nonsymmetric systems of linear equations Harrachov 2007 Peter Sonneveld en Martin van Gijzen August 24, 2007 1 Delft University of Technology

Outline • Introduction • The Induced Dimension Reduction Theorem • The IDR( s ) algorithm • Numerical experiments • Conclusions 2 August 24, 2007

Product Bi-CG methods Bi-CG solves nonsymmetric linear systems using (short) CG recursions but needs extra matvec with A H . Idea of Sonneveld: use ’wasted’ matvec in a more useful way. Result: transpose-free methods: • CGS (Sonneveld, 1989) • Bi-CGSTAB (Van der Vorst, 1992) • BiCGSTAB2 (Gutknecht, 1993) • TFQMR (Freund, 1993) • BiCGstab( ℓ ) (Sleijpen and Fokkema, 1994) • Many other variants 3 August 24, 2007

Historical remarks Sonneveld first developed IDR (1980). Analysis showed that IDR was Bi-CG combined with linear minimal residual steps. The fact that IDR is transpose free, combined with the relation with Bi-CG led to the development of a now famous algorithm: CGS. Later Van der Vorst proposed another famous method: Bi-CGSTAB, which is mathematicaly equivalent with IDR. As a result of these developments, the basic IDR idea was abandoned for the Bi-CG approach. IDR is now forgotten. 4 August 24, 2007

The IDR idea The IDR-idea is to generate a sequence of subspaces G 0 · · · G j with the following operations: - Intersect G j − 1 with fixed subspace S , - Compute G j = ( I − ω j A )( G j − 1 ∩ S ) . The subspaces G 0 · · · G j are nested and of shrinking dimension. 5 August 24, 2007

The IDR theorem Theorem 1 (IDR) Let A be any matrix in C N × N , let v 0 be any nonzero vector in C N , and let G 0 be the complete Krylov space K N ( A , v 0 ) . Let S denote any (proper) subspace of C N such that S and G 0 do not share a nontrivial invariant subspace of A , and define the sequence G j , j = 1 , 2 , . . . as G j = ( I − ω j A )( G j − 1 ∩ S ) where the ω j ’s are nonzero scalars. Then (i) G j ⊂ G j − 1 for all j > 0 . (ii) G j = { 0 } for some j ≤ N . 6 August 24, 2007

Making an IDR algorithm The IDR theorem can be used to construct solution algorithms. This is done by constructing residuals r n ∈ G j . According to the IDR theorem ultimately r n ∈ G M = { 0 } . In order to generate residuals and corresponding solution approximations we first look at the basic recursions. 7 August 24, 2007

Krylov methods: basic recursion (1) A Krylov-type solver produces iterates x n , for which the residuals r n = b − Ax n are in the Krylov spaces K n ( A , r 0 ) = r 0 ⊕ Ar 0 ⊕ A 2 r 0 ⊕ · · · ⊕ A n r 0 , The next residual r n +1 can be generated by b j � r n +1 = − α Ar n + β j r n − j . j =0 The parameters α, β j determine the specific method and must be such that x n +1 can be computed. 8 August 24, 2007

Krylov methods: basic recursion (2) By using the difference vector ∆ r k = r k +1 − r k = − A ( x n +1 − x n ) , an explicit way to satisfy this requirement is b j � r n +1 = r n − α Ar n − γ j ∆ r n − j , j =1 which leads to the following update for the x estimate: b j � x n +1 = x n + α r n − γ j ∆ x n − j , j =1 9 August 24, 2007

Computation of a new residual (1) Residuals are computed that are forced to be in the subspaces G j , by application of the IDR-theorem. The residual r n +1 is in G j +1 if r n +1 = ( I − ω j +1 A ) v with v ∈ G j ∩ S . The main problem is to find v . 10 August 24, 2007

Computation of a new residual (2) 11 August 24, 2007

Computation of a new residual (3) The vector v is a combination of the residuals r l in G j . b j � v = r n − γ j ∆ r n − j . j =1 Let the space S be the left null space of some N × s matrix P : S = N ( P H ) . P = ( p 1 p 2 . . . p s ) , Since v is also in S = N ( P H ) , it must satisfy P H v = 0 . Combining these two yields an s × � j linear system for the coefficients γ j that (normally) is uniquely solvable if � j = s . 12 August 24, 2007

Computation of a new residual (4) Hence with the residual r n , and a matrix ∆ R consisting of the last s residual differences: ∆ R = (∆ r n − 1 ∆ r n − 2 . . . ∆ r n − s ) a suitable v can be found by ( P H ∆ R ) c = P H r n Solve s × s system v = r n − ∆ Rc Calculate 13 August 24, 2007

Building G j +1 (1) Assume r n and all columns of ∆ R are in G j , and let r n +1 be calculated as r n +1 = v − ω j +1 Av Then r n +1 ∈ G j +1 . Since G j +1 ⊂ G j (theorem) we automatically have r n +1 ∈ G j Now the next ∆ R is made by repeating the calculations. In this way we find s + 1 residuals in G j +1 14 August 24, 2007

Building G j +1 (2) 15 August 24, 2007

Building G j +1 (3) 16 August 24, 2007

Building G j +1 (4) 17 August 24, 2007

A few details 1. The first s + 1 residuals, starting with r 0 can be constructed by any Krylov-based iteration, such as a local minimum residual method. 2. In our actual implementation, all steps are identical. However, in calculating the first residual in G j +1 , a new value ω j +1 may be chosen. We choose ω j +1 such that � v − ω j +1 Av � is minimal. 18 August 24, 2007

Basic IDR( s ) algorithm. while � r n � > TOL or n < MAXIT do for k = 0 to s do Solve c from P H d R n c = P H r n v = r n − d R n c ; t = Av ; if k = 0 then ω = ( t H v ) / ( t H t ) ; end if d r n = − d R n c − ω t ; d x n = − d X n c + ω v ; r n +1 = r n + d r n ; x n +1 = x n + d x n ; n = n + 1; d R n = ( d r n − 1 · · · d r n − s ); d X n = ( d x n − 1 · · · d x n − s ) ; end for end while 19 August 24, 2007

Vector operations per MATVEC Method DOT AXPY Memory Requirements IDR(1) 2 4 8 2 2 5 5 IDR(2) 11 3 6 4 2 9 7 IDR(4) 17 5 10 6 2 13 9 IDR(6) 23 7 14 n +1 n +1 Full GMRES n + 2 2 2 BiCGSTAB 2 3 7 20 August 24, 2007

Relation with other methods Although the approach is different, IDR( s ) is closely related to some Bi-CGSTAB methods: • IDR(1) and Bi-CGSTAB yield the same residuals at the even steps. • ML( k )BiCGSTAB (Yeung and Chen, 1999) seems closely related to IDR( s ), BUT • IDR( s ) is MUCH simpler (both conceptually and its implementation) • Other, more natural extensions are possible, e.g. to avoid breakdown. 21 August 24, 2007

Performance of IDR( s ) The IDR theorem states that - it is possible to generate a sequence of nested subspace G j of shrinking dimension, - but does not say how fast the dimension shrinks It can be proven that the dimension reduction is (normally) s , So dim ( G j +1 ) = dim ( G j ) − s . IDR( s ) requires at at most N + N s matrix-vector multiplications to compute the exact solution. 22 August 24, 2007

Numerical experiments We will present two typical numerical examples • A 2D Ocean Circulation Problem • A 3D Helmholtz Problem 23 August 24, 2007

A 2D Ocean Circulation Problem We compare IDR( s ) with Full GMRES, restarted GMRES and Bi-CGSTAB. This ocean example is representative for a wide class of CFD problems. We will compare: • Rate of convergence • Stagnation level (of the true residual norm) 24 August 24, 2007

Stommel’s model for ocean circulation Balance between bottom friction, wind stress and Coriolis force. − r ∆ ψ − β ∂ψ ∂x − = ( ∇ × F ) z plus circulation condition around islands k � � r ∂ψ ∂n ds = − F · s ds. Γ k Γ k • ψ : streamfunction • r : bottom friction parameter • β : Coriolis parameter • F : Wind stress 25 August 24, 2007

Discretization of the ocean problem • Discretization with linear finite elements • Results in nonsymmetric system of 42248 • Eigenvalues are (almost) real Solution parameters: • ILU(0) preconditioning • P : s − 1 random vectors plus r 0 (for comparison with Bi-CGSTAB) 26 August 24, 2007

Solution of the ocean problem 80 60 40 20 0 −20 −40 −60 −80 0 50 100 150 200 250 300 350 27 August 24, 2007

Convergence for the ocean problem Convergence 2D ocean problem 5 10 IDR1 IDR2 IDR4 IDR6 0 10 BICGSTAB GMRES Scaled residual norm −5 10 −10 10 −15 10 −20 10 0 100 200 300 400 500 600 700 Number of MATVECS 28 August 24, 2007

Some observations • Required number of MATVECS decreases if s is increased. IDR(4) and IDR(6) are close to the optimal convergence curve of full GMRES. • Convergence curves of IDR(1) and Bi -CGSTAB coincide. • Stagnation levels of IDR( s ) comparable with Bi-CGSTAB. 29 August 24, 2007

Required number of MATVECS Method Number of MATVECS Vectors Full GMRES 265 268 GMRES(20) > 10000 23 GMRES(50) 4671 53 Bi -CGSTAB 411 7 IDR(1) 420 8 IDR(2) 339 11 IDR(4) 315 17 IDR(6) 307 23 Tolerance: � b − Ax n � < 10 − 8 � b � 30 August 24, 2007

A 3D Helmholtz Problem Example models sound propagation in a room of 4 × 4 × 4 m 3 . A harmonic sound source gives acoustic pressure field p ( x ) e 2 πift . p ( x , t ) = � The pressure function � p can be determined from − (2 πf ) 2 p − ∆ � p = δ ( x − x s ) in Ω . � c 2 in which - c : the sound speed (340 m/s) - δ ( x − x s ) : the harmonic point source, in the center of the room. 31 August 24, 2007