# Inverse problem in linear ODE

I have a linear ordinary differential equation (ODE) with a system matrix with constant coefficients: $$\dot{y}(t) = \mathcal{A}\; y(t), \quad y(0) = y_0$$ with $y(t) \in \mathbb{R}^{n \times 1}$ and $\mathcal{A} \in \mathbb{R}^{n \times n}$ where $\mathcal{A}$ has constant coefficients. I have a (large) series of $r$ (with $r \gg n^2$) datapoints $y_{i}$ for time points $t_i$ (including $y_0$ for $t=t_0$). The elements of the matrix $\mathcal{A}$ are unknown and I'm interested in estimating these.

In general, this is considered as an inverse problem. Is there a more specific terminology for this example (constant coefficients, linear ODE, system matrix identifcation)?

What are the current techniques to solve this type of problem? How sensitive is the problem to noisy data?

My "brute force" approach would be

1. Write the elements of matrix $\mathcal{A}$ in a vector $a = [a_{11}, a_{12}, \ldots, a_{1n}, a_{21}, \ldots, a_{2n}, \ldots, a_{nn}]$
2. Write down the equations for each $t_i$ time point using some approximation for $\dot{y}(t) \approx \frac{y_{i} - y_{i-1}}{\Delta t}$
3. Build a linear least squares problem using this systems matrix and the vector of unknowns.
4. I can perform an SVD decomposition $\mathcal{A} = U\Sigma V^{t}$
5. Since the problem is homogeneous, the LS optimum would be the last column of $V$.

But the problem is of course, that if the number of data points gets larger, the computation of the full SVD becomes memory and time consuming.

Given the simplicity of the problem formulation (leading to a very structured matrix $\mathcal{A}$), I would be surprised if there is no specific, more optimal, method developed for this type of problem. Would an approach with the formal solution using the matrix exponential $y(t) = e^{\mathcal{A}t}y_{0}$ lead to a more optimal algorithm?

• In frequency domain this is $Y(i\omega) = (i\omega I -A)^{-1}Y(0)$. In control theory circles this is called system identification. Nov 22, 2017 at 11:54
• Note that you want $r > n^{2}$, not just $n$ because there are $n^{2}$ parameters to identify. Nov 22, 2017 at 12:37
• Note that this is not a linear inverse problem even if the ODE is linear, since the mapping $a\mapsto y$ is not. (This also means that $r>n^2$ is in general neither necessary nor sufficent.) Another keyword would be "output least squares" (a.k.a. Tikhonov regularization) together with "adjoint approach". Nov 22, 2017 at 13:22
• As @percusse already mentioned, this is (a simple case of) linear black-box system identification, about which there are of course dozens of textbooks available. You might be interested in, e.g. chapter 7 of the book of Söderström and Stoica: user.it.uu.se/%7Ets/sysidbook.pdf. Feb 5, 2019 at 19:32
• On another note, where does the "very structured matrix A" come from? I see no corresponding assumption being made above? If you assume some specific strucure (e.g. a particular transfer function model --> see the Söderström/Stoica book again), the identification problem becomes of course much easier. Feb 5, 2019 at 19:34

Consider the vector $$y_i$$ for a given time $$t_i$$. We are looking for the linear operator $$\mathcal{B}$$ that maps $$y_i$$ to $$y_{i+1}$$ by noting that

\begin{align} \dot{y}&= \mathcal{A} y\\ \frac{y_{i+1}-y_i}{\Delta t}&=\mathcal{A} y_i\\ \Rightarrow y_{i+1}&=(\mathcal{A}\Delta t+I) y_i\\ y_{i+1}&=\mathcal{B} y_i \end{align}

For a set of vectors $$y_i$$, we can construct the matrices $$\mathcal{Y}$$ and $$\mathcal{Y}'$$. \begin{align} \mathcal{Y}= \begin{bmatrix} \vert & \vert & & \vert \\ y_1 & y_2 & \dots & y_{n-1} \\ \vert & \vert & & \vert \end{bmatrix} \\ \mathcal{Y}'= \begin{bmatrix} \vert & \vert & & \vert \\ y_2 & y_3 & \dots & y_n \\ \vert & \vert & & \vert \end{bmatrix} \end{align} We can therefore write: \begin{align} \mathcal{Y}'=\mathcal{B}\mathcal{Y} \end{align}

The optimization then involves solving:

\begin{align} \Phi(\mathcal{B})=\min_{\mathcal{B}}\vert\vert\mathcal{B}\mathcal{Y}-\mathcal{Y}'\vert\vert_F \end{align}

Here, $$\vert\vert\cdot\vert\vert_F$$ denotes the Frobenius norm.

We can, as you said, minimize the objective function via an SVD, by seeing that

\begin{align} \mathcal{B}\mathcal{Y}&=\mathcal{Y}'\\ \mathcal{B}\;{U_k^{\,}S_k^{\,}V_k^{\,}}^T&=\mathcal{Y}'\\ \Rightarrow \mathcal{B}&=\mathcal{Y}'{V_k^{\,}S_k^{-1}U_k^{\,}}^T \end{align}

Here, $$k$$ denotes the number of singular vectors/values used for the construction of the inverse operator. You could e.g. use Matlabs svds function for the subset of the first $$k$$ singular vectors/values.

Alternatively, we can iterate the solution via gradient descent. The derivative of the with respect to $$\mathcal{B}$$ of the objective function is:

\begin{align} \frac{\partial \Phi}{\partial \mathcal{B}}=2(\mathcal{B}\mathcal{Y}-\mathcal{Y}')\mathcal{Y}^T \end{align}

Here, we can preallocate the products $$\mathcal{Y}'\mathcal{Y}^T$$ and $$\mathcal{Y}\mathcal{Y}^T$$. This may be faster than waiting for the SVD.

We thus iterate until we hit some stopping criterion:

\begin{align} \mathcal{B}_{k+1} = \mathcal{B}_{k}-\alpha(\mathcal{B}_{k}\mathcal{Y}-\mathcal{Y}')\mathcal{Y}^T \end{align}

With $$\alpha$$ being the step size.

Finally, we can get $$\mathcal{A}$$ by

\begin{align} \mathcal{A} = \frac{\mathcal{B}-I}{\Delta t} \end{align}

To prevent overfitting due to noise, we can limit the size of $$\mathcal{B}$$ via adding a regularizer to the objective function. Often times this is implemented via a proximal mapping onto a ball in a given metric. Then, each iteration is projected back into a feasible set. This is a form of constrained optimization.