There are continual and compelling needs for computing many eigenpairs of very large Hermitian matrix in physical simulations and data analysis. Though the Lanczos method is effective for computing a few eigenvalues, it can be expensive for computing a large number of eigenvalues. To improve the performance of the Lanczos method, in this talk, we will present a combination of explicit external deflation (EED) with an s-step variant of thick-restart Lanczos (s-step TRLan). The s-step Lanczos method can achieve an order of s reduction in data movement while the EED enables to compute eigenpairs in batches along with a number of other advantages

Bai, Zhaojun

English

ScholarWorks@UARK

University of Arkansas, Fayetteville ScholarWorks@UARK Mathematical Sciences Spring Lecture Series Mathematical Sciences 4-8-2021 Lecture 08: Partial Eigen Decomposition of Large Symmetric Matrices via Thick-Restart Lanczos with Explicit External Deflation and its Communication-Avoiding Variant Zhaojun Bai University of California, Davis Follow this and additional works at: https://scholarworks.uark.edu/mascsls  Part of the Digital Communications and Networking Commons, Dynamical Systems Commons, Numerical Analysis and Computation Commons, Numerical Analysis and Scientific Computing Commons, and the Partial Differential Equations Commons Citation Bai, Z. (2021). Lecture 08: Partial Eigen Decomposition of Large Symmetric Matrices via Thick-Restart Lanczos with Explicit External Deflation and its Communication-Avoiding Variant. Mathematical Sciences Spring Lecture Series. Retrieved from https://scholarworks.uark.edu/mascsls/14 This Video is brought to you for free and open access by the Mathematical Sciences at ScholarWorks@UARK. It has been accepted for inclusion in Mathematical Sciences Spring Lecture Series by an authorized administrator of ScholarWorks@UARK. For more information, please contact ccmiddle@uark.edu. Partial Eigen Decomposition of Large Symmetric Matricesvia Thick-Restart Lanczos with Explicit External Deflationand its Communication-Avoiding VariantZhaojun BaiUniversity of California, Davis46th Annual Spring Lecture Series, Univ. of ArkansasApril 5-9, 2021 (Virtual)1 / 39CollaboratorsI Jack J. Dongarra, Univ of TennesseeI Chao-Ping Lin, UC DavisI Ding Lu, Univ of KentuckyI Ichitaro Yamazaki, Sandia National Labs.2 / 39Introduction1. Problem statement (“many eigenpairs computation”):Let (λi, vi) be the eigenpairs of a n × n symmetric matrix A with theordering λ1 ≤ λ2 ≤ . . .. Compute the partial eigen decomposition:AVm = VmΛm,where Vm = [v1, v2, . . . , vm] and Λm = diag(λ1, λ2, . . . , λm), n is hugeand m is large.2. Emerging applications such asI electronic structure calculations of Lithium-ion electrolyte,I graphene,I dyanmics analysis of viral capsids of supramolecular systems such as Zika andWest Nile virusesI ...3 / 39Introduction3. Approaches – an incomplete list:I Full eigenvalue decomposition:I LAPACK, ScaLAPACK, PLASMA, MAGMAI ELPA (Eigenvalue Solves for Petaflop Applications) https://elpa.mpcdf.mpg.de/I EigenExa, https://www.r-ccs.riken.jp/labs/lpnctrt/projects/eigenexa/I ...I QDWHeig, [Sukkari, Ltaief, Keyes at KAUST]I SLATE, [Gates et al, U of Tennessee]Stable, but expensive, O(n2) storage and O(n3) flopsI “Spectrum slicing:”I SLEPc, https://slepc.upv.es/I EVSL, http://www.cs.umn.edu/∼saad/softwareI FEAST, http://www.ecs.umass.edu/∼polizzi/feast/I z-Pares, http://zpares.cs.tsukuba.ac.jp/I ...I SISLICE: [Williams-Young, Yang, LBNL]Scalable, but several parameters (e.g. windows) and duplicate/missingeigenvalues on the interface.4 / 39Introduction4. Using Lanczos or any other subspace projection methods for many eigenpairsare challenging – numerically and computationally:I require (internal) locking (deflation) to avoid computing the same eigenvalues(nm2 flops)I often needs large subspace (memory) (e.g., m = 2ne)5 / 39Introduction5. Challenging problem – a case demopreliminary results with EVSL (http://www.cs.umn.edu/∼saad/software)0 0.05 0.1 0.15 0.2 0.24Slices050100150200250300350# of computed eigenvalues0 0.05 0.1 0.15 0.2 0.24Slices050100150200250300350Degrees of polynomial0 0.05 0.1 0.15 0.2 0.24eigenvalue10 -1610 -1410 -1210 -1010 -8relative residue0 0.05 0.1 0.15 0.2 0.24Slices0100020003000400050006000700080009000Wall time of each slice (secs)Figure 1: Left: Image of Dengue virus from https://www.rcsb.org/structure/4cct. Thematrix dimension is 307,260 in the elastic network model. Middle: DOS sliced 5, 076 eigenvaluesinto 20 subintervals (top) and the relative residual numbers for the computed eigenpairs bypolynomial filtered TRLan (bottom). Right: The degrees of filter polynomials for each slices(top) and CPU timing for computing eigenpairs for each slices and the longest one took 8, 772secs (bottom). In addition, the orthogonality of the computed eigenvector V̂ is ‖V̂ T V̂ − I‖F =O(10−6).of-states (DOS) for spectrum slicing. The methods for computing the DOS are the classicalkernel polynomial method and Gaussian quadrature by the Lanczos procedure [46, 45]. EVSLuses a polynomial filtered thick-restart Lanczos (TRLan) to compute the eigenvalues and thecorresponding eigenvectors in each subinterval [43, 42].Case study. We applied EVSL to compute normal modes (i.e., eigenvectors corresponding thesmallest eigenvalues) of an elastic network model for dynamics analysis of viral capsids of severaldifferent viruses in the RCSB Protein Data Bank (https://www.rcsb.org/structure/4cct),also see [60, 61, 38]. We used a compute node of Edison at National Energy Research ScientificComputing Center (NERSC) with OpenMP and the multi-threaded MKL. Figure 1 reports oneof our experimental results. By these experiment results, we observed that (a) EVSL requireshigh degrees of filter polynomials for interior slices. This is also confirmed by the developersof EVSL: “There are situations that very small intervals require a degree of a few thousands.”[42]. (b) To solve multiple subproblems in parallel, there needs huge amount of memory. Forexample, in a test case called riboA of dimension 596,352, if each slice needs 1, 000 Lanczosvectors, solving all the 31 subproblems for 6,082 eigenpairs in parallel would need about 125GB.To cut down the cost of memory (and data communication), we can increase the degree offilter polynomials and/or use restarting strategy. Both strategies could significantly increase thecomputational time. (c) There are also a number of issues related to the interfaces, such as howto handle duplicate and missing eigenvalues and orthogonality between eigenvectors computedfrom different slices and how to validate. (d) More critically, it is still an open problem whether itis possible to completely avoid the factorization of B when computing the DOS of the frequentlyencountered generalized eigenvalue problems Ax = λBx. There is a recent report [66] on usingan iterative procedure to approximate the actions of B−1 and B−1/2 on vectors to address thisissue, see Section 3.1.2 for further discussion.5Left: Dengue virus https://www.rcsb.org/structure/4cct. The matrixdimension is 307,260 in a elastic network model.Middle: DOS sliced 5, 076 eigenvalues into 20 subintervals (top) and therelative residual norms of eigenpairs (bottom). ‖V̂ T V̂ − I‖F = O(10−6).Right: Degrees of filter polynomials (top) CPU timing for each slices6 / 39Introduction6. Many eigenpairs computation is challenging even when vast computationalresources are available.7. This talk is about an on-going development ofLanczos + Explicit external deflation + Communication-avoidingfor computing many eigenpairs.7 / 39Rest of the talkI. Explicit External Deflation (EED)II. Backward stability of EEDIII. Communication-avoiding algorithm for MPK (CA-MPK)IV. Eigensolver (TRLan) with EED and CA-MPKV. Concluding remarksVI. References8 / 39I. Explicit External Deflation (EED)9 / 39I. Explicit External Deflation (EED)1. EED [Hotelling’33, Wielandt’44]Let AV = V Λ be the eigen-decomposition of A, partitionV = [Vk, Vn−k] and Λ =[ΛkΛn−k],where Vk consists of eigenvectors corresponding to first k eigenvalues.DefineÂ = A+ σVkVTkThen(a) The eigenvalues of Â areλi(Â) ={λi + σ for 1 ≤ i ≤ kλi for k + 1 ≤ i ≤ n(b) A and Â have the same eigenvectorsTherefore, after the first k eigenvalues are computed, one can pick a propershift σ to move these eigenvalues away, and then compute the next keigenvalues of A using the matrix Â.10 / 39I. Explicit External Deflation (EED)2. Governing equations of j steps of exact EED process: for j = 1, 2, . . . ,Aj = Aj−1 + σjvjvTj = A+ VjΣjVTj and A0 = A.Ajvj+1 = λj+1vj+1AjVj = Vj(Λj +Σj)with initial Av1 = λ1v1, where Σj = diag(σ1, . . . , σj).3. From the governing equations, we have the desired partial eigendecompositionAVj = VjΛj .11 / 39I. Explicit External Deflation (EED)4. Advantages of EED:I No need to explicitly orthogonalize the projection subspace to the computedeigenvectors in contrast to the implicit internal deflation (“locking”)I Easily incooperated into existing eigensolvers, such as TRLan and ARPACK.I Straightforward extension to generalized symmetric eigenproblem Av = λBv:(A+ σBVkVTk B)x = λBxwithout factoring/inverting the matrix B.I Readily exploit structures of A and B, for example, see EED for the linearresponse eigenvalue problem [Bai-Li-Lin’17].I Extensible to other eigenvalue-type problems, such as SVD.I Nonlinear eigenvector problem (NEPv) analogy: level shifting12 / 39I. Explicit External Deflation (EED)5. Two key issues of EED:(a) Numerical stability with approximate eigenvectors V̂k:Âj = A+ V̂jΣj V̂Tj(b) Cost of matrix powers kernel (MPK) for generating Krylov subspace:[p0(Âj)v0, p1(Âj)v0, . . . , ps(Âj)v0]where {pk(·)} are recursively defined polynomials.13 / 39II. Backward stability of EED14 / 39II. Backward stability of EED1. Previous work: [Wilkinson’65], [Parlett’82], [Saad’89], [Jang/Lee’06]2. Governing equations of j steps inexact EED:Âj = Âj−1 + σj v̂j v̂Tj = A+ V̂jΣj V̂TjÂj v̂j+1 = λ̂j+1v̂j+1 + ηj+1Âj V̂j = V̂j(Λ̂j +Σj) + V̂jΣjΦj + Ejwhere Â0 = A,‖ηj+1‖ ≤ tol · ‖A‖Λ̂j = diag(λ̂1, . . . , λ̂j)V̂j = [v̂1, . . . , v̂j ],Σj = diag(σ1, . . . , σj),Φj = utri(V̂Tj V̂j − Ij),Ej = [η1, η2, . . . , ηj ].and tol is prescribed relative residual tolerance.15 / 39II. Backward stability of EED3. Notion of backward stability for computed eigenpairs (Λ̂j+1, V̂j+1) withrelative residual tolerance tol:I The loss of orthogonalityωj+1 = ‖V̂ Tj+1V̂j+1 − I‖F = O(tol)I The backward errorδj+1 = min∆∈H‖∆‖F = O(tol · ‖A‖),whereH ={∆ | (A+∆)Qj+1 = Qj+1Λ̂j+1, ∆ = ∆T , Qj+1 = orth(V̂j+1)}16 / 39II. Backward stability of EED4. Two key quantitiesI Spectral gapγj ≡ minλ∈Ij+1,θ∈Jj|λ− θ|,where Ij+1 = {λ̂1, . . . , λ̂j , λ̂j+1} is the set of computed eigenvalues, andJj = {λ̂1 + σ1, . . . , λ̂j + σj} is the set of computed eigenvalues with shifts.I Shift-gap ratioτj ≡1γj· max1≤i≤j|σi|.5. Theorem (Rule of Thumb)Under mild assumptions, if γ−1j ‖A‖ = O(1) and τj = O(1), thenωj+1 = O(tol) and δj+1 = O(tol · ‖A‖).6. Dynamical choice of shifts σj to satisfy the conditions γ−1j ‖A‖ = O(1) andτj = O(1).17 / 39II. Backward stability of EED7. Numerics of TRLED = TRLan + EED:I Test matrices:matrix n [λmin, λmax] [λlow, λupper] neLaplacian 40, 000 [0, 7.9995] [0, 0.07] 205worms20 20, 055 [0, 6.0450] [0, 0.05] 289SiO 33, 401 [−1.6745, 84.3139] [−1.7, 2.0] 182Si34H36 97, 569 [−1.1586, 42.9396] [−1.2, 0.4] 310Ge87H76 112, 985 [−1.214, 32.764] [−1.3,−0.0053] 318Ge99H100 112, 985 [−1.226, 32.703] [−1.3,−0.0096] 37218 / 39II. Backward stability of EED7. Numerics of TRLED = TRLan + EED, cont’dI Results:matrix n̂e jmax ωn̂e ‖Rn̂e‖F/AnormCPU time (sec.)TRLED TRLanLaplacian 205 60 1.93 · 10−8 6.33 · 10−8 66.5 86.0worms20 289 86 2.63 · 10−8 7.24 · 10−8 57.3 74.8SiO 182 41 2.33 · 10−8 4.71 · 10−8 42.4 47.1Si34H36 310 72 3.41 · 10−8 7.50 · 10−8 309.9 310.4Ge87H76 318 66 4.08 · 10−8 8.50 · 10−8 388.7 421.0Ge99H100 372 74 3.65 · 10−8 7.63 · 10−8 501.1 533.4I EED profile of Ge99H100:19 / 39II. Backward stability of EED7. Numerics of TRLED = TRLan + EED, cont’dI Observations:(1) All desired eigenvalues are successfully computed: ne = n̂e(2) computed eigenpairs are backward stable:ωne = O(tol) and δne ≈ ‖Rne‖ = O(tol · ‖A‖).(3) EED steps do not significant slow down in execution time, in fact, slightly improveddue to low “internal” memory usage.20 / 39III. Communication-avoiding algorithm for computing MPK21 / 39III. Communication-avoiding algorithm for computing MPK1. Sparse-plus-low-rank matrix powers kernel (MPK)[p0(B)v, p1(B)v, . . . , ps(B)v]withB = A+ σUkUTk = sparse + low rankand AUk = UkΛk, UTk Uk = Ik and pj(·) are polynomials defined recursively.2. For simplicity, consider polynomials {pj(·)} in the monomial basis andcompute the MPK:[p0(B)v, p1(B)v, . . . , ps(B)v] =[x,Bx,B2x, . . . , Bsx]≡[x(0), x(1), x(2), . . . , x(s)]3. Standard MPK algorithm1: x(0) = x;2: for j = 1 : s do3: x(j) = Ax(j−1) + σ(Uk(UTk x(j−1)))4: end for22 / 39III. Communication-avoiding algorithm for computing MPK4. ExampleI B = A+ σUkUTk , where A is 2D Laplacian and Uk are eigenvectorsI matrix powers kernel Vs = [p0(B)v0, p1(B)v0, . . . , ps(B)v0]I Timing of the standard algorithm in MATLAB – blue curve.n = 250 000 and k = 10000 5 10 15 20s-value00.511.522.533.544.5tic-toc in secondscpu timing of sparse+low rank matrix powers kernelstd recurrencebc alg.I Q: what’s the red curve?23 / 39III. Communication-avoiding algorithm for computing MPK5. An algorithm has two costs:arithmetic (flops) + movement of data (communication)6. Communication is the bottleneck on modern architectures.7. Q: How to exploit the sparse-plus-low-rank matrix structure to reducecommunication cost?Ans.: use a specialized communication-avoiding algorithm developed byLeiserson-Rao-Toledo’97 (“out-of-core”) and Knight-Carson-Demmel’13(“exploiting data sparsity”)24 / 39III. Communication-avoiding algorithm for computing MPK8. Communication-Avoiding (CA) algorithm for the MPK1: x(0) = x2: b0 = UTk x(0)3: W(j)k = Λjk + σ∑ji=1 Λi−1k (Λk + σ)j−i for j = 1 : s− 1,4: bj =W(j)k b0 for j = 1 : s− 15: [c0, c1, . . . , cs−1] = Uk[b0, b1, . . . , bs−1]6: for j = 1 : s do7: x(j) = Ax(j−1) + σcj−18: end for9. Benefits of CA-MPK algorithm:I Reduced flopsflopsstd = nnzA · s+ nks+ nksflopsca = nnzA · s+ nks+ nk +O(k2s)I Reduced movement of data: Uk is only accessed twice.25 / 39III. Communication-avoiding algorithm for computing MPK10. Derivation of a simplified blocking covers (communication-avoiding) MPKalgorithmI Recall the identity Cj = (E + F )j = Ej +∑ji=1Ei−1FCj−iI Let B = A + σUkUTk and AUk = UkΣk , then for j = 1, . . . , s,Bj=(A + σUkUTk)j= Aj+ σj∑i=1Ai−1UkUTk Bj−i= Aj+ σj∑i=1UkΛi−1kUTk Bj−iI Then for j = 0, 1, . . . , s − 1,UTk Bjx = UTk Ajx + σj∑i=1Λi−1kUTk Bj−ix = ΛjkUTk x + σj∑i=1Λi−1k(Λk + σ)j−iUTk x=Λjk+ σj∑i=1Λi−1k(Λk + σ)j−iUTk x ≡ W (j)k · UTk x ≡ bjI Compute [c0, c1, . . . , cs−1] = Uk[b0, b1, . . . , bs−1]I Finally for j = 1, . . . , s,x(j) ≡ Bjx = B(Bj−1x) = (A + σUkUTk )(Bj−1x)= A(Bj−1x) + σUkUTk (Bj−1x) = A(Bj−1x) + σUk(UTk Bj−1x)= Ax(j−1)+ σUkbj−1 = Ax(j−1)+ σcj−126 / 39III. Communication-avoiding algorithm for computing MPK11. Example, cont’d:I B = A+ σUkUTk , where A is 2D Laplacian and Uk are eigenvectorsI Matrix powers kernel Vs = [p0(B)v0, p1(B)v0, . . . , ps(B)v0]I Timing in MATLABI Standard algorithm – blue curveI CA algorithm – red curven = 250000 and k = 10000 5 10 15 20s-value00.511.522.533.544.5tic-toc in secondscpu timing of sparse+low rank matrix powers kernelstd recurrencebc alg.27 / 39III. Communication-avoiding algorithm for computing MPK11. Example, cont’d:I B = A+ σUkUTk , where A is 2D Laplacian and Uk are eigenvectorsI Matrix powers kernel Vs = [p0(B)v0, p1(B)v0, . . . , ps(B)v0]I Timing in one node (32 cores) of Cori (NERSC)k = 100 k = 20012. Rounding error analysis of CA-MPK – skip.28 / 39IV. Lanczos algorithm with EED and MPK29 / 39IV. Lanczos algorithm with EED and MPK1. Lanczos algorithmI Lanczos processAQj = QjTj + βjqj+1eTjwhere span{Qj} = span{q,Aq, . . . , Aj−1q} (Krylov subspace)I Rayleigh-Ritz approximationTjxi = θixi(λi, vi) ≈ (θi, Qjxi)2. Lacnzos algorithm is efficient for computing a few exterior eigenvalues (andeigenvectors).3. Two main kernelsI Matrix-Vector multiply (SpMV) for generating Krylov subspace: AqI Re-orthogonalization for maintaining orthonormal basis vectors: Qj4. Two variants of Lanczos method:I Thick-restart Lanczos (TRLan) – to reduce memory requirementsI s-step Lanczos (s-Lanczos) – replace SpMV by MPK (Matrix Powers Kernel)to reduce communication cost30 / 39IV. Lanczos algorithm with EED and MPK5. sTRLED = s-step-TRLan + EED + CA-MPK31 / 39IV. Lanczos algorithm with EED and MPK6. Preliminary results on performance ofsTRLED = s-step-TRLan + EED + CA-MPKSpeedup on single core for Si34H26 with Standard, Block-cover andSpecialized MPKs:2452/240/0 +229/03186/31+2/0-1/03623/35+1/00/04125/40-1/00/04440/440/00/04324/42-86/0+7/04746/480/00/02.2x 2.2x 1.8x 2.2x 1.9x 2.3x 1.9x 2.2x 1.9x 2.3x 1.9x 2.3x 1.9x 2.3x100 200 300 400 500 600 700# of eigenpairs050100150200250300350400450500550time (s)DeflateSpMVReOrthoRestartn = 97569, ne = 100, . . . , 700computing 100 eigenvalues at a time with m = 200 and s = 532 / 39IV. Lanczos algorithm with EED and MPK7. Preliminary results on strong-parallel scaling ofsTRLED = s-step-TRLan + EED + CA-MPKspeedup over standard MPK on one processor for Si87H76:3.1x3.1x3.0x2.8x 2.8x 2.7x 2.7x 2.5x1  2  4  8  16 32 64 128Process count0100020003000400050006000700080009000Time to solution (s)DeflateSpMVReOrthoRestart14 8 16 32 64 128Process count020406080100120140Speedup over standard on one process standardspecializedn = 240, 369, ne = 700computing 100 eigenvalues at a time with m = 200 and s = 533 / 39IV. Lanczos algorithm with EED and MPK7. Preliminary results on strong-parallel scaling ofsTRLED = s-step-TRLan + EED + CA-MPKspeedup over standard MPK 2D Laplacian:n = 250, 000, ne = 5, 000,computing 500 eigenvalues at a time with m = 1000 and s = 534 / 39V. Concluding remarks35 / 39V. Concluding remarks1. Many eigenpairs computation in emerging applications, a challengingproblem even when vast computational resources are available.2. Two techniques discussed in this talk:I Explicit external deflation (EED) for reliably moving away computedeigenpairs,I a communication-avoiding matrix powers kernel (CA-MPK) for fastsparse-plus-low-rank MPK3. The capability of being able to efficiently compute large number ofeigenvalues will not just be appealing, but also mandatory for the nextgeneration of eigensolvers.36 / 39VI. References37 / 39VI. References1. B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM Classics Edition,1998.2. Y. Saad, Numerical Methods for Large Eigenvalue Problems, Revised edition,SIAM, 2011.3. Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. A. van der Vorst, eds.Templates for the Solution of Algebraic Eigenvalue Problems: A PracticalGuide, SIAM, 2000.4. C. Leiserson, S. Rao, S. Toledo, Efficient out-of-core algorithms for linearrelaxation using blocking covers, J. Comput. Syst. Sci. Int. Vol.54,pp.332-–344, 1997.5. N. Knight, E. Carson, J. Demmel, Exploiting data sparsity in parallel matrixpowers computations, PPAM, pp.15–25, 2014.38 / 39VI. References6. S. Kim, A. T. Chronopoulos, A class of Lanczos-like algorithms implementedon parallel computers, Parallel Computing, Vol.17, pp.763–778, 1991.7. I. Yamazaki, K. Wu, A communication-avoiding thick-restart Lanczosmethod on a distributed-memory system, in: Workshop on Algorithms andProgramming Tools for next-generation high-performance scientific andsoftware (HPCC), pp.345–354, 2011.8. E. Carson and J. W. Demmel, Accuracy of the s-step Lanczos method forthe symmetric eigenproblem in finite precision, SIAM J. Matrix Anal. Appl.Vol.36, pp.783-819, 2015.9. Z. Bai, J. Dongarra, D. Lu and I. Yamazaki, Matrix powers kernels forthick-restart Lanczos with explicit external deflation, IPDPS, May 201910. C. Lin, D. Lu and Z. Bai, Backward stability analysis of explicit externaldeflation for symmetric eigenvalue problems, to appear soon39 / 39

Lecture 08: Partial Eigen Decomposition of Large Symmetric Matrices via Thick-Restart Lanczos with Explicit External Deflation and its Communication-Avoiding Variant

https://scholarworks.uark.edu/cgi/viewcontent.cgi?article=1012&context=mascsls

Lecture 08: Partial Eigen Decomposition of Large Symmetric Matrices via Thick-Restart Lanczos with Explicit External Deflation and its Communication-Avoiding Variant

Abstract

Similar works

Full text

Available Versions

ScholarWorks@UARK