Analytic Parameter Jacobians

[mjohnson541]

Forward sensitivity speeds are currently held back by the parameter jacobian evaluations, which can be enormously expensive to compute for our problems (lots of parameters) without sparse differentiation. For simple problems the fastest way to do this should be by adding analytic parameter jacobians. This should be easier to program than the variable jacobians as they should be more consistent between reactors.

[ChrisRackauckas]

You shouldn't need to calculate the Jacobian in order to compute the forward sensitivity equations. How were you doing the call?

[mjohnson541]

So I appreciate that we don't need the variable jacobians unless they're very fast since we can instead calculate the VJP, but I thought we still needed to calculate the parameter jacobian (F) in https://diffeq.sciml.ai/stable/extras/sensitivity_math/? F is quite expensive because we typically have 100x more parameters than we have variables.

Currently we're being beaten pretty badly by RMG's internal sensitivity analysis code, which uses analytic Jy and Jp rather than VJP.

For forward sensitivities we currently feed in jacp! calculated by ForwardDiff to the ODEFunction only. ReverseDiff/Tracker would certainly make more sense, but I'm still working on making the internals more compatible with those.

[ChrisRackauckas]

You don't need the full Jacobian since jacobian-vector products are the primitive of forward-mode, so you can do a single dual number f(Dual(x,v)) for J(x)*v. In fact, you can do this with numerical differencing as well, since J*v is just the directional derivative: (f(x+eps*v)-f(x))/eps. So forward sensitivities should never need the Jacobian, and that's a naturally sparse calculation. You don't want to use reverse-mode here since that's asymptotically worse for this case.

But the forward sensitivity implementation in DiffEqSensitivity.jl should already be doing this by default, with the absolute fastest way to do forward-mode though being to just use ForwardDiff over the integrator (demonstrated in https://arxiv.org/abs/1812.01892).

[mjohnson541]

I guess I don't understand how you can formulate the parameter derivatives in the forward sensitivity equations as a jacobian vector product? Don't they have to be computed separately by auto diff or finite differences?

[ChrisRackauckas]

Oooh, parameter derivatives are the parameter Jacobian. My bad, I thought you were talking about the other term. Yes, that term needs a fast Jacobian, probably sparse, and analytical will be faster than sparse even in its best case, though it should be asymptotically similar but this would be a spot to really use an analytical Jacobian if you have one.

[ChrisRackauckas]

Though ForwardDiff over the solver makes this naturally sparse and fast, so I wouldn't use the continuous sensitivity form at all and go discrete.

[mjohnson541]

I've been eyeing JVODE_BDF(), but I think currently differentiating over the solver is not possible if we stick with CVODE_BDF() and I've had trouble finding anything competitive with CVODE_BDF() on these problems.

[ChrisRackauckas]

TRBDF2 and KenCarp4 haven't done well?

[ChrisRackauckas]

I've been eyeing JVODE_BDF(), but

We have a private repo with a few things being developed, but won't be ready until December. That JVODE_BDF will likely be one of the things @YingboMa gets to do a lot of time on in the next few months.

[mjohnson541]

Cool! The test I ran was a couple years ago I'll run some benchmarks on TRBDF2 and KenCarp4!

[ChrisRackauckas]

A couple years ago was bad. We bested CVODE in https://benchmarks.sciml.ai/html/StiffODE/BCR.html only within the last year. That's what I'd point to for a 1,000 ODE performance case.

How many ODEs are you doing and how sparse?

[mjohnson541]

Our systems are usually 100-3000 ODEs. The sparsity depends a lot on the size, it should look like the sparsity in this paper: https://www.sciencedirect.com/science/article/pii/S0010218001003522. Maybe ~90% sparse for large stuff.