(Preconditioning) Chebyshev subspace iteration applied to sequences - PowerPoint PPT Presentation

Mitglied der Helmholtz-Gemeinschaft (Preconditioning) Chebyshev subspace iteration applied to sequences of dense eigenproblems in ab initio simulations NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli

Motivation and Goals Electronic Structure Mathematical Model Band energy gap + SIMULATIONS Conductivity Algorithmic Structure Forces, etc. Extracting & Exploiting Information

Motivation and Goals Electronic Structure Mathematical Model Band energy gap + SIMULATIONS Conductivity Algorithmic Structure Forces, etc. Extracting & Exploiting Information Performance More Efficiency Physics Scalability NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 2

Outline Sequences of generalized eigenproblems in all-electron computations The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) Exploiting approximate eigenvectors: numerical results Optimizations schemes NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 3

Outline Sequences of generalized eigenproblems in all-electron computations The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) Exploiting approximate eigenvectors: numerical results Optimizations schemes NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 4

Density Functional Theory scheme Self-consistent cycle Solve a set of Initial guess Compute discretized eigenproblems for charge density Kohn-Sham P ( ℓ ) k 1 ... P ( ℓ ) equations n start ( r ) k N No Compute new OUTPUT Yes Converged? Electronic charge density | n ( ℓ ) − n ( ℓ − 1 ) | < η structure, n ( ℓ ) ( r ) ... NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 5

Density Functional Theory scheme Self-consistent cycle Solve a set of Initial guess Compute discretized eigenproblems for charge density Kohn-Sham P ( ℓ ) k 1 ... P ( ℓ ) equations n start ( r ) k N No Compute new OUTPUT Yes Converged? Electronic charge density | n ( ℓ ) − n ( ℓ − 1 ) | < η structure, n ( ℓ ) ( r ) ... Observations: 1 every P ( ℓ ) : A ( ℓ ) k x = B ( ℓ ) k λ x is a generalized eigenvalue problem; k 2 A and B are DENSE and hermitian (B is also pos. def.); 3 P k s with different k index have different size and are independent from each other. 4 k = 1 : 10 − 100 ℓ = 1 : 20 − 50 ; NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 5

Sequences of eigenproblems Current solving strategy The set of generalized eigenproblems P ( 1 ) ... P ( ℓ ) P ( ℓ + 1 ) ... P ( N ) is handled as a set of disjoint problems ( P ) N ; Each problem P ( ℓ ) is solved independently using a direct solver as a black-box from a standard library (i.e. ScaLAPACK). NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 6

Sequences of eigenproblems Current solving strategy The set of generalized eigenproblems P ( 1 ) ... P ( ℓ ) P ( ℓ + 1 ) ... P ( N ) is handled as a set of disjoint problems ( P ) N ; Each problem P ( ℓ ) is solved independently using a direct solver as a black-box from a standard library (i.e. ScaLAPACK). Extracting information − → searching for a new strategy Treated the set of eigenproblems ( P ) N as a sequence � P ( ℓ ) � ; Investigated the presence of a connection between adjacent problems, numerical simulations analyzed employing a parameter-based inverse problem method; collected data on angles b/w eigenvectors of adjacent eigenproblems; discovered evolution of eigenvectors along the sequence. NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 6

Angles evolution fixed k Example: a metallic compound at fixed k Evolution of subspace angle for eigenvectors of k − point 1 and lowest 75 eigs 0 10 AuAg Angle b/w eigenvectors of adjacent iterations − 2 10 − 4 10 − 6 10 − 8 10 − 10 10 2 6 10 14 18 22 Iterations (2 − > 22) NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 7

Correlation and its exploitation ∃ correlation between successive eigenvectors x ( ℓ − 1 ) and x ( ℓ ) ; angles are small after the first few iterations. Note: Mathematical model � Correlation. Correlation ⇐ numerical analysis of the simulation . Exploiting information: S IMULATION ⇒ A LGORITHM The stage is favorable to an iterative eigensolver where the solution of P ( ℓ − 1 ) is used to solve P ( ℓ ) . 1 Approximate eigenvectors can speed-up iterative solvers (EDN, M. Berljafa [arXiv:1206.3768]) 2 Developed of a block iterative eigensolver (ChFSI) that can maximally exploit the correlation (EDN, M. Berljafa [arXiv:1305.5120]) 3 ChFSI is competitive with direct methods for dense problems in ab initio methods (EDN M. Berljafa, and in preparation). NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 8

Algorithmic digression Direct solvers. Iterative solvers. NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 9

Algorithmic digression Direct solvers. Iterative solvers. ∗ ∗ ∗ ∗ ∗ ∗   ∗ ∗ ∗ ∗ ∗ ∗   ∗ ∗ ∗ ∗ ∗ ∗     ∗ ∗ ∗ ∗ ∗ ∗     ∗ ∗ ∗ ∗ ∗ ∗   ∗ ∗ ∗ ∗ ∗ ∗  ∗ ∗  ∗ ∗ ∗   ∗ ∗ ∗     ∗ ∗ ∗     ∗ ∗ ∗   ∗ ∗ NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 9

Algorithmic digression Direct solvers. Iterative solvers. ∗ ∗ ∗ ∗ ∗ ∗ | λ 1 | > | λ 2 | > | λ 3 | > ...   ∗ ∗ ∗ ∗ ∗ ∗   Ax j = λ j x j ∗ ∗ ∗ ∗ ∗ ∗     ∗ ∗ ∗ ∗ ∗ ∗   v = ∑ j γ j x j   ∗ ∗ ∗ ∗ ∗ ∗   ∗ ∗ ∗ ∗ ∗ ∗ � λ j � Av = ∑ j λ j γ j x j ⇒ A k v = ∑ j λ k j γ j x j = λ 1 x 1 + ∑ j ≥ 2 λ 1 x j  ∗ ∗  � � ∗ ∗ ∗ λ 1 � �   ∗ ∗ ∗ Rate of convergence → magnitude of � �     λ j � � ∗ ∗ ∗   � �   ∗ ∗ ∗   ∗ ∗ NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 9

Algorithmic digression Direct solvers. Iterative solvers. Sparse matrices. Dense matrices. NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 9

Algorithmic choice Direct solvers. Iterative solvers. Sparse matrices. Dense matrices. NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 10

Outline Sequences of generalized eigenproblems in all-electron computations The algorithm: Chebyshev Filtered Sub-space Iteration method ( ChFSI ) Exploiting approximate eigenvectors: numerical results Optimizations schemes NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 11

Selecting an iterative eigensolver Two are the main properties an iterative algorithm has to comply with: 1 The ability to receive as input a sizable set Z 0 of approximate eigenvectors; 2 The capacity to solve simultaneously for a substantial portion of eigenpairs. NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 12

Selecting an iterative eigensolver Two are the main properties an iterative algorithm has to comply with: 1 The ability to receive as input a sizable set Z 0 of approximate eigenvectors; 2 The capacity to solve simultaneously for a substantial portion of eigenpairs. ChFSI constitutes the natural choice: it accepts the full set of multiple starting vectors; it avoids stalling when facing small clusters of eigenvalues; when augmented with polynomial accelerators it has a much faster convergence rate; converged eigenvectors can be easily locked; the degree of the polynomial can be opportunely optimized. NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 12

The core of the algorithm: Chebyshev filter Chebyshev polynomials The Chebyshev polynomial C m of the first kind of order m , is defined as � cos ( m arccos ( x )) , x ∈ [ − 1 , 1 ] ; C m ( x ) = cosh ( m arccosh ( x )) , | x | > 1 . Three-terms recurrence relation C m + 1 ( x ) = 2 xC m ( x ) − C m − 1 ( x ) ; m ∈ N , C 0 ( x ) = 1 , C 1 ( x ) = x Degree 5 Degree 10 6 x 10 3 500 2 400 1 300 0 − 1 200 − 2 100 − 3 0 − 3 − 2 − 1 0 1 2 3 − 3 − 2 − 1 0 1 2 3 Degree 15 Degree 20 10 14 x 10 x 10 2.5 1.5 2 1 0.5 1.5 0 1 − 0.5 − 1 0.5 − 1.5 0 − 3 − 2 − 1 0 1 2 3 − 3 − 2 − 1 0 1 2 3 NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 13

The core of the algorithm: Chebyshev filter The basic principle Theorem Let | γ | > 1 and P m denote the set of polynomials of degree smaller or equal to m . Then the extremum t ∈ [ − 1 , 1 ] | p ( t ) | p ∈ P m , p ( γ )= 1 max min is reached by = C m ( t ) p m ( t ) . C m ( γ ) . A generic vector v is very quickly aligned in the direction of the eigenvector corresponding to the extremal eigenvalue λ 1 n n v m = p m ( A ) v ∑ ∑ = s i p m ( A ) x i = s i p m ( λ i ) x i i = 1 i = 1 C m ( λ i − c n e ) ∑ = s 1 x 1 + x i ∼ s i s 1 x 1 C m ( λ 1 − c ) i = 2 e NASCA 2013 Calais, France, June 24th M. Berljafa and E. Di Napoli Folie 14