Hide keyboard shortcuts

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"""Compute the action of the matrix exponential. 

2""" 

3 

4import numpy as np 

5 

6import scipy.linalg 

7import scipy.sparse.linalg 

8from scipy.sparse.linalg import aslinearoperator 

9from scipy.sparse.sputils import is_pydata_spmatrix 

10 

11__all__ = ['expm_multiply'] 

12 

13 

14def _exact_inf_norm(A): 

15 # A compatibility function which should eventually disappear. 

16 if scipy.sparse.isspmatrix(A): 

17 return max(abs(A).sum(axis=1).flat) 

18 elif is_pydata_spmatrix(A): 

19 return max(abs(A).sum(axis=1)) 

20 else: 

21 return np.linalg.norm(A, np.inf) 

22 

23 

24def _exact_1_norm(A): 

25 # A compatibility function which should eventually disappear. 

26 if scipy.sparse.isspmatrix(A): 

27 return max(abs(A).sum(axis=0).flat) 

28 elif is_pydata_spmatrix(A): 

29 return max(abs(A).sum(axis=0)) 

30 else: 

31 return np.linalg.norm(A, 1) 

32 

33 

34def _trace(A): 

35 # A compatibility function which should eventually disappear. 

36 if scipy.sparse.isspmatrix(A): 

37 return A.diagonal().sum() 

38 elif is_pydata_spmatrix(A): 

39 return A.to_scipy_sparse().diagonal().sum() 

40 else: 

41 return np.trace(A) 

42 

43 

44def _ident_like(A): 

45 # A compatibility function which should eventually disappear. 

46 if scipy.sparse.isspmatrix(A): 

47 return scipy.sparse.construct.eye(A.shape[0], A.shape[1], 

48 dtype=A.dtype, format=A.format) 

49 elif is_pydata_spmatrix(A): 

50 import sparse 

51 return sparse.eye(A.shape[0], A.shape[1], dtype=A.dtype) 

52 else: 

53 return np.eye(A.shape[0], A.shape[1], dtype=A.dtype) 

54 

55 

56def expm_multiply(A, B, start=None, stop=None, num=None, endpoint=None): 

57 """ 

58 Compute the action of the matrix exponential of A on B. 

59 

60 Parameters 

61 ---------- 

62 A : transposable linear operator 

63 The operator whose exponential is of interest. 

64 B : ndarray 

65 The matrix or vector to be multiplied by the matrix exponential of A. 

66 start : scalar, optional 

67 The starting time point of the sequence. 

68 stop : scalar, optional 

69 The end time point of the sequence, unless `endpoint` is set to False. 

70 In that case, the sequence consists of all but the last of ``num + 1`` 

71 evenly spaced time points, so that `stop` is excluded. 

72 Note that the step size changes when `endpoint` is False. 

73 num : int, optional 

74 Number of time points to use. 

75 endpoint : bool, optional 

76 If True, `stop` is the last time point. Otherwise, it is not included. 

77 

78 Returns 

79 ------- 

80 expm_A_B : ndarray 

81 The result of the action :math:`e^{t_k A} B`. 

82 

83 Notes 

84 ----- 

85 The optional arguments defining the sequence of evenly spaced time points 

86 are compatible with the arguments of `numpy.linspace`. 

87 

88 The output ndarray shape is somewhat complicated so I explain it here. 

89 The ndim of the output could be either 1, 2, or 3. 

90 It would be 1 if you are computing the expm action on a single vector 

91 at a single time point. 

92 It would be 2 if you are computing the expm action on a vector 

93 at multiple time points, or if you are computing the expm action 

94 on a matrix at a single time point. 

95 It would be 3 if you want the action on a matrix with multiple 

96 columns at multiple time points. 

97 If multiple time points are requested, expm_A_B[0] will always 

98 be the action of the expm at the first time point, 

99 regardless of whether the action is on a vector or a matrix. 

100 

101 References 

102 ---------- 

103 .. [1] Awad H. Al-Mohy and Nicholas J. Higham (2011) 

104 "Computing the Action of the Matrix Exponential, 

105 with an Application to Exponential Integrators." 

106 SIAM Journal on Scientific Computing, 

107 33 (2). pp. 488-511. ISSN 1064-8275 

108 http://eprints.ma.man.ac.uk/1591/ 

109 

110 .. [2] Nicholas J. Higham and Awad H. Al-Mohy (2010) 

111 "Computing Matrix Functions." 

112 Acta Numerica, 

113 19. 159-208. ISSN 0962-4929 

114 http://eprints.ma.man.ac.uk/1451/ 

115 

116 Examples 

117 -------- 

118 >>> from scipy.sparse import csc_matrix 

119 >>> from scipy.sparse.linalg import expm, expm_multiply 

120 >>> A = csc_matrix([[1, 0], [0, 1]]) 

121 >>> A.todense() 

122 matrix([[1, 0], 

123 [0, 1]], dtype=int64) 

124 >>> B = np.array([np.exp(-1.), np.exp(-2.)]) 

125 >>> B 

126 array([ 0.36787944, 0.13533528]) 

127 >>> expm_multiply(A, B, start=1, stop=2, num=3, endpoint=True) 

128 array([[ 1. , 0.36787944], 

129 [ 1.64872127, 0.60653066], 

130 [ 2.71828183, 1. ]]) 

131 >>> expm(A).dot(B) # Verify 1st timestep 

132 array([ 1. , 0.36787944]) 

133 >>> expm(1.5*A).dot(B) # Verify 2nd timestep 

134 array([ 1.64872127, 0.60653066]) 

135 >>> expm(2*A).dot(B) # Verify 3rd timestep 

136 array([ 2.71828183, 1. ]) 

137 """ 

138 if all(arg is None for arg in (start, stop, num, endpoint)): 

139 X = _expm_multiply_simple(A, B) 

140 else: 

141 X, status = _expm_multiply_interval(A, B, start, stop, num, endpoint) 

142 return X 

143 

144 

145def _expm_multiply_simple(A, B, t=1.0, balance=False): 

146 """ 

147 Compute the action of the matrix exponential at a single time point. 

148 

149 Parameters 

150 ---------- 

151 A : transposable linear operator 

152 The operator whose exponential is of interest. 

153 B : ndarray 

154 The matrix to be multiplied by the matrix exponential of A. 

155 t : float 

156 A time point. 

157 balance : bool 

158 Indicates whether or not to apply balancing. 

159 

160 Returns 

161 ------- 

162 F : ndarray 

163 :math:`e^{t A} B` 

164 

165 Notes 

166 ----- 

167 This is algorithm (3.2) in Al-Mohy and Higham (2011). 

168 

169 """ 

170 if balance: 

171 raise NotImplementedError 

172 if len(A.shape) != 2 or A.shape[0] != A.shape[1]: 

173 raise ValueError('expected A to be like a square matrix') 

174 if A.shape[1] != B.shape[0]: 

175 raise ValueError('the matrices A and B have incompatible shapes') 

176 ident = _ident_like(A) 

177 n = A.shape[0] 

178 if len(B.shape) == 1: 

179 n0 = 1 

180 elif len(B.shape) == 2: 

181 n0 = B.shape[1] 

182 else: 

183 raise ValueError('expected B to be like a matrix or a vector') 

184 u_d = 2**-53 

185 tol = u_d 

186 mu = _trace(A) / float(n) 

187 A = A - mu * ident 

188 A_1_norm = _exact_1_norm(A) 

189 if t*A_1_norm == 0: 

190 m_star, s = 0, 1 

191 else: 

192 ell = 2 

193 norm_info = LazyOperatorNormInfo(t*A, A_1_norm=t*A_1_norm, ell=ell) 

194 m_star, s = _fragment_3_1(norm_info, n0, tol, ell=ell) 

195 return _expm_multiply_simple_core(A, B, t, mu, m_star, s, tol, balance) 

196 

197 

198def _expm_multiply_simple_core(A, B, t, mu, m_star, s, tol=None, balance=False): 

199 """ 

200 A helper function. 

201 """ 

202 if balance: 

203 raise NotImplementedError 

204 if tol is None: 

205 u_d = 2 ** -53 

206 tol = u_d 

207 F = B 

208 eta = np.exp(t*mu / float(s)) 

209 for i in range(s): 

210 c1 = _exact_inf_norm(B) 

211 for j in range(m_star): 

212 coeff = t / float(s*(j+1)) 

213 B = coeff * A.dot(B) 

214 c2 = _exact_inf_norm(B) 

215 F = F + B 

216 if c1 + c2 <= tol * _exact_inf_norm(F): 

217 break 

218 c1 = c2 

219 F = eta * F 

220 B = F 

221 return F 

222 

223 

224# This table helps to compute bounds. 

225# They seem to have been difficult to calculate, involving symbolic 

226# manipulation of equations, followed by numerical root finding. 

227_theta = { 

228 # The first 30 values are from table A.3 of Computing Matrix Functions. 

229 1: 2.29e-16, 

230 2: 2.58e-8, 

231 3: 1.39e-5, 

232 4: 3.40e-4, 

233 5: 2.40e-3, 

234 6: 9.07e-3, 

235 7: 2.38e-2, 

236 8: 5.00e-2, 

237 9: 8.96e-2, 

238 10: 1.44e-1, 

239 # 11 

240 11: 2.14e-1, 

241 12: 3.00e-1, 

242 13: 4.00e-1, 

243 14: 5.14e-1, 

244 15: 6.41e-1, 

245 16: 7.81e-1, 

246 17: 9.31e-1, 

247 18: 1.09, 

248 19: 1.26, 

249 20: 1.44, 

250 # 21 

251 21: 1.62, 

252 22: 1.82, 

253 23: 2.01, 

254 24: 2.22, 

255 25: 2.43, 

256 26: 2.64, 

257 27: 2.86, 

258 28: 3.08, 

259 29: 3.31, 

260 30: 3.54, 

261 # The rest are from table 3.1 of 

262 # Computing the Action of the Matrix Exponential. 

263 35: 4.7, 

264 40: 6.0, 

265 45: 7.2, 

266 50: 8.5, 

267 55: 9.9, 

268 } 

269 

270 

271def _onenormest_matrix_power(A, p, 

272 t=2, itmax=5, compute_v=False, compute_w=False): 

273 """ 

274 Efficiently estimate the 1-norm of A^p. 

275 

276 Parameters 

277 ---------- 

278 A : ndarray 

279 Matrix whose 1-norm of a power is to be computed. 

280 p : int 

281 Non-negative integer power. 

282 t : int, optional 

283 A positive parameter controlling the tradeoff between 

284 accuracy versus time and memory usage. 

285 Larger values take longer and use more memory 

286 but give more accurate output. 

287 itmax : int, optional 

288 Use at most this many iterations. 

289 compute_v : bool, optional 

290 Request a norm-maximizing linear operator input vector if True. 

291 compute_w : bool, optional 

292 Request a norm-maximizing linear operator output vector if True. 

293 

294 Returns 

295 ------- 

296 est : float 

297 An underestimate of the 1-norm of the sparse matrix. 

298 v : ndarray, optional 

299 The vector such that ||Av||_1 == est*||v||_1. 

300 It can be thought of as an input to the linear operator 

301 that gives an output with particularly large norm. 

302 w : ndarray, optional 

303 The vector Av which has relatively large 1-norm. 

304 It can be thought of as an output of the linear operator 

305 that is relatively large in norm compared to the input. 

306 

307 """ 

308 #XXX Eventually turn this into an API function in the _onenormest module, 

309 #XXX and remove its underscore, 

310 #XXX but wait until expm_multiply goes into scipy. 

311 return scipy.sparse.linalg.onenormest(aslinearoperator(A) ** p) 

312 

313class LazyOperatorNormInfo: 

314 """ 

315 Information about an operator is lazily computed. 

316 

317 The information includes the exact 1-norm of the operator, 

318 in addition to estimates of 1-norms of powers of the operator. 

319 This uses the notation of Computing the Action (2011). 

320 This class is specialized enough to probably not be of general interest 

321 outside of this module. 

322 

323 """ 

324 def __init__(self, A, A_1_norm=None, ell=2, scale=1): 

325 """ 

326 Provide the operator and some norm-related information. 

327 

328 Parameters 

329 ---------- 

330 A : linear operator 

331 The operator of interest. 

332 A_1_norm : float, optional 

333 The exact 1-norm of A. 

334 ell : int, optional 

335 A technical parameter controlling norm estimation quality. 

336 scale : int, optional 

337 If specified, return the norms of scale*A instead of A. 

338 

339 """ 

340 self._A = A 

341 self._A_1_norm = A_1_norm 

342 self._ell = ell 

343 self._d = {} 

344 self._scale = scale 

345 

346 def set_scale(self,scale): 

347 """ 

348 Set the scale parameter. 

349 """ 

350 self._scale = scale 

351 

352 def onenorm(self): 

353 """ 

354 Compute the exact 1-norm. 

355 """ 

356 if self._A_1_norm is None: 

357 self._A_1_norm = _exact_1_norm(self._A) 

358 return self._scale*self._A_1_norm 

359 

360 def d(self, p): 

361 """ 

362 Lazily estimate d_p(A) ~= || A^p ||^(1/p) where ||.|| is the 1-norm. 

363 """ 

364 if p not in self._d: 

365 est = _onenormest_matrix_power(self._A, p, self._ell) 

366 self._d[p] = est ** (1.0 / p) 

367 return self._scale*self._d[p] 

368 

369 def alpha(self, p): 

370 """ 

371 Lazily compute max(d(p), d(p+1)). 

372 """ 

373 return max(self.d(p), self.d(p+1)) 

374 

375def _compute_cost_div_m(m, p, norm_info): 

376 """ 

377 A helper function for computing bounds. 

378 

379 This is equation (3.10). 

380 It measures cost in terms of the number of required matrix products. 

381 

382 Parameters 

383 ---------- 

384 m : int 

385 A valid key of _theta. 

386 p : int 

387 A matrix power. 

388 norm_info : LazyOperatorNormInfo 

389 Information about 1-norms of related operators. 

390 

391 Returns 

392 ------- 

393 cost_div_m : int 

394 Required number of matrix products divided by m. 

395 

396 """ 

397 return int(np.ceil(norm_info.alpha(p) / _theta[m])) 

398 

399 

400def _compute_p_max(m_max): 

401 """ 

402 Compute the largest positive integer p such that p*(p-1) <= m_max + 1. 

403 

404 Do this in a slightly dumb way, but safe and not too slow. 

405 

406 Parameters 

407 ---------- 

408 m_max : int 

409 A count related to bounds. 

410 

411 """ 

412 sqrt_m_max = np.sqrt(m_max) 

413 p_low = int(np.floor(sqrt_m_max)) 

414 p_high = int(np.ceil(sqrt_m_max + 1)) 

415 return max(p for p in range(p_low, p_high+1) if p*(p-1) <= m_max + 1) 

416 

417 

418def _fragment_3_1(norm_info, n0, tol, m_max=55, ell=2): 

419 """ 

420 A helper function for the _expm_multiply_* functions. 

421 

422 Parameters 

423 ---------- 

424 norm_info : LazyOperatorNormInfo 

425 Information about norms of certain linear operators of interest. 

426 n0 : int 

427 Number of columns in the _expm_multiply_* B matrix. 

428 tol : float 

429 Expected to be 

430 :math:`2^{-24}` for single precision or 

431 :math:`2^{-53}` for double precision. 

432 m_max : int 

433 A value related to a bound. 

434 ell : int 

435 The number of columns used in the 1-norm approximation. 

436 This is usually taken to be small, maybe between 1 and 5. 

437 

438 Returns 

439 ------- 

440 best_m : int 

441 Related to bounds for error control. 

442 best_s : int 

443 Amount of scaling. 

444 

445 Notes 

446 ----- 

447 This is code fragment (3.1) in Al-Mohy and Higham (2011). 

448 The discussion of default values for m_max and ell 

449 is given between the definitions of equation (3.11) 

450 and the definition of equation (3.12). 

451 

452 """ 

453 if ell < 1: 

454 raise ValueError('expected ell to be a positive integer') 

455 best_m = None 

456 best_s = None 

457 if _condition_3_13(norm_info.onenorm(), n0, m_max, ell): 

458 for m, theta in _theta.items(): 

459 s = int(np.ceil(norm_info.onenorm() / theta)) 

460 if best_m is None or m * s < best_m * best_s: 

461 best_m = m 

462 best_s = s 

463 else: 

464 # Equation (3.11). 

465 for p in range(2, _compute_p_max(m_max) + 1): 

466 for m in range(p*(p-1)-1, m_max+1): 

467 if m in _theta: 

468 s = _compute_cost_div_m(m, p, norm_info) 

469 if best_m is None or m * s < best_m * best_s: 

470 best_m = m 

471 best_s = s 

472 best_s = max(best_s, 1) 

473 return best_m, best_s 

474 

475 

476def _condition_3_13(A_1_norm, n0, m_max, ell): 

477 """ 

478 A helper function for the _expm_multiply_* functions. 

479 

480 Parameters 

481 ---------- 

482 A_1_norm : float 

483 The precomputed 1-norm of A. 

484 n0 : int 

485 Number of columns in the _expm_multiply_* B matrix. 

486 m_max : int 

487 A value related to a bound. 

488 ell : int 

489 The number of columns used in the 1-norm approximation. 

490 This is usually taken to be small, maybe between 1 and 5. 

491 

492 Returns 

493 ------- 

494 value : bool 

495 Indicates whether or not the condition has been met. 

496 

497 Notes 

498 ----- 

499 This is condition (3.13) in Al-Mohy and Higham (2011). 

500 

501 """ 

502 

503 # This is the rhs of equation (3.12). 

504 p_max = _compute_p_max(m_max) 

505 a = 2 * ell * p_max * (p_max + 3) 

506 

507 # Evaluate the condition (3.13). 

508 b = _theta[m_max] / float(n0 * m_max) 

509 return A_1_norm <= a * b 

510 

511 

512def _expm_multiply_interval(A, B, start=None, stop=None, 

513 num=None, endpoint=None, balance=False, status_only=False): 

514 """ 

515 Compute the action of the matrix exponential at multiple time points. 

516 

517 Parameters 

518 ---------- 

519 A : transposable linear operator 

520 The operator whose exponential is of interest. 

521 B : ndarray 

522 The matrix to be multiplied by the matrix exponential of A. 

523 start : scalar, optional 

524 The starting time point of the sequence. 

525 stop : scalar, optional 

526 The end time point of the sequence, unless `endpoint` is set to False. 

527 In that case, the sequence consists of all but the last of ``num + 1`` 

528 evenly spaced time points, so that `stop` is excluded. 

529 Note that the step size changes when `endpoint` is False. 

530 num : int, optional 

531 Number of time points to use. 

532 endpoint : bool, optional 

533 If True, `stop` is the last time point. Otherwise, it is not included. 

534 balance : bool 

535 Indicates whether or not to apply balancing. 

536 status_only : bool 

537 A flag that is set to True for some debugging and testing operations. 

538 

539 Returns 

540 ------- 

541 F : ndarray 

542 :math:`e^{t_k A} B` 

543 status : int 

544 An integer status for testing and debugging. 

545 

546 Notes 

547 ----- 

548 This is algorithm (5.2) in Al-Mohy and Higham (2011). 

549 

550 There seems to be a typo, where line 15 of the algorithm should be 

551 moved to line 6.5 (between lines 6 and 7). 

552 

553 """ 

554 if balance: 

555 raise NotImplementedError 

556 if len(A.shape) != 2 or A.shape[0] != A.shape[1]: 

557 raise ValueError('expected A to be like a square matrix') 

558 if A.shape[1] != B.shape[0]: 

559 raise ValueError('the matrices A and B have incompatible shapes') 

560 ident = _ident_like(A) 

561 n = A.shape[0] 

562 if len(B.shape) == 1: 

563 n0 = 1 

564 elif len(B.shape) == 2: 

565 n0 = B.shape[1] 

566 else: 

567 raise ValueError('expected B to be like a matrix or a vector') 

568 u_d = 2**-53 

569 tol = u_d 

570 mu = _trace(A) / float(n) 

571 

572 # Get the linspace samples, attempting to preserve the linspace defaults. 

573 linspace_kwargs = {'retstep': True} 

574 if num is not None: 

575 linspace_kwargs['num'] = num 

576 if endpoint is not None: 

577 linspace_kwargs['endpoint'] = endpoint 

578 samples, step = np.linspace(start, stop, **linspace_kwargs) 

579 

580 # Convert the linspace output to the notation used by the publication. 

581 nsamples = len(samples) 

582 if nsamples < 2: 

583 raise ValueError('at least two time points are required') 

584 q = nsamples - 1 

585 h = step 

586 t_0 = samples[0] 

587 t_q = samples[q] 

588 

589 # Define the output ndarray. 

590 # Use an ndim=3 shape, such that the last two indices 

591 # are the ones that may be involved in level 3 BLAS operations. 

592 X_shape = (nsamples,) + B.shape 

593 X = np.empty(X_shape, dtype=np.result_type(A.dtype, B.dtype, float)) 

594 t = t_q - t_0 

595 A = A - mu * ident 

596 A_1_norm = _exact_1_norm(A) 

597 ell = 2 

598 norm_info = LazyOperatorNormInfo(t*A, A_1_norm=t*A_1_norm, ell=ell) 

599 if t*A_1_norm == 0: 

600 m_star, s = 0, 1 

601 else: 

602 m_star, s = _fragment_3_1(norm_info, n0, tol, ell=ell) 

603 

604 # Compute the expm action up to the initial time point. 

605 X[0] = _expm_multiply_simple_core(A, B, t_0, mu, m_star, s) 

606 

607 # Compute the expm action at the rest of the time points. 

608 if q <= s: 

609 if status_only: 

610 return 0 

611 else: 

612 return _expm_multiply_interval_core_0(A, X, 

613 h, mu, q, norm_info, tol, ell,n0) 

614 elif not (q % s): 

615 if status_only: 

616 return 1 

617 else: 

618 return _expm_multiply_interval_core_1(A, X, 

619 h, mu, m_star, s, q, tol) 

620 elif (q % s): 

621 if status_only: 

622 return 2 

623 else: 

624 return _expm_multiply_interval_core_2(A, X, 

625 h, mu, m_star, s, q, tol) 

626 else: 

627 raise Exception('internal error') 

628 

629 

630def _expm_multiply_interval_core_0(A, X, h, mu, q, norm_info, tol, ell, n0): 

631 """ 

632 A helper function, for the case q <= s. 

633 """ 

634 

635 # Compute the new values of m_star and s which should be applied 

636 # over intervals of size t/q 

637 if norm_info.onenorm() == 0: 

638 m_star, s = 0, 1 

639 else: 

640 norm_info.set_scale(1./q) 

641 m_star, s = _fragment_3_1(norm_info, n0, tol, ell=ell) 

642 norm_info.set_scale(1) 

643 

644 for k in range(q): 

645 X[k+1] = _expm_multiply_simple_core(A, X[k], h, mu, m_star, s) 

646 return X, 0 

647 

648 

649def _expm_multiply_interval_core_1(A, X, h, mu, m_star, s, q, tol): 

650 """ 

651 A helper function, for the case q > s and q % s == 0. 

652 """ 

653 d = q // s 

654 input_shape = X.shape[1:] 

655 K_shape = (m_star + 1, ) + input_shape 

656 K = np.empty(K_shape, dtype=X.dtype) 

657 for i in range(s): 

658 Z = X[i*d] 

659 K[0] = Z 

660 high_p = 0 

661 for k in range(1, d+1): 

662 F = K[0] 

663 c1 = _exact_inf_norm(F) 

664 for p in range(1, m_star+1): 

665 if p > high_p: 

666 K[p] = h * A.dot(K[p-1]) / float(p) 

667 coeff = float(pow(k, p)) 

668 F = F + coeff * K[p] 

669 inf_norm_K_p_1 = _exact_inf_norm(K[p]) 

670 c2 = coeff * inf_norm_K_p_1 

671 if c1 + c2 <= tol * _exact_inf_norm(F): 

672 break 

673 c1 = c2 

674 X[k + i*d] = np.exp(k*h*mu) * F 

675 return X, 1 

676 

677 

678def _expm_multiply_interval_core_2(A, X, h, mu, m_star, s, q, tol): 

679 """ 

680 A helper function, for the case q > s and q % s > 0. 

681 """ 

682 d = q // s 

683 j = q // d 

684 r = q - d * j 

685 input_shape = X.shape[1:] 

686 K_shape = (m_star + 1, ) + input_shape 

687 K = np.empty(K_shape, dtype=X.dtype) 

688 for i in range(j + 1): 

689 Z = X[i*d] 

690 K[0] = Z 

691 high_p = 0 

692 if i < j: 

693 effective_d = d 

694 else: 

695 effective_d = r 

696 for k in range(1, effective_d+1): 

697 F = K[0] 

698 c1 = _exact_inf_norm(F) 

699 for p in range(1, m_star+1): 

700 if p == high_p + 1: 

701 K[p] = h * A.dot(K[p-1]) / float(p) 

702 high_p = p 

703 coeff = float(pow(k, p)) 

704 F = F + coeff * K[p] 

705 inf_norm_K_p_1 = _exact_inf_norm(K[p]) 

706 c2 = coeff * inf_norm_K_p_1 

707 if c1 + c2 <= tol * _exact_inf_norm(F): 

708 break 

709 c1 = c2 

710 X[k + i*d] = np.exp(k*h*mu) * F 

711 return X, 2