Parameter estimation in ODEs. Modelling and computational issues

==Parameter estimation in ODEs. Modelling and computational issues==

'''Abstract.''' In this work we discuss a variational approach for the determination of the parameters of systems of ordinary differential equations (ODE). We construct a model for fitting observed noisy data into the given dynamical system. Also we explain in detail the advantage of using the adjoint equation method to compute the derivatives or gradients, which are needed for the application of gradient methods and quasi-Newton algorithms to find the minimum of the cost function. In particular, we consider two classic iterative algorithms: the conjugate gradient (CG) algorithm and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. For educational purposes we try to explain several numerical and computational issues with some detail and illustrate them with the Susceptible, Exposed, Infected, Recovered, and Deceased (SEIRD) epidemiological model.

'''Keywords.''' Inverse problem; noisy data; parameter determination; adjoint method; variational approach.

==1 Introduction==

Systems of ordinary (and partial) differential equations are an important tool to model the physical state of a real phenomenon that arise in many areas of applied sciences and engineering. Predicting the future behaviour or allowing control of these processes requires not only accurately describing the system but also finding or improving its parameters. Estimation of unknown or inaccurate parameters in turn requires fitting partially observed noisy data or experimental measurements to the model. More generally, in the statistics and machine learning literature various methods have been employed to fit differential equations to data, from maximum likelihood approaches, <span id='citeF-13'></span>[[#cite-13|[13]]] to Bayesian sampling algorithms, <span id='citeF-6'></span>[[#cite-6|[6]]] or traditional deterministic approaches. Thus, parameter estimation needs efficient ODE (forward) solvers, optimization routines, statistical and possible stochastic procedures. There are several stochastic, deterministic and hybrid optimization routines <span id='citeF-2'></span>[[#cite-2|[2]]]. Contrary to stochastic algorithms, deterministic ones are computationally efficient but they tend to converge to local minima. However, they are the departure point in many applications and in the design of better and efficient procedures. For these methods the gradient is combined with information from line searches or methods involving a Newton, quasi-Newton (low-rank) or Fisher information based curvature estimators to update model parameters, <span id='citeF-9'></span>[[#cite-9|[9]]]. The main computational bottleneck in these algorithms is the computation of the gradient (or the curvature) of the parametric cost function. Then, efficient methods to evaluate gradients or for parametric sensitivity analysis of differential&#8211;algebraic models is important, no only for the determination of parameters, but also in other areas of application like model simplification, data assimilation, optimal control, process sensitivity, uncertainty analysis, and experimental design, among others, for a wide range of scientific and engineering problems.

In this work we concentrate in deterministic optimization procedures, mainly those for quadratic non-linear programming and based on gradient methods and quasi-Newton algorithms. Our purpose is to explain in some detail modelling, algorithmic and computational issues to a wide, and possible non-expert, audience. In Section 2 we introduce the quadratic non-linear model that incorporates noisy data into the cost function and it is adapted to the SEIRD epidemiological deterministic model, which is employed to illustrate the algorithms and their related computational issues. Section 3 is devoted to explaining the adjoint equation method for computing gradients of the given cost function and its advantages over other methods. In section 4 we describe the CG and BFGS optimization algorithms and discuss important computational issues. Numerical results are shown in Section 5 and finally, in Section 6, we give some conclusions and perspectives for future work in order to improve the model and overcome some drawbacks and algorithmic problems.

==2 The quadratic model==

Let <math display="inline">\mathbf{x}(t) \in \mathbb{R}^d</math> be a state variable at time <math display="inline">t \in I = [t_0,t_f ]</math> of a continuous time ordinary differential equation (ODE) satisfying the following initial value problem:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\mathbf{x}}(t) = \mathbf{f}(t,\mathbf{x}(t), \boldsymbol{\theta }),  \qquad \mathbf{x}(t_0) = \mathbf{x}_0,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where the upper dot on the left hand side denotes derivation with respect to time. The vector function <math display="inline">\mathbf{f}: \mathbb{R}\times \mathbb{R}^d\times \mathbb{R}^{np} \mapsto \mathbb{R}^d</math> depends on the parameter vector <math display="inline">\boldsymbol{\theta }\in \mathbb{R}^{np}</math> with <math display="inline">np</math> the number of parameters. The solution at time <math display="inline">t</math> of this problem is denoted by <math display="inline">\mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta })</math> when its dependence of <math display="inline">\mathbf{x}_0</math> and <math display="inline">\boldsymbol{\theta }</math> is made explicit. Frequently, we use the short notation <math display="inline">\mathbf{x}(t)</math> for simplicity, like in equation ([[#eq-1|1]]) above. A set of vector measurements at times <math display="inline">t_0 \le t_1 < \ldots < t_m \le t_f</math> in <math display="inline">I</math> are available:

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}_i = \mathbf{x}(t_i;\mathbf{x}_0,\boldsymbol{\theta }) + \boldsymbol{\epsilon }_i, \quad i = 1,\ldots ,m  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

where <math display="inline">\boldsymbol{\epsilon }_i \in \mathbb{R}^d</math> are independent random vectors and represent measurement errors with a multivariate Gaussian distribution having zero mean and vector variances <math display="inline">\boldsymbol{\sigma }_i^2 \in \mathbb{R}^d</math>. In many problems not all components of the state <math display="inline">\mathbf{x}(t)</math> are observable, so this vector variable is decomposed in observable variables <math display="inline">\overline{\mathbf{x}}</math> and non-observable variables <math display="inline">\underline{\mathbf{x}}</math>, which can be regarded as orthogonal projections of <math display="inline">\mathbf{x}</math> over the coordinates of these observable and non-observable variables. A model to estimate the unknown parameter vector <math display="inline">\boldsymbol{\theta }</math> and the initial conditions <math display="inline">\mathbf{x}_0</math>, from the given measurements, relays on the minimization of the least squares objective function

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\ell (\mathbf{x}_0,\boldsymbol{\theta }) = \frac{1}{2} \left\Vert \frac{\mathbf{x}_0 - \mathbf{s}_0}{\boldsymbol{\sigma }_0} \right\Vert ^2_{d} + \frac{1}{2}\sum _{i=1}^{m}  \left\Vert \frac{\overline{\mathbf{x}}(t_i; \mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i}{\overline{\boldsymbol{\sigma }}_i} \right\Vert ^2_{no},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where the state variable <math display="inline">\mathbf{x}(t)</math> is subject to satisfy the ODE ([[#eq-1|1]]). The first norm is the euclidean norm in <math display="inline">\mathbb{R}^d</math> and the norms into the sum are the euclidean norms in <math display="inline">\mathbb{R}^{no}</math>, <math display="inline">no =</math> number of observable variables. The fixed vector <math display="inline">\mathbf{s}_0</math> denotes an experimental measurement of the initial conditions.

Remark 1: We assume that all components of the state variable are observable at the initial time <math display="inline">t_0</math>, otherwise the first term in ([[#eq-3|3]]) can be modified accordingly.  The most used norms, in the construction of objective functions <math display="inline">\ell (\mathbf{x}_0,\boldsymbol{\theta })</math>, are those of the form <math display="inline">|| \mathbf{x}||_p = (\sum _{j=1}^d |x_j|^p)^{1/p}</math> with <math display="inline">p=1</math>, <math display="inline">p=2</math>, <math display="inline">p = \infty </math>. Choosing the most appropriate norm depends on the particular application and of the properties of the state variable. As mentioned before we use the usual euclidean norm, i.e. <math display="inline">p=2</math>.

Remark 2: Observe that the quantities in ([[#eq-3|3]]) are vectors, so the quotients are computed component-wise. In general, if <math display="inline">\Sigma _i</math> is the covariance matrix at experimental time <math display="inline">t_i</math>, then

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> ||(\overline{\mathbf{x}}(t_i; \mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i) / (\overline{\boldsymbol{\sigma }}_i)||^2_{no} =  (\overline{\mathbf{x}}(t_i; \mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i)^T W_i (\overline{\mathbf{x}}(t_i; \mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i) </math>
|}
|}

where <math display="inline">W_i = \Sigma _i^{-1}</math> is the precision matrix. If the symmetric matrix <math display="inline">W</math> is positive definite, then it defines a norm in the corresponding subspace of the observable variables, as <math display="inline">||\overline{\mathbf{x}}||_{_W} = \overline{\mathbf{x}}^T W_i\,\overline{\mathbf{x}}</math>. Of course, it is possible to define the cost function ([[#eq-3|3]]) taking <math display="inline">\overline{\boldsymbol{\sigma }}_i</math> as the unitary vector for every <math display="inline">i</math>, however the minimization process using iterative gradient methods converge faster when using the weights <math display="inline">\overline{\boldsymbol{\sigma }}_i</math>. The attraction ball around the global minimum is also larger in this case.

<span id='theorem-SEIRD'></span>Example 1: To illustrate the previous concepts and notation we consider the SEIRD epidemiological deterministic model that describe the dynamics of the propagation of an infections disease, like COVID-19, <span id='citeF-12'></span>[[#cite-12|[12]]]. If we have a closed constant population of <math display="inline">N</math> individuals, at each time <math display="inline">t</math> a compartmental model separates this population in five segments: susceptible, exposed, infectious, recovered and dead, denoted by <math display="inline">S(t)</math>, <math display="inline">E(t)</math>, <math display="inline">I(t)</math>, <math display="inline">R(t)</math>, <math display="inline">D(t)</math>, respectively. The system of equations they satisfy is given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{dS}{dt} = -\frac{\alpha }{N} S\,I </math>
|-
| style="text-align: center;" | <math> \frac{dE}{dt} = \frac{\alpha }{N} S\,I  - \beta \,E </math>
|-
| style="text-align: center;" | <math> \frac{dI}{dt} =  \beta \,E - \frac{1}{T_I}\, I  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|-
| style="text-align: center;" | <math> \frac{dR}{dt} = \frac{1-f}{T_I}\, I </math>
|-
| style="text-align: center;" | <math> \frac{dD}{dt} = \frac{f}{T_I}\, I   </math>
|}
|}

which is complemented with appropriate initial conditions, <math display="inline">S(t_0)</math>, <math display="inline">E(t_0)</math>, <math display="inline">I(t_0)</math>, <math display="inline">R(t_0)</math>, <math display="inline">D(t_0)</math>. In this system <math display="inline">\alpha </math> is the infection rate, <math display="inline">\beta </math> is the incubation rate, <math display="inline">T_I</math> is the average infectious period and <math display="inline">f</math> is the fraction of individuals who die. Here the dimension is <math display="inline">d = 5</math> and the state variable is the vector <math display="inline">\mathbf{x}(t) = (S(t),E(t),I(t),R(t),D(t))^T</math>, the number of parameters is <math display="inline">np = 4</math> and <math display="inline">\boldsymbol{\theta }= (\alpha ,\beta ,T_I, f)^T \sim (\alpha ,\beta ,\gamma ,\mu )^T</math> with <math display="inline">\gamma = 1/T_I</math> and <math display="inline">\mu = f/T_1</math>. Figure [[#img-1|1]] shows the solution <math display="inline">\mathbf{x}(t)</math>, obtained with the standard RK4 solver, of  ([[#eq-1|1]])&#8211;([[#eq-5|5]]) in the time interval <math display="inline">[t_0,t_f] = [40,80]</math> (days), with exact parameters <math display="inline">\boldsymbol{\theta }=(\alpha ,\beta ,\gamma ,\mu )^T = (1,1/7,1/5,1/70)^T</math>, initial condition <math display="inline">\mathbf{x}_0 = (91647, 4853, 1755, 1620, 125)^T</math> and total population <math display="inline">N = 10^5</math>. The points are experimental measurements at times <math display="inline">t_i = 40 + i</math>, <math display="inline">1\le i\le 13</math> (<math display="inline">m=13</math>). Of course, not always all variables are observable, as shown in this figure. Most frequently the observed variables reported by the medical or government agencies are those people in the population that have been infected, recovered and died, i.e. <math display="inline">\overline{\mathbf{x}_i} = (I_i,R_i,D_i)</math> at times <math display="inline">t_i</math>, so that the non-observable variables are <math display="inline">\underline{\mathbf{x}}_i = (S_i,E_i)</math>.  <div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-exact_noisy_dat.png|300px|The solution of the SEIRD model and synthetic measurements with white noise.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' The solution of the SEIRD model and synthetic measurements with white noise.
|}
Assuming that the total population is constant or its change is negligible during the time period, then above system must be complemented by the conservation equation

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>S(t) + E(t) + I(t) + R(t) + D(t) = N, \quad \forall \, t\in [t_0,t_f],  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

so ([[#theorem-SEIRD|1]])-([[#eq-5|5]]) describe an differential-algebraic system.  Accordingly, we may modify the optimization model ([[#eq-6|6]]), adding a penalized term, obtaining the extended model:

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L(\mathbf{x}_0,\boldsymbol{\theta }) = \frac{1}{2} \left\Vert \frac{\mathbf{x}_0 - \mathbf{s}_0}{\boldsymbol{\sigma }_0} \right\Vert ^2_5 + \frac{k}{2}(\mathbf{x}_0\cdot \mathbf{1} - N)^2 + \frac{1}{2}\sum _{i=1}^{m}  \left\Vert \frac{\overline{\mathbf{x}}(t_i; \mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i}{\overline{\boldsymbol{\sigma }}_i} \right\Vert ^2_3,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

The penalty parameter <math display="inline">k</math> may be proportional to <math display="inline">N</math>.  The added penalized term in ([[#eq-6|6]]) helps stabilizing the optimization process to estimate the initial conditions. The constant vector <math display="inline">\mathbf{1} \in \mathbb{R}^d</math>  has components all equal to one, so the scalar product <math display="inline">\mathbf{x}_0\cdot \mathbf{1}</math> is equal to the sum of the components of <math display="inline">\mathbf{x}_0</math>.

==3 Variational approach==

Our goal is to find the initial conditions <math display="inline">\mathbf{x}_0 \in \mathbb{R}^d</math> and the parameters <math display="inline">\boldsymbol{\theta }\in \mathbb{R}^{np}</math> that minimize the cost function ([[#eq-6|6]]) subject to ([[#theorem-SEIRD|1]]) and ([[#eq-5|5]]), given that we have a set of noisy experimental measurements ([[#eq-2|2]]) at different times. Since the model is deterministic and all the variables and functions are both continuous and smooth, we can employ gradient descent methods or quasi-Newton algorithms and its variants. The most expensive and delicate task, when applying these methods, is the calculation of the gradient or Jacobians of the cost function at each iteration. Some options are available to compute these quantities as shown in <span id='citeF-3'></span>[[#cite-3|[3]]], <span id='citeF-14'></span>[[#cite-14|[14]]], <span id='citeF-4'></span>[[#cite-4|[4]]], and references therein. Here, we are interested on two approaches which are based on variational calculus.

Let <math display="inline">\mathbf{z}= (x_0,\boldsymbol{\theta })^T \in \mathbb{R}^{d+np}</math> be the vector which contain the unknown initial conditions and the unknown  vector of parameters of the dynamical system, then to first order

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla L(\mathbf{z}) \cdot \delta \mathbf{z}\approx L(\mathbf{z}+ \delta \mathbf{z}) - L(\mathbf{z})  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where <math display="inline">\delta \mathbf{z}= (\delta \mathbf{x}_0, \delta \boldsymbol{\theta })^T</math> is an small increment of <math display="inline">\mathbf{z}</math>. If we want to be more specific we can write

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla _{\mathbf{x}_0}\, L(\mathbf{x}_0,\boldsymbol{\theta })\cdot \delta \mathbf{x}_0 + \nabla _\theta \, L(\mathbf{x}_0,\boldsymbol{\theta })\cdot \delta \boldsymbol{\theta }\approx L(\mathbf{x}_0 + \delta \mathbf{x}_0, \boldsymbol{\theta }+ \delta \boldsymbol{\theta }) - L(\mathbf{x}_0,\boldsymbol{\theta }), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

where <math display="inline">\nabla _{\mathbf{x}_0}</math> and <math display="inline">\nabla _\theta </math> denote the gradients with respect <math display="inline">\mathbf{x}_0</math> and <math display="inline">\boldsymbol{\theta }</math>, receptively and the dot is used to denote the corresponding scalar products. Evaluating directly <math display="inline">L(\mathbf{x}_0 + \delta \mathbf{x}_0, \boldsymbol{\theta }+ \delta \boldsymbol{\theta })</math> from ([[#eq-6|6]]), and simplifying, we obtain

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla _{x_0}\, L(\mathbf{x}_0,\boldsymbol{\theta }) \cdot \delta \mathbf{x}_0 + \nabla _\theta \, L(\mathbf{x}_0,\boldsymbol{\theta })\cdot \delta \boldsymbol{\theta }\approx \frac{\mathbf{x}_0 - \mathbf{s}_0}{\boldsymbol{\sigma }_0^2}\, \cdot \delta  \mathbf{x}_0 + k(\mathbf{x}_0\cdot \mathbf{1} - N)\,\mathbf{1}\cdot \delta \mathbf{x}_0  </math>
|-
| style="text-align: center;" | <math>  + \sum _{i=1}^{m}  \frac{\overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i}{\overline{\boldsymbol{\sigma }}_i^2} \cdot \left( \frac{\partial \overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta })}{\partial \mathbf{x}_0} \delta \mathbf{x}_0 + \frac{\partial \overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta })}{\partial \boldsymbol{\theta }} \delta \boldsymbol{\theta } \right).  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

At the limit <math display="inline">\delta \mathbf{z}=(\delta \mathbf{x}_0,\delta \boldsymbol{\theta })^T \rightarrow \mathbf{0}</math> we obtain the gradient <math display="inline">\nabla L(\mathbf{z}) = (\nabla _{x_0}\, L(\mathbf{x}_0,\boldsymbol{\theta }), \nabla _\theta \, L(\mathbf{x}_0,\boldsymbol{\theta }) )^T</math> with the above derivatives distributed accordingly. In ([[#eq-9|9]]), matrices <math display="inline">\partial \overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta }) / \partial \mathbf{x}_0</math> and <math display="inline">\partial \overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta }) / \partial \boldsymbol{\theta }</math> are the Jacobians of the observable variables <math display="inline">\overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta })</math> with respect <math display="inline">\mathbf{x}_0</math> and <math display="inline">\boldsymbol{\theta }</math>, respectively. In the literature, their partial derivatives are commonly called ''sensitivities of the observables'', and their evaluation may require considerable computational effort. The first matrix is of size <math display="inline">no\times d</math> and the second matrix is of size <math display="inline">no\times np</math>, then the total number of partial derivatives  that we must compute in ([[#eq-9|9]]) is <math display="inline">m\times no (d + np)</math>. For instance, in the above example for the SEIRD model, if the number of observable variables is <math display="inline">no = 3</math>, and there are <math display="inline">m =13</math> experimental measurements, we have to compute <math display="inline">13 \times 3 (5 + 4) = 351</math> sensitivities at each iteration of a typical gradient or quasi-Newton method. The most basic method to compute these derivatives is the finite difference method (see <span id='citeF-14'></span>[[#cite-14|[14]]]), but it has some disadvantages. As mentioned before, we concentrate only in variational methods, which may be more sophisticated but they are commonly more efficient.

===3.1 Variational method to compute the sensitivities===

This method is also known as the ''forward approach'', because a forward dynamical systems is solved to compute the sensitivities. We consider first the calculation of the sensitivities with respect to the parameter <math display="inline">\boldsymbol{\theta }</math>, so let us denote the Jacobian <math display="inline">\partial \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }) / \partial \boldsymbol{\theta }</math> by <math display="inline"> J (t;\mathbf{x}_0,\boldsymbol{\theta })</math> for simplicity. Then

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  J (t;\mathbf{x}_0,\boldsymbol{\theta })\, \delta \boldsymbol{\theta }\approx \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }+ \delta \boldsymbol{\theta })  - \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }). </math>
|}
|}

Taking derivatives with respect to <math display="inline">t</math>, we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{J} (t;\boldsymbol{\theta })\,  \delta \boldsymbol{\theta }\approx \dot{\mathbf{x}}(t;\mathbf{x}_0,\boldsymbol{\theta }+ \delta \boldsymbol{\theta }) - \dot{\mathbf{x}}(t;\mathbf{x}_0,\boldsymbol{\theta }) </math>
|-
| style="text-align: center;" | <math> = \mathbf{f}(\mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }+\delta \boldsymbol{\theta }), \boldsymbol{\theta }+ \delta \boldsymbol{\theta }) - \mathbf{f}(\mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }), \boldsymbol{\theta }) </math>
|-
| style="text-align: center;" | <math> \approx \mathbf{f}_\mathbf{x}(\mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }), \boldsymbol{\theta })   J (t;\mathbf{x}_0,\boldsymbol{\theta })\, \delta \boldsymbol{\theta }+ \mathbf{f}_{\boldsymbol{\theta }} (\mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }), \boldsymbol{\theta }) \, \delta \boldsymbol{\theta }. </math>
|}
|}

Then, <math display="inline">\mathbf{x}(t) = \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta })</math> and <math display="inline">J(t) = J (t;\mathbf{x}_0,\boldsymbol{\theta })</math> satisfy the following variational system:

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t),\boldsymbol{\theta }), \quad t_0 < t \le t_f, </math>
|-
| style="text-align: center;" | <math> \mathbf{x}(t_0) = \mathbf{x}_0, </math>
|-
| style="text-align: center;" | <math> \dot{J} (t) = \mathbf{f}_\mathbf{x}(\mathbf{x}(t),\boldsymbol{\theta })  J(t) + \mathbf{f}_{\boldsymbol{\theta }} (\mathbf{x}(t),\boldsymbol{\theta }), \quad t_0 < t \le t_f,  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|-
| style="text-align: center;" | <math> J (t_0) = \mathbb{O}.  </math>
|}
|}

Last two relations in ([[#eq-10|10]]) form a system of matricial differential equations and its initial condition reflects the fact that <math display="inline">\mathbf{x}_0</math> does not depend on <math display="inline">\theta </math>, because

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> J(t_0) = J (t_0;\mathbf{x}_0,\boldsymbol{\theta }) = \partial \mathbf{x}(t_0;\mathbf{x}_0,\boldsymbol{\theta }) / \partial \boldsymbol{\theta }= \partial \mathbf{x}_0/ \partial \boldsymbol{\theta }= \mathbb{O} \quad \hbox{(the null matrix)}. </math>
|}
|}

A similar variational system is satisfied by the matrix with the sensitivities with respect <math display="inline">\mathbf{x}_0</math>. Therefore, with this approach, most of the computational work is concentrated in the solution of two matricial systems of differential equations and their evaluation at the experimental times <math display="inline">t_i</math>, <math display="inline">1\le i \le m</math>. The amount of computational work will accumulate at each new iteration of the optimization algorithm, which in some cases may be prohibitive.

===3.2 The adjoint method to compute the sensitivities===

The total variation of <math display="inline">\mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta })</math> with respect to <math display="inline">\mathbf{x}_0</math> and <math display="inline">\boldsymbol{\theta }</math> is given by

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }) = \mathbf{x}(t;\mathbf{x}_0 + \delta \mathbf{x}_0,\boldsymbol{\theta }+ \delta \boldsymbol{\theta }) - \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }) = \frac{\partial \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta })}{\partial \mathbf{x}_0} \delta \mathbf{x}_0 + \frac{\partial \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta })}{\partial \boldsymbol{\theta }} \delta \boldsymbol{\theta }.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

Our goal is to avoid the explicit calculation of the Jacobian matrices on the right hand side. Before we proceed further, let us first introduce the following Hilbert spaces:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>V = L^2([t_0,t_f];\mathbb{R}^d) = \left\{\mathbf{v}:[t_0,t_f] \rightarrow \mathbb{R}^d   \left\vert  \int _{t_0}^{t_f} ||\mathbf{v}(t)||^2\,dt < \infty \right.\right\}, </math>
|-
| style="text-align: center;" | <math> H = H^{1}([t_0,t_f];\mathbb{R}^d) = \{ \mathbf{v}\in V  |  d\mathbf{v}/ dt \in V \} , </math>
|}
|}

where the <math display="inline">||\mathbf{v}(t)||^2 = \mathbf{v}(t)^T\mathbf{v}(t) = \mathbf{v}(t) \cdot \mathbf{v}(t)</math> is the usual scalar product in <math display="inline">\mathbb{R}^d</math>. Then a natural inner product in <math display="inline">V</math> is given by <math display="inline">\langle \mathbf{u}, \mathbf{v}\rangle = \int _{t_0}^{t_f} \mathbf{u}(t) \cdot \mathbf{v}(t) \,dt</math> with induced norm given by <math display="inline">\Vert \mathbf{u}\Vert _V = \langle \mathbf{u}, \mathbf{u}\rangle ^{1/2}</math>.

Since <math display="inline">\mathbf{x}(t) = \mathbf{x}(t;\mathbf{x}_0,\boldsymbol{\theta }) \in H</math> satisfies the state equation ([[#eq-1|1]]), then the following inner product is null

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \int _{t_0}^{t_f} \left(\dot{\mathbf{x}}(t) - \mathbf{f}(\mathbf{x}(t),\mathbf{x}_0,\boldsymbol{\theta }) \right)\cdot \mathbf{p}(t)\, dt = 0, \quad \forall  \mathbf{p}\in H. </math>
|}
|}

Differentiating this expression, we get

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \int _{t_0}^{t_f} \!\! \delta \dot{\mathbf{x}}(t)\cdot \mathbf{p}(t)\, dt = \int _{t_0}^{t_f} \!\! \left[\, \mathbf{f}_\mathbf{x}(\mathbf{x}(t),\mathbf{x}_0,\boldsymbol{\theta })\,\delta \mathbf{x}(t) +  \mathbf{f}_{\mathbf{x}_0} (\mathbf{x}(t),\mathbf{x}_0,\boldsymbol{\theta })\,\delta \mathbf{x}_0 + \mathbf{f}_{\boldsymbol{\theta }} (\mathbf{x}(t),\mathbf{x}_0,\boldsymbol{\theta })\,\delta \boldsymbol{\theta }\, \right]\cdot \mathbf{p}(t)\, dt, </math>
|}
|}

where <math display="inline">\mathbf{f}_{\mathbf{x}}</math>, <math display="inline">\mathbf{f}_{\mathbf{x}_0}</math> and <math display="inline">\mathbf{f}_{\boldsymbol{\theta }}</math> are the Jacobians of <math display="inline">\mathbf{f}</math> with respect to <math display="inline">\mathbf{x}</math>, <math display="inline">\mathbf{x}_0</math> and <math display="inline">\boldsymbol{\theta }</math>, respectively. Integrating by parts the left hand side, and observing that on the right hand side <math display="inline">\left[\mathbf{f}_{\mathbf{x}}\, \delta \mathbf{x}\right]\cdot \mathbf{p}= \left[\mathbf{f}_{\mathbf{x}}^T\, \mathbf{p}\right]\cdot \delta \mathbf{x}</math>, <math display="inline">\left[\mathbf{f}_{\boldsymbol{\theta }}\, \delta \mathbf{x}_0 \right]\cdot \mathbf{p}= \left[\mathbf{f}_{\boldsymbol{\theta }}^T\, \mathbf{p}\right]\cdot \delta \mathbf{x}_0</math>, and <math display="inline">\left[\mathbf{f}_{\boldsymbol{\theta }}\, \delta \boldsymbol{\theta }\right]\cdot \mathbf{p}= \left[\mathbf{f}_{\boldsymbol{\theta }}^T\, \mathbf{p}\right]\cdot \delta \boldsymbol{\theta }</math>, we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{p}(t) \!\cdot \!  \delta \mathbf{x}(t) |_{t_0}^{t_f} - \int _{t_0}^{t_f}  \!\!\! \dot{\mathbf{p}}(t) \!\cdot \! \delta \mathbf{x}(t)\, dt  =   \int _{t_0}^{t_f} \!\!\left\{\left[\mathbf{f}_{\mathbf{x}}^T\, \mathbf{p}\right](t) \!\cdot \!  \delta \mathbf{x}(t) +  \left[\mathbf{f}_{\mathbf{x}_0}^T\, \mathbf{p}\right](t) \!\cdot \! \delta \mathbf{x}_0(t) + \left[\mathbf{f}_{\boldsymbol{\theta }}^T\, \mathbf{p}\right](t) \!\cdot \! \delta \boldsymbol{\theta }\right\}dt, </math>
|}
|}

and, assuming that <math display="inline">\mathbf{f}</math> does not depend explicitly of <math display="inline">\mathbf{x}_0</math> (it depends only  through <math display="inline">\mathbf{x}</math>, but the variation w.r.t. <math display="inline">\mathbf{x}_0</math> is already accounted in <math display="inline">\delta \mathbf{x}</math>), then

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{p}(t) \cdot \delta \mathbf{x}(t) |_{t_0}^{t_f} - \int _{t_0}^{t_f}   \left\{\dot{\mathbf{p}}(t) + \left[\mathbf{f}_{\mathbf{x}}^T\, \mathbf{p}\right](t) \right\}\cdot \delta \mathbf{x}(t)\, dt =    \int _{t_0}^{t_f}  \left[\mathbf{f}_{\boldsymbol{\theta }}^T\, \mathbf{p}\right](t) \cdot \delta \boldsymbol{\theta } \, dt,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where <math display="inline">\delta \mathbf{x}</math> is given by ([[#eq-11|11]]) and arise in the last term of ([[#eq-9|9]]). One way to avoid the explicit computation of ([[#eq-11|11]]) is forcing a relation with ([[#eq-12|12]]) by introducing the following adjoint equation

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>- \dot{\mathbf{p}}(t) = \mathbf{f}_{\mathbf{x}}(\mathbf{x}(t),\boldsymbol{\theta })^T\, \mathbf{p}(t) + \sum _{i=1}^{m} \frac{(\overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i) }{\overline{\boldsymbol{\sigma }}_i^2}  \, \boldsymbol{\delta }_\mathbf{D}(t-t_i),    t_f > t \ge t_0,   </math>
|-
| style="text-align: center;" | <math>   \mathbf{p}(t_f) = \mathbf{0},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

with <math display="inline">\boldsymbol{\delta }_\mathbf{D}(t-t_i)</math> the Dirac measure centred at <math display="inline">t_i</math>. This adjoint equation (a backward in time differential equation) contains all the information about the experimental data <math display="inline">\{ \overline{\mathbf{x}}_i\} _{i=1}^{m}</math>, and its variational formulation is obtained multiplying by a differentiable test function <math display="inline">\boldsymbol{\phi }(t)</math> and integrating:

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>- \int _{t_0}^{t_f}   \left\{\dot{\mathbf{p}}(t) + \left[\mathbf{f}_{\mathbf{x}}^T\, \mathbf{p}\right](t) \right\}\cdot \boldsymbol{\phi }(t)\, dt =   \int _{t_0}^{t_f} \sum _{i=1}^{m} \frac{(\overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i) }{\overline{\boldsymbol{\sigma }}_i^2}  \, \boldsymbol{\delta }_\mathbf{D}(t-t_i) \cdot \boldsymbol{\phi }(t) \, dt.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

The integral on the right hand side is equal to <math display="inline">\sum _{i=1}^{m} \frac{(\overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i) }{\overline{\boldsymbol{\sigma }}_i} \cdot \boldsymbol{\phi }(t_i)</math>. Choosing <math display="inline">\boldsymbol{\phi }(t) = \delta \mathbf{x}(t)</math> in ([[#eq-14|14]]) and substituting the result in ([[#eq-12|12]]), we obtain

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>-\mathbf{p}(t_0) \cdot \delta \mathbf{x}(t_0) + \sum _{i=1}^{m} \frac{(\overline{\mathbf{x}}(t_i;\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i) }{\overline{\boldsymbol{\sigma }}_i^2} \cdot \delta \mathbf{x}(t_i) =    \int _{t_0}^{t_f}  \left[\mathbf{f}_{\boldsymbol{\theta }}^T\, \mathbf{p}\right](t) \cdot \delta \boldsymbol{\theta } \, dt,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">\mathbf{p}(t)</math> is the solution of the adjoint equation. Finally, the substitution of this equation into ([[#eq-9|9]]), gives

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla _{x_0}\,  L(\mathbf{x}_0,\boldsymbol{\theta }) \cdot \delta \mathbf{x}_0 + \nabla _\theta \, L(\mathbf{x}_0,\boldsymbol{\theta })\cdot \delta \boldsymbol{\theta }</math>
|-
| style="text-align: center;" | <math>   \approx \frac{\mathbf{x}_0 - \mathbf{s}_0}{\boldsymbol{\sigma }_0^2}\, \cdot \delta  \mathbf{x}_0 + k(\mathbf{x}_0\cdot \mathbf{1} - N)\,\mathbf{1}\cdot \delta \mathbf{x}_0 + \mathbf{p}(t_0) \cdot \delta \mathbf{x}_0 + \left[\int _{t_0}^{t_f}  \left[\mathbf{f}_{\boldsymbol{\theta }}^T\, \mathbf{p}\right](t)\,dt \right] \cdot \delta \boldsymbol{\theta }.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

Therefore, the gradient of the objective function ([[#eq-6|6]]) is <math display="inline">\nabla L(\mathbf{z}) = (\nabla _{x_0}L(\mathbf{x}_0,\boldsymbol{\theta }), \nabla _{\boldsymbol{\theta }}L(\mathbf{x}_0,\boldsymbol{\theta }))^T</math>, with

<span id="eq-17"></span>
<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla _{x_0}L(\mathbf{x}_0,\boldsymbol{\theta }) = \frac{\mathbf{x}_0 - \mathbf{s}_0}{\boldsymbol{\sigma }_0^2} + k(\mathbf{x}_0\cdot \mathbf{1} - N)\,\mathbf{1}  + \mathbf{p}(t_0),   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|-
| style="text-align: center;" | <math> \nabla _{\boldsymbol{\theta }}L(\mathbf{x}_0,\boldsymbol{\theta }) = \int _{t_0}^{t_f} \mathbf{f}_{\boldsymbol{\theta }}(\mathbf{x}(t),\mathbf{x}_0,\boldsymbol{\theta })^T\, \mathbf{p}(t) \, dt,  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}
|}

where <math display="inline">\mathbf{x}(t)</math> solves the state equation ([[#eq-1|1]]) and <math display="inline">\mathbf{p}(t)</math> solves the adjoint equation ([[#eq-13|13]]).

Remark 3: Formulas ([[#eq-17|17]])-([[#eq-18|18]]) avoid the explicit calculation of the sensitivities, they only requiere the solution of the  state equation ([[#eq-1|1]]) and of the adjoint equation ([[#eq-13|13]]), regardless of the number of experimental data, <math display="inline">m</math>, the number of parameters, <math display="inline">np</math>, and of the observable variables, <math display="inline">no</math>. Furthermore, these same equations must be solved to compute either the gradient with respect <math display="inline">\mathbf{x}_0</math> or with respect <math display="inline">\boldsymbol{\theta }</math>, or both gradients simultaneously.

<span id='theorem-Lagrangian'></span>Remark 4: The solution of the adjoint equation turns out to be the Lagrange multiplier of the optimization problem, with objective function <math display="inline">L(\mathbf{x}_0,\boldsymbol{\theta })</math>, subject to the constraint ([[#eq-1|1]]). To show this property, let us introduce the Lagrangian function

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}(\mathbf{x}_0,\boldsymbol{\theta },\mathbf{p}) = L(\mathbf{x}_0,\boldsymbol{\theta }) + \int _{t_0}^{t_f} \mathbf{p}(t)^T \left[\mathbf{f}(t,\mathbf{x}(t), \boldsymbol{\theta }) - \dot{\mathbf{x}}(t) \right]\, dt,     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

where <math display="inline">\mathbf{p}</math> is the Lagrange multiplier associated to the given constraint, and the last term is the inner product of this multiplier with the restriction. Our goal is to compute <math display="inline">\partial L / \partial \boldsymbol{\theta }</math> from this expression. Formally, the second term on the right hand side of ([[#theorem-Lagrangian|4]]) is zero because the state <math display="inline">\mathbf{x}(t)</math> solves the ODE, therefore <math display="inline">\partial \mathcal{L}/ \partial \boldsymbol{\theta }= \partial L / \partial \boldsymbol{\theta }</math>. However, the differentiation of <math display="inline">\mathcal{L}</math> reveals additional information and gives more freedom. Doing integration by parts of the term <math display="inline">- \int _{t_0}^{t_f} \mathbf{p}(t)^T\dot{\mathbf{x}}(t)</math> in ([[#theorem-Lagrangian|4]]) first, and then taking the derivative with respect to <math display="inline">\boldsymbol{\theta }</math>, we obtain

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial \mathcal{L}}{\partial \boldsymbol{\theta }} = \frac{\partial L}{\partial \boldsymbol{\theta }} - \left.\mathbf{p}^T \frac{\partial \mathbf{x}}{\partial \boldsymbol{\theta }} \right|_{t_0}^{t_f} + \int _{t_0}^{t_f} \!\! \dot{\mathbf{p}}^T \frac{\partial \mathbf{x}}{\partial \boldsymbol{\theta }} dt +  \int _{t_0}^{t_f} \!\! \mathbf{p}^T \left[\mathbf{f}_{\mathbf{x}}\frac{\partial \mathbf{x}}{\partial \boldsymbol{\theta }} + \mathbf{f}_{\boldsymbol{\theta }}\right]dt </math>
|-
| style="text-align: center;" | <math>  = \sum _{i=1}^m \left(\frac{\overline{\mathbf{x}}(t_i,\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i}{\overline{\boldsymbol{\sigma }}_i^2}\right)^T \!\!\frac{\partial \overline{\mathbf{x}}}{\partial \theta } - \left.\mathbf{p}^T \frac{\partial \mathbf{x}}{\partial \boldsymbol{\theta }} \right|_{t_0}^{t_f} + \int _{t_0}^{t_f} \!\! \left(\dot{\mathbf{p}} + \mathbf{f}_\mathbf{x}^T \mathbf{p}\right)^T \frac{\partial \mathbf{x}}{\partial \boldsymbol{\theta }} dt + \int _{t_0}^{t_f} \!\! \mathbf{p}^T\mathbf{f}_{\boldsymbol{\theta }} \, dt.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

The first term on the right hand side is obtained directly from ([[#eq-6|6]]) and can be expressed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \sum _{i=1}^m \left(\frac{\overline{\mathbf{x}}(t_i,\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i}{\overline{\boldsymbol{\sigma }}_i^2}\right)^T \!\!\frac{\partial \overline{\mathbf{x}}}{\partial \theta }  = \int _{t_0}^{t_f} \sum _{i=1}^m \left(\frac{\overline{\mathbf{x}}(t_i,\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i}{\overline{\boldsymbol{\sigma }}_i^2}\right)^T \!\!\frac{\partial \overline{\mathbf{x}}}{\partial \theta } \,\delta _D(t-t_i) dt. </math>
|}
|}

Now, if <math display="inline">\mathbf{p}</math> satisfies the adjoint equation ([[#eq-13|13]]) then the sum of first and third terms in ([[#eq-20|20]]) vanish, and also the boundary term <math display="inline">\left.\mathbf{p}^T \frac{\partial \mathbf{x}}{\partial \boldsymbol{\theta }} \right|_{t_0}^{t_f} </math>, since <math display="inline">\mathbf{p}(t_f) = \mathbf{0}</math> and <math display="inline">\partial \mathbf{x}(t_0) /\partial \boldsymbol{\theta }= \mathbb{O}</math>. Therefore

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \frac{\partial \mathcal{L}}{\partial \boldsymbol{\theta }} = \frac{\partial L}{\partial \boldsymbol{\theta }} = \int _{t_0}^{t_f} \mathbf{p}^T\mathbf{f}_{\boldsymbol{\theta }} \, dt \quad \Longrightarrow \quad \nabla _{\boldsymbol{\theta }} L(\mathbf{x}_0,\boldsymbol{\theta }) = \int _{t_0}^{t_f} \mathbf{f}_{\boldsymbol{\theta }}(t,\mathbf{x}(t),\boldsymbol{\theta })^T\,\mathbf{p}(t) \, dt. </math>
|}
|}

A similar development can be applied to obtain <math display="inline">\partial L/\partial \mathbf{x}_0</math>.

===3.3 Solution of the adjoint equation and computation of the gradient===

The adjoint equation ([[#eq-13|13]]) is a system of ODE with backward in time propagation. We apply a change of variable from the symmetric relation <math display="inline">t_f - \tau = t - t_0</math>:

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p(t) = p(t_f - (\tau{-}t_0)) \equiv p_A (\tau ) \quad \Longrightarrow \quad  \dot{p}(t) = -\dot{p}_A (\tau ),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

obtaining the equivalent system with forward in time dynamics

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{\mathbf{p}}_A (\tau ) = \mathbf{f}_{\mathbf{x}} (\mathbf{x}_A(\tau ),\boldsymbol{\theta })^T\, \mathbf{p}_A(\tau ) + \sum _{i=1}^{m} \frac{\overline{\mathbf{x}}_A(\tau _i,\mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_{m+1-i}}{\overline{\boldsymbol{\sigma }}_{m+1-i}^2}  \, \delta _D(\tau _i-\tau ) ,\quad   t_0 \le \tau < t_f,  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|-
| style="text-align: center;" | <math>  \mathbf{p}_A (\tau _0) = \mathbf{0}, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}
|}

which can be solved with any usual numerical ODE solver like Runge&#8211;Kutta methods.

For the SEIRD model, if the observable variables are <math display="inline">\overline{\mathbf{x}}(t) = (I(t),R(t),D(t))^T</math>, the adjoint equation ([[#eq-22|22]]) takes the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[\!\! \begin{array}{c}\dot{p}_{_{A1}}(\tau ) \\ \dot{p}_{_{A2}}(\tau ) \\ \dot{p}_{_{A3}}(\tau ) \\ \dot{p}_{_{A4}}(\tau ) \\ \dot{p}_{_{A5}}(\tau )  \end{array} \!\! \right]= \left[\!\! \begin{array}{ccrcr}-\frac{\alpha \,I(\tau )}{N} & \frac{\alpha \,I(\tau )}{N} & 0 & 0 & 0 \\ 0 & -\beta & \beta & 0 & 0 \\ -\frac{\alpha \,S(\tau )}{N} & \frac{\alpha \,S(\tau )}{N} & -\gamma & \gamma -\mu & \mu \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{array} \!\! \right]\!\!  \left[\!\! \begin{array}{c}p_{_{A1}}(\tau ) \\ p_{_{A2}}(\tau ) \\ p_{_{A3}}(\tau ) \\ p_{_{A4}}(\tau ) \\ p_{_{A5}}(\tau ) \end{array} \!\! \right]+  \left[\!\! \begin{array}{c}0 \\ 0 \\ \sum _{i=1}^{m} \frac{I(\tau _i) - I_{m+1-i}}{\sigma _{i3}^2} \delta _D(\tau _i-\tau ) \\ \sum _{i=1}^{m} \frac{R(\tau _i) - R_{m+1-i}}{\sigma _{i4}^2} \delta _D(\tau _i-\tau ) \\  \sum _{i=1}^{m} \frac{D(\tau _i) - D_{m+1-i}}{\sigma _{i5}^2} \delta _D(\tau _i-\tau ) \end{array} \!\! \right], </math>
|}
|}

with <math display="inline">\gamma = 1/T_1</math>, <math display="inline">\mu = f/T_1</math> in ([[#theorem-SEIRD|1]]). After solving this forward adjoint equation we recover the solution <math display="inline">\mathbf{p}</math> of the backward adjoint equation ([[#eq-13|13]]) with <math display="inline">\mathbf{p}(t_i) = \mathbf{p}_A(\tau _{m+1-i})</math>, <math display="inline">i = 1,\ldots ,m</math> or by interpolation and using ([[#eq-21|21]]) for other values of <math display="inline">t</math>.

The last step to compute the gradient is the calculation of the integral in ([[#eq-18|18]]). Observe that

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{f}_{\boldsymbol{\theta }}(\mathbf{x}(t),\boldsymbol{\theta })^T\, \mathbf{p}(t) = \left[\!\! \begin{array}{rrrr}-\dfrac{S(t)\,I(t)}{N} & \dfrac{S(t)\,I(t)}{N} & 0 & 0 \\ 0 & -E(t) & E(t) & 0 \\  0 & 0 & -I(t) & I(t) \end{array} \!\! \right]\left[\!\! \begin{array}{c}p_1(t) \\ p_2(t) \\ p_3(t) \\ p_4(t) \end{array} \!\! \right]</math>
|-
| style="text-align: center;" | <math>  = \left[\begin{array}{c}\dfrac{S(t)\,I(t)}{N} \, ( p_2(t) - p_1(t) ) \\ E(t) \, (p_3(t) - p_2(t)) \\ I(t) \, (p_4(t) - p_3(t)) \end{array}  \right]\equiv \mathbf{H}(t) </math>
|}
|}

Then, using the Simpson's rule in a set of given times (nodes of the mesh or interpolated times) we have

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{t_0}^{t_f} \mathbf{f}_{\boldsymbol{\theta }}(\mathbf{x}(t),\boldsymbol{\theta })^T\, \mathbf{p}(t)\, dt = \frac{\Delta t}{6}\sum _{j=1}^{N-1} \left\{\mathbf{H}(t_{j-1}) + 4 \mathbf{H} (t_j) + \mathbf{H}(t_{j+1}) \right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

==4 Numerical algorithms for the optimization problem==

===4.1 Conjugate gradient algorithm===

This algorithm is one of the most important algorithms for quadratic optimization problems with positive definite Hessians and for unconstrained continuous convex optimization. It may be considered as a variant of gradient descent where the search directions are generated progressively based on the orthogonality of the residuals and conjugacy of the search directions. The conjugate directions are calculated at each iteration  a linear combination of the most recent negative gradient and the last conjugate direction, as indicated in step 9 of Algorithm 1 bellow. Its computational cost is comparable to steepest descent, but it has faster convergence, specially for ill conditioned problems, <span id='citeF-9'></span>[[#cite-9|[9]]]. 
{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:'''Initialization''' 1. Initial guess <math display="inline">\mathbf{z}^0 = \widehat{\mathbf{z}} = (\widehat{\mathbf{x}}_0,\widehat{\boldsymbol{\theta }})^T</math>. 2. Initial gradient <math display="inline">\mathbf{g}^0 = \nabla L(\mathbf{z}^0)</math>. 3. Initial direction <math display="inline">\mathbf{d}^0 = -\mathbf{g}^0</math>. '''Descent''' For <math display="inline">\ell \ge 0</math>, given <math display="inline">\mathbf{z}^\ell </math>, <math display="inline">\mathbf{g}^\ell </math>, <math display="inline">\mathbf{d}^\ell </math>, find <math display="inline">\mathbf{z}^{\ell{+1}}</math>, <math display="inline">\mathbf{g}^{\ell{+1}}</math>, <math display="inline">\mathbf{d}^{\ell{+1}}</math>, doing the following 4. Find <math display="inline">\rho _\ell = \displaystyle \arg \min _{\rho \ge 0} \varphi (\rho ) = L(\mathbf{z}^\ell + \rho \,\mathbf{d}^\ell )</math> 5. Update <math display="inline">\mathbf{z}^{\ell{+1}} = \mathbf{z}^\ell + \rho _\ell \,\mathbf{d}^\ell </math>. 6. Evaluate <math display="inline">\mathbf{g}^{\ell + 1} = \nabla L(\mathbf{z}^{\ell + 1})</math> '''Convergence test and new direction''' <math>\Vert \mathbf{g}^{\ell +1} \Vert \le \epsilon \Vert \mathbf{g}^0 \Vert </math>  7. Take <math>\mathbf{z}^* = \mathbf{z}^{\ell +1}</math>. Stop and exit.   8. Evaluate <math>\beta _\ell =  \dfrac{\mathbf{g}^{\ell +1}\cdot \mathbf{g}^{\ell +1}}{\mathbf{g}^\ell \cdot  \mathbf{g}^\ell }</math>9. Update <math>\mathbf{d}^{\ell + 1} = -\mathbf{g}^\ell + \beta _\ell \, \mathbf{d}^\ell </math>10. Make <math>\ell = \ell + 1</math> and go back to 4.  


|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-1'></span>'''Algorithm. 1''' Conjugate gradient algorithm
|}

If <math display="inline">\beta _\ell = 0</math> for all <math display="inline">\ell </math> we recover steepest descent. Some variants for computing <math display="inline">\beta _\ell \ne 0</math>, include:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{ll}\bullet  \hbox{Fletcher--Reeves: } \beta _\ell =  \dfrac{\mathbf{g}^{\ell +1}\cdot \mathbf{g}^{\ell +1}}{\mathbf{g}^\ell \cdot  \mathbf{g}^\ell } & \quad  \bullet  \hbox{ Polak--Ribiere: } \beta _\ell = \dfrac{\mathbf{g}^{\ell{+1}} \cdot (\mathbf{g}^{\ell{+1}} - \mathbf{g}^\ell )}{\mathbf{g}^\ell \cdot  \mathbf{g}^\ell } \\ \bullet  \hbox{ Hestenes--Stiefel: } \beta _\ell = \dfrac{\mathbf{g}^{\ell{+1}} \cdot (\mathbf{g}^{\ell{+1}} - \mathbf{g}^\ell )}{\mathbf{d}^\ell \cdot (\mathbf{g}^{\ell{+1}} - \mathbf{g}^\ell )} & \quad   \bullet  \hbox{ Dai--Yuan: } \beta _\ell = \dfrac{\mathbf{g}^{\ell +1}\cdot \mathbf{g}^{\ell +1}}{\mathbf{d}^\ell \cdot (\mathbf{g}^{\ell{+1}} - \mathbf{g}^\ell )} \end{array} </math>
|}
|}

We use the short notations F&#8211;R, P&#8211;R, H&#8211;S, D&#8211;Y, for these variants.

===4.2 BFGS algorithm===

This is one of the most popular quasi-Newton algorithms for nonlinear optimization. It is usually more effective because it involves information about curvature, besides the information about the gradient. The curvature information is incorporated by approximating the Hessian matrix during the iteration process, which formaly plays the role of  preconditioner for the gradient. The general idea about these methods comes from the second order approximation:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> L(\mathbf{z}^\ell + \Delta \mathbf{z}) \approx L(\mathbf{z}^\ell ) + \nabla L(\mathbf{z}^\ell )^T \Delta \mathbf{z}+ \frac{1}{2}   \Delta \mathbf{z}^T A\, \Delta \mathbf{z}, </math>
|}
|}

with a matrix <math display="inline">A</math> close to the Hessian. Then the so called secant condition is obtained:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \nabla L(\mathbf{z}^\ell + \Delta \mathbf{z}) \approx \nabla L(\mathbf{z}^\ell ) + A\, \Delta \mathbf{z}. </math>
|}
|}

Forcing the gradient to be zero (to look for a minimum), we obtain the Newton step

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \Delta \mathbf{z}= - H\, \nabla L(\mathbf{z}^\ell ), \quad \hbox{where} \quad H \approx A^{-1}  </math>
|}
|}

and the following search direction <math display="inline">\mathbf{d}^\ell </math> is obtained;

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{d}^\ell = -H\,\mathbf{g}^\ell , \quad \hbox{with} \quad \mathbf{d}^\ell = \Delta \mathbf{z},\quad \mathbf{g}^\ell = \nabla L(\mathbf{z}^\ell{).} </math>
|}
|}

The Hessian and its inverse are updated at every iteration adding rank&#8211;one or rank&#8211;two matrices. The BFGS algorithm updates the approximated Hessian with a rank&#8211;two matrix as showing in step 9 of Algorithm 2 bellow. 
{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:'''Initialization''' 1. Initial guess ans initial Hessian: <math display="inline">\mathbf{z}^0</math> given, and <math display="inline">H^0 = I</math>. 2. Initial gradient <math display="inline">\mathbf{g}^0 = \nabla L(\mathbf{z}^0)</math>. 3. Initial direction <math display="inline">\mathbf{d}^0 = -\mathbf{g}^0</math>. '''Descent''' For <math display="inline">\ell \ge 0</math>, given <math display="inline">\mathbf{z}^\ell </math>, <math display="inline">\mathbf{g}^\ell </math>, <math display="inline">\mathbf{d}^\ell </math>, <math display="inline">H^\ell </math> find <math display="inline">\mathbf{z}^{\ell{+1}}</math>, <math display="inline">\mathbf{g}^{\ell{+1}}</math>, <math display="inline">\mathbf{d}^{\ell{+1}}</math>, <math display="inline">H^{\ell +1}</math>, doing the following 4. Find <math display="inline">\rho _\ell = \displaystyle \arg \min _{\rho \ge 0} \varphi (\rho ) = L(\mathbf{z}^\ell + \rho \,\mathbf{d}^\ell )</math> 5. Update <math display="inline">\mathbf{z}^{\ell{+1}} = \mathbf{z}^\ell + \rho _\ell \,\mathbf{d}^\ell </math>. 6. Evaluate <math display="inline">\mathbf{g}^{\ell + 1} = \nabla L(\mathbf{z}^{\ell + 1})</math> '''Convergence test and new direction''' <math>\Vert \mathbf{g}^{\ell +1} \Vert \le \epsilon \Vert \mathbf{g}^0 \Vert </math>  7. Do <math>\mathbf{z}^* = \mathbf{z}^{\ell +1}</math>. Stop and exit.   8. Evaluate <math>\mathbf{u}= \Delta \mathbf{z}^\ell = \mathbf{z}^{\ell +1} - \mathbf{z}^\ell </math> and <math>\mathbf{v}= \Delta \mathbf{g}^\ell = \mathbf{g}^{\ell +1} - \mathbf{g}^\ell </math>  9. Update <math>\displaystyle H^{\ell{+1}} = \left(I - \frac{\mathbf{v}\mathbf{u}^T}{\mathbf{u}^T\mathbf{v}} \right)H^\ell \left( I - \frac{\mathbf{u}\mathbf{v}^T}{\mathbf{u}^T\mathbf{v}} \right)+ \frac{\mathbf{v}\mathbf{v}^T}{\mathbf{u}^T\mathbf{v}} </math>  10. Update <math>\mathbf{d}^{\ell + 1} = -H^{\ell{+1}} \mathbf{g}^{\ell +1}</math>  11. Make <math>\ell = \ell +1</math> and go back to 4  


|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-2'></span>'''Algorithm. 2''' BFGS algorithm
|}

Remark 5: Another very popular method for least squares problems is the Gauss-Newton method and it variant, the Levenverg-Marquadt method (see <span id='citeF-9'></span>[[#cite-9|[9]]]), where the Jacobian of the residual vector <math display="inline">\mathbf{r} = \left\{\left\Vert \frac{\overline{\mathbf{x}}(t_i; \mathbf{x}_0,\boldsymbol{\theta }) - \overline{\mathbf{x}}_i}{\overline{\boldsymbol{\sigma }}_i} \right\Vert ^2 \right\}_{i=1}^m</math> with respect to <math display="inline">\mathbf{x}_0</math> and <math display="inline">\boldsymbol{\theta }</math> must be computed at each iteration. These partial derivatives can also be computed with variational methods.

===4.3 Line search===

The most critical step in Algorithms [[#algorithm-1|1]] and [[#algorithm-2|2]] is the solution of the one dimensional optimization problem at step 4. This is not a trivial step and requires careful treatment. There are several algorithms for this problem, the most common are line search methods and trust region methods. Some classic references are <span id='citeF-8'></span>[[#cite-8|[8]]] and <span id='citeF-9'></span>[[#cite-9|[9]]], or the more recent one <span id='citeF-15'></span>[[#cite-15|[15]]], while some publications, e,g, <span id='citeF-1'></span>[[#cite-1|[1]]] and <span id='citeF-10'></span>[[#cite-10|[10]]], show that this topic is still under development.

Given that we have a efficient way to compute the derivative <math display="inline">\varphi ^\prime (\rho ) = \nabla L(\mathbf{z}^\ell + \rho \,\mathbf{d}^\ell ) \cdot \mathbf{d}^\ell </math>, we may use the secant method, whose iteration formula is

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho _{k+1} = \rho _k - \frac{ \varphi ^{\,\prime }(\rho _k) (\rho _k - \rho _{k-1}) }{ \varphi ^{\,\prime }(\rho _k) - \varphi ^{\,\prime }(\rho _{k-1}) } = \frac{ \rho _{k-1} \varphi ^{\,\prime }(\rho _k) - \rho _k \varphi ^{\,\prime }(\rho _{k-1})}{\varphi ^{\,\prime }(\rho _k) - \varphi ^{\,\prime }(\rho _{k-1})}, \qquad k = 1,2,\ldots  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

with two initial values <math display="inline">\rho _0 = 0</math> and <math display="inline">\rho _1 \approx \epsilon \, \frac{|| \mathbf{z}||}{|| \mathbf{d}||}</math>, <math display="inline">\epsilon < 1</math>. The initial value <math display="inline">\rho _1</math> takes into account a proper scaling (see step 5). Computing <math display="inline">\varphi ^{\,\prime } (\rho )</math> requires solving the state equation ([[#eq-1|1]]) with initial conditions <math display="inline">\mathbf{x}_0^\ell + \rho \mathbf{d}_0^\ell </math> and parameter <math display="inline">\boldsymbol{\theta }^\ell + \rho \mathbf{d}_{\boldsymbol{\theta }}^\ell </math>, obtaining a solution that we call <math display="inline">\mathbf{x}_\rho ^\ell </math>, and then solving the corresponding adjoint equation  ([[#eq-13|13]]) with <math display="inline">\mathbf{x}= \mathbf{x}_\rho ^\ell </math> and <math display="inline">\boldsymbol{\theta }= \boldsymbol{\theta }^\ell + \rho \mathbf{d}_{\boldsymbol{\theta }}^\ell </math>, obtaining the solution <math display="inline">\mathbf{p}_\rho ^\ell </math>. Thus

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varphi ^{\,\prime } (\rho ) \! = \! \left[\!\! \begin{array}{c}\mathbf{p}_\rho ^\ell (t_0) + \dfrac{\overline{\mathbf{x}}_0^\ell + \rho \,\overline{\mathbf{d}}_0^\ell - \overline{\mathbf{y}}_0}{\overline{\boldsymbol{\sigma }}_0^2} + k \left([\mathbf{x}_0^\ell + \rho \mathbf{d}_0^\ell ]\cdot \mathbf{1} - N \right)\, \mathbf{1} \\  \displaystyle \int _{t_0}^{t_f} \mathbf{f}_{\boldsymbol{\theta }} ( \mathbf{x}_\rho ^\ell , \boldsymbol{\theta }^\ell + \rho \,\mathbf{d}_\theta ^\ell )^T \mathbf{p}_\rho ^\ell \, dt \end{array} \!\!\right]\!\cdot \! \left[\!\! \begin{array}{c}\mathbf{d}_0^\ell \\ \mathbf{d}_\theta ^\ell \end{array} \!\!\right].  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

Remark 6: The secant method for line search may be improved, incorporating standard bracketing strategies that keeps track of upper and lower bounds for the location of the root <span id='citeF-1'></span>[[#cite-1|[1]]]. We will try this enhancement in a future work.

<span id='theorem-Newton_line_search'></span>Remark 7:  Newton's method for line is also possible. However, it is more costly because we must compute <math display="inline">\varphi ^{\,\prime \prime } (\rho )</math>, i.e. the derivative of ([[#eq-26|26]]) with respect to <math display="inline">\rho </math> . This operation involves the solution of another two systems of ODE at each iteration: one forward-in-time problem (related to the state equation) and a backward-in-time problem (related to the adjoint equation), where the Jacobians <math display="inline">\mathbf{f}_\mathbf{x}</math> and <math display="inline">\mathbf{f}_{\boldsymbol{\theta }}</math> must be evaluated at different values. We leave this task for a future work.

==5 Numerical results for the SEIRD model==

To validate the fitting model and the proposed optimization algorithms we consider the SEIRD model, described in Section 2. Our base or reference true solution is the one obtained numerically in the time interval is <math display="inline">[t_0,t_f] = [40,80]</math> (days), parameters <math display="inline">\boldsymbol{\theta }=(\alpha ,\beta ,\gamma ,\mu )^T = (1,1/7,1/5,1/70)^T</math>, initial condition <math display="inline">\mathbf{x}_0 = (91647, 4853, 1755, 1620, 125)^T</math> and total population <math display="inline">N = 10^5</math>, shown in Figure [[#img-1|1]]. Synthetic data is generated adding white noise to the true solution with the random Gaussian generator of Matlab, with zero mean and standard deviations <math display="inline">\mathbf{s}(t_i) = noise\_level*\mathbf{x}(t_i)</math> at the times <math display="inline">t_i</math> where we suppose to have experimental measurements. The proposed model and numerical algorithms are tested at three levels of noise, 0.05, 0.1, and 0.2. We will divide the experiments in two parts: 1) only the vector parameter <math display="inline">\boldsymbol{\theta }</math> is unknown, 2) both the initial conditions <math display="inline">\mathbf{x}_0</math> and the vector parameter <math display="inline">\boldsymbol{\theta }</math> are unknown.

<span id='theorem-tabla1'></span>Example 2: Case when <math display="inline">\mathbf{x}_0</math> is known and <math display="inline">\boldsymbol{\theta }</math> is unknown, <math display="inline">noise\_level = 0.1</math>.

In many problems we are interested in recovering the vector of unknown parameters <math display="inline">\boldsymbol{\theta }</math>, assuming that we are given the exact initial conditions <math display="inline">\mathbf{x}_0</math> and the experimental data <math display="inline">\{ \overline{\mathbf{x}}_i\} _{i=1}^m</math> at the corresponding times <math display="inline">t_i</math>. We consider synthetic experimental noisy data with <math display="inline">noise\_level = 0.1</math>, where the observable variables are <math display="inline">\overline{\mathbf{x}}(t_i) = (I_i,R_i,D_i)^T</math> in the time window <math display="inline">t = 40+i</math>, <math display="inline">1\le i \le 13</math>. Table [[#theorem-tabla1|2]] shows the numerical results obtained with the CG-algorithm (F-R variant) and the BFGS-algorithms, with initial guess <math display="inline">\theta ^0 = (1.4,0.09,0.2,0.001)^T</math> and tolerance <math display="inline">\epsilon = 10^{-8}</math> to stop the iterations. The relative error is computed component-wise.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |Table. 1 <span style="text-align: center; font-size: 75%;">Numerical results with the CG and BFGS algorithms.</span>
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |    Method                      
| style="border-left: 2px solid;border-right: 2px solid;" | CG (F-R variant) 
| style="border-left: 2px solid;" | BFGS 
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" | <math display="inline">\boldsymbol{\theta }^0</math>                
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(1.4,0.09,0.2,0.001)</math>
| style="border-left: 2px solid;" | <math>(1.4,0.09,0.2,0.001)</math>
|-
| style="border-right: 2px solid;" | Data time window              
| style="border-left: 2px solid;border-right: 2px solid;" | <math>t_i = 40 +i</math>, <math display="inline">1\le i\le 13</math>      
| style="border-left: 2px solid;" | <math>t_i = 40 + i</math>, <math display="inline">1\le i \le 13</math>  
|-
| style="border-right: 2px solid;" | <math>\epsilon </math>, no. iters. 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\epsilon = 10^{-8}</math>, 168  
| style="border-left: 2px solid;" | <math>\epsilon = 10^{-8}</math>, 14 
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |  Computed <math display="inline">\boldsymbol{\theta }</math>   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(1.0769, 0.1425,0.2180,0.0143)</math>
| style="border-left: 2px solid;" | <math>(1.0768,0.1426,0.2179,0.0143)</math>
|- style="border-bottom: 2px solid;"
| style="border-right: 2px solid;" | Relative error             
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(0.0769,0.0022,0.0898,4.8e(-5))</math>
| style="border-left: 2px solid;" | <math>(0.0768,0.0021,0.0897,0.0001) \quad </math>

|}

We obtain almost the same numerical value for the computed <math display="inline">\boldsymbol{\theta }</math> with both algorithms, the main difference is the number of iterations for each algorithm to achieve convergence to the given tolerance. Overall, the most important feature in this experiment is that the numerical computation is stable and the relative error is smaller than the noise level 0.1 (10%) for each parameter.  Figure [[#img-2|2]] (left), shows the behaviour of <math display="inline">||\nabla L(\boldsymbol{\theta })||</math> (logarithmic scale) and clearly show the faster convergence of the BFGS algorithm. Figure [[#img-2|2]] (right) illustrate the convergence of the parameters <math display="inline">\boldsymbol{\theta }^\ell = \left(\alpha ^\ell ,\beta ^\ell ,\gamma ^\ell ,\mu ^\ell \right)^T</math>. Finally, Figure [[#img-3|3]] shows the dynamics of the true solution <math display="inline">\mathbf{x}</math> (continuous lines) and the numerical solution obtained with the computed parameters (dashed lines), along with the noisy data (points).  We want to emphasize that the initial value for parameter <math display="inline">\gamma ^0 = 0.2</math> is the exact one, since this a parameter is usually assumed to be known. However the numerical algorithms converge for other initial values, like <math display="inline">\gamma ^0 = 0.1</math> or 0.3 with a similar number of iterations. <div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-gradient_10_CG_BFGS.png|210px|]]
|[[Image:Draft_Juarez Valencia_878068377-parameters_10_CG_BFGS.png|210px|<span style="text-align: center; font-size: 75%;">Left: graph of (∇L(θ)^\ell ) against iteration value \ell . Right figure: Convergence of θ^\ell = (α^\ell ,β^\ell ,γ^\ell ,μ^\ell )<sup>T</sup> with respect to iteration \ell  for CG (continuous lines) and for BFGS (dashed lines).</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 2:''' <span style="text-align: center; font-size: 75%;">Left: graph of <math>\log \left(\nabla L(\boldsymbol{\theta })^\ell \right)</math> against iteration value <math>\ell </math>. Right figure: Convergence of <math>\boldsymbol{\theta }^\ell = \left(\alpha ^\ell ,\beta ^\ell ,\gamma ^\ell ,\mu ^\ell \right)^T</math> with respect to iteration <math>\ell </math> for CG (continuous lines) and for BFGS (dashed lines).</span>
|}
<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-x_10_CG_BFGS.png|240px|<span style="text-align: center; font-size: 75%;">The dynamics of the true solution x(t) = (S(t),E(t),I(t),R(t),D(t) )<sup>T</sup> (continuous lines) and the one obtained with computed θ (dashes lines), along with noisy data for the observable variables (I<sub>i</sub>,R<sub>i</sub>,D<sub>i</sub>) (points).</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' <span style="text-align: center; font-size: 75%;">The dynamics of the true solution <math>\mathbf{x}(t) = \left(S(t),E(t),I(t),R(t),D(t) \right)^T</math> (continuous lines) and the one obtained with computed <math>\boldsymbol{\theta }</math> (dashes lines), along with noisy data for the observable variables <math>(I_i,R_i,D_i)</math> (points).</span>
|}

<span id='theorem-ajuste_a_N'></span>Example 3: Case when both <math display="inline">\mathbf{x}_0</math> and <math display="inline">\boldsymbol{\theta }</math> are unknown, <math display="inline">noise\_level = 0.1</math>.

This problem needs a more careful treatment, we now have to compute nine parameters instead of four and the scale of both unknowns is different: <math display="inline">\mathbf{x}_0</math> may have components of the order of <math display="inline">10^{5}</math> while <math display="inline">\boldsymbol{\theta }</math> has components of the order at most <math display="inline">10^0 = 1</math>. The initial guess <math display="inline">\boldsymbol{\theta }^0</math> to start iterations is the same for both algorithms (CG and BFGS). However, the election of the initial guess <math display="inline">\mathbf{x}_0^0</math> is more subtle. We first fix the value of <math display="inline">\mathbf{s}_0</math> in model ([[#eq-6|6]]) with the following formula

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{s}_0 = round(N/\widehat{N})\, \widehat{\mathbf{x}}_0 \quad \hbox{where} \quad \widehat{N} = \sum _{i=1}^5 \widehat{x}_{i\,0},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

where <math display="inline">\widehat{\mathbf{x}}_0</math> is obtained from the noisy data at <math display="inline">t =40</math>. The adjustment in ([[#theorem-ajuste_a_N|3]]) ensures that the sum of the components of <math display="inline">\mathbf{s}_0</math> be equal to <math display="inline">N = 10^5</math>. Observe that <math display="inline">\widehat{N}</math> is not guaranteed to be equal to <math display="inline">N</math>, because experimental data have inherent noise. Finally, we choose <math display="inline">\mathbf{x}_0^0 = \mathbf{s}_0</math> for the CG algorithm and <math display="inline">\mathbf{x}_0^0 = \widehat{\mathbf{x}}_0</math> for the BFGS algorithm. The CG algorithm does not always converge properly with an arbitrary initial value like <math display="inline">\mathbf{x}_0^0 = \widehat{\mathbf{x}}_0</math>. Table [[#table-2|2]] summarize  the numerical results.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 <span style="text-align: center; font-size: 75%;">Numerical results for computed <math>\boldsymbol{\theta }</math> and <math>\mathbf{x}_0</math> with CG (P&#8211;R) and BFGS algorithms.</span>
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |    Method 
| style="border-left: 2px solid;border-right: 2px solid;" | CG (P&#8211;R variant) 
| style="border-left: 2px solid;" | BFGS  
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" | <math display="inline">\boldsymbol{\theta }^0</math>                    
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(1.40, 0.09,0.20,0.001)</math>
| style="border-left: 2px solid;" | <math>(1.40, 0.09,0.20,0.001)</math>
|-
| style="border-right: 2px solid;" | <math>\mathbf{x}_0^0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(91503,4978,1872,1643,132)</math>
| style="border-left: 2px solid;" | <math>(91503,4978,1872,1643,132)</math>
|-
| style="border-right: 2px solid;" | Data time window          
| style="border-left: 2px solid;border-right: 2px solid;" | <math>t_i = 40 + (2i-1)</math>, <math display="inline">1\le i\le 8</math>     
| style="border-left: 2px solid;" | <math>t_i = 40 +(2i-1)</math>, <math display="inline">1\le i\le 8</math> 
|-
| style="border-right: 2px solid;" | <math>\epsilon </math>,  no. iters.  
| style="border-left: 2px solid;border-right: 2px solid;" | <math>10^{-5}</math>, 229 
| style="border-left: 2px solid;" | <math>10^{-6}</math>, 41  
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |  Computed <math display="inline">\boldsymbol{\theta }</math>       
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(1.0075,0.1420,0.2018,0.0134)</math>
| style="border-left: 2px solid;" | <math>(1.0077,0.1419,0.2003,0.0132)</math>
|-
| style="border-right: 2px solid;" | Relative error                 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(0.0075,0.0058,0.0090,0.0590)</math>
| style="border-left: 2px solid;" | <math>(0.0077,0.0064,0.0016,0.0758)</math>
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |  Computed <math display="inline">\mathbf{x}_0</math>    
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(91386,4972,1870,1641,132)</math>
| style="border-left: 2px solid;" | <math>(91340,4996,1844,1681,139)</math>
|- style="border-bottom: 2px solid;"
| style="border-right: 2px solid;" | Relative error                 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(0.0029,0.0245,0.0654,0.0128,0.0548)</math>
| style="border-left: 2px solid;" | <math>(0.0034,0.0294,0.0505,0.0379,0.1148)</math>

|}

This time the CG algorithm does not admit tolerances smaller than <math display="inline">\epsilon = 10^{-5}</math> and also the P-R variant to compute <math display="inline">\beta _{\ell }</math> at step 8 turn out to be more efficient than F-R. We observe that the window time for experimental data is wider but with less data points for both algorithms. The computed <math display="inline">\boldsymbol{\theta }</math> obtained with both algorithms is almost the same, but the computed <math display="inline">\mathbf{x}_0</math> exhibits more discrepancy. This lack of stability on the estimation of initial conditions from noisy data is already known by the scientific community. In fact, if <math display="inline">\mathbf{f}</math> is Lipschitz continuous with constant <math display="inline">L>0</math>, then the sensitivity of the solution  of the system ([[#eq-1|1]]) with respect to initial conditions is given by <math display="inline">||\mathbf{x}(t;\mathbf{x}_0) - \mathbf{x}(t;\mathbf{x}_0 + \delta \mathbf{x}_0)|| \le e^{L(t-t_0)} ||\delta \mathbf{x}_0||</math>. 

Figures [[#img-4|4]] and [[#img-5|5]] show information about the convergence of both methods like in the previous example. We have only added in Figure [[#img-5|5]] (left) a plot that shows convergence of <math display="inline">\mathbf{x}_0^\ell </math>. The main difference with respect to the previous example is that convergence not only is slower for both methods but also it is not as smooth as before. The convergence curves oscillate a lot more due to a destabilization effect of the unknown initial conditions. However, Figure [[#img-5|5]] (right) shows that the overall numerical reproduction of the true solution <math display="inline">\mathbf{x}(t)</math> is much better than in the previous example (all curves are very close to each other).  <div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-g2_10_CG_BFGS.png|210px|]]
|[[Image:Draft_Juarez Valencia_878068377-p2_10_CG_BFGS.png|210px|<span style="text-align: center; font-size: 75%;">Left: gradient behaviour obtained with the CG (blue line) and the BFGS (red line) algorithms. Right: convergence of θ with CG (continuous lines) and BFGS (dashed lines) algorithms.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 4:''' <span style="text-align: center; font-size: 75%;">Left: gradient behaviour obtained with the CG (blue line) and the BFGS (red line) algorithms. Right: convergence of <math>\boldsymbol{\theta }</math> with CG (continuous lines) and BFGS (dashed lines) algorithms.</span>
|}
<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-x0_10_CG_BFGS.png|240px|]]
|[[Image:Draft_Juarez Valencia_878068377-x2_10_CG_BFGS.png|240px|<span style="text-align: center; font-size: 75%;">Left: convergence of x₀ with the CG (dashed lines) and BFGS (continuous lines) algorithms. Right:  dynamics of the true solution (continuous lines) and the numerical obtained from (x₀,θ) computed with  CG (dased line) and BFGS (dash-pointed line) algorithms, and noisy data for observable variables (I<sub>i</sub>,R<sub>i</sub>,D<sub>i</sub>) (points)</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' <span style="text-align: center; font-size: 75%;">Left: convergence of <math>\mathbf{x}_0</math> with the CG (dashed lines) and BFGS (continuous lines) algorithms. Right:  dynamics of the true solution (continuous lines) and the numerical obtained from <math>(\mathbf{x}_0,\boldsymbol{\theta })</math> computed with  CG (dased line) and BFGS (dash-pointed line) algorithms, and noisy data for observable variables <math>(I_i,R_i,D_i)</math> (points)</span>
|}

<span id='theorem-tabla3'></span>Example 4: This last example includes numerical results obtained with the BFGS algorithm only, but with perturbations <math display="inline">noise\_level = 0.05, 0.2</math> for the generation of synthetic noisy measurements. Table [[#theorem-tabla3|4]] summarizes the numerical results.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |Table. 3 <span style="text-align: center; font-size: 75%;">Effect of the noisy data on the numerical results for <math>\boldsymbol{\theta }</math> and <math>\mathbf{x}_0</math> computed with the BFGS algorithm.</span>
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" | <math display="inline">noise\_level</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.05\ (5%)</math>
| style="border-left: 2px solid;" | <math>0.15\ (15%)</math>
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" | <math display="inline">\boldsymbol{\theta }^0</math>                   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(1.40, 0.09,0.20,0.001)</math>
| style="border-left: 2px solid;" | <math>(1.40, 0.09,0.20,0.001)</math>
|-
| style="border-right: 2px solid;" | <math>\mathbf{x}_0^0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(91645,4782,1844,1605,124)</math>
| style="border-left: 2px solid;" | <math>(90474,5782,1628,987,129)</math>
|-
| style="border-right: 2px solid;" | Data time window         
| style="border-left: 2px solid;border-right: 2px solid;" | <math>t_i = 40 +i</math>, <math display="inline">1\le i\le 10</math>      
| style="border-left: 2px solid;" | <math>t_i = 40 + 2i</math>, <math display="inline">1\le i \le 8</math>  
|-
| style="border-right: 2px solid;" | <math>\epsilon </math>,  no. iters.  
| style="border-left: 2px solid;border-right: 2px solid;" | <math>10^{-8}</math>, 30 
| style="border-left: 2px solid;" | <math>10^{-8}</math>, 43   
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |  Computed <math display="inline">\boldsymbol{\theta }</math>      
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(1.0256,0.1356,0.2055,0.0143)</math>
| style="border-left: 2px solid;" | <math>(0.9953,0.1385,0.2126,0.0153)</math>
|-
| style="border-right: 2px solid;" | Relative error                
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(0.0256,0.0508,0.0277,0.0004)</math>
| style="border-left: 2px solid;" | <math>(0.0047,0.0307,0.0628,0.0705)</math>
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |  Computed <math display="inline">\mathbf{x}_0</math>   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(91546,4795,1886,1646,127)</math>
| style="border-left: 2px solid;" | <math>(90759,5780,1637,1690,135)</math>
|- style="border-bottom: 2px solid;"
| style="border-right: 2px solid;" | Relative error                
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(0.0011,0.0119,0.0746,0.0160,0.0194)</math>
| style="border-left: 2px solid;" | <math>(0.0097,0.1910,0.0675,0.0433,0.0770)</math>

|}

Comparing these results with those for BFGS in Table [[#table-2|2]] we observe that the speed of convergence decreases with increasing <math display="inline">noise\_level</math> for the same tolerance. Also, since the initial guess <math display="inline">\mathbf{x}_0^0</math> is closer to exact <math display="inline">\mathbf{x}_0</math> for the 5% noisy case, then the accuracy of the computed value is better in this case. This sensitivity is not so evident in the calculation of the parameters <math display="inline">\boldsymbol{\theta }</math>, since the achieved accuracy is comparable. Another feature is that the time window of noisy data is wider the higher the <math display="inline">noise\_level</math> to achieve convergent results. Figures [[#img-6|6]] and [[#img-7|7]] illustrate the performance of the BFGS algorithm with respect to <math display="inline">noise\_level</math>. <div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-g_05_15_BFGS.png|210px|]]
|[[Image:Draft_Juarez Valencia_878068377-p_05_15_BFGS.png|210px|<span style="text-align: center; font-size: 75%;">Left: plot of (∇L(x₀,θ)) against iteration obtained with the BFGS algorithm for 5% noisy data (blue line) and 15% noisy data (red line). Right: convergence of θ for 5% noisy data (continuous lines) and 15% noisy data (dashed lines).</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 6:''' <span style="text-align: center; font-size: 75%;">Left: plot of <math>\log (\nabla L(\mathbf{x}_0,\boldsymbol{\theta }))</math> against iteration obtained with the BFGS algorithm for 5% noisy data (blue line) and 15% noisy data (red line). Right: convergence of <math>\boldsymbol{\theta }</math> for 5% noisy data (continuous lines) and 15% noisy data (dashed lines).</span>
|}
<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-x0_05_15_BFGS.png|240px|]]
|[[Image:Draft_Juarez Valencia_878068377-x_05_15_BFGS.png|240px|<span style="text-align: center; font-size: 75%;">Left: convergence of x₀ for 5% noisy data (continuous lines) and 15% noisy data (dashed lines). Right:  dynamics of the true solution (continuous lines) and of the numerical solution obtained with 5% noisy data (dashed lines) and 15% noisy data (dash-pointed lines). Here we only show the points with 15% noisy data (see Table [[#table-2|2]]).</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 7:''' <span style="text-align: center; font-size: 75%;">Left: convergence of <math>\mathbf{x}_0</math> for 5% noisy data (continuous lines) and 15% noisy data (dashed lines). Right:  dynamics of the true solution (continuous lines) and of the numerical solution obtained with 5% noisy data (dashed lines) and 15% noisy data (dash-pointed lines). Here we only show the points with 15% noisy data (see Table [[#table-2|2]]).</span>
|}
Another way to enhance convergence of the optimization algorithms is adding observable variables. For instance, adding <math display="inline">E</math> as observable variable for the case of 15% noise and using the same numerical parameters in Table [[#theorem-tabla3|4]], the BFGS algorithm converges in 27 iterations for the given tolerance. The best improvement, besides the faster convergence, is the estimation of the initial conditions, as shown in Table [[#table-4|4]]


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-4'></span>Table. 4 <span style="text-align: center; font-size: 75%;">Numerical results obtained with the BFGS algorithm, adding <math>E</math> as observable variable.</span>
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |  Parameter    
| style="border-left: 2px solid;border-right: 2px solid;" | Computed                   
| style="border-left: 2px solid;" | Relative error                 
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" | <math display="inline">\boldsymbol{\theta }</math>      
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(1.0106,0.1415,0.2131,0.0153)</math>
| style="border-left: 2px solid;" | <math>(0.0106,0.0092,0.0657,0.0731)</math>
|- style="border-bottom: 2px solid;"
| style="border-right: 2px solid;" | <math>\mathbf{x}_0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>(91154,5381,1632,1698,135)</math>
| style="border-left: 2px solid;" | <math>(0.0054,0.1088,0.0704,0.0480,0.0821)</math>

|}

Figures [[#img-6|6]] and [[#img-7|7]] show the corresponding improvements. <div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-g_15_EIRD.png|210px|]]
|[[Image:Draft_Juarez Valencia_878068377-p_15_EIRD.png|210px|<span style="text-align: center; font-size: 75%;">Left: plot of (∇L(x₀,θ)) against iteration obtained with the BFGS algorithm for 15% noisy data and observable variables (E,I,R,D). Right: convergence of θ.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 8:''' <span style="text-align: center; font-size: 75%;">Left: plot of <math>\log (\nabla L(\mathbf{x}_0,\boldsymbol{\theta }))</math> against iteration obtained with the BFGS algorithm for 15% noisy data and observable variables <math>(E,I,R,D)</math>. Right: convergence of <math>\boldsymbol{\theta }</math>.</span>
|}
<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Juarez Valencia_878068377-x0_15_EIRD.png|240px|]]
|[[Image:Draft_Juarez Valencia_878068377-x_15_EIRD.png|240px|<span style="text-align: center; font-size: 75%;">Left: convergence of x₀ for 15% noisy data and observable variables (E,I,R,D). Right: dynamics of the true solution (continuous lines) and of the numerical solution (dashed lines).</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 9:''' <span style="text-align: center; font-size: 75%;">Left: convergence of <math>\mathbf{x}_0</math> for 15% noisy data and observable variables <math>(E,I,R,D)</math>. Right: dynamics of the true solution (continuous lines) and of the numerical solution (dashed lines).</span>
|}

==6 Conclusions==

We have introduced a deterministic model for fitting observed noisy data into a given dynamical system to find initial conditions and the parameters of the associated system of ordinary differential equations. The classical CG and BFGS optimization algorithms are employed to minimize the quadratic non-linear cost function. It is shown the advantage of using the adjoint equation approach to find the derivatives or gradients. We explain with some detail the implementation of this methods and algorithms with the SEIRD epidemiological model. However, this approach can be equally applied to other problems modelled by ODEs.

Similar numerical results are obtained with both algorithms, CG and BFGS using the same tolerance to achieve a given accuracy, but as expected the BFGS algorithm has better convergence properties and it is more robust. Numerical results show that more experimental data points and more observable variables increase the convergence properties of these algorithms. On the other hand, the higher the noise of the experimental data the slower is the convergence of the optimization algorithms. The main drawback of the proposed methodology is that it is sensitive to the location of noisy data and also to the initial guesses for initial conditions. However, if the algorithms converge properly, then the numerical results obtained are more accurate when <math display="inline">\mathbf{x}_0</math> is also estimated along with <math display="inline">\boldsymbol{\theta }</math>.

For future work, we want to overcome some difficulties or deficiencies that arise with the proposed model and numerical algorithms. First, we must include explicitly into the fitting model ([[#eq-6|6]]) the positivity constraint of the unknown parameters, <math display="inline">\mathbf{x}_0</math> and <math display="inline">\boldsymbol{\theta }</math>, specially for those that are relatively small and extend the proposed algorithms accordingly, like in <span id='citeF-7'></span>[[#cite-7|[7]]]. The inherent instability and difficulty to find the initial conditions may be fixed incorporating the technique of multiple shooting, e.g. <span id='citeF-5'></span>[[#cite-5|[5]]], <span id='citeF-11'></span>[[#cite-11|[11]]]. Concerning the efficiency of optimization algorithms, we still need to test the Gauss-Newton method and if necessary its variant, the Levenberg&#8211;Marquardt algorithm. Finally, as mentioned before, the line search strategy is crucial for gradient descent algorithms. We may improve the performance of the secant method incorporating bracketing strategies like in <span id='citeF-1'></span>[[#cite-1|[1]]], or trying the Newton's method as mentioned in remark [[#theorem-Newton_line_search|7]].

==Acknowledgements==

We want to acknowledge the Department of Mathematics at Universidad Autónoma Metroprolitana &#8211; Izatapalapa and to CONACyT for the support for this research work.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  I. K. Argyros, M. A. Hernández-Verón, M. J. Rubio, M. J., ''On the Convergence of Secant-Like Methods'', in Current Trends in Mathematical Analysis and Its Interdisciplinary Applications (2019) 141&#8211;183.

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  J. R. Banga, C. G. Moles, A. A. Alonso, ''Global optimization of bioprocesses using stochastic and hybrid methods'', in: Frontiers in global optimization, Springer (2004) pp. 45&#8211;70.

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''  J. Calver and W. Enright, ''Numerical methods for computing sensitivities for ODEs and DDEs'', Numerical Algorithms 74(4) (2017) 1101&#8211;1117.

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  Y. Cao, S. Li,L. Petzold, R. Serban, ''Adjoint sensitivity analysis for differential algebraic equations: the adjoint DAE system and its numerical solution''. SIAM J. Sci. Comput. 3(24) (2003), 1076&#8211;1089.

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''   F. Carbonell, Y. Iturria-Medina, J.C. Jimenez, ''Multiple Shooting-Local Linearization method for the identification of dynamical systems'', Communications in Nonlinear Science and Numerical Simulation, 37 C (2016) 292&#8211;304.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  Calderhead, B., Girolami, M. ''Estimating Bayes factors via thermodynamic integration and population MCMC'' Comput. Stat. Data Anal. 53 (12), (2009), 4028&#8211;4045.

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  M. Victoria Chávez, L. Héctor Juárez, Yasmín A. Ríos, ''Penalization and augmented Lagrangian for OD demand matrix estimation from transit segment counts'', Transportmetrica A, Transport Science, 15(2) (2019), 915&#8211;943.

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  J. E. Dennis and Robert B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Englewood Cliffs: Prentice-Hall. (1983).

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  Jorge Nocedal, Stephen J. Wright, Numerical Optimization, New York: Springer, 1999.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  I. F. D. Oliveira, R. H. C. Takahashi, An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality, ACM Transactions on Mathematical Software. 47(1) (2020) 5:1&#8211;5:24.

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  Ozgur Aydogmus, Ali Hakan TOR, ''A Modified Multiple Shooting Algorithm for Parameter Estimation in ODEs Using Adjoint Sensitivity Analysis'', Applied Mathematics and Computation Volume 390(1), (2021) 125644.

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''  Elena L. Piccolomini and Fabiana Zama, ''Monitoring Italian COVID-19 spread by an adaptive SEIRD model'', medRxiv preprint doi: https://doi.org/10.1101/2020.04.03.20049734, April 6, 2020.

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  Ramsay, J., Hooker, H., Campbell, D., Cao, J.  ''Parameter estimation for differential equations: a generalized smoothing approach''. J. R. Stat. Soc. Ser. B 69 (5), (2007) 741&#8211;796.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  B. Sengupta, K.J. Friston, W.D. Penny, ''Efficient gradient computation for dynamical models'', NeuroImage 98 (2014) 521&#8211;527.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  Wenyu Sun, Ya-Xiang Yuan, Optimization Theory and Methods:Nonlinear Programming. New York: Springer, 2006.