sparse.linalg._expm_multiply

Compute the action of the matrix exponential.

Module Contents

Classes

LazyOperatorNormInfo(self,A,A_1_norm=None,ell=2,scale=1) Information about an operator is lazily computed.

Functions

_exact_inf_norm(A)
_exact_1_norm(A)
_trace(A)
_ident_like(A)
expm_multiply(A,B,start=None,stop=None,num=None,endpoint=None) Compute the action of the matrix exponential of A on B.
_expm_multiply_simple(A,B,t=1.0,balance=False) Compute the action of the matrix exponential at a single time point.
_expm_multiply_simple_core(A,B,t,mu,m_star,s,tol=None,balance=False) A helper function.
_onenormest_matrix_power(A,p,t=2,itmax=5,compute_v=False,compute_w=False) Efficiently estimate the 1-norm of A^p.
_compute_cost_div_m(m,p,norm_info) A helper function for computing bounds.
_compute_p_max(m_max) Compute the largest positive integer p such that p*(p-1) <= m_max + 1.
_fragment_3_1(norm_info,n0,tol,m_max=55,ell=2) A helper function for the _expm_multiply_* functions.
_condition_3_13(A_1_norm,n0,m_max,ell) A helper function for the _expm_multiply_* functions.
_expm_multiply_interval(A,B,start=None,stop=None,num=None,endpoint=None,balance=False,status_only=False) Compute the action of the matrix exponential at multiple time points.
_expm_multiply_interval_core_0(A,X,h,mu,q,norm_info,tol,ell,n0) A helper function, for the case q <= s.
_expm_multiply_interval_core_1(A,X,h,mu,m_star,s,q,tol) A helper function, for the case q > s and q % s == 0.
_expm_multiply_interval_core_2(A,X,h,mu,m_star,s,q,tol) A helper function, for the case q > s and q % s > 0.
_exact_inf_norm(A)
_exact_1_norm(A)
_trace(A)
_ident_like(A)
expm_multiply(A, B, start=None, stop=None, num=None, endpoint=None)

Compute the action of the matrix exponential of A on B.

A : transposable linear operator
The operator whose exponential is of interest.
B : ndarray
The matrix or vector to be multiplied by the matrix exponential of A.
start : scalar, optional
The starting time point of the sequence.
stop : scalar, optional
The end time point of the sequence, unless endpoint is set to False. In that case, the sequence consists of all but the last of num + 1 evenly spaced time points, so that stop is excluded. Note that the step size changes when endpoint is False.
num : int, optional
Number of time points to use.
endpoint : bool, optional
If True, stop is the last time point. Otherwise, it is not included.
expm_A_B : ndarray
The result of the action .

The optional arguments defining the sequence of evenly spaced time points are compatible with the arguments of numpy.linspace.

The output ndarray shape is somewhat complicated so I explain it here. The ndim of the output could be either 1, 2, or 3. It would be 1 if you are computing the expm action on a single vector at a single time point. It would be 2 if you are computing the expm action on a vector at multiple time points, or if you are computing the expm action on a matrix at a single time point. It would be 3 if you want the action on a matrix with multiple columns at multiple time points. If multiple time points are requested, expm_A_B[0] will always be the action of the expm at the first time point, regardless of whether the action is on a vector or a matrix.

[1]Awad H. Al-Mohy and Nicholas J. Higham (2011) “Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators.” SIAM Journal on Scientific Computing, 33 (2). pp. 488-511. ISSN 1064-8275 http://eprints.ma.man.ac.uk/1591/
[2]Nicholas J. Higham and Awad H. Al-Mohy (2010) “Computing Matrix Functions.” Acta Numerica, 19. 159-208. ISSN 0962-4929 http://eprints.ma.man.ac.uk/1451/
>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import expm, expm_multiply
>>> A = csc_matrix([[1, 0], [0, 1]])
>>> A.todense()
matrix([[1, 0],
        [0, 1]], dtype=int64)
>>> B = np.array([np.exp(-1.), np.exp(-2.)])
>>> B
array([ 0.36787944,  0.13533528])
>>> expm_multiply(A, B, start=1, stop=2, num=3, endpoint=True)
array([[ 1.        ,  0.36787944],
       [ 1.64872127,  0.60653066],
       [ 2.71828183,  1.        ]])
>>> expm(A).dot(B)                  # Verify 1st timestep
array([ 1.        ,  0.36787944])
>>> expm(1.5*A).dot(B)              # Verify 2nd timestep
array([ 1.64872127,  0.60653066])
>>> expm(2*A).dot(B)                # Verify 3rd timestep
array([ 2.71828183,  1.        ])
_expm_multiply_simple(A, B, t=1.0, balance=False)

Compute the action of the matrix exponential at a single time point.

A : transposable linear operator
The operator whose exponential is of interest.
B : ndarray
The matrix to be multiplied by the matrix exponential of A.
t : float
A time point.
balance : bool
Indicates whether or not to apply balancing.
F : ndarray

This is algorithm (3.2) in Al-Mohy and Higham (2011).

_expm_multiply_simple_core(A, B, t, mu, m_star, s, tol=None, balance=False)

A helper function.

_onenormest_matrix_power(A, p, t=2, itmax=5, compute_v=False, compute_w=False)

Efficiently estimate the 1-norm of A^p.

A : ndarray
Matrix whose 1-norm of a power is to be computed.
p : int
Non-negative integer power.
t : int, optional
A positive parameter controlling the tradeoff between accuracy versus time and memory usage. Larger values take longer and use more memory but give more accurate output.
itmax : int, optional
Use at most this many iterations.
compute_v : bool, optional
Request a norm-maximizing linear operator input vector if True.
compute_w : bool, optional
Request a norm-maximizing linear operator output vector if True.
est : float
An underestimate of the 1-norm of the sparse matrix.
v : ndarray, optional
The vector such that ||Av||_1 == est*||v||_1. It can be thought of as an input to the linear operator that gives an output with particularly large norm.
w : ndarray, optional
The vector Av which has relatively large 1-norm. It can be thought of as an output of the linear operator that is relatively large in norm compared to the input.
class LazyOperatorNormInfo(A, A_1_norm=None, ell=2, scale=1)

Information about an operator is lazily computed.

The information includes the exact 1-norm of the operator, in addition to estimates of 1-norms of powers of the operator. This uses the notation of Computing the Action (2011). This class is specialized enough to probably not be of general interest outside of this module.

__init__(A, A_1_norm=None, ell=2, scale=1)

Provide the operator and some norm-related information.

A : linear operator
The operator of interest.
A_1_norm : float, optional
The exact 1-norm of A.
ell : int, optional
A technical parameter controlling norm estimation quality.
scale : int, optional
If specified, return the norms of scale*A instead of A.
set_scale(scale)

Set the scale parameter.

onenorm()

Compute the exact 1-norm.

d(p)

Lazily estimate d_p(A) ~= || A^p ||^(1/p) where ||.|| is the 1-norm.

alpha(p)

Lazily compute max(d(p), d(p+1)).

_compute_cost_div_m(m, p, norm_info)

A helper function for computing bounds.

This is equation (3.10). It measures cost in terms of the number of required matrix products.

m : int
A valid key of _theta.
p : int
A matrix power.
norm_info : LazyOperatorNormInfo
Information about 1-norms of related operators.
cost_div_m : int
Required number of matrix products divided by m.
_compute_p_max(m_max)

Compute the largest positive integer p such that p*(p-1) <= m_max + 1.

Do this in a slightly dumb way, but safe and not too slow.

m_max : int
A count related to bounds.
_fragment_3_1(norm_info, n0, tol, m_max=55, ell=2)

A helper function for the _expm_multiply_* functions.

norm_info : LazyOperatorNormInfo
Information about norms of certain linear operators of interest.
n0 : int
Number of columns in the _expm_multiply_* B matrix.
tol : float
Expected to be for single precision or for double precision.
m_max : int
A value related to a bound.
ell : int
The number of columns used in the 1-norm approximation. This is usually taken to be small, maybe between 1 and 5.
best_m : int
Related to bounds for error control.
best_s : int
Amount of scaling.

This is code fragment (3.1) in Al-Mohy and Higham (2011). The discussion of default values for m_max and ell is given between the definitions of equation (3.11) and the definition of equation (3.12).

_condition_3_13(A_1_norm, n0, m_max, ell)

A helper function for the _expm_multiply_* functions.

A_1_norm : float
The precomputed 1-norm of A.
n0 : int
Number of columns in the _expm_multiply_* B matrix.
m_max : int
A value related to a bound.
ell : int
The number of columns used in the 1-norm approximation. This is usually taken to be small, maybe between 1 and 5.
value : bool
Indicates whether or not the condition has been met.

This is condition (3.13) in Al-Mohy and Higham (2011).

_expm_multiply_interval(A, B, start=None, stop=None, num=None, endpoint=None, balance=False, status_only=False)

Compute the action of the matrix exponential at multiple time points.

A : transposable linear operator
The operator whose exponential is of interest.
B : ndarray
The matrix to be multiplied by the matrix exponential of A.
start : scalar, optional
The starting time point of the sequence.
stop : scalar, optional
The end time point of the sequence, unless endpoint is set to False. In that case, the sequence consists of all but the last of num + 1 evenly spaced time points, so that stop is excluded. Note that the step size changes when endpoint is False.
num : int, optional
Number of time points to use.
endpoint : bool, optional
If True, stop is the last time point. Otherwise, it is not included.
balance : bool
Indicates whether or not to apply balancing.
status_only : bool
A flag that is set to True for some debugging and testing operations.
F : ndarray
status : int
An integer status for testing and debugging.

This is algorithm (5.2) in Al-Mohy and Higham (2011).

There seems to be a typo, where line 15 of the algorithm should be moved to line 6.5 (between lines 6 and 7).

_expm_multiply_interval_core_0(A, X, h, mu, q, norm_info, tol, ell, n0)

A helper function, for the case q <= s.

_expm_multiply_interval_core_1(A, X, h, mu, m_star, s, q, tol)

A helper function, for the case q > s and q % s == 0.

_expm_multiply_interval_core_2(A, X, h, mu, m_star, s, q, tol)

A helper function, for the case q > s and q % s > 0.