MInimizing cost function using iterative search for a minimum method

Question

I want to estimated the parameters $\ \hat{\theta} $ of a model using an iterative search for the minimum of a cost function. The cost function is defined as follows:

$$ V_N(\hat{\theta}) = \frac{1}{N}\sum_{i=1}^N(y(t_k)-\hat{y}(t_k|\theta))^2 $$

where $\ y $ is the output of the system and $\ \hat{y} $ is the estimated output of the system. The system is described by the following differential equations:

$$ \dot{h}(t) = -\theta\sqrt{2g}\sqrt{h(t)}+u(t) $$

$$ \hat{y}(t|\theta) = \theta\sqrt{2g}\sqrt{h(t)} $$

where $\ u $ is the input to the system. Suppose that data for both input and output of the system are collected and available. The equation for updating the estimate for the unknown parameters $\ \theta $ is:

$$ \hat{\theta}_{i+1} = \hat{\theta}_i-μ_i[\frac{d^2}{d\theta^2}(V_N(\hat{\theta}))]^{-1}\frac{d}{d\theta}(V_N(\hat{\theta})) $$ where $\ μ_i $ is a step length determined so that: $\ V_N(\hat{\theta}_{i+1}) < V_N(\hat{\theta}_i) $. The derivatives of the cost function are:

$$ \frac{d^2}{d\theta^2}(V_N(\hat{\theta})) = \frac{1}{N}\sum_{i=1}^N(\frac{d}{d\hat{\theta}}\hat{y}(t_i|\hat{\theta}))((\frac{d}{d\hat{\theta}}\hat{y}(t_i|\hat{\theta}))^T-\frac{1}{N}\sum_{i=1}^N\frac{d^2}{d\hat{\theta}^2}(\hat{y}(t_i|\hat{\theta}))(y(t_i)-\hat{y}(t_i|\hat{\theta})) $$

which by neglecting the second sum comes down to the Gauss-Newton Method. Considering this the whole problem is solved by finding the way to compute:

$$ ψ(t,\hat{\theta})=\frac{d}{d\hat{\theta}}\hat{y}(t|\hat{\theta}) $$

By working out the math, the following differential equations are obtained:

$$ z(t,\hat{\theta}) = \frac{d}{d\hat{\theta}}x(t,\hat{\theta}) $$ $$ \frac{d}{dt}z(t,\hat{\theta}) = -\frac{\hat{\theta}\sqrt{2g}}{2\sqrt{x(t,\hat{\theta})}}z(t,\hat{\theta})-\sqrt{2gx(t,\hat{\theta})} $$

$$ ψ^T(t,\hat{\theta}) = \frac{d}{d\hat{\theta}}\hat{y}(t,\hat{\theta}) = \frac{\hat{\theta}\sqrt{2g}}{2\sqrt{x(t,\hat{\theta})}}z(t,\hat{\theta})+\sqrt{2gx(t,\hat{\theta})} $$

In order to calculate $\ \frac{d}{d\hat{\theta}}{\hat{y}(t,\hat{\theta})} $ we need to first compute $\ x(t,\hat{\theta}) $ and $\ z(t,\hat{\theta}) $ since $\ g $ is the gravity constant. Suppose that we have an initial guess for the value of $\ \hat{\theta} $, now $\ x(t,\hat{\theta}) $ is obtained by the differential equation $\ \dot{x}(t,\hat{\theta}) $ since we also have the value of the input data $\ u $ and some initial condition for $\ x(t,\hat{\theta}) $. My question is how to compute $\ z(t,\hat{\theta}) $ since there are two equations which give as output $\ z(t,\hat{\theta}) $ ?

Is the sequence in which the equations should be computed the following:

1. Compute $\ x(t,\theta) $
2. Compute $\ z(t,\theta) $
3. Compute $\ \frac{d}{d\theta}\hat{y}(t,\theta) $
4. Compute $\ \frac{d^2}{d\theta^2}V_N(\theta $
5. Compute $\ \frac{d}{d\theta}V_N(\theta) $
6. Update estimation of $\ \hat{\theta} $

Would an ODE solver of MATLAB do the work in one bunch ?

Answer 1

I don't see how two equations give $z$ as output. Nevertheless, your sequence of computations looks reasonable, except I would combine steps one and two into:

1. Simultaneously solve for $x$ and the sensitivities $z$. This is a extended ODE system that you could throw at a built-in MATLAB ODE solver:

$$ \begin{bmatrix} x'(t,\hat{\theta}) \\ z'(t,\hat{\theta}) \end{bmatrix} = \begin{bmatrix} -\theta\sqrt{2g}\sqrt{x(t)}+u(t) \\ -\frac{\hat{\theta}\sqrt{2g}}{2\sqrt{x(t,\hat{\theta})}}z(t,\hat{\theta})-\sqrt{2gx(t,\hat{\theta})} \end{bmatrix} $$

If your error $V$ used a continuous $L_2$ norm instead of a discrete $\ell_2$ norm, the integral and its sensitivities could be computed as the solution to an ODE, extending the system even more. Maybe this would be preferable.

Since you're using MATLAB, you may be interested in the MATLODE package, which provides several types of time integration schemes. Its tangent linear and adjoint model integration can compute sensitivities with respect to model parameters as well as initial conditions. Note its discrete adjoint integrators may be of interest for efficiently computing the application of $\frac{d}{d \widehat{\theta}} \widehat{y}^T$ with a vector as needed in the Hessian approximation.

Answer 2

I have been thinking that it might be easier if one changes first the variables in the differential equation. That way one can bypass the function $h(t)$ and deal with fewer functions. Since $$\hat{y}(t) = \hat{y}_{\theta}(t) = \hat{y}(t \,| \, \theta) = \theta \, \sqrt{2g} \, \sqrt{h(t)}$$ change the dependent variable $$\hat{y} = \theta \, \sqrt{2g} \, \sqrt{h}$$ whose converse is $$h = \frac{1}{2g \,\theta^2} \,\, \hat{y}^2$$ Consequently \begin{align} \frac{dh}{dt} = \frac{1}{g \,\theta^2} \,\,\hat{y} \, \frac{d\hat{y}}{dt} = -\,\theta \, \sqrt{2g}\, \sqrt{h} \, +\, u(t) = -\, \hat{y} \, + \, u(t) \end{align} so the differential equation after change of variable becomes $$\frac{d\hat{y}}{dt} = g \,\theta^2 \, \frac{u(t)}{\hat{y}} \, - \, g \,\theta^2 \,$$ The first variation of $\hat{y}$ with respect to $\theta$, i.e. the first derivatives of $\hat{y}$ with respect to $\theta$, satisfy the system of linear equations $$\frac{d}{dt}\,\partial_{\theta}\hat{y} = - \, g \,\theta^2 \, \frac{u(t)}{\hat{y}^2}\, \partial_{\theta}\hat{y} \, + \,2\, g \,\theta \, \frac{u(t)}{\hat{y}} \, - \,2\, g \,\theta $$ Thus, as a whole, the functions $\big(\,\hat{y}(t \, | \, \theta),\,\,\, \partial_{\theta}\hat{y}(t \, | \, \theta)\, \big)$ satisfy the system of ODEs \begin{align} &\frac{d\hat{y}}{dt} = g \,\theta^2 \, \frac{u(t)}{\hat{y}} \, - \, g \,\theta^2 \\ &\frac{d}{dt}\,\partial_{\theta}\hat{y} = - \, g \,\theta^2 \, \frac{u(t)}{\hat{y}^2}\, \partial_{\theta}\hat{y} \, + \, 2\,g \,\theta \, \frac{u(t)}{\hat{y}} \, - \,2\, g \,\theta \end{align} The gradient of your functional $$V_N(\theta) =\frac{1}{N}\sum_{k=1}^{N} \Big(\,\hat{y}(t_k \, | \, \theta) \, - \,y(t_k) \,\Big)^2 $$ should be $$\partial_{\theta}V_N(\theta) =\frac{2}{N}\sum_{k=1}^{N} \Big(\,\hat{y}(t_k \, | \, \theta) \, - \,y(t_k) \,\Big)\, \partial_{\theta}\hat{y}(t_k \, | \, \theta)$$ I do not know how important this Newton's method is to you (it is probably faster than what I am going to propose), but if you do not feel like taking second derivatives of $\hat{y}$ with respect to $\theta$, you can try first the slower gradient descend method: $$\theta_{i+1} = \theta_{i} - \mu_i \, \partial_{\theta}V_N(\theta_i)$$ Then the parameter fitting algorithm should be like this:

You start with given data points $y(t_1), \, y(t_2), \, ..., y(t_N)$. Take a reasonable guess $\theta_1$ for the parameter. Assume that at step $n$, you have calculated an approximate parameter $\theta_n$:

Step 1: Use an ODE solver, like Runge-Kutta or something like that, to generate solutions of the initial value problem \begin{align} &\frac{d\hat{y}}{dt} = g \,\theta_n^2 \, \frac{u(t)}{\hat{y}} \, - \, g \,\theta_n^2 \\ &\frac{d}{dt}\,\partial_{\theta}\hat{y} = - \, g \,\theta_n^2 \, \frac{u(t)}{\hat{y}^2}\, \partial_{\theta}\hat{y} \, + \,2\, g \,\theta \, \frac{u(t)}{\hat{y}} \, - \,2\, g \,\theta \\ &\hat{y}(t_1) = y(t_1)\\ &\partial_{\theta}\hat{y}(t_1) = 0 \end{align} and obtain the sequence of solution data $$\Big\{ \, \big(\,\hat{y}(t_k \, | \, \theta),\,\, \partial_{\theta}\hat{y}(t_k \, | \, \theta)\, \big) \, \,\, \big{|} \,\,\, k = 1, 2, 3, ..., N \,\, \Big\}$$

Step 2: The sequence of solution data from step 1 allows us to calculate the cost function $$V_N(\theta) =\frac{1}{N}\sum_{k=1}^{N} \Big(\,\hat{y}(t_k \, | \, \theta) \, - \,y(t_k) \,\Big)^2$$

Step 3: If $V_{N}(\theta_n) < \varepsilon$, stop the algorithm and take as a solution $\theta = \theta_n$. Otherwise, if the error is $V_N(\theta_n) \geq \varepsilon$, calculate the gradient of $V_N$, using the solution data from step 1: $$\partial_{\theta}V_N(\theta_n) =\frac{2}{N}\sum_{k=1}^{N} \Big(\,\hat{y}(t_k \, | \, \theta_n) \, - \,y(t_k) \,\Big)\, \partial_{\theta}\hat{y}(t_k \, | \, \theta_n)$$

Step 4: Using the gradient calculated in step 3, update the $\theta$ parameter: $$\theta_{n+1} = \theta_{n} - \mu_n \, \partial_{\theta}V_N(\theta_n)$$

Step 5: Go back to step 1 with the newly calculated parameters $\theta_{n+1}$.

Even if you really want the Newton's method, you can add the extra differential equations for the second derivatives of $\hat{y}(t \, | \, \theta)$ with respect to $\theta$ in step 1 and add the calculation of the hessian of the cost function $V_{N}(\theta)$ in step 3.