Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_isolve/_gcrotmk.py: 5%
197 statements
« prev ^ index » next coverage.py v7.4.1, created at 2024-02-14 06:37 +0000
1# Copyright (C) 2015, Pauli Virtanen <pav@iki.fi>
2# Distributed under the same license as SciPy.
4import numpy as np
5from numpy.linalg import LinAlgError
6from scipy.linalg import (get_blas_funcs, qr, solve, svd, qr_insert, lstsq)
7from .iterative import _get_atol_rtol
8from scipy.sparse.linalg._isolve.utils import make_system
9from scipy._lib.deprecation import _NoValue, _deprecate_positional_args
12__all__ = ['gcrotmk']
15def _fgmres(matvec, v0, m, atol, lpsolve=None, rpsolve=None, cs=(), outer_v=(),
16 prepend_outer_v=False):
17 """
18 FGMRES Arnoldi process, with optional projection or augmentation
20 Parameters
21 ----------
22 matvec : callable
23 Operation A*x
24 v0 : ndarray
25 Initial vector, normalized to nrm2(v0) == 1
26 m : int
27 Number of GMRES rounds
28 atol : float
29 Absolute tolerance for early exit
30 lpsolve : callable
31 Left preconditioner L
32 rpsolve : callable
33 Right preconditioner R
34 cs : list of (ndarray, ndarray)
35 Columns of matrices C and U in GCROT
36 outer_v : list of ndarrays
37 Augmentation vectors in LGMRES
38 prepend_outer_v : bool, optional
39 Whether augmentation vectors come before or after
40 Krylov iterates
42 Raises
43 ------
44 LinAlgError
45 If nans encountered
47 Returns
48 -------
49 Q, R : ndarray
50 QR decomposition of the upper Hessenberg H=QR
51 B : ndarray
52 Projections corresponding to matrix C
53 vs : list of ndarray
54 Columns of matrix V
55 zs : list of ndarray
56 Columns of matrix Z
57 y : ndarray
58 Solution to ||H y - e_1||_2 = min!
59 res : float
60 The final (preconditioned) residual norm
62 """
64 if lpsolve is None:
65 def lpsolve(x):
66 return x
67 if rpsolve is None:
68 def rpsolve(x):
69 return x
71 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'], (v0,))
73 vs = [v0]
74 zs = []
75 y = None
76 res = np.nan
78 m = m + len(outer_v)
80 # Orthogonal projection coefficients
81 B = np.zeros((len(cs), m), dtype=v0.dtype)
83 # H is stored in QR factorized form
84 Q = np.ones((1, 1), dtype=v0.dtype)
85 R = np.zeros((1, 0), dtype=v0.dtype)
87 eps = np.finfo(v0.dtype).eps
89 breakdown = False
91 # FGMRES Arnoldi process
92 for j in range(m):
93 # L A Z = C B + V H
95 if prepend_outer_v and j < len(outer_v):
96 z, w = outer_v[j]
97 elif prepend_outer_v and j == len(outer_v):
98 z = rpsolve(v0)
99 w = None
100 elif not prepend_outer_v and j >= m - len(outer_v):
101 z, w = outer_v[j - (m - len(outer_v))]
102 else:
103 z = rpsolve(vs[-1])
104 w = None
106 if w is None:
107 w = lpsolve(matvec(z))
108 else:
109 # w is clobbered below
110 w = w.copy()
112 w_norm = nrm2(w)
114 # GCROT projection: L A -> (1 - C C^H) L A
115 # i.e. orthogonalize against C
116 for i, c in enumerate(cs):
117 alpha = dot(c, w)
118 B[i,j] = alpha
119 w = axpy(c, w, c.shape[0], -alpha) # w -= alpha*c
121 # Orthogonalize against V
122 hcur = np.zeros(j+2, dtype=Q.dtype)
123 for i, v in enumerate(vs):
124 alpha = dot(v, w)
125 hcur[i] = alpha
126 w = axpy(v, w, v.shape[0], -alpha) # w -= alpha*v
127 hcur[i+1] = nrm2(w)
129 with np.errstate(over='ignore', divide='ignore'):
130 # Careful with denormals
131 alpha = 1/hcur[-1]
133 if np.isfinite(alpha):
134 w = scal(alpha, w)
136 if not (hcur[-1] > eps * w_norm):
137 # w essentially in the span of previous vectors,
138 # or we have nans. Bail out after updating the QR
139 # solution.
140 breakdown = True
142 vs.append(w)
143 zs.append(z)
145 # Arnoldi LSQ problem
147 # Add new column to H=Q@R, padding other columns with zeros
148 Q2 = np.zeros((j+2, j+2), dtype=Q.dtype, order='F')
149 Q2[:j+1,:j+1] = Q
150 Q2[j+1,j+1] = 1
152 R2 = np.zeros((j+2, j), dtype=R.dtype, order='F')
153 R2[:j+1,:] = R
155 Q, R = qr_insert(Q2, R2, hcur, j, which='col',
156 overwrite_qru=True, check_finite=False)
158 # Transformed least squares problem
159 # || Q R y - inner_res_0 * e_1 ||_2 = min!
160 # Since R = [R'; 0], solution is y = inner_res_0 (R')^{-1} (Q^H)[:j,0]
162 # Residual is immediately known
163 res = abs(Q[0,-1])
165 # Check for termination
166 if res < atol or breakdown:
167 break
169 if not np.isfinite(R[j,j]):
170 # nans encountered, bail out
171 raise LinAlgError()
173 # -- Get the LSQ problem solution
175 # The problem is triangular, but the condition number may be
176 # bad (or in case of breakdown the last diagonal entry may be
177 # zero), so use lstsq instead of trtrs.
178 y, _, _, _, = lstsq(R[:j+1,:j+1], Q[0,:j+1].conj())
180 B = B[:,:j+1]
182 return Q, R, B, vs, zs, y, res
185@_deprecate_positional_args(version="1.14.0")
186def gcrotmk(A, b, x0=None, *, tol=_NoValue, maxiter=1000, M=None, callback=None,
187 m=20, k=None, CU=None, discard_C=False, truncate='oldest',
188 atol=None, rtol=1e-5):
189 """
190 Solve a matrix equation using flexible GCROT(m,k) algorithm.
192 Parameters
193 ----------
194 A : {sparse matrix, ndarray, LinearOperator}
195 The real or complex N-by-N matrix of the linear system.
196 Alternatively, ``A`` can be a linear operator which can
197 produce ``Ax`` using, e.g.,
198 ``scipy.sparse.linalg.LinearOperator``.
199 b : ndarray
200 Right hand side of the linear system. Has shape (N,) or (N,1).
201 x0 : ndarray
202 Starting guess for the solution.
203 rtol, atol : float, optional
204 Parameters for the convergence test. For convergence,
205 ``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
206 The default is ``rtol=1e-5``, the default for ``atol`` is ``rtol``.
208 .. warning::
210 The default value for ``atol`` will be changed to ``0.0`` in
211 SciPy 1.14.0.
212 maxiter : int, optional
213 Maximum number of iterations. Iteration will stop after maxiter
214 steps even if the specified tolerance has not been achieved.
215 M : {sparse matrix, ndarray, LinearOperator}, optional
216 Preconditioner for A. The preconditioner should approximate the
217 inverse of A. gcrotmk is a 'flexible' algorithm and the preconditioner
218 can vary from iteration to iteration. Effective preconditioning
219 dramatically improves the rate of convergence, which implies that
220 fewer iterations are needed to reach a given error tolerance.
221 callback : function, optional
222 User-supplied function to call after each iteration. It is called
223 as callback(xk), where xk is the current solution vector.
224 m : int, optional
225 Number of inner FGMRES iterations per each outer iteration.
226 Default: 20
227 k : int, optional
228 Number of vectors to carry between inner FGMRES iterations.
229 According to [2]_, good values are around m.
230 Default: m
231 CU : list of tuples, optional
232 List of tuples ``(c, u)`` which contain the columns of the matrices
233 C and U in the GCROT(m,k) algorithm. For details, see [2]_.
234 The list given and vectors contained in it are modified in-place.
235 If not given, start from empty matrices. The ``c`` elements in the
236 tuples can be ``None``, in which case the vectors are recomputed
237 via ``c = A u`` on start and orthogonalized as described in [3]_.
238 discard_C : bool, optional
239 Discard the C-vectors at the end. Useful if recycling Krylov subspaces
240 for different linear systems.
241 truncate : {'oldest', 'smallest'}, optional
242 Truncation scheme to use. Drop: oldest vectors, or vectors with
243 smallest singular values using the scheme discussed in [1,2].
244 See [2]_ for detailed comparison.
245 Default: 'oldest'
246 tol : float, optional, deprecated
248 .. deprecated:: 1.12.0
249 `gcrotmk` keyword argument ``tol`` is deprecated in favor of
250 ``rtol`` and will be removed in SciPy 1.14.0.
252 Returns
253 -------
254 x : ndarray
255 The solution found.
256 info : int
257 Provides convergence information:
259 * 0 : successful exit
260 * >0 : convergence to tolerance not achieved, number of iterations
262 Examples
263 --------
264 >>> import numpy as np
265 >>> from scipy.sparse import csc_matrix
266 >>> from scipy.sparse.linalg import gcrotmk
267 >>> R = np.random.randn(5, 5)
268 >>> A = csc_matrix(R)
269 >>> b = np.random.randn(5)
270 >>> x, exit_code = gcrotmk(A, b, atol=1e-5)
271 >>> print(exit_code)
272 0
273 >>> np.allclose(A.dot(x), b)
274 True
276 References
277 ----------
278 .. [1] E. de Sturler, ''Truncation strategies for optimal Krylov subspace
279 methods'', SIAM J. Numer. Anal. 36, 864 (1999).
280 .. [2] J.E. Hicken and D.W. Zingg, ''A simplified and flexible variant
281 of GCROT for solving nonsymmetric linear systems'',
282 SIAM J. Sci. Comput. 32, 172 (2010).
283 .. [3] M.L. Parks, E. de Sturler, G. Mackey, D.D. Johnson, S. Maiti,
284 ''Recycling Krylov subspaces for sequences of linear systems'',
285 SIAM J. Sci. Comput. 28, 1651 (2006).
287 """
288 A,M,x,b,postprocess = make_system(A,M,x0,b)
290 if not np.isfinite(b).all():
291 raise ValueError("RHS must contain only finite numbers")
293 if truncate not in ('oldest', 'smallest'):
294 raise ValueError(f"Invalid value for 'truncate': {truncate!r}")
296 matvec = A.matvec
297 psolve = M.matvec
299 if CU is None:
300 CU = []
302 if k is None:
303 k = m
305 axpy, dot, scal = None, None, None
307 if x0 is None:
308 r = b.copy()
309 else:
310 r = b - matvec(x)
312 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'], (x, r))
314 b_norm = nrm2(b)
316 # we call this to get the right atol/rtol and raise warnings as necessary
317 atol, rtol = _get_atol_rtol('gcrotmk', b_norm, tol, atol, rtol)
319 if b_norm == 0:
320 x = b
321 return (postprocess(x), 0)
323 if discard_C:
324 CU[:] = [(None, u) for c, u in CU]
326 # Reorthogonalize old vectors
327 if CU:
328 # Sort already existing vectors to the front
329 CU.sort(key=lambda cu: cu[0] is not None)
331 # Fill-in missing ones
332 C = np.empty((A.shape[0], len(CU)), dtype=r.dtype, order='F')
333 us = []
334 j = 0
335 while CU:
336 # More memory-efficient: throw away old vectors as we go
337 c, u = CU.pop(0)
338 if c is None:
339 c = matvec(u)
340 C[:,j] = c
341 j += 1
342 us.append(u)
344 # Orthogonalize
345 Q, R, P = qr(C, overwrite_a=True, mode='economic', pivoting=True)
346 del C
348 # C := Q
349 cs = list(Q.T)
351 # U := U P R^-1, back-substitution
352 new_us = []
353 for j in range(len(cs)):
354 u = us[P[j]]
355 for i in range(j):
356 u = axpy(us[P[i]], u, u.shape[0], -R[i,j])
357 if abs(R[j,j]) < 1e-12 * abs(R[0,0]):
358 # discard rest of the vectors
359 break
360 u = scal(1.0/R[j,j], u)
361 new_us.append(u)
363 # Form the new CU lists
364 CU[:] = list(zip(cs, new_us))[::-1]
366 if CU:
367 axpy, dot = get_blas_funcs(['axpy', 'dot'], (r,))
369 # Solve first the projection operation with respect to the CU
370 # vectors. This corresponds to modifying the initial guess to
371 # be
372 #
373 # x' = x + U y
374 # y = argmin_y || b - A (x + U y) ||^2
375 #
376 # The solution is y = C^H (b - A x)
377 for c, u in CU:
378 yc = dot(c, r)
379 x = axpy(u, x, x.shape[0], yc)
380 r = axpy(c, r, r.shape[0], -yc)
382 # GCROT main iteration
383 for j_outer in range(maxiter):
384 # -- callback
385 if callback is not None:
386 callback(x)
388 beta = nrm2(r)
390 # -- check stopping condition
391 beta_tol = max(atol, rtol * b_norm)
393 if beta <= beta_tol and (j_outer > 0 or CU):
394 # recompute residual to avoid rounding error
395 r = b - matvec(x)
396 beta = nrm2(r)
398 if beta <= beta_tol:
399 j_outer = -1
400 break
402 ml = m + max(k - len(CU), 0)
404 cs = [c for c, u in CU]
406 try:
407 Q, R, B, vs, zs, y, pres = _fgmres(matvec,
408 r/beta,
409 ml,
410 rpsolve=psolve,
411 atol=max(atol, rtol*b_norm)/beta,
412 cs=cs)
413 y *= beta
414 except LinAlgError:
415 # Floating point over/underflow, non-finite result from
416 # matmul etc. -- report failure.
417 break
419 #
420 # At this point,
421 #
422 # [A U, A Z] = [C, V] G; G = [ I B ]
423 # [ 0 H ]
424 #
425 # where [C, V] has orthonormal columns, and r = beta v_0. Moreover,
426 #
427 # || b - A (x + Z y + U q) ||_2 = || r - C B y - V H y - C q ||_2 = min!
428 #
429 # from which y = argmin_y || beta e_1 - H y ||_2, and q = -B y
430 #
432 #
433 # GCROT(m,k) update
434 #
436 # Define new outer vectors
438 # ux := (Z - U B) y
439 ux = zs[0]*y[0]
440 for z, yc in zip(zs[1:], y[1:]):
441 ux = axpy(z, ux, ux.shape[0], yc) # ux += z*yc
442 by = B.dot(y)
443 for cu, byc in zip(CU, by):
444 c, u = cu
445 ux = axpy(u, ux, ux.shape[0], -byc) # ux -= u*byc
447 # cx := V H y
448 hy = Q.dot(R.dot(y))
449 cx = vs[0] * hy[0]
450 for v, hyc in zip(vs[1:], hy[1:]):
451 cx = axpy(v, cx, cx.shape[0], hyc) # cx += v*hyc
453 # Normalize cx, maintaining cx = A ux
454 # This new cx is orthogonal to the previous C, by construction
455 try:
456 alpha = 1/nrm2(cx)
457 if not np.isfinite(alpha):
458 raise FloatingPointError()
459 except (FloatingPointError, ZeroDivisionError):
460 # Cannot update, so skip it
461 continue
463 cx = scal(alpha, cx)
464 ux = scal(alpha, ux)
466 # Update residual and solution
467 gamma = dot(cx, r)
468 r = axpy(cx, r, r.shape[0], -gamma) # r -= gamma*cx
469 x = axpy(ux, x, x.shape[0], gamma) # x += gamma*ux
471 # Truncate CU
472 if truncate == 'oldest':
473 while len(CU) >= k and CU:
474 del CU[0]
475 elif truncate == 'smallest':
476 if len(CU) >= k and CU:
477 # cf. [1,2]
478 D = solve(R[:-1,:].T, B.T).T
479 W, sigma, V = svd(D)
481 # C := C W[:,:k-1], U := U W[:,:k-1]
482 new_CU = []
483 for j, w in enumerate(W[:,:k-1].T):
484 c, u = CU[0]
485 c = c * w[0]
486 u = u * w[0]
487 for cup, wp in zip(CU[1:], w[1:]):
488 cp, up = cup
489 c = axpy(cp, c, c.shape[0], wp)
490 u = axpy(up, u, u.shape[0], wp)
492 # Reorthogonalize at the same time; not necessary
493 # in exact arithmetic, but floating point error
494 # tends to accumulate here
495 for cp, up in new_CU:
496 alpha = dot(cp, c)
497 c = axpy(cp, c, c.shape[0], -alpha)
498 u = axpy(up, u, u.shape[0], -alpha)
499 alpha = nrm2(c)
500 c = scal(1.0/alpha, c)
501 u = scal(1.0/alpha, u)
503 new_CU.append((c, u))
504 CU[:] = new_CU
506 # Add new vector to CU
507 CU.append((cx, ux))
509 # Include the solution vector to the span
510 CU.append((None, x.copy()))
511 if discard_C:
512 CU[:] = [(None, uz) for cz, uz in CU]
514 return postprocess(x), j_outer + 1