Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/linalg/_special_matrices.py: 10%
248 statements
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-23 06:43 +0000
1import math
2import warnings
4import numpy as np
5from numpy.lib.stride_tricks import as_strided
7__all__ = ['tri', 'tril', 'triu', 'toeplitz', 'circulant', 'hankel',
8 'hadamard', 'leslie', 'kron', 'block_diag', 'companion',
9 'helmert', 'hilbert', 'invhilbert', 'pascal', 'invpascal', 'dft',
10 'fiedler', 'fiedler_companion', 'convolution_matrix']
13# -----------------------------------------------------------------------------
14# matrix construction functions
15# -----------------------------------------------------------------------------
17#
18# *Note*: tri{,u,l} is implemented in NumPy, but an important bug was fixed in
19# 2.0.0.dev-1af2f3, the following tri{,u,l} definitions are here for backwards
20# compatibility.
22def tri(N, M=None, k=0, dtype=None):
23 """
24 .. deprecated:: 1.11.0
25 `tri` is deprecated in favour of `numpy.tri` and will be removed in
26 SciPy 1.13.0.
28 Construct (N, M) matrix filled with ones at and below the kth diagonal.
30 The matrix has A[i,j] == 1 for j <= i + k
32 Parameters
33 ----------
34 N : int
35 The size of the first dimension of the matrix.
36 M : int or None, optional
37 The size of the second dimension of the matrix. If `M` is None,
38 `M = N` is assumed.
39 k : int, optional
40 Number of subdiagonal below which matrix is filled with ones.
41 `k` = 0 is the main diagonal, `k` < 0 subdiagonal and `k` > 0
42 superdiagonal.
43 dtype : dtype, optional
44 Data type of the matrix.
46 Returns
47 -------
48 tri : (N, M) ndarray
49 Tri matrix.
51 Examples
52 --------
53 >>> from scipy.linalg import tri
54 >>> tri(3, 5, 2, dtype=int)
55 array([[1, 1, 1, 0, 0],
56 [1, 1, 1, 1, 0],
57 [1, 1, 1, 1, 1]])
58 >>> tri(3, 5, -1, dtype=int)
59 array([[0, 0, 0, 0, 0],
60 [1, 0, 0, 0, 0],
61 [1, 1, 0, 0, 0]])
63 """
64 warnings.warn("'tri'/'tril/'triu' are deprecated as of SciPy 1.11.0 and "
65 "will be removed in v1.13.0. Please use "
66 "numpy.(tri/tril/triu) instead.",
67 DeprecationWarning, stacklevel=2)
69 if M is None:
70 M = N
71 if isinstance(M, str):
72 # pearu: any objections to remove this feature?
73 # As tri(N,'d') is equivalent to tri(N,dtype='d')
74 dtype = M
75 M = N
76 m = np.greater_equal.outer(np.arange(k, N+k), np.arange(M))
77 if dtype is None:
78 return m
79 else:
80 return m.astype(dtype)
83def tril(m, k=0):
84 """
85 .. deprecated:: 1.11.0
86 `tril` is deprecated in favour of `numpy.tril` and will be removed in
87 SciPy 1.13.0.
89 Make a copy of a matrix with elements above the kth diagonal zeroed.
91 Parameters
92 ----------
93 m : array_like
94 Matrix whose elements to return
95 k : int, optional
96 Diagonal above which to zero elements.
97 `k` == 0 is the main diagonal, `k` < 0 subdiagonal and
98 `k` > 0 superdiagonal.
100 Returns
101 -------
102 tril : ndarray
103 Return is the same shape and type as `m`.
105 Examples
106 --------
107 >>> from scipy.linalg import tril
108 >>> tril([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
109 array([[ 0, 0, 0],
110 [ 4, 0, 0],
111 [ 7, 8, 0],
112 [10, 11, 12]])
114 """
115 m = np.asarray(m)
116 out = tri(m.shape[0], m.shape[1], k=k, dtype=m.dtype.char) * m
117 return out
120def triu(m, k=0):
121 """
122 .. deprecated:: 1.11.0
123 `tril` is deprecated in favour of `numpy.triu` and will be removed in
124 SciPy 1.13.0.
126 Make a copy of a matrix with elements below the kth diagonal zeroed.
128 Parameters
129 ----------
130 m : array_like
131 Matrix whose elements to return
132 k : int, optional
133 Diagonal below which to zero elements.
134 `k` == 0 is the main diagonal, `k` < 0 subdiagonal and
135 `k` > 0 superdiagonal.
137 Returns
138 -------
139 triu : ndarray
140 Return matrix with zeroed elements below the kth diagonal and has
141 same shape and type as `m`.
143 Examples
144 --------
145 >>> from scipy.linalg import triu
146 >>> triu([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
147 array([[ 1, 2, 3],
148 [ 4, 5, 6],
149 [ 0, 8, 9],
150 [ 0, 0, 12]])
152 """
153 m = np.asarray(m)
154 out = (1 - tri(m.shape[0], m.shape[1], k - 1, m.dtype.char)) * m
155 return out
158def toeplitz(c, r=None):
159 """
160 Construct a Toeplitz matrix.
162 The Toeplitz matrix has constant diagonals, with c as its first column
163 and r as its first row. If r is not given, ``r == conjugate(c)`` is
164 assumed.
166 Parameters
167 ----------
168 c : array_like
169 First column of the matrix. Whatever the actual shape of `c`, it
170 will be converted to a 1-D array.
171 r : array_like, optional
172 First row of the matrix. If None, ``r = conjugate(c)`` is assumed;
173 in this case, if c[0] is real, the result is a Hermitian matrix.
174 r[0] is ignored; the first row of the returned matrix is
175 ``[c[0], r[1:]]``. Whatever the actual shape of `r`, it will be
176 converted to a 1-D array.
178 Returns
179 -------
180 A : (len(c), len(r)) ndarray
181 The Toeplitz matrix. Dtype is the same as ``(c[0] + r[0]).dtype``.
183 See Also
184 --------
185 circulant : circulant matrix
186 hankel : Hankel matrix
187 solve_toeplitz : Solve a Toeplitz system.
189 Notes
190 -----
191 The behavior when `c` or `r` is a scalar, or when `c` is complex and
192 `r` is None, was changed in version 0.8.0. The behavior in previous
193 versions was undocumented and is no longer supported.
195 Examples
196 --------
197 >>> from scipy.linalg import toeplitz
198 >>> toeplitz([1,2,3], [1,4,5,6])
199 array([[1, 4, 5, 6],
200 [2, 1, 4, 5],
201 [3, 2, 1, 4]])
202 >>> toeplitz([1.0, 2+3j, 4-1j])
203 array([[ 1.+0.j, 2.-3.j, 4.+1.j],
204 [ 2.+3.j, 1.+0.j, 2.-3.j],
205 [ 4.-1.j, 2.+3.j, 1.+0.j]])
207 """
208 c = np.asarray(c).ravel()
209 if r is None:
210 r = c.conjugate()
211 else:
212 r = np.asarray(r).ravel()
213 # Form a 1-D array containing a reversed c followed by r[1:] that could be
214 # strided to give us toeplitz matrix.
215 vals = np.concatenate((c[::-1], r[1:]))
216 out_shp = len(c), len(r)
217 n = vals.strides[0]
218 return as_strided(vals[len(c)-1:], shape=out_shp, strides=(-n, n)).copy()
221def circulant(c):
222 """
223 Construct a circulant matrix.
225 Parameters
226 ----------
227 c : (N,) array_like
228 1-D array, the first column of the matrix.
230 Returns
231 -------
232 A : (N, N) ndarray
233 A circulant matrix whose first column is `c`.
235 See Also
236 --------
237 toeplitz : Toeplitz matrix
238 hankel : Hankel matrix
239 solve_circulant : Solve a circulant system.
241 Notes
242 -----
243 .. versionadded:: 0.8.0
245 Examples
246 --------
247 >>> from scipy.linalg import circulant
248 >>> circulant([1, 2, 3])
249 array([[1, 3, 2],
250 [2, 1, 3],
251 [3, 2, 1]])
253 """
254 c = np.asarray(c).ravel()
255 # Form an extended array that could be strided to give circulant version
256 c_ext = np.concatenate((c[::-1], c[:0:-1]))
257 L = len(c)
258 n = c_ext.strides[0]
259 return as_strided(c_ext[L-1:], shape=(L, L), strides=(-n, n)).copy()
262def hankel(c, r=None):
263 """
264 Construct a Hankel matrix.
266 The Hankel matrix has constant anti-diagonals, with `c` as its
267 first column and `r` as its last row. If `r` is not given, then
268 `r = zeros_like(c)` is assumed.
270 Parameters
271 ----------
272 c : array_like
273 First column of the matrix. Whatever the actual shape of `c`, it
274 will be converted to a 1-D array.
275 r : array_like, optional
276 Last row of the matrix. If None, ``r = zeros_like(c)`` is assumed.
277 r[0] is ignored; the last row of the returned matrix is
278 ``[c[-1], r[1:]]``. Whatever the actual shape of `r`, it will be
279 converted to a 1-D array.
281 Returns
282 -------
283 A : (len(c), len(r)) ndarray
284 The Hankel matrix. Dtype is the same as ``(c[0] + r[0]).dtype``.
286 See Also
287 --------
288 toeplitz : Toeplitz matrix
289 circulant : circulant matrix
291 Examples
292 --------
293 >>> from scipy.linalg import hankel
294 >>> hankel([1, 17, 99])
295 array([[ 1, 17, 99],
296 [17, 99, 0],
297 [99, 0, 0]])
298 >>> hankel([1,2,3,4], [4,7,7,8,9])
299 array([[1, 2, 3, 4, 7],
300 [2, 3, 4, 7, 7],
301 [3, 4, 7, 7, 8],
302 [4, 7, 7, 8, 9]])
304 """
305 c = np.asarray(c).ravel()
306 if r is None:
307 r = np.zeros_like(c)
308 else:
309 r = np.asarray(r).ravel()
310 # Form a 1-D array of values to be used in the matrix, containing `c`
311 # followed by r[1:].
312 vals = np.concatenate((c, r[1:]))
313 # Stride on concatenated array to get hankel matrix
314 out_shp = len(c), len(r)
315 n = vals.strides[0]
316 return as_strided(vals, shape=out_shp, strides=(n, n)).copy()
319def hadamard(n, dtype=int):
320 """
321 Construct an Hadamard matrix.
323 Constructs an n-by-n Hadamard matrix, using Sylvester's
324 construction. `n` must be a power of 2.
326 Parameters
327 ----------
328 n : int
329 The order of the matrix. `n` must be a power of 2.
330 dtype : dtype, optional
331 The data type of the array to be constructed.
333 Returns
334 -------
335 H : (n, n) ndarray
336 The Hadamard matrix.
338 Notes
339 -----
340 .. versionadded:: 0.8.0
342 Examples
343 --------
344 >>> from scipy.linalg import hadamard
345 >>> hadamard(2, dtype=complex)
346 array([[ 1.+0.j, 1.+0.j],
347 [ 1.+0.j, -1.-0.j]])
348 >>> hadamard(4)
349 array([[ 1, 1, 1, 1],
350 [ 1, -1, 1, -1],
351 [ 1, 1, -1, -1],
352 [ 1, -1, -1, 1]])
354 """
356 # This function is a slightly modified version of the
357 # function contributed by Ivo in ticket #675.
359 if n < 1:
360 lg2 = 0
361 else:
362 lg2 = int(math.log(n, 2))
363 if 2 ** lg2 != n:
364 raise ValueError("n must be an positive integer, and n must be "
365 "a power of 2")
367 H = np.array([[1]], dtype=dtype)
369 # Sylvester's construction
370 for i in range(0, lg2):
371 H = np.vstack((np.hstack((H, H)), np.hstack((H, -H))))
373 return H
376def leslie(f, s):
377 """
378 Create a Leslie matrix.
380 Given the length n array of fecundity coefficients `f` and the length
381 n-1 array of survival coefficients `s`, return the associated Leslie
382 matrix.
384 Parameters
385 ----------
386 f : (N,) array_like
387 The "fecundity" coefficients.
388 s : (N-1,) array_like
389 The "survival" coefficients, has to be 1-D. The length of `s`
390 must be one less than the length of `f`, and it must be at least 1.
392 Returns
393 -------
394 L : (N, N) ndarray
395 The array is zero except for the first row,
396 which is `f`, and the first sub-diagonal, which is `s`.
397 The data-type of the array will be the data-type of ``f[0]+s[0]``.
399 Notes
400 -----
401 .. versionadded:: 0.8.0
403 The Leslie matrix is used to model discrete-time, age-structured
404 population growth [1]_ [2]_. In a population with `n` age classes, two sets
405 of parameters define a Leslie matrix: the `n` "fecundity coefficients",
406 which give the number of offspring per-capita produced by each age
407 class, and the `n` - 1 "survival coefficients", which give the
408 per-capita survival rate of each age class.
410 References
411 ----------
412 .. [1] P. H. Leslie, On the use of matrices in certain population
413 mathematics, Biometrika, Vol. 33, No. 3, 183--212 (Nov. 1945)
414 .. [2] P. H. Leslie, Some further notes on the use of matrices in
415 population mathematics, Biometrika, Vol. 35, No. 3/4, 213--245
416 (Dec. 1948)
418 Examples
419 --------
420 >>> from scipy.linalg import leslie
421 >>> leslie([0.1, 2.0, 1.0, 0.1], [0.2, 0.8, 0.7])
422 array([[ 0.1, 2. , 1. , 0.1],
423 [ 0.2, 0. , 0. , 0. ],
424 [ 0. , 0.8, 0. , 0. ],
425 [ 0. , 0. , 0.7, 0. ]])
427 """
428 f = np.atleast_1d(f)
429 s = np.atleast_1d(s)
430 if f.ndim != 1:
431 raise ValueError("Incorrect shape for f. f must be 1D")
432 if s.ndim != 1:
433 raise ValueError("Incorrect shape for s. s must be 1D")
434 if f.size != s.size + 1:
435 raise ValueError("Incorrect lengths for f and s. The length"
436 " of s must be one less than the length of f.")
437 if s.size == 0:
438 raise ValueError("The length of s must be at least 1.")
440 tmp = f[0] + s[0]
441 n = f.size
442 a = np.zeros((n, n), dtype=tmp.dtype)
443 a[0] = f
444 a[list(range(1, n)), list(range(0, n - 1))] = s
445 return a
448def kron(a, b):
449 """
450 Kronecker product.
452 The result is the block matrix::
454 a[0,0]*b a[0,1]*b ... a[0,-1]*b
455 a[1,0]*b a[1,1]*b ... a[1,-1]*b
456 ...
457 a[-1,0]*b a[-1,1]*b ... a[-1,-1]*b
459 Parameters
460 ----------
461 a : (M, N) ndarray
462 Input array
463 b : (P, Q) ndarray
464 Input array
466 Returns
467 -------
468 A : (M*P, N*Q) ndarray
469 Kronecker product of `a` and `b`.
471 Examples
472 --------
473 >>> from numpy import array
474 >>> from scipy.linalg import kron
475 >>> kron(array([[1,2],[3,4]]), array([[1,1,1]]))
476 array([[1, 1, 1, 2, 2, 2],
477 [3, 3, 3, 4, 4, 4]])
479 """
480 if not a.flags['CONTIGUOUS']:
481 a = np.reshape(a, a.shape)
482 if not b.flags['CONTIGUOUS']:
483 b = np.reshape(b, b.shape)
484 o = np.outer(a, b)
485 o = o.reshape(a.shape + b.shape)
486 return np.concatenate(np.concatenate(o, axis=1), axis=1)
489def block_diag(*arrs):
490 """
491 Create a block diagonal matrix from provided arrays.
493 Given the inputs `A`, `B` and `C`, the output will have these
494 arrays arranged on the diagonal::
496 [[A, 0, 0],
497 [0, B, 0],
498 [0, 0, C]]
500 Parameters
501 ----------
502 A, B, C, ... : array_like, up to 2-D
503 Input arrays. A 1-D array or array_like sequence of length `n` is
504 treated as a 2-D array with shape ``(1,n)``.
506 Returns
507 -------
508 D : ndarray
509 Array with `A`, `B`, `C`, ... on the diagonal. `D` has the
510 same dtype as `A`.
512 Notes
513 -----
514 If all the input arrays are square, the output is known as a
515 block diagonal matrix.
517 Empty sequences (i.e., array-likes of zero size) will not be ignored.
518 Noteworthy, both [] and [[]] are treated as matrices with shape ``(1,0)``.
520 Examples
521 --------
522 >>> import numpy as np
523 >>> from scipy.linalg import block_diag
524 >>> A = [[1, 0],
525 ... [0, 1]]
526 >>> B = [[3, 4, 5],
527 ... [6, 7, 8]]
528 >>> C = [[7]]
529 >>> P = np.zeros((2, 0), dtype='int32')
530 >>> block_diag(A, B, C)
531 array([[1, 0, 0, 0, 0, 0],
532 [0, 1, 0, 0, 0, 0],
533 [0, 0, 3, 4, 5, 0],
534 [0, 0, 6, 7, 8, 0],
535 [0, 0, 0, 0, 0, 7]])
536 >>> block_diag(A, P, B, C)
537 array([[1, 0, 0, 0, 0, 0],
538 [0, 1, 0, 0, 0, 0],
539 [0, 0, 0, 0, 0, 0],
540 [0, 0, 0, 0, 0, 0],
541 [0, 0, 3, 4, 5, 0],
542 [0, 0, 6, 7, 8, 0],
543 [0, 0, 0, 0, 0, 7]])
544 >>> block_diag(1.0, [2, 3], [[4, 5], [6, 7]])
545 array([[ 1., 0., 0., 0., 0.],
546 [ 0., 2., 3., 0., 0.],
547 [ 0., 0., 0., 4., 5.],
548 [ 0., 0., 0., 6., 7.]])
550 """
551 if arrs == ():
552 arrs = ([],)
553 arrs = [np.atleast_2d(a) for a in arrs]
555 bad_args = [k for k in range(len(arrs)) if arrs[k].ndim > 2]
556 if bad_args:
557 raise ValueError("arguments in the following positions have dimension "
558 "greater than 2: %s" % bad_args)
560 shapes = np.array([a.shape for a in arrs])
561 out_dtype = np.result_type(*[arr.dtype for arr in arrs])
562 out = np.zeros(np.sum(shapes, axis=0), dtype=out_dtype)
564 r, c = 0, 0
565 for i, (rr, cc) in enumerate(shapes):
566 out[r:r + rr, c:c + cc] = arrs[i]
567 r += rr
568 c += cc
569 return out
572def companion(a):
573 """
574 Create a companion matrix.
576 Create the companion matrix [1]_ associated with the polynomial whose
577 coefficients are given in `a`.
579 Parameters
580 ----------
581 a : (N,) array_like
582 1-D array of polynomial coefficients. The length of `a` must be
583 at least two, and ``a[0]`` must not be zero.
585 Returns
586 -------
587 c : (N-1, N-1) ndarray
588 The first row of `c` is ``-a[1:]/a[0]``, and the first
589 sub-diagonal is all ones. The data-type of the array is the same
590 as the data-type of ``1.0*a[0]``.
592 Raises
593 ------
594 ValueError
595 If any of the following are true: a) ``a.ndim != 1``;
596 b) ``a.size < 2``; c) ``a[0] == 0``.
598 Notes
599 -----
600 .. versionadded:: 0.8.0
602 References
603 ----------
604 .. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK:
605 Cambridge University Press, 1999, pp. 146-7.
607 Examples
608 --------
609 >>> from scipy.linalg import companion
610 >>> companion([1, -10, 31, -30])
611 array([[ 10., -31., 30.],
612 [ 1., 0., 0.],
613 [ 0., 1., 0.]])
615 """
616 a = np.atleast_1d(a)
618 if a.ndim != 1:
619 raise ValueError("Incorrect shape for `a`. `a` must be "
620 "one-dimensional.")
622 if a.size < 2:
623 raise ValueError("The length of `a` must be at least 2.")
625 if a[0] == 0:
626 raise ValueError("The first coefficient in `a` must not be zero.")
628 first_row = -a[1:] / (1.0 * a[0])
629 n = a.size
630 c = np.zeros((n - 1, n - 1), dtype=first_row.dtype)
631 c[0] = first_row
632 c[list(range(1, n - 1)), list(range(0, n - 2))] = 1
633 return c
636def helmert(n, full=False):
637 """
638 Create an Helmert matrix of order `n`.
640 This has applications in statistics, compositional or simplicial analysis,
641 and in Aitchison geometry.
643 Parameters
644 ----------
645 n : int
646 The size of the array to create.
647 full : bool, optional
648 If True the (n, n) ndarray will be returned.
649 Otherwise the submatrix that does not include the first
650 row will be returned.
651 Default: False.
653 Returns
654 -------
655 M : ndarray
656 The Helmert matrix.
657 The shape is (n, n) or (n-1, n) depending on the `full` argument.
659 Examples
660 --------
661 >>> from scipy.linalg import helmert
662 >>> helmert(5, full=True)
663 array([[ 0.4472136 , 0.4472136 , 0.4472136 , 0.4472136 , 0.4472136 ],
664 [ 0.70710678, -0.70710678, 0. , 0. , 0. ],
665 [ 0.40824829, 0.40824829, -0.81649658, 0. , 0. ],
666 [ 0.28867513, 0.28867513, 0.28867513, -0.8660254 , 0. ],
667 [ 0.2236068 , 0.2236068 , 0.2236068 , 0.2236068 , -0.89442719]])
669 """
670 H = np.tril(np.ones((n, n)), -1) - np.diag(np.arange(n))
671 d = np.arange(n) * np.arange(1, n+1)
672 H[0] = 1
673 d[0] = n
674 H_full = H / np.sqrt(d)[:, np.newaxis]
675 if full:
676 return H_full
677 else:
678 return H_full[1:]
681def hilbert(n):
682 """
683 Create a Hilbert matrix of order `n`.
685 Returns the `n` by `n` array with entries `h[i,j] = 1 / (i + j + 1)`.
687 Parameters
688 ----------
689 n : int
690 The size of the array to create.
692 Returns
693 -------
694 h : (n, n) ndarray
695 The Hilbert matrix.
697 See Also
698 --------
699 invhilbert : Compute the inverse of a Hilbert matrix.
701 Notes
702 -----
703 .. versionadded:: 0.10.0
705 Examples
706 --------
707 >>> from scipy.linalg import hilbert
708 >>> hilbert(3)
709 array([[ 1. , 0.5 , 0.33333333],
710 [ 0.5 , 0.33333333, 0.25 ],
711 [ 0.33333333, 0.25 , 0.2 ]])
713 """
714 values = 1.0 / (1.0 + np.arange(2 * n - 1))
715 h = hankel(values[:n], r=values[n - 1:])
716 return h
719def invhilbert(n, exact=False):
720 """
721 Compute the inverse of the Hilbert matrix of order `n`.
723 The entries in the inverse of a Hilbert matrix are integers. When `n`
724 is greater than 14, some entries in the inverse exceed the upper limit
725 of 64 bit integers. The `exact` argument provides two options for
726 dealing with these large integers.
728 Parameters
729 ----------
730 n : int
731 The order of the Hilbert matrix.
732 exact : bool, optional
733 If False, the data type of the array that is returned is np.float64,
734 and the array is an approximation of the inverse.
735 If True, the array is the exact integer inverse array. To represent
736 the exact inverse when n > 14, the returned array is an object array
737 of long integers. For n <= 14, the exact inverse is returned as an
738 array with data type np.int64.
740 Returns
741 -------
742 invh : (n, n) ndarray
743 The data type of the array is np.float64 if `exact` is False.
744 If `exact` is True, the data type is either np.int64 (for n <= 14)
745 or object (for n > 14). In the latter case, the objects in the
746 array will be long integers.
748 See Also
749 --------
750 hilbert : Create a Hilbert matrix.
752 Notes
753 -----
754 .. versionadded:: 0.10.0
756 Examples
757 --------
758 >>> from scipy.linalg import invhilbert
759 >>> invhilbert(4)
760 array([[ 16., -120., 240., -140.],
761 [ -120., 1200., -2700., 1680.],
762 [ 240., -2700., 6480., -4200.],
763 [ -140., 1680., -4200., 2800.]])
764 >>> invhilbert(4, exact=True)
765 array([[ 16, -120, 240, -140],
766 [ -120, 1200, -2700, 1680],
767 [ 240, -2700, 6480, -4200],
768 [ -140, 1680, -4200, 2800]], dtype=int64)
769 >>> invhilbert(16)[7,7]
770 4.2475099528537506e+19
771 >>> invhilbert(16, exact=True)[7,7]
772 42475099528537378560
774 """
775 from scipy.special import comb
776 if exact:
777 if n > 14:
778 dtype = object
779 else:
780 dtype = np.int64
781 else:
782 dtype = np.float64
783 invh = np.empty((n, n), dtype=dtype)
784 for i in range(n):
785 for j in range(0, i + 1):
786 s = i + j
787 invh[i, j] = ((-1) ** s * (s + 1) *
788 comb(n + i, n - j - 1, exact=exact) *
789 comb(n + j, n - i - 1, exact=exact) *
790 comb(s, i, exact=exact) ** 2)
791 if i != j:
792 invh[j, i] = invh[i, j]
793 return invh
796def pascal(n, kind='symmetric', exact=True):
797 """
798 Returns the n x n Pascal matrix.
800 The Pascal matrix is a matrix containing the binomial coefficients as
801 its elements.
803 Parameters
804 ----------
805 n : int
806 The size of the matrix to create; that is, the result is an n x n
807 matrix.
808 kind : str, optional
809 Must be one of 'symmetric', 'lower', or 'upper'.
810 Default is 'symmetric'.
811 exact : bool, optional
812 If `exact` is True, the result is either an array of type
813 numpy.uint64 (if n < 35) or an object array of Python long integers.
814 If `exact` is False, the coefficients in the matrix are computed using
815 `scipy.special.comb` with `exact=False`. The result will be a floating
816 point array, and the values in the array will not be the exact
817 coefficients, but this version is much faster than `exact=True`.
819 Returns
820 -------
821 p : (n, n) ndarray
822 The Pascal matrix.
824 See Also
825 --------
826 invpascal
828 Notes
829 -----
830 See https://en.wikipedia.org/wiki/Pascal_matrix for more information
831 about Pascal matrices.
833 .. versionadded:: 0.11.0
835 Examples
836 --------
837 >>> from scipy.linalg import pascal
838 >>> pascal(4)
839 array([[ 1, 1, 1, 1],
840 [ 1, 2, 3, 4],
841 [ 1, 3, 6, 10],
842 [ 1, 4, 10, 20]], dtype=uint64)
843 >>> pascal(4, kind='lower')
844 array([[1, 0, 0, 0],
845 [1, 1, 0, 0],
846 [1, 2, 1, 0],
847 [1, 3, 3, 1]], dtype=uint64)
848 >>> pascal(50)[-1, -1]
849 25477612258980856902730428600
850 >>> from scipy.special import comb
851 >>> comb(98, 49, exact=True)
852 25477612258980856902730428600
854 """
856 from scipy.special import comb
857 if kind not in ['symmetric', 'lower', 'upper']:
858 raise ValueError("kind must be 'symmetric', 'lower', or 'upper'")
860 if exact:
861 if n >= 35:
862 L_n = np.empty((n, n), dtype=object)
863 L_n.fill(0)
864 else:
865 L_n = np.zeros((n, n), dtype=np.uint64)
866 for i in range(n):
867 for j in range(i + 1):
868 L_n[i, j] = comb(i, j, exact=True)
869 else:
870 L_n = comb(*np.ogrid[:n, :n])
872 if kind == 'lower':
873 p = L_n
874 elif kind == 'upper':
875 p = L_n.T
876 else:
877 p = np.dot(L_n, L_n.T)
879 return p
882def invpascal(n, kind='symmetric', exact=True):
883 """
884 Returns the inverse of the n x n Pascal matrix.
886 The Pascal matrix is a matrix containing the binomial coefficients as
887 its elements.
889 Parameters
890 ----------
891 n : int
892 The size of the matrix to create; that is, the result is an n x n
893 matrix.
894 kind : str, optional
895 Must be one of 'symmetric', 'lower', or 'upper'.
896 Default is 'symmetric'.
897 exact : bool, optional
898 If `exact` is True, the result is either an array of type
899 ``numpy.int64`` (if `n` <= 35) or an object array of Python integers.
900 If `exact` is False, the coefficients in the matrix are computed using
901 `scipy.special.comb` with `exact=False`. The result will be a floating
902 point array, and for large `n`, the values in the array will not be the
903 exact coefficients.
905 Returns
906 -------
907 invp : (n, n) ndarray
908 The inverse of the Pascal matrix.
910 See Also
911 --------
912 pascal
914 Notes
915 -----
917 .. versionadded:: 0.16.0
919 References
920 ----------
921 .. [1] "Pascal matrix", https://en.wikipedia.org/wiki/Pascal_matrix
922 .. [2] Cohen, A. M., "The inverse of a Pascal matrix", Mathematical
923 Gazette, 59(408), pp. 111-112, 1975.
925 Examples
926 --------
927 >>> from scipy.linalg import invpascal, pascal
928 >>> invp = invpascal(5)
929 >>> invp
930 array([[ 5, -10, 10, -5, 1],
931 [-10, 30, -35, 19, -4],
932 [ 10, -35, 46, -27, 6],
933 [ -5, 19, -27, 17, -4],
934 [ 1, -4, 6, -4, 1]])
936 >>> p = pascal(5)
937 >>> p.dot(invp)
938 array([[ 1., 0., 0., 0., 0.],
939 [ 0., 1., 0., 0., 0.],
940 [ 0., 0., 1., 0., 0.],
941 [ 0., 0., 0., 1., 0.],
942 [ 0., 0., 0., 0., 1.]])
944 An example of the use of `kind` and `exact`:
946 >>> invpascal(5, kind='lower', exact=False)
947 array([[ 1., -0., 0., -0., 0.],
948 [-1., 1., -0., 0., -0.],
949 [ 1., -2., 1., -0., 0.],
950 [-1., 3., -3., 1., -0.],
951 [ 1., -4., 6., -4., 1.]])
953 """
954 from scipy.special import comb
956 if kind not in ['symmetric', 'lower', 'upper']:
957 raise ValueError("'kind' must be 'symmetric', 'lower' or 'upper'.")
959 if kind == 'symmetric':
960 if exact:
961 if n > 34:
962 dt = object
963 else:
964 dt = np.int64
965 else:
966 dt = np.float64
967 invp = np.empty((n, n), dtype=dt)
968 for i in range(n):
969 for j in range(0, i + 1):
970 v = 0
971 for k in range(n - i):
972 v += comb(i + k, k, exact=exact) * comb(i + k, i + k - j,
973 exact=exact)
974 invp[i, j] = (-1)**(i - j) * v
975 if i != j:
976 invp[j, i] = invp[i, j]
977 else:
978 # For the 'lower' and 'upper' cases, we computer the inverse by
979 # changing the sign of every other diagonal of the pascal matrix.
980 invp = pascal(n, kind=kind, exact=exact)
981 if invp.dtype == np.uint64:
982 # This cast from np.uint64 to int64 OK, because if `kind` is not
983 # "symmetric", the values in invp are all much less than 2**63.
984 invp = invp.view(np.int64)
986 # The toeplitz matrix has alternating bands of 1 and -1.
987 invp *= toeplitz((-1)**np.arange(n)).astype(invp.dtype)
989 return invp
992def dft(n, scale=None):
993 """
994 Discrete Fourier transform matrix.
996 Create the matrix that computes the discrete Fourier transform of a
997 sequence [1]_. The nth primitive root of unity used to generate the
998 matrix is exp(-2*pi*i/n), where i = sqrt(-1).
1000 Parameters
1001 ----------
1002 n : int
1003 Size the matrix to create.
1004 scale : str, optional
1005 Must be None, 'sqrtn', or 'n'.
1006 If `scale` is 'sqrtn', the matrix is divided by `sqrt(n)`.
1007 If `scale` is 'n', the matrix is divided by `n`.
1008 If `scale` is None (the default), the matrix is not normalized, and the
1009 return value is simply the Vandermonde matrix of the roots of unity.
1011 Returns
1012 -------
1013 m : (n, n) ndarray
1014 The DFT matrix.
1016 Notes
1017 -----
1018 When `scale` is None, multiplying a vector by the matrix returned by
1019 `dft` is mathematically equivalent to (but much less efficient than)
1020 the calculation performed by `scipy.fft.fft`.
1022 .. versionadded:: 0.14.0
1024 References
1025 ----------
1026 .. [1] "DFT matrix", https://en.wikipedia.org/wiki/DFT_matrix
1028 Examples
1029 --------
1030 >>> import numpy as np
1031 >>> from scipy.linalg import dft
1032 >>> np.set_printoptions(precision=2, suppress=True) # for compact output
1033 >>> m = dft(5)
1034 >>> m
1035 array([[ 1. +0.j , 1. +0.j , 1. +0.j , 1. +0.j , 1. +0.j ],
1036 [ 1. +0.j , 0.31-0.95j, -0.81-0.59j, -0.81+0.59j, 0.31+0.95j],
1037 [ 1. +0.j , -0.81-0.59j, 0.31+0.95j, 0.31-0.95j, -0.81+0.59j],
1038 [ 1. +0.j , -0.81+0.59j, 0.31-0.95j, 0.31+0.95j, -0.81-0.59j],
1039 [ 1. +0.j , 0.31+0.95j, -0.81+0.59j, -0.81-0.59j, 0.31-0.95j]])
1040 >>> x = np.array([1, 2, 3, 0, 3])
1041 >>> m @ x # Compute the DFT of x
1042 array([ 9. +0.j , 0.12-0.81j, -2.12+3.44j, -2.12-3.44j, 0.12+0.81j])
1044 Verify that ``m @ x`` is the same as ``fft(x)``.
1046 >>> from scipy.fft import fft
1047 >>> fft(x) # Same result as m @ x
1048 array([ 9. +0.j , 0.12-0.81j, -2.12+3.44j, -2.12-3.44j, 0.12+0.81j])
1049 """
1050 if scale not in [None, 'sqrtn', 'n']:
1051 raise ValueError("scale must be None, 'sqrtn', or 'n'; "
1052 "{!r} is not valid.".format(scale))
1054 omegas = np.exp(-2j * np.pi * np.arange(n) / n).reshape(-1, 1)
1055 m = omegas ** np.arange(n)
1056 if scale == 'sqrtn':
1057 m /= math.sqrt(n)
1058 elif scale == 'n':
1059 m /= n
1060 return m
1063def fiedler(a):
1064 """Returns a symmetric Fiedler matrix
1066 Given an sequence of numbers `a`, Fiedler matrices have the structure
1067 ``F[i, j] = np.abs(a[i] - a[j])``, and hence zero diagonals and nonnegative
1068 entries. A Fiedler matrix has a dominant positive eigenvalue and other
1069 eigenvalues are negative. Although not valid generally, for certain inputs,
1070 the inverse and the determinant can be derived explicitly as given in [1]_.
1072 Parameters
1073 ----------
1074 a : (n,) array_like
1075 coefficient array
1077 Returns
1078 -------
1079 F : (n, n) ndarray
1081 See Also
1082 --------
1083 circulant, toeplitz
1085 Notes
1086 -----
1088 .. versionadded:: 1.3.0
1090 References
1091 ----------
1092 .. [1] J. Todd, "Basic Numerical Mathematics: Vol.2 : Numerical Algebra",
1093 1977, Birkhauser, :doi:`10.1007/978-3-0348-7286-7`
1095 Examples
1096 --------
1097 >>> import numpy as np
1098 >>> from scipy.linalg import det, inv, fiedler
1099 >>> a = [1, 4, 12, 45, 77]
1100 >>> n = len(a)
1101 >>> A = fiedler(a)
1102 >>> A
1103 array([[ 0, 3, 11, 44, 76],
1104 [ 3, 0, 8, 41, 73],
1105 [11, 8, 0, 33, 65],
1106 [44, 41, 33, 0, 32],
1107 [76, 73, 65, 32, 0]])
1109 The explicit formulas for determinant and inverse seem to hold only for
1110 monotonically increasing/decreasing arrays. Note the tridiagonal structure
1111 and the corners.
1113 >>> Ai = inv(A)
1114 >>> Ai[np.abs(Ai) < 1e-12] = 0. # cleanup the numerical noise for display
1115 >>> Ai
1116 array([[-0.16008772, 0.16666667, 0. , 0. , 0.00657895],
1117 [ 0.16666667, -0.22916667, 0.0625 , 0. , 0. ],
1118 [ 0. , 0.0625 , -0.07765152, 0.01515152, 0. ],
1119 [ 0. , 0. , 0.01515152, -0.03077652, 0.015625 ],
1120 [ 0.00657895, 0. , 0. , 0.015625 , -0.00904605]])
1121 >>> det(A)
1122 15409151.999999998
1123 >>> (-1)**(n-1) * 2**(n-2) * np.diff(a).prod() * (a[-1] - a[0])
1124 15409152
1126 """
1127 a = np.atleast_1d(a)
1129 if a.ndim != 1:
1130 raise ValueError("Input 'a' must be a 1D array.")
1132 if a.size == 0:
1133 return np.array([], dtype=float)
1134 elif a.size == 1:
1135 return np.array([[0.]])
1136 else:
1137 return np.abs(a[:, None] - a)
1140def fiedler_companion(a):
1141 """ Returns a Fiedler companion matrix
1143 Given a polynomial coefficient array ``a``, this function forms a
1144 pentadiagonal matrix with a special structure whose eigenvalues coincides
1145 with the roots of ``a``.
1147 Parameters
1148 ----------
1149 a : (N,) array_like
1150 1-D array of polynomial coefficients in descending order with a nonzero
1151 leading coefficient. For ``N < 2``, an empty array is returned.
1153 Returns
1154 -------
1155 c : (N-1, N-1) ndarray
1156 Resulting companion matrix
1158 See Also
1159 --------
1160 companion
1162 Notes
1163 -----
1164 Similar to `companion` the leading coefficient should be nonzero. In the case
1165 the leading coefficient is not 1, other coefficients are rescaled before
1166 the array generation. To avoid numerical issues, it is best to provide a
1167 monic polynomial.
1169 .. versionadded:: 1.3.0
1171 References
1172 ----------
1173 .. [1] M. Fiedler, " A note on companion matrices", Linear Algebra and its
1174 Applications, 2003, :doi:`10.1016/S0024-3795(03)00548-2`
1176 Examples
1177 --------
1178 >>> import numpy as np
1179 >>> from scipy.linalg import fiedler_companion, eigvals
1180 >>> p = np.poly(np.arange(1, 9, 2)) # [1., -16., 86., -176., 105.]
1181 >>> fc = fiedler_companion(p)
1182 >>> fc
1183 array([[ 16., -86., 1., 0.],
1184 [ 1., 0., 0., 0.],
1185 [ 0., 176., 0., -105.],
1186 [ 0., 1., 0., 0.]])
1187 >>> eigvals(fc)
1188 array([7.+0.j, 5.+0.j, 3.+0.j, 1.+0.j])
1190 """
1191 a = np.atleast_1d(a)
1193 if a.ndim != 1:
1194 raise ValueError("Input 'a' must be a 1-D array.")
1196 if a.size <= 2:
1197 if a.size == 2:
1198 return np.array([[-(a/a[0])[-1]]])
1199 return np.array([], dtype=a.dtype)
1201 if a[0] == 0.:
1202 raise ValueError('Leading coefficient is zero.')
1204 a = a/a[0]
1205 n = a.size - 1
1206 c = np.zeros((n, n), dtype=a.dtype)
1207 # subdiagonals
1208 c[range(3, n, 2), range(1, n-2, 2)] = 1.
1209 c[range(2, n, 2), range(1, n-1, 2)] = -a[3::2]
1210 # superdiagonals
1211 c[range(0, n-2, 2), range(2, n, 2)] = 1.
1212 c[range(0, n-1, 2), range(1, n, 2)] = -a[2::2]
1213 c[[0, 1], 0] = [-a[1], 1]
1215 return c
1218def convolution_matrix(a, n, mode='full'):
1219 """
1220 Construct a convolution matrix.
1222 Constructs the Toeplitz matrix representing one-dimensional
1223 convolution [1]_. See the notes below for details.
1225 Parameters
1226 ----------
1227 a : (m,) array_like
1228 The 1-D array to convolve.
1229 n : int
1230 The number of columns in the resulting matrix. It gives the length
1231 of the input to be convolved with `a`. This is analogous to the
1232 length of `v` in ``numpy.convolve(a, v)``.
1233 mode : str
1234 This is analogous to `mode` in ``numpy.convolve(v, a, mode)``.
1235 It must be one of ('full', 'valid', 'same').
1236 See below for how `mode` determines the shape of the result.
1238 Returns
1239 -------
1240 A : (k, n) ndarray
1241 The convolution matrix whose row count `k` depends on `mode`::
1243 ======= =========================
1244 mode k
1245 ======= =========================
1246 'full' m + n -1
1247 'same' max(m, n)
1248 'valid' max(m, n) - min(m, n) + 1
1249 ======= =========================
1251 See Also
1252 --------
1253 toeplitz : Toeplitz matrix
1255 Notes
1256 -----
1257 The code::
1259 A = convolution_matrix(a, n, mode)
1261 creates a Toeplitz matrix `A` such that ``A @ v`` is equivalent to
1262 using ``convolve(a, v, mode)``. The returned array always has `n`
1263 columns. The number of rows depends on the specified `mode`, as
1264 explained above.
1266 In the default 'full' mode, the entries of `A` are given by::
1268 A[i, j] == (a[i-j] if (0 <= (i-j) < m) else 0)
1270 where ``m = len(a)``. Suppose, for example, the input array is
1271 ``[x, y, z]``. The convolution matrix has the form::
1273 [x, 0, 0, ..., 0, 0]
1274 [y, x, 0, ..., 0, 0]
1275 [z, y, x, ..., 0, 0]
1276 ...
1277 [0, 0, 0, ..., x, 0]
1278 [0, 0, 0, ..., y, x]
1279 [0, 0, 0, ..., z, y]
1280 [0, 0, 0, ..., 0, z]
1282 In 'valid' mode, the entries of `A` are given by::
1284 A[i, j] == (a[i-j+m-1] if (0 <= (i-j+m-1) < m) else 0)
1286 This corresponds to a matrix whose rows are the subset of those from
1287 the 'full' case where all the coefficients in `a` are contained in the
1288 row. For input ``[x, y, z]``, this array looks like::
1290 [z, y, x, 0, 0, ..., 0, 0, 0]
1291 [0, z, y, x, 0, ..., 0, 0, 0]
1292 [0, 0, z, y, x, ..., 0, 0, 0]
1293 ...
1294 [0, 0, 0, 0, 0, ..., x, 0, 0]
1295 [0, 0, 0, 0, 0, ..., y, x, 0]
1296 [0, 0, 0, 0, 0, ..., z, y, x]
1298 In the 'same' mode, the entries of `A` are given by::
1300 d = (m - 1) // 2
1301 A[i, j] == (a[i-j+d] if (0 <= (i-j+d) < m) else 0)
1303 The typical application of the 'same' mode is when one has a signal of
1304 length `n` (with `n` greater than ``len(a)``), and the desired output
1305 is a filtered signal that is still of length `n`.
1307 For input ``[x, y, z]``, this array looks like::
1309 [y, x, 0, 0, ..., 0, 0, 0]
1310 [z, y, x, 0, ..., 0, 0, 0]
1311 [0, z, y, x, ..., 0, 0, 0]
1312 [0, 0, z, y, ..., 0, 0, 0]
1313 ...
1314 [0, 0, 0, 0, ..., y, x, 0]
1315 [0, 0, 0, 0, ..., z, y, x]
1316 [0, 0, 0, 0, ..., 0, z, y]
1318 .. versionadded:: 1.5.0
1320 References
1321 ----------
1322 .. [1] "Convolution", https://en.wikipedia.org/wiki/Convolution
1324 Examples
1325 --------
1326 >>> import numpy as np
1327 >>> from scipy.linalg import convolution_matrix
1328 >>> A = convolution_matrix([-1, 4, -2], 5, mode='same')
1329 >>> A
1330 array([[ 4, -1, 0, 0, 0],
1331 [-2, 4, -1, 0, 0],
1332 [ 0, -2, 4, -1, 0],
1333 [ 0, 0, -2, 4, -1],
1334 [ 0, 0, 0, -2, 4]])
1336 Compare multiplication by `A` with the use of `numpy.convolve`.
1338 >>> x = np.array([1, 2, 0, -3, 0.5])
1339 >>> A @ x
1340 array([ 2. , 6. , -1. , -12.5, 8. ])
1342 Verify that ``A @ x`` produced the same result as applying the
1343 convolution function.
1345 >>> np.convolve([-1, 4, -2], x, mode='same')
1346 array([ 2. , 6. , -1. , -12.5, 8. ])
1348 For comparison to the case ``mode='same'`` shown above, here are the
1349 matrices produced by ``mode='full'`` and ``mode='valid'`` for the
1350 same coefficients and size.
1352 >>> convolution_matrix([-1, 4, -2], 5, mode='full')
1353 array([[-1, 0, 0, 0, 0],
1354 [ 4, -1, 0, 0, 0],
1355 [-2, 4, -1, 0, 0],
1356 [ 0, -2, 4, -1, 0],
1357 [ 0, 0, -2, 4, -1],
1358 [ 0, 0, 0, -2, 4],
1359 [ 0, 0, 0, 0, -2]])
1361 >>> convolution_matrix([-1, 4, -2], 5, mode='valid')
1362 array([[-2, 4, -1, 0, 0],
1363 [ 0, -2, 4, -1, 0],
1364 [ 0, 0, -2, 4, -1]])
1365 """
1366 if n <= 0:
1367 raise ValueError('n must be a positive integer.')
1369 a = np.asarray(a)
1370 if a.ndim != 1:
1371 raise ValueError('convolution_matrix expects a one-dimensional '
1372 'array as input')
1373 if a.size == 0:
1374 raise ValueError('len(a) must be at least 1.')
1376 if mode not in ('full', 'valid', 'same'):
1377 raise ValueError(
1378 "'mode' argument must be one of ('full', 'valid', 'same')")
1380 # create zero padded versions of the array
1381 az = np.pad(a, (0, n-1), 'constant')
1382 raz = np.pad(a[::-1], (0, n-1), 'constant')
1384 if mode == 'same':
1385 trim = min(n, len(a)) - 1
1386 tb = trim//2
1387 te = trim - tb
1388 col0 = az[tb:len(az)-te]
1389 row0 = raz[-n-tb:len(raz)-tb]
1390 elif mode == 'valid':
1391 tb = min(n, len(a)) - 1
1392 te = tb
1393 col0 = az[tb:len(az)-te]
1394 row0 = raz[-n-tb:len(raz)-tb]
1395 else: # 'full'
1396 col0 = az
1397 row0 = raz[-n:]
1398 return toeplitz(col0, row0)