# # Author: Travis Oliphant, March 2002 # from __future__ import division, print_function, absolute_import __all__ = ['expm','expm2','expm3','cosm','sinm','tanm','coshm','sinhm', 'tanhm','logm','funm','signm','sqrtm', 'expm_frechet', 'expm_cond', 'fractional_matrix_power'] from numpy import (Inf, dot, diag, exp, product, logical_not, cast, ravel, transpose, conjugate, absolute, amax, sign, isfinite, sqrt, single) import numpy as np import warnings # Local imports from .misc import norm from .basic import solve, inv from .special_matrices import triu from .decomp import eig from .decomp_svd import svd from .decomp_schur import schur, rsf2csf from ._expm_frechet import expm_frechet, expm_cond from ._matfuncs_sqrtm import sqrtm eps = np.finfo(float).eps feps = np.finfo(single).eps _array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1} ############################################################################### # Utility functions. def _asarray_square(A): """ Wraps asarray with the extra requirement that the input be a square matrix. The motivation is that the matfuncs module has real functions that have been lifted to square matrix functions. Parameters ---------- A : array_like A square matrix. Returns ------- out : ndarray An ndarray copy or view or other representation of A. """ A = np.asarray(A) if len(A.shape) != 2 or A.shape[0] != A.shape[1]: raise ValueError('expected square array_like input') return A def _maybe_real(A, B, tol=None): """ Return either B or the real part of B, depending on properties of A and B. The motivation is that B has been computed as a complicated function of A, and B may be perturbed by negligible imaginary components. If A is real and B is complex with small imaginary components, then return a real copy of B. The assumption in that case would be that the imaginary components of B are numerical artifacts. Parameters ---------- A : ndarray Input array whose type is to be checked as real vs. complex. B : ndarray Array to be returned, possibly without its imaginary part. tol : float Absolute tolerance. Returns ------- out : real or complex array Either the input array B or only the real part of the input array B. """ # Note that booleans and integers compare as real. if np.isrealobj(A) and np.iscomplexobj(B): if tol is None: tol = {0:feps*1e3, 1:eps*1e6}[_array_precision[B.dtype.char]] if np.allclose(B.imag, 0.0, atol=tol): B = B.real return B ############################################################################### # Matrix functions. def fractional_matrix_power(A, t): # This fixes some issue with imports; # this function calls onenormest which is in scipy.sparse. A = _asarray_square(A) import scipy.linalg._matfuncs_inv_ssq return scipy.linalg._matfuncs_inv_ssq.fractional_matrix_power(A, t) def logm(A, disp=True): """ Compute matrix logarithm. The matrix logarithm is the inverse of expm: expm(logm(`A`)) == `A` Parameters ---------- A : (N, N) array_like Matrix whose logarithm to evaluate disp : bool, optional Print warning if error in the result is estimated large instead of returning estimated error. (Default: True) Returns ------- logm : (N, N) ndarray Matrix logarithm of `A` errest : float (if disp == False) 1-norm of the estimated error, ||err||_1 / ||A||_1 """ A = _asarray_square(A) # Avoid circular import ... this is OK, right? import scipy.linalg._matfuncs_inv_ssq F = scipy.linalg._matfuncs_inv_ssq.logm(A) errtol = 1000*eps #TODO use a better error approximation errest = norm(expm(F)-A,1) / norm(A,1) if disp: if not isfinite(errest) or errest >= errtol: print("logm result may be inaccurate, approximate err =", errest) return F else: return F, errest def expm(A, q=None): """ Compute the matrix exponential using Pade approximation. Parameters ---------- A : (N, N) array_like or sparse matrix Matrix to be exponentiated. Returns ------- expm : (N, N) ndarray Matrix exponential of `A`. References ---------- .. [1] Awad H. Al-Mohy and Nicholas J. Higham (2009) "A New Scaling and Squaring Algorithm for the Matrix Exponential." SIAM Journal on Matrix Analysis and Applications. 31 (3). pp. 970-989. ISSN 1095-7162 """ if q is not None: msg = "argument q=... in scipy.linalg.expm is deprecated." warnings.warn(msg, DeprecationWarning) # Input checking and conversion is provided by sparse.linalg.expm(). import scipy.sparse.linalg return scipy.sparse.linalg.expm(A) # deprecated, but probably should be left there in the long term @np.deprecate(new_name="expm") def expm2(A): """ Compute the matrix exponential using eigenvalue decomposition. Parameters ---------- A : (N, N) array_like Matrix to be exponentiated Returns ------- expm2 : (N, N) ndarray Matrix exponential of `A` """ A = _asarray_square(A) t = A.dtype.char if t not in ['f','F','d','D']: A = A.astype('d') t = 'd' s, vr = eig(A) vri = inv(vr) r = dot(dot(vr, diag(exp(s))), vri) if t in ['f', 'd']: return r.real.astype(t) else: return r.astype(t) # deprecated, but probably should be left there in the long term @np.deprecate(new_name="expm") def expm3(A, q=20): """ Compute the matrix exponential using Taylor series. Parameters ---------- A : (N, N) array_like Matrix to be exponentiated q : int Order of the Taylor series used is `q-1` Returns ------- expm3 : (N, N) ndarray Matrix exponential of `A` """ A = _asarray_square(A) n = A.shape[0] t = A.dtype.char if t not in ['f','F','d','D']: A = A.astype('d') t = 'd' eA = np.identity(n, dtype=t) trm = np.identity(n, dtype=t) castfunc = cast[t] for k in range(1, q): trm[:] = trm.dot(A) / castfunc(k) eA += trm return eA def cosm(A): """ Compute the matrix cosine. This routine uses expm to compute the matrix exponentials. Parameters ---------- A : (N, N) array_like Input array Returns ------- cosm : (N, N) ndarray Matrix cosine of A """ A = _asarray_square(A) if np.iscomplexobj(A): return 0.5*(expm(1j*A) + expm(-1j*A)) else: return expm(1j*A).real def sinm(A): """ Compute the matrix sine. This routine uses expm to compute the matrix exponentials. Parameters ---------- A : (N, N) array_like Input array. Returns ------- sinm : (N, N) ndarray Matrix cosine of `A` """ A = _asarray_square(A) if np.iscomplexobj(A): return -0.5j*(expm(1j*A) - expm(-1j*A)) else: return expm(1j*A).imag def tanm(A): """ Compute the matrix tangent. This routine uses expm to compute the matrix exponentials. Parameters ---------- A : (N, N) array_like Input array. Returns ------- tanm : (N, N) ndarray Matrix tangent of `A` """ A = _asarray_square(A) return _maybe_real(A, solve(cosm(A), sinm(A))) def coshm(A): """ Compute the hyperbolic matrix cosine. This routine uses expm to compute the matrix exponentials. Parameters ---------- A : (N, N) array_like Input array. Returns ------- coshm : (N, N) ndarray Hyperbolic matrix cosine of `A` """ A = _asarray_square(A) return _maybe_real(A, 0.5 * (expm(A) + expm(-A))) def sinhm(A): """ Compute the hyperbolic matrix sine. This routine uses expm to compute the matrix exponentials. Parameters ---------- A : (N, N) array_like Input array. Returns ------- sinhm : (N, N) ndarray Hyperbolic matrix sine of `A` """ A = _asarray_square(A) return _maybe_real(A, 0.5 * (expm(A) - expm(-A))) def tanhm(A): """ Compute the hyperbolic matrix tangent. This routine uses expm to compute the matrix exponentials. Parameters ---------- A : (N, N) array_like Input array Returns ------- tanhm : (N, N) ndarray Hyperbolic matrix tangent of `A` """ A = _asarray_square(A) return _maybe_real(A, solve(coshm(A), sinhm(A))) def funm(A, func, disp=True): """ Evaluate a matrix function specified by a callable. Returns the value of matrix-valued function ``f`` at `A`. The function ``f`` is an extension of the scalar-valued function `func` to matrices. Parameters ---------- A : (N, N) array_like Matrix at which to evaluate the function func : callable Callable object that evaluates a scalar function f. Must be vectorized (eg. using vectorize). disp : bool, optional Print warning if error in the result is estimated large instead of returning estimated error. (Default: True) Returns ------- funm : (N, N) ndarray Value of the matrix function specified by func evaluated at `A` errest : float (if disp == False) 1-norm of the estimated error, ||err||_1 / ||A||_1 """ A = _asarray_square(A) # Perform Shur decomposition (lapack ?gees) T, Z = schur(A) T, Z = rsf2csf(T,Z) n,n = T.shape F = diag(func(diag(T))) # apply function to diagonal elements F = F.astype(T.dtype.char) # e.g. when F is real but T is complex minden = abs(T[0,0]) # implement Algorithm 11.1.1 from Golub and Van Loan # "matrix Computations." for p in range(1,n): for i in range(1,n-p+1): j = i + p s = T[i-1,j-1] * (F[j-1,j-1] - F[i-1,i-1]) ksl = slice(i,j-1) val = dot(T[i-1,ksl],F[ksl,j-1]) - dot(F[i-1,ksl],T[ksl,j-1]) s = s + val den = T[j-1,j-1] - T[i-1,i-1] if den != 0.0: s = s / den F[i-1,j-1] = s minden = min(minden,abs(den)) F = dot(dot(Z, F), transpose(conjugate(Z))) F = _maybe_real(A, F) tol = {0:feps, 1:eps}[_array_precision[F.dtype.char]] if minden == 0.0: minden = tol err = min(1, max(tol,(tol/minden)*norm(triu(T,1),1))) if product(ravel(logical_not(isfinite(F))),axis=0): err = Inf if disp: if err > 1000*tol: print("funm result may be inaccurate, approximate err =", err) return F else: return F, err def signm(A, disp=True): """ Matrix sign function. Extension of the scalar sign(x) to matrices. Parameters ---------- A : (N, N) array_like Matrix at which to evaluate the sign function disp : bool, optional Print warning if error in the result is estimated large instead of returning estimated error. (Default: True) Returns ------- signm : (N, N) ndarray Value of the sign function at `A` errest : float (if disp == False) 1-norm of the estimated error, ||err||_1 / ||A||_1 Examples -------- >>> from scipy.linalg import signm, eigvals >>> a = [[1,2,3], [1,2,1], [1,1,1]] >>> eigvals(a) array([ 4.12488542+0.j, -0.76155718+0.j, 0.63667176+0.j]) >>> eigvals(signm(a)) array([-1.+0.j, 1.+0.j, 1.+0.j]) """ A = _asarray_square(A) def rounded_sign(x): rx = np.real(x) if rx.dtype.char == 'f': c = 1e3*feps*amax(x) else: c = 1e3*eps*amax(x) return sign((absolute(rx) > c) * rx) result, errest = funm(A, rounded_sign, disp=0) errtol = {0:1e3*feps, 1:1e3*eps}[_array_precision[result.dtype.char]] if errest < errtol: return result # Handle signm of defective matrices: # See "E.D.Denman and J.Leyva-Ramos, Appl.Math.Comp., # 8:237-250,1981" for how to improve the following (currently a # rather naive) iteration process: # a = result # sometimes iteration converges faster but where?? # Shifting to avoid zero eigenvalues. How to ensure that shifting does # not change the spectrum too much? vals = svd(A, compute_uv=0) max_sv = np.amax(vals) # min_nonzero_sv = vals[(vals>max_sv*errtol).tolist().count(1)-1] # c = 0.5/min_nonzero_sv c = 0.5/max_sv S0 = A + c*np.identity(A.shape[0]) prev_errest = errest for i in range(100): iS0 = inv(S0) S0 = 0.5*(S0 + iS0) Pp = 0.5*(dot(S0,S0)+S0) errest = norm(dot(Pp,Pp)-Pp,1) if errest < errtol or prev_errest == errest: break prev_errest = errest if disp: if not isfinite(errest) or errest >= errtol: print("signm result may be inaccurate, approximate err =", errest) return S0 else: return S0, errest