Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/sparse/linalg/isolve/_gcrotmk.py : 4%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Copyright (C) 2015, Pauli Virtanen <pav@iki.fi>
2# Distributed under the same license as SciPy.
4import warnings
5import numpy as np
6from numpy.linalg import LinAlgError
7from scipy.linalg import (get_blas_funcs, qr, solve, svd, qr_insert, lstsq)
8from scipy.sparse.linalg.isolve.utils import make_system
11__all__ = ['gcrotmk']
14def _fgmres(matvec, v0, m, atol, lpsolve=None, rpsolve=None, cs=(), outer_v=(),
15 prepend_outer_v=False):
16 """
17 FGMRES Arnoldi process, with optional projection or augmentation
19 Parameters
20 ----------
21 matvec : callable
22 Operation A*x
23 v0 : ndarray
24 Initial vector, normalized to nrm2(v0) == 1
25 m : int
26 Number of GMRES rounds
27 atol : float
28 Absolute tolerance for early exit
29 lpsolve : callable
30 Left preconditioner L
31 rpsolve : callable
32 Right preconditioner R
33 CU : list of (ndarray, ndarray)
34 Columns of matrices C and U in GCROT
35 outer_v : list of ndarrays
36 Augmentation vectors in LGMRES
37 prepend_outer_v : bool, optional
38 Whether augmentation vectors come before or after
39 Krylov iterates
41 Raises
42 ------
43 LinAlgError
44 If nans encountered
46 Returns
47 -------
48 Q, R : ndarray
49 QR decomposition of the upper Hessenberg H=QR
50 B : ndarray
51 Projections corresponding to matrix C
52 vs : list of ndarray
53 Columns of matrix V
54 zs : list of ndarray
55 Columns of matrix Z
56 y : ndarray
57 Solution to ||H y - e_1||_2 = min!
58 res : float
59 The final (preconditioned) residual norm
61 """
63 if lpsolve is None:
64 lpsolve = lambda x: x
65 if rpsolve is None:
66 rpsolve = lambda x: x
68 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'], (v0,))
70 vs = [v0]
71 zs = []
72 y = None
73 res = np.nan
75 m = m + len(outer_v)
77 # Orthogonal projection coefficients
78 B = np.zeros((len(cs), m), dtype=v0.dtype)
80 # H is stored in QR factorized form
81 Q = np.ones((1, 1), dtype=v0.dtype)
82 R = np.zeros((1, 0), dtype=v0.dtype)
84 eps = np.finfo(v0.dtype).eps
86 breakdown = False
88 # FGMRES Arnoldi process
89 for j in range(m):
90 # L A Z = C B + V H
92 if prepend_outer_v and j < len(outer_v):
93 z, w = outer_v[j]
94 elif prepend_outer_v and j == len(outer_v):
95 z = rpsolve(v0)
96 w = None
97 elif not prepend_outer_v and j >= m - len(outer_v):
98 z, w = outer_v[j - (m - len(outer_v))]
99 else:
100 z = rpsolve(vs[-1])
101 w = None
103 if w is None:
104 w = lpsolve(matvec(z))
105 else:
106 # w is clobbered below
107 w = w.copy()
109 w_norm = nrm2(w)
111 # GCROT projection: L A -> (1 - C C^H) L A
112 # i.e. orthogonalize against C
113 for i, c in enumerate(cs):
114 alpha = dot(c, w)
115 B[i,j] = alpha
116 w = axpy(c, w, c.shape[0], -alpha) # w -= alpha*c
118 # Orthogonalize against V
119 hcur = np.zeros(j+2, dtype=Q.dtype)
120 for i, v in enumerate(vs):
121 alpha = dot(v, w)
122 hcur[i] = alpha
123 w = axpy(v, w, v.shape[0], -alpha) # w -= alpha*v
124 hcur[i+1] = nrm2(w)
126 with np.errstate(over='ignore', divide='ignore'):
127 # Careful with denormals
128 alpha = 1/hcur[-1]
130 if np.isfinite(alpha):
131 w = scal(alpha, w)
133 if not (hcur[-1] > eps * w_norm):
134 # w essentially in the span of previous vectors,
135 # or we have nans. Bail out after updating the QR
136 # solution.
137 breakdown = True
139 vs.append(w)
140 zs.append(z)
142 # Arnoldi LSQ problem
144 # Add new column to H=Q*R, padding other columns with zeros
145 Q2 = np.zeros((j+2, j+2), dtype=Q.dtype, order='F')
146 Q2[:j+1,:j+1] = Q
147 Q2[j+1,j+1] = 1
149 R2 = np.zeros((j+2, j), dtype=R.dtype, order='F')
150 R2[:j+1,:] = R
152 Q, R = qr_insert(Q2, R2, hcur, j, which='col',
153 overwrite_qru=True, check_finite=False)
155 # Transformed least squares problem
156 # || Q R y - inner_res_0 * e_1 ||_2 = min!
157 # Since R = [R'; 0], solution is y = inner_res_0 (R')^{-1} (Q^H)[:j,0]
159 # Residual is immediately known
160 res = abs(Q[0,-1])
162 # Check for termination
163 if res < atol or breakdown:
164 break
166 if not np.isfinite(R[j,j]):
167 # nans encountered, bail out
168 raise LinAlgError()
170 # -- Get the LSQ problem solution
172 # The problem is triangular, but the condition number may be
173 # bad (or in case of breakdown the last diagonal entry may be
174 # zero), so use lstsq instead of trtrs.
175 y, _, _, _, = lstsq(R[:j+1,:j+1], Q[0,:j+1].conj())
177 B = B[:,:j+1]
179 return Q, R, B, vs, zs, y, res
182def gcrotmk(A, b, x0=None, tol=1e-5, maxiter=1000, M=None, callback=None,
183 m=20, k=None, CU=None, discard_C=False, truncate='oldest',
184 atol=None):
185 """
186 Solve a matrix equation using flexible GCROT(m,k) algorithm.
188 Parameters
189 ----------
190 A : {sparse matrix, dense matrix, LinearOperator}
191 The real or complex N-by-N matrix of the linear system.
192 Alternatively, ``A`` can be a linear operator which can
193 produce ``Ax`` using, e.g.,
194 ``scipy.sparse.linalg.LinearOperator``.
195 b : {array, matrix}
196 Right hand side of the linear system. Has shape (N,) or (N,1).
197 x0 : {array, matrix}
198 Starting guess for the solution.
199 tol, atol : float, optional
200 Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``.
201 The default for ``atol`` is `tol`.
203 .. warning::
205 The default value for `atol` will be changed in a future release.
206 For future compatibility, specify `atol` explicitly.
207 maxiter : int, optional
208 Maximum number of iterations. Iteration will stop after maxiter
209 steps even if the specified tolerance has not been achieved.
210 M : {sparse matrix, dense matrix, LinearOperator}, optional
211 Preconditioner for A. The preconditioner should approximate the
212 inverse of A. gcrotmk is a 'flexible' algorithm and the preconditioner
213 can vary from iteration to iteration. Effective preconditioning
214 dramatically improves the rate of convergence, which implies that
215 fewer iterations are needed to reach a given error tolerance.
216 callback : function, optional
217 User-supplied function to call after each iteration. It is called
218 as callback(xk), where xk is the current solution vector.
219 m : int, optional
220 Number of inner FGMRES iterations per each outer iteration.
221 Default: 20
222 k : int, optional
223 Number of vectors to carry between inner FGMRES iterations.
224 According to [2]_, good values are around m.
225 Default: m
226 CU : list of tuples, optional
227 List of tuples ``(c, u)`` which contain the columns of the matrices
228 C and U in the GCROT(m,k) algorithm. For details, see [2]_.
229 The list given and vectors contained in it are modified in-place.
230 If not given, start from empty matrices. The ``c`` elements in the
231 tuples can be ``None``, in which case the vectors are recomputed
232 via ``c = A u`` on start and orthogonalized as described in [3]_.
233 discard_C : bool, optional
234 Discard the C-vectors at the end. Useful if recycling Krylov subspaces
235 for different linear systems.
236 truncate : {'oldest', 'smallest'}, optional
237 Truncation scheme to use. Drop: oldest vectors, or vectors with
238 smallest singular values using the scheme discussed in [1,2].
239 See [2]_ for detailed comparison.
240 Default: 'oldest'
242 Returns
243 -------
244 x : array or matrix
245 The solution found.
246 info : int
247 Provides convergence information:
249 * 0 : successful exit
250 * >0 : convergence to tolerance not achieved, number of iterations
252 References
253 ----------
254 .. [1] E. de Sturler, ''Truncation strategies for optimal Krylov subspace
255 methods'', SIAM J. Numer. Anal. 36, 864 (1999).
256 .. [2] J.E. Hicken and D.W. Zingg, ''A simplified and flexible variant
257 of GCROT for solving nonsymmetric linear systems'',
258 SIAM J. Sci. Comput. 32, 172 (2010).
259 .. [3] M.L. Parks, E. de Sturler, G. Mackey, D.D. Johnson, S. Maiti,
260 ''Recycling Krylov subspaces for sequences of linear systems'',
261 SIAM J. Sci. Comput. 28, 1651 (2006).
263 """
264 A,M,x,b,postprocess = make_system(A,M,x0,b)
266 if not np.isfinite(b).all():
267 raise ValueError("RHS must contain only finite numbers")
269 if truncate not in ('oldest', 'smallest'):
270 raise ValueError("Invalid value for 'truncate': %r" % (truncate,))
272 if atol is None:
273 warnings.warn("scipy.sparse.linalg.gcrotmk called without specifying `atol`. "
274 "The default value will change in the future. To preserve "
275 "current behavior, set ``atol=tol``.",
276 category=DeprecationWarning, stacklevel=2)
277 atol = tol
279 matvec = A.matvec
280 psolve = M.matvec
282 if CU is None:
283 CU = []
285 if k is None:
286 k = m
288 axpy, dot, scal = None, None, None
290 r = b - matvec(x)
292 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'], (x, r))
294 b_norm = nrm2(b)
296 if discard_C:
297 CU[:] = [(None, u) for c, u in CU]
299 # Reorthogonalize old vectors
300 if CU:
301 # Sort already existing vectors to the front
302 CU.sort(key=lambda cu: cu[0] is not None)
304 # Fill-in missing ones
305 C = np.empty((A.shape[0], len(CU)), dtype=r.dtype, order='F')
306 us = []
307 j = 0
308 while CU:
309 # More memory-efficient: throw away old vectors as we go
310 c, u = CU.pop(0)
311 if c is None:
312 c = matvec(u)
313 C[:,j] = c
314 j += 1
315 us.append(u)
317 # Orthogonalize
318 Q, R, P = qr(C, overwrite_a=True, mode='economic', pivoting=True)
319 del C
321 # C := Q
322 cs = list(Q.T)
324 # U := U P R^-1, back-substitution
325 new_us = []
326 for j in range(len(cs)):
327 u = us[P[j]]
328 for i in range(j):
329 u = axpy(us[P[i]], u, u.shape[0], -R[i,j])
330 if abs(R[j,j]) < 1e-12 * abs(R[0,0]):
331 # discard rest of the vectors
332 break
333 u = scal(1.0/R[j,j], u)
334 new_us.append(u)
336 # Form the new CU lists
337 CU[:] = list(zip(cs, new_us))[::-1]
339 if CU:
340 axpy, dot = get_blas_funcs(['axpy', 'dot'], (r,))
342 # Solve first the projection operation with respect to the CU
343 # vectors. This corresponds to modifying the initial guess to
344 # be
345 #
346 # x' = x + U y
347 # y = argmin_y || b - A (x + U y) ||^2
348 #
349 # The solution is y = C^H (b - A x)
350 for c, u in CU:
351 yc = dot(c, r)
352 x = axpy(u, x, x.shape[0], yc)
353 r = axpy(c, r, r.shape[0], -yc)
355 # GCROT main iteration
356 for j_outer in range(maxiter):
357 # -- callback
358 if callback is not None:
359 callback(x)
361 beta = nrm2(r)
363 # -- check stopping condition
364 beta_tol = max(atol, tol * b_norm)
366 if beta <= beta_tol and (j_outer > 0 or CU):
367 # recompute residual to avoid rounding error
368 r = b - matvec(x)
369 beta = nrm2(r)
371 if beta <= beta_tol:
372 j_outer = -1
373 break
375 ml = m + max(k - len(CU), 0)
377 cs = [c for c, u in CU]
379 try:
380 Q, R, B, vs, zs, y, pres = _fgmres(matvec,
381 r/beta,
382 ml,
383 rpsolve=psolve,
384 atol=max(atol, tol*b_norm)/beta,
385 cs=cs)
386 y *= beta
387 except LinAlgError:
388 # Floating point over/underflow, non-finite result from
389 # matmul etc. -- report failure.
390 break
392 #
393 # At this point,
394 #
395 # [A U, A Z] = [C, V] G; G = [ I B ]
396 # [ 0 H ]
397 #
398 # where [C, V] has orthonormal columns, and r = beta v_0. Moreover,
399 #
400 # || b - A (x + Z y + U q) ||_2 = || r - C B y - V H y - C q ||_2 = min!
401 #
402 # from which y = argmin_y || beta e_1 - H y ||_2, and q = -B y
403 #
405 #
406 # GCROT(m,k) update
407 #
409 # Define new outer vectors
411 # ux := (Z - U B) y
412 ux = zs[0]*y[0]
413 for z, yc in zip(zs[1:], y[1:]):
414 ux = axpy(z, ux, ux.shape[0], yc) # ux += z*yc
415 by = B.dot(y)
416 for cu, byc in zip(CU, by):
417 c, u = cu
418 ux = axpy(u, ux, ux.shape[0], -byc) # ux -= u*byc
420 # cx := V H y
421 hy = Q.dot(R.dot(y))
422 cx = vs[0] * hy[0]
423 for v, hyc in zip(vs[1:], hy[1:]):
424 cx = axpy(v, cx, cx.shape[0], hyc) # cx += v*hyc
426 # Normalize cx, maintaining cx = A ux
427 # This new cx is orthogonal to the previous C, by construction
428 try:
429 alpha = 1/nrm2(cx)
430 if not np.isfinite(alpha):
431 raise FloatingPointError()
432 except (FloatingPointError, ZeroDivisionError):
433 # Cannot update, so skip it
434 continue
436 cx = scal(alpha, cx)
437 ux = scal(alpha, ux)
439 # Update residual and solution
440 gamma = dot(cx, r)
441 r = axpy(cx, r, r.shape[0], -gamma) # r -= gamma*cx
442 x = axpy(ux, x, x.shape[0], gamma) # x += gamma*ux
444 # Truncate CU
445 if truncate == 'oldest':
446 while len(CU) >= k and CU:
447 del CU[0]
448 elif truncate == 'smallest':
449 if len(CU) >= k and CU:
450 # cf. [1,2]
451 D = solve(R[:-1,:].T, B.T).T
452 W, sigma, V = svd(D)
454 # C := C W[:,:k-1], U := U W[:,:k-1]
455 new_CU = []
456 for j, w in enumerate(W[:,:k-1].T):
457 c, u = CU[0]
458 c = c * w[0]
459 u = u * w[0]
460 for cup, wp in zip(CU[1:], w[1:]):
461 cp, up = cup
462 c = axpy(cp, c, c.shape[0], wp)
463 u = axpy(up, u, u.shape[0], wp)
465 # Reorthogonalize at the same time; not necessary
466 # in exact arithmetic, but floating point error
467 # tends to accumulate here
468 for cp, up in new_CU:
469 alpha = dot(cp, c)
470 c = axpy(cp, c, c.shape[0], -alpha)
471 u = axpy(up, u, u.shape[0], -alpha)
472 alpha = nrm2(c)
473 c = scal(1.0/alpha, c)
474 u = scal(1.0/alpha, u)
476 new_CU.append((c, u))
477 CU[:] = new_CU
479 # Add new vector to CU
480 CU.append((cx, ux))
482 # Include the solution vector to the span
483 CU.append((None, x.copy()))
484 if discard_C:
485 CU[:] = [(None, uz) for cz, uz in CU]
487 return postprocess(x), j_outer + 1