Numerical stability analysis of phi functions

Phi functions.md

include("Phi.jl")
using Phi
using PyPlot

Utilizing the Sidje formula

For scalar $z$:

$$ \exp\Big( \begin{bmatrix} z & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ 0 & 0 & 0 & \cdots & 0 \end{bmatrix}\Big) = \begin{bmatrix} \exp(z) & \phi_1(z) & \phi_2(z) & \cdots & \phi_p(z)\\ 0 & 1 & \frac{1}{1!} & \cdots & \frac{1}{(p-1)!} \\ 0 & 0 & 1 & \cdots & \frac{1}{(p-2)!} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix}$$

For matrix case, see Expokit paper.

# Scalar
z = 0.1
K = 4
P = phi(z,K)
Palt = phi_alt(z,K)
P ≈ Palt

true

# Matrix
A = randn(2,2)
K = 4
P = phim(A,K)
Palt = phim_alt(A,K)
P ≈ Palt

true

Numerical stability analysis

Scalar case

# Float64
K = 6
zero_points = phi(0.0, K)
zs = logspace(0, -12, 13)
residuals = zeros(K+1, length(zs))
for i in 1:length(zs)
#     phis = phi(zs[i], K)
    phis = phi_alt(zs[i], K)
    @. residuals[:,i] = abs(phis - zero_points)
end

for k = 0:K
    plot(zs, residuals[k+1,:], "-o", label="|ϕ$k(z) - ϕ$k(0)|")
end
loglog()
legend()
xlim(xlim()[2], xlim()[1])
xlabel("z")
title("Float64, scalar, Sidje");

Matrix case

# Float64
K = 4
dim = 5
zero_points = phim(zeros(dim,dim), K)
srand(0); A = randn(dim,dim)
zs = logspace(0, -12, 13)
residuals = zeros(K+1, length(zs))
for i in 1:length(zs)
#     phims = phim(zs[i]*A, K)
    phims = phim_alt(zs[i]*A, K)
    for k in 1:K+1
        residuals[k,i] = norm(phims[k] - zero_points[k], Inf)
    end
end

for k = 0:K
    plot(zs, residuals[k+1,:], "-o", label="|ϕ$k(zA) - ϕ$k(0)|")
end
loglog()
legend()
xlim(xlim()[2], xlim()[1])
xlabel("z")
title("Float64, $dim x $dim matrix, Sidje");

Computational cost

Scalar phi vs phi_alt:

z = 0.1
@time begin
    for i = 1:1000
        phi(z, 4)
    end
end

  0.000130 seconds (2.00 k allocations: 250.000 KiB)

@time begin
    for i = 1:1000
        phi_alt(z, 4)
    end
end

  0.043846 seconds (41.00 k allocations: 10.269 MiB, 5.13% gc time)

z = big(0.1)
@time begin
    for i = 1:1000
        phi(z, 4)
    end
end

  0.015232 seconds (74.00 k allocations: 2.869 MiB)

Matrix phim vs phim_alt

dim = 10
A = randn(dim,dim);

@time begin
    for i = 1:1000
        phim(A, 4)
    end
end

  0.202246 seconds (87.00 k allocations: 49.683 MiB, 8.80% gc time)

@time begin
    for i = 1:1000
        phim_alt(A, 4)
    end
end

  1.011265 seconds (877.00 k allocations: 737.366 MiB, 18.17% gc time)

A = big.(A);

@time begin
    for i = 1:1000
        phim(A, 4)
    end
end

  9.409701 seconds (94.35 M allocations: 3.829 GiB, 29.62% gc time)

Phi.jl

module Phi
using Expokit.padm
export phi, phi_alt, phim, phimc, phim_alt

"""
    phi(z,k) -> [phi_0(z),phi_1(z),...,phi_k(z)]

Compute the value of the phi functions of z for all orders up to k.
"""
function phi(z::T, k) where {T <: Number}
    # Compute zero points
    zero_points = Vector{T}(k+1)
    zero_points[1] = one(T)
    for i = 1:k
        zero_points[i+1] = zero_points[i] / i
    end
    
    if iszero(z)
        return zero_points
    else
        # Compute phi(z) using the recurrence relation
        out = Vector{T}(k+1)
        out[1] = exp(z) # phi_0(z) = exp(z)
        for i = 1:k
            out[i+1] = (out[i] - zero_points[i]) / z
        end
        return out
    end
end

"""
    phi_alt(z,k) -> [phi_0(z),phi_1(z),...,phi_k(z)]

Compute the value of the phi functions of z for all orders up to k.

This alternative implemenation uses the matrix formula and computes all 
phi functions via one computation of `expm`.
"""
function phi_alt(z::T, k) where {T <: Number}
    # Construct the matrix M
    M = zeros(T, k+1, k+1)
    M[1,1] = z
    for i = 1:k
        M[i,i+1] = one(T)
    end
    # Compute expm(M)
    if T == BigFloat || T == Complex{BigFloat}
        P = padm(M)
    else
        P = expm(M)
    end
    return P[1,:]
end

"""
    phim(A,k,[p]) -> [phi_0(A),phi_1(A),...,phi_k(A)]

Compute the value of the matrix phi functions of A for all orders up to k.

Depending on the dtype of A, its matrix exponential is either calculated using 
`expm` (for regular floats) or `Expokit.padm` (for BigFloats). In the latter 
case p will determine the degree of the Pade approximation used.
"""
function phim(A::AbstractMatrix{T}, k; p=6) where {T <: Number}
    @assert size(A,1) == size(A,2) "Dimension mismatch"
    # Compute zero points
    zero_points = Vector{Matrix{T}}(k+1)
    zero_points[1] = eye(A)
    for i = 1:k
        zero_points[i+1] = zero_points[i] ./ i
    end
    
    if iszero(A)
        return zero_points
    else
        # Compute phi(A) using the recurrence relation
        out = Vector{Matrix{T}}(k+1)
        # phi_0(A) = exp(A)
        if T == BigFloat || T == Complex{BigFloat}
            out[1] = padm(A, p=p)
        else
            out[1] = expm(A)
        end
        for i = 1:k
            out[i+1] = A \ (out[i] - zero_points[i])
        end
        return out
    end
end

"""
    phimc(A,c,k) -> [phi_0(A)*c phi_1(A)*c ... phi_k(A)*c]

Compute dense matrix-phi-vector products using the Sidje formula. 
[Sidje, R. B. (1998). Expokit: a software package for computing matrix 
exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 
130-156.]
"""
function phimc(A::AbstractMatrix{T}, c::AbstractVector{T}, k) where {T <: Number}
    @assert size(A,1) == size(A,2) == length(c) "Dimension mismatch"
    m = size(A,1)
    # Construct the extended matrix
    M = zeros(T, m+k, m+k)
    M[1:m, 1:m] = A
    M[1:m, m+1] = c
    for i = m+1:m+k-1
        M[i, i+1] = one(T)
    end
    # Compute expm(M)
    if T == BigFloat || T == Complex{BigFloat}
        P = padm(M)
    else
        P = expm(M)
    end
    # Extract results
    out = zeros(T, m, k+1)
    out[:, 1] = P[1:m, 1:m] * c # expm(A)
    for i = 1:k
        out[:, i+1] = P[1:m, m+i]
    end
    return out
end

"""
    phim_alt(A,k) -> [phi_0(A),phi_1(A),...,phi_k(A)]

Compute phim using phimc with the basis vectors e_j as c.
"""
function phim_alt(A::AbstractMatrix{T}, k) where {T <: Number}
    m = size(A,1)
    out = [Matrix{T}(m, m) for i = 1:k+1]
    E = eye(A)
    for i = 1:m
        e = E[:, i]
        P = phimc(A, e, k) # [phi_0(A)*e phi_1(A)*e ... phi_k(A)*e]
        for j = 1:k+1
            out[j][:, i] = P[:, j]
        end
    end
    return out
end

end