Skip to content

add support to linear maps #30

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
acroy merged 3 commits into from
Aug 9, 2018
Merged

add support to linear maps #30

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f8a53a3
test linear operator
GiggleLiu Jul 27, 2018
9d6f284
update anorm, throw warning.
GiggleLiu Aug 7, 2018
8254cb4
minor change
GiggleLiu Aug 9, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 21 additions & 7 deletions src/expmv.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,22 @@
export expmv, expmv!

function default_anorm(A)
try
Compat.opnorm(A, Inf)
catch err
if err isa MethodError
warn("opnorm($(typeof(A)), Inf) is not defined, fall back to using `anorm = 1.0`.
To suppress this warning, please specify the anorm parameter manually.")
1.0
else
throw(err)
end
end
end


"""
expmv{T}(t, A, vec; [tol], [m], [norm])
expmv{T}(t, A, vec; [tol], [m], [norm], [anorm])

Calculate matrix exponential acting on some vector, ``w = e^{tA}v``,
using the Krylov subspace approximation.
Expand All @@ -13,20 +28,20 @@ function expmv{T}( t::Number,
A, vec::Vector{T};
tol::Real=1e-7,
m::Int=min(30, size(A, 1)),
norm=Base.norm)
norm=Base.norm, anorm=default_anorm(A))
result = convert(Vector{promote_type(eltype(A), T, typeof(t))}, copy(vec))
expmv!(t, A, result; tol=tol, m=m, norm=norm)
expmv!(t, A, result; tol=tol, m=m, norm=norm, anorm=anorm)
return result
end

expmv!{T}( t::Number,
A, vec::Vector{T};
tol::Real=1e-7,
m::Int=min(30,size(A,1)),
norm=Base.norm) = expmv!(vec, t, A, vec; tol=tol, m=m, norm=norm)
norm=Base.norm, anorm=default_anorm(A)) = expmv!(vec, t, A, vec; tol=tol, m=m, norm=norm, anorm=anorm)

function expmv!{T}( w::Vector{T}, t::Number, A, vec::Vector{T};
tol::Real=1e-7, m::Int=min(30,size(A,1)), norm=Base.norm)
function expmv!{T}(w::Vector{T}, t::Number, A, vec::Vector{T};
tol::Real=1e-7, m::Int=min(30,size(A,1)), norm=Base.norm, anorm=default_anorm(A))

if size(vec,1) != size(A,2)
error("dimension mismatch")
Expand All @@ -39,7 +54,6 @@ function expmv!{T}( w::Vector{T}, t::Number, A, vec::Vector{T};
btol = 1e-7 # tolerance for "happy-breakdown"
maxiter = 10 # max number of time-step refinements

anorm = norm(A, Inf)
rndoff= anorm*eps()

# estimate first time-step and round to two significant digits
Expand Down
11 changes: 5 additions & 6 deletions src/phimv.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,20 +35,20 @@ function phimv{T}( t::Number,
A, u::Vector{T}, vec::Vector{T};
tol::Real=1e-7,
m::Int=min(30, size(A, 1)),
norm=Base.norm)
norm=Base.norm, anorm=default_anorm(A))
result = convert(Vector{promote_type(eltype(A), T, typeof(t))}, copy(vec))
phimv!(t, A, u, result; tol=tol, m=m, norm=norm)
phimv!(t, A, u, result; tol=tol, m=m, norm=norm, anorm=anorm)
return result
end

phimv!{T}( t::Number,
A, u::Vector{T}, vec::Vector{T};
tol::Real=1e-7,
m::Int=min(30, size(A, 1)),
norm=Base.norm) = phimv!(vec, t, A, u, vec; tol=tol, m=m, norm=norm)
norm=Base.norm, anorm=default_anorm(A)) = phimv!(vec, t, A, u, vec; tol=tol, m=m, norm=norm, anorm=anorm)

function phimv!{T}( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
tol::Real=1e-7, m::Int=min(30, size(A, 1)), norm=Base.norm)
tol::Real=1e-7, m::Int=min(30, size(A, 1)), norm=Base.norm, anorm=default_anorm(A))

if size(vec, 1) != size(A, 2)
error("dimension mismatch")
Expand All @@ -61,7 +61,6 @@ function phimv!{T}( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
btol = 1e-7 # tolerance for "happy-breakdown"
maxiter = 10 # max number of time-step refinements

anorm = norm(A, Inf)
rndoff = anorm*eps()

# estimate first time-step and round to two significant digits
Expand Down Expand Up @@ -94,7 +93,7 @@ function phimv!{T}( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
scale!(copy!(vm[1], A*w+u), 1/beta)

for j = 1:m
A_mul_B!(p, A, vm[j])
Base.A_mul_B!(p, A, vm[j])

for i = 1:j
hm[i, j] = dot(vm[i], p)
Expand Down
74 changes: 74 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,22 @@
using Expokit
using Base.Test

struct LinearOp
m
end
Copy link
Contributor

@garrison garrison Aug 10, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

wouldn't this be better as:

struct LinearOp{T}
    m::T
end

Unless there is a recent development I am missing, I believe m will be of type Any without this change.

Copy link
Contributor Author

@GiggleLiu GiggleLiu Aug 10, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right, thanks for pointing this out, add a type parameter can increase type stability a lot.
Since this LinearOp is only used for testing and this branch has been merged, I probability won't fix it right now. Users are suggested to use the LinearMaps.jl package. LinearMap type in this package is much more powerful and supports many abstract operations.

Copy link
Contributor

@garrison garrison Aug 10, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. Let me explain a bit more why I asked. Since there is no documentation in this pull request, I looked at the code because I was a bit confused how to use this functionality. Although the pull request mentions LinearMaps, I actually mistakenly started to think that support for LinearMaps had been removed during the review process of this pull request in favor of LinearOp.

That said, why not add LinearMaps to test/REQUIRE and use it in the tests if the goal is to support it?

Copy link
Owner

@acroy acroy Aug 11, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the title of the PR was a bit misleading. Basically, a new keyword anormwas added which also lifts the need for Ato support norm/opnorm and hence allows using LinearMaps. As @GiggleLiu says the test is not really an example for using LinearMaps.

Maybe it would be a good a idea to add such an example to README or somewhere else, though?

Copy link
Contributor Author

@GiggleLiu GiggleLiu Aug 11, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is a good idea to add some examples. Also, it is the time to maintain a document using Documenter, rather than adding more examples into README.

Copy link
Owner

@acroy acroy Aug 11, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@GiggleLiu Sure, it is on the list. OTOH there isn't that much to document ...


@static if VERSION < v"0.7-"
Base.A_mul_B!(y, lo::LinearOp, x) = A_mul_B!(y, lo.m, x)
else
import LinearAlgebra: mul!, A_mul_B!
mul!(y, lo::LinearOp, x) = mul!(y, lo.m, x)
A_mul_B!(y, lo::LinearOp, x) = A_mul_B!(y, lo.m, x)
end
Base.size(lo::LinearOp, i::Int) = size(lo.m, i)
Base.eltype(lo::LinearOp) = eltype(lo.m)
import Base: *
*(lo::LinearOp, v::Vector) = Base.A_mul_B!(similar(v), lo, v)

function test_expmv(n::Int)

A = sprand(n,n,0.4)
Expand Down Expand Up @@ -40,6 +56,30 @@ function test_expmv3()
return e1, e2
end

function test_expmv_linop(n::Int)
A = LinearOp(sprand(n,n,0.2) + 1im*sprand(n,n,0.2))
v = eye(n,1)[:]+0im*eye(n,1)[:]

tic()
w1 = expmv(1.0, A, v, anorm=norm(A.m, Inf))
t1 = toc()

tic()
w2 = expmv(1.0, A, v)
t2 = toc()

tic()
w3 = expm(full(A.m))*v
t3 = toc()

return max(norm(w1-w2)/norm(w2), norm(w1-w3)/norm(w3)), t1, t2, t3
end

println("testing linear operator n=1000 (first expmv, then expm)")
res, t1, t2, t3 = test_expmv_linop(1000)
println("residuum: $res\n")
@test res < 1e-6

println("testing real n=100 (first expmv, then expm)")
res, t1, t2 = test_expmv(100)
println("residuum: $res\n")
Expand Down Expand Up @@ -247,3 +287,37 @@ println("testing complex n=1000 (first phimv, then expm)")
res, t1, t2 = test_phimv2(1000)
println("residuum: $res\n")
@test res < 1e-6

function test_phimv_linop(n::Int)
p = 0.1
found = false
A = LinearOp(sprand(n,n,p))
u = rand(n)
x = similar(u)
while !found
try
copy!(x, A.m\u)
found = true
catch
A = LinearOp(sprand(n,n,p))
end
end
vec = eye(n, 1)[:]

w1 = phimv(1.0, A, u, vec) # warmup
tic()
w1 = phimv(1.0, A, u, vec)
t1 = toc()

w2 = expm(full(A.m))*(vec+x)-x # warmup
tic()
w2 = expm(full(A.m))*(vec+x)-x
t2 = toc()

return norm(w1-w2)/norm(w2), t1, t2
end

println("testing real n=1000 (first phimv, then expm), the linear operator version.")
res, t1, t2 = test_phimv_linop(1000)
println("residuum: $res\n")
@test res < 1e-6