4 ms·

> I don't think this really backs up your statement though. I'm not sure which statement you're reading, but everything you said so far is in agreement with my

> I don't think this really backs up your statement though. I'm not sure which statement you're reading, but everything you said so far is in agreement with my statement: >> Intel has an ODE solver that seems unbelievably good at solving certain kinds of differential equations. I'm surprised/disappointed it wasn't even mentioned. I would've loved to see an evaluation of it against MATLAB (or other solvers) more than any other pair in the list. Note that I have specifically stated that this method seems better on (at least) certain kinds of equations, and that, therefore, I would have loved to see a comparison of it with everything else. In particular, note that I never claimed the method Intel uses is better, at least (especially?) not in general. If I knew that to be the case, I wouldn't have even had a reason to hope for a comparison! The entire point of my comment was that I would be interested in an evaluation of its performance, and I'm glad to see that you feel similarly. :)

StefanKarpinski 9y ago
Gotta say, "long-abandoned implementation of some unclear method that was considerably faster than Matlab's known-to-be-slow implementation on some problems" does not sound particularly interesting. I suspect that Chris is only interested at all because he is spectacularly thorough (as evinced his comparisons with other ODE solvers and his code) and the method might not yet be included in DifferentialEquations.jl. There's little reason to believe that Intel was onto something "magical" here; it seems that they just had a good implementation of a decent method.
ChrisRackauckas 9y ago
Finally got on campus to be able to read the paper in detail. It's an interesting method but it doesn't really seem to have a clear "best use case". It is a Rosenbrock-W (or ROW) type of method, which means that it's a Rosenbrock method where the order of accuracy is not dependent on having an "exact Jacobian". It then uses this fact to allow for "freezing", i.e. storing the same Jacobian across multiple steps. The reason this is nice is because then you can store the factorization and re-use it. (m,k) is with m being the number of linear solves, and k being the number of RHS function evaluations. The paper creates a bunch, and what Intel used was a (5,2) method. However, the (5,2) method, in order to have an embedded method for adaptivity, has its coefficients (for the Rosenbrock tableau) almost completely determined. What this means is that the paper is not able to "optimize" the method for low error. Normally what's done is that a method is written in terms of tableaus, and the numbers for the tableau are chosen such that the principle terms in the error, the coefficients of the leading error terms, are minimized. Sufficiently good RK and Rosenbrock methods can get like 10^-8, meaning that for large enough dt they almost act like they are an order higher. Of course, the large amount of assumptions to get the given method means that this couldn't be done. But does the algorithm have a place? This is the issue... it doesn't necessarily have a good fit anywhere. Standard Rosenbrock methods require a new Jacobian each timestep, but that's fine if factorization isn't sufficiently expensive when is true for smaller ODE systems. Also if the problem is "very nonlinear", i.e. not semi-stiff but solidly stiff, Jacobian re-usage doesn't help anyways. So standard coefficient-optimized Rosenbrock methods do really well when the equations are very stiff and the number of ODEs is small enough. (coefficient-optimized) Modern (E)SDIRK methods naturally allow for Jacobian re-usage but are otherwise similar to the Rosenbrock methods, making them do well in the small number of ODEs but semi-stiff case. The issue with this method is that it tries to sit in the middle here but it's not really optimized for either case. But the issue really goes back to the 2 in the (5,2). Early on in ODE literature, people believed that reducing the number of steps was always a good idea. You can do 5th order RK methods in 7 steps (thus 7 f evaluations, where u'=f with f being the user given ODE function), so why would you ever do 8? However, in the early 90's when we finally got our hands on PCs, people found out that the assumption to keep stages minimized doesn't empirically do well. Instead, what you can do is have more stages, but then use that to design the method such that the local truncation error is orders of magnitude lower. The result is that you can take much larger steps in order to solve the equation but get the same error, making everything quicker. To compensate for the fact that you are now taking large steps you have to have "dense output", i.e. a continuous extension of the method which can use the internal stage data to turn the solution into a continuous function sol(t) which doesn't require any evaluations of f. While this is more cumbersome in to implement, these "higher stage with dense output" methods, the engineering is very much worth it and if you take a look at pretty much any production-quality method these days you'll see they've ended up with these Rosenbrock, (E)SDIRK, or FIRK methods... except for multistep methods. Multistep methods though use history to achieve high order with just a single function evaluation in its implicit function per step. This is thus what's best when f is expensive, and is another reason why stage reduction isn't necessarily helpful because then you're ending up in an area where other methods are already pretty optimal. The classic is the BDF method. It can do the Jacobian re-usage as well. While these aren't L-stable above order 2, these can be adaptive order and do really well in mildly stiff situations. So that said, there really isn't a good niche for these types of methods. I have other low-stage Rosenbrock tableaus implemented and it comes out pretty much as I say here: more error for the same amount of time as others. Rosenbrock-W (ROW) methods are not important anymore because with autodifferentiation we don't need to do numerical differentiation and thus we do get exact Jacobians easily, which negates the purpose of the ROW method (the reason why they were invented because of numerical errors with numerical differentiation causing order loss in Rosenbrock methods. The tradeoff is a lot more assumptions and a lot less coefficient optimization). So the method Intel chose, extrapolating from other low-stage Rosenbrock methods and its other properties, would do best in a cases where the equation is not too nonlinear (you want to make use of the Jacobian freezing, otherwise there's no point since the standard Rosenbrock methods would simply get less error for the same steps), but also it's very stiff (otherwise you're in the regime where BDF dominates), but also not too stiff where you'd want to use the IMEX properties of the modern (E)SDIRK methods, and you'd also want f sufficiently costly that the higher-stage (E)SDIRK methods don't do better. This puts it into a very very small niche, which is probably why this method was lost to the literature and not really adopted. But then again, since it does sit right in the middle of a lot of things, if you're to pick one method to implement, this is probably a pretty good choice. So while Intel probably did a very good implementation of it, I wouldn't see a reason to go looking for it somewhere else. If this method did well, then either a Rosenbrock, (E)SDIRK, or BDF method will probably do better depending on what regime it more falls into.
- throwaway613834 9y ago
  Wow! Thank you for the comprehensive review. Leaves me with nothing else to say but thanks :)

StefanKarpinski 9y ago

Gotta say, "long-abandoned implementation of some unclear method that was considerably faster than Matlab's known-to-be-slow implementation on some problems" does not sound particularly interesting. I suspect that Chris is only interested at all because he is spectacularly thorough (as evinced his comparisons with other ODE solvers and his code) and the method might not yet be included in DifferentialEquations.jl. There's little reason to believe that Intel was onto something "magical" here; it seems that they just had a good implementation of a decent method.

ChrisRackauckas 9y ago

Finally got on campus to be able to read the paper in detail. It's an interesting method but it doesn't really seem to have a clear "best use case". It is a Rosenbrock-W (or ROW) type of method, which means that it's a Rosenbrock method where the order of accuracy is not dependent on having an "exact Jacobian". It then uses this fact to allow for "freezing", i.e. storing the same Jacobian across multiple steps. The reason this is nice is because then you can store the factorization and re-use it. (m,k) is with m being the number of linear solves, and k being the number of RHS function evaluations. The paper creates a bunch, and what Intel used was a (5,2) method. However, the (5,2) method, in order to have an embedded method for adaptivity, has its coefficients (for the Rosenbrock tableau) almost completely determined. What this means is that the paper is not able to "optimize" the method for low error. Normally what's done is that a method is written in terms of tableaus, and the numbers for the tableau are chosen such that the principle terms in the error, the coefficients of the leading error terms, are minimized. Sufficiently good RK and Rosenbrock methods can get like 10^-8, meaning that for large enough dt they almost act like they are an order higher. Of course, the large amount of assumptions to get the given method means that this couldn't be done. But does the algorithm have a place? This is the issue... it doesn't necessarily have a good fit anywhere. Standard Rosenbrock methods require a new Jacobian each timestep, but that's fine if factorization isn't sufficiently expensive when is true for smaller ODE systems. Also if the problem is "very nonlinear", i.e. not semi-stiff but solidly stiff, Jacobian re-usage doesn't help anyways. So standard coefficient-optimized Rosenbrock methods do really well when the equations are very stiff and the number of ODEs is small enough. (coefficient-optimized) Modern (E)SDIRK methods naturally allow for Jacobian re-usage but are otherwise similar to the Rosenbrock methods, making them do well in the small number of ODEs but semi-stiff case. The issue with this method is that it tries to sit in the middle here but it's not really optimized for either case. But the issue really goes back to the 2 in the (5,2). Early on in ODE literature, people believed that reducing the number of steps was always a good idea. You can do 5th order RK methods in 7 steps (thus 7 f evaluations, where u'=f with f being the user given ODE function), so why would you ever do 8? However, in the early 90's when we finally got our hands on PCs, people found out that the assumption to keep stages minimized doesn't empirically do well. Instead, what you can do is have more stages, but then use that to design the method such that the local truncation error is orders of magnitude lower. The result is that you can take much larger steps in order to solve the equation but get the same error, making everything quicker. To compensate for the fact that you are now taking large steps you have to have "dense output", i.e. a continuous extension of the method which can use the internal stage data to turn the solution into a continuous function sol(t) which doesn't require any evaluations of f. While this is more cumbersome in to implement, these "higher stage with dense output" methods, the engineering is very much worth it and if you take a look at pretty much any production-quality method these days you'll see they've ended up with these Rosenbrock, (E)SDIRK, or FIRK methods... except for multistep methods. Multistep methods though use history to achieve high order with just a single function evaluation in its implicit function per step. This is thus what's best when f is expensive, and is another reason why stage reduction isn't necessarily helpful because then you're ending up in an area where other methods are already pretty optimal. The classic is the BDF method. It can do the Jacobian re-usage as well. While these aren't L-stable above order 2, these can be adaptive order and do really well in mildly stiff situations. So that said, there really isn't a good niche for these types of methods. I have other low-stage Rosenbrock tableaus implemented and it comes out pretty much as I say here: more error for the same amount of time as others. Rosenbrock-W (ROW) methods are not important anymore because with autodifferentiation we don't need to do numerical differentiation and thus we do get exact Jacobians easily, which negates the purpose of the ROW method (the reason why they were invented because of numerical errors with numerical differentiation causing order loss in Rosenbrock methods. The tradeoff is a lot more assumptions and a lot less coefficient optimization). So the method Intel chose, extrapolating from other low-stage Rosenbrock methods and its other properties, would do best in a cases where the equation is not too nonlinear (you want to make use of the Jacobian freezing, otherwise there's no point since the standard Rosenbrock methods would simply get less error for the same steps), but also it's very stiff (otherwise you're in the regime where BDF dominates), but also not too stiff where you'd want to use the IMEX properties of the modern (E)SDIRK methods, and you'd also want f sufficiently costly that the higher-stage (E)SDIRK methods don't do better. This puts it into a very very small niche, which is probably why this method was lost to the literature and not really adopted. But then again, since it does sit right in the middle of a lot of things, if you're to pick one method to implement, this is probably a pretty good choice. So while Intel probably did a very good implementation of it, I wouldn't see a reason to go looking for it somewhere else. If this method did well, then either a Rosenbrock, (E)SDIRK, or BDF method will probably do better depending on what regime it more falls into.

throwaway613834 9y ago
Wow! Thank you for the comprehensive review. Leaves me with nothing else to say but thanks :)