9.4 Krylov subspace methods

Krylov subspace methods are powerful tools for solving large, sparse linear systems. They work by iteratively building up a solution space, using clever math tricks to make each step count. These methods are a game-changer for tackling problems that would be too big or slow to solve directly.

In this section, we'll dive into the two main Krylov methods: Conjugate Gradient (CG) for symmetric matrices and Generalized Minimal Residual (GMRES) for non-symmetric ones. We'll explore how they work, their pros and cons, and when to use each one.

Krylov Subspace Methods

Fundamentals and Derivation

Top images from around the web for Fundamentals and Derivation

• 

  Submitted: Inexact Krylov iterations and relaxation strategies with fast-multipole boundary ... View original

  Is this image relevant?

• 

  Submitted: Inexact Krylov iterations and relaxation strategies with fast-multipole boundary ... View original

  Is this image relevant?

1 of 1

Top images from around the web for Fundamentals and Derivation

• 

  Submitted: Inexact Krylov iterations and relaxation strategies with fast-multipole boundary ... View original

  Is this image relevant?

• 

  Submitted: Inexact Krylov iterations and relaxation strategies with fast-multipole boundary ... View original

  Is this image relevant?

1 of 1

• Krylov subspace methods solve large, sparse linear systems Ax = b iteratively, where A represents an n × n matrix and b a vector of length n
• Define Krylov subspace as $K_m(A, v) = span\{v, Av, A^2v, ..., A^{m-1}v\}$, typically using initial residual $r_0 = b - Ax_0$ as v
• Approximate solution x within Krylov subspace grows in dimension each iteration
• Employ Arnoldi iteration to generate orthonormal basis for Krylov subspace
• Classify into two main categories
  • Lanczos process-based methods for symmetric matrices
  • Arnoldi process-based methods for non-symmetric matrices
• Minimize residual norm $||b - Ax_k||$ over Krylov subspace as fundamental principle
• Apply preconditioning techniques to improve convergence by transforming original system
  • Examples of preconditioning methods (Jacobi, Gauss-Seidel, incomplete LU factorization)

Advanced Concepts and Techniques

• Utilize matrix-vector products as primary operations
  • Suitable for large, sparse systems where A storage occurs implicitly
• Implement restarting for GMRES (GMRES(m)) to limit memory requirements
  • Example: Restart every 20 iterations (GMRES(20))
• Employ preconditioning for both CG and GMRES to enhance convergence
  • Left preconditioning: Solve $M^{-1}Ax = M^{-1}b$
  • Right preconditioning: Solve $AM^{-1}y = b$, then $x = M^{-1}y$
• Consider hybrid methods combining multiple Krylov techniques
  • Example: GMRES-DR (Deflated Restarting) for challenging eigenvalue distributions

CG and GMRES for Linear Systems

Conjugate Gradient Method

• Design CG for symmetric positive definite matrices
• Minimize A-norm of error at each iteration: $||x - x_k||_A = \sqrt{(x - x_k)^T A (x - x_k)}$
• Generate sequence of conjugate (A-orthogonal) vectors
  • Two vectors u and v are A-orthogonal if $u^T A v = 0$
• Update solution and residual vectors using simple recurrence relations
  • $x_{k+1} = x_k + \alpha_k p_k$
  • $r_{k+1} = r_k - \alpha_k A p_k$
• Calculate step size $\alpha_k$ and update direction $p_{k+1}$ using dot products
  • $\alpha_k = \frac{r_k^T r_k}{p_k^T A p_k}$
  • $\beta_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$
  • $p_{k+1} = r_{k+1} + \beta_k p_k$

Generalized Minimal Residual Method

• Apply GMRES to general non-symmetric matrices
• Minimize Euclidean norm of residual over Krylov subspace: $\min_{x \in K_m} ||b - Ax||_2$
• Use Arnoldi iteration to build orthonormal basis for Krylov subspace
  • Generate orthonormal vectors $v_1, v_2, ..., v_m$ spanning $K_m(A, r_0)$
• Solve least-squares problem to find optimal solution in subspace
  • Minimize $||c - H_m y||_2$ where $H_m$ upper Hessenberg matrix from Arnoldi process
• Require more storage than CG due to growing orthonormal basis
  • Storage requirement increases linearly with iteration count

Convergence and Complexity of Krylov Methods

Convergence Analysis

• Relate convergence rate to spectral properties of coefficient matrix A
  • Condition number: ratio of largest to smallest singular value
  • Eigenvalue distribution: clustering and spread
• Bound CG convergence using ratio of largest to smallest eigenvalue of A
  • $||x - x_k||_A \leq 2 \left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^k ||x - x_0||_A$
  • $\kappa$ condition number of A
• Analyze GMRES convergence based on eigenvalue clustering and non-normality of A
  • Convergence bounds more complex than CG
  • Non-normal matrices can exhibit superlinear convergence behavior

Computational Complexity

• Calculate CG complexity per iteration
  • Dense matrices: $O(n^2)$ operations
  • Sparse matrices: $O(nnz)$ operations (nnz number of non-zero elements)
• Determine GMRES complexity for m iterations
  • Dense matrices: $O(mn^2)$ operations
  • Sparse matrices: $O(m^2n + m \cdot nnz)$ operations
• Compare storage requirements
  • CG: $O(n)$ storage
  • GMRES: $O(mn)$ storage, grows with iterations
• Balance preconditioning benefits against added computational cost
  • Example: Incomplete Cholesky preconditioning for CG adds $O(n)$ operations per iteration

Implementation and Comparison of Krylov Methods

Implementation Considerations

• Develop efficient routines for core operations
  • Sparse matrix-vector multiplication
  • Vector operations (dot products, vector additions)
  • Orthogonalization procedures (Gram-Schmidt, Householder)
• Choose initial guess $x_0$ to impact convergence rate
  • Zero vector
  • Solution from previous similar system
• Define stopping criteria
  • Relative residual norm: $||r_k|| / ||b|| < \epsilon$
  • Absolute residual norm: $||r_k|| < \epsilon$
• Address numerical stability issues
  • Implement iterative refinement for GMRES
  • Utilize mixed-precision arithmetic for improved accuracy

Performance Comparison and Optimization

• Evaluate Krylov methods using multiple metrics
  • Convergence rate: iterations required for desired accuracy
  • Computational time: wall clock time for solution
  • Memory usage: peak memory consumption during execution
  • Robustness: performance across varied problem types
• Implement parallel versions to reduce computation time
  • Parallelize matrix-vector products
  • Distribute vector operations across multiple processors
• Develop hybrid solvers for challenging problems
  • Combine direct and iterative methods
  • Example: Use direct solver for preconditioner, Krylov method for main system