SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
521 stars 198 forks source link

Use `ExpMethodDiagonalization` rather than `Generic` #2234

Closed oscardssmith closed 1 month ago

oscardssmith commented 1 month ago

It's ~2x faster.

Checklist

Additional context

Add any other context about the problem here.

oscardssmith commented 1 month ago

I think this isn't the way forward. Timing indicates that diagonalization is only faster for matrices bigger than ~512x512 (at least on my computer). I think the better solution needs to come from the ExponentialUtilities side which needs a ExpMethodDefault method which uses a good method based on the problem sparsity and size.

julia> for i in 1:10
           t_generic = @belapsed exponential!(A, ExpMethodGeneric()) setup=(A=rand(2^$i, 2^$i)) evals=1
           t_diag = @belapsed exponential!(A, ExpMethodDiagonalization()) setup=(A=rand(2^$i, 2^$i)) evals=1
           @show i, t_diag, t_generic
       end
(i, t_diag, t_generic) = (1, 2.582e-6, 5.84e-7)
(i, t_diag, t_generic) = (2, 5.613e-6, 1.448e-6)
(i, t_diag, t_generic) = (3, 1.5839e-5, 2.579e-6)
(i, t_diag, t_generic) = (4, 5.5466e-5, 9.875e-6)
(i, t_diag, t_generic) = (5, 0.000238276, 5.8513e-5)
(i, t_diag, t_generic) = (6, 0.00131012, 0.000436813)
(i, t_diag, t_generic) = (7, 0.009235047, 0.003087585)
(i, t_diag, t_generic) = (8, 0.039572306, 0.023596429)
(i, t_diag, t_generic) = (9, 0.24952722, 0.260065999)
(i, t_diag, t_generic) = (10, 0.99962065, 2.685850028)
oscardssmith commented 1 month ago

On the other hand, expv! seems to dominate beyond 10x10, so we possibly should just always use that.

julia> for i in 1:10
           t_generic = @belapsed exponential!(A, ExpMethodGeneric()) setup=(A=rand(2^$i, 2^$i)) evals=1
           t_diag = @belapsed exponential!(A, ExpMethodDiagonalization()) setup=(A=rand(2^$i, 2^$i)) evals=1
           t_expv = @belapsed expv(1.0, A, b) setup=(A,b) = (rand(2^$i, 2^$i), rand(2^$i)) evals=1
           @show i, t_diag, t_generic, t_expv
       end
(i, t_diag, t_generic, t_expv) = (1, 2.676e-6, 6.02e-7, 1.604e-6)
(i, t_diag, t_generic, t_expv) = (2, 5.685e-6, 1.432e-6, 2.822e-6)
(i, t_diag, t_generic, t_expv) = (3, 1.5297e-5, 2.552e-6, 1.1585e-5)
(i, t_diag, t_generic, t_expv) = (4, 5.5605e-5, 9.664e-6, 3.1198e-5)
(i, t_diag, t_generic, t_expv) = (5, 0.000244218, 5.7904e-5, 7.4628e-5)
(i, t_diag, t_generic, t_expv) = (6, 0.00130707, 0.000436641, 8.3198e-5)
(i, t_diag, t_generic, t_expv) = (7, 0.009175025, 0.003052096, 0.000115425)
(i, t_diag, t_generic, t_expv) = (8, 0.039352406, 0.023071895, 0.000252655)
(i, t_diag, t_generic, t_expv) = (9, 0.243172999, 0.25509279, 0.001513738)
(i, t_diag, t_generic, t_expv) = (10, 0.889689215, 2.469432628, 0.003373091)