search

acroy / Expokit.jl

Julia implementation of EXPOKIT routines
Other
26 stars 11 forks source link

Benchmarks #21

Open acroy opened 6 years ago

acroy commented 6 years ago

As briefly discussed in #9 and #18 having benchmarks of a couple of "real world problems" would be nice. Some suitable examples are provided in this paper, among them the 3D quantum harmonic oscillator. Additionally, @mforets suggested to "use Hurwitz-stable matrices arising from control systems".

New benchmarks should go into the folder benchmarks/ and this issue can be used for discussions.

mforets commented 6 years ago

Do you have any preferred data format? I have been using MAT.jl. There is also HDF5.jl.

mforets commented 6 years ago

I'll try PkgBenchmark.jl. Do you have experience with it? It allows to define benchmark groups and benchmark between different branches.

acroy commented 6 years ago

Under the hood MATLABs mat format also uses hdf5, so from that perspective it probably doesn’t matter. If we stick to Julia, we could also use JLD. For sparse matrices there is also HarwellRutherfordBoeing.jl, which should be able to read the matrices provided with EXPOKIT.

PkgBenchmark.jl sounds very nice - I haven’t tried it yet.

mforets commented 6 years ago

See PR #28 for a benchmark script using PkgBenchmark. I used JLD to store the matrices as you suggested.

matteoacrossi commented 6 years ago

As we disccused in #30, I'd be interested in having a common benchmark suite for the package ExpmV.jl, that uses a different algorithm. I'm not really an expert of matrix exponentiation and I don't know how their speed and accuracy depends on the features of the matrix, so it would be interesting to have a comprehensive suite, including real-life matrices from different fields.

I don't know what would be the best way to implement this (a new package? or a new repository?) so that both packages can take advantage of the benchmark suite.

GiggleLiu commented 6 years ago

As a real world application, please consider the householder reflection 2.*v*v' - eye(N). Normally, one should use its low rank decomposition like

H(x) = 2 .* v*(v'*x) - x

which is a linear operator.

I really hope ExpmV can also support a general linear operator as the input too.

acroy commented 6 years ago

@matteoacrossi I think a separate package would make sense - possibly starting from #28. Everyone interested could add their method(s), run the benchmark suite and push the results. We would have to find a way to visualize the results somehow, but that is a separate issue.