Choosing the best options for a PEtab problem

PEtab.jl provides several gradient and hessian methods that can be used with the ODE solvers in the DifferentialEquations.jl package. You can choose from a variety of options when creating a PEtabODEProblem. If you do not specify any of these options, appropriate options will be selected automatically based on an extensive benchmark study. These default options usually work well for specific problem types. In the following section, we will discuss the main findings of the benchmark study.

Small models ($\leq 20$ parameters and $\leq 15$ ODE:s)

ODE solver: For small stiff models, the Rosenbrock Rodas5P() solver is often the fastest and most accurate option. While Julia bdf-solvers such as QNDF() can also perform well, they are often less reliable and less accurate than Rodas5P(). If the model is "mildly" stiff, composite solvers such as AutoVern7(Rodas5P()) often perform best. Regardless of solver, we recommend using low tolerances (around abstol, reltol = 1e-8, 1e-8) to obtain accurate gradients.

Gradient method: For small models, forward-mode automatic differentiation via ForwardDiff.jl tends to be the best performing option, and is often twice as fast as the forward-sensitivity equations approach in AMICI. Therefore, we recommend using gradient_method=:ForwardDiff.

Hessian method: For small models, it is computationally feasible to compute an accurate full Hessian via ForwardDiff.jl. For most models we benchmarked, using a provided Hessian improved convergence. Therefore, we recommend using hessian_method=:ForwardDiff.

Overall, for a small model, a good setup often is:

petab_problem = PEtabODEProblem(petab_model,
                                ode_solver=ODESolver(Rodas5P(), abstol=1e-8, reltol=1e-8),
                                gradient_method=:ForwardDiff,
                                hessian_method=:ForwardDiff)

Medium-sized models ($\leq 75$ parameters and $\leq 50$ ODE:s)

ODE solver: For medium-sized stiff models, bdf-solvers like QNDF() are often fast enough and accurate. However, they can fail for certain models with many events when low tolerances are used. In such cases, KenCarp4() is a good alternative. Another option is Sundials' CVODE_BDF(), but it is written in C++ and not compatible with forward-mode automatic differentiation. To obtain accurate gradients, we recommend using low tolerances (around abstol, reltol = 1e-8, 1e-8) regardless of solver.

Gradient method: For medium-sized models, when using the Gauss-Newton method to approximate the Hessian, we recommend computing the gradient via the forward sensitivities (gradient_method=:ForwardEquations), where the sensitivities are computed via forward-mode automatic differentiation (sensealg=:ForwardDiff). This way, the sensitivities can be reused when computing the Hessian if the optimizer always computes the gradient first. Otherwise, if a BFGS Hessian-approximation is used, gradient_method=:ForwardDiff often performs best.

Hessian method: For medium-sized models, it is often computationally infeasible to compute an accurate full Hessian via ForwardDiff.jl. Instead, we recommend the Gauss-Newton Hessian approximation, which often performs better than the commonly used (L)-BFGS approximation. Thus, we recommend hessian_method=:GaussNewton.

Overall, when the gradient is always computed before the Hessian in the optimizer, a good setup is often:

petab_problem = PEtabODEProblem(petab_model,
                                ode_solver=ODESolver(QNDF(), abstol=1e-8, reltol=1e-8),
                                gradient_method=:ForwardEquations,
                                hessian_method=:GaussNewton,
                                reuse_sensitivities=true)

Otherwise, a good setup is:

petab_problem = PEtabODEProblem(petab_model,
                                ode_solver=ODESolver(QNDF(), abstol=1e-8, reltol=1e-8),
                                gradient_method=:ForwardDiff,
                                hessian_method=:GaussNewton)

Large models ($\geq 75$ parameters and $\geq 50$ ODE:s)

ODE solver: To efficiently solve large models, we recommend benchmarking different ODE solvers that are suitble for large problems; such as QNDF(), FBDF(), KenCarp4(), and CVODE_BDF(). You can also try providing the ODE solver with a sparse Jacobian (sparse_jacobian::Bool=false) and testing different linear solvers such as CVODE_BDF(linsolve=:KLU). Check out this link for more information on solving large stiff models in Julia.

Gradient method: For large models, the most scalable approach is adjoint sensitivity analysis (gradient_method=:Adjoint). We support InterpolatingAdjoint() and QuadratureAdjoint() from SciMLSensitivity (see their documentation for info), but we recommend InterpolatingAdjoint() because it's more reliable.

Hessian method: For large models, computing the sensitives (Gauss-Newton) or a full hessian is not computationally feasible. Thus, the best option is often to use an L-(BFGS) approximation. BFGS support is built into most available optimizers such as Optim.jl, Ipopt.jl, and Fides.py.

All in all, for a large model, a good setup often is:

petab_problem = PEtabODEProblem(petab_model,
                                ode_solver=ODESolver(CVODE_BDF(), abstol=1e-8, reltol=1e-8),
                                ode_solver_gradient=ODESolver(CVODE_BDF(), abstol=1e-8, reltol=1e-8),
                                gradient_method=:Adjoint,
                                sensealg=InterpolatingAdjoint())