add support to linear maps #30

GiggleLiu · 2018-07-27T19:09:19Z

I changed the interface slightly to weaken the dependancy of norm(A, Inf), in order to

add support to LinearMaps.jl package. LinearMaps maps a matrix vector multiplication to a linear operator.
provide an option to remove the cost of calculating norm of input matrix, which can be expensive for some user defined matrices.

The value of norm will not affect the correctness of result, thus it is proper to take a default value, e.g. 1.0.

Now a linear operator requires the following 3 interfaces

A_mul_B!(y, A, x)
size(A, i::Int)
eltype(A)

acroy · 2018-07-29T13:35:33Z

Thanks! Looks reasonable and useful. Did you do any benchmarks by any chance?

GiggleLiu · 2018-07-29T14:16:14Z

@acroy From the output time in test for a 1000 x 1000 sparse matrix, the time consumption decreased by an order when the program fall back to using default anorm = 1.

Do you know what may explain this counter intuitive result?

acroy · 2018-07-29T16:50:59Z

Hard to say. How did you do the test?

GiggleLiu · 2018-07-29T17:55:09Z

It is in runtests.jl

function test_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

The second run uses anorm = 1

acroy · 2018-07-30T21:12:55Z

You have to run those tests at least twice to get a realistic timing.

Following #29 it might be good to use Compat.opnorm instead for the norm of the operator?

GiggleLiu · 2018-08-02T10:57:51Z

using LinearMaps
using Compat, BenchmarkTools
using Expokit

n = 1000
sp = sprand(n,n,0.01)
A = LinearMap{ComplexF64}(x->sp*x, n)
v = eye(n,1)[:]+0im*eye(n,1)[:]

@benchmark Expokit.expmv(1.0, $A, $v, anorm=$(Compat.opnorm(sp, Inf)))

BenchmarkTools.Trial: 
  memory estimate:  2.65 MiB
  allocs estimate:  1389
  --------------
  minimum time:     7.952 ms (0.00% GC)
  median time:      9.418 ms (0.00% GC)
  mean time:        11.178 ms (0.95% GC)
  maximum time:     69.372 ms (0.00% GC)
  --------------
  samples:          447
  evals/sample:     1

@benchmark Expokit.expmv(1.0, $A, $v, anorm=1.0)

BenchmarkTools.Trial: 
  memory estimate:  1.62 MiB
  allocs estimate:  722
  --------------
  minimum time:     4.034 ms (0.00% GC)
  median time:      4.658 ms (0.00% GC)
  mean time:        5.517 ms (1.32% GC)
  maximum time:     61.518 ms (0.00% GC)
  --------------
  samples:          905
  evals/sample:     1

CPU info: Intel(R) Core(TM) i5-7200U CPU @ 2.50GHz

acroy · 2018-08-02T15:53:29Z

Interesting. So it is more a factor of 2. I guess the algorithm takes different numbers of iterations depending on anorm. What do you get for the norm of the matrix?

It is probably a “feature” of the original algorithm ... Maybe you could open a separate issue for this?

GiggleLiu · 2018-08-03T16:18:56Z

It is somehow related to the stopping criteria, but I haven't thuroully read the reference yet and can not have a conclusion.

I have opened an issue here #31 .

acroy · 2018-08-03T19:24:44Z

Great! I will have a look when I have access to a computer again.

I thought about the default behavior of the new keyword. Wouldn’t it be easier to take norm(A, Inf) as a default value for anorm? If the type of A doesn’t support this method one gets a method error and the user can supply some other value. Moreover, this behavior would be consistent with the original behavior...

GiggleLiu · 2018-08-04T04:00:34Z

I agree that the program should have consistent behavior.
But here, if anorm does not change the result at all, throwing an error is not nessesary. To user, it is still consistent.

We should find out whether anorm is important enough (at least in some specific cases) to let user set a specific value.

GiggleLiu · 2018-08-04T05:29:00Z

@acroy
I am considering using LinearMaps.jl + Expokit.jl to realize time evolution in a project of quantum circuit simulation

https://github.com/QuantumBFS/Yao.jl

It will be helpful if you would consider revising and merging this PR at the time when Julia 1.0 releases.

By the way, do you know the bechmark in https://github.com/matteoacrossi/ExpmV.jl? Although their project is less active, the algorithm looks more advanced.

acroy · 2018-08-06T07:29:15Z

Don't worry I am certainly going to merge this feature. However, I think you misunderstood me: In your PR the default value of anorm is 1 which changes the behavior compared to the original implementation. Moreover, checking if the type supports norm is IMO not necessary since Julia will throw a method error anyways (and it is preferable to throw an error than silently changing something under the hood). Therefore, I suggest to use anorm=Compat.opnorm(A, Inf) as a default. If A doesn't support norm (or opnorm) or it is too expensive the user can supply a value for anorm otherwise everything works as before.

acroy · 2018-08-06T07:37:01Z

Regarding the benchmark of ExpmV: I know about it and we discussed it in #1 and #9. I think one should do more benchmarking with "real life" matrices to get a full picture. For example, Expokit has several parameters which can be tweaked to optimize the performance.

GiggleLiu · 2018-08-07T03:07:33Z

OK, got it. I think we should provide more meaningful error information.
e.g.

ERROR: MethodError: the norm of given matrix can not be obtained,
try to define opnorm(YOUR_TYPE, Inf) for given matrix or provide anorm manually.

is better than

ERROR: MethodError: no method matching opnorm(YOUR_TYPE)

because, with the second error, user still does not know what to do.

Or, we can

WARNING: opnorm(YOUR_TYPE, Inf) is not defined for input matrix,
 fall back to using `anorm = 1.0`, to suppress this warning,
please specify the anorm parameter manually.

Warning or error, Which do you prefer?

acroy · 2018-08-07T09:06:22Z

Fair enough, instead of silently setting anorm to 1.0 you would throw an error with one of the messages you proposed. Please add the same code also for the function phimv.

acroy

Looks fine. I just have few minor remarks.

acroy · 2018-08-07T17:29:23Z

src/expmv.jl

@@ -1,5 +1,20 @@
 export expmv, expmv!

+function default_anorm(A)
+    try
+        norm(A, Inf)


I would suggest Compat.opnorm here.

acroy · 2018-08-07T17:31:36Z

src/expmv.jl

@@ -13,7 +28,7 @@ function expmv{T}( t::Number,
                   A, vec::Vector{T};
                   tol::Real=1e-7,
                   m::Int=min(30, size(A, 1)),
-                   norm=Base.norm, anorm::Real=0)
+                   norm=Base.norm, anorm=nothing)


And here just anorm=default_anorm(A) (and everywhere else as well).

matteoacrossi · 2018-08-08T11:10:00Z

Hi, I'm maintaining https://github.com/matteoacrossi/ExpmV.jl that @GiggleLiu mentioned above. I actually forked from https://github.com/marcusps/ExpmV.jl which was inactive, because I needed to extend the expmv function for a span of the parameter t (btw, does Expokit.jl support this?). Then I thought it would have been a nice idea to release the package for other users as well.

ExpmV is a translation of a more recent algorithm published by Al-Mohy and Higham in 2010, (the original MATLAB implementation is here https://github.com/higham/expmv), which is not based on Krylov subspace projection. In the paper they compare it to Expokit and discuss pros and cons of both algorithms.

I think it would be interesting to implement the same real-life benchmarks for both packages to have a detailed comparison (I've seen that in PR #28 @mforets has added some).

One could even think of merging the two packages (or create a wrapper package) and then give the user the option to choose which algorithm to use.

ChrisRackauckas · 2018-08-08T12:37:43Z

It would also be nice to test the numerical accuracy, and time, on different functions. @MSeeker1340 did a bunch of accuracy tests with ExponentialUtilities.jl during its development in OrdinaryDiffEq.jl and it would be nice to have a comparison to help with further development.

acroy · 2018-08-09T06:57:53Z

@matteoacrossi Great to hear that ExpmV is alive and being take care of. As I already wrote we had some discussions about the Krylov method and the AH approach in #1. Basically, I created Expokit.jl because I needed the functionality back then and I wanted to provide a port of expv to help people to try Julia. Then @mforets stepped up an helped to finish the port of the other functions in Expokit. Anyways, I think there is nothing wrong with having different packages for now - maybe everything converges into ExponentialUtilities.jl? Another idea to unify different approaches was to introduce an Expm type (as a linear operator)...

I like the idea of setting up a (sparate) package/repo to collect benchmarks and compare different approaches and packages. This could be quite interesting and helpful. Maybe we should continue to discuss this in a separate issue (#21) though?

garrison · 2018-08-10T03:48:04Z

test/runtests.jl

@@ -1,6 +1,22 @@
 using Expokit
 using Base.Test

+struct LinearOp
+    m
+end


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.

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.

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?

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?

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.

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

garrison · 2018-08-11T12:08:40Z

Thanks 💯 for the clarifications. Does the merging of this mean that #29 can be closed?

test linear operator

f8a53a3

GiggleLiu mentioned this pull request Jul 27, 2018

Add support to LinearMaps.jl marcusps/ExpmV.jl#12

Open

GiggleLiu mentioned this pull request Aug 3, 2018

Influence of anorm on performance of expmv and phimv? #31

Open

update anorm, throw warning.

9d6f284

acroy reviewed Aug 7, 2018

View reviewed changes

matteoacrossi mentioned this pull request Aug 9, 2018

Benchmarks #21

Open

minor change

8254cb4

acroy merged commit b5e020e into acroy:master Aug 9, 2018

garrison reviewed Aug 10, 2018

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

add support to linear maps #30

add support to linear maps #30

GiggleLiu commented Jul 27, 2018 •

edited

Loading

acroy commented Jul 29, 2018

GiggleLiu commented Jul 29, 2018

acroy commented Jul 29, 2018

GiggleLiu commented Jul 29, 2018

acroy commented Jul 30, 2018 •

edited

Loading

GiggleLiu commented Aug 2, 2018 •

edited

Loading

acroy commented Aug 2, 2018

GiggleLiu commented Aug 3, 2018

acroy commented Aug 3, 2018

GiggleLiu commented Aug 4, 2018

GiggleLiu commented Aug 4, 2018 •

edited

Loading

acroy commented Aug 6, 2018

acroy commented Aug 6, 2018

GiggleLiu commented Aug 7, 2018 •

edited

Loading

acroy commented Aug 7, 2018

acroy left a comment

acroy Aug 7, 2018

acroy Aug 7, 2018

matteoacrossi commented Aug 8, 2018 •

edited

Loading

ChrisRackauckas commented Aug 8, 2018

acroy commented Aug 9, 2018

garrison Aug 10, 2018 •

edited

Loading

GiggleLiu Aug 10, 2018 •

edited

Loading

garrison Aug 10, 2018

acroy Aug 11, 2018

GiggleLiu Aug 11, 2018

acroy Aug 11, 2018

garrison commented Aug 11, 2018

add support to linear maps #30

add support to linear maps #30

Conversation

GiggleLiu commented Jul 27, 2018 • edited Loading

Now a linear operator requires the following 3 interfaces

acroy commented Jul 29, 2018

GiggleLiu commented Jul 29, 2018

acroy commented Jul 29, 2018

GiggleLiu commented Jul 29, 2018

acroy commented Jul 30, 2018 • edited Loading

GiggleLiu commented Aug 2, 2018 • edited Loading

acroy commented Aug 2, 2018

GiggleLiu commented Aug 3, 2018

acroy commented Aug 3, 2018

GiggleLiu commented Aug 4, 2018

GiggleLiu commented Aug 4, 2018 • edited Loading

acroy commented Aug 6, 2018

acroy commented Aug 6, 2018

GiggleLiu commented Aug 7, 2018 • edited Loading

acroy commented Aug 7, 2018

acroy left a comment

Choose a reason for hiding this comment

acroy Aug 7, 2018

Choose a reason for hiding this comment

acroy Aug 7, 2018

Choose a reason for hiding this comment

matteoacrossi commented Aug 8, 2018 • edited Loading

ChrisRackauckas commented Aug 8, 2018

acroy commented Aug 9, 2018

garrison Aug 10, 2018 • edited Loading

Choose a reason for hiding this comment

GiggleLiu Aug 10, 2018 • edited Loading

Choose a reason for hiding this comment

garrison Aug 10, 2018

Choose a reason for hiding this comment

acroy Aug 11, 2018

Choose a reason for hiding this comment

GiggleLiu Aug 11, 2018

Choose a reason for hiding this comment

acroy Aug 11, 2018

Choose a reason for hiding this comment

garrison commented Aug 11, 2018

GiggleLiu commented Jul 27, 2018 •

edited

Loading

acroy commented Jul 30, 2018 •

edited

Loading

GiggleLiu commented Aug 2, 2018 •

edited

Loading

GiggleLiu commented Aug 4, 2018 •

edited

Loading

GiggleLiu commented Aug 7, 2018 •

edited

Loading

matteoacrossi commented Aug 8, 2018 •

edited

Loading

garrison Aug 10, 2018 •

edited

Loading

GiggleLiu Aug 10, 2018 •

edited

Loading