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

• 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