Thanks for the detailed remarks, and for pointing out that my statement above about the novelty of our implementation was wrong - CasADi uses the same algorithm as was released in 2013. My apologies. Joel Andersson pointed this out to us and we added a cite to his thesis several months ago.
As for Sundials and FatODE, my understanding was that they used finite differences, forward mode, or differentiating the solver operations for at least some aspect of their sensitivity analysis.
On another topic, I think you might be misunderstanding how we're running our adjoint sensitivity analysis. You say
> 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.
This isn't true - when we do the reverse solve, we get all gradients using a _single_ solve going backwards in time. I'm now realizing that the misunderstanding might be caused by our Fig. 2, which shows that multiple backward solves are necessary when the loss depends on the state at multiple time points.
I agree that our statement about the stability of implicit over explicit was overly broad, thanks for pointing that out. Can you suggest a more accurate statement about the advantages of implicit over explicit methods?
I also agree there is a lot of numerical work to be done in this area, and I'm glad that people more knowledgable than us (such as yourself) are looking at it too!
>As for Sundials and FatODE, my understanding was that they used finite differences, forward mode, or differentiating the solver operations for at least some aspect of their sensitivity analysis.
No, their sensitivity analysis doesn't use finite differences. They perform the sensitivity analysis as described in their documentation. If methods for the Jacobian calculation are not provided, then they utilize finite differences on the Jacobian calculation of course, using the same routine as for the stuff solver. But with a given Jacobian function there's no finite differences or forward mode.
>This isn't true - when we do the reverse solve, we get all gradients using a _single_ solve going backwards in time. I'm now realizing that the misunderstanding might be caused by our Fig. 2, which shows that multiple backward solves are necessary when the loss depends on the state at multiple time points.
No. Of course it's a single solve going backwards in time. I see what my misread was, but I'm surprised you'd attempt to solve the equation backwards like that because it's known to not be stable. Without a reversible integrator (implicit Adams is not reversible), it's a well-known result that the method drifts from the true solution doing a backwards integration, so the values so z(t) needs to computed with forward passes in order to be correct. A good test equation for this is probably the Lornez equation with standard parameters over a time like [0,300]. The backwards pass will diverge and not necessarily be on the same butterfly wing as the forwards pass. The Julia code using CVODE is shown here ( https://gist.github.com/ChrisRackauckas/fef4ae7778320530d44b... ) and you see that starting from [1.0,1.0,1.0] the result of going backwards is then off by [-17.5445, -14.7706, 39.7985]. So there you go, CVODE's Adams method ends up on a different "wing" of the butterfly when integrated backwards, ending up not even close to the actual initial point (CVODE is the successor to LSODE, both by Alan Hindmarsh, but utilizes constant leading coefficient forms to reduce computations. So not exactly the same as the paper, but very close). Thus to ensure correctness, existing sensitivity analysis packages only get Jacobians of f using data from forward passes. The Neural ODEs may have had a small enough Lyopunov coefficient or a short enough integration that this wasn't an issue, but it is in general something to note. Of course, if you are only doing this on Hamiltonian systems...
>I agree that our statement about the stability of implicit over explicit was overly broad, thanks for pointing that out. Can you suggest a more accurate statement about the advantages of implicit over explicit methods?
Any statement on it is too broad to be useful. Runge-Kutta Chebyshev methods are explicit methods for stiff systems. Implicit Adams is an implicit method for non-stiff systems. And there's many more examples. The stiffness handling also depends on implementation details. Using functional iteration on a BDF method reduces the region of stability, which is why BDF needs to use Newton's method for solving the implicit equation in order to be applicable to stiff ODEs. It's best to just talk about the stability of individual methods and their implementation.
> But with a given Jacobian function there's no finite differences or forward mode.
Right, but instantiating an entire Jacobian is always going to scale at least quadratically with time. The point I was trying to make is that the existing non-adjoint approaches were never going to scale to large systems with millions of parameters. This is the main attraction of reverse-mode, and it appeared to me that this was a major obstacle for fitting large models using existing packages (excepting CasADi).
> I'm surprised you'd attempt to solve the equation backwards like that because it's known to not be stable.
> it's a well-known result that the method drifts from the true solution doing a backwards integration, so the values so z(t) needs to computed with forward passes in order to be correct.
I agree that a purely reverse-mode gradient solve will diverge from the forward trajectory to some degree. But to say the gradients are 'correct' or not seems a bit strange to me. Every numerical solve introduces some degree of error, and I think the most useful discussion to have is the tradeoff between computation cost and numerical error. Re-solving the system forwards is one strategy to reduce error at the cost of computation. Another strategy would be reducing the error tolerance of the reverse solve. There are situations where our strategy might give worse precision wrt the parameter gradients for a given computational budget, but I wouldn't dismiss it out of hand. Especially since it's about as computationally cheap as one could hope for - O(1) memory and similar time cost as the forward solve. Also, it worked for our applications.
> Any statement on it is too broad to be useful.
I appreciate the detailed reply. But is there anything you can say about when to try implicit methods over explicit? What was the motivation for developing implicit methods in the first place?
As for Sundials and FatODE, my understanding was that they used finite differences, forward mode, or differentiating the solver operations for at least some aspect of their sensitivity analysis.
On another topic, I think you might be misunderstanding how we're running our adjoint sensitivity analysis. You say
> 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.
This isn't true - when we do the reverse solve, we get all gradients using a _single_ solve going backwards in time. I'm now realizing that the misunderstanding might be caused by our Fig. 2, which shows that multiple backward solves are necessary when the loss depends on the state at multiple time points.
I agree that our statement about the stability of implicit over explicit was overly broad, thanks for pointing that out. Can you suggest a more accurate statement about the advantages of implicit over explicit methods?
I also agree there is a lot of numerical work to be done in this area, and I'm glad that people more knowledgable than us (such as yourself) are looking at it too!