[oivey]

Yeah, this is what I meant. It’s virtually impossible for Python-based libraries to catch up to the huge number of differential equation solvers that are available in Julia, and that work on a GPU, and that can be automatically differentiated through, and that all work through a common interface, and that all work around a well defined ecosystem.

More broadly, I’d say Julia has more potential than Python when it comes to automatic differentiation. JAX is really cool, but Julia’s automatic differentiation tools are targeting the whole language without modification and restriction (outside those imposed by calculus). Potentially some really fancy stuff could get embedded in neural nets beyond differential equations and with relatively little effort.

[orbifold]

I find that claim highly dubious, solving differential equations rarely yields to completely generic methods and there are plenty of solver suites in other languages, among them C++ and wrappers for Python. To use scientific machine learning you don't necessarily implement all methods under the sun, just a method appropriate for the domain at hand.

[oivey]

Implementing automatic differentiation through C++ solvers or Python wrappers is going to be exceptionally difficult. If you want a generic toolbox like torch for scientific machine learning you would want every method under the sun.

That’s the moat Julia has here. A full rewrite of the differential equation ecosystem wasn’t necessary in order to support automatic differentiation.

[krastanov]

I do not know of another differential equations framework nearly as polished, performant, and feautere rich as the one in Julia, in any language. I actually know of no other framework that has both modern solvers and autodiff working together. The only other autodiff solvers I know of are ancient RK4-like solvers implemented in things like Tensorflow. It is hard to overstate how monumentally enabling such tools are for old-school science and they currently exist only in Julia.

[orbifold]

Sundials (https://computing.llnl.gov/projects/sundials) is arguably more feature complete, has many optimised implementations scaling to super computers and supports both adjoint and forward sensitivity analysis. Btw. several "research" papers at machine learning conferences on "Neural Differential Equations" presented "discoveries" covered in the documentation of Sundials https://sundials.readthedocs.io/en/latest/cvodes/Mathematics....

[ChrisRackauckas]

Let's talk details then. On the small stiff ODE and non-stiff end SUNDIALS is bested quite easily (for example https://benchmarks.sciml.ai/html/StiffODE/Hires.html), but that's pretty well known so not much to say there. For the larger ODE end, it's missing Runge-Kutta-Chebyschev methods, exponential integrators, etc. which do amazing on a lot of equations, and in direct head to head FBDF has been outperforming CVODE BDF for over a year now on most equations (as demonstrated in the online benchmarks, and now there are some applications papers starting to show that). Where SUNDIALS does do well is for distributed NVector is really good, along with it's heterogeneous compute. Its new multirate methods are good too. But then it's lacking methods for SDEs, DDEs, RODEs, higher precision, probabistic numerics, etc. so it's hard to call it feature filled. It's also lacking a complete symbolic front end like ModelingToolkit for handling higher index DAEs, so while IDA currently outperforms DFBDF head-to-head, the optimal thing to do is still to use a Julia solver just post MTK optimized (for example see https://benchmarks.sciml.ai/html/DAE/ChemicalAkzoNobel.html and for the large equation side see some of the stuff on HVAC models). Sundials has deep customizability in the linear and nonlinear solvers but that's also seen on the Julia side.

And while SUNDIALS does have an adjoint implementation, it does not construct the vjp through reverse mode so direct (naive) usage is really lacking in performance. Yes I know that there is a way to define the vjp function in the docs, but it's buried enough in there that all of the examples I've looked at of SUNDIALS in applications in the wild (like in AMIGO2) don't make use of it so there's a real performance loss seen by end users for many cases. DiffEq of course automates this.

Don't get me wrong, SUNDIALS is really good, still bests the Julia solvers on the distributed case, and there is a small performance advantage to CVODE in the one doc example on the 2D Brusselator after adding iLU preconditioning (https://diffeq.sciml.ai/dev/tutorials/advanced_ode_example/ and the low tolerance case of the Filament benchmark (https://benchmarks.sciml.ai/html/MOLPDE/Filament.html, though the medium tolerance case is bested by ROCK4). For years DifferentialEquations.jl used wrappers to SUNDIALS as a crutch for the cases where SUNDIALS still outperformed it, so nowadays it both has a fairly complete wrapper of SUNDIALS but also has benchmarks of different algorithms outperforming it.

We've been benchmarking against it weekly for years now, and within the last two years we've really seen a reversal so now it's a rare sight to see SUNDIALS perform the best. And there's a lot of momentum too: new parallel solvers, new GPU solvers, etc. are coming out from both crews so there's a fun arm's race going on right now, so we will see how it evolves (I definitely want to wrap those multirate methods ASAP and play with the theory a bit)

> Btw. several "research" papers at machine learning conferences on "Neural Differential Equations" presented "discoveries" covered in the documentation of Sundials https://sundials.readthedocs.io/en/latest/cvodes/Mathematics....

Of course the Julia folks know this because we were the ones mentioning this at conferences, since we would routinely choose to benchmark against SUNDIALS implementations instead of random unstable ML ones.

[krastanov]

Independently on our views on the comparisons of Sundials and DifferentialEquations.jl, I want to wholeheartedly agree with you that modern research in Neural Differential Equations and ML in general is rediscovering a ton of stuff that is obvious of everyone that has done serious numerical work.

[adgjlsfhk1]

while there are c++ solvers, they are missing the flexibility needed to solve many types of problems efficiently. they almost all use finite differencing, which is inaccurate and slow, and they rarely let users customize the linear solve and nonlinear solve algorithms. What the Julia diffeq suite is working towards is complete control over solver options, while picking good defaults when not customized. They let you pick forward and backward mode AD (or a combination of the two), problem specific preconditioning methods, custom root finding algorithms (including gpu based linear solvers), arbitrary precision methods for verification (or if you're just paranoid), interval arithmetic, linear uncertainty propegation, and more.

[orbifold]

As far as I'm aware there is considerable prior art, since solving this problem has been relevant since at least the 70s. Granted the barrier to entry was far higher then, then now, but we are talking about scientific / engineering applications and there the prior investment in specialised tools has always been justified. A contemporary example is http://www.dolfin-adjoint.org/en/latest/, but there are others. Sundials for example also both supports completely custom linear and non-linear solvers in conjunction with adjoint sensitivity analysis.

[adgjlsfhk1]

there is lots of prior art but it's still not ideal. the one you posted, for example only has reverse mode AD which sometimes is good, but forward mode is sometimes asymptotically faster (and has lower overhead, and some numerical stability issues)

[ChrisRackauckas]

Sundials does have forward sensitivities, though the implementation does not SIMD the calculations of the primal with the dual parts like direct AD of the solver does. This makes it a bit worse on forward mode than direct differentiation of a solver (with an optimizing compiler) would do, which was something that took us a few years to really pin down (we had to implement our own to really benchmark it, and got the same phenomena demonstrated in https://arxiv.org/abs/1812.01892). The nice thing about the direct forward mode is that it is easier to opt the extra equations out of the linear solver, something that AD systems won't do by default (but it can be done)