A Robust and Efficient Parallel SVD Solver Based on Restarted - PowerPoint PPT Presentation

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation A Robust and Efficient Parallel SVD Solver Based on Restarted Lanczos Bidiagonalization Jose E. Roman Universidad Polit´ ecnica de Valencia, Spain (joint work with V. Hernandez and A. Tomas) Harrachov 2007

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Outline Overview of SLEPc 1 Summary of Functionality The SVD in SLEPc Restarted Lanczos Bidiagonalization 2 Lanczos Bidiagonalization Dealing with Loss of Orthogonality Restarted Bidiagonalization Enhancements for Parallel Efficiency Evaluation 3 Numerical Robustness Parallel Performance

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation What Users Need Provided by PETSc ◮ Abstraction of mathematical objects: vectors and matrices ◮ Efficient linear solvers (direct or iterative) ◮ Easy programming interface ◮ Run-time flexibility, full control over the solution process ◮ Parallel computing, mostly transparent to the user Provided by SLEPc ◮ State-of-the-art eigensolvers and SVD solvers ◮ Spectral transformations

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Summary PETSc : Portable, Extensible Toolkit for Scientific Computation Software for the scalable (parallel) solution of algebraic systems arising from partial differential equation (PDE) simulations ◮ Developed at Argonne National Lab since 1991 ◮ Usable from C, C++, Fortran77/90 ◮ Focus on abstraction, portability, interoperability ◮ Extensive documentation and examples ◮ Freely available and supported through email Current version: 2.3.3 (released May 2007) http://www.mcs.anl.gov/petsc

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Summary SLEPc : Scalable Library for Eigenvalue Problem Computations A general library for solving large-scale sparse eigenproblems on parallel computers ◮ For standard and generalized eigenproblems ◮ For real and complex arithmetic ◮ For Hermitian or non-Hermitian problems Current version: 2.3.3 (released June 2007) http://www.grycap.upv.es/slepc

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation PETSc/SLEPc Numerical Components PETSc SLEPc Nonlinear Systems Time Steppers SVD Solvers Trust Pseudo Line Backward Cross Matrix Lanczos Thick Res. Cyclic Other Euler Other Region Time Step Search Euler Product Lanczos Krylov Subspace Methods Eigensolvers GMRES CG CGS Bi-CGStab TFQMR Richardson Chebychev Other Krylov-Schur Arnoldi Lanczos Other Preconditioners Spectral Transform Additive Block Jacobi ILU ICC LU Other Shift Shift-and-invert Cayley Fold Schwarz Jacobi Matrices Compressed Block Compressed Block Dense Other Sparse Row Sparse Row Diagonal Index Sets Vectors Indices Block Indices Stride Other

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Singular Value Decomposition (SVD) A = U Σ V ∗ where ◮ A is an m × n rectangular matrix ◮ U = [ u 1 , u 2 , . . . , u m ] is a m × m unitary matrix ◮ V = [ v 1 , v 2 , . . . , v n ] is a n × n unitary matrix ◮ Σ is a m × n diagonal matrix with entries Σ ii = σ i ◮ σ 1 ≥ σ 2 ≥ · · · ≥ σ n ≥ 0 ◮ If A is real, U and V are real and orthogonal ◮ Each ( σ i , u i , v i ) is called a singular triplet

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Thin SVD V ∗ n = U n Σ n A In SLEPc we compute a partial SVD, that is, only a subset of the singular triplets

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation SVD Solvers Based on EPS Cross product matrix A ∗ Ax = λx AA ∗ y = λy The eigenvalues are λ i = σ 2 i and the eigenvectors x i = v i or y i = u i Cyclic matrix � 0 � A H ( A ) x = λx H ( A ) = A ∗ 0 � ± u i � 1 The eigenvalues are ± σ i and the eigenvectors √ v i 2

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Basic Usage Usual steps for solving an SVD problem with SLEPc: 1. Create an SVD object 2. Define the SVD problem 3. (Optionally) Specify options for the solution 4. Run the SVD solver 5. Retrieve the computed solution 6. Destroy the SVD object All these operations are done via a generic interface, common to all the SVD solvers

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Simple Example SVD svd; /* singular value solver context */ Mat A; /* matrix */ Vec u, v; /* singular vectors */ PetscReal sigma; /* singular value */ SVDCreate(PETSC_COMM_WORLD, &svd); SVDSetOperator(svd, A); SVDSetFromOptions(svd); SVDSolve(svd); SVDGetConverged(svd, &nconv); for (i=0; i<nconv; i++) { SVDGetSingularTriplet(svd, i, &sigma, u, v); } SVDDestroy(svd);

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Run-Time Examples % program -svd_type lanczos -svd_tol 1e-12 -svd_max_it 200 % program -svd_type trlanczos -svd_nsv 4 % program -svd_type cross -svd_eps_type krylovschur -svd_ncv 30 -svd_smallest -svd_monitor_draw % program -svd_type cyclic -svd_eps_type arpack -svd_st_type sinvert -svd_st_shift 1 % mpirun -np 16 program ...

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Bidiagonalization Compute the SVD in two stages [Golub and Kahan, 1965]: Q ∗ 1. A = PBQ ∗ = A P B 2. B = X Σ Y ∗ , with U = PX and V = QY Lanczos bidiagonalization computes part of the info: P k , B k , Q k → Ritz approximations: ˜ σ i , u i = P k x i , ˜ ˜ v i = Q k y i

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Lanczos Bidiagonalization Equating the first k columns AQ k = P k B k A ∗ P k = Q k B ∗ k + β k q k +1 e ∗ k   α 1 β 1 α 2 β 2     α 3 β 3 α j = p ∗ j Aq j   B k =   ... ... β j = p ∗ j Aq j +1       α k − 1 β k − 1   α k Double recursion: α j p j = Aq j − β j − 1 p j − 1 , β j q j +1 = A ∗ p j − α j q j

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Golub-Kahan-Lanczos Bidiagonalization Golub-Kahan-Lanczos algorithm Choose a unit-norm vector q 1 Set β 0 = 0 For j = 1 , 2 , . . . , k p j = Aq j − β j − 1 p j − 1 α j = � p j � 2 p j = p j /α j q j +1 = A ∗ p j − α j q j β j = � q j +1 � 2 q j +1 = q j +1 /β j end Loss of orthogonality is an issue

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Algorithm with Orthogonalization Lanczos bidiagonalization with orthogonalization Choose a unit-norm vector q 1 For j = 1 , 2 , . . . , k p j = Aq j Orthogonalize p j with respect to P j − 1 α j = � p j � 2 p j = p j /α j q j +1 = A ∗ p j Orthogonalize q j +1 with respect to Q j β j = � q j +1 � 2 q j +1 = q j +1 /β j end

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation One-Sided Orthogonalization Orthogonalizing right vectors is enough [Simon and Zha, 2000] One-Sided Lanczos bidiagonalization Choose a unit-norm vector q 1 Set β 0 = 0 For j = 1 , 2 , . . . , k p j = Aq j − β j − 1 p j − 1 α j = � p j � 2 p j = p j /α j q j +1 = A ∗ p j Orthogonalize q j +1 with respect to Q j β j = � q j +1 � 2 q j +1 = q j +1 /β j end

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Restarted Bidiagonalization Required k can be arbitrarily large (slow convergence, many singular triplets) Problems: storage and computational effort Solution: restart the computation when a certain k is reached Explicit restart: re-run with a “better” q 1 (e.g. use Ritz vector associated to the first value) Thick restart: a better alternative that avoids to explicitly compute a new initial vector Idea: keep ℓ -dimensional subspace with relevant spectral information [Baglama and Reichel, 2005]

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Thick-Restart Lanczos Bidiagonalization Compact Lanczos Bidiagonalization of order ℓ + 1 : A ˜ Q ℓ +1 = ˜ P ℓ +1 ˜ B ℓ +1 A ∗ ˜ P ℓ +1 = ˜ Q ℓ +1 ˜ ℓ +1 + ˜ B ∗ β ℓ +1 ˜ q ℓ +2 e ∗ k +1 ˜ Q ℓ +1 = [˜ v 1 , ˜ v 2 , . . . , ˜ v ℓ , q k +1 ] ˜ v i = Q k y i right Ritz vectors Residual of full decomposition ˜ P ℓ +1 = [˜ u 1 , ˜ u 2 , . . . , ˜ u ℓ , ˜ p ℓ +1 ] u i = P k x i left Ritz vectors ˜ New left initial vector f = Aq k +1 − � ℓ p ℓ +1 = f/ � f � , ˜ i =1 ˜ ρ i ˜ u i , α ℓ +1 = � f � ˜ ˜ ˜ q ℓ +2 = g/ � g � , g = A ∗ ˜ p ℓ +1 − ˜ α ℓ +1 q k +1 , β ℓ +1 = � g �

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Thick-Restart Lanczos Bidiagonalization Compact Lanczos Bidiagonalization of order ℓ + 1 : A ˜ Q ℓ +1 = ˜ P ℓ +1 ˜ B ℓ +1 A ∗ ˜ P ℓ +1 = ˜ Q ℓ +1 ˜ ℓ +1 + ˜ B ∗ β ℓ +1 ˜ q ℓ +2 e ∗ k +1  ˜ σ 1 ρ 1 ˜  σ 2 ˜ ρ 2 ˜    .  ... ˜ . B ℓ +1 =   .     σ ℓ ˜ ρ ℓ ˜   α ℓ +1 ˜ σ i are Ritz values ˜ f = Aq k +1 − � ℓ p ℓ +1 = f/ � f � , ˜ i =1 ˜ ρ i ˜ u i , α ℓ +1 = � f � ˜

Overview of SLEPc Restarted Lanczos Bidiagonalization Evaluation Thick-Restart Lanczos Bidiagonalization The decomposition can be extended to order k and the Lanczos relations are maintained   ˜ σ 1 ρ 1 ˜ σ 2 ˜ ρ 2 ˜   .  ...  .  .      σ ℓ ˜ ρ ℓ ˜ ˜   B k =   α ℓ +1 β ℓ +1 ˜    ... ...        α k − 1 β k − 1   α k