Numerical behavior of stationary and two-step splitting iteration - PowerPoint PPT Presentation

Numerical behavior of stationary and two-step splitting iteration methods Miro Rozloˇ zn´ ık joint results with Zhong-zhi Bai and Pavel Jir´ anek Institute of Computer Science, Czech Academy of Sciences, Prague, Czech Republic Numerical Analysis and Scientific Computation with Applications (NASCA13), June 24-25-26, 2013, Universite du Littoral Cˆ ote d’Opale (ULCO), Calais, France

Delay of convergence and maximum attainable accuracy

Stationary iterative methods ◮ A x = b , A = M − N , G = M − 1 N , F = NM − 1 ◮ A: M x k +1 = N x k + b B: x k +1 = x k + M − 1 ( b − A x k ) ◮ Inexact solution of systems with M : every computed solution y of M y = z is interpreted as an exact solution of a system with perturbed data and relative perturbation bounded by parameter τ such that ( M + ∆ M ) y = z, � ∆ M� ≤ τ �M� , τk ( M ) ≪ 1 ◮ Higham, Knight 1993: M triangular, τ = O ( u )

Accuracy of the computed approximate solution A: M x k +1 = N x k + b ≤ τ �M − 1 |� ( �M� + �N� ) � ˆ x k +1 − x � max i =0 ,...,k {� ˆ x i �} � x � 1 − �G� � x � � b − A ˆ x k +1 � x k +1 � ≤ τ �M� � I − F� max i =0 ,...,k {� ˆ x i �} � b � + �A�� ˆ �A� 1 − �F� � x � x k +1 = x k + M − 1 ( b − A x k ) B: �M − 1 � ( �M� + �N� ) � ˆ x k +1 − x � max i =0 ,...,k {� ˆ x i �} ≤ O ( u ) � x � 1 − �G� − 2 τ �M − 1 ��M� � x � � b − A ˆ x k +1 � x k +1 � ≤ O ( u ) �M� + �N� � I − F� max i =0 ,...,k {� ˆ x i �} � b � + �A�� ˆ �A� 1 − �F� − 2 τ �M − 1 ��M� � x � Higham, Knight 1993, Bai, R, 2012

Numerical experiments: small model example A = tridiag(1 , 4 , 1) ∈ R 100 × 100 , b = A · ones(100 , 1) , κ ( A ) = � A � · � A − 1 � = 5 . 9990 · 0 . 4998 ≈ 2 . 9983 A = M − N , M = D − L, N = U

0 10 τ = 10 −2 −2 10 −4 10 τ = 10 −6 relative error norms −6 10 −8 10 τ = 10 −10 −10 10 −12 10 −14 10 τ = O(u) −16 10 0 5 10 15 20 25 30 35 40 45 50 iteration number k

0 10 τ = 10 −2 −2 10 −4 10 normwise backward error τ = 10 −6 −6 10 −8 10 τ = 10 −10 −10 10 −12 10 τ = O(u) −14 10 −16 10 0 5 10 15 20 25 30 35 40 45 50 iteration number k

0 10 −5 10 relative error norms τ = O(u), τ =10 −10 , τ =10 −6 , τ =10 −2 −10 10 −15 10 0 5 10 15 20 25 30 35 40 45 50 iteration number k

0 10 −2 10 −4 10 normwise backward errors −6 10 −8 10 τ = O(u), τ =10 −10 , τ =10 −6 , τ =10 −2 −10 10 −12 10 −14 10 −16 10 0 5 10 15 20 25 30 35 40 45 50 iteration number k

Two-step splitting iteration methods M 1 x k +1 / 2 = N 1 x k + b, A = M 1 − N 1 M 2 x k +1 = N 2 x k +1 / 2 + b, A = M 2 − N 2 Numerous solution schemes: Hermitian/skew-Hermitian (HSS) splitting, modified Hermitian/skew-Hermitian (MHSS) splitting, normal Hermitian/skew-Hermitian (NSS) splitting, preconditioned variant of modified Hermitian/skew-Hermitian (PMHSS) splitting and other splittings, ... Bai, Golub, Ng 2003, 2007, 2008; Bai 2009 Bai, Benzi, Chen 2010, 2011; Bai, Benzi, Chen, Wang 2012 � ˆ x k +1 − x � τ 1 �M − 1 2 N 2 ��M − 1 1 � ( �M 1 � + �N 1 � ) + τ 2 �M − 1 � � 2 � ( �M 2 � + �N 2 � ) � � x � max i =0 , 1 / 2 ,...,k +1 / 2 {� ˆ x i �} � x �

Two-step splitting iteration methods x k +1 / 2 = x k + M − 1 1 ( b − A x k ) x k +1 = x k +1 / 2 + M − 1 2 ( b − A x k +1 / 2 ) ⇔ x k +1 = x k + ( M − 1 + M − 1 − M − 1 2 AM − 1 1 )( b − A x k ) 1 2 = x k + ( I + M − 1 2 N 1 ) M − 1 1 ( b − A x k ) = x k + M − 1 2 ( I + N 2 M − 1 1 )( b − A x k ) � ˆ x k +1 − x � 1 ) � ( �M� + �N� )max i =0 ,...,k {� ˆ x i �} � O ( u ) �M − 1 2 ( I + N 2 M − 1 � x � � x �

Numerical experiments: small model example A = tridiag(2 , 4 , 1) ∈ R 100 × 100 , b = A · ones(100 , 1) , κ ( A ) = � A � · � A − 1 � = 5 . 9990 · 0 . 4998 ≈ 2 . 9983 H = 1 S = 1 2( A + A T ) , 2( A − A T ) A = H + S , H = tridiag(3 2 , 4 , 3 2) , S = tridiag(1 2 , 0 , − 1 2)

0 10 −2 10 τ 1 = 10 −2 , τ 2 =10 −10 −4 10 τ 1 = 10 −6 , τ 2 =10 −10 relative error norms −6 10 −8 10 τ 1 = 10 −10 , τ 2 =10 −10 −10 10 τ 1 = O(u), τ 2 =10 −10 −12 10 −14 10 τ 1 =O(u), 10 −10 , 10 −6 , 10 −2 , τ 2 = 10 −10 −16 10 0 10 20 30 40 50 60 70 80 90 100 iteration number k

0 10 −2 10 τ 1 = 10 −10 , τ 2 =10 −2 −4 10 τ 1 = 10 −10 , τ 2 =10 −6 relative error norms −6 10 −8 10 τ 1 = 10 −10 , τ 2 =10 −10 −10 10 τ 1 = 10 −10 , τ 2 =O(u) −12 10 −14 10 τ 1 = 10 −10 , τ 2 =O(u), 10 −10 , 10 −6 , 10 −2 −16 10 0 10 20 30 40 50 60 70 80 90 100 iteration number k

Saddle point problems We consider a saddle point problem with the symmetric 2 × 2 block form � A � � x � � f � B = . B T 0 y 0 ◮ A is a square n × n nonsingular (symmetric positive definite) matrix, ◮ B is a rectangular n × m matrix of (full column) rank m .

Schur complement reduction method ◮ Compute y as a solution of the Schur complement system B T A − 1 By = B T A − 1 f, ◮ compute x as a solution of Ax = f − By. ◮ Segregated vs. coupled approach: x k and y k approximate solutions to x and y , respectively. ◮ Inexact solution of systems with A : every computed solution ˆ u of Au = b is interpreted as an exact solution of a perturbed system ( A + ∆ A )ˆ u = b + ∆ b, � ∆ A � ≤ τ � A � , � ∆ b � ≤ τ � b � , τκ ( A ) ≪ 1 .

Iterative solution of the Schur complement system choose y 0 , solve Ax 0 = f − By 0    compute α k and p ( y )    k     y k +1 = y k + α k p ( y )   k     solve Ap ( x ) = − Bp ( y ) �    � k k    �    back-substitution:  outer �    �   �  iteration A: x k +1 = x k + α k p ( x )  k , inner �    �   iteration �  B: solve Ax k +1 = f − By k +1 ,  �     �     C: solve Au k = f − Ax k − By k +1 , �     �     �   x k +1 = x k + u k .   �     − α k B T p ( x ) r ( y ) k +1 = r ( y )   k k

Accuracy in the saddle point system � f − Ax k − By k � ≤ O ( α 1 ) κ ( A ) 1 − τκ ( A ) ( � f � + � B � Y k ) , k � ≤ O ( α 2 ) κ ( A ) � − B T x k − r ( y ) 1 − τκ ( A ) � A − 1 �� B � ( � f � + � B � Y k ) , Y k ≡ max {� y i � | i = 0 , 1 , . . . , k } . Back-substitution scheme α 1 α 2 A : Generic update τ u x k +1 = x k + α k p ( x ) k B : Direct substitution τ τ � x k +1 = A − 1 ( f − By k +1 ) additional C : Corrected dir. subst. system with A u τ x k +1 = x k + A − 1 ( f − Ax k − By k +1 ) − B T A − 1 f + B T A − 1 By k = − B T x k − B T A − 1 ( f − Ax k − By k )

Numerical experiments: a small model example A = tridiag(1 , 4 , 1) ∈ R 100 × 100 , B = rand(100 , 20) , f = rand(100 , 1) , κ ( A ) = � A � · � A − 1 � = 5 . 9990 · 0 . 4998 ≈ 2 . 9983 , κ ( B ) = � B � · � B † � = 7 . 1695 · 0 . 4603 ≈ 3 . 3001 . [R, Simoncini, 2002]

Generic update: x k +1 = x k + α k p ( x ) k 0 0 10 10 τ = 10 −2 −2 10 −2 10 (y) || −4 (y) ||/||r 0 10 −4 10 τ = 10 −6 relative residual norms ||−B T x k ||/||−B T x 0 ||, ||r k −6 −6 10 residual norm ||f−Ax k −By k || 10 −8 −8 10 10 τ = 10 −10 −10 10 −10 10 −12 10 −12 10 τ = O(u) −14 −14 10 10 τ = O(u), τ = 10 −10 , τ =10 −6 , τ =10 −2 −16 −16 10 10 −18 −18 10 10 0 50 100 150 200 250 300 0 50 100 150 200 250 300 iteration number k iteration number k

Direct substitution: x k +1 = A − 1 ( f − By k +1 ) 0 0 10 10 τ = 10 −2 −2 τ = 10 −2 10 −2 10 (y) || −4 (y) ||/||r 0 10 −4 10 τ = 10 −6 relative residual norms ||−B T x k ||/||−B T x 0 ||, ||r k τ = 10 −6 −6 −6 10 residual norm ||f−Ax k −By k || 10 −8 −8 10 10 τ = 10 −10 τ = 10 −10 −10 10 −10 10 −12 10 −12 10 τ = O(u) −14 −14 10 10 τ = O(u) −16 −16 10 10 −18 −18 10 10 0 50 100 150 200 250 300 0 50 100 150 200 250 300 iteration number k iteration number k

Corrected direct substitution: x k +1 = x k + A − 1 ( f − Ax k − By k +1 ) 0 0 10 10 −2 τ = 10 −2 10 −2 10 (y) || −4 (y) ||/||r 0 10 −4 10 relative residual norms ||−B T x k ||/||−B T x 0 ||, ||r k τ = 10 −6 −6 −6 10 residual norm ||f−Ax k −By k || 10 −8 −8 10 10 τ = 10 −10 −10 10 −10 10 −12 10 −12 10 −14 −14 10 10 τ = O(u), τ = 10 −10 , τ =10 −6 , τ =10 −2 τ = O(u) −16 −16 10 10 −18 −18 10 10 0 50 100 150 200 250 300 0 50 100 150 200 250 300 iteration number k iteration number k