# Speed up solution of a very large system of ODEs

I need to solve many very large systems of first order ODEs, which describe some chemical reactions. The number of variables (in each system) is on the order of $$n \sim 10^5$$. I am using ALGLIB vector ODE solver, which under the hood implements Cash-Karp adaptive ODE solver.

Chemical equations are often considered as stiff and it is recommended to use implicit methods. Unfortunately, implicit methods cannot work with such number of variables at they will require inverting huge matrices. This, in turn, requires on the order of $$\sim n^3$$ operations just to start from and that's before rounding errors are taken into account. Last time when I "played" around with an order smaller matrices ($$n \sim 10^4$$), I needed around 50 digits of precision in order to get reasonably correct results. Such precision is not natively supported by current processors and all together this effectively rules out implicit methods for such number of variables.

From another side, adaptive explicit methods, like the one used in ALGLIB, calculate the errors and decrease the step to a very small number. That often makes the solution time insanely large (like months).

My guess is that this happens when some of the variables start to approach zero. As all variables must be strictly non-negative numbers (they are concentrations of some substances), once some of them get close to zero, then the relative errors immediately become very large and the algorithm then decreases the step. In fact, the algorithm would routinely overshoot zero thus making some variables negative. And after that everything would blow up. So, I treat all variables less than zero as exact zeros, when calculating the derivative. That makes the algorithm stable but still does not delete the fact that it likely decreases the step too much.

The parallelization for a single ODE system does not work well, or, shall I say, does not work at all. This is due to the fact that the ODE solver repeatedly calls a fairly quick function, which calculates the derivative. Such functions are based on just pure math and for such functions (even with the input array of $$n \sim 10^5$$ size) the parallelization overhead just kills the purpose. Rather, the parallelization is achieved by spawning multiple models at the same time where each model runs under a single thread. This can be done spawning multiple threads or multiple external processes (with each thread / process running a single "model"). Both have some benefits and drawbacks and this is irrelevant to the question. The bottom line is that parallelizing a single model kills performance.

The system has an integral of motion: the total number or "atoms" (a sum over all variables multiplied by some weights) in the system must be constant. And since all variables are strictly non-negative, this results in the upper bound for a derivative.

I wonder what can be done to solve such systems of ODEs [much] faster.

• I assume your system is sparse? Are you exploiting this sparsity in an implicit solver? I suggest taking a look at CVODE in SUNDIALS, computing.llnl.gov/projects/sundials – Bill Greene Feb 27 at 23:54
• Does the formulation lend itself to parallelization if any kind? It would of interest to me to know more about how you’re doing the computations. – A rural reader Feb 28 at 0:38
• @KonstantinKonstantinov if explicit methods decrease their time step such that the solution visually seems far too precisely resolved, that is an indicator that your problem is stiff, as is very often the case with chemistry-related models. What you said about implicit solvers seems rather wrong, as implicit solvers are routinely used on very large scale systems (e.g. eactive fluid dynamics). Also, I guess you have a mesh that discretizes a physical space. Then in each cell you have reactions occurring. If that's the case, you can look into operator splitting (continued in next comment). – Laurent90 Feb 28 at 18:18
• @KonstantinKonstantinov Operator splitting allows you to solve the reaction operator on its own with specialized integrators, and the other operators (diffusion, convection...) with other integrators (explicit). That way, you can deal with the stiff chemistry with an implicit method (e.g Radau). Moreover, if you indeed have a mesh, each cell "reacts" independently of its neighbours, therefore you actually have multiple 0D reactors which you can solve sequentially/in parallel. Otherwise, you can solve them together easily, as the Jacobian will be very sparse (bandwith~number of species). – Laurent90 Feb 28 at 18:22
• Also non-negativity of your variables is a very specific problem that has already been discussed in this forum before I believe. Usually, lowering the absolute error tolerance specified to your adaptive integrator allows for a better handling of the near zero values. – Laurent90 Feb 28 at 18:26

If you want a parallel linear algebra library that handles arbitrary precision, I have had good luck with Elemental. That might allow you to use an implicit method.

I use my own fork which replaces MPFR with GMP (about 2x faster, but not all operations work). There is another fork by LLNL, but I do not know how well it works.

Elemental has support for a very large number of operations. There is some documentation here. I also have a copy of the documentation repository, but you have to build it yourself. The original developer stopped working on it a few years ago, so support is basically self-service.

• I am not going to down vote [yet] but the question was about solving very large ODE system and not about linear algebra. So, I am kindly asking to remove this "answer". I will down vote it in a week if it is not removed. Thanks. – Konstantin Konstantinov Mar 4 at 3:32
• @KonstantinKonstantinov: Though shit. What do you lose by just ignoring this answer, and maybe (i) take that as an indication that you could be more compact/precise in your OP, and maybe (ii) appreciate the effort put in this answer even though it might not completely satisfy you? After all, some guy is trying to help you for free... – davidhigh Mar 4 at 6:21
• @KonstantinKonstantinov One strategy you mentioned was to use an implicit method. You thought it was infeasible because It would require inverting a large matrix with high precision arithmetic. This library might make it feasible. – Damascus Steel Mar 4 at 7:19
• Using an arbitrary precision scalar can improve convergence behavior, especially when the solvers are unable to produce correct step directions when using very small numbers. This may be a better approach than intervening and zeroing out a variable explicitly. ALGLIB has an older build with MPFR – Charlie S Mar 4 at 11:37
• @DamascusSteel I am afraid that you did not pay enough attention to the size of the problem. The number of variables that I have is around 1.25*10^5, which makes a square matrix of doubles to take (1.25*10^5)^2 / (1024^3) GB or around 115GB. Double the precision and this will become 230 GB and that's before anything is done with such a matrix. So, nothing matrix related will work for that kind of problem. – Konstantin Konstantinov Mar 5 at 12:10

Following the ideas presented by @BillGreene, I put a small C# NET5 interop to FORTRAN ODE Solver DLSODE: https://github.com/kkkmail/OdePackInterop. The interop includes C# tests and identical F# tests for a chemical-like system with 100,001 variables. The number of variables can be adjusted via command line. The results there are self-explanatory.

The approach that worked was to treat all negative values as exact zeros (hereinafter called non-negativity constraint) when calculating derivative and use "functional" corrector iterator (as it is called in DLSODE). There seems to be no big difference between "Adams" and "Bdf" solver methods, though "Adams" seems a little bit faster.

Another corrector iterator candidate "Chord with diagonal Jacobian" did not conserve the integral of motion in all tests. However, it came close and the integral of motion deviates only in the third digit (in contract to functional solver where it is conserved with more than about 10 digits of precision).

All other corrector iterators in DLSODE solver and all other solvers in that package require full or banded Jacobian and that will increase calculation time dramatically, not to mention the challenges of dealing with row / column swapping in multidimensional arrays when performing an interop between NET and FORTRAN. Subsequently, they were ruled out.

2021-03-27 Update

The things changed when mentioned above solvers were presented with the real chemical problems, which had on the order of 10^5 variables.

The estimated average run time (or actual run time when the solvers did complete) per [some specific] model were as follows (assuming that a model is run in a single thread):

AlgLib / Cash-Carp / Use non-negativity constraint: about 8,000 computer years. The step size became outrageously small and that explained why it stalled.

DLSODE / (Adams or BDF) / Functional / Use non-negativity constraint: about 2.5 computer years. And it would blow up without non-negativity constraint.

DLSODE / BDF / "Chord with diagonal Jacobian"/ Use non-negativity constraint: somewhere between half a day to two days (depending on the model). That was the only one that managed to complete within reasonable time. And it would also blow up without non-negativity constraint.