Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_special_matrices.py: 10%

1import math 

2import numpy as np 

3from numpy.lib.stride_tricks import as_strided 

4 

5__all__ = ['tri', 'tril', 'triu', 'toeplitz', 'circulant', 'hankel', 

6 'hadamard', 'leslie', 'kron', 'block_diag', 'companion', 

7 'helmert', 'hilbert', 'invhilbert', 'pascal', 'invpascal', 'dft', 

8 'fiedler', 'fiedler_companion', 'convolution_matrix'] 

9 

10 

11# ----------------------------------------------------------------------------- 

12# matrix construction functions 

13# ----------------------------------------------------------------------------- 

14 

15# 

16# *Note*: tri{,u,l} is implemented in NumPy, but an important bug was fixed in 

17# 2.0.0.dev-1af2f3, the following tri{,u,l} definitions are here for backwards 

18# compatibility. 

19 

20def tri(N, M=None, k=0, dtype=None): 

21 """ 

22 Construct (N, M) matrix filled with ones at and below the kth diagonal. 

23 

24 The matrix has A[i,j] == 1 for j <= i + k 

25 

26 Parameters 

27 ---------- 

28 N : int 

29 The size of the first dimension of the matrix. 

30 M : int or None, optional 

31 The size of the second dimension of the matrix. If `M` is None, 

32 `M = N` is assumed. 

33 k : int, optional 

34 Number of subdiagonal below which matrix is filled with ones. 

35 `k` = 0 is the main diagonal, `k` < 0 subdiagonal and `k` > 0 

36 superdiagonal. 

37 dtype : dtype, optional 

38 Data type of the matrix. 

39 

40 Returns 

41 ------- 

42 tri : (N, M) ndarray 

43 Tri matrix. 

44 

45 Examples 

46 -------- 

47 >>> from scipy.linalg import tri 

48 >>> tri(3, 5, 2, dtype=int) 

49 array([[1, 1, 1, 0, 0], 

50 [1, 1, 1, 1, 0], 

51 [1, 1, 1, 1, 1]]) 

52 >>> tri(3, 5, -1, dtype=int) 

53 array([[0, 0, 0, 0, 0], 

54 [1, 0, 0, 0, 0], 

55 [1, 1, 0, 0, 0]]) 

56 

57 """ 

58 if M is None: 

59 M = N 

60 if isinstance(M, str): 

61 # pearu: any objections to remove this feature? 

62 # As tri(N,'d') is equivalent to tri(N,dtype='d') 

63 dtype = M 

64 M = N 

65 m = np.greater_equal.outer(np.arange(k, N+k), np.arange(M)) 

66 if dtype is None: 

67 return m 

68 else: 

69 return m.astype(dtype) 

70 

71 

72def tril(m, k=0): 

73 """ 

74 Make a copy of a matrix with elements above the kth diagonal zeroed. 

75 

76 Parameters 

77 ---------- 

78 m : array_like 

79 Matrix whose elements to return 

80 k : int, optional 

81 Diagonal above which to zero elements. 

82 `k` == 0 is the main diagonal, `k` < 0 subdiagonal and 

83 `k` > 0 superdiagonal. 

84 

85 Returns 

86 ------- 

87 tril : ndarray 

88 Return is the same shape and type as `m`. 

89 

90 Examples 

91 -------- 

92 >>> from scipy.linalg import tril 

93 >>> tril([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1) 

94 array([[ 0, 0, 0], 

95 [ 4, 0, 0], 

96 [ 7, 8, 0], 

97 [10, 11, 12]]) 

98 

99 """ 

100 m = np.asarray(m) 

101 out = tri(m.shape[0], m.shape[1], k=k, dtype=m.dtype.char) * m 

102 return out 

103 

104 

105def triu(m, k=0): 

106 """ 

107 Make a copy of a matrix with elements below the kth diagonal zeroed. 

108 

109 Parameters 

110 ---------- 

111 m : array_like 

112 Matrix whose elements to return 

113 k : int, optional 

114 Diagonal below which to zero elements. 

115 `k` == 0 is the main diagonal, `k` < 0 subdiagonal and 

116 `k` > 0 superdiagonal. 

117 

118 Returns 

119 ------- 

120 triu : ndarray 

121 Return matrix with zeroed elements below the kth diagonal and has 

122 same shape and type as `m`. 

123 

124 Examples 

125 -------- 

126 >>> from scipy.linalg import triu 

127 >>> triu([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1) 

128 array([[ 1, 2, 3], 

129 [ 4, 5, 6], 

130 [ 0, 8, 9], 

131 [ 0, 0, 12]]) 

132 

133 """ 

134 m = np.asarray(m) 

135 out = (1 - tri(m.shape[0], m.shape[1], k - 1, m.dtype.char)) * m 

136 return out 

137 

138 

139def toeplitz(c, r=None): 

140 """ 

141 Construct a Toeplitz matrix. 

142 

143 The Toeplitz matrix has constant diagonals, with c as its first column 

144 and r as its first row. If r is not given, ``r == conjugate(c)`` is 

145 assumed. 

146 

147 Parameters 

148 ---------- 

149 c : array_like 

150 First column of the matrix. Whatever the actual shape of `c`, it 

151 will be converted to a 1-D array. 

152 r : array_like, optional 

153 First row of the matrix. If None, ``r = conjugate(c)`` is assumed; 

154 in this case, if c[0] is real, the result is a Hermitian matrix. 

155 r[0] is ignored; the first row of the returned matrix is 

156 ``[c[0], r[1:]]``. Whatever the actual shape of `r`, it will be 

157 converted to a 1-D array. 

158 

159 Returns 

160 ------- 

161 A : (len(c), len(r)) ndarray 

162 The Toeplitz matrix. Dtype is the same as ``(c[0] + r[0]).dtype``. 

163 

164 See Also 

165 -------- 

166 circulant : circulant matrix 

167 hankel : Hankel matrix 

168 solve_toeplitz : Solve a Toeplitz system. 

169 

170 Notes 

171 ----- 

172 The behavior when `c` or `r` is a scalar, or when `c` is complex and 

173 `r` is None, was changed in version 0.8.0. The behavior in previous 

174 versions was undocumented and is no longer supported. 

175 

176 Examples 

177 -------- 

178 >>> from scipy.linalg import toeplitz 

179 >>> toeplitz([1,2,3], [1,4,5,6]) 

180 array([[1, 4, 5, 6], 

181 [2, 1, 4, 5], 

182 [3, 2, 1, 4]]) 

183 >>> toeplitz([1.0, 2+3j, 4-1j]) 

184 array([[ 1.+0.j, 2.-3.j, 4.+1.j], 

185 [ 2.+3.j, 1.+0.j, 2.-3.j], 

186 [ 4.-1.j, 2.+3.j, 1.+0.j]]) 

187 

188 """ 

189 c = np.asarray(c).ravel() 

190 if r is None: 

191 r = c.conjugate() 

192 else: 

193 r = np.asarray(r).ravel() 

194 # Form a 1-D array containing a reversed c followed by r[1:] that could be 

195 # strided to give us toeplitz matrix. 

196 vals = np.concatenate((c[::-1], r[1:])) 

197 out_shp = len(c), len(r) 

198 n = vals.strides[0] 

199 return as_strided(vals[len(c)-1:], shape=out_shp, strides=(-n, n)).copy() 

200 

201 

202def circulant(c): 

203 """ 

204 Construct a circulant matrix. 

205 

206 Parameters 

207 ---------- 

208 c : (N,) array_like 

209 1-D array, the first column of the matrix. 

210 

211 Returns 

212 ------- 

213 A : (N, N) ndarray 

214 A circulant matrix whose first column is `c`. 

215 

216 See Also 

217 -------- 

218 toeplitz : Toeplitz matrix 

219 hankel : Hankel matrix 

220 solve_circulant : Solve a circulant system. 

221 

222 Notes 

223 ----- 

224 .. versionadded:: 0.8.0 

225 

226 Examples 

227 -------- 

228 >>> from scipy.linalg import circulant 

229 >>> circulant([1, 2, 3]) 

230 array([[1, 3, 2], 

231 [2, 1, 3], 

232 [3, 2, 1]]) 

233 

234 """ 

235 c = np.asarray(c).ravel() 

236 # Form an extended array that could be strided to give circulant version 

237 c_ext = np.concatenate((c[::-1], c[:0:-1])) 

238 L = len(c) 

239 n = c_ext.strides[0] 

240 return as_strided(c_ext[L-1:], shape=(L, L), strides=(-n, n)).copy() 

241 

242 

243def hankel(c, r=None): 

244 """ 

245 Construct a Hankel matrix. 

246 

247 The Hankel matrix has constant anti-diagonals, with `c` as its 

248 first column and `r` as its last row. If `r` is not given, then 

249 `r = zeros_like(c)` is assumed. 

250 

251 Parameters 

252 ---------- 

253 c : array_like 

254 First column of the matrix. Whatever the actual shape of `c`, it 

255 will be converted to a 1-D array. 

256 r : array_like, optional 

257 Last row of the matrix. If None, ``r = zeros_like(c)`` is assumed. 

258 r[0] is ignored; the last row of the returned matrix is 

259 ``[c[-1], r[1:]]``. Whatever the actual shape of `r`, it will be 

260 converted to a 1-D array. 

261 

262 Returns 

263 ------- 

264 A : (len(c), len(r)) ndarray 

265 The Hankel matrix. Dtype is the same as ``(c[0] + r[0]).dtype``. 

266 

267 See Also 

268 -------- 

269 toeplitz : Toeplitz matrix 

270 circulant : circulant matrix 

271 

272 Examples 

273 -------- 

274 >>> from scipy.linalg import hankel 

275 >>> hankel([1, 17, 99]) 

276 array([[ 1, 17, 99], 

277 [17, 99, 0], 

278 [99, 0, 0]]) 

279 >>> hankel([1,2,3,4], [4,7,7,8,9]) 

280 array([[1, 2, 3, 4, 7], 

281 [2, 3, 4, 7, 7], 

282 [3, 4, 7, 7, 8], 

283 [4, 7, 7, 8, 9]]) 

284 

285 """ 

286 c = np.asarray(c).ravel() 

287 if r is None: 

288 r = np.zeros_like(c) 

289 else: 

290 r = np.asarray(r).ravel() 

291 # Form a 1-D array of values to be used in the matrix, containing `c` 

292 # followed by r[1:]. 

293 vals = np.concatenate((c, r[1:])) 

294 # Stride on concatenated array to get hankel matrix 

295 out_shp = len(c), len(r) 

296 n = vals.strides[0] 

297 return as_strided(vals, shape=out_shp, strides=(n, n)).copy() 

298 

299 

300def hadamard(n, dtype=int): 

301 """ 

302 Construct an Hadamard matrix. 

303 

304 Constructs an n-by-n Hadamard matrix, using Sylvester's 

305 construction. `n` must be a power of 2. 

306 

307 Parameters 

308 ---------- 

309 n : int 

310 The order of the matrix. `n` must be a power of 2. 

311 dtype : dtype, optional 

312 The data type of the array to be constructed. 

313 

314 Returns 

315 ------- 

316 H : (n, n) ndarray 

317 The Hadamard matrix. 

318 

319 Notes 

320 ----- 

321 .. versionadded:: 0.8.0 

322 

323 Examples 

324 -------- 

325 >>> from scipy.linalg import hadamard 

326 >>> hadamard(2, dtype=complex) 

327 array([[ 1.+0.j, 1.+0.j], 

328 [ 1.+0.j, -1.-0.j]]) 

329 >>> hadamard(4) 

330 array([[ 1, 1, 1, 1], 

331 [ 1, -1, 1, -1], 

332 [ 1, 1, -1, -1], 

333 [ 1, -1, -1, 1]]) 

334 

335 """ 

336 

337 # This function is a slightly modified version of the 

338 # function contributed by Ivo in ticket #675. 

339 

340 if n < 1: 

341 lg2 = 0 

342 else: 

343 lg2 = int(math.log(n, 2)) 

344 if 2 ** lg2 != n: 

345 raise ValueError("n must be an positive integer, and n must be " 

346 "a power of 2") 

347 

348 H = np.array([[1]], dtype=dtype) 

349 

350 # Sylvester's construction 

351 for i in range(0, lg2): 

352 H = np.vstack((np.hstack((H, H)), np.hstack((H, -H)))) 

353 

354 return H 

355 

356 

357def leslie(f, s): 

358 """ 

359 Create a Leslie matrix. 

360 

361 Given the length n array of fecundity coefficients `f` and the length 

362 n-1 array of survival coefficients `s`, return the associated Leslie 

363 matrix. 

364 

365 Parameters 

366 ---------- 

367 f : (N,) array_like 

368 The "fecundity" coefficients. 

369 s : (N-1,) array_like 

370 The "survival" coefficients, has to be 1-D. The length of `s` 

371 must be one less than the length of `f`, and it must be at least 1. 

372 

373 Returns 

374 ------- 

375 L : (N, N) ndarray 

376 The array is zero except for the first row, 

377 which is `f`, and the first sub-diagonal, which is `s`. 

378 The data-type of the array will be the data-type of ``f[0]+s[0]``. 

379 

380 Notes 

381 ----- 

382 .. versionadded:: 0.8.0 

383 

384 The Leslie matrix is used to model discrete-time, age-structured 

385 population growth [1]_ [2]_. In a population with `n` age classes, two sets 

386 of parameters define a Leslie matrix: the `n` "fecundity coefficients", 

387 which give the number of offspring per-capita produced by each age 

388 class, and the `n` - 1 "survival coefficients", which give the 

389 per-capita survival rate of each age class. 

390 

391 References 

392 ---------- 

393 .. [1] P. H. Leslie, On the use of matrices in certain population 

394 mathematics, Biometrika, Vol. 33, No. 3, 183--212 (Nov. 1945) 

395 .. [2] P. H. Leslie, Some further notes on the use of matrices in 

396 population mathematics, Biometrika, Vol. 35, No. 3/4, 213--245 

397 (Dec. 1948) 

398 

399 Examples 

400 -------- 

401 >>> from scipy.linalg import leslie 

402 >>> leslie([0.1, 2.0, 1.0, 0.1], [0.2, 0.8, 0.7]) 

403 array([[ 0.1, 2. , 1. , 0.1], 

404 [ 0.2, 0. , 0. , 0. ], 

405 [ 0. , 0.8, 0. , 0. ], 

406 [ 0. , 0. , 0.7, 0. ]]) 

407 

408 """ 

409 f = np.atleast_1d(f) 

410 s = np.atleast_1d(s) 

411 if f.ndim != 1: 

412 raise ValueError("Incorrect shape for f. f must be 1D") 

413 if s.ndim != 1: 

414 raise ValueError("Incorrect shape for s. s must be 1D") 

415 if f.size != s.size + 1: 

416 raise ValueError("Incorrect lengths for f and s. The length" 

417 " of s must be one less than the length of f.") 

418 if s.size == 0: 

419 raise ValueError("The length of s must be at least 1.") 

420 

421 tmp = f[0] + s[0] 

422 n = f.size 

423 a = np.zeros((n, n), dtype=tmp.dtype) 

424 a[0] = f 

425 a[list(range(1, n)), list(range(0, n - 1))] = s 

426 return a 

427 

428 

429def kron(a, b): 

430 """ 

431 Kronecker product. 

432 

433 The result is the block matrix:: 

434 

435 a[0,0]*b a[0,1]*b ... a[0,-1]*b 

436 a[1,0]*b a[1,1]*b ... a[1,-1]*b 

437 ... 

438 a[-1,0]*b a[-1,1]*b ... a[-1,-1]*b 

439 

440 Parameters 

441 ---------- 

442 a : (M, N) ndarray 

443 Input array 

444 b : (P, Q) ndarray 

445 Input array 

446 

447 Returns 

448 ------- 

449 A : (M*P, N*Q) ndarray 

450 Kronecker product of `a` and `b`. 

451 

452 Examples 

453 -------- 

454 >>> from numpy import array 

455 >>> from scipy.linalg import kron 

456 >>> kron(array([[1,2],[3,4]]), array([[1,1,1]])) 

457 array([[1, 1, 1, 2, 2, 2], 

458 [3, 3, 3, 4, 4, 4]]) 

459 

460 """ 

461 if not a.flags['CONTIGUOUS']: 

462 a = np.reshape(a, a.shape) 

463 if not b.flags['CONTIGUOUS']: 

464 b = np.reshape(b, b.shape) 

465 o = np.outer(a, b) 

466 o = o.reshape(a.shape + b.shape) 

467 return np.concatenate(np.concatenate(o, axis=1), axis=1) 

468 

469 

470def block_diag(*arrs): 

471 """ 

472 Create a block diagonal matrix from provided arrays. 

473 

474 Given the inputs `A`, `B` and `C`, the output will have these 

475 arrays arranged on the diagonal:: 

476 

477 [[A, 0, 0], 

478 [0, B, 0], 

479 [0, 0, C]] 

480 

481 Parameters 

482 ---------- 

483 A, B, C, ... : array_like, up to 2-D 

484 Input arrays. A 1-D array or array_like sequence of length `n` is 

485 treated as a 2-D array with shape ``(1,n)``. 

486 

487 Returns 

488 ------- 

489 D : ndarray 

490 Array with `A`, `B`, `C`, ... on the diagonal. `D` has the 

491 same dtype as `A`. 

492 

493 Notes 

494 ----- 

495 If all the input arrays are square, the output is known as a 

496 block diagonal matrix. 

497 

498 Empty sequences (i.e., array-likes of zero size) will not be ignored. 

499 Noteworthy, both [] and [[]] are treated as matrices with shape ``(1,0)``. 

500 

501 Examples 

502 -------- 

503 >>> import numpy as np 

504 >>> from scipy.linalg import block_diag 

505 >>> A = [[1, 0], 

506 ... [0, 1]] 

507 >>> B = [[3, 4, 5], 

508 ... [6, 7, 8]] 

509 >>> C = [[7]] 

510 >>> P = np.zeros((2, 0), dtype='int32') 

511 >>> block_diag(A, B, C) 

512 array([[1, 0, 0, 0, 0, 0], 

513 [0, 1, 0, 0, 0, 0], 

514 [0, 0, 3, 4, 5, 0], 

515 [0, 0, 6, 7, 8, 0], 

516 [0, 0, 0, 0, 0, 7]]) 

517 >>> block_diag(A, P, B, C) 

518 array([[1, 0, 0, 0, 0, 0], 

519 [0, 1, 0, 0, 0, 0], 

520 [0, 0, 0, 0, 0, 0], 

521 [0, 0, 0, 0, 0, 0], 

522 [0, 0, 3, 4, 5, 0], 

523 [0, 0, 6, 7, 8, 0], 

524 [0, 0, 0, 0, 0, 7]]) 

525 >>> block_diag(1.0, [2, 3], [[4, 5], [6, 7]]) 

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

527 [ 0., 2., 3., 0., 0.], 

528 [ 0., 0., 0., 4., 5.], 

529 [ 0., 0., 0., 6., 7.]]) 

530 

531 """ 

532 if arrs == (): 

533 arrs = ([],) 

534 arrs = [np.atleast_2d(a) for a in arrs] 

535 

536 bad_args = [k for k in range(len(arrs)) if arrs[k].ndim > 2] 

537 if bad_args: 

538 raise ValueError("arguments in the following positions have dimension " 

539 "greater than 2: %s" % bad_args) 

540 

541 shapes = np.array([a.shape for a in arrs]) 

542 out_dtype = np.result_type(*[arr.dtype for arr in arrs]) 

543 out = np.zeros(np.sum(shapes, axis=0), dtype=out_dtype) 

544 

545 r, c = 0, 0 

546 for i, (rr, cc) in enumerate(shapes): 

547 out[r:r + rr, c:c + cc] = arrs[i] 

548 r += rr 

549 c += cc 

550 return out 

551 

552 

553def companion(a): 

554 """ 

555 Create a companion matrix. 

556 

557 Create the companion matrix [1]_ associated with the polynomial whose 

558 coefficients are given in `a`. 

559 

560 Parameters 

561 ---------- 

562 a : (N,) array_like 

563 1-D array of polynomial coefficients. The length of `a` must be 

564 at least two, and ``a[0]`` must not be zero. 

565 

566 Returns 

567 ------- 

568 c : (N-1, N-1) ndarray 

569 The first row of `c` is ``-a[1:]/a[0]``, and the first 

570 sub-diagonal is all ones. The data-type of the array is the same 

571 as the data-type of ``1.0*a[0]``. 

572 

573 Raises 

574 ------ 

575 ValueError 

576 If any of the following are true: a) ``a.ndim != 1``; 

577 b) ``a.size < 2``; c) ``a[0] == 0``. 

578 

579 Notes 

580 ----- 

581 .. versionadded:: 0.8.0 

582 

583 References 

584 ---------- 

585 .. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK: 

586 Cambridge University Press, 1999, pp. 146-7. 

587 

588 Examples 

589 -------- 

590 >>> from scipy.linalg import companion 

591 >>> companion([1, -10, 31, -30]) 

592 array([[ 10., -31., 30.], 

593 [ 1., 0., 0.], 

594 [ 0., 1., 0.]]) 

595 

596 """ 

597 a = np.atleast_1d(a) 

598 

599 if a.ndim != 1: 

600 raise ValueError("Incorrect shape for `a`. `a` must be " 

601 "one-dimensional.") 

602 

603 if a.size < 2: 

604 raise ValueError("The length of `a` must be at least 2.") 

605 

606 if a[0] == 0: 

607 raise ValueError("The first coefficient in `a` must not be zero.") 

608 

609 first_row = -a[1:] / (1.0 * a[0]) 

610 n = a.size 

611 c = np.zeros((n - 1, n - 1), dtype=first_row.dtype) 

612 c[0] = first_row 

613 c[list(range(1, n - 1)), list(range(0, n - 2))] = 1 

614 return c 

615 

616 

617def helmert(n, full=False): 

618 """ 

619 Create an Helmert matrix of order `n`. 

620 

621 This has applications in statistics, compositional or simplicial analysis, 

622 and in Aitchison geometry. 

623 

624 Parameters 

625 ---------- 

626 n : int 

627 The size of the array to create. 

628 full : bool, optional 

629 If True the (n, n) ndarray will be returned. 

630 Otherwise the submatrix that does not include the first 

631 row will be returned. 

632 Default: False. 

633 

634 Returns 

635 ------- 

636 M : ndarray 

637 The Helmert matrix. 

638 The shape is (n, n) or (n-1, n) depending on the `full` argument. 

639 

640 Examples 

641 -------- 

642 >>> from scipy.linalg import helmert 

643 >>> helmert(5, full=True) 

644 array([[ 0.4472136 , 0.4472136 , 0.4472136 , 0.4472136 , 0.4472136 ], 

645 [ 0.70710678, -0.70710678, 0. , 0. , 0. ], 

646 [ 0.40824829, 0.40824829, -0.81649658, 0. , 0. ], 

647 [ 0.28867513, 0.28867513, 0.28867513, -0.8660254 , 0. ], 

648 [ 0.2236068 , 0.2236068 , 0.2236068 , 0.2236068 , -0.89442719]]) 

649 

650 """ 

651 H = np.tril(np.ones((n, n)), -1) - np.diag(np.arange(n)) 

652 d = np.arange(n) * np.arange(1, n+1) 

653 H[0] = 1 

654 d[0] = n 

655 H_full = H / np.sqrt(d)[:, np.newaxis] 

656 if full: 

657 return H_full 

658 else: 

659 return H_full[1:] 

660 

661 

662def hilbert(n): 

663 """ 

664 Create a Hilbert matrix of order `n`. 

665 

666 Returns the `n` by `n` array with entries `h[i,j] = 1 / (i + j + 1)`. 

667 

668 Parameters 

669 ---------- 

670 n : int 

671 The size of the array to create. 

672 

673 Returns 

674 ------- 

675 h : (n, n) ndarray 

676 The Hilbert matrix. 

677 

678 See Also 

679 -------- 

680 invhilbert : Compute the inverse of a Hilbert matrix. 

681 

682 Notes 

683 ----- 

684 .. versionadded:: 0.10.0 

685 

686 Examples 

687 -------- 

688 >>> from scipy.linalg import hilbert 

689 >>> hilbert(3) 

690 array([[ 1. , 0.5 , 0.33333333], 

691 [ 0.5 , 0.33333333, 0.25 ], 

692 [ 0.33333333, 0.25 , 0.2 ]]) 

693 

694 """ 

695 values = 1.0 / (1.0 + np.arange(2 * n - 1)) 

696 h = hankel(values[:n], r=values[n - 1:]) 

697 return h 

698 

699 

700def invhilbert(n, exact=False): 

701 """ 

702 Compute the inverse of the Hilbert matrix of order `n`. 

703 

704 The entries in the inverse of a Hilbert matrix are integers. When `n` 

705 is greater than 14, some entries in the inverse exceed the upper limit 

706 of 64 bit integers. The `exact` argument provides two options for 

707 dealing with these large integers. 

708 

709 Parameters 

710 ---------- 

711 n : int 

712 The order of the Hilbert matrix. 

713 exact : bool, optional 

714 If False, the data type of the array that is returned is np.float64, 

715 and the array is an approximation of the inverse. 

716 If True, the array is the exact integer inverse array. To represent 

717 the exact inverse when n > 14, the returned array is an object array 

718 of long integers. For n <= 14, the exact inverse is returned as an 

719 array with data type np.int64. 

720 

721 Returns 

722 ------- 

723 invh : (n, n) ndarray 

724 The data type of the array is np.float64 if `exact` is False. 

725 If `exact` is True, the data type is either np.int64 (for n <= 14) 

726 or object (for n > 14). In the latter case, the objects in the 

727 array will be long integers. 

728 

729 See Also 

730 -------- 

731 hilbert : Create a Hilbert matrix. 

732 

733 Notes 

734 ----- 

735 .. versionadded:: 0.10.0 

736 

737 Examples 

738 -------- 

739 >>> from scipy.linalg import invhilbert 

740 >>> invhilbert(4) 

741 array([[ 16., -120., 240., -140.], 

742 [ -120., 1200., -2700., 1680.], 

743 [ 240., -2700., 6480., -4200.], 

744 [ -140., 1680., -4200., 2800.]]) 

745 >>> invhilbert(4, exact=True) 

746 array([[ 16, -120, 240, -140], 

747 [ -120, 1200, -2700, 1680], 

748 [ 240, -2700, 6480, -4200], 

749 [ -140, 1680, -4200, 2800]], dtype=int64) 

750 >>> invhilbert(16)[7,7] 

751 4.2475099528537506e+19 

752 >>> invhilbert(16, exact=True)[7,7] 

753 42475099528537378560 

754 

755 """ 

756 from scipy.special import comb 

757 if exact: 

758 if n > 14: 

759 dtype = object 

760 else: 

761 dtype = np.int64 

762 else: 

763 dtype = np.float64 

764 invh = np.empty((n, n), dtype=dtype) 

765 for i in range(n): 

766 for j in range(0, i + 1): 

767 s = i + j 

768 invh[i, j] = ((-1) ** s * (s + 1) * 

769 comb(n + i, n - j - 1, exact) * 

770 comb(n + j, n - i - 1, exact) * 

771 comb(s, i, exact) ** 2) 

772 if i != j: 

773 invh[j, i] = invh[i, j] 

774 return invh 

775 

776 

777def pascal(n, kind='symmetric', exact=True): 

778 """ 

779 Returns the n x n Pascal matrix. 

780 

781 The Pascal matrix is a matrix containing the binomial coefficients as 

782 its elements. 

783 

784 Parameters 

785 ---------- 

786 n : int 

787 The size of the matrix to create; that is, the result is an n x n 

788 matrix. 

789 kind : str, optional 

790 Must be one of 'symmetric', 'lower', or 'upper'. 

791 Default is 'symmetric'. 

792 exact : bool, optional 

793 If `exact` is True, the result is either an array of type 

794 numpy.uint64 (if n < 35) or an object array of Python long integers. 

795 If `exact` is False, the coefficients in the matrix are computed using 

796 `scipy.special.comb` with `exact=False`. The result will be a floating 

797 point array, and the values in the array will not be the exact 

798 coefficients, but this version is much faster than `exact=True`. 

799 

800 Returns 

801 ------- 

802 p : (n, n) ndarray 

803 The Pascal matrix. 

804 

805 See Also 

806 -------- 

807 invpascal 

808 

809 Notes 

810 ----- 

811 See https://en.wikipedia.org/wiki/Pascal_matrix for more information 

812 about Pascal matrices. 

813 

814 .. versionadded:: 0.11.0 

815 

816 Examples 

817 -------- 

818 >>> from scipy.linalg import pascal 

819 >>> pascal(4) 

820 array([[ 1, 1, 1, 1], 

821 [ 1, 2, 3, 4], 

822 [ 1, 3, 6, 10], 

823 [ 1, 4, 10, 20]], dtype=uint64) 

824 >>> pascal(4, kind='lower') 

825 array([[1, 0, 0, 0], 

826 [1, 1, 0, 0], 

827 [1, 2, 1, 0], 

828 [1, 3, 3, 1]], dtype=uint64) 

829 >>> pascal(50)[-1, -1] 

830 25477612258980856902730428600 

831 >>> from scipy.special import comb 

832 >>> comb(98, 49, exact=True) 

833 25477612258980856902730428600 

834 

835 """ 

836 

837 from scipy.special import comb 

838 if kind not in ['symmetric', 'lower', 'upper']: 

839 raise ValueError("kind must be 'symmetric', 'lower', or 'upper'") 

840 

841 if exact: 

842 if n >= 35: 

843 L_n = np.empty((n, n), dtype=object) 

844 L_n.fill(0) 

845 else: 

846 L_n = np.zeros((n, n), dtype=np.uint64) 

847 for i in range(n): 

848 for j in range(i + 1): 

849 L_n[i, j] = comb(i, j, exact=True) 

850 else: 

851 L_n = comb(*np.ogrid[:n, :n]) 

852 

853 if kind == 'lower': 

854 p = L_n 

855 elif kind == 'upper': 

856 p = L_n.T 

857 else: 

858 p = np.dot(L_n, L_n.T) 

859 

860 return p 

861 

862 

863def invpascal(n, kind='symmetric', exact=True): 

864 """ 

865 Returns the inverse of the n x n Pascal matrix. 

866 

867 The Pascal matrix is a matrix containing the binomial coefficients as 

868 its elements. 

869 

870 Parameters 

871 ---------- 

872 n : int 

873 The size of the matrix to create; that is, the result is an n x n 

874 matrix. 

875 kind : str, optional 

876 Must be one of 'symmetric', 'lower', or 'upper'. 

877 Default is 'symmetric'. 

878 exact : bool, optional 

879 If `exact` is True, the result is either an array of type 

880 ``numpy.int64`` (if `n` <= 35) or an object array of Python integers. 

881 If `exact` is False, the coefficients in the matrix are computed using 

882 `scipy.special.comb` with `exact=False`. The result will be a floating 

883 point array, and for large `n`, the values in the array will not be the 

884 exact coefficients. 

885 

886 Returns 

887 ------- 

888 invp : (n, n) ndarray 

889 The inverse of the Pascal matrix. 

890 

891 See Also 

892 -------- 

893 pascal 

894 

895 Notes 

896 ----- 

897 

898 .. versionadded:: 0.16.0 

899 

900 References 

901 ---------- 

902 .. [1] "Pascal matrix", https://en.wikipedia.org/wiki/Pascal_matrix 

903 .. [2] Cohen, A. M., "The inverse of a Pascal matrix", Mathematical 

904 Gazette, 59(408), pp. 111-112, 1975. 

905 

906 Examples 

907 -------- 

908 >>> from scipy.linalg import invpascal, pascal 

909 >>> invp = invpascal(5) 

910 >>> invp 

911 array([[ 5, -10, 10, -5, 1], 

912 [-10, 30, -35, 19, -4], 

913 [ 10, -35, 46, -27, 6], 

914 [ -5, 19, -27, 17, -4], 

915 [ 1, -4, 6, -4, 1]]) 

916 

917 >>> p = pascal(5) 

918 >>> p.dot(invp) 

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

920 [ 0., 1., 0., 0., 0.], 

921 [ 0., 0., 1., 0., 0.], 

922 [ 0., 0., 0., 1., 0.], 

923 [ 0., 0., 0., 0., 1.]]) 

924 

925 An example of the use of `kind` and `exact`: 

926 

927 >>> invpascal(5, kind='lower', exact=False) 

928 array([[ 1., -0., 0., -0., 0.], 

929 [-1., 1., -0., 0., -0.], 

930 [ 1., -2., 1., -0., 0.], 

931 [-1., 3., -3., 1., -0.], 

932 [ 1., -4., 6., -4., 1.]]) 

933 

934 """ 

935 from scipy.special import comb 

936 

937 if kind not in ['symmetric', 'lower', 'upper']: 

938 raise ValueError("'kind' must be 'symmetric', 'lower' or 'upper'.") 

939 

940 if kind == 'symmetric': 

941 if exact: 

942 if n > 34: 

943 dt = object 

944 else: 

945 dt = np.int64 

946 else: 

947 dt = np.float64 

948 invp = np.empty((n, n), dtype=dt) 

949 for i in range(n): 

950 for j in range(0, i + 1): 

951 v = 0 

952 for k in range(n - i): 

953 v += comb(i + k, k, exact=exact) * comb(i + k, i + k - j, 

954 exact=exact) 

955 invp[i, j] = (-1)**(i - j) * v 

956 if i != j: 

957 invp[j, i] = invp[i, j] 

958 else: 

959 # For the 'lower' and 'upper' cases, we computer the inverse by 

960 # changing the sign of every other diagonal of the pascal matrix. 

961 invp = pascal(n, kind=kind, exact=exact) 

962 if invp.dtype == np.uint64: 

963 # This cast from np.uint64 to int64 OK, because if `kind` is not 

964 # "symmetric", the values in invp are all much less than 2**63. 

965 invp = invp.view(np.int64) 

966 

967 # The toeplitz matrix has alternating bands of 1 and -1. 

968 invp *= toeplitz((-1)**np.arange(n)).astype(invp.dtype) 

969 

970 return invp 

971 

972 

973def dft(n, scale=None): 

974 """ 

975 Discrete Fourier transform matrix. 

976 

977 Create the matrix that computes the discrete Fourier transform of a 

978 sequence [1]_. The nth primitive root of unity used to generate the 

979 matrix is exp(-2*pi*i/n), where i = sqrt(-1). 

980 

981 Parameters 

982 ---------- 

983 n : int 

984 Size the matrix to create. 

985 scale : str, optional 

986 Must be None, 'sqrtn', or 'n'. 

987 If `scale` is 'sqrtn', the matrix is divided by `sqrt(n)`. 

988 If `scale` is 'n', the matrix is divided by `n`. 

989 If `scale` is None (the default), the matrix is not normalized, and the 

990 return value is simply the Vandermonde matrix of the roots of unity. 

991 

992 Returns 

993 ------- 

994 m : (n, n) ndarray 

995 The DFT matrix. 

996 

997 Notes 

998 ----- 

999 When `scale` is None, multiplying a vector by the matrix returned by 

1000 `dft` is mathematically equivalent to (but much less efficient than) 

1001 the calculation performed by `scipy.fft.fft`. 

1002 

1003 .. versionadded:: 0.14.0 

1004 

1005 References 

1006 ---------- 

1007 .. [1] "DFT matrix", https://en.wikipedia.org/wiki/DFT_matrix 

1008 

1009 Examples 

1010 -------- 

1011 >>> import numpy as np 

1012 >>> from scipy.linalg import dft 

1013 >>> np.set_printoptions(precision=2, suppress=True) # for compact output 

1014 >>> m = dft(5) 

1015 >>> m 

1016 array([[ 1. +0.j , 1. +0.j , 1. +0.j , 1. +0.j , 1. +0.j ], 

1017 [ 1. +0.j , 0.31-0.95j, -0.81-0.59j, -0.81+0.59j, 0.31+0.95j], 

1018 [ 1. +0.j , -0.81-0.59j, 0.31+0.95j, 0.31-0.95j, -0.81+0.59j], 

1019 [ 1. +0.j , -0.81+0.59j, 0.31-0.95j, 0.31+0.95j, -0.81-0.59j], 

1020 [ 1. +0.j , 0.31+0.95j, -0.81+0.59j, -0.81-0.59j, 0.31-0.95j]]) 

1021 >>> x = np.array([1, 2, 3, 0, 3]) 

1022 >>> m @ x # Compute the DFT of x 

1023 array([ 9. +0.j , 0.12-0.81j, -2.12+3.44j, -2.12-3.44j, 0.12+0.81j]) 

1024 

1025 Verify that ``m @ x`` is the same as ``fft(x)``. 

1026 

1027 >>> from scipy.fft import fft 

1028 >>> fft(x) # Same result as m @ x 

1029 array([ 9. +0.j , 0.12-0.81j, -2.12+3.44j, -2.12-3.44j, 0.12+0.81j]) 

1030 """ 

1031 if scale not in [None, 'sqrtn', 'n']: 

1032 raise ValueError("scale must be None, 'sqrtn', or 'n'; " 

1033 "%r is not valid." % (scale,)) 

1034 

1035 omegas = np.exp(-2j * np.pi * np.arange(n) / n).reshape(-1, 1) 

1036 m = omegas ** np.arange(n) 

1037 if scale == 'sqrtn': 

1038 m /= math.sqrt(n) 

1039 elif scale == 'n': 

1040 m /= n 

1041 return m 

1042 

1043 

1044def fiedler(a): 

1045 """Returns a symmetric Fiedler matrix 

1046 

1047 Given an sequence of numbers `a`, Fiedler matrices have the structure 

1048 ``F[i, j] = np.abs(a[i] - a[j])``, and hence zero diagonals and nonnegative 

1049 entries. A Fiedler matrix has a dominant positive eigenvalue and other 

1050 eigenvalues are negative. Although not valid generally, for certain inputs, 

1051 the inverse and the determinant can be derived explicitly as given in [1]_. 

1052 

1053 Parameters 

1054 ---------- 

1055 a : (n,) array_like 

1056 coefficient array 

1057 

1058 Returns 

1059 ------- 

1060 F : (n, n) ndarray