# Least square estimation using Euler discretization

I have a set of differential equations $\dot{x}_i = f_i(x, \theta)$, with $x = [x_1, \ldots, x_n]^\top \in \mathbb{R}^n$ and $\theta \in \mathbb{R}^m$. Measurements of the variables $x_i$ are available at certain equispaced time instants, let say $\bar{x}_{i,k} = x_i(t_k)$, with $t_k = t_0 + (k-1)\Delta$ and $\Delta >0$ for all $k = 0, \ldots, T$.

It turns out that $f_i(x, \theta) = m_i(x)^\top \theta$, where $m_i(x)$ is a vector in $\mathbb{R}^m$. The linear dependence of the system on the parameters suggests me that I can try to estimate $\theta$ using least square approach on the basis of the observations $\bar{x}_{i,k}$.

To this aim, I employ Euler discretization:

$$x_{i,k+1} = x_{i,k} + \Delta m_i(x_k)^\top \theta,$$

where $x_{i,k} = x_i(t_k)$ and $x_k = [x_{1,k}, \ldots, x_{n,k}]^\top$.

At this point, I build the following quantities:

$$Y(x) = \left[ \begin{array}{c} x_{1,1} - x_{1,0} \\ x_{2,1} - x_{2,0}\\ \vdots \\ x_{n,1} - x_{n,0}\\ x_{1,2} - x_{1,1}\\ x_{2,2} - x_{2,1}\\ \vdots \\ x_{n,2} - x_{n,1}\\ \vdots \\ x_{1,T} - x_{1,T-1}\\ x_{2,T} - x_{2,T-1}\\ \vdots \\ x_{n,T} - x_{n,T}\\ \end{array} \right] \text{and}~ U(x) = \Delta\left[ \begin{array}{c} m_1(x_0)^\top \\ m_2(x_0)^\top\\ \vdots \\ m_n(x_0)^\top\\ m_1(x_1)^\top\\ m_2(x_1)^\top\\ \vdots \\ m_n(x_1)^\top\\ \vdots \\ m_1(x_{T-1})^\top\\ m_2(x_{T-1})^\top\\ \vdots \\ m_n(x_{T-1})^\top\\ \end{array} \right],$$

in order to write the whole stuff as:

$$Y(x) = U(x) \theta.$$

The least square estimation of $\theta$ on the basis of observations $\hat{x}$ is:

$$\hat\theta = (U(\hat{x})^\top U(\hat{x}))^{-1}U(\hat{x})^\top Y(\hat{x}).$$

Here comes my question.

The objective function I'm minimizing is:

$$\frac{1}{n(T-1)}\sum_{i=1}^n\sum_{k=0}^{T-1} \left(\hat{x}_{i,k+1}- \hat{x}_{i,k}- \Delta m_i(\hat{x}_k)^\top \theta\right)^2.$$

Anyway, this is not so intuitive, because in the "standard" case one has:

$$\frac{1}{n(T-1)}\sum_{i=1}^n\sum_{k=0}^{T-1} \left(\text{DATA}_{i,k} - \text{MODEL}_{i,k} \cdot \text{PARAMETERS} \right)^2.$$

Is there another way to write down my problem as a "standard" one? Moreover, are there numerical advantages to write the problem in some particular form?

• Since you have a dynamical system for which measurements arrive over time it seems to me more appropriate to use a recursive least squares algorithm. Other than that, I do not see any mismatch between your objective function and the "standard" one. – RTJ Dec 15 '15 at 12:15
• @CTNT thanks a lot for the hints. At this point, why not a Kalman Filter? – the_candyman Dec 15 '15 at 17:15
• Yes, you can also use a Kalman filter. – RTJ Dec 15 '15 at 17:32
• @CTNT in your opinion, which approach is better? RLS or KF? – the_candyman Dec 15 '15 at 17:37
• If you know the measurement noise covariance matrix then go for the Kalman filter on the parameters(i.e. considering dynamics $\dot{\theta}=0$). Otherwise I would probably use RLS. – RTJ Dec 15 '15 at 17:42