Inverse problem in linear ODE

Question

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?

Answer 1

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.

A good book to read about this is:

Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems

Author(s): J. Nathan Kutz , Steven L. Brunton , Bingni W. Brunton and Joshua L. Proctor