Coverage for /usr/lib/python3/dist-packages/scipy/linalg/_matfuncs.py: 15%
212 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-06-14 15:25 +0200
1#
2# Author: Travis Oliphant, March 2002
3#
4from itertools import product
6import numpy as np
7from numpy import (Inf, dot, diag, prod, logical_not, ravel, transpose,
8 conjugate, absolute, amax, sign, isfinite, triu)
9from numpy.lib.scimath import sqrt as csqrt
11# Local imports
12from scipy.linalg import LinAlgError, bandwidth
13from ._misc import norm
14from ._basic import solve, inv
15from ._decomp_svd import svd
16from ._decomp_schur import schur, rsf2csf
17from ._expm_frechet import expm_frechet, expm_cond
18from ._matfuncs_sqrtm import sqrtm
19from ._matfuncs_expm import pick_pade_structure, pade_UV_calc
21__all__ = ['expm', 'cosm', 'sinm', 'tanm', 'coshm', 'sinhm', 'tanhm', 'logm',
22 'funm', 'signm', 'sqrtm', 'fractional_matrix_power', 'expm_frechet',
23 'expm_cond', 'khatri_rao']
25eps = np.finfo('d').eps
26feps = np.finfo('f').eps
28_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
31###############################################################################
32# Utility functions.
35def _asarray_square(A):
36 """
37 Wraps asarray with the extra requirement that the input be a square matrix.
39 The motivation is that the matfuncs module has real functions that have
40 been lifted to square matrix functions.
42 Parameters
43 ----------
44 A : array_like
45 A square matrix.
47 Returns
48 -------
49 out : ndarray
50 An ndarray copy or view or other representation of A.
52 """
53 A = np.asarray(A)
54 if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
55 raise ValueError('expected square array_like input')
56 return A
59def _maybe_real(A, B, tol=None):
60 """
61 Return either B or the real part of B, depending on properties of A and B.
63 The motivation is that B has been computed as a complicated function of A,
64 and B may be perturbed by negligible imaginary components.
65 If A is real and B is complex with small imaginary components,
66 then return a real copy of B. The assumption in that case would be that
67 the imaginary components of B are numerical artifacts.
69 Parameters
70 ----------
71 A : ndarray
72 Input array whose type is to be checked as real vs. complex.
73 B : ndarray
74 Array to be returned, possibly without its imaginary part.
75 tol : float
76 Absolute tolerance.
78 Returns
79 -------
80 out : real or complex array
81 Either the input array B or only the real part of the input array B.
83 """
84 # Note that booleans and integers compare as real.
85 if np.isrealobj(A) and np.iscomplexobj(B):
86 if tol is None:
87 tol = {0:feps*1e3, 1:eps*1e6}[_array_precision[B.dtype.char]]
88 if np.allclose(B.imag, 0.0, atol=tol):
89 B = B.real
90 return B
93###############################################################################
94# Matrix functions.
97def fractional_matrix_power(A, t):
98 """
99 Compute the fractional power of a matrix.
101 Proceeds according to the discussion in section (6) of [1]_.
103 Parameters
104 ----------
105 A : (N, N) array_like
106 Matrix whose fractional power to evaluate.
107 t : float
108 Fractional power.
110 Returns
111 -------
112 X : (N, N) array_like
113 The fractional power of the matrix.
115 References
116 ----------
117 .. [1] Nicholas J. Higham and Lijing lin (2011)
118 "A Schur-Pade Algorithm for Fractional Powers of a Matrix."
119 SIAM Journal on Matrix Analysis and Applications,
120 32 (3). pp. 1056-1078. ISSN 0895-4798
122 Examples
123 --------
124 >>> import numpy as np
125 >>> from scipy.linalg import fractional_matrix_power
126 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
127 >>> b = fractional_matrix_power(a, 0.5)
128 >>> b
129 array([[ 0.75592895, 1.13389342],
130 [ 0.37796447, 1.88982237]])
131 >>> np.dot(b, b) # Verify square root
132 array([[ 1., 3.],
133 [ 1., 4.]])
135 """
136 # This fixes some issue with imports;
137 # this function calls onenormest which is in scipy.sparse.
138 A = _asarray_square(A)
139 import scipy.linalg._matfuncs_inv_ssq
140 return scipy.linalg._matfuncs_inv_ssq._fractional_matrix_power(A, t)
143def logm(A, disp=True):
144 """
145 Compute matrix logarithm.
147 The matrix logarithm is the inverse of
148 expm: expm(logm(`A`)) == `A`
150 Parameters
151 ----------
152 A : (N, N) array_like
153 Matrix whose logarithm to evaluate
154 disp : bool, optional
155 Print warning if error in the result is estimated large
156 instead of returning estimated error. (Default: True)
158 Returns
159 -------
160 logm : (N, N) ndarray
161 Matrix logarithm of `A`
162 errest : float
163 (if disp == False)
165 1-norm of the estimated error, ||err||_1 / ||A||_1
167 References
168 ----------
169 .. [1] Awad H. Al-Mohy and Nicholas J. Higham (2012)
170 "Improved Inverse Scaling and Squaring Algorithms
171 for the Matrix Logarithm."
172 SIAM Journal on Scientific Computing, 34 (4). C152-C169.
173 ISSN 1095-7197
175 .. [2] Nicholas J. Higham (2008)
176 "Functions of Matrices: Theory and Computation"
177 ISBN 978-0-898716-46-7
179 .. [3] Nicholas J. Higham and Lijing lin (2011)
180 "A Schur-Pade Algorithm for Fractional Powers of a Matrix."
181 SIAM Journal on Matrix Analysis and Applications,
182 32 (3). pp. 1056-1078. ISSN 0895-4798
184 Examples
185 --------
186 >>> import numpy as np
187 >>> from scipy.linalg import logm, expm
188 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
189 >>> b = logm(a)
190 >>> b
191 array([[-1.02571087, 2.05142174],
192 [ 0.68380725, 1.02571087]])
193 >>> expm(b) # Verify expm(logm(a)) returns a
194 array([[ 1., 3.],
195 [ 1., 4.]])
197 """
198 A = _asarray_square(A)
199 # Avoid circular import ... this is OK, right?
200 import scipy.linalg._matfuncs_inv_ssq
201 F = scipy.linalg._matfuncs_inv_ssq._logm(A)
202 F = _maybe_real(A, F)
203 errtol = 1000*eps
204 #TODO use a better error approximation
205 errest = norm(expm(F)-A,1) / norm(A,1)
206 if disp:
207 if not isfinite(errest) or errest >= errtol:
208 print("logm result may be inaccurate, approximate err =", errest)
209 return F
210 else:
211 return F, errest
214def expm(A):
215 """Compute the matrix exponential of an array.
217 Parameters
218 ----------
219 A : ndarray
220 Input with last two dimensions are square ``(..., n, n)``.
222 Returns
223 -------
224 eA : ndarray
225 The resulting matrix exponential with the same shape of ``A``
227 Notes
228 -----
229 Implements the algorithm given in [1], which is essentially a Pade
230 approximation with a variable order that is decided based on the array
231 data.
233 For input with size ``n``, the memory usage is in the worst case in the
234 order of ``8*(n**2)``. If the input data is not of single and double
235 precision of real and complex dtypes, it is copied to a new array.
237 For cases ``n >= 400``, the exact 1-norm computation cost, breaks even with
238 1-norm estimation and from that point on the estimation scheme given in
239 [2] is used to decide on the approximation order.
241 References
242 ----------
243 .. [1] Awad H. Al-Mohy and Nicholas J. Higham, (2009), "A New Scaling
244 and Squaring Algorithm for the Matrix Exponential", SIAM J. Matrix
245 Anal. Appl. 31(3):970-989, :doi:`10.1137/09074721X`
247 .. [2] Nicholas J. Higham and Francoise Tisseur (2000), "A Block Algorithm
248 for Matrix 1-Norm Estimation, with an Application to 1-Norm
249 Pseudospectra." SIAM J. Matrix Anal. Appl. 21(4):1185-1201,
250 :doi:`10.1137/S0895479899356080`
252 Examples
253 --------
254 >>> import numpy as np
255 >>> from scipy.linalg import expm, sinm, cosm
257 Matrix version of the formula exp(0) = 1:
259 >>> expm(np.zeros((3, 2, 2)))
260 array([[[1., 0.],
261 [0., 1.]],
262 <BLANKLINE>
263 [[1., 0.],
264 [0., 1.]],
265 <BLANKLINE>
266 [[1., 0.],
267 [0., 1.]]])
269 Euler's identity (exp(i*theta) = cos(theta) + i*sin(theta))
270 applied to a matrix:
272 >>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
273 >>> expm(1j*a)
274 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
275 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
276 >>> cosm(a) + 1j*sinm(a)
277 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
278 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
280 """
281 a = np.asarray(A)
282 if a.size == 1 and a.ndim < 2:
283 return np.array([[np.exp(a.item())]])
285 if a.ndim < 2:
286 raise LinAlgError('The input array must be at least two-dimensional')
287 if a.shape[-1] != a.shape[-2]:
288 raise LinAlgError('Last 2 dimensions of the array must be square')
289 n = a.shape[-1]
290 # Empty array
291 if min(*a.shape) == 0:
292 return np.empty_like(a)
294 # Scalar case
295 if a.shape[-2:] == (1, 1):
296 return np.exp(a)
298 if not np.issubdtype(a.dtype, np.inexact):
299 a = a.astype(float)
300 elif a.dtype == np.float16:
301 a = a.astype(np.float32)
303 # Explicit formula for 2x2 case, formula (2.2) in [1]
304 # without Kahan's method numerical instabilities can occur.
305 if a.shape[-2:] == (2, 2):
306 a1, a2, a3, a4 = (a[..., [0], [0]],
307 a[..., [0], [1]],
308 a[..., [1], [0]],
309 a[..., [1], [1]])
310 mu = csqrt((a1-a4)**2 + 4*a2*a3)/2. # csqrt slow but handles neg.vals
312 eApD2 = np.exp((a1+a4)/2.)
313 AmD2 = (a1 - a4)/2.
314 coshMu = np.cosh(mu)
315 sinchMu = np.ones_like(coshMu)
316 mask = mu != 0
317 sinchMu[mask] = np.sinh(mu[mask]) / mu[mask]
318 eA = np.empty((a.shape), dtype=mu.dtype)
319 eA[..., [0], [0]] = eApD2 * (coshMu + AmD2*sinchMu)
320 eA[..., [0], [1]] = eApD2 * a2 * sinchMu
321 eA[..., [1], [0]] = eApD2 * a3 * sinchMu
322 eA[..., [1], [1]] = eApD2 * (coshMu - AmD2*sinchMu)
323 if np.isrealobj(a):
324 return eA.real
325 return eA
327 # larger problem with unspecified stacked dimensions.
328 n = a.shape[-1]
329 eA = np.empty(a.shape, dtype=a.dtype)
330 # working memory to hold intermediate arrays
331 Am = np.empty((5, n, n), dtype=a.dtype)
333 # Main loop to go through the slices of an ndarray and passing to expm
334 for ind in product(*[range(x) for x in a.shape[:-2]]):
335 aw = a[ind]
337 lu = bandwidth(aw)
338 if not any(lu): # a is diagonal?
339 eA[ind] = np.diag(np.exp(np.diag(aw)))
340 continue
342 # Generic/triangular case; copy the slice into scratch and send.
343 # Am will be mutated by pick_pade_structure
344 Am[0, :, :] = aw
345 m, s = pick_pade_structure(Am)
347 if s != 0: # scaling needed
348 Am[:4] *= [[[2**(-s)]], [[4**(-s)]], [[16**(-s)]], [[64**(-s)]]]
350 pade_UV_calc(Am, n, m)
351 eAw = Am[0]
353 if s != 0: # squaring needed
355 if (lu[1] == 0) or (lu[0] == 0): # lower/upper triangular
356 # This branch implements Code Fragment 2.1 of [1]
358 diag_aw = np.diag(aw)
359 # einsum returns a writable view
360 np.einsum('ii->i', eAw)[:] = np.exp(diag_aw * 2**(-s))
361 # super/sub diagonal
362 sd = np.diag(aw, k=-1 if lu[1] == 0 else 1)
364 for i in range(s-1, -1, -1):
365 eAw = eAw @ eAw
367 # diagonal
368 np.einsum('ii->i', eAw)[:] = np.exp(diag_aw * 2.**(-i))
369 exp_sd = _exp_sinch(diag_aw * (2.**(-i))) * (sd * 2**(-i))
370 if lu[1] == 0: # lower
371 np.einsum('ii->i', eAw[1:, :-1])[:] = exp_sd
372 else: # upper
373 np.einsum('ii->i', eAw[:-1, 1:])[:] = exp_sd
375 else: # generic
376 for _ in range(s):
377 eAw = eAw @ eAw
379 # Zero out the entries from np.empty in case of triangular input
380 if (lu[0] == 0) or (lu[1] == 0):
381 eA[ind] = np.triu(eAw) if lu[0] == 0 else np.tril(eAw)
382 else:
383 eA[ind] = eAw
385 return eA
388def _exp_sinch(x):
389 # Higham's formula (10.42), might overflow, see GH-11839
390 lexp_diff = np.diff(np.exp(x))
391 l_diff = np.diff(x)
392 mask_z = l_diff == 0.
393 lexp_diff[~mask_z] /= l_diff[~mask_z]
394 lexp_diff[mask_z] = np.exp(x[:-1][mask_z])
395 return lexp_diff
398def cosm(A):
399 """
400 Compute the matrix cosine.
402 This routine uses expm to compute the matrix exponentials.
404 Parameters
405 ----------
406 A : (N, N) array_like
407 Input array
409 Returns
410 -------
411 cosm : (N, N) ndarray
412 Matrix cosine of A
414 Examples
415 --------
416 >>> import numpy as np
417 >>> from scipy.linalg import expm, sinm, cosm
419 Euler's identity (exp(i*theta) = cos(theta) + i*sin(theta))
420 applied to a matrix:
422 >>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
423 >>> expm(1j*a)
424 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
425 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
426 >>> cosm(a) + 1j*sinm(a)
427 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
428 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
430 """
431 A = _asarray_square(A)
432 if np.iscomplexobj(A):
433 return 0.5*(expm(1j*A) + expm(-1j*A))
434 else:
435 return expm(1j*A).real
438def sinm(A):
439 """
440 Compute the matrix sine.
442 This routine uses expm to compute the matrix exponentials.
444 Parameters
445 ----------
446 A : (N, N) array_like
447 Input array.
449 Returns
450 -------
451 sinm : (N, N) ndarray
452 Matrix sine of `A`
454 Examples
455 --------
456 >>> import numpy as np
457 >>> from scipy.linalg import expm, sinm, cosm
459 Euler's identity (exp(i*theta) = cos(theta) + i*sin(theta))
460 applied to a matrix:
462 >>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
463 >>> expm(1j*a)
464 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
465 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
466 >>> cosm(a) + 1j*sinm(a)
467 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
468 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
470 """
471 A = _asarray_square(A)
472 if np.iscomplexobj(A):
473 return -0.5j*(expm(1j*A) - expm(-1j*A))
474 else:
475 return expm(1j*A).imag
478def tanm(A):
479 """
480 Compute the matrix tangent.
482 This routine uses expm to compute the matrix exponentials.
484 Parameters
485 ----------
486 A : (N, N) array_like
487 Input array.
489 Returns
490 -------
491 tanm : (N, N) ndarray
492 Matrix tangent of `A`
494 Examples
495 --------
496 >>> import numpy as np
497 >>> from scipy.linalg import tanm, sinm, cosm
498 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
499 >>> t = tanm(a)
500 >>> t
501 array([[ -2.00876993, -8.41880636],
502 [ -2.80626879, -10.42757629]])
504 Verify tanm(a) = sinm(a).dot(inv(cosm(a)))
506 >>> s = sinm(a)
507 >>> c = cosm(a)
508 >>> s.dot(np.linalg.inv(c))
509 array([[ -2.00876993, -8.41880636],
510 [ -2.80626879, -10.42757629]])
512 """
513 A = _asarray_square(A)
514 return _maybe_real(A, solve(cosm(A), sinm(A)))
517def coshm(A):
518 """
519 Compute the hyperbolic matrix cosine.
521 This routine uses expm to compute the matrix exponentials.
523 Parameters
524 ----------
525 A : (N, N) array_like
526 Input array.
528 Returns
529 -------
530 coshm : (N, N) ndarray
531 Hyperbolic matrix cosine of `A`
533 Examples
534 --------
535 >>> import numpy as np
536 >>> from scipy.linalg import tanhm, sinhm, coshm
537 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
538 >>> c = coshm(a)
539 >>> c
540 array([[ 11.24592233, 38.76236492],
541 [ 12.92078831, 50.00828725]])
543 Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))
545 >>> t = tanhm(a)
546 >>> s = sinhm(a)
547 >>> t - s.dot(np.linalg.inv(c))
548 array([[ 2.72004641e-15, 4.55191440e-15],
549 [ 0.00000000e+00, -5.55111512e-16]])
551 """
552 A = _asarray_square(A)
553 return _maybe_real(A, 0.5 * (expm(A) + expm(-A)))
556def sinhm(A):
557 """
558 Compute the hyperbolic matrix sine.
560 This routine uses expm to compute the matrix exponentials.
562 Parameters
563 ----------
564 A : (N, N) array_like
565 Input array.
567 Returns
568 -------
569 sinhm : (N, N) ndarray
570 Hyperbolic matrix sine of `A`
572 Examples
573 --------
574 >>> import numpy as np
575 >>> from scipy.linalg import tanhm, sinhm, coshm
576 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
577 >>> s = sinhm(a)
578 >>> s
579 array([[ 10.57300653, 39.28826594],
580 [ 13.09608865, 49.86127247]])
582 Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))
584 >>> t = tanhm(a)
585 >>> c = coshm(a)
586 >>> t - s.dot(np.linalg.inv(c))
587 array([[ 2.72004641e-15, 4.55191440e-15],
588 [ 0.00000000e+00, -5.55111512e-16]])
590 """
591 A = _asarray_square(A)
592 return _maybe_real(A, 0.5 * (expm(A) - expm(-A)))
595def tanhm(A):
596 """
597 Compute the hyperbolic matrix tangent.
599 This routine uses expm to compute the matrix exponentials.
601 Parameters
602 ----------
603 A : (N, N) array_like
604 Input array
606 Returns
607 -------
608 tanhm : (N, N) ndarray
609 Hyperbolic matrix tangent of `A`
611 Examples
612 --------
613 >>> import numpy as np
614 >>> from scipy.linalg import tanhm, sinhm, coshm
615 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
616 >>> t = tanhm(a)
617 >>> t
618 array([[ 0.3428582 , 0.51987926],
619 [ 0.17329309, 0.86273746]])
621 Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))
623 >>> s = sinhm(a)
624 >>> c = coshm(a)
625 >>> t - s.dot(np.linalg.inv(c))
626 array([[ 2.72004641e-15, 4.55191440e-15],
627 [ 0.00000000e+00, -5.55111512e-16]])
629 """
630 A = _asarray_square(A)
631 return _maybe_real(A, solve(coshm(A), sinhm(A)))
634def funm(A, func, disp=True):
635 """
636 Evaluate a matrix function specified by a callable.
638 Returns the value of matrix-valued function ``f`` at `A`. The
639 function ``f`` is an extension of the scalar-valued function `func`
640 to matrices.
642 Parameters
643 ----------
644 A : (N, N) array_like
645 Matrix at which to evaluate the function
646 func : callable
647 Callable object that evaluates a scalar function f.
648 Must be vectorized (eg. using vectorize).
649 disp : bool, optional
650 Print warning if error in the result is estimated large
651 instead of returning estimated error. (Default: True)
653 Returns
654 -------
655 funm : (N, N) ndarray
656 Value of the matrix function specified by func evaluated at `A`
657 errest : float
658 (if disp == False)
660 1-norm of the estimated error, ||err||_1 / ||A||_1
662 Notes
663 -----
664 This function implements the general algorithm based on Schur decomposition
665 (Algorithm 9.1.1. in [1]_).
667 If the input matrix is known to be diagonalizable, then relying on the
668 eigendecomposition is likely to be faster. For example, if your matrix is
669 Hermitian, you can do
671 >>> from scipy.linalg import eigh
672 >>> def funm_herm(a, func, check_finite=False):
673 ... w, v = eigh(a, check_finite=check_finite)
674 ... ## if you further know that your matrix is positive semidefinite,
675 ... ## you can optionally guard against precision errors by doing
676 ... # w = np.maximum(w, 0)
677 ... w = func(w)
678 ... return (v * w).dot(v.conj().T)
680 References
681 ----------
682 .. [1] Gene H. Golub, Charles F. van Loan, Matrix Computations 4th ed.
684 Examples
685 --------
686 >>> import numpy as np
687 >>> from scipy.linalg import funm
688 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
689 >>> funm(a, lambda x: x*x)
690 array([[ 4., 15.],
691 [ 5., 19.]])
692 >>> a.dot(a)
693 array([[ 4., 15.],
694 [ 5., 19.]])
696 """
697 A = _asarray_square(A)
698 # Perform Shur decomposition (lapack ?gees)
699 T, Z = schur(A)
700 T, Z = rsf2csf(T,Z)
701 n,n = T.shape
702 F = diag(func(diag(T))) # apply function to diagonal elements
703 F = F.astype(T.dtype.char) # e.g., when F is real but T is complex
705 minden = abs(T[0,0])
707 # implement Algorithm 11.1.1 from Golub and Van Loan
708 # "matrix Computations."
709 for p in range(1,n):
710 for i in range(1,n-p+1):
711 j = i + p
712 s = T[i-1,j-1] * (F[j-1,j-1] - F[i-1,i-1])
713 ksl = slice(i,j-1)
714 val = dot(T[i-1,ksl],F[ksl,j-1]) - dot(F[i-1,ksl],T[ksl,j-1])
715 s = s + val
716 den = T[j-1,j-1] - T[i-1,i-1]
717 if den != 0.0:
718 s = s / den
719 F[i-1,j-1] = s
720 minden = min(minden,abs(den))
722 F = dot(dot(Z, F), transpose(conjugate(Z)))
723 F = _maybe_real(A, F)
725 tol = {0:feps, 1:eps}[_array_precision[F.dtype.char]]
726 if minden == 0.0:
727 minden = tol
728 err = min(1, max(tol,(tol/minden)*norm(triu(T,1),1)))
729 if prod(ravel(logical_not(isfinite(F))),axis=0):
730 err = Inf
731 if disp:
732 if err > 1000*tol:
733 print("funm result may be inaccurate, approximate err =", err)
734 return F
735 else:
736 return F, err
739def signm(A, disp=True):
740 """
741 Matrix sign function.
743 Extension of the scalar sign(x) to matrices.
745 Parameters
746 ----------
747 A : (N, N) array_like
748 Matrix at which to evaluate the sign function
749 disp : bool, optional
750 Print warning if error in the result is estimated large
751 instead of returning estimated error. (Default: True)
753 Returns
754 -------
755 signm : (N, N) ndarray
756 Value of the sign function at `A`
757 errest : float
758 (if disp == False)
760 1-norm of the estimated error, ||err||_1 / ||A||_1
762 Examples
763 --------
764 >>> from scipy.linalg import signm, eigvals
765 >>> a = [[1,2,3], [1,2,1], [1,1,1]]
766 >>> eigvals(a)
767 array([ 4.12488542+0.j, -0.76155718+0.j, 0.63667176+0.j])
768 >>> eigvals(signm(a))
769 array([-1.+0.j, 1.+0.j, 1.+0.j])
771 """
772 A = _asarray_square(A)
774 def rounded_sign(x):
775 rx = np.real(x)
776 if rx.dtype.char == 'f':
777 c = 1e3*feps*amax(x)
778 else:
779 c = 1e3*eps*amax(x)
780 return sign((absolute(rx) > c) * rx)
781 result, errest = funm(A, rounded_sign, disp=0)
782 errtol = {0:1e3*feps, 1:1e3*eps}[_array_precision[result.dtype.char]]
783 if errest < errtol:
784 return result
786 # Handle signm of defective matrices:
788 # See "E.D.Denman and J.Leyva-Ramos, Appl.Math.Comp.,
789 # 8:237-250,1981" for how to improve the following (currently a
790 # rather naive) iteration process:
792 # a = result # sometimes iteration converges faster but where??
794 # Shifting to avoid zero eigenvalues. How to ensure that shifting does
795 # not change the spectrum too much?
796 vals = svd(A, compute_uv=False)
797 max_sv = np.amax(vals)
798 # min_nonzero_sv = vals[(vals>max_sv*errtol).tolist().count(1)-1]
799 # c = 0.5/min_nonzero_sv
800 c = 0.5/max_sv
801 S0 = A + c*np.identity(A.shape[0])
802 prev_errest = errest
803 for i in range(100):
804 iS0 = inv(S0)
805 S0 = 0.5*(S0 + iS0)
806 Pp = 0.5*(dot(S0,S0)+S0)
807 errest = norm(dot(Pp,Pp)-Pp,1)
808 if errest < errtol or prev_errest == errest:
809 break
810 prev_errest = errest
811 if disp:
812 if not isfinite(errest) or errest >= errtol:
813 print("signm result may be inaccurate, approximate err =", errest)
814 return S0
815 else:
816 return S0, errest
819def khatri_rao(a, b):
820 r"""
821 Khatri-rao product
823 A column-wise Kronecker product of two matrices
825 Parameters
826 ----------
827 a : (n, k) array_like
828 Input array
829 b : (m, k) array_like
830 Input array
832 Returns
833 -------
834 c: (n*m, k) ndarray
835 Khatri-rao product of `a` and `b`.
837 See Also
838 --------
839 kron : Kronecker product
841 Notes
842 -----
843 The mathematical definition of the Khatri-Rao product is:
845 .. math::
847 (A_{ij} \bigotimes B_{ij})_{ij}
849 which is the Kronecker product of every column of A and B, e.g.::
851 c = np.vstack([np.kron(a[:, k], b[:, k]) for k in range(b.shape[1])]).T
853 Examples
854 --------
855 >>> import numpy as np
856 >>> from scipy import linalg
857 >>> a = np.array([[1, 2, 3], [4, 5, 6]])
858 >>> b = np.array([[3, 4, 5], [6, 7, 8], [2, 3, 9]])
859 >>> linalg.khatri_rao(a, b)
860 array([[ 3, 8, 15],
861 [ 6, 14, 24],
862 [ 2, 6, 27],
863 [12, 20, 30],
864 [24, 35, 48],
865 [ 8, 15, 54]])
867 """
868 a = np.asarray(a)
869 b = np.asarray(b)
871 if not (a.ndim == 2 and b.ndim == 2):
872 raise ValueError("The both arrays should be 2-dimensional.")
874 if not a.shape[1] == b.shape[1]:
875 raise ValueError("The number of columns for both arrays "
876 "should be equal.")
878 # c = np.vstack([np.kron(a[:, k], b[:, k]) for k in range(b.shape[1])]).T
879 c = a[..., :, np.newaxis, :] * b[..., np.newaxis, :, :]
880 return c.reshape((-1,) + c.shape[2:])