Challenges in Analysis of Algebraic Iterative Solvers Jörg Liesen Technical University of Berlin and Zdenˇ ek Strakoš Charles University in Prague and Czech Academy of Sciences http://www.karlin.mff.cuni.cz/˜strakos Workshop in honor of K. R. Rajagopal Prague, March 2012
Cornelius Lanczos, March 9, 1947 “The reason why I am strongly drawn to such approximation mathematics problems is not the practical applicability of the solution, but rather the fact that a very “economical” solution is possible only when it is very “adequate”. To obtain a solution in very few steps means nearly always that one has found a way that does justice to the inner nature of the problem.” Z. Strakoš 2
Albert Einstein, March 18, 1947 “Your remark on the importance of adapted approximation methods makes very good sense to me, and I am convinced that this is a fruitful mathematical aspect, and not just a utilitarian one.” Z. Strakoš 3
Algebraic iterative computations ● In iterative methods applied to linear algebraic problems, computational cost of finding sufficiently accurate approximation to the exact solution heavily depends on the particular data, i.e., ❋ on the underlying real world problem, ❋ on the mathematical model, ❋ on its discretisation. ● Any evaluation of cost in iterative computations must take into account effects of rounding errors. ● In mathematical modeling of real world phenomena, the accuracy of the computed approximation must be related to the underlying phenomena. Its evaluation can not be restricted to algebra. Z. Strakoš 4
Is there any algebraic error worth consideration? Knupp and Salari, 2003: “There may be incomplete iterative convergence (IICE) or round-off-error that is polluting the results. If the code uses an iterative solver, then one must be sure that the iterative stopping criteria is sufficiently tight so that the numerical and discrete solutions are close to one another. Usually in order-verification tests, one sets the iterative stopping criterion to just above the level of machine precision to circumvent this possibility.” Why do we care? Is not all these algebraic stuff linear and simple? 5 Z. Strakoš
HPD, x 0 , r 0 , p 0 = r 0 Conjugate Gradients: A � x − x n � A = u ∈ x 0 + K n ( A,r 0 ) � x − u � A min K n ( A, r 0 ) ≡ span { r 0 , Ar 0 , · · · , A n − 1 r 0 } γ n − 1 = ( r n − 1 , r n − 1 ) / ( p n − 1 , Ap n − 1 ) x n = x n − 1 + γ n − 1 p n − 1 r n = r n − 1 − γ n − 1 Ap n − 1 δ n = ( r n , r n ) / ( r n − 1 , r n − 1 ) p n = r n + δ n p n − 1 . Hestenes and Stiefel (1952), Lanczos (1950, 1952) This algebraic stuff is nothing but linear! Z. Strakoš 6
CG is the Gauss-Christoffel Quadrature � ξ Ax = b , x 0 ← → ω ( λ ) , f ( λ ) dω ( λ ) ζ ↑ ↑ n � � � ω ( n ) θ ( n ) ω ( n ) ( λ ) , T n y n = � r 0 � e 1 ← → f i i i =1 x n = x 0 + W n y n ω ( n ) ( λ ) − → ω ( λ ) 7 Z. Strakoš
ω ( λ ) Distribution function ω i = | ( s i , w 1 ) | 2 , w 1 = r 0 / � r 0 � λ i , s i A , are the eigenpairs of ω N 1 . . ω 4 . ω 3 ω 2 ω 1 0 . . . . . . L λ 1 λ 2 λ 3 λ N U Hestenes and Stiefel (1952), Lanczos (1952, almost unknown) Z. Strakoš 8
CG does model reduction matching 2 n moments � U n � � − 1 � λ − 1 dω ( λ ) ω ( n ) θ ( n ) = + R n ( f ) i i L i =1 � x − x 0 � 2 n -th Gauss quadrature + � x − x n � 2 A A = � r 0 � 2 � r 0 � 2 n − 1 � b ∗ A − 1 b γ j � r j � 2 n A − 1 r n . + r ∗ With x 0 = 0 , = j =0 Golub, Meurant, Reichel, Boley, Gutknecht, Saylor, Smolarski, ......... , Meurant and S (2006), Golub and Meurant (2010), S and Tichý (2011), Liesen, S, Krylov subspace methods, OUP (2012) Z. Strakoš 9
Outline 1. CG convergence bounds based on Chebyshev polynomials 2. Sensitivity of the Gauss-Christoffel quadrature 3. PDE discretizations and matrix computations Z. Strakoš 10
1 Linear bounds for the nonlinear method? � A 1 / 2 p ( A )( x − x 0 ) � � x − x n � A = min p (0)=1 deg( p ) ≤ n � Y p (Λ) Y ∗ A 1 / 2 ( x − x 0 ) � = min p (0)=1 deg( p ) ≤ n � � ≤ min 1 ≤ j ≤ N | p ( λ j ) | max � x − x 0 � A p (0)=1 deg( p ) ≤ n Using the shifted Chebyshev polynomials on the interval [ λ 1 , λ N ] , �� � n κ ( A ) − 1 � x − x n � A ≤ 2 � x − x 0 � A . � κ ( A ) + 1 Z. Strakoš 11
1 Minimization property and the bound This bound has a remarkably wiggling history: ● Markov (1890) ● Flanders and Shortley (1950) ● Lanczos (1953), Kincaid (1947), Young (1954, ... ) ● Stiefel (1958), Rutishauser (1959) ● Meinardus (1963), Kaniel (1966) ● Daniel (1967a, 1967b) ● Luenberger (1969) It is relevant to the Chebyshev method! Z. Strakoš 12
1 Composite bounds considering large outliers? This bound should not be used in connection with the behaviour of CG unless κ ( A ) = λ N /λ 1 is really small or unless the (very special) distribution of eigenvalues makes it relevant. In particular, one should be very careful while using it as a part of a composite bound in the presence of the large outlying eigenvalues � � T n − s ( λ j ) | � � min 1 ≤ j ≤ N | q s ( λ j ) p ( λ j ) | max ≤ 1 ≤ j ≤ N | q s ( λ j ) | max � � T n − s (0) � � p (0)=1 deg( p ) ≤ n − s � � T n − s ( λ j ) � � < max � . � � T n − s (0) � 1 ≤ j ≤ N − s This Chebyshev method bound on the interval [ λ 1 , λ N − s ] is then valid after s initial steps. Z. Strakoš 13
1 Quote (2009, ... ): the desired accuracy ǫ Theorem. After � � � ln(2 /ǫ ) λ N − s k = s + 2 λ 1 iteration steps the CG will produce the approximate solution x n satisfying � x − x n � A ≤ ǫ � x − x 0 � A . This recently republished and used statement is in finite precision arithmetic not true at all. Z. Strakoš 14
1 Axelsson (1976), Jennings (1977) p. 72: ... it may be inferred that rounding errors ... affects the convergence rate when large outlying eigenvalues are present. Z. Strakoš 15
1 The composite bounds completely fail 0 0 10 10 −5 −5 10 10 −10 −10 10 10 −15 −15 10 10 0 20 40 60 80 100 0 20 40 60 80 100 Composite bounds with varying number of outliers: Exact CG (left) and FP CG (right), Gergelits (2011). Z. Strakoš 16
2 CG and Gauss-Christoffel quadrature errors � U n � � − 1 � λ − 1 dω ( λ ) ω ( n ) θ ( n ) = + R n ( f ) i i L i =1 � x − x 0 � 2 n -th Gauss quadrature + � x − x n � 2 A A = � r 0 � 2 � r 0 � 2 Consider two slightly different distribution functions with � U λ − 1 dω ( λ ) ≈ I n I ω = ω L � U λ − 1 d ˜ ω ( λ ) ≈ I n I ˜ ω = ω ˜ L Z. Strakoš 17
2 Sensitivity of the Gauss-Christoffel Q. 0 10 −5 10 quadrature error − perturbed integral quadrature error − original integral −10 10 0 5 10 15 20 iteration n difference − estimates 0 10 difference − integrals −5 10 −10 10 0 5 10 15 20 iteration n Z. Strakoš 18
2 The point goes back to 1814 1. Gauss-Christoffel quadrature for a small number of quadrature nodes can be highly sensitive to small changes in the distribution function that enlarge its support. In particular, the difference between the corresponding quadrature approximations (using the same number of quadrature nodes) can be many orders of magnitude larger than the difference between the integrals being approximated. 2. This sensitivity in Gauss-Christoffel quadrature can be observed for discontinuous, continuous, and even analytic distribution functions, and for analytic integrands uncorrelated with changes in the distribution functions, with no singularity close to the interval of integration. Z. Strakoš 19
2 Theorem - O’Leary, S, Tichý (2007) Consider distribution functions ω ( λ ) and ω ( λ ) ˜ on [ L, U ] . Let p n ( λ ) = ( λ − ˜ λ 1 ) . . . ( λ − ˜ p n ( λ ) = ( λ − λ 1 ) . . . ( λ − λ n ) and ˜ λ n ) be the n th orthogonal polynomials corresponding to ω and ω ˜ respectively, with p s ( λ ) = ( λ − ξ 1 ) . . . ( λ − ξ s ) ˆ their least common multiple. ∆ n f ′′ If is continuous on [ L, U ] , then the difference between ω, ˜ ω I n I n I ω I ˜ the approximation to and the approximation to ω , ω ω ˜ obtained from the n -point Gauss-Christoffel quadrature, is bounded as � � � U � U � � | ∆ n � � ω | ≤ p s ( λ ) f [ ξ 1 , . . . , ξ s , λ ] dω ( λ ) − ˆ p s ( λ ) f [ ξ 1 , . . . , ξ s , λ ] d ˜ ˆ ω ( λ ) ω, ˜ � � � � L L � � � U � U � � � � + f ( λ ) dω ( λ ) − f ( λ ) d ˜ ω ( λ ) � . � � � L L Z. Strakoš 20
3 Take very simple model boundary value problem − ∆ u = 16 η 1 η 2 (1 − η 1 )(1 − η 2 ) on the unit square with zero Dirichlet boundary conditions. Galerkin finite element method (FEM) discretization with linear basis functions on the regular triangular grid with the mesh size h = 1 / ( m + 1) , where m is the number of inner nodes in each direction. Discrete (piecewise linear) solution N � u h = ζ j φ j ( η 1 , η 2 ) . j =1 Computational error u − u ( n ) u h − u ( n ) = u − u h + . h h � �� � � �� � � �� � discretisation error total error algebraic error Z. Strakoš 21
Recommend
More recommend
