Finding Rightmost Eigenvalues of Large, Sparse, Nonsymmetric Parameterized Eigenvalue ProblemsMotivationProblem StatementMethodArnoldi Algorithm (continue)Matrix TransformationMatlab Code of Arnoldi MethodTest ProblemTest Problem (continue)Slide 10Implicitly Restarted Arnoldi (IRA)IRA (continue)Slide 13Slide 14Slide 15Future Work (AMSC 664)Finding Rightmost Eigenvalues of Large, Sparse, Nonsymmetric Parameterized Eigenvalue ProblemsMinghao WuAMSC [email protected]: Dr. Howard ElmanDepartment of Computer [email protected]•To determine the stability of the linearized system of the form:AxxB • The steady state solution x* is - stable, if all the eigenvalues of Ax = λBx have negative real parts; - unstable, otherwise.3Problem StatementTo find the rightmost eigenvalues of:BxAxwhere matrices A and B are • Real N-by-N• Large• Sparse• Nonsymmetric• Depend on one or several parameters4Method•Eigensolver: Arnoldi Algorithm - Iterative method - Based on Krylov subspace:   1112111, uAuAAuuuAKkk  - Computes Arnoldi decomposition:TkkkkkkeuHUAU1where:[Uk uk+1]: an orthonormal basis of Hk: k-by-k upper Hessenberg matrix, k « Nβk: scalar, ek: k-by-1 vector [0 0 … 0 1]’  111uAAuuk5Arnoldi Algorithm (continue)•Eigenvalues of Hk approximate eigenvalues of A Premultiply previous equation by transpose(Uk):.1 kTkkTkkkkTkkTkHeuUHUUAUU Let (λ,z) be an eigenpair of Hk, then:   .zUλzUAkkTkkTkkTkUzUUzzHU As k increases, ||Residual|| will decrease.When k = N, Residual = 0.    0 zUzUAUkkTkResidual6Matrix Transformation•Motivation - Arnoldi Algorithm cannot solve generalized eigenvalue problem - It converges to well – separated extremal eigenvalues, not rightmost eigenvalues•Shift – Invert Transformation   1;,1xxTBxAxBBABATSISI7Matlab Code of Arnoldi Method•Arnoldi Algorithm (with Shift – Invert matrix transformation) routine: [v,X,U,H]=SI_Arnoldi(A,B,k,sigma) Input: A, B: matrix A and B in Ax = λBx k: number of eigenpairs wanted sigma: the shift σ in shift – invert matrix transformation Output: v: a vector of k computed eigenvalues X: k eigenvectors associated with the eigenvalues U: the Krylov basis Uk+1 H: the upper Hessenberg matrix Hk8Test ProblemOlmstead Model (see Olmstead et al (1986)): SucbSuRucuSutxxxxt13with boundary conditions: ,0,0  xSuThis model represents the flow of a layer of viscoelastic fluid heated from below. u: the speed of the fluidS: related to viscoelastic forcesb,c: scalars, R: scalar, Rayleigh number9Test Problem (continue)•Discretize the model with finite differences - grid – size: h = 1/(N/2) - (4) can be written as dy/dt = f(y) with] .[2/2/2211 NNTSuSuSuy • Evaluate the Jacobian matrix A = df/dy at steady state solution y* - N = 1000, b = 2, c = 0.1, R = 0.6 - y* = 0 - A = df/dy(y*) is a nonsymmetric sparse matrix with bandwidth 610Test Problem (continue)Computational Result:•Rightmost eigenvalues: λ1,2 = 0 ± 0.4472i•Residual: ||Axi - λi xi|| = 8.4504e-012, i=1,2•The result agrees with the literature.11Implicitly Restarted Arnoldi (IRA)•Motivation - Large k is not practical Example: Size of A: 10,000Value of k: 100Memory required tostore U100 in double precision:10 MB- When B is singular, Arnoldi algorithm may give rise to spurious eigenvalues12IRA (continue)•Basic idea of Implicitly Restarted Arnoldi Filter out the unwanted eigendirections from the starting vector by using the most recent spectrum information and a clever filtering technique•IRA steps 1. Compute m eigenpairs (k<m«N) by Arnoldi method with starting vector u1 2. Filter out the m-k unwanted eigendirections from u1 (Key Technique: shifted QR algorithm) 3. Restart the process with filtered starting vector till the k eigenvalues of interest converge13Test ProblemxMxCCKT0000K: 200-by-200 matrix, full rank;C: 200-by-100 matrix, full rank;M: 200-by-200 matrix, full rank.Eigenvalue problem with this kind of block structure appears in the stability analysis of steady state solution of Navier – Stokes equations for incompressible flow.14Test Problem (continue)•Use Matlab function “rand” to generate K, C, M •-2.7377 ≤ Re(λ) ≤ 49.9129•Find out 10 rightmost eigenvalues•Use the IRA code written by Fei Xue15Test Problem (continue)Exact Eigenavalues49.9129 + 0i2.9112 + 1.1256i2.9112 - 1.1256i2.5036 + 0.0624i2.5036 - 0.0624i2.3792 + 0i2.1318 + 0.9356i2.1318 - 0.9356i2.1081 + 1.3539i2.1081 - 1.3539iComputed Eigenavalues(Arnoldi)49.9129 + 0i193.8412 - 7113830.9524i193.8412 + 7113830.9524i3.0891 + 0i2.6112 + 0i-0.5752 - 2.7079i-0.5752 + 2.7079i-1.0901 - 0.7532i-1.0901 + 0.7532i-47.3106 + 0iComputed Eigenavalues(IRA)49.9129 + 0i2.9112 - 1.1256i2.9112 + 1.1256i2.5036 - 0.0624i2.5036 + 0.0624i2.3792 + 0i2.1318 - 0.9356i2.1318 + 0.9356i2.1081 - 1.3539i2.1081 + 1.3539iComputational Result: (shift σ = 60)16Future Work (AMSC 664)•Solve the third test problem•Implement iterative solvers for linear