Coverage for /usr/lib/python3/dist-packages/scipy/sparse/linalg/_isolve/

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 cs : 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 def lpsolve(x):

65 return x

66 if rpsolve is None:

67 def rpsolve(x):

68 return x

70 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'], (v0,))

72 vs = [v0]

73 zs = []

74 y = None

75 res = np.nan

77 m = m + len(outer_v)

79 # Orthogonal projection coefficients

80 B = np.zeros((len(cs), m), dtype=v0.dtype)

82 # H is stored in QR factorized form

83 Q = np.ones((1, 1), dtype=v0.dtype)

84 R = np.zeros((1, 0), dtype=v0.dtype)

86 eps = np.finfo(v0.dtype).eps

88 breakdown = False

90 # FGMRES Arnoldi process

91 for j in range(m):

92 # L A Z = C B + V H

94 if prepend_outer_v and j < len(outer_v):

95 z, w = outer_v[j]

96 elif prepend_outer_v and j == len(outer_v):

97 z = rpsolve(v0)

98 w = None

99 elif not prepend_outer_v and j >= m - len(outer_v):

100 z, w = outer_v[j - (m - len(outer_v))]

101 else:

102 z = rpsolve(vs[-1])

103 w = None

104

105 if w is None:

106 w = lpsolve(matvec(z))

107 else:

108 # w is clobbered below

109 w = w.copy()

110

111 w_norm = nrm2(w)

112

113 # GCROT projection: L A -> (1 - C C^H) L A

114 # i.e. orthogonalize against C

115 for i, c in enumerate(cs):

116 alpha = dot(c, w)

117 B[i,j] = alpha

118 w = axpy(c, w, c.shape[0], -alpha) # w -= alpha*c

119

120 # Orthogonalize against V

121 hcur = np.zeros(j+2, dtype=Q.dtype)

122 for i, v in enumerate(vs):

123 alpha = dot(v, w)

124 hcur[i] = alpha

125 w = axpy(v, w, v.shape[0], -alpha) # w -= alpha*v

126 hcur[i+1] = nrm2(w)

127

128 with np.errstate(over='ignore', divide='ignore'):

129 # Careful with denormals

130 alpha = 1/hcur[-1]

131

132 if np.isfinite(alpha):

133 w = scal(alpha, w)

134

135 if not (hcur[-1] > eps * w_norm):

136 # w essentially in the span of previous vectors,

137 # or we have nans. Bail out after updating the QR

138 # solution.

139 breakdown = True

140

141 vs.append(w)

142 zs.append(z)

143

144 # Arnoldi LSQ problem

145

146 # Add new column to H=Q@R, padding other columns with zeros

147 Q2 = np.zeros((j+2, j+2), dtype=Q.dtype, order='F')

148 Q2[:j+1,:j+1] = Q

149 Q2[j+1,j+1] = 1

150

151 R2 = np.zeros((j+2, j), dtype=R.dtype, order='F')

152 R2[:j+1,:] = R

153

154 Q, R = qr_insert(Q2, R2, hcur, j, which='col',

155 overwrite_qru=True, check_finite=False)

156

157 # Transformed least squares problem

158 # || Q R y - inner_res_0 * e_1 ||_2 = min!

159 # Since R = [R'; 0], solution is y = inner_res_0 (R')^{-1} (Q^H)[:j,0]

160

161 # Residual is immediately known

162 res = abs(Q[0,-1])

163

164 # Check for termination

165 if res < atol or breakdown:

166 break

167

168 if not np.isfinite(R[j,j]):

169 # nans encountered, bail out

170 raise LinAlgError()

171

172 # -- Get the LSQ problem solution

173

174 # The problem is triangular, but the condition number may be

175 # bad (or in case of breakdown the last diagonal entry may be

176 # zero), so use lstsq instead of trtrs.

177 y, _, _, _, = lstsq(R[:j+1,:j+1], Q[0,:j+1].conj())

178

179 B = B[:,:j+1]

180

181 return Q, R, B, vs, zs, y, res

182

183

184def gcrotmk(A, b, x0=None, tol=1e-5, maxiter=1000, M=None, callback=None,

185 m=20, k=None, CU=None, discard_C=False, truncate='oldest',

186 atol=None):

187 """

188 Solve a matrix equation using flexible GCROT(m,k) algorithm.

189

190 Parameters

191 ----------

192 A : {sparse matrix, ndarray, LinearOperator}

193 The real or complex N-by-N matrix of the linear system.

194 Alternatively, ``A`` can be a linear operator which can

195 produce ``Ax`` using, e.g.,

196 ``scipy.sparse.linalg.LinearOperator``.

197 b : ndarray

198 Right hand side of the linear system. Has shape (N,) or (N,1).

199 x0 : ndarray

200 Starting guess for the solution.

201 tol, atol : float, optional

202 Tolerances for convergence, ``norm(residual) <= max(tol*norm(b), atol)``.

203 The default for ``atol`` is `tol`.

204

205 .. warning::

206

207 The default value for `atol` will be changed in a future release.

208 For future compatibility, specify `atol` explicitly.

209 maxiter : int, optional

210 Maximum number of iterations. Iteration will stop after maxiter

211 steps even if the specified tolerance has not been achieved.

212 M : {sparse matrix, ndarray, LinearOperator}, optional

213 Preconditioner for A. The preconditioner should approximate the

214 inverse of A. gcrotmk is a 'flexible' algorithm and the preconditioner

215 can vary from iteration to iteration. Effective preconditioning

216 dramatically improves the rate of convergence, which implies that

217 fewer iterations are needed to reach a given error tolerance.

218 callback : function, optional

219 User-supplied function to call after each iteration. It is called

220 as callback(xk), where xk is the current solution vector.

221 m : int, optional

222 Number of inner FGMRES iterations per each outer iteration.

223 Default: 20

224 k : int, optional

225 Number of vectors to carry between inner FGMRES iterations.

226 According to [2]_, good values are around m.

227 Default: m

228 CU : list of tuples, optional

229 List of tuples ``(c, u)`` which contain the columns of the matrices

230 C and U in the GCROT(m,k) algorithm. For details, see [2]_.

231 The list given and vectors contained in it are modified in-place.

232 If not given, start from empty matrices. The ``c`` elements in the

233 tuples can be ``None``, in which case the vectors are recomputed

234 via ``c = A u`` on start and orthogonalized as described in [3]_.

235 discard_C : bool, optional

236 Discard the C-vectors at the end. Useful if recycling Krylov subspaces

237 for different linear systems.

238 truncate : {'oldest', 'smallest'}, optional

239 Truncation scheme to use. Drop: oldest vectors, or vectors with

240 smallest singular values using the scheme discussed in [1,2].

241 See [2]_ for detailed comparison.

242 Default: 'oldest'

243

244 Returns

245 -------

246 x : ndarray

247 The solution found.

248 info : int

249 Provides convergence information:

250

251 * 0 : successful exit

252 * >0 : convergence to tolerance not achieved, number of iterations

253

254 Examples

255 --------

256 >>> import numpy as np

257 >>> from scipy.sparse import csc_matrix

258 >>> from scipy.sparse.linalg import gcrotmk

259 >>> R = np.random.randn(5, 5)

260 >>> A = csc_matrix(R)

261 >>> b = np.random.randn(5)

262 >>> x, exit_code = gcrotmk(A, b)

263 >>> print(exit_code)

264 0

265 >>> np.allclose(A.dot(x), b)

266 True

267

268 References

269 ----------

270 .. [1] E. de Sturler, ''Truncation strategies for optimal Krylov subspace

271 methods'', SIAM J. Numer. Anal. 36, 864 (1999).

272 .. [2] J.E. Hicken and D.W. Zingg, ''A simplified and flexible variant

273 of GCROT for solving nonsymmetric linear systems'',

274 SIAM J. Sci. Comput. 32, 172 (2010).

275 .. [3] M.L. Parks, E. de Sturler, G. Mackey, D.D. Johnson, S. Maiti,

276 ''Recycling Krylov subspaces for sequences of linear systems'',

277 SIAM J. Sci. Comput. 28, 1651 (2006).

278

279 """

280 A,M,x,b,postprocess = make_system(A,M,x0,b)

281

282 if not np.isfinite(b).all():

283 raise ValueError("RHS must contain only finite numbers")

284

285 if truncate not in ('oldest', 'smallest'):

286 raise ValueError(f"Invalid value for 'truncate': {truncate!r}")

287

288 if atol is None:

289 warnings.warn("scipy.sparse.linalg.gcrotmk called without specifying `atol`. "

290 "The default value will change in the future. To preserve "

291 "current behavior, set ``atol=tol``.",

292 category=DeprecationWarning, stacklevel=2)

293 atol = tol

294

295 matvec = A.matvec

296 psolve = M.matvec

297

298 if CU is None:

299 CU = []

300

301 if k is None:

302 k = m

303

304 axpy, dot, scal = None, None, None

305

306 if x0 is None:

307 r = b.copy()

308 else:

309 r = b - matvec(x)

310

311 axpy, dot, scal, nrm2 = get_blas_funcs(['axpy', 'dot', 'scal', 'nrm2'], (x, r))

312

313 b_norm = nrm2(b)

314 if b_norm == 0:

315 x = b

316 return (postprocess(x), 0)

317

318 if discard_C:

319 CU[:] = [(None, u) for c, u in CU]

320

321 # Reorthogonalize old vectors

322 if CU:

323 # Sort already existing vectors to the front

324 CU.sort(key=lambda cu: cu[0] is not None)

325

326 # Fill-in missing ones

327 C = np.empty((A.shape[0], len(CU)), dtype=r.dtype, order='F')

328 us = []

329 j = 0

330 while CU:

331 # More memory-efficient: throw away old vectors as we go

332 c, u = CU.pop(0)

333 if c is None:

334 c = matvec(u)

335 C[:,j] = c

336 j += 1

337 us.append(u)

338

339 # Orthogonalize

340 Q, R, P = qr(C, overwrite_a=True, mode='economic', pivoting=True)

341 del C

342

343 # C := Q

344 cs = list(Q.T)

345

346 # U := U P R^-1, back-substitution

347 new_us = []

348 for j in range(len(cs)):

349 u = us[P[j]]

350 for i in range(j):

351 u = axpy(us[P[i]], u, u.shape[0], -R[i,j])

352 if abs(R[j,j]) < 1e-12 * abs(R[0,0]):

353 # discard rest of the vectors

354 break

355 u = scal(1.0/R[j,j], u)

356 new_us.append(u)

357

358 # Form the new CU lists

359 CU[:] = list(zip(cs, new_us))[::-1]

360

361 if CU:

362 axpy, dot = get_blas_funcs(['axpy', 'dot'], (r,))

363

364 # Solve first the projection operation with respect to the CU

365 # vectors. This corresponds to modifying the initial guess to

366 # be

367 #

368 # x' = x + U y

369 # y = argmin_y || b - A (x + U y) ||^2

370 #

371 # The solution is y = C^H (b - A x)

372 for c, u in CU:

373 yc = dot(c, r)

374 x = axpy(u, x, x.shape[0], yc)

375 r = axpy(c, r, r.shape[0], -yc)

376

377 # GCROT main iteration

378 for j_outer in range(maxiter):

379 # -- callback

380 if callback is not None:

381 callback(x)

382

383 beta = nrm2(r)

384

385 # -- check stopping condition

386 beta_tol = max(atol, tol * b_norm)

387

388 if beta <= beta_tol and (j_outer > 0 or CU):

389 # recompute residual to avoid rounding error

390 r = b - matvec(x)

391 beta = nrm2(r)

392

393 if beta <= beta_tol:

394 j_outer = -1

395 break

396

397 ml = m + max(k - len(CU), 0)

398

399 cs = [c for c, u in CU]

400

401 try:

402 Q, R, B, vs, zs, y, pres = _fgmres(matvec,

403 r/beta,

404 ml,

405 rpsolve=psolve,

406 atol=max(atol, tol*b_norm)/beta,

407 cs=cs)

408 y *= beta

409 except LinAlgError:

410 # Floating point over/underflow, non-finite result from

411 # matmul etc. -- report failure.

412 break

413

414 #

415 # At this point,

416 #

417 # [A U, A Z] = [C, V] G; G = [ I B ]

418 # [ 0 H ]

419 #

420 # where [C, V] has orthonormal columns, and r = beta v_0. Moreover,

421 #

422 # || b - A (x + Z y + U q) ||_2 = || r - C B y - V H y - C q ||_2 = min!

423 #

424 # from which y = argmin_y || beta e_1 - H y ||_2, and q = -B y

425 #

426

427 #

428 # GCROT(m,k) update

429 #

430

431 # Define new outer vectors

432

433 # ux := (Z - U B) y

434 ux = zs[0]*y[0]

435 for z, yc in zip(zs[1:], y[1:]):

436 ux = axpy(z, ux, ux.shape[0], yc) # ux += z*yc

437 by = B.dot(y)

438 for cu, byc in zip(CU, by):

439 c, u = cu

440 ux = axpy(u, ux, ux.shape[0], -byc) # ux -= u*byc

441

442 # cx := V H y

443 hy = Q.dot(R.dot(y))

444 cx = vs[0] * hy[0]

445 for v, hyc in zip(vs[1:], hy[1:]):

446 cx = axpy(v, cx, cx.shape[0], hyc) # cx += v*hyc

447

448 # Normalize cx, maintaining cx = A ux

449 # This new cx is orthogonal to the previous C, by construction

450 try:

451 alpha = 1/nrm2(cx)

452 if not np.isfinite(alpha):

453 raise FloatingPointError()

454 except (FloatingPointError, ZeroDivisionError):

455 # Cannot update, so skip it

456 continue

457

458 cx = scal(alpha, cx)

459 ux = scal(alpha, ux)

460

461 # Update residual and solution

462 gamma = dot(cx, r)

463 r = axpy(cx, r, r.shape[0], -gamma) # r -= gamma*cx

464 x = axpy(ux, x, x.shape[0], gamma) # x += gamma*ux

465

466 # Truncate CU

467 if truncate == 'oldest':

468 while len(CU) >= k and CU:

469 del CU[0]

470 elif truncate == 'smallest':

471 if len(CU) >= k and CU:

472 # cf. [1,2]

473 D = solve(R[:-1,:].T, B.T).T

474 W, sigma, V = svd(D)

475

476 # C := C W[:,:k-1], U := U W[:,:k-1]

477 new_CU = []

478 for j, w in enumerate(W[:,:k-1].T):

479 c, u = CU[0]

480 c = c * w[0]

481 u = u * w[0]

482 for cup, wp in zip(CU[1:], w[1:]):

483 cp, up = cup

484 c = axpy(cp, c, c.shape[0], wp)

485 u = axpy(up, u, u.shape[0], wp)

486

487 # Reorthogonalize at the same time; not necessary

488 # in exact arithmetic, but floating point error

489 # tends to accumulate here

490 for cp, up in new_CU:

491 alpha = dot(cp, c)

492 c = axpy(cp, c, c.shape[0], -alpha)

493 u = axpy(up, u, u.shape[0], -alpha)

494 alpha = nrm2(c)

495 c = scal(1.0/alpha, c)

496 u = scal(1.0/alpha, u)

497

498 new_CU.append((c, u))

499 CU[:] = new_CU

500

501 # Add new vector to CU

502 CU.append((cx, ux))

503

504 # Include the solution vector to the span

505 CU.append((None, x.copy()))

506 if discard_C:

507 CU[:] = [(None, uz) for cz, uz in CU]

508

509 return postprocess(x), j_outer + 1

Coverage for /usr/lib/python3/dist-packages/scipy/sparse/linalg/_isolve/_gcrotmk.py: 4%

197 statements