Preconditioned Locally Harmonic Residual Methods for Interior - PowerPoint PPT Presentation

Preconditioned Locally Harmonic Residual Methods for Interior Eigenvalue Computations Eugene Vecharynski Lawrence Berkeley National Laboratory Computational Research Division 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.7 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.34 0.79 0.34 0.79 0.00 0.00 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.62 0.27 0.00 0.23 0.12 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.62 0.27 0.00 0.23 0.12 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.65 0.29 0.15 0.00 0.39 0.24 0.17 0.11 0.04 0.00 0.00 0.00 0.00 0.00 0.29 0.15 0.00 0.39 0.24 0.17 0.11 0.04 0.00 0.00 0.00 0.00 0.00 0.09 0.06 0.00 0,44 0.81 0.63 0.49 0.19 0.14 0.02 0.00 0.00 0.00 0.09 0.06 0.00 0,44 0.81 0.63 0.49 0.19 0.14 0.02 0.00 0.00 0.00 0.01 0.00 0.00 0.21 0.00 0.12 0.35 0.41 0.23 0.15 0.01 0.02 0.00 0.01 0.00 0.00 0.21 0.00 0.12 0.35 0.41 0.23 0.15 0.01 0.02 0.00 0.6 0.00 0.01 0.00 0.03 0.02 0.04 0.11 0.22 0.51 0.39 0.86 0.91 0.13 0.00 0.01 0.00 0.03 0.02 0.04 0.11 0.22 0.51 0.39 0.86 0.91 0.13 0.00 0.00 0.00 0.00 0.00 0.01 0.05 0.09 0.31 0.61 0.09 0.01 0.89 0.00 0.00 0.00 0.00 0.00 0.01 0.05 0.09 0.31 0.61 0.09 0.01 0.89 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.08 0.01 0.00 0.00 0.03 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.08 0.01 0.00 0.00 0.03 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.55 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Im 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.5 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.45 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.4 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.34 0.34 0.79 0.79 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.62 0.27 0.00 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.62 0.27 0.00 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.35 0.29 0.00 0.92 0.29 0.00 0.92 0.00 0.17 0.11 0.04 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.17 0.11 0.04 0.00 0.00 0.00 0.00 0.00 0.00 0.09 0.00 0.09 0.00 0.00 0.00 0.00 0.63 0.49 0.19 0.14 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.63 0.49 0.19 0.14 0.02 0.00 0.00 0.00 0.01 0.00 0.01 0.00 0.00 0.00 0.00 0.12 0.35 0.41 0.23 0.15 0.01 0.02 0.00 0.00 0.00 0.00 0.12 0.35 0.41 0.23 0.15 0.01 0.02 0.00 0.3 0.00 0.01 0.00 0.00 0.00 0.04 0.11 0.22 0.51 0.39 0.86 0.91 0.13 0.00 0.01 0.00 0.00 0.00 0.04 0.11 0.22 0.51 0.39 0.86 0.91 0.13 − 0.4 − 0.3 − 0.2 − 0.1 0 0.00 0.00 0.00 0.00 0.00 0.01 0.05 0.09 0.31 0.61 0.09 0.01 0.89 0.00 0.00 0.00 0.00 0.00 0.01 0.05 0.09 0.31 0.61 0.09 0.01 0.89 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.08 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.08 0.01 0.00 0.00 0.00 Re 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 SIAM Conference on Applied Linear Algebra, October 29, 2015 Atlanta, GA

The Preconditioned Locally Harmonic Residual methods A new class of robust block preconditioned iterative eigensolvers for computing interior eigenpairs I Hermitian eigenproblems The PLHR algorithm (EV, A. Knyazev) I Non-Hermitian eigenproblems The Generalized PLHR (GPLHR) algorithm (EV, C. Yang, F. Xue)

Problem Compute a subset of k right eigenpairs ( λ , x ) of a non-Hermitian matrix pair ( A , B ) that are closest to a given shift σ ∈ C A , B ∈ C n × n Ax = λ Bx , A , B can be extremely large 0.7 0.65 A , B may not be explicitly given 0.6 0.55 Im ( A , B ) is regular 0.5 0.45 0.4 σ can point to the spectrum’s 0.35 interior or exterior 0.3 − 0.4 − 0.3 − 0.2 − 0.1 0 Re

Many origins of the problem Study of excited states of molecules (complex scaling, EOM-CC) Chemical reactions, stability analysis in fluid dynamics, crystal growth simulations, power systems, Markov chains, plasmas, ETC.

Available eigensolvers I (Inexact) Inverse subspace iteration Su ffi ciently accurate ( A − σ B ) − 1 is needed I Arnoldi (ARPACK) Su ffi ciently accurate ( A − σ B ) − 1 is needed I Block Generalized Davidson (Morgan-92) Robustness issues, may need a larger search subspace I JDQR/JDQZ (Fokkema et. al. SISC-98) “One-by-one” eigenpair computation, lack of BLAS3

Need an eigensolver that ... I Avoids “shift-and-invert” I Takes advantage of preconditioning, T ≈ ( A − σ B ) − 1 I Performs block iteration (parallel performance, BLAS3) I Is robust and reliable if memory is limited I Handles standard and generalized, interior and exterior, eigenproblems in a uniform manner I Is simple (to the extent a non-Hermitian solver can be ;-))

An example from the “Hermitian world”: LOBPCG The LOBPCG algorithm (Knyazev SISC-01): n X ( i ) , T ( AX ( i ) − BX ( i ) Λ ( i ) ) , X ( i − 1) o X ( i +1) ← col I Block iterations I No inversions I Short-term recurrence I The Rayleigh–Ritz procedure over low-dimensional subspaces

An example from the “Hermitian world”: LOBPCG The LOBPCG algorithm (Knyazev SISC-01): n X ( i ) , T ( AX ( i ) − BX ( i ) Λ ( i ) ) , X ( i − 1) o X ( i +1) ← col I Block iterations I No inversions I Short-term recurrence I The Rayleigh–Ritz procedure over low-dimensional subspaces Does NOT extend to non-Hermitian eigenproblems

Di ffi culties 1. The eigendecomposition AX = BX Λ may not exist or X can be ill-conditioned. Partial generalized Schur form should be computed ⇢ = AV QR A , λ j = R A ( j , j ) / R B ( j , j ) = BV QR B 2. A correct definition of a block residual is not clear 3. The LOBPCG-like trial subspace is “too small” 4. The standard Rayleigh–Ritz may not be appropriate for interior eigenvalue computations

The “ Q -free” partial generalized Schur form I Partial generalized Schur form ⇢ AV = QR A , λ j = R A ( j , j ) / R B ( j , j ) BV = QR B I The “ Q -free” partial generalized Schur form AVM B = BVM A , λ j = M A ( j , j ) / M B ( j , j ) , where M A , M B are upper triangular. – Gives a numerically stable analogue of AX Λ B = BX Λ A – Allows defining a block residual AVM B − BVM A

The “ Q -free” partial generalized Schur form I Partial generalized Schur form ⇢ AV = QR A , λ j = R A ( j , j ) / R B ( j , j ) BV = QR B I The “ Q -free” partial generalized Schur form AVM B = BVM A , λ j = M A ( j , j ) / M B ( j , j ) , where M A , M B are upper triangular. – Gives a numerically stable analogue of AX Λ B = BX Λ A – Allows defining a block residual AVM B − BVM A Does the “Q-free” form exist for any regular pair ( A , B )?

The “ Q -free” partial generalized Schur form For any regular pair ( A , B ), there exist AV M B = BV M A , λ j = M A ( j , j ) / M B ( j , j ) where M A = G 2 G − 1 R A , M B = I − G 1 G − 1 R A , with an upper triangular G = R A G 1 + R B G 2 , G ( j , j ) = 1 , and G 1 , G 2 diagonal 8 0 , | R A ( j , j ) | < | R B ( j , j ) | < 1 − R B ( j , j ) G 1 ( j , j ) = , otherwise : R A ( j , j ) ⇢ 1 / R B ( j , j ) , | R A ( j , j ) | < | R B ( j , j ) | G 2 ( j , j ) = 1 , otherwise

The simplest GPLHR trial subspace Given approximations V ( i ) , R ( i ) A , and R ( i ) B , the “ Q -free” partial generalized Schur form allows I Defining a preconditioned residual W ( i ) = T ( AV ( i ) M ( i ) B − BV ( i ) M ( i ) A ) , where T is a preconditioning operator. I Constructing a LOBPCG-like trial subspace Z = col { V ( i ) , W ( i ) , P ( i ) } where P ( i ) is a block of additional search directions.

The simplest GPLHR trial subspace Given approximations V ( i ) , R ( i ) A , and R ( i ) B , the “ Q -free” partial generalized Schur form allows I Defining a preconditioned residual W ( i ) = T ( AV ( i ) M ( i ) B − BV ( i ) M ( i ) A ) , where T is a preconditioning operator. I Constructing a LOBPCG-like trial subspace Z = col { V ( i ) , W ( i ) , P ( i ) } where P ( i ) is a block of additional search directions. How to construct a larger trial subspace?