Expv: Matrix Exponentials Times Vectors

The main functionality of ExponentialUtilities is the computation of matrix-phi-vector products. The phi functions are defined as

ϕ_0(z) = exp(z)
ϕ_(k+1)(z) = (ϕ_k(z) - 1) / z

In exponential algorithms, products in the form of ϕ_m(tA)b is frequently encountered. Instead of computing the matrix function first and then computing the matrix-vector product, the common alternative is to construct a Krylov subspace K_m(A,b) and then approximate the matrix-phi-vector product.

Support for matrix-free operators

You can use any object as the "matrix" A as long as it implements the following linear operator interface:

  • Base.eltype(A)
  • Base.size(A, dim)
  • LinearAlgebra.mul!(y, A, x) (for computing y = A * x in place).
  • LinearAlgebra.opnorm(A, p=Inf). If this is not implemented or the default implementation can be slow, you can manually pass in the operator norm (a rough estimate is fine) using the keyword argument opnorm.
  • LinearAlgebra.ishermitian(A). If this is not implemented or the default implementation can be slow, you can manually pass in the value using the keyword argument ishermitian.

Core API

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

Compute the matrix-exponential-vector product using Krylov.

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.

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

Compute the matrix-phi-vector products using Krylov. k >= 1.

The phi functions are defined as

\[\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1}{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.

the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631. Formula (10).

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

Non-allocating version of expv that uses precomputed Krylov subspace Ks.

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

Alternative interface for calculating the action of exp(t*A) on the vector b, storing the result in w. The Krylov iteration is terminated when an error estimate for the matrix exponential in the generated subspace is below the requested tolerance. Ks is a KrylovSubspace and typeof(cache)<:HermitianSubspaceCache, the exact type decides which algorithm is used to compute the subspace exponential.

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

Non-allocating version of 'phiv' that uses precomputed Krylov subspace Ks.

Missing docstring.

Missing docstring for exp_timestep. Check Documenter's build log for details.

ExponentialUtilities.phiv_timestepFunction
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 Neisen & Wright. If adaptive==true, the time step and Krylov subsapce 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 parameter 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.

evaluating the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.

Missing docstring.

Missing docstring for exp_timestep!. Check Documenter's build log for details.

ExponentialUtilities.phiFunction
phi(z,k[;cache]) -> [phi_0(z),phi_1(z),...,phi_k(z)]

Compute the scalar phi functions for all orders up to k.

The phi functions are defined as

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

Instead of using the recurrence relation, which is numerically unstable, a formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-156. Theorem 1).

phi(A,k[;cache]) -> [phi_0(A),phi_1(A),...,phi_k(A)]

Compute the matrix phi functions for all orders up to k. k >= 1.

The phi functions are defined as

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

Calls phiv_dense on each of the basis vectors to obtain the answer. If A is Diagonal, instead calls the scalar phi on each diagonal element and the return values are also Diagonals

Caches

Missing docstring.

Missing docstring for ExpvCache. Check Documenter's build log for details.

Missing docstring.

Missing docstring for PhivCache. Check Documenter's build log for details.

  • 1Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating
  • 1Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for