Jacobian for system of ODEs


I may have some code ready somewhere which does what you want… but apart from that it sounds that you are solving the large ODE system for multiple patients, right? If that is the case, then you should by all means use the new map_rect function and parallelize the evaluation of the log-likelihood. It is a bit of a pain to do the packing and unpacking of the data, yes, but in your case you will get almost for sure a linear speedup of your computation time with the number of cores per chain.


Yes, hierarchical model for multiple patients. If there is a way to supply jacobian for ODEs in cmdstan 2.18 I would be grateful. I would like to use map_rect - maybe you could share some example?

In general, to speed up things I was thinking about approximated inference (variational inference) and transforming ODEs to algebraic equations (via orthogonal collocation or some other techniques, perhaps you could suggest). Since stan already has AE solver it would be great to supply jacobian for system of AEs.

In reality I think it would be nice to expose C++ interface for folks who love supplying jacobians for ODEs or AEs. It is OK that it is hard to use - but sometimes one has to bite a bullet.


The manual has an example and there is a StanCon contribution forthcoming from my end which will demonstrate things. For your models you really should use map_rect. These type of problems are ideal for threading or MPI and I would expect that you get linear performance increases with the number of cores up to the point where you have 1 core per patient.

The need to pass in the Jacobian will be much reduced (and BTW, the stiff ODE solver got ~15% faster in 2.18 in examples I tried).


Looking forward for your StanCon contribution while Bob’s logistic regression example is quite self explanatory (http://discourse.mc-stan.org/t/map-rect-in-the-stan-language-mpi-is-coming/3547/2). I still suspect that even with map_rect supplying jacobian for ODEs would reduce time even more (simply because calculating jacobian numerically requires multiple solutions of ODEs). 15% from 5 days is significant but I would like to push the boundaries.


Of course, supplying the Jacobian will help on top of the parallelization… but I can tell you that with brute-force 80 cores I was able to get a 2.7 days runtime down to 1h for an ODE problem. At that point I don’t bother with Jacobians any more… but this was running on a high-performance computing cluster not all of us have available.


Cool, I don’t have such luxury. BTW, if interested I can share model/data to test on the cluster.