A divide-and-conquer algorithm for a symmetric eigenproblem Binh T. - PowerPoint PPT Presentation

A divide-and-conquer algorithm for a symmetric eigenproblem A divide-and-conquer algorithm for a symmetric eigenproblem Binh T. Nguyen Anh-Duc Luong-Thanh Nguyen Dinh Thuc Bui Van Thach March 5, 2012

A divide-and-conquer algorithm for a symmetric eigenproblem Outline Introduction 1 Proposed algorithm 2 Numerical experiments 3 Conclusion 4

A divide-and-conquer algorithm for a symmetric eigenproblem Introduction Problem Statement Compute the eigensystem of the symmetric tri-block-diagonal matrix:  A 11 A 12 0 0 ... 0  0 ... 0 A 21 A 22 A 23     A = 0 A 32 A 33 A 34 ... 0      ... ... ... ... ... ...    0 0 0 ... A N , N − 1 A NN where A ij ∈ R ( 2 M + 1 ) × ( 2 M + 1 ) , A i , i + 1 has only one non-zero entry, A i + 1 , i = A T i , i + 1 , and P = ( 2 M + 1 ) N is the size of matrix A

A divide-and-conquer algorithm for a symmetric eigenproblem Introduction Objectives The proposed algorithm needs to satisfy: Accuracy Fast (faster than O ( P 3 ) where P is the size of matrix)

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - Dividing step We represent matrix A as following � � B 11 B 12 A = B 21 B 22 Using rank-one decomposition, the matrix A can be decomposed as following � ˜ � 0 B 11 + ρ zz T A = ˜ 0 B 22 Both ˜ B 11 and ˜ B 22 have the form of original matrix with half of the size. We recursively split the matrix until the submatrices contain only one primitive block.

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - Primitive block Since computing the eigensystem of primitive block is a general eigenproblem, QR algorithm will be used The complexity for this step is O ( M 3 ) because the size of primitive block is ( 2 M + 1 ) × ( 2 M + 1 )

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - Conquering step Suppose that Q 1 D 1 Q T 1 and Q 2 D 2 Q T 2 are eigendecompositions of ˜ B 11 and ˜ B 22 respectively. We have: � � Q 1 D 1 Q T 0 + ρ zz T = QDQ T + ρ zz T 1 A = Q 2 D 2 Q T 0 2 � � � � 0 0 Q 1 D 1 where Q = and D = 0 Q 2 0 D 2 Computing this eigendecomposition is just an rank-one modification of symmetric eigenproblem.

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - Rank-one modification of symmetric eigenproblem Bunch et al proposed an algorithm for rank-one update eigenproblem in 1978 However, its computational cost is very expensive So we use an extension of Fast Multipole Method (FMM) to accelerate the algorithm

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - FMM (a) Straight forward scheme. (b) Middle man scheme.

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - FMM Suppose that we want to evaluate the function: n � Φ( y i ) = c j ϕ ( y i − x j ) j = 1 at n points { y i } n i = 1 , where { c j } n j = 1 are constants and ϕ ( x ) is log ( x ) . The direct computation takes O ( n 2 ) . But FMM only takes O ( n + n ) = O ( n ) .

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - Extension of FMM We extend idea of FMM to 2 other functions 1 / x and 1 / x 2

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - Complexity analysis

A divide-and-conquer algorithm for a symmetric eigenproblem Proposed algorithm Proposed algorithm - Complexity analysis Let T ( N ) be complexity of our proposed algorithm for matrix A . We have the recurrence as following: � T ( N ) = 2 T ( N 2 ) + O ( N 2 M 2 ) T ( 1 ) = O ( M 3 ) The solution for this system is O ( NM 3 ) + O ( P 2 ) .

A divide-and-conquer algorithm for a symmetric eigenproblem Numerical experiments Numerical experiments Compare our proposed algorithm (denotes PA) with 3 other algorithms (ADC, CDC, QR) Execution Time Residual Error Orthogonality of eigenvectors

A divide-and-conquer algorithm for a symmetric eigenproblem Numerical experiments Numerical experiments We use 4 sets of test matrix in our numerical experiments: Test 1: The test matrix is generated from the discretization process of Laplace operator in 2D lattices Test 2: Replace all the diagonal block A ii in test 1 by Wilkinson matrix W 2 M + 1 Test 3: Replace the diagonal block A ii in test 1 by by Wilkinson matrix W 2 M + 1 where i is even Test 4: Random symmetric matrix for each diagonal block Since none of our test matrices is particularly, FMM is not used in this comparison.

A divide-and-conquer algorithm for a symmetric eigenproblem Numerical experiments Numerical experiments - Execution Time

A divide-and-conquer algorithm for a symmetric eigenproblem Numerical experiments Numerical experiments - Execution Time

A divide-and-conquer algorithm for a symmetric eigenproblem Numerical experiments Numerical experiments - Residual Error

A divide-and-conquer algorithm for a symmetric eigenproblem Numerical experiments Numerical experiments - Orthogonality Error

A divide-and-conquer algorithm for a symmetric eigenproblem Numerical experiments Numerical experiments Large residual errors occur with PA, ADC and CDC in comparison with QR. We suspect that is because of different strategies in algorithms.

A divide-and-conquer algorithm for a symmetric eigenproblem Numerical experiments Numerical experiments

A divide-and-conquer algorithm for a symmetric eigenproblem Conclusion Conclusion Achievements: Come up with an algorithm with complexity O ( NM 3 ) + O ( P 2 ) for a symmetric tri-block-diagonal problem when FMM is used When FMM is not used, our algorithm still achieves faster running time than other algorithms. Existing problems: Large residual errors in numerical experiments

A divide-and-conquer algorithm for a symmetric eigenproblem Conclusion Q & A

A divide-and-conquer algorithm for a symmetric eigenproblem Conclusion Thank you for your attention!