Yes, we've had this working in Julia for a few years now (with Flux.jl and ReverseDiff.jl tracing-based reverse-mode) and it's surprising that wasn't know because it's been pretty widely publicized (and mentioned to the authors). In fact, we even have a preprint showing that reverse-mode adjoint isn't what you want to use in most scenarios:
https://arxiv.org/abs/1812.01892
We identified that tracing-based reverse-mode adjoint + solving without mutation is not competitive with the other sensitivity analysis methods (this is the method which they have implemented in this repository) when comparing it to methods which are fully compiled and utilize mutation. Since the methods in SciPy don't utilize mutation (and add quite a bit of their own overhead, and don't necessarily work well with the JIT compilation being of context switching, as previously demonstrated in http://juliadiffeq.org/2018/04/30/Jupyter.html ), I wouldn't be surprised that they saw a speedup anyways over the odeint+adjoint method described in the paper, but that's mostly because the method from the paper is highly inefficient.
In terms of efficient implementations, the first ones which are equivalent to a reverse-mode adjoint were actually those of FATODE ( https://dl.acm.org/citation.cfm?id=2048596&dl=ACM&coll=DL ) which implemented quite a few Runge-Kutta and Rosenbrock methods with forward-mode and reverse-mode transforms of the solvers manually written in there (described as "discrete sensitivity analysis") for hooking into Fortran AD implementations. JuliaDiffEq's methods are based along this existing heritage and testing.
But sensitivity analysis goes back a lot further. In their paper though they utilize the adjoint sensitivity analysis, which the most common software one would point to for this is SUNDIALS CVODES (https://computation.llnl.gov/sites/default/files/public/cvs_... ). If you compare page 26 of the SUNDIALS manual to the Neural ODE paper, you'll notice that it's eerily similar because it's almost the same method. And actually, 2.7.1 describes an improvement to the method in the Neural ODE paper by using "checkpointing". In the Neural ODE paper, to do a reverse solve of the adjoint ODE it solve the forward ODE from the beginning time point until the point. Clearly, this is really slow because it requires a lot of forward solves over long intervals. You can instead save a few time points (checkpoints) along the forward solve and use those as the starting point, and this is the method implemented in CVODES and IDAS. In JuliaDiffEq, since we were focused on less memory intensive applications (CVODES did it for PDEs), we did it in a way where you get a continuous solution in one forward solve, and then just use the interpolation for an essentially free calculation for the Jacobian at a given time point. This of course goes from many ODE solves to a single ODE solve, at the cost of using more memory, and we proposed implementing checkpointing in a GSoC. But anyways, the derivation of the paper for the adjoint method wasn't new and it was actually already improved by 2005 in the TOMS paper on CVODES ( https://computation.llnl.gov/projects/sundials/toms_cvodes_w... ), with other methods for improving runtime done as well.
In a previous Twitter thread, the senior author also claimed that using the autodifferentiation to do the vector-Jacobian products was a first. For reference, what I mean is that there is a vjp in the adjoint ODE and this can be solved without explicitly building the Jacobian and seeding the backsolve of the derivative function f appropriately. I will admit that in JuliaDiffEq we weren't doing this, but we have since improved our implementation to make use of this trick (and the advantage of this method is shown in the Arxiv paper linked above). That said, the Neural ODE paper isn't the first group to do this either, with earlier work actually being traced back to CasADi (https://link.springer.com/chapter/10.1007/978-3-642-30023-3_...), and of course the author of CasADi was confused as well in the same Twitter thread about the novelty claim here.
And this is just the tip of the iceberg for schemes utilizing adjoint methods. There's entire pharmacometric software (NONMEM) and engineering software (Modelica) built back in the 80's for building the adjoint equations in ways to get these sensitivities for parameter estimation of methods.
So I honestly cannot find a single part of the adjoint method either as mentioned or implemented where this PyTorch implementation is the first. I think a lot of this comes from a lack of technical expertise in the area of numerical differential equations and just an honest mistake though. There are other glaring examples of this in the paper as well. For example, in the paper they mention:
>To solve ODE initial value problems numerically, we use the implicit Adams method implemented in LSODE and VODE and interfaced through the scipy.integrate package. Being an implicit method, it has better guarantees than explicit methods such as Runge-Kutta but requires solving a nonlinear optimization problem at every step.
However, it's well known that not all implicit methods have better stability guarantees than explicit methods. Specifically, the implicit Adams methods they mention are Adams-Bashforth-Moulton methods are not unconditionally stable and have stability similar to explicit methods (you can find this in a lot of numerical analysis course lecture notes, like this one: http://www.math.vt.edu/people/embree/math5466/lecture37_38.p... . The best resource on this is probably Hairer Solving Ordinary Differential Equations I: Non-stiff Problems). Because of this lack of stability, LSODE with Adams coefficients is only recommended non-stiff equations and common stiff test examples like the ROBER will cause it to fail. Instead, the backwards differentiation formula (BDF) is the similar method which is used for stiff equations (there are lots of sources on this. For example, ode113 is the Adams and ode15s is the BDF, and go look at MATLAB's documentation for how they do the recommendations. Another example is CVODE which implements both the Adams and BDF methods, and only recommends the Adams method for non-stiff equations). This also means that the current set of PyTorch differential equation solvers is only applicable to (some) non-stiff ODEs BTW.
I do want to end this critique though by saying that the paper itself describes a very new and novel application of differential equation solvers as part of a neural net in a way that seems very beneficial. So even though what was demonstrated is not technically sophisticated compared to current software and techniques, their application is a very interesting idea and a new research direction which is deserving of a best paper. I don't think we should lose track that the idea of an RNN as an Euler discretization and thus building the continuous one to get better training times is a very neat solution to a machine learning problem and I very much respect the authors for this advance.
I think that what we can learn from this case is that numerical differential equation and machine learning communities need to come together to do this in a way that pushes the research and the software forward. We are starting to do this in Julia, and that's why in our Arxiv paper you can see the JuliaDiffEq crew has teamed up with the machine learning and AD people (Jerret Revels and Mike Innes) to make sure we can get the full adaptive + events + stiffness handling etc. solves working with AD, and get these utilized in new applications for machine learning. With the help of one of the paper's authors we hope to put out a package in January for implementing said neural ODEs in Julia. From there we can all probe the method our own ways. Since my background is in numerical differential equation solvers, I plan to look at this to see what kinds of solvers are optimal in this application (I think some EPIRK type methods might be the right deal). The ML folks will likely want to actually use it as a model to train and predict. The AD folk see this as a good test case for source-to-source transformations, since building a tape a la Flux or PyTorch can be pretty inefficient for large nonlinear operations like an ODE solver (so we plan on seeing what happens with Zygote.jl on this. Right now it errors, but we'll see!).
I see a very nice future merging differential equations and machine learning in various ways, and this paper is a fantastic start in that direction. Let's just make sure we cite the relevant existing work from both communities.