Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/linalg/special_matrices.py : 10%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import math
2import numpy as np
3from numpy.lib.stride_tricks import as_strided
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']
11# -----------------------------------------------------------------------------
12# matrix construction functions
13# -----------------------------------------------------------------------------
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.
20def tri(N, M=None, k=0, dtype=None):
21 """
22 Construct (N, M) matrix filled with ones at and below the kth diagonal.
24 The matrix has A[i,j] == 1 for i <= j + k
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.
40 Returns
41 -------
42 tri : (N, M) ndarray
43 Tri matrix.
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]])
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)
72def tril(m, k=0):
73 """
74 Make a copy of a matrix with elements above the kth diagonal zeroed.
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.
85 Returns
86 -------
87 tril : ndarray
88 Return is the same shape and type as `m`.
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]])
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
105def triu(m, k=0):
106 """
107 Make a copy of a matrix with elements below the kth diagonal zeroed.
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.
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`.
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]])
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
139def toeplitz(c, r=None):
140 """
141 Construct a Toeplitz matrix.
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.
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.
159 Returns
160 -------
161 A : (len(c), len(r)) ndarray
162 The Toeplitz matrix. Dtype is the same as ``(c[0] + r[0]).dtype``.
164 See Also
165 --------
166 circulant : circulant matrix
167 hankel : Hankel matrix
168 solve_toeplitz : Solve a Toeplitz system.
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.
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]])
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()
202def circulant(c):
203 """
204 Construct a circulant matrix.
206 Parameters
207 ----------
208 c : (N,) array_like
209 1-D array, the first column of the matrix.
211 Returns
212 -------
213 A : (N, N) ndarray
214 A circulant matrix whose first column is `c`.
216 See Also
217 --------
218 toeplitz : Toeplitz matrix
219 hankel : Hankel matrix
220 solve_circulant : Solve a circulant system.
222 Notes
223 -----
224 .. versionadded:: 0.8.0
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]])
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()
243def hankel(c, r=None):
244 """
245 Construct a Hankel matrix.
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.
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.
262 Returns
263 -------
264 A : (len(c), len(r)) ndarray
265 The Hankel matrix. Dtype is the same as ``(c[0] + r[0]).dtype``.
267 See Also
268 --------
269 toeplitz : Toeplitz matrix
270 circulant : circulant matrix
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]])
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()
300def hadamard(n, dtype=int):
301 """
302 Construct an Hadamard matrix.
304 Constructs an n-by-n Hadamard matrix, using Sylvester's
305 construction. `n` must be a power of 2.
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.
314 Returns
315 -------
316 H : (n, n) ndarray
317 The Hadamard matrix.
319 Notes
320 -----
321 .. versionadded:: 0.8.0
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]])
335 """
337 # This function is a slightly modified version of the
338 # function contributed by Ivo in ticket #675.
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")
348 H = np.array([[1]], dtype=dtype)
350 # Sylvester's construction
351 for i in range(0, lg2):
352 H = np.vstack((np.hstack((H, H)), np.hstack((H, -H))))
354 return H
357def leslie(f, s):
358 """
359 Create a Leslie matrix.
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.
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.
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]``.
380 Notes
381 -----
382 .. versionadded:: 0.8.0
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.
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)
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. ]])
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.")
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
429def kron(a, b):
430 """
431 Kronecker product.
433 The result is the block matrix::
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
440 Parameters
441 ----------
442 a : (M, N) ndarray
443 Input array
444 b : (P, Q) ndarray
445 Input array
447 Returns
448 -------
449 A : (M*P, N*Q) ndarray
450 Kronecker product of `a` and `b`.
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]])
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)
470def block_diag(*arrs):
471 """
472 Create a block diagonal matrix from provided arrays.
474 Given the inputs `A`, `B` and `C`, the output will have these
475 arrays arranged on the diagonal::
477 [[A, 0, 0],
478 [0, B, 0],
479 [0, 0, C]]
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)``.
487 Returns
488 -------
489 D : ndarray
490 Array with `A`, `B`, `C`, ... on the diagonal. `D` has the
491 same dtype as `A`.
493 Notes
494 -----
495 If all the input arrays are square, the output is known as a
496 block diagonal matrix.
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)``.
501 Examples
502 --------
503 >>> from scipy.linalg import block_diag
504 >>> A = [[1, 0],
505 ... [0, 1]]
506 >>> B = [[3, 4, 5],
507 ... [6, 7, 8]]
508 >>> C = [[7]]
509 >>> P = np.zeros((2, 0), dtype='int32')
510 >>> block_diag(A, B, C)
511 array([[1, 0, 0, 0, 0, 0],
512 [0, 1, 0, 0, 0, 0],
513 [0, 0, 3, 4, 5, 0],
514 [0, 0, 6, 7, 8, 0],
515 [0, 0, 0, 0, 0, 7]])
516 >>> block_diag(A, P, B, C)
517 array([[1, 0, 0, 0, 0, 0],
518 [0, 1, 0, 0, 0, 0],
519 [0, 0, 0, 0, 0, 0],
520 [0, 0, 0, 0, 0, 0],
521 [0, 0, 3, 4, 5, 0],
522 [0, 0, 6, 7, 8, 0],
523 [0, 0, 0, 0, 0, 7]])
524 >>> block_diag(1.0, [2, 3], [[4, 5], [6, 7]])
525 array([[ 1., 0., 0., 0., 0.],
526 [ 0., 2., 3., 0., 0.],
527 [ 0., 0., 0., 4., 5.],
528 [ 0., 0., 0., 6., 7.]])
530 """
531 if arrs == ():
532 arrs = ([],)
533 arrs = [np.atleast_2d(a) for a in arrs]
535 bad_args = [k for k in range(len(arrs)) if arrs[k].ndim > 2]
536 if bad_args:
537 raise ValueError("arguments in the following positions have dimension "
538 "greater than 2: %s" % bad_args)
540 shapes = np.array([a.shape for a in arrs])
541 out_dtype = np.find_common_type([arr.dtype for arr in arrs], [])
542 out = np.zeros(np.sum(shapes, axis=0), dtype=out_dtype)
544 r, c = 0, 0
545 for i, (rr, cc) in enumerate(shapes):
546 out[r:r + rr, c:c + cc] = arrs[i]
547 r += rr
548 c += cc
549 return out
552def companion(a):
553 """
554 Create a companion matrix.
556 Create the companion matrix [1]_ associated with the polynomial whose
557 coefficients are given in `a`.
559 Parameters
560 ----------
561 a : (N,) array_like
562 1-D array of polynomial coefficients. The length of `a` must be
563 at least two, and ``a[0]`` must not be zero.
565 Returns
566 -------
567 c : (N-1, N-1) ndarray
568 The first row of `c` is ``-a[1:]/a[0]``, and the first
569 sub-diagonal is all ones. The data-type of the array is the same
570 as the data-type of ``1.0*a[0]``.
572 Raises
573 ------
574 ValueError
575 If any of the following are true: a) ``a.ndim != 1``;
576 b) ``a.size < 2``; c) ``a[0] == 0``.
578 Notes
579 -----
580 .. versionadded:: 0.8.0
582 References
583 ----------
584 .. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK:
585 Cambridge University Press, 1999, pp. 146-7.
587 Examples
588 --------
589 >>> from scipy.linalg import companion
590 >>> companion([1, -10, 31, -30])
591 array([[ 10., -31., 30.],
592 [ 1., 0., 0.],
593 [ 0., 1., 0.]])
595 """
596 a = np.atleast_1d(a)
598 if a.ndim != 1:
599 raise ValueError("Incorrect shape for `a`. `a` must be "
600 "one-dimensional.")
602 if a.size < 2:
603 raise ValueError("The length of `a` must be at least 2.")
605 if a[0] == 0:
606 raise ValueError("The first coefficient in `a` must not be zero.")
608 first_row = -a[1:] / (1.0 * a[0])
609 n = a.size
610 c = np.zeros((n - 1, n - 1), dtype=first_row.dtype)
611 c[0] = first_row
612 c[list(range(1, n - 1)), list(range(0, n - 2))] = 1
613 return c
616def helmert(n, full=False):
617 """
618 Create an Helmert matrix of order `n`.
620 This has applications in statistics, compositional or simplicial analysis,
621 and in Aitchison geometry.
623 Parameters
624 ----------
625 n : int
626 The size of the array to create.
627 full : bool, optional
628 If True the (n, n) ndarray will be returned.
629 Otherwise the submatrix that does not include the first
630 row will be returned.
631 Default: False.
633 Returns
634 -------
635 M : ndarray
636 The Helmert matrix.
637 The shape is (n, n) or (n-1, n) depending on the `full` argument.
639 Examples
640 --------
641 >>> from scipy.linalg import helmert
642 >>> helmert(5, full=True)
643 array([[ 0.4472136 , 0.4472136 , 0.4472136 , 0.4472136 , 0.4472136 ],
644 [ 0.70710678, -0.70710678, 0. , 0. , 0. ],
645 [ 0.40824829, 0.40824829, -0.81649658, 0. , 0. ],
646 [ 0.28867513, 0.28867513, 0.28867513, -0.8660254 , 0. ],
647 [ 0.2236068 , 0.2236068 , 0.2236068 , 0.2236068 , -0.89442719]])
649 """
650 H = np.tril(np.ones((n, n)), -1) - np.diag(np.arange(n))
651 d = np.arange(n) * np.arange(1, n+1)
652 H[0] = 1
653 d[0] = n
654 H_full = H / np.sqrt(d)[:, np.newaxis]
655 if full:
656 return H_full
657 else:
658 return H_full[1:]
661def hilbert(n):
662 """
663 Create a Hilbert matrix of order `n`.
665 Returns the `n` by `n` array with entries `h[i,j] = 1 / (i + j + 1)`.
667 Parameters
668 ----------
669 n : int
670 The size of the array to create.
672 Returns
673 -------
674 h : (n, n) ndarray
675 The Hilbert matrix.
677 See Also
678 --------
679 invhilbert : Compute the inverse of a Hilbert matrix.
681 Notes
682 -----
683 .. versionadded:: 0.10.0
685 Examples
686 --------
687 >>> from scipy.linalg import hilbert
688 >>> hilbert(3)
689 array([[ 1. , 0.5 , 0.33333333],
690 [ 0.5 , 0.33333333, 0.25 ],
691 [ 0.33333333, 0.25 , 0.2 ]])
693 """
694 values = 1.0 / (1.0 + np.arange(2 * n - 1))
695 h = hankel(values[:n], r=values[n - 1:])
696 return h
699def invhilbert(n, exact=False):
700 """
701 Compute the inverse of the Hilbert matrix of order `n`.
703 The entries in the inverse of a Hilbert matrix are integers. When `n`
704 is greater than 14, some entries in the inverse exceed the upper limit
705 of 64 bit integers. The `exact` argument provides two options for
706 dealing with these large integers.
708 Parameters
709 ----------
710 n : int
711 The order of the Hilbert matrix.
712 exact : bool, optional
713 If False, the data type of the array that is returned is np.float64,
714 and the array is an approximation of the inverse.
715 If True, the array is the exact integer inverse array. To represent
716 the exact inverse when n > 14, the returned array is an object array
717 of long integers. For n <= 14, the exact inverse is returned as an
718 array with data type np.int64.
720 Returns
721 -------
722 invh : (n, n) ndarray
723 The data type of the array is np.float64 if `exact` is False.
724 If `exact` is True, the data type is either np.int64 (for n <= 14)
725 or object (for n > 14). In the latter case, the objects in the
726 array will be long integers.
728 See Also
729 --------
730 hilbert : Create a Hilbert matrix.
732 Notes
733 -----
734 .. versionadded:: 0.10.0
736 Examples
737 --------
738 >>> from scipy.linalg import invhilbert
739 >>> invhilbert(4)
740 array([[ 16., -120., 240., -140.],
741 [ -120., 1200., -2700., 1680.],
742 [ 240., -2700., 6480., -4200.],
743 [ -140., 1680., -4200., 2800.]])
744 >>> invhilbert(4, exact=True)
745 array([[ 16, -120, 240, -140],
746 [ -120, 1200, -2700, 1680],
747 [ 240, -2700, 6480, -4200],
748 [ -140, 1680, -4200, 2800]], dtype=int64)
749 >>> invhilbert(16)[7,7]
750 4.2475099528537506e+19
751 >>> invhilbert(16, exact=True)[7,7]
752 42475099528537378560
754 """
755 from scipy.special import comb
756 if exact:
757 if n > 14:
758 dtype = object
759 else:
760 dtype = np.int64
761 else:
762 dtype = np.float64
763 invh = np.empty((n, n), dtype=dtype)
764 for i in range(n):
765 for j in range(0, i + 1):
766 s = i + j
767 invh[i, j] = ((-1) ** s * (s + 1) *
768 comb(n + i, n - j - 1, exact) *
769 comb(n + j, n - i - 1, exact) *
770 comb(s, i, exact) ** 2)
771 if i != j:
772 invh[j, i] = invh[i, j]
773 return invh
776def pascal(n, kind='symmetric', exact=True):
777 """
778 Returns the n x n Pascal matrix.
780 The Pascal matrix is a matrix containing the binomial coefficients as
781 its elements.
783 Parameters
784 ----------
785 n : int
786 The size of the matrix to create; that is, the result is an n x n
787 matrix.
788 kind : str, optional
789 Must be one of 'symmetric', 'lower', or 'upper'.
790 Default is 'symmetric'.
791 exact : bool, optional
792 If `exact` is True, the result is either an array of type
793 numpy.uint64 (if n < 35) or an object array of Python long integers.
794 If `exact` is False, the coefficients in the matrix are computed using
795 `scipy.special.comb` with `exact=False`. The result will be a floating
796 point array, and the values in the array will not be the exact
797 coefficients, but this version is much faster than `exact=True`.
799 Returns
800 -------
801 p : (n, n) ndarray
802 The Pascal matrix.
804 See Also
805 --------
806 invpascal
808 Notes
809 -----
810 See https://en.wikipedia.org/wiki/Pascal_matrix for more information
811 about Pascal matrices.
813 .. versionadded:: 0.11.0
815 Examples
816 --------
817 >>> from scipy.linalg import pascal
818 >>> pascal(4)
819 array([[ 1, 1, 1, 1],
820 [ 1, 2, 3, 4],
821 [ 1, 3, 6, 10],
822 [ 1, 4, 10, 20]], dtype=uint64)
823 >>> pascal(4, kind='lower')
824 array([[1, 0, 0, 0],
825 [1, 1, 0, 0],
826 [1, 2, 1, 0],
827 [1, 3, 3, 1]], dtype=uint64)
828 >>> pascal(50)[-1, -1]
829 25477612258980856902730428600
830 >>> from scipy.special import comb
831 >>> comb(98, 49, exact=True)
832 25477612258980856902730428600
834 """
836 from scipy.special import comb
837 if kind not in ['symmetric', 'lower', 'upper']:
838 raise ValueError("kind must be 'symmetric', 'lower', or 'upper'")
840 if exact:
841 if n >= 35:
842 L_n = np.empty((n, n), dtype=object)
843 L_n.fill(0)
844 else:
845 L_n = np.zeros((n, n), dtype=np.uint64)
846 for i in range(n):
847 for j in range(i + 1):
848 L_n[i, j] = comb(i, j, exact=True)
849 else:
850 L_n = comb(*np.ogrid[:n, :n])
852 if kind == 'lower':
853 p = L_n
854 elif kind == 'upper':
855 p = L_n.T
856 else:
857 p = np.dot(L_n, L_n.T)
859 return p
862def invpascal(n, kind='symmetric', exact=True):
863 """
864 Returns the inverse of the n x n Pascal matrix.
866 The Pascal matrix is a matrix containing the binomial coefficients as
867 its elements.
869 Parameters
870 ----------
871 n : int
872 The size of the matrix to create; that is, the result is an n x n
873 matrix.
874 kind : str, optional
875 Must be one of 'symmetric', 'lower', or 'upper'.
876 Default is 'symmetric'.
877 exact : bool, optional
878 If `exact` is True, the result is either an array of type
879 ``numpy.int64`` (if `n` <= 35) or an object array of Python integers.
880 If `exact` is False, the coefficients in the matrix are computed using
881 `scipy.special.comb` with `exact=False`. The result will be a floating
882 point array, and for large `n`, the values in the array will not be the
883 exact coefficients.
885 Returns
886 -------
887 invp : (n, n) ndarray
888 The inverse of the Pascal matrix.
890 See Also
891 --------
892 pascal
894 Notes
895 -----
897 .. versionadded:: 0.16.0
899 References
900 ----------
901 .. [1] "Pascal matrix", https://en.wikipedia.org/wiki/Pascal_matrix
902 .. [2] Cohen, A. M., "The inverse of a Pascal matrix", Mathematical
903 Gazette, 59(408), pp. 111-112, 1975.
905 Examples
906 --------
907 >>> from scipy.linalg import invpascal, pascal
908 >>> invp = invpascal(5)
909 >>> invp
910 array([[ 5, -10, 10, -5, 1],
911 [-10, 30, -35, 19, -4],
912 [ 10, -35, 46, -27, 6],
913 [ -5, 19, -27, 17, -4],
914 [ 1, -4, 6, -4, 1]])
916 >>> p = pascal(5)
917 >>> p.dot(invp)
918 array([[ 1., 0., 0., 0., 0.],
919 [ 0., 1., 0., 0., 0.],
920 [ 0., 0., 1., 0., 0.],
921 [ 0., 0., 0., 1., 0.],
922 [ 0., 0., 0., 0., 1.]])
924 An example of the use of `kind` and `exact`:
926 >>> invpascal(5, kind='lower', exact=False)
927 array([[ 1., -0., 0., -0., 0.],
928 [-1., 1., -0., 0., -0.],
929 [ 1., -2., 1., -0., 0.],
930 [-1., 3., -3., 1., -0.],
931 [ 1., -4., 6., -4., 1.]])
933 """
934 from scipy.special import comb
936 if kind not in ['symmetric', 'lower', 'upper']:
937 raise ValueError("'kind' must be 'symmetric', 'lower' or 'upper'.")
939 if kind == 'symmetric':
940 if exact:
941 if n > 34:
942 dt = object
943 else:
944 dt = np.int64
945 else:
946 dt = np.float64
947 invp = np.empty((n, n), dtype=dt)
948 for i in range(n):
949 for j in range(0, i + 1):
950 v = 0
951 for k in range(n - i):
952 v += comb(i + k, k, exact=exact) * comb(i + k, i + k - j,
953 exact=exact)
954 invp[i, j] = (-1)**(i - j) * v
955 if i != j:
956 invp[j, i] = invp[i, j]
957 else:
958 # For the 'lower' and 'upper' cases, we computer the inverse by
959 # changing the sign of every other diagonal of the pascal matrix.
960 invp = pascal(n, kind=kind, exact=exact)
961 if invp.dtype == np.uint64:
962 # This cast from np.uint64 to int64 OK, because if `kind` is not
963 # "symmetric", the values in invp are all much less than 2**63.
964 invp = invp.view(np.int64)
966 # The toeplitz matrix has alternating bands of 1 and -1.
967 invp *= toeplitz((-1)**np.arange(n)).astype(invp.dtype)
969 return invp
972def dft(n, scale=None):
973 """
974 Discrete Fourier transform matrix.
976 Create the matrix that computes the discrete Fourier transform of a
977 sequence [1]_. The nth primitive root of unity used to generate the
978 matrix is exp(-2*pi*i/n), where i = sqrt(-1).
980 Parameters
981 ----------
982 n : int
983 Size the matrix to create.
984 scale : str, optional
985 Must be None, 'sqrtn', or 'n'.
986 If `scale` is 'sqrtn', the matrix is divided by `sqrt(n)`.
987 If `scale` is 'n', the matrix is divided by `n`.
988 If `scale` is None (the default), the matrix is not normalized, and the
989 return value is simply the Vandermonde matrix of the roots of unity.
991 Returns
992 -------
993 m : (n, n) ndarray
994 The DFT matrix.
996 Notes
997 -----
998 When `scale` is None, multiplying a vector by the matrix returned by
999 `dft` is mathematically equivalent to (but much less efficient than)
1000 the calculation performed by `scipy.fft.fft`.
1002 .. versionadded:: 0.14.0
1004 References
1005 ----------
1006 .. [1] "DFT matrix", https://en.wikipedia.org/wiki/DFT_matrix
1008 Examples
1009 --------
1010 >>> from scipy.linalg import dft
1011 >>> np.set_printoptions(precision=2, suppress=True) # for compact output
1012 >>> m = dft(5)
1013 >>> m
1014 array([[ 1. +0.j , 1. +0.j , 1. +0.j , 1. +0.j , 1. +0.j ],
1015 [ 1. +0.j , 0.31-0.95j, -0.81-0.59j, -0.81+0.59j, 0.31+0.95j],
1016 [ 1. +0.j , -0.81-0.59j, 0.31+0.95j, 0.31-0.95j, -0.81+0.59j],
1017 [ 1. +0.j , -0.81+0.59j, 0.31-0.95j, 0.31+0.95j, -0.81-0.59j],
1018 [ 1. +0.j , 0.31+0.95j, -0.81+0.59j, -0.81-0.59j, 0.31-0.95j]])
1019 >>> x = np.array([1, 2, 3, 0, 3])
1020 >>> m @ x # Compute the DFT of x
1021 array([ 9. +0.j , 0.12-0.81j, -2.12+3.44j, -2.12-3.44j, 0.12+0.81j])
1023 Verify that ``m @ x`` is the same as ``fft(x)``.
1025 >>> from scipy.fft import fft
1026 >>> fft(x) # Same result as m @ x
1027 array([ 9. +0.j , 0.12-0.81j, -2.12+3.44j, -2.12-3.44j, 0.12+0.81j])
1028 """
1029 if scale not in [None, 'sqrtn', 'n']:
1030 raise ValueError("scale must be None, 'sqrtn', or 'n'; "
1031 "%r is not valid." % (scale,))
1033 omegas = np.exp(-2j * np.pi * np.arange(n) / n).reshape(-1, 1)
1034 m = omegas ** np.arange(n)
1035 if scale == 'sqrtn':
1036 m /= math.sqrt(n)
1037 elif scale == 'n':
1038 m /= n
1039 return m
1042def fiedler(a):
1043 """Returns a symmetric Fiedler matrix
1045 Given an sequence of numbers `a`, Fiedler matrices have the structure
1046 ``F[i, j] = np.abs(a[i] - a[j])``, and hence zero diagonals and nonnegative
1047 entries. A Fiedler matrix has a dominant positive eigenvalue and other
1048 eigenvalues are negative. Although not valid generally, for certain inputs,
1049 the inverse and the determinant can be derived explicitly as given in [1]_.
1051 Parameters
1052 ----------
1053 a : (n,) array_like
1054 coefficient array
1056 Returns
1057 -------
1058 F : (n, n) ndarray
1060 See Also
1061 --------
1062 circulant, toeplitz
1064 Notes
1065 -----
1067 .. versionadded:: 1.3.0
1069 References
1070 ----------
1071 .. [1] J. Todd, "Basic Numerical Mathematics: Vol.2 : Numerical Algebra",
1072 1977, Birkhauser, :doi:`10.1007/978-3-0348-7286-7`
1074 Examples
1075 --------
1076 >>> from scipy.linalg import det, inv, fiedler
1077 >>> a = [1, 4, 12, 45, 77]
1078 >>> n = len(a)
1079 >>> A = fiedler(a)
1080 >>> A
1081 array([[ 0, 3, 11, 44, 76],
1082 [ 3, 0, 8, 41, 73],
1083 [11, 8, 0, 33, 65],
1084 [44, 41, 33, 0, 32],
1085 [76, 73, 65, 32, 0]])
1087 The explicit formulas for determinant and inverse seem to hold only for
1088 monotonically increasing/decreasing arrays. Note the tridiagonal structure
1089 and the corners.
1091 >>> Ai = inv(A)
1092 >>> Ai[np.abs(Ai) < 1e-12] = 0. # cleanup the numerical noise for display
1093 >>> Ai
1094 array([[-0.16008772, 0.16666667, 0. , 0. , 0.00657895],
1095 [ 0.16666667, -0.22916667, 0.0625 , 0. , 0. ],
1096 [ 0. , 0.0625 , -0.07765152, 0.01515152, 0. ],
1097 [ 0. , 0. , 0.01515152, -0.03077652, 0.015625 ],
1098 [ 0.00657895, 0. , 0. , 0.015625 , -0.00904605]])
1099 >>> det(A)
1100 15409151.999999998
1101 >>> (-1)**(n-1) * 2**(n-2) * np.diff(a).prod() * (a[-1] - a[0])
1102 15409152
1104 """
1105 a = np.atleast_1d(a)
1107 if a.ndim != 1:
1108 raise ValueError("Input 'a' must be a 1D array.")
1110 if a.size == 0:
1111 return np.array([], dtype=float)
1112 elif a.size == 1:
1113 return np.array([[0.]])
1114 else:
1115 return np.abs(a[:, None] - a)
1118def fiedler_companion(a):
1119 """ Returns a Fiedler companion matrix
1121 Given a polynomial coefficient array ``a``, this function forms a
1122 pentadiagonal matrix with a special structure whose eigenvalues coincides
1123 with the roots of ``a``.
1125 Parameters
1126 ----------
1127 a : (N,) array_like
1128 1-D array of polynomial coefficients in descending order with a nonzero
1129 leading coefficient. For ``N < 2``, an empty array is returned.
1131 Returns
1132 -------
1133 c : (N-1, N-1) ndarray
1134 Resulting companion matrix
1136 Notes
1137 -----
1138 Similar to `companion` the leading coefficient should be nonzero. In the case
1139 the leading coefficient is not 1, other coefficients are rescaled before
1140 the array generation. To avoid numerical issues, it is best to provide a
1141 monic polynomial.
1143 .. versionadded:: 1.3.0
1145 See Also
1146 --------
1147 companion
1149 References
1150 ----------
1151 .. [1] M. Fiedler, " A note on companion matrices", Linear Algebra and its
1152 Applications, 2003, :doi:`10.1016/S0024-3795(03)00548-2`
1154 Examples
1155 --------
1156 >>> from scipy.linalg import fiedler_companion, eigvals
1157 >>> p = np.poly(np.arange(1, 9, 2)) # [1., -16., 86., -176., 105.]
1158 >>> fc = fiedler_companion(p)
1159 >>> fc
1160 array([[ 16., -86., 1., 0.],
1161 [ 1., 0., 0., 0.],
1162 [ 0., 176., 0., -105.],
1163 [ 0., 1., 0., 0.]])
1164 >>> eigvals(fc)
1165 array([7.+0.j, 5.+0.j, 3.+0.j, 1.+0.j])
1167 """
1168 a = np.atleast_1d(a)
1170 if a.ndim != 1:
1171 raise ValueError("Input 'a' must be a 1-D array.")
1173 if a.size <= 2:
1174 if a.size == 2:
1175 return np.array([[-(a/a[0])[-1]]])
1176 return np.array([], dtype=a.dtype)
1178 if a[0] == 0.:
1179 raise ValueError('Leading coefficient is zero.')
1181 a = a/a[0]
1182 n = a.size - 1
1183 c = np.zeros((n, n), dtype=a.dtype)
1184 # subdiagonals
1185 c[range(3, n, 2), range(1, n-2, 2)] = 1.
1186 c[range(2, n, 2), range(1, n-1, 2)] = -a[3::2]
1187 # superdiagonals
1188 c[range(0, n-2, 2), range(2, n, 2)] = 1.
1189 c[range(0, n-1, 2), range(1, n, 2)] = -a[2::2]
1190 c[[0, 1], 0] = [-a[1], 1]
1192 return c
1195def convolution_matrix(a, n, mode='full'):
1196 """
1197 Construct a convolution matrix.
1199 Constructs the Toeplitz matrix representing one-dimensional
1200 convolution [1]_. See the notes below for details.
1202 Parameters
1203 ----------
1204 a : (m,) array_like
1205 The 1-D array to convolve.
1206 n : int
1207 The number of columns in the resulting matrix. It gives the length
1208 of the input to be convolved with `a`. This is analogous to the
1209 length of `v` in ``numpy.convolve(a, v)``.
1210 mode : str
1211 This is analogous to `mode` in ``numpy.convolve(v, a, mode)``.
1212 It must be one of ('full', 'valid', 'same').
1213 See below for how `mode` determines the shape of the result.
1215 Returns
1216 -------
1217 A : (k, n) ndarray
1218 The convolution matrix whose row count `k` depends on `mode`::
1220 ======= =========================
1221 mode k
1222 ======= =========================
1223 'full' m + n -1
1224 'same' max(m, n)
1225 'valid' max(m, n) - min(m, n) + 1
1226 ======= =========================
1228 See Also
1229 --------
1230 toeplitz : Toeplitz matrix
1232 Notes
1233 -----
1234 The code::
1236 A = convolution_matrix(a, n, mode)
1238 creates a Toeplitz matrix `A` such that ``A @ v`` is equivalent to
1239 using ``convolve(a, v, mode)``. The returned array always has `n`
1240 columns. The number of rows depends on the specified `mode`, as
1241 explained above.
1243 In the default 'full' mode, the entries of `A` are given by::
1245 A[i, j] == (a[i-j] if (0 <= (i-j) < m) else 0)
1247 where ``m = len(a)``. Suppose, for example, the input array is
1248 ``[x, y, z]``. The convolution matrix has the form::
1250 [x, 0, 0, ..., 0, 0]
1251 [y, x, 0, ..., 0, 0]
1252 [z, y, x, ..., 0, 0]
1253 ...
1254 [0, 0, 0, ..., x, 0]
1255 [0, 0, 0, ..., y, x]
1256 [0, 0, 0, ..., z, y]
1257 [0, 0, 0, ..., 0, z]
1259 In 'valid' mode, the entries of `A` are given by::
1261 A[i, j] == (a[i-j+m-1] if (0 <= (i-j+m-1) < m) else 0)
1263 This corresponds to a matrix whose rows are the subset of those from
1264 the 'full' case where all the coefficients in `a` are contained in the
1265 row. For input ``[x, y, z]``, this array looks like::
1267 [z, y, x, 0, 0, ..., 0, 0, 0]
1268 [0, z, y, x, 0, ..., 0, 0, 0]
1269 [0, 0, z, y, x, ..., 0, 0, 0]
1270 ...
1271 [0, 0, 0, 0, 0, ..., x, 0, 0]
1272 [0, 0, 0, 0, 0, ..., y, x, 0]
1273 [0, 0, 0, 0, 0, ..., z, y, x]
1275 In the 'same' mode, the entries of `A` are given by::
1277 d = (m - 1) // 2
1278 A[i, j] == (a[i-j+d] if (0 <= (i-j+d) < m) else 0)
1280 The typical application of the 'same' mode is when one has a signal of
1281 length `n` (with `n` greater than ``len(a)``), and the desired output
1282 is a filtered signal that is still of length `n`.
1284 For input ``[x, y, z]``, this array looks like::
1286 [y, x, 0, 0, ..., 0, 0, 0]
1287 [z, y, x, 0, ..., 0, 0, 0]
1288 [0, z, y, x, ..., 0, 0, 0]
1289 [0, 0, z, y, ..., 0, 0, 0]
1290 ...
1291 [0, 0, 0, 0, ..., y, x, 0]
1292 [0, 0, 0, 0, ..., z, y, x]
1293 [0, 0, 0, 0, ..., 0, z, y]
1295 .. versionadded:: 1.5.0
1297 References
1298 ----------
1299 .. [1] "Convolution", https://en.wikipedia.org/wiki/Convolution
1301 Examples
1302 --------
1303 >>> from scipy.linalg import convolution_matrix
1304 >>> A = convolution_matrix([-1, 4, -2], 5, mode='same')
1305 >>> A
1306 array([[ 4, -1, 0, 0, 0],
1307 [-2, 4, -1, 0, 0],
1308 [ 0, -2, 4, -1, 0],
1309 [ 0, 0, -2, 4, -1],
1310 [ 0, 0, 0, -2, 4]])
1312 Compare multiplication by `A` with the use of `numpy.convolve`.
1314 >>> x = np.array([1, 2, 0, -3, 0.5])
1315 >>> A @ x
1316 array([ 2. , 6. , -1. , -12.5, 8. ])
1318 Verify that ``A @ x`` produced the same result as applying the
1319 convolution function.
1321 >>> np.convolve([-1, 4, -2], x, mode='same')
1322 array([ 2. , 6. , -1. , -12.5, 8. ])
1324 For comparison to the case ``mode='same'`` shown above, here are the
1325 matrices produced by ``mode='full'`` and ``mode='valid'`` for the
1326 same coefficients and size.
1328 >>> convolution_matrix([-1, 4, -2], 5, mode='full')
1329 array([[-1, 0, 0, 0, 0],
1330 [ 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],
1335 [ 0, 0, 0, 0, -2]])
1337 >>> convolution_matrix([-1, 4, -2], 5, mode='valid')
1338 array([[-2, 4, -1, 0, 0],
1339 [ 0, -2, 4, -1, 0],
1340 [ 0, 0, -2, 4, -1]])
1341 """
1342 if n <= 0:
1343 raise ValueError('n must be a positive integer.')
1345 a = np.asarray(a)
1346 if a.ndim != 1:
1347 raise ValueError('convolution_matrix expects a one-dimensional '
1348 'array as input')
1349 if a.size == 0:
1350 raise ValueError('len(a) must be at least 1.')
1352 if mode not in ('full', 'valid', 'same'):
1353 raise ValueError(
1354 "'mode' argument must be one of ('full', 'valid', 'same')")
1356 # create zero padded versions of the array
1357 az = np.pad(a, (0, n-1), 'constant')
1358 raz = np.pad(a[::-1], (0, n-1), 'constant')
1360 if mode == 'same':
1361 trim = min(n, len(a)) - 1
1362 tb = trim//2
1363 te = trim - tb
1364 col0 = az[tb:len(az)-te]
1365 row0 = raz[-n-tb:len(raz)-tb]
1366 elif mode == 'valid':
1367 tb = min(n, len(a)) - 1
1368 te = tb
1369 col0 = az[tb:len(az)-te]
1370 row0 = raz[-n-tb:len(raz)-tb]
1371 else: # 'full'
1372 col0 = az
1373 row0 = raz[-n:]
1374 return toeplitz(col0, row0)