Expv: Matrix Exponentials Times Vectors

expv(t, A, b; kwargs) -> exp(tA)b

Compute the matrix-exponential-vector product with a Krylov approximation.

Arguments

• t: scalar time or scale factor.
• A: square matrix or matrix-free operator satisfying the Matrix-Free Operator Interface page in the manual.
• b: input vector.
• Ks: a precomputed KrylovSubspace for the (A, b) pair.

Keywords

• mode: :happy_breakdown (default) builds a regular Arnoldi basis; :error_estimate uses the Hermitian error-estimate method.
• cache: optional ExpvCache for the reduced exponential.
• expmethod: matrix-exponential implementation for the reduced problem.
• Remaining keywords for the (t, A, b) form are forwarded to arnoldi.

Returns

The vector approximating $\exp(t A)b$.

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
expv(0.1, A, b; m = 2)

A Krylov subspace is constructed using arnoldi and exp! is called on the Hessenberg matrix. Consult arnoldi for the values of the keyword arguments. An alternative algorithm, where an error estimate generated on-the-fly is used to terminate the Krylov iteration, can be employed by setting the kwarg mode=:error_estimate.

expv(t, Ks; cache) -> exp(tA)b

Compute the expv product using a pre-constructed Krylov subspace.

phiv(t,A,b,k;correct,kwargs) -> [phi_0(tA)b phi_1(tA)b ... phi_k(tA)b][, errest]

Compute matrix-phi-vector products with a Krylov approximation. k >= 1.

Arguments

• t: scalar time or scale factor.
• A: square matrix or matrix-free operator satisfying the Matrix-Free Operator Interface page in the manual.
• b: input vector.
• k: highest phi-function order.
• Ks: precomputed KrylovSubspace for the (A, b) pair.

Keywords

• cache: optional PhivCache for reduced phi-function work.
• correct: apply the last-Arnoldi-vector correction to orders 0:k-1.
• errest: return (w, estimate) instead of just w.
• Remaining keywords for the (t, A, b, k) form are forwarded to arnoldi.

Returns

An n by k + 1 matrix whose columns are $\varphi_j(tA)b$ for j = 0:k, or that matrix paired with an error estimate when errest=true.

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
phiv(0.1, A, b, 2; m = 2)

The phi functions are defined as

\[\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - \varphi_k(0)}{z}\]

A Krylov subspace is constructed using arnoldi and phiv_dense is called on the Hessenberg matrix. If correct=true, then phi0 through phik-1 are updated using the last Arnoldi vector v_m+1 ^[1]. If errest=true then an additional error estimate for the second-to-last phi is also returned. For the additional keyword arguments, consult arnoldi.

phiv(t,Ks,k;correct,kwargs) -> [phi0(tA)b phi1(tA)b ... phi_k(tA)b][, errest]

Compute the matrix-phi-vector products using a pre-constructed Krylov subspace.

expv!(w,t,Ks[;cache]) -> w

Compute $\exp(t A)b$ from a precomputed Krylov basis without allocating the output vector.

Arguments

• w: output vector, overwritten in place.
• t: scalar time or scale factor.
• Ks: populated KrylovSubspace.

Keywords

• cache: nothing or a reusable ExpvCache.
• expmethod: reduced matrix-exponential implementation.

Returns

The mutated w.

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
Ks = arnoldi(A, b; m = 2)
w = similar(b)
expv!(w, 0.1, Ks)

expv!(w, t, A, b, Ks, cache)

Calculate the action of exp(t*A) on b, storing the result in w. The Krylov iteration terminates when its subspace error estimate is below the requested tolerance. This method currently supports Hermitian operators and a cache returned by get_subspace_cache.

Arguments

• w: output vector, overwritten in place.
• t: scalar time or scale factor.
• A: Hermitian matrix or matrix-free operator.
• b: input vector.
• Ks: reusable KrylovSubspace with real recurrence coefficients.
• cache: subspace-exponential workspace from get_subspace_cache.

Keywords

• atol = 1e-8: absolute stopping tolerance.
• rtol = 1e-4: relative stopping tolerance, scaled by norm(b).
• m: maximum Lanczos dimension, defaulting to min(Ks.maxiter, size(A, 1)).
• ishermitian = LinearAlgebra.ishermitian(A): must be true; non-Hermitian error estimation is not implemented.
• verbose = false: print each error estimate when true.
• expmethod: reduced matrix-exponential method. Reserved for compatibility.

Returns

The mutated w.

Examples

using LinearAlgebra

A = Hermitian([-2.0 1.0; 1.0 -2.0])
b = ComplexF64[1, 0]
Ks = KrylovSubspace{ComplexF64, Float64}(length(b), 2)
cache = get_subspace_cache(Ks)
w = similar(b)
expv!(w, 0.1, A, b, Ks, cache; atol = 1.0e-10, rtol = 1.0e-8)

phiv!(w,t,Ks,k[;cache,correct,errest]) -> w[,errest]

Compute phi-vector products from a precomputed Krylov basis without allocating the output matrix.

Arguments

• w: output matrix with size(w, 2) == k + 1, overwritten in place.
• t: scalar time or scale factor.
• Ks: populated KrylovSubspace.
• k: highest phi-function order.

Keywords

cache, correct, and errest have the same meaning as for phiv.

Returns

The mutated w, or (w, estimate) when errest=true.

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
Ks = arnoldi(A, b; m = 2)
w = similar(b, length(b), 3)
phiv!(w, 0.1, Ks, 2)

expv_timestep(ts,A,b[;adaptive,tol,kwargs...]) -> U

Evaluates the matrix exponentiation-vector product using time stepping

\[u = \exp(tA)b\]

ts is an array of time snapshots for u, with U[:,j] ≈ u(ts[j]). ts can also be just one value, in which case only the end result is returned and U is a vector.

The time stepping formula of Niesen & Wright is used ^[1]. If the time step tau is not specified, it is chosen according to (17) of Niesen & Wright. If adaptive==true, the time step and Krylov subspace size adaptation scheme of Niesen & Wright is used, the relative tolerance of which can be set using the keyword parameter tol. The delta and gamma parameters of the adaptation scheme can also be adjusted.

Set verbose=true to print out the internal steps (for debugging). For the other keyword arguments, consult arnoldi and phiv, which are used internally.

Note that this function is just a special case of phiv_timestep with a more intuitive interface (vector b instead of a n-by-1 matrix B).

Arguments

• ts or t: requested output time(s).
• A: square matrix or matrix-free operator satisfying the Matrix-Free Operator Interface page in the manual.
• b: initial vector.

Keywords

adaptive, tol, tau, m, iop, opnorm, ishermitian, delta, gamma, NA, and verbose are forwarded to phiv_timestep.

Returns

A vector for scalar t, or a matrix whose column j corresponds to ts[j].

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
expv_timestep(0.1, A, b; tol = 1.0e-8)

phiv_timestep(ts,A,B[;adaptive,tol,kwargs...]) -> U

Evaluates the linear combination of phi-vector products using time stepping

\[u = \varphi_0(tA)b_0 + t\varphi_1(tA)b_1 + \cdots + t^p\varphi_p(tA)b_p\]

ts is an array of time snapshots for u, with U[:,j] ≈ u(ts[j]). ts can also be just one value, in which case only the end result is returned and U is a vector.

The time stepping formula of Niesen & Wright is used ^[1]. If the time step tau is not specified, it is chosen according to (17) of Niesen & Wright. If adaptive==true, the time step and Krylov subspace size adaptation scheme of Niesen & Wright is used, the relative tolerance of which can be set using the keyword parameter tol. The delta and gamma parameters of the adaptation scheme can also be adjusted.

The adaptation scheme drives the internal error estimate below tol times a scalar operator-norm scale of A. By default (opnorm === nothing) that scale is estimated matrix-free from the Arnoldi Hessenberg built on the first Krylov step, so opnorm(A, Inf) is never evaluated: this needs no operator method beyond mul! – suitable for matrix types (e.g. sparse GPU arrays) that do not support opnorm – and reflects A's action on the Krylov subspace the exponential actually explores. The first Arnoldi step supplies the scale used to seed the initial tau. Pass opnorm as a precomputed scalar (a norm or bound) or a function opnorm(A, Inf) to override the estimate with an explicit operator-norm scale, which additionally seeds the initial tau.

When encountering a happy breakdown in the Krylov subspace construction, the time step is set to the remainder of the time interval since time stepping is no longer necessary.

Set verbose=true to print out the internal steps (for debugging). For the other keyword arguments, consult arnoldi and phiv, which are used internally.

Arguments

• ts or t: requested output time(s).
• A: square matrix or matrix-free operator satisfying the Matrix-Free Operator Interface page in the manual.
• B: matrix whose columns are the coefficients b_0, ..., b_p.

Keywords

• tau: initial internal step size. 0 selects an estimate.
• m: maximum Krylov dimension per internal step.
• tol: requested relative tolerance for adaptive stepping.
• opnorm: nothing, a scalar bound, or callable (A, Inf) -> bound.
• adaptive: enable joint step-size and Krylov-dimension adaptation.
• iop, correct, caches, delta, ishermitian, gamma, NA, and verbose: advanced Arnoldi, correction, workspace, and adaptation controls.

Returns

A vector for scalar t, or a matrix whose column j is the requested linear combination at ts[j].

Examples

A = [-2.0 1.0; 0.0 -1.0]
B = [1.0 0.5; 0.0 0.0]
phiv_timestep(0.1, A, B; tol = 1.0e-8)

expv_timestep!(u,t,A,b[;kwargs]) -> u

Non-allocating version of expv_timestep.

Arguments

• u or U: output vector or output matrix, overwritten in place.
• t or ts: requested output time(s).
• A: square matrix or matrix-free operator.
• b: initial vector.

Keywords

• caches: optional reusable workspace tuple from the internal timestep-cache constructor.
• Remaining keywords have the same meaning as for expv_timestep.

Returns

The mutated u or U.

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
u = similar(b)
expv_timestep!(u, 0.1, A, b; tol = 1.0e-8)

phiv_timestep!(U,ts,A,B[;kwargs]) -> U

Non-allocating version of phiv_timestep.

Arguments

• u or U: output vector or matrix, overwritten in place.
• t or ts: requested output time(s).
• A: square matrix or matrix-free operator.
• B: coefficient matrix for the phi-function combination.

Keywords

The keywords are the same as for phiv_timestep.

Returns

The mutated u or U.

Examples

A = [-2.0 1.0; 0.0 -1.0]
B = [1.0 0.5; 0.0 0.0]
u = similar(B, size(B, 1))
phiv_timestep!(u, 0.1, A, B; tol = 1.0e-8)

kiops(tau_out, A, u; kwargs...) -> (w, stats)

Evaluate a linear combination of phi functions at the requested output times using the adaptive KIOPS Krylov algorithm:

\[ w_i = φ_0(t_i A) u_0 + φ_1(t_i A) u_1 + φ_2(t_i A) u_2 + \cdots.\]

The size of the Krylov subspace is changed dynamically during the integration. The Krylov subspace is computed using the incomplete orthogonalization method.

Arguments

• tau_out: scalar or ordered collection of output times. All entries must have the same sign.
• A: square matrix or linear operator used as the phi-function argument.
• u: vector u_0, or matrix whose columns are u_0, u_1, ... in the linear combination above.

Keywords

• tol::Real = 1e-7: requested local convergence tolerance.
• mmin::Int = 10: minimum adaptive Krylov dimension.
• mmax::Int = 128: maximum adaptive Krylov dimension.
• m::Int = min(mmin, mmax): initial Krylov dimension estimate.
• iop::Int = 2: incomplete-orthogonalization length.
• ishermitian::Bool = LinearAlgebra.ishermitian(A): whether A is Hermitian. Only set this to true when the operator obeys that contract.
• opnorm = LinearAlgebra.opnorm(A, Inf): operator-norm scale used by the adaptive controller.
• task1::Bool = false: apply the KIOPS task-one final scaling when true.

Returns

Returns (w, stats), where w contains the requested linear combinations and stats is the tuple (substeps, rejected_steps, krylov_steps, matrix_exponentials, final_krylov_dimension).

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
w, stats = kiops(0.1, A, b)

References

• Gaudreault, S., Rainwater, G. and Tokman, M., 2018. KIOPS: A fast adaptive Krylov subspace solver for exponential integrators. Journal of Computational Physics. Based on the PHIPM and EXPMVP codes (http://www1.maths.leeds.ac.uk/~jitse/software.html). https://gitlab.com/stephane.gaudreault/kiops.
• Niesen, J. and Wright, W.M., 2011. A Krylov subspace method for option pricing. SSRN 1799124
• Niesen, J. and Wright, W.M., 2012. Algorithm 919: A Krylov subspace algorithm for evaluating the $φ$-functions appearing in exponential integrators. ACM Transactions on Mathematical Software (TOMS), 38(3), p.22

ExpvCache{T}(maxiter::Int)

Reusable scratch workspace for the in-place expv! Krylov matrix-exponential-vector product with element type T. It holds a single flat buffer sized for an maxiter×maxiter Hessenberg matrix; pass the cache as the cache keyword to expv! to avoid reallocating that buffer on repeated calls. The buffer is grown automatically (via resize!) if a later call requests a larger subspace than the one it was allocated for.

Arguments

• T: element type of the Krylov Hessenberg workspace.
• maxiter: largest Krylov dimension expected for repeated calls.

Returns

An ExpvCache sized for Krylov dimensions through maxiter.

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
cache = ExpvCache{Float64}(30)
expv!(similar(b), 0.1, arnoldi(A, b); cache)

Fields

• mem::Vector{T}: flat storage of length maxiter^2 reshaped on demand into the m×m working copy of the Hessenberg matrix.
• expcache: typed, size-keyed reduced matrix-exponential workspaces.
• expcol::Vector{T}: contiguous storage for the first column of the reduced matrix exponential.

PhivCache(w, maxiter::Int, p::Int)

Reusable scratch workspace for the in-place phiv! Krylov matrix-phi-vector product. It packs all the intermediate buffers needed to evaluate ϕ_0(tA)b … ϕ_p(tA)b over a Krylov subspace of dimension up to maxiter into a single flat vector, whose element type matches eltype(w) (the output array w). Pass the cache as the cache keyword to phiv! to avoid reallocating these buffers on repeated calls; the buffer is grown automatically (via resize!) if a later call requests a larger subspace or higher order than the one it was allocated for.

The useview type parameter records whether views into the flat buffer are usable (true for ordinary Arrays, false for GPU arrays, which get freshly allocated reshaped copies instead).

Arguments

• w: representative output array whose element type and device determine the workspace layout.
• maxiter: largest Krylov dimension expected for repeated calls.
• p: highest phi-function order expected for repeated calls.

Returns

A PhivCache sized for Krylov dimensions through maxiter and phi orders through p.

Examples

A = [-2.0 1.0; 0.0 -1.0]
b = [1.0, 0.0]
cache = PhivCache(similar(b, length(b), 3), 30, 2)
phiv!(similar(b, length(b), 3), 0.1, arnoldi(A, b), 2; cache)

Fields

• mem::Vector{T}: flat storage that is carved (by get_caches) into the subspace vector, a Hessenberg working copy, and the two augmented matrices used by the phi-function recurrence.
• expcache: exponential! workspaces (see alloc_mem) keyed by extended-matrix size, reused across calls. W is a single concrete workspace type; see ExpvCache for why one type covers all sizes.
• coeffs::Vector{T}: coefficient scratch used by timestep evaluations.
• ts1::Vector{Float64}: one-element time buffer used by scalar-time wrappers.

StegrCache(T, n::Integer)

Subspace-exponential cache for the error-estimate variant of expv! (the mode = :error_estimate path) on Hermitian operators with element type T. It is a concrete HermitianSubspaceCache sized for a Lanczos subspace of dimension up to n. Every iteration diagonalizes the tridiagonal subspace matrix with a preallocated implicit-QL eigensolver, so applying the subspace exponential allocates nothing. Construct one directly, or let get_subspace_cache build the right cache for a given KrylovSubspace.

Arguments

• T: element type of the propagated subspace vector.
• n::Integer: maximum Lanczos subspace dimension.

Returns

A StegrCache sized for subspaces through dimension n.

Fields

• v::Vector{T}: the subspace-propagated vector (length n) that holds the result of applying the subspace exponential.
• w::Vector{T}: scratch vector (length n) for intermediate values.
• d::Vector{R}, e::Vector{R}: working copies (length n) of the tridiagonal diagonal and off-diagonal, overwritten by the eigensolver, where R = real(T).
• Z::Matrix{R}: n × n eigenvector workspace.

Examples

cache = StegrCache(ComplexF64, 30)

get_subspace_cache(Ks::KrylovSubspace) -> SubspaceCache

Construct the subspace-exponential cache appropriate for the Krylov subspace Ks, for use with the error-estimate variant of expv!. For a real (Hermitian) subspace this returns a StegrCache sized to Ks.maxiter, which diagonalizes the tridiagonal subspace matrix with a preallocated implicit-QL eigensolver. Non-Hermitian (complex) subspaces are not yet supported and raise an error.

Arguments

• Ks: a populated or reusable KrylovSubspace with real-valued Hessenberg coefficients.

Returns

A StegrCache that can be supplied to the error-estimate expv! method. This function is a cache constructor, not an extension point.

Examples

Ks = KrylovSubspace{ComplexF64, Float64}(100, 30)
cache = get_subspace_cache(Ks)