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

1import math 

2 

3import numpy as np 

4from numpy.lib.stride_tricks import as_strided 

5 

6__all__ = ['toeplitz', 'circulant', 'hankel', 

7 'hadamard', 'leslie', 'kron', 'block_diag', 'companion', 

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

9 'fiedler', 'fiedler_companion', 'convolution_matrix'] 

10 

11 

12# ----------------------------------------------------------------------------- 

13# matrix construction functions 

14# ----------------------------------------------------------------------------- 

15 

16 

17def toeplitz(c, r=None): 

18 """ 

19 Construct a Toeplitz matrix. 

20 

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

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

23 assumed. 

24 

25 Parameters 

26 ---------- 

27 c : array_like 

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

29 will be converted to a 1-D array. 

30 r : array_like, optional 

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

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

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

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

35 converted to a 1-D array. 

36 

37 Returns 

38 ------- 

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

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

41 

42 See Also 

43 -------- 

44 circulant : circulant matrix 

45 hankel : Hankel matrix 

46 solve_toeplitz : Solve a Toeplitz system. 

47 

48 Notes 

49 ----- 

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

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

52 versions was undocumented and is no longer supported. 

53 

54 Examples 

55 -------- 

56 >>> from scipy.linalg import toeplitz 

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

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

59 [2, 1, 4, 5], 

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

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

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

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

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

65 

66 """ 

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

68 if r is None: 

69 r = c.conjugate() 

70 else: 

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

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

73 # strided to give us toeplitz matrix. 

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

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

76 n = vals.strides[0] 

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

78 

79 

80def circulant(c): 

81 """ 

82 Construct a circulant matrix. 

83 

84 Parameters 

85 ---------- 

86 c : (N,) array_like 

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

88 

89 Returns 

90 ------- 

91 A : (N, N) ndarray 

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

93 

94 See Also 

95 -------- 

96 toeplitz : Toeplitz matrix 

97 hankel : Hankel matrix 

98 solve_circulant : Solve a circulant system. 

99 

100 Notes 

101 ----- 

102 .. versionadded:: 0.8.0 

103 

104 Examples 

105 -------- 

106 >>> from scipy.linalg import circulant 

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

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

109 [2, 1, 3], 

110 [3, 2, 1]]) 

111 

112 """ 

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

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

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

116 L = len(c) 

117 n = c_ext.strides[0] 

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

119 

120 

121def hankel(c, r=None): 

122 """ 

123 Construct a Hankel matrix. 

124 

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

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

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

128 

129 Parameters 

130 ---------- 

131 c : array_like 

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

133 will be converted to a 1-D array. 

134 r : array_like, optional 

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

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

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

138 converted to a 1-D array. 

139 

140 Returns 

141 ------- 

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

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

144 

145 See Also 

146 -------- 

147 toeplitz : Toeplitz matrix 

148 circulant : circulant matrix 

149 

150 Examples 

151 -------- 

152 >>> from scipy.linalg import hankel 

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

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

155 [17, 99, 0], 

156 [99, 0, 0]]) 

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

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

159 [2, 3, 4, 7, 7], 

160 [3, 4, 7, 7, 8], 

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

162 

163 """ 

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

165 if r is None: 

166 r = np.zeros_like(c) 

167 else: 

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

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

170 # followed by r[1:]. 

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

172 # Stride on concatenated array to get hankel matrix 

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

174 n = vals.strides[0] 

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

176 

177 

178def hadamard(n, dtype=int): 

179 """ 

180 Construct an Hadamard matrix. 

181 

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

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

184 

185 Parameters 

186 ---------- 

187 n : int 

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

189 dtype : dtype, optional 

190 The data type of the array to be constructed. 

191 

192 Returns 

193 ------- 

194 H : (n, n) ndarray 

195 The Hadamard matrix. 

196 

197 Notes 

198 ----- 

199 .. versionadded:: 0.8.0 

200 

201 Examples 

202 -------- 

203 >>> from scipy.linalg import hadamard 

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

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

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

207 >>> hadamard(4) 

208 array([[ 1, 1, 1, 1], 

209 [ 1, -1, 1, -1], 

210 [ 1, 1, -1, -1], 

211 [ 1, -1, -1, 1]]) 

212 

213 """ 

214 

215 # This function is a slightly modified version of the 

216 # function contributed by Ivo in ticket #675. 

217 

218 if n < 1: 

219 lg2 = 0 

220 else: 

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

222 if 2 ** lg2 != n: 

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

224 "a power of 2") 

225 

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

227 

228 # Sylvester's construction 

229 for i in range(0, lg2): 

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

231 

232 return H 

233 

234 

235def leslie(f, s): 

236 """ 

237 Create a Leslie matrix. 

238 

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

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

241 matrix. 

242 

243 Parameters 

244 ---------- 

245 f : (N,) array_like 

246 The "fecundity" coefficients. 

247 s : (N-1,) array_like 

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

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

250 

251 Returns 

252 ------- 

253 L : (N, N) ndarray 

254 The array is zero except for the first row, 

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

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

257 

258 Notes 

259 ----- 

260 .. versionadded:: 0.8.0 

261 

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

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

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

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

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

267 per-capita survival rate of each age class. 

268 

269 References 

270 ---------- 

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

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

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

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

275 (Dec. 1948) 

276 

277 Examples 

278 -------- 

279 >>> from scipy.linalg import leslie 

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

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

282 [ 0.2, 0. , 0. , 0. ], 

283 [ 0. , 0.8, 0. , 0. ], 

284 [ 0. , 0. , 0.7, 0. ]]) 

285 

286 """ 

287 f = np.atleast_1d(f) 

288 s = np.atleast_1d(s) 

289 if f.ndim != 1: 

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

291 if s.ndim != 1: 

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

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

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

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

296 if s.size == 0: 

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

298 

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

300 n = f.size 

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

302 a[0] = f 

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

304 return a 

305 

306 

307def kron(a, b): 

308 """ 

309 Kronecker product. 

310 

311 The result is the block matrix:: 

312 

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

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

315 ... 

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

317 

318 Parameters 

319 ---------- 

320 a : (M, N) ndarray 

321 Input array 

322 b : (P, Q) ndarray 

323 Input array 

324 

325 Returns 

326 ------- 

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

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

329 

330 Examples 

331 -------- 

332 >>> from numpy import array 

333 >>> from scipy.linalg import kron 

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

335 array([[1, 1, 1, 2, 2, 2], 

336 [3, 3, 3, 4, 4, 4]]) 

337 

338 """ 

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

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

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

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

343 o = np.outer(a, b) 

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

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

346 

347 

348def block_diag(*arrs): 

349 """ 

350 Create a block diagonal matrix from provided arrays. 

351 

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

353 arrays arranged on the diagonal:: 

354 

355 [[A, 0, 0], 

356 [0, B, 0], 

357 [0, 0, C]] 

358 

359 Parameters 

360 ---------- 

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

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

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

364 

365 Returns 

366 ------- 

367 D : ndarray 

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

369 same dtype as `A`. 

370 

371 Notes 

372 ----- 

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

374 block diagonal matrix. 

375 

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

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

378 

379 Examples 

380 -------- 

381 >>> import numpy as np 

382 >>> from scipy.linalg import block_diag 

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

384 ... [0, 1]] 

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

386 ... [6, 7, 8]] 

387 >>> C = [[7]] 

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

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

390 array([[1, 0, 0, 0, 0, 0], 

391 [0, 1, 0, 0, 0, 0], 

392 [0, 0, 3, 4, 5, 0], 

393 [0, 0, 6, 7, 8, 0], 

394 [0, 0, 0, 0, 0, 7]]) 

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

396 array([[1, 0, 0, 0, 0, 0], 

397 [0, 1, 0, 0, 0, 0], 

398 [0, 0, 0, 0, 0, 0], 

399 [0, 0, 0, 0, 0, 0], 

400 [0, 0, 3, 4, 5, 0], 

401 [0, 0, 6, 7, 8, 0], 

402 [0, 0, 0, 0, 0, 7]]) 

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

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

405 [ 0., 2., 3., 0., 0.], 

406 [ 0., 0., 0., 4., 5.], 

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

408 

409 """ 

410 if arrs == (): 

411 arrs = ([],) 

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

413 

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

415 if bad_args: 

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

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

418 

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

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

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

422 

423 r, c = 0, 0 

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

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

426 r += rr 

427 c += cc 

428 return out 

429 

430 

431def companion(a): 

432 """ 

433 Create a companion matrix. 

434 

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

436 coefficients are given in `a`. 

437 

438 Parameters 

439 ---------- 

440 a : (N,) array_like 

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

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

443 

444 Returns 

445 ------- 

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

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

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

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

450 

451 Raises 

452 ------ 

453 ValueError 

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

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

456 

457 Notes 

458 ----- 

459 .. versionadded:: 0.8.0 

460 

461 References 

462 ---------- 

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

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

465 

466 Examples 

467 -------- 

468 >>> from scipy.linalg import companion 

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

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

471 [ 1., 0., 0.], 

472 [ 0., 1., 0.]]) 

473 

474 """ 

475 a = np.atleast_1d(a) 

476 

477 if a.ndim != 1: 

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

479 "one-dimensional.") 

480 

481 if a.size < 2: 

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

483 

484 if a[0] == 0: 

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

486 

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

488 n = a.size 

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

490 c[0] = first_row 

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

492 return c 

493 

494 

495def helmert(n, full=False): 

496 """ 

497 Create an Helmert matrix of order `n`. 

498 

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

500 and in Aitchison geometry. 

501 

502 Parameters 

503 ---------- 

504 n : int 

505 The size of the array to create. 

506 full : bool, optional 

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

508 Otherwise the submatrix that does not include the first 

509 row will be returned. 

510 Default: False. 

511 

512 Returns 

513 ------- 

514 M : ndarray 

515 The Helmert matrix. 

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

517 

518 Examples 

519 -------- 

520 >>> from scipy.linalg import helmert 

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

522 array([[ 0.4472136 , 0.4472136 , 0.4472136 , 0.4472136 , 0.4472136 ], 

523 [ 0.70710678, -0.70710678, 0. , 0. , 0. ], 

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

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

526 [ 0.2236068 , 0.2236068 , 0.2236068 , 0.2236068 , -0.89442719]]) 

527 

528 """ 

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

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

531 H[0] = 1 

532 d[0] = n 

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

534 if full: 

535 return H_full 

536 else: 

537 return H_full[1:] 

538 

539 

540def hilbert(n): 

541 """ 

542 Create a Hilbert matrix of order `n`. 

543 

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

545 

546 Parameters 

547 ---------- 

548 n : int 

549 The size of the array to create. 

550 

551 Returns 

552 ------- 

553 h : (n, n) ndarray 

554 The Hilbert matrix. 

555 

556 See Also 

557 -------- 

558 invhilbert : Compute the inverse of a Hilbert matrix. 

559 

560 Notes 

561 ----- 

562 .. versionadded:: 0.10.0 

563 

564 Examples 

565 -------- 

566 >>> from scipy.linalg import hilbert 

567 >>> hilbert(3) 

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

569 [ 0.5 , 0.33333333, 0.25 ], 

570 [ 0.33333333, 0.25 , 0.2 ]]) 

571 

572 """ 

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

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

575 return h 

576 

577 

578def invhilbert(n, exact=False): 

579 """ 

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

581 

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

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

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

585 dealing with these large integers. 

586 

587 Parameters 

588 ---------- 

589 n : int 

590 The order of the Hilbert matrix. 

591 exact : bool, optional 

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

593 and the array is an approximation of the inverse. 

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

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

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

597 array with data type np.int64. 

598 

599 Returns 

600 ------- 

601 invh : (n, n) ndarray 

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

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

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

605 array will be long integers. 

606 

607 See Also 

608 -------- 

609 hilbert : Create a Hilbert matrix. 

610 

611 Notes 

612 ----- 

613 .. versionadded:: 0.10.0 

614 

615 Examples 

616 -------- 

617 >>> from scipy.linalg import invhilbert 

618 >>> invhilbert(4) 

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

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

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

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

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

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

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

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

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

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

629 4.2475099528537506e+19 

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

631 42475099528537378560 

632 

633 """ 

634 from scipy.special import comb 

635 if exact: 

636 if n > 14: 

637 dtype = object 

638 else: 

639 dtype = np.int64 

640 else: 

641 dtype = np.float64 

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

643 for i in range(n): 

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

645 s = i + j 

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

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

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

649 comb(s, i, exact=exact) ** 2) 

650 if i != j: 

651 invh[j, i] = invh[i, j] 

652 return invh 

653 

654 

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

656 """ 

657 Returns the n x n Pascal matrix. 

658 

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

660 its elements. 

661 

662 Parameters 

663 ---------- 

664 n : int 

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

666 matrix. 

667 kind : str, optional 

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

669 Default is 'symmetric'. 

670 exact : bool, optional 

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

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

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

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

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

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

677 

678 Returns 

679 ------- 

680 p : (n, n) ndarray 

681 The Pascal matrix. 

682 

683 See Also 

684 -------- 

685 invpascal 

686 

687 Notes 

688 ----- 

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

690 about Pascal matrices. 

691 

692 .. versionadded:: 0.11.0 

693 

694 Examples 

695 -------- 

696 >>> from scipy.linalg import pascal 

697 >>> pascal(4) 

698 array([[ 1, 1, 1, 1], 

699 [ 1, 2, 3, 4], 

700 [ 1, 3, 6, 10], 

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

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

703 array([[1, 0, 0, 0], 

704 [1, 1, 0, 0], 

705 [1, 2, 1, 0], 

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

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

708 25477612258980856902730428600 

709 >>> from scipy.special import comb 

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

711 25477612258980856902730428600 

712 

713 """ 

714 

715 from scipy.special import comb 

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

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

718 

719 if exact: 

720 if n >= 35: 

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

722 L_n.fill(0) 

723 else: 

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

725 for i in range(n): 

726 for j in range(i + 1): 

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

728 else: 

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

730 

731 if kind == 'lower': 

732 p = L_n 

733 elif kind == 'upper': 

734 p = L_n.T 

735 else: 

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

737 

738 return p 

739 

740 

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

742 """ 

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

744 

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

746 its elements. 

747 

748 Parameters 

749 ---------- 

750 n : int 

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

752 matrix. 

753 kind : str, optional 

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

755 Default is 'symmetric'. 

756 exact : bool, optional 

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

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

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

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

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

762 exact coefficients. 

763 

764 Returns 

765 ------- 

766 invp : (n, n) ndarray 

767 The inverse of the Pascal matrix. 

768 

769 See Also 

770 -------- 

771 pascal 

772 

773 Notes 

774 ----- 

775 

776 .. versionadded:: 0.16.0 

777 

778 References 

779 ---------- 

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

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

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

783 

784 Examples 

785 -------- 

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

787 >>> invp = invpascal(5) 

788 >>> invp 

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

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

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

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

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

794 

795 >>> p = pascal(5) 

796 >>> p.dot(invp) 

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

798 [ 0., 1., 0., 0., 0.], 

799 [ 0., 0., 1., 0., 0.], 

800 [ 0., 0., 0., 1., 0.], 

801 [ 0., 0., 0., 0., 1.]]) 

802 

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

804 

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

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

807 [-1., 1., -0., 0., -0.], 

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

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

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

811 

812 """ 

813 from scipy.special import comb 

814 

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

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

817 

818 if kind == 'symmetric': 

819 if exact: 

820 if n > 34: 

821 dt = object 

822 else: 

823 dt = np.int64 

824 else: 

825 dt = np.float64 

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

827 for i in range(n): 

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

829 v = 0 

830 for k in range(n - i): 

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

832 exact=exact) 

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

834 if i != j: 

835 invp[j, i] = invp[i, j] 

836 else: 

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

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

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

840 if invp.dtype == np.uint64: 

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

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

843 invp = invp.view(np.int64) 

844 

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

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

847 

848 return invp 

849 

850 

851def dft(n, scale=None): 

852 """ 

853 Discrete Fourier transform matrix. 

854 

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

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

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

858 

859 Parameters 

860 ---------- 

861 n : int 

862 Size the matrix to create. 

863 scale : str, optional 

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

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

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

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

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

869 

870 Returns 

871 ------- 

872 m : (n, n) ndarray 

873 The DFT matrix. 

874 

875 Notes 

876 ----- 

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

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

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

880 

881 .. versionadded:: 0.14.0 

882 

883 References 

884 ---------- 

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

886 

887 Examples 

888 -------- 

889 >>> import numpy as np 

890 >>> from scipy.linalg import dft 

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

892 >>> m = dft(5) 

893 >>> m 

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

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

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

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

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

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

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

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

902 

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

904 

905 >>> from scipy.fft import fft 

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

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

908 """ 

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

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

911 f"{scale!r} is not valid.") 

912 

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

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

915 if scale == 'sqrtn': 

916 m /= math.sqrt(n) 

917 elif scale == 'n': 

918 m /= n 

919 return m 

920 

921 

922def fiedler(a): 

923 """Returns a symmetric Fiedler matrix 

924 

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

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

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

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

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

930 

931 Parameters 

932 ---------- 

933 a : (n,) array_like 

934 coefficient array 

935 

936 Returns 

937 ------- 

938 F : (n, n) ndarray 

939 

940 See Also 

941 -------- 

942 circulant, toeplitz 

943 

944 Notes 

945 ----- 

946 

947 .. versionadded:: 1.3.0 

948 

949 References 

950 ---------- 

951 .. [1] J. Todd, "Basic Numerical Mathematics: Vol.2 : Numerical Algebra", 

952 1977, Birkhauser, :doi:`10.1007/978-3-0348-7286-7` 

953 

954 Examples 

955 -------- 

956 >>> import numpy as np 

957 >>> from scipy.linalg import det, inv, fiedler 

958 >>> a = [1, 4, 12, 45, 77] 

959 >>> n = len(a) 

960 >>> A = fiedler(a) 

961 >>> A 

962 array([[ 0, 3, 11, 44, 76], 

963 [ 3, 0, 8, 41, 73], 

964 [11, 8, 0, 33, 65], 

965 [44, 41, 33, 0, 32], 

966 [76, 73, 65, 32, 0]]) 

967 

968 The explicit formulas for determinant and inverse seem to hold only for 

969 monotonically increasing/decreasing arrays. Note the tridiagonal structure 

970 and the corners. 

971 

972 >>> Ai = inv(A) 

973 >>> Ai[np.abs(Ai) < 1e-12] = 0. # cleanup the numerical noise for display 

974 >>> Ai 

975 array([[-0.16008772, 0.16666667, 0. , 0. , 0.00657895], 

976 [ 0.16666667, -0.22916667, 0.0625 , 0. , 0. ], 

977 [ 0. , 0.0625 , -0.07765152, 0.01515152, 0. ], 

978 [ 0. , 0. , 0.01515152, -0.03077652, 0.015625 ], 

979 [ 0.00657895, 0. , 0. , 0.015625 , -0.00904605]]) 

980 >>> det(A) 

981 15409151.999999998 

982 >>> (-1)**(n-1) * 2**(n-2) * np.diff(a).prod() * (a[-1] - a[0]) 

983 15409152 

984 

985 """ 

986 a = np.atleast_1d(a) 

987 

988 if a.ndim != 1: 

989 raise ValueError("Input 'a' must be a 1D array.") 

990 

991 if a.size == 0: 

992 return np.array([], dtype=float) 

993 elif a.size == 1: 

994 return np.array([[0.]]) 

995 else: 

996 return np.abs(a[:, None] - a) 

997 

998 

999def fiedler_companion(a): 

1000 """ Returns a Fiedler companion matrix 

1001 

1002 Given a polynomial coefficient array ``a``, this function forms a 

1003 pentadiagonal matrix with a special structure whose eigenvalues coincides 

1004 with the roots of ``a``. 

1005 

1006 Parameters 

1007 ---------- 

1008 a : (N,) array_like 

1009 1-D array of polynomial coefficients in descending order with a nonzero 

1010 leading coefficient. For ``N < 2``, an empty array is returned. 

1011 

1012 Returns 

1013 ------- 

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

1015 Resulting companion matrix 

1016 

1017 See Also 

1018 -------- 

1019 companion 

1020 

1021 Notes 

1022 ----- 

1023 Similar to `companion` the leading coefficient should be nonzero. In the case 

1024 the leading coefficient is not 1, other coefficients are rescaled before 

1025 the array generation. To avoid numerical issues, it is best to provide a 

1026 monic polynomial. 

1027 

1028 .. versionadded:: 1.3.0 

1029 

1030 References 

1031 ---------- 

1032 .. [1] M. Fiedler, " A note on companion matrices", Linear Algebra and its 

1033 Applications, 2003, :doi:`10.1016/S0024-3795(03)00548-2` 

1034 

1035 Examples 

1036 -------- 

1037 >>> import numpy as np 

1038 >>> from scipy.linalg import fiedler_companion, eigvals 

1039 >>> p = np.poly(np.arange(1, 9, 2)) # [1., -16., 86., -176., 105.] 

1040 >>> fc = fiedler_companion(p) 

1041 >>> fc 

1042 array([[ 16., -86., 1., 0.], 

1043 [ 1., 0., 0., 0.], 

1044 [ 0., 176., 0., -105.], 

1045 [ 0., 1., 0., 0.]]) 

1046 >>> eigvals(fc) 

1047 array([7.+0.j, 5.+0.j, 3.+0.j, 1.+0.j]) 

1048 

1049 """ 

1050 a = np.atleast_1d(a) 

1051 

1052 if a.ndim != 1: 

1053 raise ValueError("Input 'a' must be a 1-D array.") 

1054 

1055 if a.size <= 2: 

1056 if a.size == 2: 

1057 return np.array([[-(a/a[0])[-1]]]) 

1058 return np.array([], dtype=a.dtype) 

1059 

1060 if a[0] == 0.: 

1061 raise ValueError('Leading coefficient is zero.') 

1062 

1063 a = a/a[0] 

1064 n = a.size - 1