Inexact saddle point solvers and their limiting accuracy anek 1 , 2 , - - PowerPoint PPT Presentation

Inexact saddle point solvers and their limiting accuracy

Pavel Jir´ anek1,2, Miroslav Rozloˇ zn´ ık1,2

Faculty of Mechatronics and Interdisciplinary Engineering Studies, Technical University of Liberec, Czech Republic1 and Institute of Computer Science, Czech Academy of Sciences, Prague, Czech Republic2

Gene Golub Day at TU Berlin, February 29, 2008

1 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Workshop on Solution Methods for Saddle Point Systems, Hong Kong Baptist University, October 31, 2007

2 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Workshop on Solution Methods for Saddle Point Systems, Hong Kong Baptist University, October 31, 2007

3 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Workshop on Solution Methods for Saddle Point Systems, Hong Kong Baptist University, October 31, 2007

4 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Saddle point problems

We consider a saddle point problem with the symmetric 2 × 2 block form „ A B BT « „x y « = „f « . A is a square n × n nonsingular (symmetric positive definite) matrix, B is a rectangular n × m matrix of (full column) rank m. Applications: mixed finite element approximations, weighted least squares, constrained optimization etc. [Benzi, Golub, and Liesen, 2005].

5 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

inexact solutions of inner systems + rounding errors → inexact saddle point solver

6 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Schur complement reduction method

Compute y as a solution of the Schur complement system BT A−1By = BT A−1f, compute x as a solution of Ax = f − By. Systems with A are solved inexactly, the computed solution ¯ u of Au = b is interpreted an exact solution of a perturbed system (A + ∆A)¯ u = b + ∆b, ∆A ≤ τA, ∆b ≤ τb, τκ(A) ≪ 1.

7 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Iterative solution of the Schur complement system

choose y0, solve Ax0 = f − By0 compute αk and p(y)

k

yk+1 = yk + αkp(y)

k

˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ solve Ap(x)

k

= −Bp(y)

k

back-substitution: A: xk+1 = xk + αkp(x)

k ,

B: solve Axk+1 = f − Byk+1, C: solve Auk = f − Axk − Byk+1, xk+1 = xk + uk. 9 > > > > > > > > > = > > > > > > > > > ; inner iteration r(y)

k+1 = r(y) k

− αkBT p(x)

k

9 > > > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > > > ;

• uter

iteration

8 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Measure of the limiting accuracy

The limiting (maximum attainable) accuracy is measured by the ultimate (asymptotic) values of:

1

the Schur complement residual: BT A−1f − BT A−1Byk;

2

the residuals in the saddle point system: f − Axk − Byk and −BT xk;

3

the forward errors: x − xk and y − yk. Numerical example: A = tridiag(1, 4, 1) ∈ R100×100, B = rand(100, 20), f = rand(100, 1), κ(A) = A · A−1 = 7.1695 · 0.4603 ≈ 3.3001, κ(B) = B · B† = 5.9990 · 0.4998 ≈ 2.9983.

9 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Accuracy in the outer iteration process

BT (A + ∆A)−1Bˆ y = BT (A + ∆A)−1f, BT A−1f − BT A−1Bˆ y ≤ τκ(A) 1 − τκ(A)A−1B2ˆ y. − BT A−1f + BT A−1Byk − r(y)

k ≤ O(τ)κ(A)

1 − τκ(A)A−1B(f + BYk).

50 100 150 200 250 300 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number k relative residual norms ||BTA−1f−BTA−1Byk||/||BTA−1f−BTA−1By0||, ||r(y)

k ||/||r(y) 0 ||

10 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Accuracy in the saddle point system

−BT A−1f + BT A−1Byk = −BT xk − BT A−1(f − Axk − Byk) f − Axk − Byk ≤ O(α1)κ(A) 1 − τκ(A) (f + BYk), − BT xk − r(y)

k ≤ O(α2)κ(A)

1 − τκ(A) A−1B(f + BYk), Yk ≡ max{yi | i = 0, 1, . . . , k}. Back-substitution scheme α1 α2 A: Generic update xk+1 = xk + αkp(x)

k

τ u B: Direct substitution xk+1 = A−1(f − Byk+1) τ τ C: Corrected dir. subst. xk+1 = xk + A−1(f − Axk − Byk+1) u τ

• additional

system with A

11 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Generic update: xk+1 = xk + αkp(x)

k

50 100 150 200 250 300 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number k residual norm ||f−Axk−Byk|| 50 100 150 200 250 300 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u), τ = 10−10, τ=10−6, τ=10−2 iteration number k relative residual norms ||−BTxk||/||−BTx0||, ||rk

(y)||/||r0 (y)||

12 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Direct substitution: xk+1 = A−1(f − Byk+1)

50 100 150 200 250 300 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number k residual norm ||f−Axk−Byk|| 50 100 150 200 250 300 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number k relative residual norms ||−BTxk||/||−BTx0||, ||rk

(y)||/||r0 (y)||

13 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Corrected direct substitution: xk+1 = xk + A−1(f − Axk − Byk+1)

50 100 150 200 250 300 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u), τ = 10−10, τ=10−6, τ=10−2 iteration number k residual norm ||f−Axk−Byk|| 50 100 150 200 250 300 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number k relative residual norms ||−BTxk||/||−BTx0||, ||rk

(y)||/||r0 (y)||

14 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Forward error of computed approximate solution

x − xk ≤ γ1f − Axk − Byk + γ2 − BT xk, y − yk ≤ γ2f − Axk − Byk + γ3 − BT xk, γ1 = σ−1

min(A), γ2 = σ−1 min(B), γ3 = σ−1 min(BT A−1B).

50 100 150 200 250 300 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number k relative error norms ||x−xk||A/||x−x0||A, ||y−yk||B

TA −1B/||y−y0||B TA −1B

15 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Null-space projection method

compute x ∈ N(BT ) as a solution of the projected system (I − Π)A(I − Π)x = (I − Π)f, compute y as a solution of the least squares problem By ≈ f − Ax, Π is the orthogonal projector onto R(B). The least squares with B are solved inexactly, i.e. the computed solution ¯ v of Bv ≈ c is an exact solution of a perturbed least squares problem (B + ∆B)¯ v ≈ c + ∆c, ∆B ≤ τB, ∆c ≤ τc, τκ(B) ≪ 1.

16 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Iterative solution of the null-space projected system

choose x0, solve By0 ≈ f − Ax0 compute αk and p(x)

k

∈ N(BT ) xk+1 = xk + αkp(x)

k

˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ solve Bp(y)

k

≈ r(x)

k

− αkAp(x)

k

back-substitution: A: yk+1 = yk + p(y)

k ,

B: solve Byk+1 ≈ f − Axk+1, C: solve Bvk ≈ f − Axk+1 − Byk, yk+1 = yk + vk. 9 > > > > > > > > > = > > > > > > > > > ; inner iteration r(x)

k+1 = r(x) k

− αkAp(x)

k

− Bp(y)

k

9 > > > > > > > > > > > > > > > > > > > > > = > > > > > > > > > > > > > > > > > > > > > ;

• uter

iteration

17 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Accuracy in the saddle point system

f − Axk − Byk − r(x)

k ≤ O(α3)κ(B)

1 − τκ(B) (f + AXk), − BT xk ≤ O(τ)κ(B) 1 − τκ(B)BXk, Xk ≡ max{xi | i = 0, 1, . . . , k}. Back-substitution scheme α3 A: Generic update yk+1 = yk + p(y)

k

u B: Direct substitution yk+1 = B†(f − Axk+1) τ C: Corrected dir. subst. yk+1 = yk + B†(f − Axk+1 − Byk) u

• additional least

square with B

18 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Generic update: yk+1 = yk + p(y)

k

10 20 30 40 50 60 70 80 90 100 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u), τ = 10−10, τ=10−6, τ=10−2 iteration number relative residual norms ||f−Axk−Byk||/||f−Ax0−By0||, ||r(x)

k ||/||r(x) 0 ||

10 20 30 40 50 60 70 80 90 100 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number residual norm ||−BTxk||

19 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Direct substitution: yk+1 = B†(f − Axk+1)

10 20 30 40 50 60 70 80 90 100 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number relative residual norms ||f−Axk−Byk||/||f−Ax0−By0||, ||r(x)

k ||/||r(x) 0 ||

10 20 30 40 50 60 70 80 90 100 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number residual norm ||−BTxk||

20 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Corrected direct substitution: yk+1 = yk + B†(f − Axk+1 − Byk)

10 20 30 40 50 60 70 80 90 100 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u), τ = 10−10, τ=10−6, τ=10−2 iteration number relative residual norms ||f−Axk−Byk||/||f−Ax0−By0||, ||r(x)

k ||/||r(x) 0 ||

10 20 30 40 50 60 70 80 90 100 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 τ = O(u) τ = 10−2 τ = 10−6 τ = 10−10 iteration number residual norm ||−BTxk||

21 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Conclusions

All bounds of the limiting accuracy depend on the maximum norm of computed iterates, cf. [Greenbaum, 1997]. The accuracy measured by the residuals of the saddle point problem depends on the choice of the back-substitution scheme [J, R, 2008]. Care must be taken when solving nonsymmetric systems [J, R, 2007].

50 100 150 200 250 300 350 400 450 10

−15

10

−10

10

−5

10 10

5

iteration number k relative residual norms ||BTA−1f−BTA−1Byk||/||BTA−1f|| and ||r(y)

k ||/||BTA−1f||

The residuals in the outer iteration process and the forward errors of computed approximations are proportional to the backward error in solution of inner systems.

22 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

Thank you for your attention.

http://www.cs.cas.cz/∼miro

• P. Jir´

anek and M. Rozloˇ zn´ ık. Maximum attainable accuracy of inexact saddle point solvers. SIAM J. Matrix Anal. Appl., 29(4):1297–1321, 2008.

• P. Jir´

anek and M. Rozloˇ zn´ ık. Limiting accuracy of segregated solution methods for nonsymmetric saddle point problems. J. Comput. Appl. Math., to appear.

23 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy

References

• M. Benzi, G. H. Golub, and J. Liesen. Numerical solution of saddle point
• problems. Acta Numer., 14:1–137, 2005.
• A. Greenbaum. Estimating the attainable accuracy of recursively computed

residual methods. SIAM J. Matrix Anal. Appl., 18(3):535–551, 1997.

• P. Jir´

anek and M. Rozloˇ zn´ ık. Limiting accuracy of segregated solution methods for nonsymmetric saddle point problems. J. Comput. Appl. Math., 2007. to appear.

• P. Jir´

anek and M. Rozloˇ zn´ ık. Maximum attainable accuracy of inexact saddle point solvers. SIAM J. Matrix Anal. Appl., 29(4):1297–1321, 2008.

24 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Inexact saddle point solvers and their limiting accuracy