Symbolics.jl: Speed of derivative computations vs. SymEngine.jl

[ranocha]

I benchmarked some derivative computations using Symbolics.jl, SymEngine.jl, and SymPy.jl. As expected, SymPy.jl was the slowest solution in my case. However, SymEngine.jl was significantly faster than Symbolics.jl. Here is an example:

julia> using Pkg; Pkg.activate(temp=true); Pkg.add(["SymEngine", "Symbolics", "BenchmarkTools"])

julia> using SymEngine, Symbolics, BenchmarkTools

julia> f(u) = u[2] / (u[1]^2 + u[2]^2)
f (generic function with 1 method)

julia> u_symbolics = @variables u1 u2

julia> expand_derivatives(( Differential(u_symbolics[1])(f(u_symbolics)) ))
-2u1*(u2 / ((u1^2 + u2^2)^2))

julia> @benchmark expand_derivatives($( Differential(u_symbolics[1])(f(u_symbolics)) ))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  19.477 μs …  4.431 ms  ┊ GC (min … max): 0.00% … 99.14%
 Time  (median):     21.305 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   22.915 μs ± 61.455 μs  ┊ GC (mean ± σ):  3.78% ±  1.40%

     ▃█▇▄                                                      
  ▂▃▅█████▆▆▅▅▄▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▂▂▁▁▁▁▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  19.5 μs         Histogram: frequency by time        37.5 μs <

 Memory estimate: 10.97 KiB, allocs estimate: 233.

julia> u_symengine = @vars u1 u2
(u1, u2)

julia> diff((f(u_symengine)), (u_symengine[1]))
-2*u2*u1/(u1^2 + u2^2)^2

julia> @benchmark diff($(f(u_symengine)), $(u_symengine[1]))
BenchmarkTools.Trial: 10000 samples with 7 evaluations.
 Range (min … max):  4.198 μs …  35.084 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.268 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.320 μs ± 510.445 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄█▇▆▅▂▁▁▁▁                                                  ▂
  ████████████▆▆▅▃▅▃▃▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▃▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆▅ █
  4.2 μs       Histogram: log(frequency) by time       6.1 μs <

 Memory estimate: 336 bytes, allocs estimate: 20.

1. Am I doing something wrong with Symbolics.jl?
2. Is there a way for me to get better performance of such derivative calculations using Symbolics.jl?

Originally posted in https://discourse.julialang.org/t/symbolics-jl-speed-of-derivative-computations-vs-symengine-jl/68803; ChrisRackauckas told me to file an issue here

┆Issue is synchronized with this Trello card by Unito

[ranocha]

Even for a much simpler example such as f(u) = u[1] * u[2], I get

julia> @benchmark expand_derivatives($( Differential(u_symbolics[1])(f(u_symbolics)) ))
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.557 μs … 456.999 μs  ┊ GC (min … max): 0.00% … 99.02%
 Time  (median):     2.693 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.885 μs ±   4.612 μs  ┊ GC (mean ± σ):  1.57% ±  0.99%

  ▆██▇▅▅▅▄▃▁                     ▁                            ▂
  ██████████▇▆▆▆▇▇▅▄▄▁▁▁▁▁▃▆██▇▇▇██▇▅▅▄▃▅▄▃▅▁▃▃▃▁▄▃▃▃▅▆▅▅▇▇▇▇ █
  2.56 μs      Histogram: log(frequency) by time      5.94 μs <

 Memory estimate: 1.27 KiB, allocs estimate: 36.

julia> @benchmark diff($(f(u_symengine)), $(u_symengine[1]))
BenchmarkTools.Trial: 10000 samples with 39 evaluations.
 Range (min … max):  907.000 ns …   3.526 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     936.462 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   961.282 ns ± 149.500 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃█▃▄                                                           
  ████▅▄▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▂▂▂▂▂▂▂▂▂▂▁▁▁▂▂▂▂▁▂▂▁▂▁▂▂▂▂▁▂▂▂▂▂▂▂▂▂ ▂
  907 ns           Histogram: frequency by time         1.75 μs <

 Memory estimate: 96 bytes, allocs estimate: 5.

[ranocha]

See also #397: I have an example of a combination of computing many derivatives, their weighted sum, and substituting some variables in the result. For the complete operation, Symbolics.jl is ca. 70% slower than SymPy.jl and ca. 60x slower than SymEngine.jl.

[ChrisRackauckas]

I wonder if it's essentially the same issue.

[YingboMa]

Nothing is fundamentally slow. It's just that everything is not fast due to the amount of allocations.

[ranocha]

Well, it's so slow that SymPy.jl can be much faster in the end. This benchmark uses differentiation of symbolic expressions and substitution of rational numbers for two variables afterwards. There,

• SymPy.jl needs 2.2 seconds for the differentiation stuff and 0.8 seconds for the substitution
• Symbolics.jl needs 1.4 seconds for the differentiation stuff and 3.7 seconds for the substitution

Thus, it's much better to use SymPy.jl in this case (ignoring SymEngine.jl, which needs only 0.03 seconds for both tasks together). I don't want to say that Symbolics.jl is bad or anything like that - it's a fantastic piece of software and a lot of effort went into it. It's just a case that might help to debug performance issues such that Symbolics.jl can reach the goals set in its official announcement

It is built from the ground up with performance in mind. We use specialized structures for automatic simplification to match the performance of the most fully optimized C++ libraries.

[ChrisRackauckas]

There will always be operations to keep improving, and this seems to be one.

[ChrisRackauckas]

https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/main-65.pdf might be an interesting read.

[YingboMa]

Interesting!

[YingboMa]

I am not sure if that has much relevance in just the symbolic differentiation context.

[JTaets]

Disclaimer: Im not an expert and I might be a 100% wrong here. (I don't even know if something like this is already in symbolics.jl)
I just want to point out that afaik working with graphs is the key to make symbolics.jl fast like how casadi does it (https://www.youtube.com/watch?v=mYOkLkS5yqc especially at 8:09)

For example if you have:
a = long function of x
b = a + a^2 + a^3 + .... + a^100
and you do derivative(b,x), then afaik symbolics.jl will differentiate through a a 100 times
Look at the youtube video for how casadi avoids this, just like forward and backward diff doesn't do this differentiation a 100 times.

Advantages:

• You can easily calculate the symbolic derivative through very complex functions much faster
• If you also keep this graph structure for build_function, it will be MUCH faster! With casadi you can do symbolic differentiation through full ODE solves for example (because it boils down to the same complexity as forward or backwards differentiation)

Casadi doesn't allow simplifications of the symbolics, just so they can keep this efficient graph structure afaik.
In my Opinion, there should be 2 options:

• simplify equations for length (like symbolics.jl does now)
• simplify for speed (like casadi does: don't simplify anything, just keep the graph structure and generate/differentiate functions with this graph structure like shown in the video)

[ChrisRackauckas]

It has a natural graph structure. The issue is more dispatching. Yingbo has a way to fix the dispatching that would improve some speeds by 100x by reducing types.

[JTaets]

Since in my code, doing differentiation was very slow, I assumed it didn't have a graph backend. My apologies!

I looked deeper and my problem lies in build_function generating something very slow (without graph structure), a will be calculated a 100 times in the build function if i'm not wrong again?

For example in the following code, if a graph backend is used for generating the functions, then the following should only do 10 sums, 10 squares and 10 products. But the generated function is very long.
It even crashed my computer when putting a higher bound than 10.

using Symbolics

function integrate(x)
   for i in 1:10
      x+=0.1*x^2
   end
   x
end

x, = @variables x
integrate_symbolic = build_function(integrate(x),x)