Add Expm module

[schillic]

This can be used in other packages to abstract from how exp is computed.

Whether the FastExpm package is really useful is another question. From experiments it only pays off for very sparse matrices. See the benchmark below.

We should discuss the interface. Currently everything is automatic, but after loading FastExpm, code can get slower. There should be an option to choose the backend even after loading the package. This could be done as in LazySets by calling a function set_expm_backend!. (Maybe we should outsource that functionality from LazySets here as well.)

julia> using ReachabilityBase.Expm, SparseArrays

julia> n = 50;

julia> M = rand(n, n);

julia> M2 = sprand(n, n, 0.1)
50×50 SparseMatrixCSC{Float64, Int64} with 245 stored entries:
⎡⠀⢀⠀⠁⠢⠀⠐⠀⢁⠀⡀⢄⠌⠈⠢⠀⠀⢐⠀⡌⠌⡀⠄⠀⠀⎤
⎢⠀⡀⠀⡂⠀⠀⡀⠀⠀⢠⢀⢀⡀⠀⠀⠀⠂⠀⠀⠀⠀⠢⢀⢐⠈⎥
⎢⠄⠀⠄⠡⠁⠈⠁⠈⠨⠀⢄⠀⢀⠀⠀⠀⠀⡠⠁⠀⠀⠀⠐⠀⠀⎥
⎢⠁⢀⠂⠀⢀⠊⠀⢈⡀⠒⠂⠊⠄⢀⠄⠀⠡⠀⠁⡒⠀⠄⠂⢠⠄⎥
⎢⢀⠁⠊⠀⠃⢂⠀⢀⠀⠤⠀⠀⠁⠠⠀⠀⠠⠀⠀⠀⠁⠀⠀⠀⠄⎥
⎢⠊⠃⠀⠀⢀⠀⠘⢈⠀⡄⠤⠀⠀⢀⠈⣄⠠⢈⠀⠀⠀⠀⠀⠀⡀⎥
⎢⠠⠀⡠⠁⢈⠀⠄⠀⠀⠁⠔⠔⢌⠁⠀⠀⠀⠀⠄⡀⠀⠀⠒⢈⠀⎥
⎢⠀⠀⢀⠀⠀⠀⠀⠀⠈⡀⠦⠄⠀⠌⠀⠀⠀⠀⠁⠀⠅⠰⠀⠰⢂⎥
⎢⢀⢈⠄⠀⢂⣠⠅⠠⠄⠀⠀⠀⠈⠀⠀⠀⠁⠐⠀⠤⢂⠀⠐⠊⠁⎥
⎢⠀⠄⢀⠀⠀⠡⠀⢀⠀⠂⡀⠀⠀⠀⠀⠨⠀⠄⡀⠀⠀⠠⡐⡈⠠⎥
⎢⠀⠀⠀⠢⠮⠀⠀⢀⠉⠀⠑⠀⠀⠀⠀⠀⠠⠀⢁⡀⡂⡀⠂⠀⠀⎥
⎢⠀⠀⡀⠀⠄⢀⠃⠀⠀⠈⠈⠀⠀⠨⠂⠀⠉⠀⠀⠂⠀⠀⠀⠀⠀⎥
⎣⠀⠃⠀⠀⠊⠀⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠐⠀⠀⠀⠀⠂⠀⠈⎦

julia> @time expm(M);  # precompile
  3.402406 seconds (3.49 M allocations: 232.121 MiB, 3.12% gc time, 99.95% compilation time)

julia> @time expm(Matrix(M2));  # precompile
  0.071427 seconds (74.84 k allocations: 5.205 MiB, 99.28% compilation time)

julia> @time expm(M);
  0.000330 seconds (15 allocations: 118.797 KiB)

julia> @time expm(Matrix(M2));
  0.000333 seconds (17 allocations: 138.406 KiB)

julia> import FastExpm

julia> @time expm(M2);  # precompile
  4.967181 seconds (5.09 M allocations: 343.023 MiB, 2.58% gc time, 99.98% compilation time)

julia> @time expm(M2);
  0.000766 seconds (53 allocations: 630.875 KiB)

julia> using ReachabilityBase.Expm, SparseArrays

julia> n = 50
50

julia> M = rand(n, n);

julia> M2 = sprand(n, n, 0.01)
50×50 SparseMatrixCSC{Float64, Int64} with 30 stored entries:
⎡⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠐⠀⠐⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⢀⠀⠂⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⎥
⎢⠀⠀⠀⠐⠀⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠀⠂⠀⠠⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠠⠀⠀⠀⠀⠐⠀⠀⎥
⎣⠀⠀⠈⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠁⠀⠀⠀⠀⠀⠀⠀⎦

julia> @time expm(M);  # precompile
  3.597239 seconds (3.49 M allocations: 232.151 MiB, 3.55% gc time, 99.98% compilation time)

julia> @time expm(Matrix(M2));  # precompile
  0.099280 seconds (74.84 k allocations: 5.205 MiB, 26.58% gc time, 99.53% compilation time)

julia> @time expm(M);
  0.000321 seconds (15 allocations: 118.797 KiB)

julia> @time expm(Matrix(M2));
  0.000413 seconds (17 allocations: 138.406 KiB)

julia> import FastExpm

julia> @time expm(M2);  # precompile
  5.022449 seconds (5.09 M allocations: 342.826 MiB, 4.33% gc time, 99.99% compilation time)

julia> @time expm(M2);
  0.000342 seconds (32 allocations: 353.422 KiB)

[mforets]

Keeping this on hold for now as we investigate how the exp interface should look like.