linalg.matfuncs

Module Contents

Functions

_asarray_square(A) Wraps asarray with the extra requirement that the input be a square matrix.
_maybe_real(A,B,tol=None) Return either B or the real part of B, depending on properties of A and B.
fractional_matrix_power(A,t) Compute the fractional power of a matrix.
logm(A,disp=True) Compute matrix logarithm.
expm(A) Compute the matrix exponential using Pade approximation.
cosm(A) Compute the matrix cosine.
sinm(A) Compute the matrix sine.
tanm(A) Compute the matrix tangent.
coshm(A) Compute the hyperbolic matrix cosine.
sinhm(A) Compute the hyperbolic matrix sine.
tanhm(A) Compute the hyperbolic matrix tangent.
funm(A,func,disp=True) Evaluate a matrix function specified by a callable.
signm(A,disp=True) Matrix sign function.
_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.

A : array_like
A square matrix.
out : ndarray
An ndarray copy or view or other representation of A.
_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.

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.
out : real or complex array
Either the input array B or only the real part of the input array B.
fractional_matrix_power(A, t)

Compute the fractional power of a matrix.

Proceeds according to the discussion in section (6) of [1]_.

A : (N, N) array_like
Matrix whose fractional power to evaluate.
t : float
Fractional power.
X : (N, N) array_like
The fractional power of the matrix.
[1]Nicholas J. Higham and Lijing lin (2011) “A Schur-Pade Algorithm for Fractional Powers of a Matrix.” SIAM Journal on Matrix Analysis and Applications, 32 (3). pp. 1056-1078. ISSN 0895-4798
>>> from scipy.linalg import fractional_matrix_power
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> b = fractional_matrix_power(a, 0.5)
>>> b
array([[ 0.75592895,  1.13389342],
       [ 0.37796447,  1.88982237]])
>>> np.dot(b, b)      # Verify square root
array([[ 1.,  3.],
       [ 1.,  4.]])
logm(A, disp=True)

Compute matrix logarithm.

The matrix logarithm is the inverse of expm: expm(logm(A)) == A

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)
logm : (N, N) ndarray
Matrix logarithm of A
errest : float

(if disp == False)

1-norm of the estimated error, ||err||_1 / ||A||_1

[1]Awad H. Al-Mohy and Nicholas J. Higham (2012) “Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm.” SIAM Journal on Scientific Computing, 34 (4). C152-C169. ISSN 1095-7197
[2]Nicholas J. Higham (2008) “Functions of Matrices: Theory and Computation” ISBN 978-0-898716-46-7
[3]Nicholas J. Higham and Lijing lin (2011) “A Schur-Pade Algorithm for Fractional Powers of a Matrix.” SIAM Journal on Matrix Analysis and Applications, 32 (3). pp. 1056-1078. ISSN 0895-4798
>>> from scipy.linalg import logm, expm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> b = logm(a)
>>> b
array([[-1.02571087,  2.05142174],
       [ 0.68380725,  1.02571087]])
>>> expm(b)         # Verify expm(logm(a)) returns a
array([[ 1.,  3.],
       [ 1.,  4.]])
expm(A)

Compute the matrix exponential using Pade approximation.

A : (N, N) array_like or sparse matrix
Matrix to be exponentiated.
expm : (N, N) ndarray
Matrix exponential of A.
[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
>>> from scipy.linalg import expm, sinm, cosm

Matrix version of the formula exp(0) = 1:

>>> expm(np.zeros((2,2)))
array([[ 1.,  0.],
       [ 0.,  1.]])

Euler’s identity (exp(i*theta) = cos(theta) + i*sin(theta)) applied to a matrix:

>>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
>>> expm(1j*a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
       [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
>>> cosm(a) + 1j*sinm(a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
       [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
cosm(A)

Compute the matrix cosine.

This routine uses expm to compute the matrix exponentials.

A : (N, N) array_like
Input array
cosm : (N, N) ndarray
Matrix cosine of A
>>> from scipy.linalg import expm, sinm, cosm

Euler’s identity (exp(i*theta) = cos(theta) + i*sin(theta)) applied to a matrix:

>>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
>>> expm(1j*a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
       [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
>>> cosm(a) + 1j*sinm(a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
       [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
sinm(A)

Compute the matrix sine.

This routine uses expm to compute the matrix exponentials.

A : (N, N) array_like
Input array.
sinm : (N, N) ndarray
Matrix sine of A
>>> from scipy.linalg import expm, sinm, cosm

Euler’s identity (exp(i*theta) = cos(theta) + i*sin(theta)) applied to a matrix:

>>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
>>> expm(1j*a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
       [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
>>> cosm(a) + 1j*sinm(a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
       [ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
tanm(A)

Compute the matrix tangent.

This routine uses expm to compute the matrix exponentials.

A : (N, N) array_like
Input array.
tanm : (N, N) ndarray
Matrix tangent of A
>>> from scipy.linalg import tanm, sinm, cosm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> t = tanm(a)
>>> t
array([[ -2.00876993,  -8.41880636],
       [ -2.80626879, -10.42757629]])

Verify tanm(a) = sinm(a).dot(inv(cosm(a)))

>>> s = sinm(a)
>>> c = cosm(a)
>>> s.dot(np.linalg.inv(c))
array([[ -2.00876993,  -8.41880636],
       [ -2.80626879, -10.42757629]])
coshm(A)

Compute the hyperbolic matrix cosine.

This routine uses expm to compute the matrix exponentials.

A : (N, N) array_like
Input array.
coshm : (N, N) ndarray
Hyperbolic matrix cosine of A
>>> from scipy.linalg import tanhm, sinhm, coshm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> c = coshm(a)
>>> c
array([[ 11.24592233,  38.76236492],
       [ 12.92078831,  50.00828725]])

Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))

>>> t = tanhm(a)
>>> s = sinhm(a)
>>> t - s.dot(np.linalg.inv(c))
array([[  2.72004641e-15,   4.55191440e-15],
       [  0.00000000e+00,  -5.55111512e-16]])
sinhm(A)

Compute the hyperbolic matrix sine.

This routine uses expm to compute the matrix exponentials.

A : (N, N) array_like
Input array.
sinhm : (N, N) ndarray
Hyperbolic matrix sine of A
>>> from scipy.linalg import tanhm, sinhm, coshm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> s = sinhm(a)
>>> s
array([[ 10.57300653,  39.28826594],
       [ 13.09608865,  49.86127247]])

Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))

>>> t = tanhm(a)
>>> c = coshm(a)
>>> t - s.dot(np.linalg.inv(c))
array([[  2.72004641e-15,   4.55191440e-15],
       [  0.00000000e+00,  -5.55111512e-16]])
tanhm(A)

Compute the hyperbolic matrix tangent.

This routine uses expm to compute the matrix exponentials.

A : (N, N) array_like
Input array
tanhm : (N, N) ndarray
Hyperbolic matrix tangent of A
>>> from scipy.linalg import tanhm, sinhm, coshm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> t = tanhm(a)
>>> t
array([[ 0.3428582 ,  0.51987926],
       [ 0.17329309,  0.86273746]])

Verify tanhm(a) = sinhm(a).dot(inv(coshm(a)))

>>> s = sinhm(a)
>>> c = coshm(a)
>>> t - s.dot(np.linalg.inv(c))
array([[  2.72004641e-15,   4.55191440e-15],
       [  0.00000000e+00,  -5.55111512e-16]])
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.

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)
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

>>> from scipy.linalg import funm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> funm(a, lambda x: x*x)
array([[  4.,  15.],
       [  5.,  19.]])
>>> a.dot(a)
array([[  4.,  15.],
       [  5.,  19.]])

This function implements the general algorithm based on Schur decomposition (Algorithm 9.1.1. in [1]_).

If the input matrix is known to be diagonalizable, then relying on the eigendecomposition is likely to be faster. For example, if your matrix is Hermitian, you can do

>>> from scipy.linalg import eigh
>>> def funm_herm(a, func, check_finite=False):
...     w, v = eigh(a, check_finite=check_finite)
...     ## if you further know that your matrix is positive semidefinite,
...     ## you can optionally guard against precision errors by doing
...     # w = np.maximum(w, 0)
...     w = func(w)
...     return (v * w).dot(v.conj().T)
[1]Gene H. Golub, Charles F. van Loan, Matrix Computations 4th ed.
signm(A, disp=True)

Matrix sign function.

Extension of the scalar sign(x) to matrices.

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)
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

>>> 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])