Matching moments and matrix computations 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 SC2011 in honor of Claude Brezinski and Sebastiano Seatzu Saardinia, October 2011
CG ≡ matrix form of the Gauss-Christoffel Q. � ξ ( λ ) − 1 dω ( λ ) Ax = b , x 0 ← → ω ( λ ) , ζ ↑ ↑ n � � − 1 � ω ( n ) θ ( n ) ω ( n ) ( λ ) , T n y n = � r 0 � e 1 ← → i i i =1 x n = x 0 + W n y n ω ( n ) ( λ ) − → ω ( λ ) Z. Strakoš 2
ω ( λ ) 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) Z. Strakoš 3
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 With x 0 = 0 , n − 1 � b ∗ A − 1 b γ j � r j � 2 n A − 1 r n . + r ∗ = j =0 CG : model reduction matching 2 n moments; Golub, Meurant, Reichel, Boley, Gutknecht, Saylor, Smolarski, ......... , Meurant and S (2006), Golub and Meurant (2010), S and Tichý (2011) Z. Strakoš 4
Outline 1. CG convergence bounds based on Chebyshev polynomials 2. Sensitivity of the Gauss-Christoffel quadrature 3. PDE discretizations and matrix computations 5 Z. Strakoš
1 Beauty of Chebyshev polynomials ● Flanders and Shortley, Numerical determination of fundamental modes (1950) ● Lanczos, Chebyshev polynomials in the solution of large scale linear systems (1953) ● Stiefel, Kernel polynomials in linear algebra and their numerical applications (1958) ● Rutishauser, Theory of gradient methods (1959) For the state of the art demonstration of the beauty of Chebyshev polynomials we refer to the work of Nick Trefethen and his collaborators. Z. Strakoš 6
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 7 Z. Strakoš
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 represents nothing but the bound for the Chebyshev method, Liesen and S (2012?) Z. Strakoš 8
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š 9
q s ( λ ) 1 The polynomial has desired roots 5 4 3 2 1 0 −1 −2 −3 T 4 (x) −4 T 5 (x) −5 0 10 20 30 40 50 60 70 80 90 The Chebyshev polynomials T 4 ( λ ) , T 5 ( λ ) , and the polynomial q 1 ( λ ) , q 1 (0) = 1 having the single root at the large outlying eigenvalue. Z. Strakoš 10
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š 11
1 Mathematical model of FP CG CG in finite precision arithmetic can be seen as the exact arithmetic CG for the problem with the slightly modified distribution function with larger support, i.e., with single eigenvalues replaced by tight clusters. Paige (1971-80), Greenbaum (1989), Parlett (1990), S (1991), Greenbaum and S (1992), Notay (1993), ... , Druskin, Knizhnermann, Zemke, Wülling, Meurant, ... Recent review and update in Meurant and S, Acta Numerica (2006). Fundamental consequence: In FP computations, the composite convergence bounds eliminating in exact arithmetic large outlying eigenvalues at the cost of one iteration per eigenvalue do not, in general, work. Z. Strakoš 12
1 Axelsson (1976), quote Jennings (1977) p. 72: ... it may be inferred that rounding errors ... affects the convergence rate when large outlying eigenvalues are present. Z. Strakoš 13
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š 14
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š 15
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š 16
2 : A point going 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š 17
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š 18
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š 19
3 Local discretization and global computation Discrete (piecewise linear) solution N � u h = ζ j φ j ( η 1 , η 2 ) . j =1 ● If ζ j is known exactly, then u ( n ) = u h , and the global information is h approximated as the linear combination of the local basis functions. ● Apart from trivial cases, ζ j , which supplies the global information, is not known exactly. Z. Strakoš 20
3 Energy norm of the error Theorem Up to a small inaccuracy proportional to machine precision, �∇ ( u − u h ) � 2 + �∇ ( u h − u ( n ) �∇ ( u − u ( n ) h ) � 2 h ) � 2 = �∇ ( u − u h ) � 2 + � x − x n � 2 = A . Using zero Dirichlet boundary conditions, �∇ ( u − u h ) � 2 = �∇ u � 2 − �∇ u h � 2 . Z. Strakoš 21
Recommend
More recommend
