Krylov Subspaces

1. Krylov Subspaces Final Presentation ME 535 Computation Techniques in Mechanical Engineering Spring Quarter 2014 Chad C. Beck

2. References • Lloyd N. Trefethen and David Bau, III. Numerical Linear Algebra. SIAM, 1997. p243-283,293-301. • David S. Watkins. The Matrix Eigenvalue Problem. SIAM, 2007. p351-421. • Ulrich Hetmaniuk.Applied Linear Algebra and Intro to Numerical Analysis. lecture notes. University of Washington. Autumn 2013. • Nordsveen, Magnus, and Randi Moe. "Preconditioned Krylov Subspace Methods Used in Solving Two-Dimensional Transient Two-Phase Flows." International Journal for Numerical Methods in Fluids. 31.7 (1999): 1141-1156.

3. Krylov Sequence and Subspace • The Krylov sequence was created by Russian mathematician and engineer Alexei Krylov. • A given matrix A and vector b are multiplied to form a Krylov sequence as b, Ab, A(Ab), A(A(Ab)), … Akb. • The Krylov subspaces are the spaces spanned by the successively larger groups of these vectors

4. The need for Krylov Subspaces • Direct algorithms require on the order of m3 operations for a matrix of m-dimension. This significantly limits the ability to work with larger matrices. • Most matrices in engineering applications are sparsely populated, such as those produced from Finite Element Method (FEM) or Finite Difference Method (FDM). When direct methods are used on dense matrices, sparsity is often lost, as zero elements above/below a diagonal become non-zero. • Iterative solutions that use Krylov Subspaces can reduce the order of operation to m, a drastic improvement over m3. They also can maintain whatever sparsity exists within a matrix.

5. Krylov Subspaces • The idea is that it is less costly to look for approximations to the solution space x or the eigenvalue λ in the Krylov subspace that minimize the residual than perform a direct method, such as QR factorization. • However, the vectors b, Ab, A2b, … Akbcan quickly become linearly dependant, forming a poor basis. Orthonormalization is usually required to form an orthogonal basis

6. Krylov Subspace Iterative Methods Given a matrix A and a vector b, the associated Krylov sequence is the set of vectors b, Ab, A2b, A3b,… In each of the below methods, the idea is to project an m-dimensional problem into a lower-dimensional Krylov subspace. A* is the Hermitian conjugate of A. Table reproduced from “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

7. Arnoldi Iteration - Theory The matrices Qk generated by the Arnoldi iteration are reduced QR-factors of the Krylov matrix: The Hessenberg matrices Hk are the corresponding projections And the successive iterates are related by the formula Source: “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

8. Arnoldi Iteration - Eigenvalues • The eigenvalues of Hk are called Ritz values and are approximations of the eigenvalues of the matrix A. • Qk y are called Ritz vectors, where yis eigenvectors of Hk, and approximate eigenvectors of A. • The Ritz values converge to the extreme eigenvalues of A. These are usually the eigenvalues of interest for most Physical applications.

9. Arnoldi Iteration – Pseudo Code Arbitrary non-zero vector x0 q1 = x0/(||x0||2) for k = 1,2,… uk = Aqk for j = 1 to k hjk = qk*uk uk = uk– hjkqj end hk+1,k= ||uk||2 qk+1 = uk / hk+1,k end Adapted from “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

10. Lanczos Iteration - Theory The matrices Qk generated by the Arnoldi iteration are reduced QR-factors of the Krylov matrix: The tridiagonal matrices Tk are the corresponding projections And the successive iterates are related by the formula Source: “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

11. Lanczos Iteration - Theory • Analogous to Arnoldi The eigenvalues of Tk approximations of the eigenvalues of the matrix A, and Qk y are approximate eigenvectors of A, where yis eigenvectors of Hk,. • The Lanczos iteration builds a series of tridiagonal matrices Adapted from “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

12. Lanczos Iteration – Pseudo Code β0 = 0, q0= 0, q1 = b/(||b||2) (b is arbitrary) for k = 1,2,… v = Aqk αk= qk*v v = v – βk-1qk -1- αkqk βk = ||v||2 qk+1 = v / βk end Adapted from “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

13. GMRES Iteration - Theory • The GMRES iteration approximates the solution to Ax = b by finding a vector xk within the Krylov subspace that minimizes the norm of the residual. • The vector xk is found by solving the least squares problem, where x0is usually zero. Adapted from Prof. Hetmaniuk’s lecture notes for AMATH 584, Univ. of Washington

14. GMRES – Graphical Interpretation The least squares polynomial approximation problem of GMRES Figure from “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

15. GMRES Method • In the k-th step of the Arnoldi iteration, the norm of the residual is minimized • The solution is found with the singular value decomposition (SVD) of Hk and Adapted from Prof. Hetmaniuk’s lecture notes for AMATH 584, Univ. of Washington

16. GMRES Iteration – Pseudo Code x0 = 0, h10 = ||b||2 , k = 0 q1 = b / h10 for k = 1,2,… u = Aqk for j = 1 to k hjk = qk*u u = u – hi,kqi end hk+1,k= ||u||2 xk = Qkykwhere ||h10e1- Hkyk||2 is minimized. end Adapted from “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

17. Conjugate Gradient (CG) - Theory • When A is hermitianpositvedefinate, the CG iteration approximates the solution to Ax = b by minimizes the norm of the residual over the Krylov subspace. • The method can be summarized as follows: Minimizing the norm of the residual is accomplished by minimizing an associated quadratic equation. From an initial guess for x, a residual is found, then a direction to search, and then a new approximation for x.The CG method actually builds two bases for the Krylov space; one for the residual, and one for the search direction.

18. CG Iteration – Pseudo Code r0 = b – Ax0 for k = 1,2,… until convergence βk= r*k-1rk-1/ r*k-2rk-2 [β0 = 0] pk = rk-1 + βkpk-1 [p1= r0] αk = r*k-1rk-1/ p*kApk xk = xk-1 + αkpk rk= rk-1 - αkApk end Adapted from “Numerical Linear Algebra” by L. N. Trefethen & D. Bau

19. Application of Krylov Subspace Method M. Nordsveen and R. Moeof Norway used Krylov Subspaces in their study of a two-fluid flow model. From abstract: “A two-fluid model was developed for the simulation of separated and intermittent gas–liquid flows. The two-fluid model had momentum and mass balances for each phase. The equations comprising this model were solved numerically by applying a two-step semi-implicit time integration procedure. A finite difference numerical scheme with a staggered mesh was used. Previously, the resulting linear algebraic equations were solved by a Gaussian band solver.” They then used Krylov subspace methods (GMRES, Bi-CGSTAB) to solve the linear algebraic equations and noticed significant improvement in computation time versus the direct Gaussian method. Source: Nordsveen M, Moe R. "Preconditioned Krylov Subspace Methods Used in Solving Two-Dimensional Transient Two-Phase Flows.”

20. Application of Krylov Subspace Method • “The computational time using these iterative methods was from 2 to 6 times shorter than when using the Gaussian band solver. The speed-up in computational time increased with the number of unknowns.” Source: Nordsveen M, Moe R. "Preconditioned Krylov Subspace Methods Used in Solving Two-Dimensional Transient Two-Phase Flows.”