Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/linalg/_matfuncs.py: 16%
195 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-03-22 06:44 +0000
1#
2# Author: Travis Oliphant, March 2002
3#
4from itertools import product
6import numpy as np
7from numpy import (dot, diag, prod, logical_not, ravel, transpose,
8 conjugate, absolute, amax, sign, isfinite, triu)
10# Local imports
11from scipy.linalg import LinAlgError, bandwidth
12from ._misc import norm
13from ._basic import solve, inv
14from ._decomp_svd import svd
15from ._decomp_schur import schur, rsf2csf
16from ._expm_frechet import expm_frechet, expm_cond
17from ._matfuncs_sqrtm import sqrtm
18from ._matfuncs_expm import pick_pade_structure, pade_UV_calc
20# deprecated imports to be removed in SciPy 1.13.0
21from numpy import single # noqa: F401
23__all__ = ['expm', 'cosm', 'sinm', 'tanm', 'coshm', 'sinhm', 'tanhm', 'logm',
24 'funm', 'signm', 'sqrtm', 'fractional_matrix_power', 'expm_frechet',
25 'expm_cond', 'khatri_rao']
27eps = np.finfo('d').eps
28feps = np.finfo('f').eps
30_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
33###############################################################################
34# Utility functions.
37def _asarray_square(A):
38 """
39 Wraps asarray with the extra requirement that the input be a square matrix.
41 The motivation is that the matfuncs module has real functions that have
42 been lifted to square matrix functions.
44 Parameters
45 ----------
46 A : array_like
47 A square matrix.
49 Returns
50 -------
51 out : ndarray
52 An ndarray copy or view or other representation of A.
54 """
55 A = np.asarray(A)
56 if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
57 raise ValueError('expected square array_like input')
58 return A
61def _maybe_real(A, B, tol=None):
62 """
63 Return either B or the real part of B, depending on properties of A and B.
65 The motivation is that B has been computed as a complicated function of A,
66 and B may be perturbed by negligible imaginary components.
67 If A is real and B is complex with small imaginary components,
68 then return a real copy of B. The assumption in that case would be that
69 the imaginary components of B are numerical artifacts.
71 Parameters
72 ----------
73 A : ndarray
74 Input array whose type is to be checked as real vs. complex.
75 B : ndarray
76 Array to be returned, possibly without its imaginary part.
77 tol : float
78 Absolute tolerance.
80 Returns
81 -------
82 out : real or complex array
83 Either the input array B or only the real part of the input array B.
85 """
86 # Note that booleans and integers compare as real.
87 if np.isrealobj(A) and np.iscomplexobj(B):
88 if tol is None:
89 tol = {0: feps*1e3, 1: eps*1e6}[_array_precision[B.dtype.char]]
90 if np.allclose(B.imag, 0.0, atol=tol):
91 B = B.real
92 return B
95###############################################################################
96# Matrix functions.
99def fractional_matrix_power(A, t):
100 """
101 Compute the fractional power of a matrix.
103 Proceeds according to the discussion in section (6) of [1]_.
105 Parameters
106 ----------
107 A : (N, N) array_like
108 Matrix whose fractional power to evaluate.
109 t : float
110 Fractional power.
112 Returns
113 -------
114 X : (N, N) array_like
115 The fractional power of the matrix.
117 References
118 ----------
119 .. [1] Nicholas J. Higham and Lijing lin (2011)
120 "A Schur-Pade Algorithm for Fractional Powers of a Matrix."
121 SIAM Journal on Matrix Analysis and Applications,
122 32 (3). pp. 1056-1078. ISSN 0895-4798
124 Examples
125 --------
126 >>> import numpy as np
127 >>> from scipy.linalg import fractional_matrix_power
128 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
129 >>> b = fractional_matrix_power(a, 0.5)
130 >>> b
131 array([[ 0.75592895, 1.13389342],
132 [ 0.37796447, 1.88982237]])
133 >>> np.dot(b, b) # Verify square root
134 array([[ 1., 3.],
135 [ 1., 4.]])
137 """
138 # This fixes some issue with imports;
139 # this function calls onenormest which is in scipy.sparse.
140 A = _asarray_square(A)
141 import scipy.linalg._matfuncs_inv_ssq
142 return scipy.linalg._matfuncs_inv_ssq._fractional_matrix_power(A, t)
145def logm(A, disp=True):
146 """
147 Compute matrix logarithm.
149 The matrix logarithm is the inverse of
150 expm: expm(logm(`A`)) == `A`
152 Parameters
153 ----------
154 A : (N, N) array_like
155 Matrix whose logarithm to evaluate
156 disp : bool, optional
157 Print warning if error in the result is estimated large
158 instead of returning estimated error. (Default: True)
160 Returns
161 -------
162 logm : (N, N) ndarray
163 Matrix logarithm of `A`
164 errest : float
165 (if disp == False)
167 1-norm of the estimated error, ||err||_1 / ||A||_1
169 References
170 ----------
171 .. [1] Awad H. Al-Mohy and Nicholas J. Higham (2012)
172 "Improved Inverse Scaling and Squaring Algorithms
173 for the Matrix Logarithm."
174 SIAM Journal on Scientific Computing, 34 (4). C152-C169.
175 ISSN 1095-7197
177 .. [2] Nicholas J. Higham (2008)
178 "Functions of Matrices: Theory and Computation"
179 ISBN 978-0-898716-46-7
181 .. [3] Nicholas J. Higham and Lijing lin (2011)
182 "A Schur-Pade Algorithm for Fractional Powers of a Matrix."
183 SIAM Journal on Matrix Analysis and Applications,
184 32 (3). pp. 1056-1078. ISSN 0895-4798
186 Examples
187 --------
188 >>> import numpy as np
189 >>> from scipy.linalg import logm, expm
190 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
191 >>> b = logm(a)
192 >>> b
193 array([[-1.02571087, 2.05142174],
194 [ 0.68380725, 1.02571087]])
195 >>> expm(b) # Verify expm(logm(a)) returns a
196 array([[ 1., 3.],
197 [ 1., 4.]])
199 """
200 A = _asarray_square(A)
201 # Avoid circular import ... this is OK, right?
202 import scipy.linalg._matfuncs_inv_ssq
203 F = scipy.linalg._matfuncs_inv_ssq._logm(A)
204 F = _maybe_real(A, F)
205 errtol = 1000*eps
206 # TODO use a better error approximation
207 errest = norm(expm(F)-A, 1) / norm(A, 1)
208 if disp:
209 if not isfinite(errest) or errest >= errtol:
210 print("logm result may be inaccurate, approximate err =", errest)
211 return F
212 else:
213 return F, errest
216def expm(A):
217 """Compute the matrix exponential of an array.
219 Parameters
220 ----------
221 A : ndarray
222 Input with last two dimensions are square ``(..., n, n)``.
224 Returns
225 -------
226 eA : ndarray
227 The resulting matrix exponential with the same shape of ``A``
229 Notes
230 -----
231 Implements the algorithm given in [1], which is essentially a Pade
232 approximation with a variable order that is decided based on the array
233 data.
235 For input with size ``n``, the memory usage is in the worst case in the
236 order of ``8*(n**2)``. If the input data is not of single and double
237 precision of real and complex dtypes, it is copied to a new array.
239 For cases ``n >= 400``, the exact 1-norm computation cost, breaks even with
240 1-norm estimation and from that point on the estimation scheme given in
241 [2] is used to decide on the approximation order.
243 References
244 ----------
245 .. [1] Awad H. Al-Mohy and Nicholas J. Higham, (2009), "A New Scaling
246 and Squaring Algorithm for the Matrix Exponential", SIAM J. Matrix
247 Anal. Appl. 31(3):970-989, :doi:`10.1137/09074721X`
249 .. [2] Nicholas J. Higham and Francoise Tisseur (2000), "A Block Algorithm
250 for Matrix 1-Norm Estimation, with an Application to 1-Norm
251 Pseudospectra." SIAM J. Matrix Anal. Appl. 21(4):1185-1201,
252 :doi:`10.1137/S0895479899356080`
254 Examples
255 --------
256 >>> import numpy as np
257 >>> from scipy.linalg import expm, sinm, cosm
259 Matrix version of the formula exp(0) = 1:
261 >>> expm(np.zeros((3, 2, 2)))
262 array([[[1., 0.],
263 [0., 1.]],
264 <BLANKLINE>
265 [[1., 0.],
266 [0., 1.]],
267 <BLANKLINE>
268 [[1., 0.],
269 [0., 1.]]])
271 Euler's identity (exp(i*theta) = cos(theta) + i*sin(theta))
272 applied to a matrix:
274 >>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
275 >>> expm(1j*a)
276 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
277 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
278 >>> cosm(a) + 1j*sinm(a)
279 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
280 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
282 """
283 a = np.asarray(A)
284 if a.size == 1 and a.ndim < 2:
285 return np.array([[np.exp(a.item())]])
287 if a.ndim < 2:
288 raise LinAlgError('The input array must be at least two-dimensional')
289 if a.shape[-1] != a.shape[-2]:
290 raise LinAlgError('Last 2 dimensions of the array must be square')
291 n = a.shape[-1]
292 # Empty array
293 if min(*a.shape) == 0:
294 return np.empty_like(a)
296 # Scalar case
297 if a.shape[-2:] == (1, 1):
298 return np.exp(a)
300 if not np.issubdtype(a.dtype, np.inexact):
301 a = a.astype(np.float64)
302 elif a.dtype == np.float16:
303 a = a.astype(np.float32)
305 # An explicit formula for 2x2 case exists (formula (2.2) in [1]). However, without
306 # Kahan's method, numerical instabilities can occur (See gh-19584). Hence removed
307 # here until we have a more stable implementation.
309 n = a.shape[-1]
310 eA = np.empty(a.shape, dtype=a.dtype)
311 # working memory to hold intermediate arrays
312 Am = np.empty((5, n, n), dtype=a.dtype)
314 # Main loop to go through the slices of an ndarray and passing to expm
315 for ind in product(*[range(x) for x in a.shape[:-2]]):
316 aw = a[ind]
318 lu = bandwidth(aw)
319 if not any(lu): # a is diagonal?
320 eA[ind] = np.diag(np.exp(np.diag(aw)))
321 continue
323 # Generic/triangular case; copy the slice into scratch and send.
324 # Am will be mutated by pick_pade_structure
325 Am[0, :, :] = aw
326 m, s = pick_pade_structure(Am)
328 if s != 0: # scaling needed
329 Am[:4] *= [[[2**(-s)]], [[4**(-s)]], [[16**(-s)]], [[64**(-s)]]]
331 pade_UV_calc(Am, n, m)
332 eAw = Am[0]
334 if s != 0: # squaring needed
336 if (lu[1] == 0) or (lu[0] == 0): # lower/upper triangular
337 # This branch implements Code Fragment 2.1 of [1]
339 diag_aw = np.diag(aw)
340 # einsum returns a writable view
341 np.einsum('ii->i', eAw)[:] = np.exp(diag_aw * 2**(-s))
342 # super/sub diagonal
343 sd = np.diag(aw, k=-1 if lu[1] == 0 else 1)
345 for i in range(s-1, -1, -1):
346 eAw = eAw @ eAw
348 # diagonal
349 np.einsum('ii->i', eAw)[:] = np.exp(diag_aw * 2.**(-i))
350 exp_sd = _exp_sinch(diag_aw * (2.**(-i))) * (sd * 2**(-i))
351 if lu[1] == 0: # lower
352 np.einsum('ii->i', eAw[1:, :-1])[:] = exp_sd
353 else: # upper
354 np.einsum('ii->i', eAw[:-1, 1:])[:] = exp_sd
356 else: # generic
357 for _ in range(s):
358 eAw = eAw @ eAw
360 # Zero out the entries from np.empty in case of triangular input
361 if (lu[0] == 0) or (lu[1] == 0):
362 eA[ind] = np.triu(eAw) if lu[0] == 0 else np.tril(eAw)
363 else:
364 eA[ind] = eAw
366 return eA
369def _exp_sinch(x):
370 # Higham's formula (10.42), might overflow, see GH-11839
371 lexp_diff = np.diff(np.exp(x))
372 l_diff = np.diff(x)
373 mask_z = l_diff == 0.
374 lexp_diff[~mask_z] /= l_diff[~mask_z]
375 lexp_diff[mask_z] = np.exp(x[:-1][mask_z])
376 return lexp_diff
379def cosm(A):
380 """
381 Compute the matrix cosine.
383 This routine uses expm to compute the matrix exponentials.
385 Parameters
386 ----------
387 A : (N, N) array_like
388 Input array
390 Returns
391 -------
392 cosm : (N, N) ndarray
393 Matrix cosine of A
395 Examples
396 --------
397 >>> import numpy as np
398 >>> from scipy.linalg import expm, sinm, cosm
400 Euler's identity (exp(i*theta) = cos(theta) + i*sin(theta))
401 applied to a matrix:
403 >>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
404 >>> expm(1j*a)
405 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
406 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
407 >>> cosm(a) + 1j*sinm(a)
408 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
409 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
411 """
412 A = _asarray_square(A)
413 if np.iscomplexobj(A):
414 return 0.5*(expm(1j*A) + expm(-1j*A))
415 else:
416 return expm(1j*A).real
419def sinm(A):
420 """
421 Compute the matrix sine.
423 This routine uses expm to compute the matrix exponentials.
425 Parameters
426 ----------
427 A : (N, N) array_like
428 Input array.
430 Returns
431 -------
432 sinm : (N, N) ndarray
433 Matrix sine of `A`
435 Examples
436 --------
437 >>> import numpy as np
438 >>> from scipy.linalg import expm, sinm, cosm
440 Euler's identity (exp(i*theta) = cos(theta) + i*sin(theta))
441 applied to a matrix:
443 >>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
444 >>> expm(1j*a)
445 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
446 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
447 >>> cosm(a) + 1j*sinm(a)
448 array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
449 [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
451 """
452 A = _asarray_square(A)
453 if np.iscomplexobj(A):
454 return -0.5j*(expm(1j*A) - expm(-1j*A))
455 else:
456 return expm(1j*A).imag
459def tanm(A):
460 """
461 Compute the matrix tangent.
463 This routine uses expm to compute the matrix exponentials.
465 Parameters
466 ----------
467 A : (N, N) array_like
468 Input array.
470 Returns
471 -------
472 tanm : (N, N) ndarray
473 Matrix tangent of `A`
475 Examples
476 --------
477 >>> import numpy as np
478 >>> from scipy.linalg import tanm, sinm, cosm
479 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
480 >>> t = tanm(a)
481 >>> t
482 array([[ -2.00876993, -8.41880636],
483 [ -2.80626879, -10.42757629]])
485 Verify tanm(a) = sinm(a).dot(inv(cosm(a)))
487 >>> s = sinm(a)
488 >>> c = cosm(a)
489 >>> s.dot(np.linalg.inv(c))
490 array([[ -2.00876993, -8.41880636],
491 [ -2.80626879, -10.42757629]])
493 """
494 A = _asarray_square(A)
495 return _maybe_real(A, solve(cosm(A), sinm(A)))
498def coshm(A):
499 """
500 Compute the hyperbolic matrix cosine.
502 This routine uses expm to compute the matrix exponentials.
504 Parameters
505 ----------
506 A : (N, N) array_like
507 Input array.
509 Returns
510 -------
511 coshm : (N, N) ndarray
512 Hyperbolic matrix cosine of `A`
514 Examples
515 --------
516 >>> import numpy as np
517 >>> from scipy.linalg import tanhm, sinhm, coshm
518 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
519 >>> c = coshm(a)
520 >>> c
521 array([[ 11.24592233, 38.76236492],
522 [ 12.92078831, 50.00828725]])
524 Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))
526 >>> t = tanhm(a)
527 >>> s = sinhm(a)
528 >>> t - s.dot(np.linalg.inv(c))
529 array([[ 2.72004641e-15, 4.55191440e-15],
530 [ 0.00000000e+00, -5.55111512e-16]])
532 """
533 A = _asarray_square(A)
534 return _maybe_real(A, 0.5 * (expm(A) + expm(-A)))
537def sinhm(A):
538 """
539 Compute the hyperbolic matrix sine.
541 This routine uses expm to compute the matrix exponentials.
543 Parameters
544 ----------
545 A : (N, N) array_like
546 Input array.
548 Returns
549 -------
550 sinhm : (N, N) ndarray
551 Hyperbolic matrix sine of `A`
553 Examples
554 --------
555 >>> import numpy as np
556 >>> from scipy.linalg import tanhm, sinhm, coshm
557 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
558 >>> s = sinhm(a)
559 >>> s
560 array([[ 10.57300653, 39.28826594],
561 [ 13.09608865, 49.86127247]])
563 Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))
565 >>> t = tanhm(a)
566 >>> c = coshm(a)
567 >>> t - s.dot(np.linalg.inv(c))
568 array([[ 2.72004641e-15, 4.55191440e-15],
569 [ 0.00000000e+00, -5.55111512e-16]])
571 """
572 A = _asarray_square(A)
573 return _maybe_real(A, 0.5 * (expm(A) - expm(-A)))
576def tanhm(A):
577 """
578 Compute the hyperbolic matrix tangent.
580 This routine uses expm to compute the matrix exponentials.
582 Parameters
583 ----------
584 A : (N, N) array_like
585 Input array
587 Returns
588 -------
589 tanhm : (N, N) ndarray
590 Hyperbolic matrix tangent of `A`
592 Examples
593 --------
594 >>> import numpy as np
595 >>> from scipy.linalg import tanhm, sinhm, coshm
596 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
597 >>> t = tanhm(a)
598 >>> t
599 array([[ 0.3428582 , 0.51987926],
600 [ 0.17329309, 0.86273746]])
602 Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))
604 >>> s = sinhm(a)
605 >>> c = coshm(a)
606 >>> t - s.dot(np.linalg.inv(c))
607 array([[ 2.72004641e-15, 4.55191440e-15],
608 [ 0.00000000e+00, -5.55111512e-16]])
610 """
611 A = _asarray_square(A)
612 return _maybe_real(A, solve(coshm(A), sinhm(A)))
615def funm(A, func, disp=True):
616 """
617 Evaluate a matrix function specified by a callable.
619 Returns the value of matrix-valued function ``f`` at `A`. The
620 function ``f`` is an extension of the scalar-valued function `func`
621 to matrices.
623 Parameters
624 ----------
625 A : (N, N) array_like
626 Matrix at which to evaluate the function
627 func : callable
628 Callable object that evaluates a scalar function f.
629 Must be vectorized (eg. using vectorize).
630 disp : bool, optional
631 Print warning if error in the result is estimated large
632 instead of returning estimated error. (Default: True)
634 Returns
635 -------
636 funm : (N, N) ndarray
637 Value of the matrix function specified by func evaluated at `A`
638 errest : float
639 (if disp == False)
641 1-norm of the estimated error, ||err||_1 / ||A||_1
643 Notes
644 -----
645 This function implements the general algorithm based on Schur decomposition
646 (Algorithm 9.1.1. in [1]_).
648 If the input matrix is known to be diagonalizable, then relying on the
649 eigendecomposition is likely to be faster. For example, if your matrix is
650 Hermitian, you can do
652 >>> from scipy.linalg import eigh
653 >>> def funm_herm(a, func, check_finite=False):
654 ... w, v = eigh(a, check_finite=check_finite)
655 ... ## if you further know that your matrix is positive semidefinite,
656 ... ## you can optionally guard against precision errors by doing
657 ... # w = np.maximum(w, 0)
658 ... w = func(w)
659 ... return (v * w).dot(v.conj().T)
661 References
662 ----------
663 .. [1] Gene H. Golub, Charles F. van Loan, Matrix Computations 4th ed.
665 Examples
666 --------
667 >>> import numpy as np
668 >>> from scipy.linalg import funm
669 >>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
670 >>> funm(a, lambda x: x*x)
671 array([[ 4., 15.],
672 [ 5., 19.]])
673 >>> a.dot(a)
674 array([[ 4., 15.],
675 [ 5., 19.]])
677 """
678 A = _asarray_square(A)
679 # Perform Shur decomposition (lapack ?gees)
680 T, Z = schur(A)
681 T, Z = rsf2csf(T, Z)
682 n, n = T.shape
683 F = diag(func(diag(T))) # apply function to diagonal elements
684 F = F.astype(T.dtype.char) # e.g., when F is real but T is complex
686 minden = abs(T[0, 0])
688 # implement Algorithm 11.1.1 from Golub and Van Loan
689 # "matrix Computations."
690 for p in range(1, n):
691 for i in range(1, n-p+1):
692 j = i + p
693 s = T[i-1, j-1] * (F[j-1, j-1] - F[i-1, i-1])
694 ksl = slice(i, j-1)
695 val = dot(T[i-1, ksl], F[ksl, j-1]) - dot(F[i-1, ksl], T[ksl, j-1])
696 s = s + val
697 den = T[j-1, j-1] - T[i-1, i-1]
698 if den != 0.0:
699 s = s / den
700 F[i-1, j-1] = s
701 minden = min(minden, abs(den))
703 F = dot(dot(Z, F), transpose(conjugate(Z)))
704 F = _maybe_real(A, F)
706 tol = {0: feps, 1: eps}[_array_precision[F.dtype.char]]
707 if minden == 0.0:
708 minden = tol
709 err = min(1, max(tol, (tol/minden)*norm(triu(T, 1), 1)))
710 if prod(ravel(logical_not(isfinite(F))), axis=0):
711 err = np.inf
712 if disp:
713 if err > 1000*tol:
714 print("funm result may be inaccurate, approximate err =", err)
715 return F
716 else:
717 return F, err
720def signm(A, disp=True):
721 """
722 Matrix sign function.
724 Extension of the scalar sign(x) to matrices.
726 Parameters
727 ----------
728 A : (N, N) array_like
729 Matrix at which to evaluate the sign function
730 disp : bool, optional
731 Print warning if error in the result is estimated large
732 instead of returning estimated error. (Default: True)
734 Returns
735 -------
736 signm : (N, N) ndarray
737 Value of the sign function at `A`
738 errest : float
739 (if disp == False)
741 1-norm of the estimated error, ||err||_1 / ||A||_1
743 Examples
744 --------
745 >>> from scipy.linalg import signm, eigvals
746 >>> a = [[1,2,3], [1,2,1], [1,1,1]]
747 >>> eigvals(a)
748 array([ 4.12488542+0.j, -0.76155718+0.j, 0.63667176+0.j])
749 >>> eigvals(signm(a))
750 array([-1.+0.j, 1.+0.j, 1.+0.j])
752 """
753 A = _asarray_square(A)
755 def rounded_sign(x):
756 rx = np.real(x)
757 if rx.dtype.char == 'f':
758 c = 1e3*feps*amax(x)
759 else:
760 c = 1e3*eps*amax(x)
761 return sign((absolute(rx) > c) * rx)
762 result, errest = funm(A, rounded_sign, disp=0)
763 errtol = {0: 1e3*feps, 1: 1e3*eps}[_array_precision[result.dtype.char]]
764 if errest < errtol:
765 return result
767 # Handle signm of defective matrices:
769 # See "E.D.Denman and J.Leyva-Ramos, Appl.Math.Comp.,
770 # 8:237-250,1981" for how to improve the following (currently a
771 # rather naive) iteration process:
773 # a = result # sometimes iteration converges faster but where??
775 # Shifting to avoid zero eigenvalues. How to ensure that shifting does
776 # not change the spectrum too much?
777 vals = svd(A, compute_uv=False)
778 max_sv = np.amax(vals)
779 # min_nonzero_sv = vals[(vals>max_sv*errtol).tolist().count(1)-1]
780 # c = 0.5/min_nonzero_sv
781 c = 0.5/max_sv
782 S0 = A + c*np.identity(A.shape[0])
783 prev_errest = errest
784 for i in range(100):
785 iS0 = inv(S0)
786 S0 = 0.5*(S0 + iS0)
787 Pp = 0.5*(dot(S0, S0)+S0)
788 errest = norm(dot(Pp, Pp)-Pp, 1)
789 if errest < errtol or prev_errest == errest:
790 break
791 prev_errest = errest
792 if disp:
793 if not isfinite(errest) or errest >= errtol:
794 print("signm result may be inaccurate, approximate err =", errest)
795 return S0
796 else:
797 return S0, errest
800def khatri_rao(a, b):
801 r"""
802 Khatri-rao product
804 A column-wise Kronecker product of two matrices
806 Parameters
807 ----------
808 a : (n, k) array_like
809 Input array
810 b : (m, k) array_like
811 Input array
813 Returns
814 -------
815 c: (n*m, k) ndarray
816 Khatri-rao product of `a` and `b`.
818 See Also
819 --------
820 kron : Kronecker product
822 Notes
823 -----
824 The mathematical definition of the Khatri-Rao product is:
826 .. math::
828 (A_{ij} \bigotimes B_{ij})_{ij}
830 which is the Kronecker product of every column of A and B, e.g.::
832 c = np.vstack([np.kron(a[:, k], b[:, k]) for k in range(b.shape[1])]).T
834 Examples
835 --------
836 >>> import numpy as np
837 >>> from scipy import linalg
838 >>> a = np.array([[1, 2, 3], [4, 5, 6]])
839 >>> b = np.array([[3, 4, 5], [6, 7, 8], [2, 3, 9]])
840 >>> linalg.khatri_rao(a, b)
841 array([[ 3, 8, 15],
842 [ 6, 14, 24],
843 [ 2, 6, 27],
844 [12, 20, 30],
845 [24, 35, 48],
846 [ 8, 15, 54]])
848 """
849 a = np.asarray(a)
850 b = np.asarray(b)
852 if not (a.ndim == 2 and b.ndim == 2):
853 raise ValueError("The both arrays should be 2-dimensional.")
855 if not a.shape[1] == b.shape[1]:
856 raise ValueError("The number of columns for both arrays "
857 "should be equal.")
859 # c = np.vstack([np.kron(a[:, k], b[:, k]) for k in range(b.shape[1])]).T
860 c = a[..., :, np.newaxis, :] * b[..., np.newaxis, :, :]
861 return c.reshape((-1,) + c.shape[2:])