Numerical behavior of inexact saddle point solvers anek 1 , 2 , - - PowerPoint PPT Presentation

Numerical behavior of inexact saddle point solvers

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

9th IMACS International Symposium

• n Iterative Methods in Scientific Computings

March 17–20, 2008, Lille, France

1 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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].

2 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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

3 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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.

4 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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

5 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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.

6 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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 ||

7 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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

8 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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)||

9 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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)||

10 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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)||

11 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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

12 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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, 2008b].

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.

13 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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., 215:28–37, 2008.

14 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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.

15 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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

16 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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

17 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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||

18 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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||

19 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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||

20 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers

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. Maximum attainable accuracy of inexact saddle point solvers. SIAM J. Matrix Anal. Appl., 29(4):1297–1321, 2008a.

• P. Jir´

anek and M. Rozloˇ zn´ ık. Limiting accuracy of segregated solution methods for nonsymmetric saddle point problems. J. Comput. Appl. Math., 215: 28–37, 2008b.

21 Pavel Jir´ anek, Miroslav Rozloˇ zn´ ık Numerical behavior of inexact saddle point solvers