0
# Krylov Iterative Methods
1

2
Complete suite of Krylov subspace methods for iterative solution of linear systems. These methods are optimized for use with AMG preconditioning and provide the acceleration layer for PyAMG solvers.
3

4
## Capabilities
5

6
### GMRES Family
7

8
Generalized Minimal Residual methods for general (nonsymmetric) linear systems.
9

10
```python { .api }
11
def gmres(A, b, x0=None, tol=1e-5, restart=None, maxiter=None, 
12
          M=None, callback=None, residuals=None, **kwargs):
13
    """
14
    Generalized Minimal Residual method.
15
    
16
    Solves linear system Ax = b using GMRES with optional
17
    restarting and preconditioning. Most robust Krylov method
18
    for general matrices.
19
    
20
    Parameters:
21
    - A: sparse matrix or LinearOperator, coefficient matrix
22
    - b: array-like, right-hand side vector
23
    - x0: array-like, initial guess (default: zero vector)
24
    - tol: float, convergence tolerance (relative residual)
25
    - restart: int, restart parameter (default: no restart)
26
    - maxiter: int, maximum iterations (default: n)
27
    - M: sparse matrix or LinearOperator, preconditioner
28
    - callback: callable, function called after each iteration
29
    - residuals: list, stores residual norms if provided
30
    - orthog: str, orthogonalization method ('mgs', 'householder')
31
    
32
    Returns:
33
    tuple: (x, info) where:
34
        - x: array, solution vector
35
        - info: int, convergence flag:
36
            * 0: successful convergence
37
            * >0: convergence not achieved in maxiter iterations
38
    
39
    Raises:
40
    ValueError: if A and b have incompatible dimensions
41
    """
42

43
def gmres_householder(A, b, x0=None, tol=1e-5, restart=None, 
44
                     maxiter=None, M=None, **kwargs):
45
    """
46
    GMRES with Householder orthogonalization.
47
    
48
    Variant of GMRES using Householder reflections for
49
    orthogonalization, providing better numerical stability
50
    at higher computational cost.
51
    
52
    Parameters: (same as gmres)
53
    
54
    Returns:
55
    tuple: (x, info)
56
    """
57

58
def gmres_mgs(A, b, x0=None, tol=1e-5, restart=None,
59
              maxiter=None, M=None, **kwargs):
60
    """
61
    GMRES with Modified Gram-Schmidt orthogonalization.
62
    
63
    GMRES variant using Modified Gram-Schmidt for 
64
    orthogonalization, balancing stability and efficiency.
65
    
66
    Parameters: (same as gmres)
67
    
68
    Returns:
69
    tuple: (x, info)
70
    """
71

72
def fgmres(A, b, x0=None, tol=1e-5, restart=None, maxiter=None,
73
           M=None, **kwargs):
74
    """
75
    Flexible GMRES method.
76
    
77
    Variant of GMRES that allows variable preconditioning,
78
    useful with inner-outer iteration schemes and adaptive
79
    preconditioning strategies.
80
    
81
    Parameters: (same as gmres)
82
    - M: can vary between iterations for flexible preconditioning
83
    
84
    Returns:
85
    tuple: (x, info)
86
    """
87
```
88

89
**Usage Examples:**
90

91
```python
92
import pyamg
93
import numpy as np
94

95
# Basic GMRES solve
96
A = pyamg.gallery.poisson((100, 100))
97
b = np.random.rand(A.shape[0])
98
x, info = pyamg.krylov.gmres(A, b, tol=1e-8)
99
print(f"GMRES converged: {info == 0}")
100

101
# GMRES with AMG preconditioning
102
ml = pyamg.ruge_stuben_solver(A)
103
M = ml.aspreconditioner()
104
x, info = pyamg.krylov.gmres(A, b, M=M, tol=1e-10)
105

106
# GMRES with restart
107
x, info = pyamg.krylov.gmres(A, b, restart=30, maxiter=100)
108

109
# Monitor convergence
110
residuals = []
111
x, info = pyamg.krylov.gmres(A, b, residuals=residuals)
112
print(f"Residuals: {len(residuals)} iterations")
113

114
# Flexible GMRES for variable preconditioning
115
x, info = pyamg.krylov.fgmres(A, b, M=M)
116
```
117

118
### Conjugate Gradient Family
119

120
Methods for symmetric positive definite systems with optimal convergence properties.
121

122
```python { .api }
123
def cg(A, b, x0=None, tol=1e-5, maxiter=None, M=None, 
124
       callback=None, residuals=None, **kwargs):
125
    """
126
    Conjugate Gradient method.
127
    
128
    Optimal Krylov method for symmetric positive definite
129
    systems. Requires fewer operations per iteration than
130
    GMRES and has theoretical convergence guarantees.
131
    
132
    Parameters:
133
    - A: sparse matrix or LinearOperator (must be SPD)
134
    - b: array-like, right-hand side vector
135
    - x0: array-like, initial guess
136
    - tol: float, convergence tolerance  
137
    - maxiter: int, maximum iterations (default: n)
138
    - M: sparse matrix or LinearOperator, SPD preconditioner
139
    - callback: callable, iteration callback function
140
    - residuals: list, stores residual norms
141
    
142
    Returns:
143
    tuple: (x, info)
144
    
145
    Raises:
146
    ValueError: if A is not SPD or dimensions are incompatible
147
    """
148

149
def cr(A, b, x0=None, tol=1e-5, maxiter=None, M=None, **kwargs):
150
    """
151
    Conjugate Residual method.
152
    
153
    Variant of CG that minimizes A-norm of residual rather
154
    than Euclidean norm. Can be more robust for indefinite
155
    systems but requires A to be symmetric.
156
    
157
    Parameters: (same as cg, but A need only be symmetric)
158
    
159
    Returns:
160
    tuple: (x, info)
161
    """
162

163
def cgnr(A, b, x0=None, tol=1e-5, maxiter=None, M=None, **kwargs):
164
    """
165
    Conjugate Gradient Normal Residual method.
166
    
167
    Applies CG to normal equations A^T A x = A^T b.
168
    Works for rectangular or singular systems but can
169
    have poor conditioning.
170
    
171
    Parameters:
172
    - A: sparse matrix (can be rectangular or singular)
173
    - b: array-like, right-hand side
174
    - M: preconditioner for A^T A system
175
    
176
    Returns:
177
    tuple: (x, info)
178
    """
179

180
def cgne(A, b, x0=None, tol=1e-5, maxiter=None, M=None, **kwargs):
181
    """
182
    Conjugate Gradient Normal Equation method.
183
    
184
    Applies CG to normal equations A A^T y = b, x = A^T y.
185
    Alternative formulation for least squares problems.
186
    
187
    Parameters: (same as cgnr)
188
    
189
    Returns:
190
    tuple: (x, info)
191
    """
192
```
193

194
**Usage Examples:**
195

196
```python
197
# CG for SPD system
198
A = pyamg.gallery.poisson((80, 80))
199
b = np.random.rand(A.shape[0])
200
x, info = pyamg.krylov.cg(A, b, tol=1e-10)
201

202
# CG with AMG preconditioning (most efficient combination)
203
ml = pyamg.smoothed_aggregation_solver(A)
204
M = ml.aspreconditioner() 
205
x, info = pyamg.krylov.cg(A, b, M=M)
206

207
# Conjugate residual for symmetric indefinite
208
A_indef = pyamg.gallery.gauge_laplacian((50, 50))  # Singular
209
x, info = pyamg.krylov.cr(A_indef, b)
210

211
# Normal equations for rectangular system
212
A_rect = np.random.rand(100, 80)
213
b_rect = np.random.rand(100)
214
x, info = pyamg.krylov.cgnr(A_rect, b_rect)
215
```
216

217
### BiConjugate Methods
218

219
Methods for nonsymmetric systems using bi-orthogonalization.
220

221
```python { .api }
222
def bicgstab(A, b, x0=None, tol=1e-5, maxiter=None, M=None,
223
             callback=None, **kwargs):
224
    """
225
    Biconjugate Gradient Stabilized method.
226
    
227
    Stabilized variant of BiCG that avoids erratic convergence
228
    behavior. Good general-purpose method for nonsymmetric
229
    systems when GMRES memory requirements are prohibitive.
230
    
231
    Parameters:
232
    - A: sparse matrix or LinearOperator (general matrix)
233
    - b: array-like, right-hand side vector
234
    - x0: array-like, initial guess
235
    - tol: float, convergence tolerance
236
    - maxiter: int, maximum iterations  
237
    - M: sparse matrix or LinearOperator, preconditioner
238
    - callback: callable, iteration callback
239
    
240
    Returns:
241
    tuple: (x, info)
242
    
243
    Notes:
244
    - Memory requirements: O(n) vs O(kn) for GMRES
245
    - Can experience irregular convergence behavior
246
    - May break down for certain matrices
247
    """
248
```
249

250
**Usage Examples:**
251

252
```python
253
# BiCGSTAB for nonsymmetric system
254
A = pyamg.gallery.poisson((60, 60)) + 0.1 * sp.random(3600, 3600, density=0.01)
255
b = np.random.rand(A.shape[0])
256
x, info = pyamg.krylov.bicgstab(A, b, tol=1e-8)
257

258
# With preconditioning
259
ml = pyamg.ruge_stuben_solver(A)
260
M = ml.aspreconditioner()
261
x, info = pyamg.krylov.bicgstab(A, b, M=M)
262
```
263

264
### Other Iterative Methods
265

266
Additional Krylov and related methods for special cases.
267

268
```python { .api }
269
def steepest_descent(A, b, x0=None, tol=1e-5, maxiter=None, **kwargs):
270
    """
271
    Steepest descent method.
272
    
273
    Simple gradient descent method, mainly for educational
274
    purposes. Very slow convergence but robust and simple.
275
    
276
    Parameters:
277
    - A: sparse matrix (should be SPD for convergence)
278
    - b: array-like, right-hand side
279
    - x0: array-like, initial guess
280
    - tol: float, convergence tolerance
281
    - maxiter: int, maximum iterations
282
    
283
    Returns:
284
    tuple: (x, info)
285
    
286
    Notes:
287
    - Convergence rate depends on condition number
288
    - Not recommended for practical use
289
    - Useful for teaching and algorithm comparison
290
    """
291

292
def minimal_residual(A, b, x0=None, tol=1e-5, maxiter=None, 
293
                    M=None, **kwargs):
294
    """
295
    Minimal Residual method.
296
    
297
    Minimizes residual norm at each iteration without
298
    orthogonality constraints. Can be useful for certain
299
    problem classes.
300
    
301
    Parameters: (same as other Krylov methods)
302
    
303
    Returns:
304
    tuple: (x, info)
305
    """
306
```
307

308
**Usage Examples:**
309

310
```python
311
# Steepest descent (educational)
312
A = pyamg.gallery.poisson((30, 30))
313
b = np.random.rand(A.shape[0])
314
x, info = pyamg.krylov.steepest_descent(A, b, maxiter=1000)
315

316
# Minimal residual
317
x, info = pyamg.krylov.minimal_residual(A, b)
318
```
319

320
## Method Selection Guidelines
321

322
### Problem Type Recommendations
323

324
- **SPD Systems**: Conjugate Gradient with AMG preconditioning
325
- **General Symmetric**: Conjugate Residual or GMRES
326
- **Nonsymmetric**: GMRES (small-medium) or BiCGSTAB (large)
327
- **Indefinite**: GMRES with appropriate preconditioning
328
- **Singular/Rectangular**: CGNR or CGNE
329

330
### Preconditioning Strategy
331

332
```python
333
# Optimal combination: CG + SA AMG for SPD
334
A = pyamg.gallery.linear_elasticity((40, 40))
335
ml = pyamg.smoothed_aggregation_solver(A)
336
x, info = pyamg.krylov.cg(A, b, M=ml.aspreconditioner())
337

338
# General problems: GMRES + Classical AMG  
339
A = pyamg.gallery.poisson((50, 50))
340
ml = pyamg.ruge_stuben_solver(A)
341
x, info = pyamg.krylov.gmres(A, b, M=ml.aspreconditioner())
342
```
343

344
### Performance Considerations
345

346
- **Memory**: CG < BiCGSTAB < GMRES(restart) < GMRES(no restart)
347
- **Robustness**: GMRES > BiCGSTAB > CG
348
- **Speed per iteration**: CG > BiCGSTAB > GMRES
349
- **Preconditioning**: Essential for practical performance
350

351
### Convergence Monitoring
352

353
```python
354
def monitor_convergence(x):
355
    residual = np.linalg.norm(b - A*x)
356
    print(f"Iteration residual: {residual:.2e}")
357

358
x, info = pyamg.krylov.gmres(A, b, callback=monitor_convergence)
359
```
360

361
### Common Issues
362

363
- **Stagnation**: Try different method or better preconditioner
364
- **Breakdown**: Switch from BiCGSTAB to GMRES  
365
- **Memory**: Use restarted GMRES or BiCGSTAB
366
- **Slow convergence**: Improve preconditioning or initial guess