Fast derivatives of likelihood functionals for ODE based models using adjoint-state method

Abstract

We consider time series data modeled by ordinary differential equations (ODEs), widespread models in physics, chemistry, biology and science in general. The sensitivity analysis of such dynamical systems usually requires calculation of various derivatives with respect to the model parameters. We employ the adjoint state method (ASM) for efficient computation of the first and the second derivatives of likelihood functionals constrained by ODEs with respect to the parameters of the underlying ODE model. Essentially, the gradient can be computed with a cost (measured by model evaluations) that is independent of the number of the ODE model parameters and the Hessian with a linear cost in the number of the parameters instead of the quadratic one. The sensitivity analysis becomes feasible even if the parametric space is high-dimensional. The main contributions are derivation and rigorous analysis of the ASM in the statistical context, when the discrete data are coupled with the continuous ODE model. Further, we present a highly optimized implementation of the results and its benchmarks on a number of problems. The results are directly applicable in (e.g.) maximum-likelihood estimation or Bayesian sampling of ODE based statistical models, allowing for faster, more stable estimation of parameters of the underlying ODE model.

Appendix: 1 Proofs

1.1 Proof of Theorem 1

First, when compared with Theorem 4.D from Zeidler (1985) we work with $X={\mathbb {R}}^m$ and $P={\mathbb {R}}^p.$ Those are complete normed vector spaces i.e. they are Banach spaces. Second, the initial condition is dependent only on parameter $\varvec{\phi }$, not on any free parameter y as in Theorem 4.D.

Set $J = [-1, 1].$ Let us do the following rescaling: $t = sa,$ ${\varvec{z}}(s) := {{\varvec{u}}}(as) - {{\varvec{u}}}_0(\varvec{\phi })$ for all $s \in J.$ Then (1) is equivalent to

$$\begin{aligned} {\varvec{z}}'(s) - a{\varvec{f}}(as, {\varvec{z}}(s) + {{\varvec{u}}}_0(\varvec{\phi }), \varvec{\phi }) = \varvec{0}\quad \text{ for } \text{ all } s \in J, {\varvec{z}}(0) = \varvec{0}. \end{aligned}$$

This can be written as an operator equation $F({\varvec{z}}, a, \varvec{\phi }) = 0$ with the operator $F:{\varvec{Z}}\times {\varvec{A}}\rightarrow {\varvec{W}}$ and spaces ${\varvec{Z}}= \{{\varvec{z}}\in {\varvec{C}}^1(J,{\mathbb {R}}^m): {\varvec{z}}(0) = \varvec{0}\}$, ${\varvec{W}}= C(J,{\mathbb {R}}^m).$ The space ${\varvec{A}}$ contains all the parameters $(a, \varvec{\phi })$, i.e. ${\varvec{A}}= {\mathbb {R}}\times {\mathbb {R}}^p.$

Set ${\varvec{q}}= (\varvec{0}, 0, \varvec{\phi }).$ Both F and $F_z$ are continuous at ${\varvec{q}}.$ Obviously, $F({\varvec{q}})=\varvec{0}$ and $F_{{\varvec{z}}}({\varvec{q}}){\varvec{z}}= {\varvec{z}}'.$ The crucial observation is that for every ${\varvec{w}}\in {\varvec{W}},$ there exists exactly one ${\varvec{z}}\in {\varvec{Z}}$ with ${\varvec{z}}'={\varvec{w}}$, namely ${\varvec{z}}(s)=\int _0^s {\varvec{w}}(t) \,\mathrm {d}t$. Hence $F_{{\varvec{z}}}({\varvec{q}}):{\varvec{Z}}\rightarrow {\varvec{W}}$ is bijective. The implicit function theorem yields the conclusions [see e.g. Theorem 4.B in Zeidler (1985)]. $\square$

1.2 Proof of Lemma 1

After realizing that $l$ depends on $\varvec{\phi }$ only through ${{\varvec{u}}}$, (12) is formally a direct application of the chain rule ($d_{{{\varvec{u}}}}$ denotes the derivative of the metric with respect to the model state ${{\varvec{u}}}$).

The r.h.s. of (12) is a well-posed finite expression. First, we have assumed that the metric is sufficiently smooth, thus $d_{{{\varvec{u}}}}$ is continuous. Second, ${\varvec{s}}$ is continuous as well owing to Theorem 1. A distribution can be rescaled by any at least continuous function, as here $\{\delta \}$ by $d_{{{\varvec{u}}}}$.

Thus, the first differential on the l.h.s of (12) exists as well. It is moreover continuous, i.e. it is Fréchet.⁵ $\square$

1.3 Proof of Theorem 2

First, let us assume that we have already constructed a unique solution ${\varvec{v}}$ to (15) up to a certain measurement point $t_i$. The adjoint problem is solved backwards in time. Consequently, we will construct its prolongation on $[t_i, t_{i-1})$.

Let ${\varvec{v}}_i^+$ be the ODE solution just before integrating the measurement at time $t_i$, i.e. at time $t_i^+.$ We simply stop the integration at $t_i^+,$ add $d_{{{\varvec{u}}}}({\varvec{y}}_i,{\varvec{g}}(t_i, \varvec{\phi }))$ to ${\varvec{v}}_i^+$ and solve

$$\begin{aligned} \begin{aligned} d_{t}{\varvec{v}}&= -J^t_{{{\varvec{u}}}}({\varvec{f}}){\varvec{v}}, \quad t\in (t_i,t_{i-1}), \\ {\varvec{v}}(t_{i})&= {\varvec{v}}_i^+ + d_{{{\varvec{u}}}}({\varvec{y}}_i,{\varvec{g}}(t_i, \varvec{\phi })). \end{aligned} \end{aligned}$$

(32)

This is a simple linear ODE with a continuous coefficient $J^t_{{{\varvec{u}}}}({\varvec{f}}),$ since $f\in {\varvec{C}}^1$. The classical results yield the global solution on $(t_i,t_{i-1})$ [see e.g. Theorem 5.1 and Theorem 5.2 from Coddington and Levinson (1955)]. This concludes the proof of the existence and uniqueness.

Now, we prove (14). Let us without a loss of generality assume that there are no measurements in times 0 and T. It is a well-known result of theory of distributions (in the sense of functional analysis), that the classical integration by part formula

$$\begin{aligned} \int _0^T d_{t}{\varvec{v}} {\varvec{w}}\ dt = [{\varvec{v}}{\varvec{w}}]_{0}^{T} -\int _0^T {\varvec{v}}d_{t}{\varvec{w}} \ dt \end{aligned}$$

(33)

is valid for ${\varvec{w}}\in {\varvec{C}}^1$ even if the derivative $d_{t}{\varvec{v}}$ exists on [0, T] only in a weak sense, i.e. almost everywhere. Actually, (33) is the definition of the weak derivative of ${\varvec{v}}$ taking only ${\varvec{w}}\in {\varvec{C}}^1_0([0,T]).$ Consequently, since ${\varvec{s}}\in {\varvec{C}}^1([0,T]),$ we can safely proceed as follows

$$\begin{aligned} \begin{aligned} \left<{\varvec{s}},\{\delta \}d_{{{\varvec{u}}}}({\varvec{y}},{\varvec{g}}(t, \varvec{\phi }))\right>&\mathop {=}\limits ^{(15)} \left<{\varvec{s}},d_{t}{\varvec{v}} + J^t_{{{\varvec{u}}}}({\varvec{f}}){\varvec{v}}\right>\\&\mathop {=}\limits ^{(8)} -{\varvec{v}}^t(0)J_{\varvec{\phi }}({{\varvec{u}}}_0){\varvec{h}}-\left( d_{t}{\varvec{s}} - J_{{{\varvec{u}}}}({\varvec{f}}){\varvec{s}},{\varvec{v}}\right) \\&\mathop {=}\limits ^{(8)} -{\varvec{v}}^t(0)J_{\varvec{\phi }}({{\varvec{u}}}_0){\varvec{h}}-\left( J_{\varvec{\phi }}({\varvec{f}}){{\varvec{h}}},{\varvec{v}}\right) . \end{aligned} \end{aligned}$$

(34)

$\square$

1.4 Proof of Lemma 2

The existence and uniqueness of ${\varvec{\varsigma }}$ is a direct results of Theorem 1. Now, (20) is derived as follows:

$$\begin{aligned} \begin{aligned} \left<{\varvec{\varsigma }},\{\delta \}d_{{{\varvec{u}}}}({\varvec{y}},{\varvec{g}}(t, \varvec{\phi }))\right>&\mathop {=}\limits ^{(15)} \left<{\varvec{\varsigma }},d_{t}{\varvec{v}} + J^t_{{{\varvec{u}}}}({\varvec{f}}){\varvec{v}}\right>\\&\mathop {=}\limits ^{(15), (18)} -{({{\varvec{u}}}_0)}_{\varvec{\phi }\varvec{\phi }}{\varvec{h}}_1{\varvec{h}}_2\cdot {\varvec{v}}(0)\\&\quad \quad -\left( d_{t}{\varvec{\varsigma }},{\varvec{v}}\right) +\left( J_{{{\varvec{u}}}}({\varvec{f}}){\varvec{\varsigma }}, {\varvec{v}}\right) . \end{aligned} \end{aligned}$$

(35)

This after substituting for $J_{{{\varvec{u}}}}({\varvec{f}}){\varvec{\varsigma }}$ from (18) directly yields (20). Analogically to the proof of Theorem 2, we needed ${\varvec{\varsigma }}\in {\varvec{C}}^1([0,T])$ to be able to integrate by parts.$\square$