Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_matfuncs.py: 20%

1""" 

2Sparse matrix functions 

3""" 

4 

5# 

6# Authors: Travis Oliphant, March 2002 

7# Anthony Scopatz, August 2012 (Sparse Updates) 

8# Jake Vanderplas, August 2012 (Sparse Updates) 

9# 

10 

11__all__ = ['expm', 'inv', 'matrix_power'] 

12 

13import numpy as np 

14from scipy.linalg._basic import solve, solve_triangular 

15 

16from scipy.sparse._base import issparse 

17from scipy.sparse.linalg import spsolve 

18from scipy.sparse._sputils import is_pydata_spmatrix, isintlike 

19 

20import scipy.sparse 

21import scipy.sparse.linalg 

22from scipy.sparse.linalg._interface import LinearOperator 

23from scipy.sparse._construct import eye 

24 

25from ._expm_multiply import _ident_like, _exact_1_norm as _onenorm 

26 

27 

28UPPER_TRIANGULAR = 'upper_triangular' 

29 

30 

31def inv(A): 

32 """ 

33 Compute the inverse of a sparse matrix 

34 

35 Parameters 

36 ---------- 

37 A : (M, M) sparse matrix 

38 square matrix to be inverted 

39 

40 Returns 

41 ------- 

42 Ainv : (M, M) sparse matrix 

43 inverse of `A` 

44 

45 Notes 

46 ----- 

47 This computes the sparse inverse of `A`. If the inverse of `A` is expected 

48 to be non-sparse, it will likely be faster to convert `A` to dense and use 

49 `scipy.linalg.inv`. 

50 

51 Examples 

52 -------- 

53 >>> from scipy.sparse import csc_matrix 

54 >>> from scipy.sparse.linalg import inv 

55 >>> A = csc_matrix([[1., 0.], [1., 2.]]) 

56 >>> Ainv = inv(A) 

57 >>> Ainv 

58 <2x2 sparse matrix of type '<class 'numpy.float64'>' 

59 with 3 stored elements in Compressed Sparse Column format> 

60 >>> A.dot(Ainv) 

61 <2x2 sparse matrix of type '<class 'numpy.float64'>' 

62 with 2 stored elements in Compressed Sparse Column format> 

63 >>> A.dot(Ainv).toarray() 

64 array([[ 1., 0.], 

65 [ 0., 1.]]) 

66 

67 .. versionadded:: 0.12.0 

68 

69 """ 

70 # Check input 

71 if not (scipy.sparse.issparse(A) or is_pydata_spmatrix(A)): 

72 raise TypeError('Input must be a sparse matrix') 

73 

74 # Use sparse direct solver to solve "AX = I" accurately 

75 I = _ident_like(A) 

76 Ainv = spsolve(A, I) 

77 return Ainv 

78 

79 

80def _onenorm_matrix_power_nnm(A, p): 

81 """ 

82 Compute the 1-norm of a non-negative integer power of a non-negative matrix. 

83 

84 Parameters 

85 ---------- 

86 A : a square ndarray or matrix or sparse matrix 

87 Input matrix with non-negative entries. 

88 p : non-negative integer 

89 The power to which the matrix is to be raised. 

90 

91 Returns 

92 ------- 

93 out : float 

94 The 1-norm of the matrix power p of A. 

95 

96 """ 

97 # Check input 

98 if int(p) != p or p < 0: 

99 raise ValueError('expected non-negative integer p') 

100 p = int(p) 

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

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

103 

104 # Explicitly make a column vector so that this works when A is a 

105 # numpy matrix (in addition to ndarray and sparse matrix). 

106 v = np.ones((A.shape[0], 1), dtype=float) 

107 M = A.T 

108 for i in range(p): 

109 v = M.dot(v) 

110 return np.max(v) 

111 

112 

113def _is_upper_triangular(A): 

114 # This function could possibly be of wider interest. 

115 if issparse(A): 

116 lower_part = scipy.sparse.tril(A, -1) 

117 # Check structural upper triangularity, 

118 # then coincidental upper triangularity if needed. 

119 return lower_part.nnz == 0 or lower_part.count_nonzero() == 0 

120 elif is_pydata_spmatrix(A): 

121 import sparse 

122 lower_part = sparse.tril(A, -1) 

123 return lower_part.nnz == 0 

124 else: 

125 return not np.tril(A, -1).any() 

126 

127 

128def _smart_matrix_product(A, B, alpha=None, structure=None): 

129 """ 

130 A matrix product that knows about sparse and structured matrices. 

131 

132 Parameters 

133 ---------- 

134 A : 2d ndarray 

135 First matrix. 

136 B : 2d ndarray 

137 Second matrix. 

138 alpha : float 

139 The matrix product will be scaled by this constant. 

140 structure : str, optional 

141 A string describing the structure of both matrices `A` and `B`. 

142 Only `upper_triangular` is currently supported. 

143 

144 Returns 

145 ------- 

146 M : 2d ndarray 

147 Matrix product of A and B. 

148 

149 """ 

150 if len(A.shape) != 2: 

151 raise ValueError('expected A to be a rectangular matrix') 

152 if len(B.shape) != 2: 

153 raise ValueError('expected B to be a rectangular matrix') 

154 f = None 

155 if structure == UPPER_TRIANGULAR: 

156 if (not issparse(A) and not issparse(B) 

157 and not is_pydata_spmatrix(A) and not is_pydata_spmatrix(B)): 

158 f, = scipy.linalg.get_blas_funcs(('trmm',), (A, B)) 

159 if f is not None: 

160 if alpha is None: 

161 alpha = 1. 

162 out = f(alpha, A, B) 

163 else: 

164 if alpha is None: 

165 out = A.dot(B) 

166 else: 

167 out = alpha * A.dot(B) 

168 return out 

169 

170 

171class MatrixPowerOperator(LinearOperator): 

172 

173 def __init__(self, A, p, structure=None): 

174 if A.ndim != 2 or A.shape[0] != A.shape[1]: 

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

176 if p < 0: 

177 raise ValueError('expected p to be a non-negative integer') 

178 self._A = A 

179 self._p = p 

180 self._structure = structure 

181 self.dtype = A.dtype 

182 self.ndim = A.ndim 

183 self.shape = A.shape 

184 

185 def _matvec(self, x): 

186 for i in range(self._p): 

187 x = self._A.dot(x) 

188 return x 

189 

190 def _rmatvec(self, x): 

191 A_T = self._A.T 

192 x = x.ravel() 

193 for i in range(self._p): 

194 x = A_T.dot(x) 

195 return x 

196 

197 def _matmat(self, X): 

198 for i in range(self._p): 

199 X = _smart_matrix_product(self._A, X, structure=self._structure) 

200 return X 

201 

202 @property 

203 def T(self): 

204 return MatrixPowerOperator(self._A.T, self._p) 

205 

206 

207class ProductOperator(LinearOperator): 

208 """ 

209 For now, this is limited to products of multiple square matrices. 

210 """ 

211 

212 def __init__(self, *args, **kwargs): 

213 self._structure = kwargs.get('structure', None) 

214 for A in args: 

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

216 raise ValueError( 

217 'For now, the ProductOperator implementation is ' 

218 'limited to the product of multiple square matrices.') 

219 if args: 

220 n = args[0].shape[0] 

221 for A in args: 

222 for d in A.shape: 

223 if d != n: 

224 raise ValueError( 

225 'The square matrices of the ProductOperator ' 

226 'must all have the same shape.') 

227 self.shape = (n, n) 

228 self.ndim = len(self.shape) 

229 self.dtype = np.result_type(*[x.dtype for x in args]) 

230 self._operator_sequence = args 

231 

232 def _matvec(self, x): 

233 for A in reversed(self._operator_sequence): 

234 x = A.dot(x) 

235 return x 

236 

237 def _rmatvec(self, x): 

238 x = x.ravel() 

239 for A in self._operator_sequence: 

240 x = A.T.dot(x) 

241 return x 

242 

243 def _matmat(self, X): 

244 for A in reversed(self._operator_sequence): 

245 X = _smart_matrix_product(A, X, structure=self._structure) 

246 return X 

247 

248 @property 

249 def T(self): 

250 T_args = [A.T for A in reversed(self._operator_sequence)] 

251 return ProductOperator(*T_args) 

252 

253 

254def _onenormest_matrix_power(A, p, 

255 t=2, itmax=5, compute_v=False, compute_w=False, structure=None): 

256 """ 

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

258 

259 Parameters 

260 ---------- 

261 A : ndarray 

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

263 p : int 

264 Non-negative integer power. 

265 t : int, optional 

266 A positive parameter controlling the tradeoff between 

267 accuracy versus time and memory usage. 

268 Larger values take longer and use more memory 

269 but give more accurate output. 

270 itmax : int, optional 

271 Use at most this many iterations. 

272 compute_v : bool, optional 

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

274 compute_w : bool, optional 

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

276 

277 Returns 

278 ------- 

279 est : float 

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

281 v : ndarray, optional 

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

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

284 that gives an output with particularly large norm. 

285 w : ndarray, optional 

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

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

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

289 

290 """ 

291 return scipy.sparse.linalg.onenormest( 

292 MatrixPowerOperator(A, p, structure=structure)) 

293 

294 

295def _onenormest_product(operator_seq, 

296 t=2, itmax=5, compute_v=False, compute_w=False, structure=None): 

297 """ 

298 Efficiently estimate the 1-norm of the matrix product of the args. 

299 

300 Parameters 

301 ---------- 

302 operator_seq : linear operator sequence 

303 Matrices whose 1-norm of product is to be computed. 

304 t : int, optional 

305 A positive parameter controlling the tradeoff between 

306 accuracy versus time and memory usage. 

307 Larger values take longer and use more memory 

308 but give more accurate output. 

309 itmax : int, optional 

310 Use at most this many iterations. 

311 compute_v : bool, optional 

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

313 compute_w : bool, optional 

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

315 structure : str, optional 

316 A string describing the structure of all operators. 

317 Only `upper_triangular` is currently supported. 

318 

319 Returns 

320 ------- 

321 est : float 

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

323 v : ndarray, optional 

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

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

326 that gives an output with particularly large norm. 

327 w : ndarray, optional 

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

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

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

331 

332 """ 

333 return scipy.sparse.linalg.onenormest( 

334 ProductOperator(*operator_seq, structure=structure)) 

335 

336 

337class _ExpmPadeHelper: 

338 """ 

339 Help lazily evaluate a matrix exponential. 

340 

341 The idea is to not do more work than we need for high expm precision, 

342 so we lazily compute matrix powers and store or precompute 

343 other properties of the matrix. 

344 

345 """ 

346 

347 def __init__(self, A, structure=None, use_exact_onenorm=False): 

348 """ 

349 Initialize the object. 

350 

351 Parameters 

352 ---------- 

353 A : a dense or sparse square numpy matrix or ndarray 

354 The matrix to be exponentiated. 

355 structure : str, optional 

356 A string describing the structure of matrix `A`. 

357 Only `upper_triangular` is currently supported. 

358 use_exact_onenorm : bool, optional 

359 If True then only the exact one-norm of matrix powers and products 

360 will be used. Otherwise, the one-norm of powers and products 

361 may initially be estimated. 

362 """ 

363 self.A = A 

364 self._A2 = None 

365 self._A4 = None 

366 self._A6 = None 

367 self._A8 = None 

368 self._A10 = None 

369 self._d4_exact = None 

370 self._d6_exact = None 

371 self._d8_exact = None 

372 self._d10_exact = None 

373 self._d4_approx = None 

374 self._d6_approx = None 

375 self._d8_approx = None 

376 self._d10_approx = None 

377 self.ident = _ident_like(A) 

378 self.structure = structure 

379 self.use_exact_onenorm = use_exact_onenorm 

380 

381 @property 

382 def A2(self): 

383 if self._A2 is None: 

384 self._A2 = _smart_matrix_product( 

385 self.A, self.A, structure=self.structure) 

386 return self._A2 

387 

388 @property 

389 def A4(self): 

390 if self._A4 is None: 

391 self._A4 = _smart_matrix_product( 

392 self.A2, self.A2, structure=self.structure) 

393 return self._A4 

394 

395 @property 

396 def A6(self): 

397 if self._A6 is None: 

398 self._A6 = _smart_matrix_product( 

399 self.A4, self.A2, structure=self.structure) 

400 return self._A6 

401 

402 @property 

403 def A8(self): 

404 if self._A8 is None: 

405 self._A8 = _smart_matrix_product( 

406 self.A6, self.A2, structure=self.structure) 

407 return self._A8 

408 

409 @property 

410 def A10(self): 

411 if self._A10 is None: 

412 self._A10 = _smart_matrix_product( 

413 self.A4, self.A6, structure=self.structure) 

414 return self._A10 

415 

416 @property 

417 def d4_tight(self): 

418 if self._d4_exact is None: 

419 self._d4_exact = _onenorm(self.A4)**(1/4.) 

420 return self._d4_exact 

421 

422 @property 

423 def d6_tight(self): 

424 if self._d6_exact is None: 

425 self._d6_exact = _onenorm(self.A6)**(1/6.) 

426 return self._d6_exact 

427 

428 @property 

429 def d8_tight(self): 

430 if self._d8_exact is None: 

431 self._d8_exact = _onenorm(self.A8)**(1/8.) 

432 return self._d8_exact 

433 

434 @property 

435 def d10_tight(self): 

436 if self._d10_exact is None: 

437 self._d10_exact = _onenorm(self.A10)**(1/10.) 

438 return self._d10_exact 

439 

440 @property 

441 def d4_loose(self): 

442 if self.use_exact_onenorm: 

443 return self.d4_tight 

444 if self._d4_exact is not None: 

445 return self._d4_exact 

446 else: 

447 if self._d4_approx is None: 

448 self._d4_approx = _onenormest_matrix_power(self.A2, 2, 

449 structure=self.structure)**(1/4.) 

450 return self._d4_approx 

451 

452 @property 

453 def d6_loose(self): 

454 if self.use_exact_onenorm: 

455 return self.d6_tight 

456 if self._d6_exact is not None: 

457 return self._d6_exact 

458 else: 

459 if self._d6_approx is None: 

460 self._d6_approx = _onenormest_matrix_power(self.A2, 3, 

461 structure=self.structure)**(1/6.) 

462 return self._d6_approx 

463 

464 @property 

465 def d8_loose(self): 

466 if self.use_exact_onenorm: 

467 return self.d8_tight 

468 if self._d8_exact is not None: 

469 return self._d8_exact 

470 else: 

471 if self._d8_approx is None: 

472 self._d8_approx = _onenormest_matrix_power(self.A4, 2, 

473 structure=self.structure)**(1/8.) 

474 return self._d8_approx 

475 

476 @property 

477 def d10_loose(self): 

478 if self.use_exact_onenorm: 

479 return self.d10_tight 

480 if self._d10_exact is not None: 

481 return self._d10_exact 

482 else: 

483 if self._d10_approx is None: 

484 self._d10_approx = _onenormest_product((self.A4, self.A6), 

485 structure=self.structure)**(1/10.) 

486 return self._d10_approx 

487 

488 def pade3(self): 

489 b = (120., 60., 12., 1.) 

490 U = _smart_matrix_product(self.A, 

491 b[3]*self.A2 + b[1]*self.ident, 

492 structure=self.structure) 

493 V = b[2]*self.A2 + b[0]*self.ident 

494 return U, V 

495 

496 def pade5(self): 

497 b = (30240., 15120., 3360., 420., 30., 1.) 

498 U = _smart_matrix_product(self.A, 

499 b[5]*self.A4 + b[3]*self.A2 + b[1]*self.ident, 

500 structure=self.structure) 

501 V = b[4]*self.A4 + b[2]*self.A2 + b[0]*self.ident 

502 return U, V 

503 

504 def pade7(self): 

505 b = (17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.) 

506 U = _smart_matrix_product(self.A, 

507 b[7]*self.A6 + b[5]*self.A4 + b[3]*self.A2 + b[1]*self.ident, 

508 structure=self.structure) 

509 V = b[6]*self.A6 + b[4]*self.A4 + b[2]*self.A2 + b[0]*self.ident 

510 return U, V 

511 

512 def pade9(self): 

513 b = (17643225600., 8821612800., 2075673600., 302702400., 30270240., 

514 2162160., 110880., 3960., 90., 1.) 

515 U = _smart_matrix_product(self.A, 

516 (b[9]*self.A8 + b[7]*self.A6 + b[5]*self.A4 + 

517 b[3]*self.A2 + b[1]*self.ident), 

518 structure=self.structure) 

519 V = (b[8]*self.A8 + b[6]*self.A6 + b[4]*self.A4 + 

520 b[2]*self.A2 + b[0]*self.ident) 

521 return U, V 

522 

523 def pade13_scaled(self, s): 

524 b = (64764752532480000., 32382376266240000., 7771770303897600., 

525 1187353796428800., 129060195264000., 10559470521600., 

526 670442572800., 33522128640., 1323241920., 40840800., 960960., 

527 16380., 182., 1.) 

528 B = self.A * 2**-s 

529 B2 = self.A2 * 2**(-2*s) 

530 B4 = self.A4 * 2**(-4*s) 

531 B6 = self.A6 * 2**(-6*s) 

532 U2 = _smart_matrix_product(B6, 

533 b[13]*B6 + b[11]*B4 + b[9]*B2, 

534 structure=self.structure) 

535 U = _smart_matrix_product(B, 

536 (U2 + b[7]*B6 + b[5]*B4 + 

537 b[3]*B2 + b[1]*self.ident), 

538 structure=self.structure) 

539 V2 = _smart_matrix_product(B6, 

540 b[12]*B6 + b[10]*B4 + b[8]*B2, 

541 structure=self.structure) 

542 V = V2 + b[6]*B6 + b[4]*B4 + b[2]*B2 + b[0]*self.ident 

543 return U, V 

544 

545 

546def expm(A): 

547 """ 

548 Compute the matrix exponential using Pade approximation. 

549 

550 Parameters 

551 ---------- 

552 A : (M,M) array_like or sparse matrix 

553 2D Array or Matrix (sparse or dense) to be exponentiated 

554 

555 Returns 

556 ------- 

557 expA : (M,M) ndarray 

558 Matrix exponential of `A` 

559 

560 Notes 

561 ----- 

562 This is algorithm (6.1) which is a simplification of algorithm (5.1). 

563 

564 .. versionadded:: 0.12.0 

565 

566 References 

567 ---------- 

568 .. [1] Awad H. Al-Mohy and Nicholas J. Higham (2009) 

569 "A New Scaling and Squaring Algorithm for the Matrix Exponential." 

570 SIAM Journal on Matrix Analysis and Applications. 

571 31 (3). pp. 970-989. ISSN 1095-7162 

572 

573 Examples 

574 -------- 

575 >>> from scipy.sparse import csc_matrix 

576 >>> from scipy.sparse.linalg import expm 

577 >>> A = csc_matrix([[1, 0, 0], [0, 2, 0], [0, 0, 3]]) 

578 >>> A.toarray() 

579 array([[1, 0, 0], 

580 [0, 2, 0], 

581 [0, 0, 3]], dtype=int64) 

582 >>> Aexp = expm(A) 

583 >>> Aexp 

584 <3x3 sparse matrix of type '<class 'numpy.float64'>' 

585 with 3 stored elements in Compressed Sparse Column format> 

586 >>> Aexp.toarray() 

587 array([[ 2.71828183, 0. , 0. ], 

588 [ 0. , 7.3890561 , 0. ], 

589 [ 0. , 0. , 20.08553692]]) 

590 """ 

591 return _expm(A, use_exact_onenorm='auto') 

592 

593 

594def _expm(A, use_exact_onenorm): 

595 # Core of expm, separated to allow testing exact and approximate 

596 # algorithms. 

597 

598 # Avoid indiscriminate asarray() to allow sparse or other strange arrays. 

599 if isinstance(A, (list, tuple, np.matrix)): 

600 A = np.asarray(A) 

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

602 raise ValueError('expected a square matrix') 

603 

604 # gracefully handle size-0 input, 

605 # carefully handling sparse scenario 

606 if A.shape == (0, 0): 

607 out = np.zeros([0, 0], dtype=A.dtype) 

608 if issparse(A) or is_pydata_spmatrix(A): 

609 return A.__class__(out) 

610 return out 

611 

612 # Trivial case 

613 if A.shape == (1, 1): 

614 out = [[np.exp(A[0, 0])]] 

615 

616 # Avoid indiscriminate casting to ndarray to 

617 # allow for sparse or other strange arrays 

618 if issparse(A) or is_pydata_spmatrix(A): 

619 return A.__class__(out) 

620 

621 return np.array(out) 

622 

623 # Ensure input is of float type, to avoid integer overflows etc. 

624 if ((isinstance(A, np.ndarray) or issparse(A) or is_pydata_spmatrix(A)) 

625 and not np.issubdtype(A.dtype, np.inexact)): 

626 A = A.astype(float) 

627 

628 # Detect upper triangularity. 

629 structure = UPPER_TRIANGULAR if _is_upper_triangular(A) else None 

630 

631 if use_exact_onenorm == "auto": 

632 # Hardcode a matrix order threshold for exact vs. estimated one-norms. 

633 use_exact_onenorm = A.shape[0] < 200 

634 

635 # Track functions of A to help compute the matrix exponential. 

636 h = _ExpmPadeHelper( 

637 A, structure=structure, use_exact_onenorm=use_exact_onenorm) 

638 

639 # Try Pade order 3. 

640 eta_1 = max(h.d4_loose, h.d6_loose) 

641 if eta_1 < 1.495585217958292e-002 and _ell(h.A, 3) == 0: 

642 U, V = h.pade3() 

643 return _solve_P_Q(U, V, structure=structure) 

644 

645 # Try Pade order 5. 

646 eta_2 = max(h.d4_tight, h.d6_loose) 

647 if eta_2 < 2.539398330063230e-001 and _ell(h.A, 5) == 0: 

648 U, V = h.pade5() 

649 return _solve_P_Q(U, V, structure=structure) 

650 

651 # Try Pade orders 7 and 9. 

652 eta_3 = max(h.d6_tight, h.d8_loose) 

653 if eta_3 < 9.504178996162932e-001 and _ell(h.A, 7) == 0: 

654 U, V = h.pade7() 

655 return _solve_P_Q(U, V, structure=structure) 

656 if eta_3 < 2.097847961257068e+000 and _ell(h.A, 9) == 0: 

657 U, V = h.pade9() 

658 return _solve_P_Q(U, V, structure=structure) 

659 

660 # Use Pade order 13. 

661 eta_4 = max(h.d8_loose, h.d10_loose) 

662 eta_5 = min(eta_3, eta_4) 

663 theta_13 = 4.25 

664 

665 # Choose smallest s>=0 such that 2**(-s) eta_5 <= theta_13 

666 if eta_5 == 0: 

667 # Nilpotent special case 

668 s = 0 

669 else: 

670 s = max(int(np.ceil(np.log2(eta_5 / theta_13))), 0) 

671 s = s + _ell(2**-s * h.A, 13) 

672 U, V = h.pade13_scaled(s) 

673 X = _solve_P_Q(U, V, structure=structure) 

674 if structure == UPPER_TRIANGULAR: 

675 # Invoke Code Fragment 2.1. 

676 X = _fragment_2_1(X, h.A, s) 

677 else: 

678 # X = r_13(A)^(2^s) by repeated squaring. 

679 for i in range(s): 

680 X = X.dot(X) 

681 return X 

682 

683 

684def _solve_P_Q(U, V, structure=None): 

685 """ 

686 A helper function for expm_2009. 

687 

688 Parameters 

689 ---------- 

690 U : ndarray 

691 Pade numerator. 

692 V : ndarray 

693 Pade denominator. 

694 structure : str, optional 

695 A string describing the structure of both matrices `U` and `V`. 

696 Only `upper_triangular` is currently supported. 

697 

698 Notes 

699 ----- 

700 The `structure` argument is inspired by similar args 

701 for theano and cvxopt functions. 

702 

703 """ 

704 P = U + V 

705 Q = -U + V 

706 if issparse(U) or is_pydata_spmatrix(U): 

707 return spsolve(Q, P) 

708 elif structure is None: 

709 return solve(Q, P) 

710 elif structure == UPPER_TRIANGULAR: 

711 return solve_triangular(Q, P) 

712 else: 

713 raise ValueError('unsupported matrix structure: ' + str(structure)) 

714 

715 

716def _exp_sinch(a, x): 

717 """ 

718 Stably evaluate exp(a)*sinh(x)/x 

719 

720 Notes 

721 ----- 

722 The strategy of falling back to a sixth order Taylor expansion 

723 was suggested by the Spallation Neutron Source docs 

724 which was found on the internet by google search. 

725 http://www.ornl.gov/~t6p/resources/xal/javadoc/gov/sns/tools/math/ElementaryFunction.html 

726 The details of the cutoff point and the Horner-like evaluation 

727 was picked without reference to anything in particular. 

728 

729 Note that sinch is not currently implemented in scipy.special, 

730 whereas the "engineer's" definition of sinc is implemented. 

731 The implementation of sinc involves a scaling factor of pi 

732 that distinguishes it from the "mathematician's" version of sinc. 

733 

734 """ 

735 

736 # If x is small then use sixth order Taylor expansion. 

737 # How small is small? I am using the point where the relative error 

738 # of the approximation is less than 1e-14. 

739 # If x is large then directly evaluate sinh(x) / x. 

740 if abs(x) < 0.0135: 

741 x2 = x*x 

742 return np.exp(a) * (1 + (x2/6.)*(1 + (x2/20.)*(1 + (x2/42.)))) 

743 else: 

744 return (np.exp(a + x) - np.exp(a - x)) / (2*x) 

745 

746 

747def _eq_10_42(lam_1, lam_2, t_12): 

748 """ 

749 Equation (10.42) of Functions of Matrices: Theory and Computation. 

750 

751 Notes 

752 ----- 

753 This is a helper function for _fragment_2_1 of expm_2009. 

754 Equation (10.42) is on page 251 in the section on Schur algorithms. 

755 In particular, section 10.4.3 explains the Schur-Parlett algorithm. 

756 expm([[lam_1, t_12], [0, lam_1]) 

757 = 

758 [[exp(lam_1), t_12*exp((lam_1 + lam_2)/2)*sinch((lam_1 - lam_2)/2)], 

759 [0, exp(lam_2)] 

760 """ 

761 

762 # The plain formula t_12 * (exp(lam_2) - exp(lam_2)) / (lam_2 - lam_1) 

763 # apparently suffers from cancellation, according to Higham's textbook. 

764 # A nice implementation of sinch, defined as sinh(x)/x, 

765 # will apparently work around the cancellation. 

766 a = 0.5 * (lam_1 + lam_2) 

767 b = 0.5 * (lam_1 - lam_2) 

768 return t_12 * _exp_sinch(a, b) 

769 

770 

771def _fragment_2_1(X, T, s): 

772 """ 

773 A helper function for expm_2009. 

774 

775 Notes 

776 ----- 

777 The argument X is modified in-place, but this modification is not the same 

778 as the returned value of the function. 

779 This function also takes pains to do things in ways that are compatible 

780 with sparse matrices, for example by avoiding fancy indexing 

781 and by using methods of the matrices whenever possible instead of 

782 using functions of the numpy or scipy libraries themselves. 

783 

784 """ 

785 # Form X = r_m(2^-s T) 

786 # Replace diag(X) by exp(2^-s diag(T)). 

787 n = X.shape[0] 

788 diag_T = np.ravel(T.diagonal().copy()) 

789 

790 # Replace diag(X) by exp(2^-s diag(T)). 

791 scale = 2 ** -s 

792 exp_diag = np.exp(scale * diag_T) 

793 for k in range(n): 

794 X[k, k] = exp_diag[k] 

795 

796 for i in range(s-1, -1, -1): 

797 X = X.dot(X) 

798 

799 # Replace diag(X) by exp(2^-i diag(T)). 

800 scale = 2 ** -i 

801 exp_diag = np.exp(scale * diag_T) 

802 for k in range(n): 

803 X[k, k] = exp_diag[k] 

804 

805 # Replace (first) superdiagonal of X by explicit formula 

806 # for superdiagonal of exp(2^-i T) from Eq (10.42) of 

807 # the author's 2008 textbook 

808 # Functions of Matrices: Theory and Computation. 

809 for k in range(n-1): 

810 lam_1 = scale * diag_T[k] 

811 lam_2 = scale * diag_T[k+1] 

812 t_12 = scale * T[k, k+1] 

813 value = _eq_10_42(lam_1, lam_2, t_12) 

814 X[k, k+1] = value 

815 

816 # Return the updated X matrix. 

817 return X 

818 

819 

820def _ell(A, m): 

821 """ 

822 A helper function for expm_2009. 

823 

824 Parameters 

825 ---------- 

826 A : linear operator 

827 A linear operator whose norm of power we care about. 

828 m : int 

829 The power of the linear operator 

830 

831 Returns 

832 ------- 

833 value : int 

834 A value related to a bound. 

835 

836 """ 

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

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

839 

840 # The c_i are explained in (2.2) and (2.6) of the 2005 expm paper. 

841 # They are coefficients of terms of a generating function series expansion. 

842 c_i = {3: 100800., 

843 5: 10059033600., 

844 7: 4487938430976000., 

845 9: 5914384781877411840000., 

846 13: 113250775606021113483283660800000000.