GitHub

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 are frequently encountered. Instead of computing the matrix function first and then computing the matrix-vector product, the common alternative is to construct a Krylov subspaceK_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.expv — Function
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.

source
ExponentialUtilities.phiv — Function
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.

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

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

source
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.

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

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

source
ExponentialUtilities.expv_timestep — Function
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).

source
ExponentialUtilities.phiv_timestep — Function
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.

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.

source

Caches

Missing docstring.

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

Missing docstring.

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

  • 1Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631. Formula (10).
  • 1Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.