Krylov subspace methods and exascale computations: good match or lost case? Zdenˇ ek Strakoš Charles University in Prague and Czech Academy of Sciences http://www.karlin.mff.cuni.cz/˜strakos SPPEXA Symposium, Münich, January 2016
Personal prehistory Strakos, Z., Efficiency and Optimizing of Algorithms and Programs on the Host Computer / Array Processor System, Parallel Computing, 4, 1987, pp. 189-209. ● Host Computer (0.2 MFlops) / Array Processor (up to 10 MFlops). ● Large instruction overhead and slow data transfers. ● Pipelining, several arithmetic units. ● Possible overlap of data transfers and arithmetic. ● Slow scalar operations. Basic problems and principles are not even after thirty years that much different. Z. Strakoš 2
Preconditioned algebraic CG r 0 = b − Ax 0 , solve Mz 0 = r 0 , p 0 = z 0 For n = 1 , . . . , n max z ∗ n − 1 r n − 1 = α n − 1 p ∗ n − 1 Ap n − 1 = x n − 1 + α n − 1 p n − 1 , stop when the stopping criterion is satisfied x n = r n − 1 − α n − 1 Ap n − 1 r n = solve for z n Mz n r n , z ∗ n r n β n = z ∗ n − 1 r n − 1 = z n + β n p n − 1 p n End Z. Strakoš 3
Obstacles for parallelization ● Synchronized recursion. ● Matrix-vector multiplication and vector updates are linear and (possibly) fast. Preconditioning is expensive (substantial global communication). ● Scalar coefficients bring in nonlinearity and require inner products. However, for the approximation power of the methods, nonlinearity is essential! ● Parallelization can lead to numerical instabilities. Z. Strakoš 4
Parallel (communication sensitive) algorithms? ● Block recursion in order to increase arithmetic/communication ratio. ● Numerical stability is crucial. ● Stopping criteria can save the case. Size of the blocks? ● Preconditioning means an approximate solution of a part of the problem. State-of-the-art in the algorithmic developments: E. Carson, Communication-Avoiding Krylov Subspace Methods in Theory and Practice, PhD Thesis, UC at Berkeley, CA, 2015. 5 Z. Strakoš
Outline 1. Philosophy of using Krylov subspace methods 2. Nonlinear model reduction 3. Inexact Krylov? 4. Operator and algebraic preconditioning 5. Krylov subspaces and discretization 6. Stopping criteria? Z. Strakoš 6
1 Plethora of Krylov subspace methods ● Thorough analysis and fair comparison of several important methods should be given priority to overproduction of algorithmic variations. ● Krylov subspace methods are efficient providing that they “do justice to the inner nature of the problem.” (C. Lanczos, 1947). Infinite dimensional considerations are very useful. ● Oversimplification is dangerous. Widespread worst scenario analysis restricted to the operator only, universal contraction-based bounds, asymptotic considerations, unjustified or hidden restrictive assumptions. ● Results pointing out difficulties should be taken as an inspiration. They are instead unwanted and often labeled as “negative.” 7 Z. Strakoš
1 Málek and S, SIAM Spotlight, 2015 ⇒ Z. Strakoš 8
Outline 1. Philosophy of using Krylov subspace methods 2. Nonlinear model reduction 3. Inexact Krylov? 4. Operator and algebraic preconditioning 5. Krylov subspaces and discretization 6. Stopping criteria? Z. Strakoš 9
2 Operator form of the BVP and preconditioning Let be a real (infinite dimensional) Hilbert space with the inner V let V # product ( · , · ) V : V × V → R , be the dual space of bounded linear functionals on V . Consider a bounded and coercive operator A : V → V # V # and the equation in b ∈ V # . A : V → V # , A x = b , x ∈ V, Using the Riesz map, ( τ A x − τb, v ) V = 0 for all v ∈ V . The Riesz map τ can be interpreted as transformation of the original V # A x = b problem in into the equation in V τ A x = τb, τ A : V → V, x ∈ V, τb ∈ V , which is (unfortunately) called preconditioning. Z. Strakoš 10
2 Model reduction using Krylov subspaces B (= τ A ) Let be a bounded linear operator on the Hilbert space V . Choosing z 0 (= τb − τ A x 0 ) ∈ V . Consider the Krylov sequence z 0 , z 1 = B z 0 , z 2 = B z 1 = B 2 z 0 , . . . , z n = B z n − 1 = B n z n − 1 , . . . Determine a sequence of operators B n defined on the sequence of nested subspaces V n = span { z 0 , . . . , z n − 1 } , with the projector E n onto V n , such that (Vorobyev (1958, 1965)) z 1 = B z 0 = B n z 0 , z 2 = B 2 z 0 = ( B n ) 2 z 0 , . . . z n − 1 = B n − 1 z 0 = ( B n ) n − 1 z 0 , E n z n = E n B n z 0 = ( B n ) n z 0 . Z. Strakoš 11
2 Bounded self-adjoint operators in V � B x = f ← → ω ( λ ) , F ( λ ) dω ( λ ) ↑ ↑ n � � ω ( n ) θ ( n ) � ω ( n ) ( λ ) , T n y n = � f � V e 1 ← → F i i i =1 Using F ( λ ) = λ − 1 gives (assuming coercivity) � λ U n + � u − u n � 2 � − 1 � λ − 1 dω ( λ ) ω ( n ) θ ( n ) � a = i i � f � 2 λ L V i =1 Stieltjes (1894) and Vorobyev (1958) moment problems for self-adjoint bounded operators reduce to the Gauss-Christoffel quadrature (1814). No one would consider describing it by contraction. Z. Strakoš 12
2 CG in Hilbert spaces r 0 = b − A x 0 ∈ V # , p 0 = τr 0 ∈ V For n = 1 , 2 , . . . , n max α n − 1 = � r n − 1 , τr n − 1 � �A p n − 1 , p n − 1 � = ( τr n − 1 , τr n − 1 ) V ( τ A p n − 1 , p n − 1 ) V x n = x n − 1 + α n − 1 p n − 1 , stop when the stopping criterion is satisfied r n = r n − 1 − α n − 1 A p n − 1 � r n , τr n � ( τr n , τr n ) V β n = � r n − 1 , τr n − 1 � = ( τr n − 1 , τr n − 1 ) V p n = τr n + β n p n − 1 End Hayes (1954); Vorobyev (1958, 1965); Karush (1952); Stesin (1954) Superlinear convergence for (identity + compact) operators. Z. Strakoš 13
Outline 1. Philosophy of using Krylov subspace methods 2. Matching moments model reduction 3. Inexact Krylov? 4. Operator and algebraic preconditioning 5. Krylov subspaces and discretization 6. Stopping criteria? Z. Strakoš 14
3 Delay of convergence due to inexactness 0 10 0 10 −5 10 −5 10 ? −10 −10 10 10 residual smooth ubound −15 backward error 10 loss of orthogonality approximate solution −15 10 error 0 100 200 300 400 500 600 700 800 0 20 40 60 80 100 iteration number Here numerical inexactness due to roundoff. How much may we relax accuracy of the most costly operations without causing an unwanted delay and/or affecting the maximal attainable accuracy? Z. Strakoš 15
Outline 1. Philosophy of using Krylov subspace methods 2. Nonlinear model reduction 3. Inexact Krylov? 4. Operator and algebraic preconditioning 5. Krylov subspaces and discretization 6. Stopping criteria? Z. Strakoš 16
4 Restriction to finite dimensional subspace V h Φ h = ( φ ( h ) 1 , . . . , φ ( h ) N ) be a basis of the subspace V h ⊂ V , Let h = ( φ ( h )# , . . . , φ ( h )# Φ # be the canonical basis of its dual V # ) let h , 1 N ( V # Φ # h = A V h ) . Using the coordinates in Φ h and in h , � f, v � → v ∗ f , ( u, v ) V → v ∗ Mu , ( M ij ) = (( φ j , φ i ) V ) i,j =1 ,...,N , A u = A Φ h u = Φ # A u → Au , h Au ; ( A ij ) = ( a ( φ j , φ i )) i,j =1 ,...,N , τf = τ Φ # τf → M − 1 f , h f = Φ h M − 1 f ; b = Φ # h b , x n = Φ h x n , p n = Φ h p n , r n = Φ # we get with the h r n algebraic CG formulation. Z. Strakoš 17
4 Galerkin discretization gives matrix CG in V h r 0 = b − Ax 0 , solve Mz 0 = r 0 , p 0 = z 0 For n = 1 , . . . , n max z ∗ n − 1 r n − 1 α n − 1 = p ∗ n − 1 Ap n − 1 = x n − 1 + α n − 1 p n − 1 , stop when the stopping criterion is satisfied x n = r n − 1 − α n − 1 Ap n − 1 r n = solve for z n Mz n r n , z ∗ n r n = β n z ∗ n − 1 r n − 1 = z n + β n p n − 1 p n End Günnel, Herzog, Sachs (2014); Málek, S (2015) Z. Strakoš 18
4 Observations ● Unpreconditioned CG, i.e. M = I , corresponds to the discretization basis Φ orthonormal wrt ( · , · ) V . ● Orthogonalization of the discretization basis will result in the unpreconditioned algebraic CG that is applied to the preconditioned algebraic system. The resulting matrix of this preconditioned algebraic system is not sparse! ● Any algebraic preconditioning applied to the algebraic system that arise from discretization can be interpreted within this operator preconditioning framework. Z. Strakoš 19
Outline 1. Philosophy of using Krylov subspace methods 2. Nonlinear model reduction 3. Inexact Krylov? 4. Operator and algebraic preconditioning 5. Krylov subspaces and discretization 6. Stopping criteria? Z. Strakoš 20
5 Conjugate gradient method - first n steps α 1 β 2 ... ... β 2 ... ... ... T n = ... ... β n β n α n is the Jacobi matrix of the orthogonalization coefficients and the CG method is formulated by T n y n = � τr 0 � V e 1 , x n = x 0 + Q n y n , x n ∈ V n . Infinite dimensional Krylov subspace methods perform discretization via model reduction. Z. Strakoš 21
Outline 1. Philosophy of using Krylov subspace methods 2. Nonlinear model reduction 3. Inexact Krylov? 4. Operator and algebraic preconditioning 5. Krylov subspaces and discretization 6. Stopping criteria? Z. Strakoš 22
Recommend
More recommend
