Skip to content

Instantly share code, notes, and snippets.

@MSeeker1340

MSeeker1340/expRK.md Secret

Created May 25, 2018 19:54
Show Gist options
  • Save MSeeker1340/5f5f9f1333824ee7aa7052b98a0e6dc2 to your computer and use it in GitHub Desktop.
Save MSeeker1340/5f5f9f1333824ee7aa7052b98a0e6dc2 to your computer and use it in GitHub Desktop.
Download ZIP
Accuracy and performance comparison with Expokit
Raw
expRK.md 
import Expokit
import OrdinaryDiffEq
expmv_xk, expmv!_xk = Expokit.expmv, Expokit.expmv!
phimv_xk, phimv!_xk = Expokit.phimv, Expokit.phimv!
expmv_od, expmv!_od = OrdinaryDiffEq.expmv, OrdinaryDiffEq.expmv!
phimv_od, phimv!_od = OrdinaryDiffEq.phimv, OrdinaryDiffEq.phimv!;

Accuracy/Numerical Stability

n = 100
srand(0)
Q, _ = qr(randn(n,n))
lambda_A = linspace(-1.0, -0.1, n) # negative real
lambda_B = 1im * linspace(-1.0, 1.0, n) # imaginary
A = Q' * diagm(lambda_A) * Q
B = Q' * diagm(lambda_B) * Q
b = normalize(randn(n));
function test_error(L, b, h, m)
    exphL = expm(h*L)
    phihL = eltype(L).(big.(h*L) \ (big.(exphL) - I))
    w_exact = [exphL*b phihL*b];
    w_xk = [expmv_xk(h, L, b; m=m) phimv_xk(h, L, b, zeros(eltype(L), length(b)); m=m)/h]
    w_od = phimv_od(h, L, b, 1; m=m)
    errs = zeros(2,2)
    println("exp(hA)*b:")
    println("Expokit: ", norm(w_xk[:,1] - w_exact[:,1]) / norm(w_exact[:,1]), 
        "; OrdinaryDiffEq: ", norm(w_od[:,1] - w_exact[:,1]) / norm(w_exact[:,1]))
    println("phi(hA)*b:")
    println("Expokit: ", norm(w_xk[:,2] - w_exact[:,2]) / norm(w_exact[:,2]), 
        "; OrdinaryDiffEq: ", norm(w_od[:,2] - w_exact[:,2]) / norm(w_exact[:,2]))
end;

Real, negative spectrum

# Exact case (m = n)
for h in [1e-4, 1e-2, 1.0, 1e2, 1e4]
    println("# h = $h")
    test_error(A, b, h, n)
end
# h = 0.0001
exp(hA)*b:
Expokit: 2.2837617949354605e-16; OrdinaryDiffEq: 2.2837617949354605e-16
phi(hA)*b:
Expokit: 3.2923478551037227e-12; OrdinaryDiffEq: 3.2923455719465234e-12
# h = 0.01
exp(hA)*b:
Expokit: 2.54921914652591e-16; OrdinaryDiffEq: 2.54921914652591e-16
phi(hA)*b:
Expokit: 4.343341047142426e-14; OrdinaryDiffEq: 4.344095644189768e-14
# h = 1.0
exp(hA)*b:
Expokit: 3.332671455733256e-16; OrdinaryDiffEq: 4.375790537888373e-16
phi(hA)*b:
Expokit: 7.076343297022666e-16; OrdinaryDiffEq: 7.126561963869579e-16
# h = 100.0
exp(hA)*b:
Expokit: 8.297226303687405e-15; OrdinaryDiffEq: 8.103911446486427e-15
phi(hA)*b:
Expokit: 1.3365491527176778e-15; OrdinaryDiffEq: 6.197558339535075e-15
# h = 10000.0
exp(hA)*b:
Expokit: Inf; OrdinaryDiffEq: Inf
phi(hA)*b:
Expokit: 3.676650809868504e-16; OrdinaryDiffEq: 9.090183355668557e-13

m = 20
for h in [1e-4, 1e-2, 1.0, 1e2, 1e4]
    println("# h = $h")
    test_error(A, b, h, m)
end
# h = 0.0001
exp(hA)*b:
Expokit: 2.2837617949354605e-16; OrdinaryDiffEq: 2.2837617949354605e-16
phi(hA)*b:
Expokit: 3.2923478551037227e-12; OrdinaryDiffEq: 3.2923455719465234e-12
# h = 0.01
exp(hA)*b:
Expokit: 2.54921914652591e-16; OrdinaryDiffEq: 2.54921914652591e-16
phi(hA)*b:
Expokit: 4.3434058969675706e-14; OrdinaryDiffEq: 4.344095644189768e-14
# h = 1.0
exp(hA)*b:
Expokit: 3.530172695758857e-16; OrdinaryDiffEq: 3.453365244190219e-16
phi(hA)*b:
Expokit: 7.010321258908902e-16; OrdinaryDiffEq: 7.063207972186391e-16
# h = 100.0
exp(hA)*b:
Expokit: 3.3567087968792793e-9; OrdinaryDiffEq: 0.01020547120944026
phi(hA)*b:
Expokit: 2.4310063674268327e-14; OrdinaryDiffEq: 2.122741533158148e-6
# h = 10000.0
exp(hA)*b:
Expokit: Inf; OrdinaryDiffEq: NaN
phi(hA)*b:
Expokit: 2.891055892711585e-16; OrdinaryDiffEq: 2.2790997530279345e-6

m = 5
for h in [1e-4, 1e-2, 1.0, 1e2, 1e4]
    println("# h = $h")
    test_error(A, b, h, m)
end
# h = 0.0001
exp(hA)*b:
Expokit: 2.2837617949354605e-16; OrdinaryDiffEq: 2.2837617949354605e-16
phi(hA)*b:
Expokit: 3.2923478551037227e-12; OrdinaryDiffEq: 3.2923455719465234e-12
# h = 0.01
exp(hA)*b:
Expokit: 2.5489950597043657e-16; OrdinaryDiffEq: 6.77778150554718e-16
phi(hA)*b:
Expokit: 4.343454527111225e-14; OrdinaryDiffEq: 4.3441440069826416e-14
# h = 1.0
exp(hA)*b:
Expokit: 3.6079660557090645e-8; OrdinaryDiffEq: 5.996964880805403e-6
phi(hA)*b:
Expokit: 6.467693173042081e-8; OrdinaryDiffEq: 8.339329832456894e-7
# h = 100.0
exp(hA)*b:
Expokit: 0.014751062777650174; OrdinaryDiffEq: 0.9531202611175952
phi(hA)*b:
Expokit: 4.082717420602153e-9; OrdinaryDiffEq: 0.0526574814482458
# h = 10000.0



Number of iterations exceeded 10. Requested tolerance might be too high.



Stacktrace:

 [1] #phimv!#13(::Float64, ::Int64, ::Base.LinAlg.#norm, ::Function, ::Array{Float64,1}, ::Float64, ::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\MSeeker\.julia\v0.6\Expokit\src\phimv.jl:156

 [2] (::Expokit.#kw##phimv!)(::Array{Any,1}, ::Expokit.#phimv!, ::Array{Float64,1}, ::Float64, ::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,1}) at .\<missing>:0 (repeats 2 times)

 [3] #phimv#11 at C:\Users\MSeeker\.julia\v0.6\Expokit\src\phimv.jl:40 [inlined]

 [4] (::Expokit.#kw##phimv)(::Array{Any,1}, ::Expokit.#phimv, ::Float64, ::Array{Float64,2}, ::Array{Float64,1}, ::Array{Float64,1}) at .\<missing>:0

 [5] test_error(::Array{Float64,2}, ::Array{Float64,1}, ::Float64, ::Int64) at .\In[3]:5

 [6] macro expansion at .\In[6]:4 [inlined]

 [7] anonymous at .\<missing>:?

Imaginary spectrum

# Exact case (m = n)
for h in [1e-4, 1e-2, 1.0, 1e2, 1e4]
    println("# h = $h")
    test_error(B, Complex128.(b), h, n)
end
# h = 0.0001
exp(hA)*b:
Expokit: 2.792221209481551e-16; OrdinaryDiffEq: 2.354579633760155e-16
phi(hA)*b:
Expokit: 1.7839221788194472e-11; OrdinaryDiffEq: 1.7839221785417142e-11
# h = 0.01
exp(hA)*b:
Expokit: 2.525338692636447e-16; OrdinaryDiffEq: 2.676619123343307e-16
phi(hA)*b:
Expokit: 2.2380182067960868e-13; OrdinaryDiffEq: 2.238014548289598e-13
# h = 1.0
exp(hA)*b:
Expokit: 4.501308712811953e-16; OrdinaryDiffEq: 4.151099202133057e-16
phi(hA)*b:
Expokit: 4.294090524879023e-15; OrdinaryDiffEq: 4.289949067543781e-15
# h = 100.0
exp(hA)*b:
Expokit: 9.195264096758606e-15; OrdinaryDiffEq: 9.533594460301521e-15
phi(hA)*b:
Expokit: 1.0748472046238915e-14; OrdinaryDiffEq: 1.2786141607239076e-14
# h = 10000.0
exp(hA)*b:
Expokit: 8.20421675712255e-13; OrdinaryDiffEq: 8.081649652678786e-13
phi(hA)*b:
Expokit: 7.544302271453979e-13; OrdinaryDiffEq: 6.385399897989036e-13

m = 20
for h in [1e-4, 1e-2, 1.0, 1e2, 1e4]
    println("# h = $h")
    test_error(B, Complex128.(b), h, m)
end
# h = 0.0001
exp(hA)*b:
Expokit: 2.792221209481551e-16; OrdinaryDiffEq: 2.354579633760155e-16
phi(hA)*b:
Expokit: 1.7839221788194472e-11; OrdinaryDiffEq: 1.7839221785417142e-11
# h = 0.01
exp(hA)*b:
Expokit: 2.525338692636447e-16; OrdinaryDiffEq: 2.676619123343307e-16
phi(hA)*b:
Expokit: 2.2380182067960868e-13; OrdinaryDiffEq: 2.238014548289598e-13
# h = 1.0
exp(hA)*b:
Expokit: 5.023067318850921e-16; OrdinaryDiffEq: 4.151099202133057e-16
phi(hA)*b:
Expokit: 4.319018619118212e-15; OrdinaryDiffEq: 4.289949067543781e-15
# h = 100.0
exp(hA)*b:
Expokit: 7.559134210593958e-7; OrdinaryDiffEq: 1.3025847446408563
phi(hA)*b:
Expokit: 7.13373834571079e-8; OrdinaryDiffEq: 0.9833435856996274
# h = 10000.0
exp(hA)*b:
Expokit: 8.745618196360159e-5; OrdinaryDiffEq: 1.475495687044899
phi(hA)*b:
Expokit: 1.0471092637417113e-5; OrdinaryDiffEq: 1.0690155336180696

m = 5
for h in [1e-4, 1e-2, 1.0, 1e2, 1e4]
    println("# h = $h")
    test_error(B, Complex128.(b), h, m)
end
# h = 0.0001
exp(hA)*b:
Expokit: 2.792221209481551e-16; OrdinaryDiffEq: 2.354579633760155e-16
phi(hA)*b:
Expokit: 1.7839221788194472e-11; OrdinaryDiffEq: 1.7839221785417142e-11
# h = 0.01
exp(hA)*b:
Expokit: 2.506713212623583e-16; OrdinaryDiffEq: 3.2940492252317043e-14
phi(hA)*b:
Expokit: 2.238018218343006e-13; OrdinaryDiffEq: 2.2389920348710448e-13
# h = 1.0
exp(hA)*b:
Expokit: 4.7176255434700905e-8; OrdinaryDiffEq: 0.0003190647411724159
phi(hA)*b:
Expokit: 4.490530850720507e-8; OrdinaryDiffEq: 5.4549948948981035e-5
# h = 100.0
exp(hA)*b:
Expokit: 5.678315712300729e-6; OrdinaryDiffEq: 1.3822415340887466
phi(hA)*b:
Expokit: 5.684502484896983e-7; OrdinaryDiffEq: 4.445702950026552
# h = 10000.0
exp(hA)*b:
Expokit: 0.0005685779311483509; OrdinaryDiffEq: 1.3660917341845122
phi(hA)*b:
Expokit: 7.233929384412605e-5; OrdinaryDiffEq: 7.4119441316515875

Performance analysis

using OrdinaryDiffEq: KrylovSubspace, arnoldi!
m = 20
h = 1e-3
w = zeros(n)
println("Expokit: ")
@time expmv!_xk(w, h, A, b; m=m)
println("OrdinaryDiffEq: ")
Ks = KrylovSubspace{Float64}(n, m)
cache = zeros(m, m)
@time begin
    arnoldi!(Ks, A, b; m=m)
    expmv!_od(w, h, Ks; cache=cache)
end;
Expokit: 
  0.001109 seconds (53 allocations: 79.563 KiB)
OrdinaryDiffEq: 
  0.000666 seconds (481 allocations: 66.063 KiB)

m = 20
h = 1e3
w = zeros(n)
println("Expokit: ")
@time expmv!_xk(w, h, A, b; m=m)
println("OrdinaryDiffEq: ")
Ks = KrylovSubspace{Float64}(n, m)
cache = zeros(m, m)
@time begin
    arnoldi!(Ks, A, b; m=m)
    expmv!_od(w, h, Ks; cache=cache)
end;
Expokit: 
  0.005979 seconds (372 allocations: 1.101 MiB)
OrdinaryDiffEq: 
  0.000951 seconds (512 allocations: 166.891 KiB)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment