A Method for Estimating a Distribution of Eigenvalues using the AS - - PowerPoint PPT Presentation

A Method for Estimating a Distribution of Eigenvalues using the AS Method

Kenta SENZAKI 1) Hiroto TADANO 2) Tetsuya SAKURAI 1) Zhaojun BAI 3)

1) Graduate School of Systems and Information Engineering, University of Tsukuba 2) Graduate School of Informatics, Kyoto University 3) Department of Computer Science, University of California, Davis

RANMEP2008 Jan. 8, 2008 at NCTS −1−

Table of Contents

Motivation and Background

Earlier studies Technical issues What's new in this study

Combination of AS and CIRR

Flow of AS + CIRR Estimating a distribution of eigenvalue Speed up estimating a distribution

Numerical examples

Model of 8 DNA base pairs Epidermal Growth Factor Receptor (EGFR)

Summary and future works

RANMEP2008 Jan. 8, 2008 at NCTS −2−

Table of Contents

Motivation and Background

Earlier studies Technical issues What's new in this study

Combination of AS and CIRR

Flow of AS + CIRR Estimating a distribution of eigenvalue Speed up estimating a distribution

Numerical examples

Model of 8 DNA base pairs Epidermal Growth Factor Receptor (EGFR)

Summary and future works

RANMEP2008 Jan. 8, 2008 at NCTS −3−

Earlier studies

We consider to solve the large-scale eigenvalue problem Ax=λBx, where A, B ∈ Rn×n are symmetric, and B is positive definite. These eigenvalue problems appear in many applications such as

Electronic structure calculation Vibration analysis Molecular orbital (MO) computation …etc.

Our target problem

RANMEP2008 Jan. 8, 2008 at NCTS −4−

Earlier studies

Problem description

The matrix size n = 2K ~ 200K. The number of nonzero elements nnz = 100K ~ 400M.

• Sparse matrices, however, relatively dense.

About 100 interior eigenpairs (λj, xj) which is related to chemical reaction are needed.

• MOs around Highest Occupied MO (HOMO) and Lowest

Unoccupied MO (LUMO) are important.

Which solver should we use?

RANMEP2008 Jan. 8, 2008 at NCTS −5−

Earlier studies

Dence type: LAPACK, ScaLAPACK

Impractical in terms of computational time and memory.

Sparse type:

Iterative methods

• Shift-invert Lanczos method
• Jacobi-Davidson method … etc.

Projection methods

• Algebraic Sub-structuring method
• Contour Integral Rayleigh-Ritz method …etc.

RANMEP2008 Jan. 8, 2008 at NCTS −6−

Technical issues in iterative method

Shift-invert Lanczos (SIL) method

We need tens of eigenvalues, hence we may decompose the matrix A − σB as A − σB = LLT to reuse the matrix L in the iterative loop. However, in the near future, it will become unrealistic to decompose the matrix A − σB because the matrix appears in MO computation tends to have many nonzero elements and a complicated sparsity pattern. Furthermore, it is difficult to set the appropriate shift value σ without the knowledge of eigenvalue distribution.

RANMEP2008 Jan. 8, 2008 at NCTS −7−

Technical issues in projection method

Features

AS method uses the structure of matrix which is reordered and partitioned by Nested Dissection (ND) algorithm. AS can compute hundreds of approximate eigenpairs fast.

Issues

The accuracy of AS tends to be lower than SIL. If there is a large substructure after ND ordering, it takes long time to complete AS.

Algebraic Sub-structuring (AS) method

RANMEP2008 Jan. 8, 2008 at NCTS −8−

Technical issues in projection method

Features

CIRR method uses contour integration to construct the Rayleigh- Ritz subspace. CIRR is designed to be better adapted to parallel computational environment.

Issue

If the appropriate circular domain is put, the solutions are

• accurate. However, it is difficult to put such circle without the

knowledge of eigenvalue distribution.

Contour Integral Rayleigh-Ritz (CIRR) method

RANMEP2008 Jan. 8, 2008 at NCTS −9−

Technical issues

Feature of solvers for MO computation

△ ○ ○ Time ◎ ○ △ Parallel △ △ △ Accuracy

(Bad shift)

○ △ ○ Accuracy

(Good shift)

No No Yes Complete decomposition CIRR AS SIL △ ○ ○ Time ◎ ○ △ Parallel △ △ △ Accuracy

(Bad shift)

○ △ ○ Accuracy

(Good shift)

No No Yes Complete decomposition CIRR AS SIL

RANMEP2008 Jan. 8, 2008 at NCTS −10−

What’s new in this work

We use the AS method to obtain the rough approximate eigenvalues. From the result of AS, we estimate the eigenvalue distribution and set the circular domains for CIRR method. The accurate solution is computed by the CIRR method. Fully parallelized solver for generalized eigenvale problems

Both AS and CIRR can be implemented to be better adapted to parallel environment.

Goal

RANMEP2008 Jan. 8, 2008 at NCTS −11−

Table of Contents

Motivation and Background

Earlier studies Technical issues What's new in this study

Combination of AS and CIRR

Flow of AS + CIRR Estimating a distribution of eigenvalue Speed up estimating a distribution

Numerical examples

Model of 8 DNA base pairs Epidermal Growth Factor Receptor (EGFR)

Summary and future works

RANMEP2008 Jan. 8, 2008 at NCTS −12−

Flow of AS + CIRR

AS method Estimating a distribution of eigenvalue CIRR method A, B, σ (λj, xj)

approximate eigenvalues θ1, θ2, …, θN circular domain

RANMEP2008 Jan. 8, 2008 at NCTS −13−

Estimating a distribution of eigenvalue

Let θ1, θ2, …, θN be the approximate eigenvalues calculated by AS, which are included in the interval [t1, t2]. Apply the Gauss function where w is some weight. G(t) implies the distribution of eigenvalues, dense or sparse, at the point t. | : eigenvalue | : G(t)

RANMEP2008 Jan. 8, 2008 at NCTS −14−

How to put the circle

By the CIRR feature, we set the circular domain in the following two strategies.

Large circles are put at the sparse area of eigenvalues, and small circles are put at the dense area. The bound of circle is set at the sparse point of eigenvalue.

RANMEP2008 Jan. 8, 2008 at NCTS −15−

Speed up for estimating a distribution

If the target problem has the many nonzero elements, it tends to take a long time to complete AS.

In this case, the quite large substructures tend to appear after partitioning and reordering.

We use AS to obtain the rough eigenvalue distribution; we don’t need highly accurate solution from AS. We apply “Cutoff” to the target problem in order to reduce the number of nonzero elements, and to execute the AS method

• faster. We define Cutoff as follows with small positive value δ.

RANMEP2008 Jan. 8, 2008 at NCTS −16−

Speed up for estimating a distribution

Perturbation of eigenvalue by Cutoff value δ.

Let Aδ, Bδ be perturbation matrices which satisfy Aδ = A − Ac, Bδ = B − Bc and The perturbation error is estimated by following equation. In practice, the perturbation error by Cutoff tends not to exceed the Cutoff value δ.

Example

Model of 8 DNA base pairs

• 1,980 × 1,980 symmetric matrices.

Compute all eigenvalues in each Cutoff value by LAPACK (DSYGV)

† Kravanja, P., Sakurai, T., Sugiura, H., Van Barel, M., A Perturbation Result for Generalized Eigenvalue Problems and its Application to Error Estimation in a Quadrature Method for Computing Zeros of Analytic Functions, J. Comput. Appl. Math. 161(2):339-347, 2003.

RANMEP2008 Jan. 8, 2008 at NCTS −17−

Speed up for estimating a distribution

Numerical example of the perturbation

1e-7 1e-7 1e-6 1e-6 1e-5 1e-5 1e-4 1e-4 1e-3 1e-3 1e-2 1e-2 1e-1 1e-1 exact exact

10 −10 −20 −30 −40 −50 −60 −70 −80

1e-7 1e-7 1e-6 1e-6 1e-5 1e-5 1e-4 1e-4 1e-3 1e-3 1e-2 1e-2 1e-1 1e-1 exact exact

0.0 0.1 −0.1 −0.2 −0.3 −0.4 −0.5 0.4 0.3 0.2 0.5

RANMEP2008 Jan. 8, 2008 at NCTS −18−

Speed up for estimating a distribution

Numerical example of the perturbation

The number of eigenvalues which satisfy |λj − λj’| < δ. The number of eigenvalues which satisfy |λj − λj’| < δ under HOMO (λ1320).

72% 84% 93% 90% 91% 94% 95% 1441 1665 1851 1784 1800 1862 1886 1e-1 1e-2 1e-3 1e-4 1e-5 1e-6 1e-7 δ 100% 100% 100% 99% 100% 100% 100% 1320 1320 1320 1318 1320 1320 1320 1e-1 1e-2 1e-3 1e-4 1e-5 1e-6 1e-7 δ

RANMEP2008 Jan. 8, 2008 at NCTS −19−

Flow of AS with Cutoff + CIRR

Ac, Bc, σ

approximate eigenvalues θ1, θ2, …, θN A, B, circular domain

AS method Estimating a distribution of eigenvalue CIRR method (λj, xj) Cutoff A, B, σ

RANMEP2008 Jan. 8, 2008 at NCTS −20−

Table of Contents

Motivation and Background

Earlier studies Technical issues What's new in this study

Combination of AS and CIRR

Flow of AS + CIRR Estimating a distribution of eigenvalue Speed up estimating a distribution

Numerical examples

Model of 8 DNA base pairs Epidermal Growth Factor Receptor (EGFR)

Summary and future works

RANMEP2008 Jan. 8, 2008 at NCTS −21−

Numerical examples

 We applied the AS with Cutoff + CIRR to the generalized eigenvalue problem derived from the molecular orbital (MO) computation.

 Example 1

• Model of 8 DNA base pairs

 Example 2

• Epidermal Growth Factor Receptor (EGFR)

 We computed the eigenpairs around Highest Occupied MO (HOMO) and Lowest Unoccupied MO (LUMO). In the area of the chemistry, these MOs seem to be important to analyze the chemical reactions.

RANMEP2008 Jan. 8, 2008 at NCTS −22−

Numerical examples

Experimental environment

AS

• CPU : Intel Xeon 3.2GHz
• Memory : 2GByte
• Compiler : GNU C 3.4.3, GNU Fortran 3.4.3
• Package : ATLAS, LAPACK, MeTiS

CIRR

• CPU : Intel Core2Duo 2.2GHz
• Memory : 2GByte
• Compiler : icc 10.1, ifort 10.1
• Package : MKL

RANMEP2008 Jan. 8, 2008 at NCTS −23−

Numerical examples

Experimental condition

AS

• Projection size : 20% (example 1), 10% (example 2)
• ND level : 3 (example 1), 6 (example 2)

CIRR

• #Nodes : 24
• Block size : 16
• Solver : COCG
• Preconditioner : Sparse Direct Solver with Cutoff†

†M. Okada, T. Sakurai, K. Teranishi, A preconditioning using sparse direct solvers for approximate coefficient matrices,

• Trans. Japan SIAM Vol. 17, 2007, pp. 319−329 (in Japanese).

RANMEP2008 Jan. 8, 2008 at NCTS −24−

Numerical examples

Model of 8 DNA base pairs

Matrices properties:

• 1,980 × 1,980 symmetric matrices.
• 728,080 entries (18.6%) in matrix A.
• 592,458 entries (15.1%) in matrix B.
• 728,080 entries (18.6%) in matrix C = |A| + |B|.

Natural order After ND ordering Large Large Small Small

The sparsity pattern of the matrix C.

RANMEP2008 Jan. 8, 2008 at NCTS −25−

Numerical examples

The relation between computational time of AS with several Cutoff values and the number of nonzero elements.

1 105 2 105 3 105 4 105 5 105 6 105 7 105 8 105

0.0 0.5 1.0 1.5 2.0 2.5 10-7 10-6 10-5 10-4 10-3 10-2 nnz time #nonzero elements

Computational time [s]

Cutoff

2.2 [s] 0.8 [s]

RANMEP2008 Jan. 8, 2008 at NCTS −26−

Numerical examples

| : Eigenvalue (LAPACK), | : Eigenvalue (AS with Cutoff ).

1e-7 1e-7 1e-6 1e-6 1e-5 1e-5 1e-4 1e-4 1e-3 1e-3 1e-2 1e-2 − −0.1 0.1 0.0 0.0 0.2 0.2 0.1 0.1 0.3 0.3 0.4 0.4 0.5 0.5 − −0.5 0.5 − −0.4 0.4 − −0.3 0.3 − −0.2 0.2 LAPACK − −0.36 0.36

AS(1e-4)

LAPACK − −0.22 0.22 − −0.20 0.20 − −0.26 0.26 − −0.24 0.24 − −0.18 0.18 − −0.16 0.16 − −0.28 0.28 − −0.32 0.32 − −0.30 0.30 − −0.34 0.34

RANMEP2008 Jan. 8, 2008 at NCTS −27−

Numerical examples

The square ■ shows the relative error |λj − λj’|/|λj| and the bullet

• shows the residual norm ||Axj’−λj’Bxj’||2.

− −0.22 0.22 − −0.20 0.20 − −0.26 0.26 − −0.24 0.24 − −0.18 0.18 − −0.16 0.16 − −0.28 0.28 − −0.32 0.32 − −0.30 0.30 − −0.34 0.34 − −0.36 0.36 10−17 10−16 10−15 10−14 10−13 10−12 10−11

22 8 24 10 29 21 24 31 15 23 18 29 22 9 30 6 27 19 30 20 24 22 CIRR AS(1e-4)

RANMEP2008 Jan. 8, 2008 at NCTS −28−

Numerical examples

EGFR

Matrices properties:

• 26,461 × 26,461 symmetric matrices.
• 14,175,935 entries (2.0%) in matrix C.

Compute the eigenpairs around HOMO-LUMO.

Natural order After ND ordering Large Large Small Small

The sparsity pattern of the matrix C.

RANMEP2008 Jan. 8, 2008 at NCTS −29−

Numerical examples

The relation between computational time of AS with several Cutoff values and the number of nonzero elements.

5 106 1 107 1.5 107 200 400 600 800 1000 1200 10-7 10-6 10-5 10-4 10-3 10-2 10-1 nnz time #nonzero elements Computational time Cutoff 1074.5[s] 109.3[s]

RANMEP2008 Jan. 8, 2008 at NCTS −30−

Numerical examples

Eigenvalue distribution with several Cutoff values

1e-7 1e-7 1e-6 1e-6 1e-5 1e-5 1e-4 1e-4 1e-3 1e-3 1e-2 1e-2 1e-1 1e-1

− −0.1 0.1 0.0 0.0 0.2 0.2 0.1 0.1 0.3 0.3 0.4 0.4 0.5 0.5 − −0.5 0.5 − −0.4 0.4 − −0.3 0.3 − −0.2 0.2

• 0.08
• 0.08
• 0.04
• 0.04

0.04 0.04 0.00 0.00 0.08 0.08 0.12 0.12 0.16 0.16

• 0.12
• 0.12

AS(1e-4) AS(1e-4)

RANMEP2008 Jan. 8, 2008 at NCTS −31−

Numerical examples

The bullet ● shows the residual norm ||Axj’−λj’Bxj’||2 by CIRR. Residuals Residuals

AS(1e-4) AS(1e-4)

eigenvalue residual norm ||Axj’−λj’Bxj’||2

10-16 10-15 10-14 10-13 10-12

• 0.12
• 0.08
• 0.04

0.00 0.04 0.08 0.12 0.16

17 33 26 10 17 21 30 26 10 28

RANMEP2008 Jan. 8, 2008 at NCTS −32−

Table of Contents

Motivation and Background

Earlier studies Technical issues What's new in this study

Combination of AS and CIRR

Flow of AS + CIRR Estimating a distribution of eigenvalue Speed up estimating a distribution

Numerical examples

Model of 8 DNA base pairs Epidermal Growth Factor Receptor (EGFR)

Summary and future works

RANMEP2008 Jan. 8, 2008 at NCTS −33−

Summary and future works

Summary

We proposed the combination of AS and CIRR. It became faster to profile the eigenvalue distribution by combining Cutoff.

Future works

Apply this method to other applications. Speed up and parallelize this method. Compare this method with conventional methods.

Thank you for attention!