# Identifiability Analysis of Linear Ordinary Differential Equation Systems
with a Single Trajectory111This work was funded in part by the National
Institutes of Health under award number R01 AI087135 and the University of
Rochester CTSA award number UL1 TR002001 from the National Center for
Advancing Translational Sciences of the National Institutes of Health.

Xing Qiu<EMAIL_ADDRESS>Tao Xu<EMAIL_ADDRESS>Babak
Soltanalizadeh<EMAIL_ADDRESS>Hulin Wu
<EMAIL_ADDRESS>Department of Biostatistics and Computational Biology
University of Rochester, Rochester, NY, U.S.A.  Department of Biostatistics &
Data Science, School of Public Health
University of Texas Health Science Center at Houston, Houston, TX, U.S.A.

###### Abstract

Ordinary differential equations (ODEs) are widely used to model dynamical
behavior of systems. It is important to perform identifiability analysis prior
to estimating unknown parameters in ODEs (a.k.a. inverse problem), because if
a system is unidentifiable, the estimation procedure may fail or produce
erroneous and misleading results.

Although several qualitative identifiability measures have been proposed, much
less effort has been given to developing _quantitative_ (continuous) scores
that are robust to uncertainties in the data, especially for those cases in
which the data are presented as a single trajectory beginning with one initial
value.

In this paper, we first derived a closed-form representation of linear ODE
systems that are not identifiable based on a single trajectory. This
representation helps researchers design practical systems and choose the right
prior structural information in practice. Next, we proposed several
quantitative scores for identifiability analysis in practice. In simulation
studies, the proposed measures outperformed the main competing method
significantly, especially when noise was presented in the data. We also
discussed the asymptotic properties of practical identifiability for high-
dimensional ODE systems and conclude that, without additional prior
information, many random ODE systems are practically unidentifiable when the
dimension approaches infinity.

###### keywords:

linear ordinary differential equations , structural identifiability ,
practical identifiability , inverse problem , parameter estimation

††journal: Applied Mathematics and Computation

##  1 Background and Introduction

Ordinary differential equations (ODE) can be used to model complex dynamic
systems in a wide variety of disciplines including economics, physics,
engineering, chemistry, and biology [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. Such
systems are usually represented as

$\displaystyle\begin{cases}D\mathbf{x}(t)=f(\mathbf{x}(t),\mathbf{u}(t),\theta),&t\in(0,T),\\\
\mathbf{x}(0)=\mathbf{x}_{0}.\end{cases}$ (1)
$\displaystyle\mathbf{y}(t)=h(\mathbf{x}(t),\mathbf{u}(t),\theta).$ (2)

Here $\mathbf{x}(t)=(x_{1}(t),\dots,x_{d}(t))^{\prime}\in\mathbb{R}^{d}$ is
the state vector, $D=\frac{\mathrm{d}}{\mathrm{d}t}$ is the first order
differential operator222To avoid confusion, we reserve symbol ′ (apostrophe)
for matrix transpose, not the derivative with respect to $t$.,
$\mathbf{y}(t)\in R^{d}$ is the output vector, $\mathbf{u}(t)$ is a known
system input vector, $f$, $g$ are known families of linear or nonlinear
functions indexed by $\theta\in\mathbb{R}^{p}$, which is the vector of unknown
parameters to be estimated. Equation (1) is called the state equation and
Equation (2) the output or observation equation. In this paper, we focus on an
important special case of the above general ODE system: homogeneous linear ODE
system with complete observation

$D\mathbf{x}(t)=A\mathbf{x}(t),\qquad\mathbf{x}(0)=\mathbf{x}_{0}=(x_{1,0},\dots,x_{d,0})^{\prime},\qquad\mathbf{y}(t)=\mathbf{x}(t).$
(3)

Here $A\in M_{d\times d}$ is a matrix (called the system matrix) that
characterizes the mechanistic relationship between $x_{i}(t)$;
$\mathbf{y}(t)=\mathbf{x}(t)$ means that we can directly observe
$\mathbf{x}(t)$ in all dimensions.

Let $\mathbf{x}(t|A,\mathbf{x}_{0})$ be the solution curve (a.k.a. the
trajectories) of Equation (3) initiated at $\mathbf{x}_{0}$ and governed by
system matrix $A$. It is well known that $\mathbf{x}(t|A,\mathbf{x}_{0})$ can
be represented as a unique matrix exponential

$\mathbf{x}(t|A,\mathbf{x}_{0})=e^{tA}\mathbf{x}_{0}.$ (4)

As such, the _forward problem_ of Equation (3), defined as solving the ODE
system with given $A$ and $\mathbf{x}_{0}$, has been resolved in the
mathematical sense – despite of several known numerical issues in matrix
exponentials for high-dimensional data [11].

In practical applications, the parameters that characterize the ODE system,
such as $A$ and $\mathbf{x}_{0}$ in Equation (3), must be estimated from the
real data. This is known as the _inverse problem_. Over the years, many
parameter estimation methods have been developed for ODE systems [2, 6, 12,
13, 7, 14, 15, 16, 17].

In principle, before performing the parameter estimation, we need to address
an important question: are the parameters in a particular ODE model
_identifiable_ from the data? In this context, “identifiability” loosely means
that there is a _unique_ mapping between the trajectories and the parameters
of a family of ODE systems.

By now, a rich literature on the identifiability of both linear and nonlinear
ODE systems is available, see [18] for a thorough review of these methods.
Unfortunately, most of them only consider the identifiability of the system
matrix ($A$), and assume that one can choose an arbitrary initial condition
$\mathbf{x}_{0}\in\mathbb{R}^{d}$. For example, the _global identifiability_
used in some literature on nonlinear ODE identifiability (e.g., [19]) reduces
to the following definition for Equation (3), as pointed out in [20]:

######  Definition 1.1.

Linear ODE system (3) is globally identifiable in a subset $\Omega\subset
M_{d\times d}$ iff for all $A,B\in\Omega$, $A\neq B$, there exists
$\mathbf{x}_{0}\in\mathbb{R}^{d}$, such that
$\mathbf{x}(t|A,\mathbf{x}_{0})\neq\mathbf{x}(t|B,\mathbf{x}_{0})$.

However, such definition is of little use for linear ODE systems because
Stanhope and colleagues proved in [20] that, due to the linearity of Model
(3), such ODE system is _always globally identifiable_ in the entire parameter
space $M_{d\times d}$. Consequently, no computationally intensive symbolic
computation on global identifiability is needed for linear ODE systems. The
above definition of identifiability is also impractical because in many real
world applications, data are only available in the form of one trajectory
starting with a single $\mathbf{x}_{0}$. For example, influenza infection
affects the state of transcriptome of a patient, which can be modeled by
Equation (3) ([21, 22, 23, 24]). However, it is currently impossible for a
researcher to select an arbitrary $\mathbf{x}_{0}$ even in an animal study,
because not only we do not have the technology to alter whole transcriptome
globally, but also not all transcriptome states are biologically feasible.
Furthermore, we cannot repeat the same $\mathbf{x}_{0}$ for a subject either,
because the infection can have long-lasting effects to the immune system of
that subject [25].

As a response to this weakness, several researchers developed a concept known
as locally strong identifiability [26] or $\mathbf{x}_{0}$-identifiability
[27], that involves data with only one trajectory. For Equation (3), it can be
stated as follows.

######  Definition 1.2.

Equation (3) is $\mathbf{x}_{0}$-identifiable w.r.t. a given $\mathbf{x}_{0}$
iff there exists an open and dense subset $\Omega\subset M_{d\times d}$, such
that for all $A,B\in\Omega$, $A\neq B$, we have
$\mathbf{x}(t|A,\mathbf{x}_{0})\neq\mathbf{x}(t|B,\mathbf{x}_{0})$ on
$(0,\delta t)$, for some $0<\delta t<T$.

Of note, the following natural extension to $\mathbf{x}_{0}$-identifiability
was proposed in [27]:

######  Definition 1.3.

System (3) is structurally identifiable iff there exist open and dense subsets
$\Omega\subset M_{d\times d}$, $M^{0}\in\mathbb{R}^{d}$, such that for all
$A,B\in\Omega$, $A\neq B$ and all $\mathbf{x}_{0}\in M^{0}$, we have
$\mathbf{x}(t|A,\mathbf{x}_{0})\neq\mathbf{x}(t|B,\mathbf{x}_{0})$ on
$(0,\delta t)$, for some $0<\delta t<T$.

In other words, Definition 1.3 is $\mathbf{x}_{0}$-identifiability that
applies to not one $\mathbf{x}_{0}$, but an open and dense set
$M^{0}\in\mathbb{R}^{d}$. This definition is consistent with the _structural
identifiability_ [28] and _geometrical identifiability_ [26] for nonlinear ODE
systems.

As a remark, the _open and dense_ condition of $\Omega$ was designed to rule
out a set of certain “inconvenient” parameters that has zero-measure. For
example, it can be shown that if $A$ has repeated eigenvalues, it is not
identifiable with a class of other system matrices (see Section S5 and example
S7 in Supplementary Text for more details). A workaround is to simply define
$\Omega$ to be those matrices with no repeated eigenvalues, which is clearly a
dense open set in $M_{d\times d}$. However, it is theoretically possible that
the said open and dense set $\Omega$ may be “small” compared with $M_{d\times
d}$ in terms of a measure such as $\lambda_{d\times d}$, the Lebesgue measure.
See Example S9 in Supplementary Text, Section S8 for such an example. This can
be seen as a weakness because in most real world applications, there is
uncertainty in $A$, so we want to ensure that the identifiability applies to
_almost every_ $A\in M_{d\times d}$, not just a dense set with small measure
or probability. As a concrete example, $A$ may be modeled as $M+E$, where $M$
is a deterministic matrix and $E$ a perturbation term sampled from a random
matrix distribution such as the real Ginibre ensemble (GinOE, [29]). By
definition, if $E\sim\mathrm{GinOE}$, $A_{ij}$ are $i.i.d.$ standard normal
random variables, therefore the probability measure associated with GinOE and
the Lebesgue measure are absolutely continuous with each other. Therefore, the
condition that almost every $A$ is identifiable is equivalent to requiring
$\lambda_{d\times d}\left(\Omega^{c}\right)=0$, where $\Omega^{c}:=M_{d\times
d}\setminus\Omega$ is the complement of $\Omega$, the collection of all
identifiable $A$.

To the best of our knowledge, the most systematic study of linear ODE
identifiability from a single observed trajectory is provided in [20]. In this
seminal work, Stanhope and colleagues derived several necessary and sufficient
conditions of identifiability that applies to a single trajectory.
Specifically, they proposed two concepts, called “identifiability for a single
trajectory” and “unconditional identifiability”, defined as follows.

######  Definition 1.4 (identifiability for a single trajectory).

System (3) is identifiable for a single trajectory in $\Omega\subset
M_{d\times d}$ and a given $\mathbf{x}_{0}$ iff for all $A,B\in\Omega$, $A\neq
B$, we have
$\mathbf{x}(t|A,\mathbf{x}_{0})\neq\mathbf{x}(t|B,\mathbf{x}_{0})$, for some
$0<t<T$.

######  Definition 1.5 (unconditional identifiability).

System (3) is unconditionally identifiable in $\Omega\in M_{d\times d}$ iff
for all $A,B\in\Omega$, $A\neq B$ implies that for each nonzero
$\mathbf{x}_{0}\in\mathbb{R}^{d}$,
$\mathbf{x}(t|A,\mathbf{x}_{0})\neq\mathbf{x}(t|B,\mathbf{x}_{0})$, for some
$0<t<T$.

Between these two definitions, Definition 1.5 adheres more to the traditional
definition of structural identifiability for nonlinear ODE systems. Roughly
speaking, it means that System (3) is identifiable from a single trajectory
initiated from _every_ $\mathbf{x}_{0}\in\mathbb{R}^{d}$. Unfortunately, it is
of little practical use for linear ODE system because no system satisfies this
condition for an unconstrained parameter estimation problem, namely,
$\Omega=M_{d\times d}$. In fact, we showed that (Supplementary Text, Section
S1): (a) when the dimension $d$ is odd, unconditional identifiability is not
attainable for all $\Omega\subseteq M_{d\times d}$, and (b) when $d$ is even,
unconditional identifiability is not attainable for all $\Omega\subseteq
M_{d\times d}$ such that $\lambda_{d\times d}\left(M_{d\times
d}\setminus\Omega\right)=0$. In summary, a large body of prior work in
identifiability analysis are geared towards nonlinear ODEs with arbitrarily
many observed trajectories, which is of little utility to linear ODEs, and
this issue cannot be fixed by simply removing a zero-measure set from their
definitions.

One major contribution of Stanhope and colleagues is that they established a
beautiful connection between the algebraic and geometric aspects of linear ODE
systems in [20, Theorem (3.4)]. We find it easier to state this important
result by first define the following minimalist definition of identifiability.

######  Definition 1.6 ($(A,\mathbf{x}_{0})$-identifiability).

For system (3), we call $A$ is identifiable at $\mathbf{x}_{0}$ if for all
$B\in M_{d\times d}$,
$\mathbf{x}(t|A,\mathbf{x}_{0})\neq\mathbf{x}(t|B,\mathbf{x}_{0})$, for some
$0<t<T$.

######  Remarks 1.

Definition 1.6 is not equivalent to Definition 1.4 applied to
$\Omega:=M_{d\times d}$ because in Definition 1.6, $A$ is fixed and $B$ is an
arbitrary matrix in $\Omega$, while in Definition 1.4, both $A$ and $B$ are
arbitrary matrices in $\Omega$. In short, Definition 1.6 is an intrinsic
property of a _single_ system, not a _collective_ property of a set of system
matrices.

Using this definition, [20, Theorem (3.4)] can be restated as follows: the
$(A,\mathbf{x}_{0})$-identifiability holds if and only if the solution curve
$\mathbf{x}(t|A)$ is not contained in a _proper invariant subspace_ of $A$.
Based on this powerful theoretical result, they proposed to use
$\kappa(X_{1})$, the condition number of the matrix of a subset of discrete
observations (see Section 3.3.1 for more details), to test the identifiability
for discrete data with noise in practice.

However, their study is not without shortcomings. First, they did not derive
the explicit structure of the largest subset $\Omega\subseteq M_{d\times d}$
for a give $\mathbf{x}_{0}$ in identifiability analysis for a single
trajectory, nor the equivalent class of all $B\in M_{d\times d}$ such that
$\mathbf{x}(t|B,t)=\mathbf{x}(t|A,t)$ when the system $A$ is deemed
unidentifiable at a given $\mathbf{x}_{0}$. Secondly, while using
$\kappa(X_{1})$ to check the practical identifiability of an ODE system is a
clever heuristic, it has much room for improvement because: a) not all data
are used in $\kappa(X_{1})$, therefore it does not utilize data efficiently;
b) measurement errors are not directly reflected in this score and there is no
analysis of the asymptotic properties of $\kappa(X_{1})$ from the statistical
perspective; and c) by definition, $\kappa(X_{1})$ depends on the availability
of data at multiple time points, so it requires solving the ODE numerically in
simulation studies, which can be time consuming for high-dimensional systems
and/or when a large set of systems are considered.

In this study, we first derive a closed-form representation of
$(A,\mathbf{x}_{0})$-unidentifiable class, which is defined in Definition 2.1
as the collection of system matrices that are not identifiable for a given
pair of $A$ and $\mathbf{x}_{0}$. We also provide explicit structures of the
equivalent class of unidentifiable systems due to repeated eigenvalues in $A$
in Supplementary Text, Section S5. We believe these results will be valuable
for future studies that combine a priori topological constraints (e.g.,
knowing which entries in $A$ are zero in advance) and identifiability. In
light this, we give a brief discussion of the best practice of using prior
information to resolve the identifiability issues in Supplementary Text,
Section S6. More systematic studies in this direction warrant a future study.

Secondly, we specify explicit, computable principles of
$(A,\mathbf{x}_{0})$-identifiability based on either $\mathbf{x}_{0}$ or the
entire solution trajectory. These results are presented in our Theorems 2.5
and 3.2. To assist practical identifiability analyses, we propose three
continuous scores: the initial condition-based identifiability score (ICIS,
denoted as $w_{0}^{*}$ in Equation (10)), the smoothed condition number (SCN,
denoted as $\tau$ in Equation (32)), and the practical identifiability score
(PIS, denoted as $w^{*}$ in Equation (37)), to solve the aforementioned
problems. ICIS only uses $A$ and $\mathbf{x}_{0}$, therefore it does not
require numerically solving the ODE before the identifiability analysis. We
think ICIS is most suitable for designing simulated ODE systems independent of
a specific set of real data. SCN and PIS use data from all time points, which
are more suitable for practical identifiability analysis with real data. Using
extensive simulation studies, we showed that SCN and PIS correlated with
practical identifiability significantly better than $\kappa(X_{1})$ when there
was noise in the data.

In addition, we studied the asymptotic properties of practical identifiability
for high-dimensional systems with randomly generated $A$ and $\mathbf{x}_{0}$.
We reached the following interesting conclusions: a) almost every system is
$(A,\mathbf{x}_{0})$-identifiable in the sense that $\mathrm{ICIS}>0$; and b)
when $d\to\infty$, almost all systems are practically unidentifiable in the
sense that $\mathrm{ICIS}\to 0$. These two seemly contradictory conclusions
suggest that the practical identifiability of high-dimensional ODE systems is
very different from that of low-dimensional systems, and classical
mathematical identifiability analyses are insufficient for analyzing _high-
dimensional_ real world applications. The focus must be shifted towards
practical identifiability analyses characterized by _continuous_ scores,
especially with the considerations from the stochastic perspective.

Last but not the least, we provide a user-friendly R package ode.ident, with
full documentation and examples, so practitioners with minimum programming
skills can analyze the identifiability of linear ODE systems. This R package
is available at https://github.com/qiuxing/ode.ident.

##  2 $(A,\mathbf{x}_{0})$-identifiability

In this section, we focus on the mathematical inverse problem for one fully
observed trajectory. Namely, we assume that we have the complete observation
of _one_ solution curve
$\mathbf{x}(t)=(x_{1}(t),\dots,x_{d}(t))^{\prime}\in\mathbb{R}^{d}$) governed
by Equation (3) and its derivative on $[0,T]$, with no measurement error.

First, let us define the $(A,\mathbf{x}_{0})$-unidentifiable class as follows.

######  Definition 2.1 ($(A,\mathbf{x}_{0})$-unidentifiable class).

For a given system matrix $A\in M_{d\times d}$ and initial condition
$\mathbf{x}_{0}\in\mathbb{R}^{d}$, the $(A,\mathbf{x}_{0})$-unidentifiable
class, denoted by $[A]_{\mathbf{x}_{0}}$, is a subset of matrices in
$M_{d\times d}$ such that

$B\in[A]_{\mathbf{x}_{0}}\quad\text{iff}\quad\mathbf{x}(t|A,\mathbf{x}_{0})=\mathbf{x}(t|B,\mathbf{x}_{0}).$
(5)

In other words, two system matrices $A,B$ are in the same unidentifiable class
if and only if they produce the same solution trajectory at $\mathbf{x}_{0}$.

The overarching goal of this section is to understand the structure of
$[A]_{\mathbf{x}_{0}}$, and the conditions under which this class contains
only one member, therefore $A$ can be uniquely determined by the trajectory
$\mathbf{x}(t|A,\mathbf{x}_{0})$. To this end, we need to introduce an
important geometric concept called invariant subspace, which is a
generalization of eigenvectors, and its connection to the Jordan decomposition
of $A$ in Section 2.1333These concepts and results can be found in many
graduate level matrix analysis textbooks, e.g., [30]..

###  2.1 Jordan Decomposition and Invariant subspaces

######  Definition 2.2 (Invariant subspace).

An invariant subspace of a square matrix $A_{d\times d}$ is a linear subspace
$L\subseteq R^{d}$ such that for all $\mathbf{x}\in L$, $A\mathbf{x}\in L$. We
say $L$ is a _proper_ invariant subspace if $L\neq\mathbb{R}^{d}$.

By definition, we see that if a vector $\mathbf{x}$ is in a proper invariant
subspace $L$ of $A$, $A\mathbf{x}$ must also stay in $L$. Using mathematical
induction, we see that $A^{n}\mathbf{x}\in L$ for every positive integer $n$.
With a little more work, it can be proven that $e^{tA}\mathbf{x}\in L$ for
$t\in[0,T]$, where $e^{tA}$ is the matrix exponential of $tA$.

The following proposition states that the intersection and linear span (the
combination) of two invariant subspaces are invariant subspaces.

######  Proposition 2.1.

If $L_{1}$ and $L_{2}$ are invariant subspaces of $A$, then

  1. 1.

$L_{1}\cap L_{2}$ is an invariant subspace of $A$;

  2. 2.

$\mathrm{span}(L_{1},L_{2})$ is an invariant subspace of $A$.

In other words, the collection of invariant subspaces of $A$ forms a
_lattice_.

Based on random matrix theory [29, 31, 32], we know that almost every (w.r.t.
the Lebesgue measure on $M_{d\times d}$) $A\in M_{d\times d}$ has $d$ distinct
eigenvalues. This conclusion also holds for probability measures associated
with most random matrix ensembles such as Ginibre ensemble, Gaussian
orthogonal ensemble, Wishart ensemble, etc.[32]

Consequently, almost every $A\in M_{d\times d}$ has the following Jordan
decomposition

$\begin{gathered}A=Q\Lambda Q^{-1},\qquad\Lambda=\begin{pmatrix}J_{1}&&\\\
&\ddots&\\\ &&J_{K}\end{pmatrix},\qquad
Q=\begin{pmatrix}[c|c|c|c]Q_{1}&Q_{2}&\dots&Q_{K}\end{pmatrix}.\\\
J_{k}=\begin{cases}c_{k},&k=1,2,\dots,K_{1},\\\
\left(\begin{smallmatrix}a_{k}&-b_{k}\\\
b_{k}&a_{k}\end{smallmatrix}\right),&k=K_{1}+1,K_{1}+2,\dots,K.\end{cases}\\\
\dim Q_{k}=\begin{cases}1,&k=1,\dots,K_{1},\\\
2,&k=K_{1}+1,\dots,K.\end{cases}\end{gathered}$ (6)

In other words, $A$ can be decomposed into $K=K_{1}+K_{2}$ Jordan blocks, the
first $K_{1}$ such blocks are $1\times 1$ blocks corresponding with real
eigenvalues (those $c_{k}$ in Equation (6)); and the rest $K_{2}$ blocks are
$2\times 2$ blocks corresponding with complex eigenvalues $a_{k}\pm b_{k}i$.
There is a corresponding column-wise decomposition of matrix $Q$, such that
each $Q_{k}$ contains: (a) a single column vector of $Q$ which is the
eigenvector of $c_{k}$, or (b) two column vectors in $Q$ such that
$Q_{k}:=(\mathbf{v}_{k1}|\mathbf{v}_{k2})$, which are the “eigenvectors”
associated with $a_{k}\pm b_{k}i$.

Note that the word “eigenvector” in case (b) refers to a generalization of
true eigenvectors. In fact, those $2\times 2$ Jordan blocks do not have real
eigenvectors; instead, each of them is associated with a 2-dimensional
invariant subspace of $A$ and $Q_{k}=(\mathbf{v}_{k1}|\mathbf{v}_{k2})$ is a
basis of this 2-dimensional invariant subspace.

We would like to point out that based on simple enumeration of dimensions, we
have $d=K_{1}+2K_{2}$, and

$Q_{k}=\begin{cases}Q_{\cdot k},&k=1,2,\dots,K_{1},\\\
(\mathbf{v}_{k1}|\mathbf{v}_{k2}),\quad\mathbf{v}_{k1}=Q_{\cdot
2k-K_{1}-1},\;\mathbf{v}_{k2}=Q_{\cdot
2k-K_{1}},&k=K_{1}+1,\dots,K.\end{cases}$ (7)

In other words, $\mathbf{v}_{k1}$ and $\mathbf{v}_{k2}$ in $Q_{k}$ are the
($2k-K_{1}-1$)-th and ($2k-K_{1}$)-th column vectors of $Q$, respectively. For
convenience, we define the following correspondences between $i$ (the original
dimension in $J$) and $k$ (the index of invariant subspaces):

$i(k):=2k-K_{1}-1,\qquad k(i):=\begin{cases}i,&i=1,2,\dots,K_{1},\\\
K_{1}+\lceil\frac{i-K_{1}}{2}\rceil,&i=K_{1}+1,\dots,K.\end{cases}$ (8)

Using the above notation, $\mathbf{v}_{k1}=Q_{\cdot
i(k)},\;\mathbf{v}_{k2}=Q_{\cdot i(k)+1}$.

######  Theorem 2.2.

Let $L_{k}:=\mathrm{span}(Q_{k})$. Each $L_{k}$ is an invariant subspace of
$A$. Furthermore, if $L$ is an invariant subspace of $A$, it can always be
decomposed as

$L=\mathrm{span}\bigcup_{i\in S}L_{k},\qquad
S\subseteq\left\\{1,2,\dots,K\right\\}.$ (9)

From now on, we assume that $A$ has $d$ distinct eigenvalues and can be
decomposed as $A=Q\Lambda Q^{-1}$ in Equation (6). The case in which $A$ has
repeated eigenvalues will be discussed in Supplementary Text, Section S5. For
convenience, we will also denote $L_{0}:=\left\\{0_{d}\right\\}$, the trivial
proper invariant subspace of $A$ that contains only the origin.

###  2.2 Initial Condition-based Identifiability Score (ICIS)

One of our main conclusion is that the $(A,\mathbf{x}_{0})$-identifiability
defined in Definition 1.6 can be determined by the initial condition-based
identifiability score (ICIS) defined as follows.

######  Definition 2.3.

Let $\mathbf{x}_{0}\in\mathbb{R}^{d}$ and
$\tilde{\mathbf{x}}_{0}:=Q^{-1}\mathbf{x}_{0}\in\mathbb{R}^{d}$. We define the
$w_{0}^{*}$ statistic in the following equation as the Initial Condition-based
Identifiability Score (ICIS):

$\begin{gathered}w_{0,k}:=\begin{cases}\tilde{\mathbf{x}}_{0,k}\in\mathbb{R}^{1},&k=1,2,\dots,K_{1},\\\
(\tilde{\mathbf{x}}_{0,i(k)},\tilde{\mathbf{x}}_{0,i(k)+1})^{\prime}\in\mathbb{R}^{2},&k=K_{1}+1,\dots,K.\end{cases}.\\\
w_{0}^{*}:=\min_{k}|w_{0,k}|.\end{gathered}$ (10)

Here $|w_{0,k}|$ is the absolute value of $w_{0,k}$ for $k=1,2,\dots,K_{1}$,
and the Euclidean norm of $w_{0,k}$ for $k=K_{1}+1,K_{1}+2,\dots,K$.

From the geometric perspective, $\mathbf{x}_{0}$ can be decomposed into a
linear combination (oblique projections) of $Q_{k}$, and $w_{0,k}$ are the
linear coefficients of such a decomposition

$\mathbf{x}_{0}=Q\tilde{\mathbf{x}}_{0}=\sum_{k=1}^{K_{1}}w_{0,k}Q_{k}.$ (11)

Heuristically speaking, if $w_{0,k}=0$ (in $\mathbb{R}^{1}$ or
$\mathbb{R}^{2}$), $\mathbf{x}_{0}$ does not contain any information from
$L_{k}$. This is because in this case, $\mathbf{x}_{0}\in L_{-k}$, where
$L_{-k}$ is the invariant subspace of $A$ that _excludes_ $L_{k}$, which
implies that the _entire trajectory_ , $\mathbf{x}(t|A,\mathbf{x}_{0})$, is in
$L_{-k}$ (see the discussion in the beginning of Section 2.1). These ideas are
summarized in Lemma 2.3 below.

######  Lemma 2.3.

The following two statements are equivalent

  1. 1.

There exists a proper invariant subspace $L\subsetneq\mathbb{R}^{d}$ of $A$,
such that $\mathbf{x}_{0}\in L$.

  2. 2.

There exists $k\in\left\\{1,2,\dots,K\right\\}$, such that $|w_{0,k}|=0$, or
equivalently, $w_{0}^{*}=0$.

Now we are ready to present the following theorem:

######  Theorem 2.4 (Computational criterion for
$(A,\mathbf{x}_{0})$-identifiability).

Assuming that an ODE system $A$ has $d$ distinct eigenvalues. This system is
identifiable at $\mathbf{x}_{0}$ if and only if the ICIS is nonzero.

###### Proof.

Based on Lemma 3.2 and Theorem (3.4) in [20], we know that system $A$ is
identifiable at $\mathbf{x}_{0}$ if and only if $\mathbf{x}_{0}$ is not
contained in a proper invariant subspace of $A$, which is equivalent to
$w_{0}^{*}\neq 0$ based on our Lemma 2.3. ∎

Theorem 2.4 implies that, when $A$ is not identifiable at $\mathbf{x}_{0}$,
there must be a nonempty subset $S_{0}\subset\left\\{1,2,\dots,K\right\\}$
such that $|w_{0,k}|=0$, for $k\in S_{0}$. WLOG, we assume that
$\mathbf{x}_{0}\neq 0_{d}$, so its complement set
$S_{+}:=\left\\{1,2,\dots,K\right\\}\setminus S_{0}$ must be nonempty. By
construction, $L_{+}:=\mathrm{span}\bigcup_{k\in S_{+}}L_{k}$ is a proper
invariant subspace, and $\mathbf{x}_{0}\in L_{+}$.

The two index sets $S_{0}$ and $S_{+}$ induce the following diagonal binary
matrices

$I_{0}:=\mathrm{diag}(a_{i}),\quad a_{i}=\begin{cases}1,&k(i)\in S_{0},\\\
0,&k(i)\neq S_{0}.\end{cases}\qquad I_{+}:=I_{d\times d}-I_{0}.$ (12)

They can be used to construct the following decomposition of the eigenvector
matrix $Q$ and the Jordan block matrix $J$

$\begin{gathered}Q_{0}:=QI_{0},\qquad Q_{+}:=QI_{+},\qquad Q=Q_{0}+Q_{+}.\\\
J_{0}:=I_{0}JI_{0}=JI_{0},\qquad J_{+}:=I_{+}JI_{+}=JI_{+},\qquad
J=J_{0}+J_{+}.\end{gathered}$ (13)

Intuitively, $Q_{0}$ and $J_{0}$ replace column vectors in $Q$ and blocks in
$J$ into zeros if they belong to $L_{+}$. $Q_{+}$ and $J_{+}$ are defined in
exactly the opposite way.

Using these notation, we describe the explicit structure of
$(A,\mathbf{x}_{0})$-unidentifiable class in the following Theorem.

######  Theorem 2.5 (Structure of the unidentifiable classes).

$(A,\mathbf{x}_{0})$-unidentifiable class has the following explicit structure

$\begin{split}[A]_{\mathbf{x}_{0}}&=Q\left(J_{+}+I_{0}(J_{0}+D)I_{0}\right)Q^{-1}=A+Q(I_{0}DI_{0})Q^{-1},\qquad
D\in M_{d\times d}.\end{split}$ (14)

In other words, two matrices $A_{1}$ and $A_{2}$ are in the same
$(A,\mathbf{x}_{0})$-unidentifiable class (written as $A_{1}\sim A_{2}$) iff
there exist $D\in M_{d\times d}$, such that
$A_{1}-A_{2}=Q(I_{0}DI_{0})Q^{-1}$.

###### Proof.

See Section S4, Supplementary Text. ∎

######  Remarks 2.

The degrees of freedom in $[A]_{\mathbf{x}_{0}}$ is controlled by
$D_{0}:=I_{0}DI_{0}$, which has $d_{0}^{2}$ degrees of freedom (not $d^{2}$).
This is because by construction, $D_{0}$ is a sparse matrix such that its
$ij$th element satisifies

$D_{0,ij}=0,\qquad\text{if }k(i),k(j)\in S_{0}.$ (15)

###  2.3 A 3-Dimensional Example

######  Example 1.

In this example, the system matrix $A$ and its Jordan canonical form are given
as follows:

$\begin{gathered}A=\begin{pmatrix}0&1&-1\\\ 2&0&0\\\
3&1&0\end{pmatrix}=QJQ^{-1},\qquad J=\begin{pmatrix}-1&0&0\\\
0&1/2&\sqrt{7}/2\\\ 0&-\sqrt{7}/2&1/2\end{pmatrix}.\\\
Q\approx\begin{pmatrix}0.408&0&0.316\\\ -0.816&0.418&0.158\\\
-0.408&0.837&0\end{pmatrix},\qquad
Q^{-1}\approx\begin{pmatrix}0.612&-1.225&0.612\\\ 0.299&-0.598&1.494\\\
2.372&1.581&-0.791\end{pmatrix}.\end{gathered}$ (16)

Based on earlier discussions, $A$ has two proper invariant subspaces.
$L_{1}:=\mathrm{span}(Q_{\cdot 1})$ is a one-dimensional space corresponding
with the real eigenvalue $\lambda_{1}=-1$, and
$L_{2}:=\mathrm{span}\left(Q_{\cdot 2},Q_{\cdot 3}\right)$ is a two-
dimensional space corresponding with
$\lambda_{2},\lambda_{3}=1/2\pm\sqrt{7}i/2$. Here $Q_{\cdot j}$ is the $j$th
column vector of matrix $Q$.

Let us consider the following two initial conditions:

$\mathbf{x}_{0}^{(a)}:=Q\begin{pmatrix}2\\\ -1\\\
0\end{pmatrix}\approx\begin{pmatrix}0.816\\\ -2.051\\\
-1.653\end{pmatrix},\qquad\mathbf{x}_{0}^{(b)}:=Q\begin{pmatrix}0\\\ -2\\\
3\end{pmatrix}\approx\begin{pmatrix}0.949\\\ -0.362\\\ -1.673\end{pmatrix}.$

Notice that $\mathbf{x}_{0}^{(a)}$ contains information from both $L_{1}$ and
$L_{2}$, but $\mathbf{x}_{0}^{(b)}$ only contains information from $L_{2}$.
Based on earlier discussions, $A$ is identifiable at $\mathbf{x}_{0}^{(a)}$
but not $\mathbf{x}_{0}^{(b)}$. Using Equation (14), the
$(A,\mathbf{x}_{0})$-unidentifiable class in the latter case can be
represented as follows.

$\begin{gathered}I_{+}=\begin{pmatrix}0&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix},\quad I_{0}=\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&0\end{pmatrix},\\\ I_{0}DI_{0}=\begin{pmatrix}b&0&0\\\ 0&0&0\\\
0&0&0\end{pmatrix},\quad b\in\mathbb{R}.\\\ J_{+}=\begin{pmatrix}0&0&0\\\
0&1/2&\sqrt{7}/2\\\ 0&-\sqrt{7}/2&1/2\end{pmatrix},\qquad
Q(I_{0}DI_{0})Q^{-1}=bQ_{0}Q^{-1},\\\
\begin{split}[A]_{\mathbf{x}_{0}^{(b)}}&=QJ_{+}Q^{-1}+bQ_{0}Q^{-1}\\\
&=\dfrac{1}{4}\begin{pmatrix}1&2&-3\\\ 6&4&-2\\\
11&6&-1\end{pmatrix}+\dfrac{b}{4}\begin{pmatrix}1&-2&1\\\ -2&4&-2\\\
-1&2&-1\end{pmatrix}.\end{split}\end{gathered}$ (17)

Here $b\in\mathbb{R}$ is an arbitrary parameter, and the
$(A,\mathbf{x}_{0})$-unidentifiable class is characterized by $b$ times a
_full matrix_ , therefore different choices of $b$ affects all nine elements
in $A$. Consequently, we cannot determine the value of _any_ entry in $A$
without additional information. In fact, we cannot even determine whether a
particular $A_{ij}$ is zero or nonzero, which is sometimes referred to as the
network topology of $A$. For example, if we set $b=-1$, we get the original
system $A$ specified in Equation (16) with four zero entries. When we set
$b=3$, we obtain an equivalent matrix with completely different topology and
values than $A$

$\tilde{A}=\begin{pmatrix}1&-1&0\\\ 0&4&-2\\\ 2&3&-1\end{pmatrix}.$ (18)

It is easy to check that:

$e^{t\tilde{A}}\mathbf{x}_{0}^{(a)}\neq e^{tA}\mathbf{x}_{0}^{(a)},\qquad
e^{t\tilde{A}}\mathbf{x}_{0}^{(b)}=e^{tA}\mathbf{x}_{0}^{(b)}.$

Before we move on to the next topic (practical identifiability), we would like
to present two auxiliary results that are useful in practice.

  1. 1.

In Supplementary Text, Section S5, we provide a detailed analysis of
identifiability issues induced by repeated eigenvalues in $A$, and provided
the closed-form structure of unidentifiable class for these matrices in
Equation (S.23). Based on these results, we recommend researchers avoid
systems that have nearly identical eigenvalues in designing of simulation
studies.

  2. 2.

in Supplementary Text, Section S6, we show that while it is possible to use
prior information in the form of structural constraints to resolve the
identifiability issue of a problematic system $(A,\mathbf{x}_{0})$, such
constraints must be compatible with the said system, otherwise: (a) the system
may still suffer from the identifiability issue; or (b) under these
constraints, no system can generate the observed solution curve.

##  3 Data-based Identifiability Scores

We now focus on the following practical problem: to quantify the
$(A,\mathbf{x}_{0})$-identifiability from imperfect observations in real world
applications. To this end, we assume that the observed data is a set of
discrete and noisy observations on a time grid
$\left\\{t_{1},\dots,t_{J}\right\\}$:

$\begin{gathered}y_{ij}:=x_{i}(t_{j})+\epsilon_{ij},\qquad\epsilon_{ij}\sim
F_{\epsilon},\qquad i=1,2,\dots,d,\quad j=1,2,\dots,n.\end{gathered}$ (19)

In the above equation, $F_{\epsilon}$, the probability distribution of
measurement error, is assumed to be absolutely continuous w.r.t. the Lebesgue
measure on $M_{d\times J}$. For convenience, we will also use collective
notations $X=\left\\{x_{i}(t_{j})\right\\}\in M_{d\times n}$,
$Y=\\{y_{ij}\\}\in M_{d\times n}$, and $\bm{\epsilon}=\\{\epsilon_{ij}\\}\in
M_{d\times n}$. With these matrix notations, Equation (19) can be simplified
as $Y=X+\bm{\epsilon}$.

###  3.1 Minimal Signals for Reconstructing $A$

In this section, we demonstrate that even for a theoretically identifiable
system, if the “signal” in a subspace is too small, we are still not able to
reconstruct $A$ in practice.

######  Example 2.

We consider a two-dimensional system

$\begin{gathered}A=Q\Lambda Q^{-1}=\begin{pmatrix}-1.875&2.382\\\
2.382&-4.625\end{pmatrix},\\\
Q=\mathrm{Rot}_{\pi/6}=\begin{pmatrix}\frac{\sqrt{3}}{2}&-\frac{1}{2}\\\
\frac{1}{2}&\frac{\sqrt{3}}{2}\end{pmatrix},\quad
Q^{-1}=\mathrm{Rot}_{-\pi/6}=\begin{pmatrix}\frac{\sqrt{3}}{2}&\frac{1}{2}\\\
-\frac{1}{2}&\frac{\sqrt{3}}{2}\end{pmatrix},\quad\Lambda=\begin{pmatrix}-1/2&0\\\
0&-6\end{pmatrix}.\\\ \end{gathered}$ (20)

We generate the solution curve $\mathbf{x}(t)$ from this system and record its
values $x_{i}(t_{j})$ at $n=101$ equally spaced time points on $[0,1]$,
$t_{1}=0,t_{2}=0.01,\dots,t_{101}=1$. A small normal measurement error,
$\epsilon_{ij}\sim N(0,0.01^{2})$, is added to each observation.

It is easy to see that $A$ has two one-dimensional proper invariant subspaces,
$L_{1}=\mathrm{span}((\frac{\sqrt{3}}{2},\frac{1}{2})^{\prime})$,
$L_{2}=\mathrm{span}((-\frac{1}{2},\frac{\sqrt{3}}{2})^{\prime})$. Let us
consider two initial conditions

$\mathbf{x}_{0}^{(A)}=\begin{pmatrix}1\\\
1\end{pmatrix},\qquad\mathbf{x}_{0}^{(B)}=\begin{pmatrix}1.72\\\
1\end{pmatrix}.$

It can be shown that $A$ is identifiable at both $\mathbf{x}_{0}^{(A)}$ and
$\mathbf{x}_{0}^{(B)}$. Of note, we would like to mention that this analysis
can be done by applying the ICISAnalysis() function in our R package.

Using the functional two-stage method (see Section S7.2, Supplementary Text),
we are able to estimate $A$ and produce two fitted curves for both cases. The
fitted curves, denoted by $(\hat{x}_{1}(t),\hat{x}_{2}(t))^{\prime}$, look
reasonable in Figure 1 for both cases.

Figure 1: An illustration of the fitted solution trajectories in Example 2.
The black curve is $\hat{x}_{1}(t)$ and the red curve is $\hat{x}_{2}(t)$ in
both sub-figures. The discrete data (illustrated by dots in this figure) are
observed on a total number of $n=101$ time points evenly assigned on $[0,1]$.
A small normal measurement error, $\epsilon_{ij}\sim N(0,0.01^{2})$, is added
to each observation. (A): $\mathbf{x}_{0}^{(A)}=(1,1)^{\prime}$. (B):
$\mathbf{x}_{0}^{(B)}=(1.72,1)^{\prime}$. Case (A) is practically identifiable
with the functional two-stage method while Case (B) is not.

However, the two reconstructed system matrices are quite different:

$\begin{gathered}\hat{A}^{(A)}=\begin{pmatrix}-1.58&1.93\\\
1.93&-3.91\end{pmatrix},\qquad\|\hat{A}^{(A)}-A\|_{F}^{2}=1.01,\\\
\hat{A}^{(B)}=\begin{pmatrix}-0.69&0.33\\\
-1.81&2.62\end{pmatrix},\qquad\|\hat{A}^{(B)}-A\|_{F}^{2}=75.65.\end{gathered}$

Here $\|\cdot\|_{F}$ is the Frobenius norm of a matrix. It is clear that $A$
is practically identifiable at $\mathbf{x}_{0}^{(A)}$ only, not at
$\mathbf{x}_{0}^{(B)}$. From this example, we see that even when the ODE
system is mathematical identifiable, its practical identifiability may still
be an issue.

In-depth analysis shows that the unidentifiability issue in case (B) is due to
the fact that $\mathbf{x}_{0}^{(B)}$ is “almost” contained in $L_{1}$, so that
$L_{2}$ had only a tiny bit of information. Denote the basis in $L_{1}$ and
$L_{2}$ as $Q_{1}=(\frac{\sqrt{3}}{2},\frac{1}{2})^{\prime}$ and
$Q_{2}=(-\frac{1}{2},\frac{\sqrt{3}}{2})^{\prime}$, respectively. We have

$\begin{gathered}w_{0,1}^{(A)}=\langle\mathbf{x}_{0}^{(A)},\;Q_{1}\rangle=1.366,\qquad\mathrm{ICIS}^{(A)}=w_{0,2}^{(A)}=\langle\mathbf{x}_{0}^{(A)},\;Q_{2}\rangle=0.366.\\\
w_{01}^{(B)}=\langle\mathbf{x}_{0}^{(B)},\;Q_{1}\rangle=1.990,\qquad\mathrm{ICIS}^{(B)}=w_{02}^{(B)}=\langle\mathbf{x}_{0}^{(B)},\;Q_{2}\rangle=\mathbf{0.006}.\end{gathered}$

Based on the above analysis, it is easy to see that the small value of
$\mathrm{ICIS}^{(B)}=0.006$ causes the numerical problem in estimating $A$ in
case (B).

###  3.2 Identifiability for High-dimensional Systems

Based on Theorem 2.4, $A$ is identifiable at $\mathbf{x}_{0}$ if and only if
$\mathbf{x}_{0}$ is not located in a proper invariant subspace of $A$. Because
there are only finitely many ($2^{K}-1$ of them, to be more precise) proper
invariant subspaces of $A$, and each of them has dimension strictly less than
$d$ (the “proper” part of the definition), the union of all proper invariant
subspaces is only a zero-measure set of $\mathbb{R}^{d}$. In this regard, as
long as $A$ does not have repeated eigenvalues (which is true for almost every
$A\in M_{d\times d}$), $A$ is mathematically identifiable at almost every
$\mathbf{x}_{0}\in\mathbb{R}^{d}$. This fact is probably the main reason why
not many mathematicians have paid much attention to the identifiability
problem of linear ODE systems.

However, as is shown in Example 2, to have a reliable estimate of $A$ requires
more than just a _qualitative_ statement that $\mathbf{x}_{0}$ does not lie in
any proper invariant subspace of $A$. We need to ensure that when we decompose
$\mathbf{x}_{0}$ into a linear combination of components from $L_{k}$, each
one of them has _enough_ information, so that we can reconstruct the
corresponding sub-system on $L_{k}$ with noisy observations. This is the main
motivation for us to propose ICIS, a _quantitative_ measure of
identifiability.

Knowing that the collection of all proper invariant subspace has measure zero
in $\mathbb{R}^{d}$, the readers may think that while practical
identifiability issues do exist, they must be rare in practice. Unfortunately,
these issues are not that unusual when $d$ is large, in which case those
practically identifiable systems are the exceptions instead. In Supplementary
Text, Section S2, we proved that a large class of symmetric random ODE systems
are practically unidentifiable when $d\to\infty$, as stated in the following
theorem

######  Theorem 3.1.

Let us assume that:

  1. (a)

The system matrix $A\in M_{d\times d}$ is sampled from a symmetric, real-
valued random matrix ensemble with probability measure $p(A)$ that is
statistically invariant to orthogonal transformations, namely,

$p(A)=p(TAT^{\prime}),\qquad\forall T\in O(d).$ (21)

  2. (b)

The initial condition $\mathbf{x}_{0}\in\mathbb{R}^{d}$ is sampled from a
random distribution that is independent of $A$ and satisfies

$\lim_{d\to\infty}\frac{E|\mathbf{x}_{0}|^{2}}{d^{3}}=0.$ (22)

Based on the above two assumptions, the ICIS converges to zero in $L^{2}$,
namely,

$E\left(w_{0}^{*}(A,\mathbf{x}_{0})\right)^{2}\longrightarrow
0,\quad\text{when $d\to\infty$.}$ (23)

###### Proof.

The proof is provided in Section S2, Supplementary Text. ∎

######  Remarks 3.

Perhaps the most well known random matrix ensemble that satisfies Assumption
(a) is the Gaussian Orthogonal Ensemble (GOE, [32]). Many other ensembles also
satisfy this condition, such as the Wishart ensemble, Jacobi orthogonal
ensemble, etc. In fact, according to Weyl’s lemma [33, 34], a random matrix
ensemble is orthogonally invariant as long as its distribution function has
the following trace representation

$p(H)=\phi\left(\mathrm{tr}H,\mathrm{tr}H^{2},\dots,\mathrm{tr}H^{n}\right).$
(24)

Assumption (b) is a very weak condition that should be satisfied in almost all
practical applications. If $\mathbf{x}_{0,i}$ has finite second order moments,
and

$\sup_{i}Ex_{0,i}^{2}=\mu_{2}<\infty,$

we have

$E|\mathbf{x}_{0}|^{2}:=\sum_{i=1}^{d}Ex_{0,i}^{2}\leqslant
d\cdot\mu_{2}=O(d^{1}),\qquad\frac{E|\mathbf{x}_{0}|^{2}}{d^{3}}=O(d^{-2})\to
0.$

In this case, $x_{0,i}$ do not have to be independent nor identically
distributed.

We use the following simple and concrete example to illustrate the issue of
practical identifiability described by Theorem 3.1.

######  Example 3.

Let $d=100$ and assume that $A$ is an arbitrary diagonal matrix in $M_{d\times
d}$. By construction, all eigenvalues are real and $Q=I_{d}$, therefore
$w_{0,i}=x_{0,i}$, $i=1,2,\dots,d$. Let $\mathbf{x}_{0}\sim N(0_{d},I_{d})$,
in other words, $x_{0,i}$ are generated from $i.i.d.$ $N(0,1)$. For
simplicity, we write $r_{i}=|w_{0,i}|$. Apparently, $r_{i}$ are standard half
normals with relatively large expectations $Er_{i}=\sqrt{\frac{2}{\pi}}\approx
0.8$. This fact seems to suggest that, as long as the measurement error is
small ($\sigma\ll 0.8$), we would have enough information to reconstruct $A$.

In reality, the _smallest_ member of $r_{i}$ (denoted by $r_{(1)}$), has a
distribution that is statistically much smaller than a standard half normal.
Using numerical integration, we found that $Er_{(1)}\approx 0.012$, which is
66 times smaller than $Er_{i}$. Based on the lessons we learned from Example
2, we anticipate that it is almost impossible to estimate $A$ accurately in
this case.

Finally, we conduct a mini-simulation to illustrate that ODEs with random
_asymmetric_ system matrices also suffer from the identifiability issues
stated in Theorem 3.1. Specifically, we randomly generate 50 $A$ from the
standard GinOE with $d=3,5,100$ dimensions, and pair them with 50
$\mathbf{x}_{0}$ sampled from $N(0_{d},I_{d})$. We compute the ICIS for those
$(A,\mathbf{x}_{0})$ and plot them in Figure 2. We see that larger dimension
is associated with smaller ICIS, which implies that these ODE systems are more
difficult to be numerically identified.

Figure 2: Larger dimension is associated with smaller ICIS ($w_{0}^{*}$),
which implies that high-dimensional ODE systems are more difficult to be
reconstructed numerically.

In light of the above discussions, to design a well-behaved, identifiable
high-dimensional ODE system in simulations, one needs to ensure that:

  1. 1.

$A$ has nice mathematical properties, such as distinct eigenvalues; stability;
no high-frequency components, etc.

  2. 2.

$\mathbf{x}_{0}$ should not be “randomly” generated; instead, it should be
generated in a way such that $\mathrm{ICIS}(A,\mathbf{x}_{0})$ is not too
small.

###  3.3 Practical Identifiability Score (PIS)

Recall that ICIS does not depend on the full trajectory $\mathbf{x}(t)$,
therefore it is most useful in designing simulation studies. In real world
applications, it is preferable to define identifiability scores that use the
entirety of $Y$ (the discrete data measured at all time points) to quantify
the practical identifiability of the system. Let $\hat{A}(Y)$ be an estimator
of $A$ given the discrete observation. One way to quantify the practical
identifiability is to use the numerical sensitivity of $\hat{A}$, which can be
defined as the mean squared error, $\mathrm{MSE}:=E\|\hat{A}-A\|_{F}^{2}$.
Unfortunately, there are many different ways to estimate $A$, thus it is
impossible to develop a _universal_ quantity that works for all estimators. In
this section, we propose two scores based on a class of the two-stage methods,
under the assumption that $\sigma^{2}$ is small enough so that $\hat{A}-A$ is
small. These proposed scores are compared to $\kappa(Y_{1})$, the practical
identifiability measure proposed by Stanhope and colleagues in [20] in our
simulation studies.

####  3.3.1 Stanhope’s Condition Number

Stanhope and colleagues proposed to use $\kappa(Y_{1})$ as a practical
identifiability measure in [20]. Here $Y_{1}:=[Y_{\cdot 1}|Y_{\cdot
2}|\dots|Y_{\cdot d}]\in M_{d\times d}$ is the matrix of discrete data
evaluated at the first $d$ time points;
$\kappa(Y_{1}):=\|Y_{1}\|_{F}\|Y_{1}^{-1}\|_{F}$ is the _condition number_ of
$Y_{1}$, which is a quantitative measure of numerical stability of
$Y_{1}^{-1}$. It was motivated by the fact that if $Y_{1}$ is invertible and
there is no noise in the discrete observations,

$e^{\Delta tA}=Y_{2}Y_{1}^{-1},\qquad Y_{2}:=[Y_{\cdot 2}|Y_{\cdot
3}|\dots|Y_{\cdot d+1}]\in M_{d\times d}.$ (25)

####  3.3.2 Functional Two-stage Methods

######  Definition 3.1 (pairwise $L^{2}$-inner product matrix).

Let $\mathbf{x}(t)\in\mathbb{R}^{d}$ and
$\mathbf{y}(t)\in\mathbb{R}^{d^{\prime}}$ be two multidimensional functions
defined on $[0,T]$. We use notation
$\langle\mathbf{x}(t),\;\mathbf{y}(t)\rangle$ to refer to the pairwise
$L^{2}$-inner product matrix between them, which is a $d\times d^{\prime}$
matrix in which

$\langle\mathbf{x}(t),\;\mathbf{y}(t)\rangle_{i,j}:=\langle
x_{i}(t),\;y_{j}(t)\rangle_{L^{2}}=\int_{0}^{T}x_{i}(t)y_{j}(t)\mathrm{d}t.$

For convenience, we denote $\langle\mathbf{x}(t),\;\mathbf{x}(t)\rangle$ by
$\Sigma_{\mathbf{x}\mathbf{x}}$ when there is no confusion. Obviously,
$\Sigma_{\mathbf{x}\mathbf{x}}$ is a symmetric positive semi-definite matrix.
It is singular if and only if zero is one of its eigenvalues, in which case we
know

$\mathbf{v}_{0}^{\prime}\Sigma_{\mathbf{x}\mathbf{x}}\mathbf{v}_{0}=0,\qquad\mathbf{v}_{0}\in\mathbb{R}^{d},\quad\mathbf{v}_{0}\neq
0_{d}.$ (26)

where $\mathbf{v}_{0}$ is an eigenvector associated with the zero eigenvalue.

######  Theorem 3.2.

Let $\mathbf{x}(t)$ be an observed solution trajectory governed by ODE system
$A$ and initiated at $\mathbf{x}_{0}:=\mathbf{x}(0)$. This ODE system is
identifiable at $\mathbf{x}_{0}$ if and only if the pairwise inner product
matrix $\Sigma_{\mathbf{x}\mathbf{x}}$ is nonsingular (invertible).

###### Proof.

See Section S4, Supplementary Text. ∎

Notice that the ODE $D\mathbf{x}(t)=A\mathbf{x}(t)$ implies

$\langle
D\mathbf{x}(t),\;\mathbf{x}(t)\rangle=A\langle\mathbf{x}(t),\;\mathbf{x}(t)\rangle\Longrightarrow
A=\Sigma_{D\mathbf{x},\mathbf{x}}\cdot\Sigma_{\mathbf{x}\mathbf{x}}^{-1}.$
(27)

Here $\Sigma_{D\mathbf{x},\mathbf{x}}:=\langle
D\mathbf{x}(t),\;\mathbf{x}(t)\rangle$ is the pairwise inner product matrix
between $D\mathbf{x}(t)$ and $\mathbf{x}(t)$. This fact motivated the
following two-stage methods, which is a class of estimators of $A$:

$\hat{A}:=\hat{\Sigma}_{D\mathbf{x},\mathbf{x}}\cdot\hat{\Sigma}_{\mathbf{x}\mathbf{x}}^{-1}.$
(28)

Here $\hat{\Sigma}_{D\mathbf{x},\mathbf{x}}$ and
$\hat{\Sigma}_{\mathbf{x}\mathbf{x}}$ are estimates of
$\Sigma_{D\mathbf{x},\mathbf{x}}$ and $\Sigma_{\mathbf{x}\mathbf{x}}$,
respectively. There are many choices of these estimators, some of which
include tuning parameter(s). For example, one can use discretizing techniques
and finite differences to estimate them. We call this approach the _simple
two-stage method_. Alternatively, we could use roughness penalized basis
splines to estimate $\hat{\mathbf{x}}(t)$, then apply the differential
operator and integral operator to estimate those terms. We call the latter
approach the _functional two-stage method_. These two methods are described in
detail in Section S7, Supplementary Text. In either case,
$\hat{\Sigma}_{D\mathbf{x},\mathbf{x}}$ and
$\hat{\Sigma}_{\mathbf{x}\mathbf{x}}$ could be represented by the following
matrix operations

$\hat{\Sigma}_{\mathbf{x}\mathbf{x}}=YSY^{\prime},\quad\hat{\Sigma}_{D\mathbf{x},\mathbf{x}}=YLY^{\prime},\quad\hat{A}=YLY^{\prime}(YSY^{\prime})^{-1}.$
(29)

Here $S$ and $L$ were two $(n\times n)$-dimensional matrices obtained from the
particular estimating procedure, such as the smoothing step in the functional
two-stage methods.

The following two theorems state that: a) when there is no error in estimating
$\mathbf{x}(t)$, $\hat{A}$ defined in Equation (28) is exact, and b) small
errors in estimating $\mathbf{x}(t)$ and $D\mathbf{x}(t)$ only induce a small
error in $\hat{A}$ for an identifiable $(A,\mathbf{x}_{0})$.

######  Theorem 3.3.

Assume that a linear ODE system with constant coefficient is identifiable at
$\mathbf{x}_{0}$. A matrix $A\in M_{d\times d}$ is its system matrix if and
only if it satisfies

$A=\Sigma_{D\mathbf{x},\mathbf{x}}\cdot\Sigma_{\mathbf{x}\mathbf{x}}^{-1}.$
(30)

######  Theorem 3.4.

Let $\hat{\mathbf{x}}(t)$ and $D\hat{\mathbf{x}}(t)$ be the estimates of the
solution trajectory and its derivative used in Equation (28). Let
$\delta_{1}:=\|\hat{\mathbf{x}}(t)-\mathbf{x}(t)\|_{L^{2}}$ and
$\delta_{2}:=\|D\hat{\mathbf{x}}(t)-D\mathbf{x}(t)\|_{L^{2}}$ be the
estimation errors measured in $L^{2}$ norm of $\hat{\mathbf{x}}(t)$ and
$D\hat{\mathbf{x}}(t)$, respectively; and define
$\delta:=\max(\delta_{1},\delta_{2})$. We have

$\|\hat{A}-A\|_{F}\leqslant C\delta.$ (31)

Here $C$ is a multiplicative constant that depends on
$\|\hat{\mathbf{x}}(t)\|$, $\|D\hat{\mathbf{x}}(t)\|$, and the condition
number of $\Sigma_{\mathbf{x}\mathbf{x}}$.

###### Proof.

The proofs of the above two theorems are provided in Section S4, Supplementary
Text. ∎

It is well known that, with reasonable knot placement and design points
($t_{j}$s), $\hat{x}(t)$ and $D\hat{x}(t)$ obtained by roughness penalized
smoothing splines converge to $\mathbf{x}(t)$ and $D\mathbf{x}(t)$ in $L^{2}$
norm. Given Theorem 3.3, it is reasonable to assume that $\hat{A}-A$ is small
in our subsequent analyses.

####  3.3.3 Smoothed Condition Number

We first propose a straightforward generalization of Stanhope’s condition
number, called the smoothed condition number (SCN), to measure practical
identifiability:

$\tau(Y,S):=\kappa(\hat{\Sigma}_{\mathbf{x}\mathbf{x}})=\kappa(YSY^{\prime}).$
(32)

Apparently, Equation (28) is the main motivation of this generalization.
Compared with Stanhope’s $\kappa$ statistic, SCN incorporates the information
contained in the smoothing operator $S$, therefore captures more information
of the parameter estimation procedure.

####  3.3.4 Practical Identifiability Score

Based on Equation (29), a small perturbation of data could induce a small
$\|\hat{A}-A\|_{F}^{2}$. To conduct a formal sensitivity analysis, we need to
make the following additional assumptions:

  1. 1.

Measurement errors are uncorrelated:
$\mathrm{cor}(\epsilon_{ij},\epsilon_{i^{\prime}j^{\prime}})=0$ if $i\neq
i^{\prime}$ or $j\neq j^{\prime}$.

  2. 2.

These errors are relatively small, namely,
$\mathrm{var}(\epsilon_{ij})\leqslant\sigma^{2}\ll 1$ for all $i,j$.

  3. 3.

$\|\hat{A}(X)-A\|_{F}^{2}\ll\|\hat{A}(Y)-A\|_{F}^{2}$, namely, the _numerical
error_ in $\hat{A}$ due to the use of discrete data is much smaller than the
_variance_ of $\hat{A}$ caused by measurement error. This assumption can also
be expressed as $A\approx XLX^{\prime}(XSX^{\prime})^{-1}$, which is a
reasonable assumption for cases in which $n$ is large based on Theorem 3.3.

Denote $\bm{\epsilon}_{S}=\bm{\epsilon}SX^{\prime}+XS\bm{\epsilon}^{\prime}$,
$\bm{\epsilon}_{L}=\bm{\epsilon}LX^{\prime}+XL\bm{\epsilon}^{\prime}$, and
$N=\left(XSX^{\prime}\right)^{-1}$. Based on the above assumptions, we have

$\begin{split}\hat{A}(X+\bm{\epsilon})-A&=(X+\bm{\epsilon})L(X+\bm{\epsilon})^{\prime}\left((X+\bm{\epsilon})S(X+\bm{\epsilon})^{\prime}\right)^{-1}-A\\\
&\approx(XLX^{\prime}+\bm{\epsilon}LX^{\prime}+XL\bm{\epsilon}^{\prime})\left(XSX^{\prime}+\bm{\epsilon}SX^{\prime}+XS\bm{\epsilon}^{\prime}\right)^{-1}-A\\\
&\approx(XLX^{\prime}+\bm{\epsilon}_{L})\left(N-N\bm{\epsilon}_{S}N\right)-A\\\
&\approx-XLX^{\prime}N(\bm{\epsilon}_{S})N+\bm{\epsilon}_{L}N\\\
&\approx(\bm{\epsilon}_{L}-A\bm{\epsilon}_{S})N.\end{split}$ (33)

Using Proposition S1, we have

$\small\begin{gathered}\begin{split}E(\bm{\epsilon}_{L}^{\prime}\bm{\epsilon}_{L})&=E((\bm{\epsilon}L^{\prime}X^{\prime}+XL^{\prime}\bm{\epsilon}^{\prime})(\bm{\epsilon}LX^{\prime}+XL\bm{\epsilon}^{\prime}))\\\
&=E\left(\bm{\epsilon}L^{\prime}X^{\prime}\bm{\epsilon}LX^{\prime}+\bm{\epsilon}L^{\prime}X^{\prime}XL\bm{\epsilon}^{\prime}+XL^{\prime}\bm{\epsilon}^{\prime}\bm{\epsilon}LX^{\prime}+XL^{\prime}\bm{\epsilon}^{\prime}XL\bm{\epsilon}^{\prime}\right)\\\
&=\sigma^{2}\left(XL^{2}X^{\prime}+\mathrm{tr}(L^{\prime}X^{\prime}XL)\cdot
I_{d}+d\cdot XL^{\prime}LX^{\prime}+X(L^{\prime})^{2}X^{\prime}\right).\\\
\end{split}\\\
\begin{split}E\left(\bm{\epsilon}_{S}^{\prime}A^{\prime}\bm{\epsilon}_{L}\right)&=E((\bm{\epsilon}S^{\prime}X^{\prime}+XS^{\prime}\bm{\epsilon}^{\prime})A^{\prime}(\bm{\epsilon}LX^{\prime}+XL\bm{\epsilon}^{\prime}))\\\
&=E\left(\bm{\epsilon}S^{\prime}X^{\prime}A^{\prime}\bm{\epsilon}LX^{\prime}+XS^{\prime}\bm{\epsilon}^{\prime}A^{\prime}\bm{\epsilon}LX^{\prime}+\bm{\epsilon}S^{\prime}X^{\prime}A^{\prime}XL\bm{\epsilon}^{\prime}+XS^{\prime}\bm{\epsilon}^{\prime}A^{\prime}XL\bm{\epsilon}^{\prime}\right)\\\
&\leqslant\sigma^{2}\left(AXSLX^{\prime}+\mathrm{tr}(A)\cdot
XS^{\prime}LX^{\prime}+\mathrm{tr}(S^{\prime}X^{\prime}A^{\prime}XL)\cdot
I_{d}+XS^{\prime}L^{\prime}X^{\prime}A\right).\end{split}\\\
\begin{split}E\left(\bm{\epsilon}_{S}^{\prime}A^{\prime}A\bm{\epsilon}_{S}\right)&=E((\bm{\epsilon}S^{\prime}X^{\prime}+XS^{\prime}\bm{\epsilon}^{\prime})A^{\prime}A(\bm{\epsilon}SX^{\prime}+XS\bm{\epsilon}^{\prime}))\\\
&=E\big{(}\bm{\epsilon}S^{\prime}X^{\prime}A^{\prime}A\bm{\epsilon}SX^{\prime}+XS^{\prime}\bm{\epsilon}^{\prime}A^{\prime}A\bm{\epsilon}SX^{\prime}\\\
&\quad+\bm{\epsilon}S^{\prime}X^{\prime}A^{\prime}AXS\bm{\epsilon}^{\prime}+XS^{\prime}\bm{\epsilon}^{\prime}A^{\prime}AXS\bm{\epsilon}^{\prime}\big{)}\\\
&\leqslant\sigma^{2}\big{(}A^{\prime}AXS^{2}X^{\prime}+\mathrm{tr}(A^{\prime}A)\cdot
XS^{\prime}SX^{\prime}\\\
&\quad\mathrm{tr}\left(S^{\prime}X^{\prime}A^{\prime}AXS\right)\cdot
I_{d}+X(S^{\prime})^{2}X^{\prime}A^{\prime}A\big{)}.\end{split}\end{gathered}$
(34)

Therefore

$\begin{split}\|\hat{A}(X+\bm{\epsilon})-A\|_{F}^{2}&\approx\mathrm{tr}\left((\bm{\epsilon}_{L}-A\bm{\epsilon}_{S})N^{2}(\bm{\epsilon}_{L}^{\prime}-\bm{\epsilon}_{S}^{\prime}A^{\prime})\right)\\\
&=\mathrm{tr}\left((\bm{\epsilon}_{L}^{\prime}-\bm{\epsilon}_{S}^{\prime}A^{\prime})(\bm{\epsilon}_{L}-A\bm{\epsilon}_{S})N^{2}\right)\\\
&=\mathrm{tr}\left((\bm{\epsilon}_{L}^{\prime}\bm{\epsilon}_{L}-\bm{\epsilon}_{S}^{\prime}A^{\prime}\bm{\epsilon}_{L}-\bm{\epsilon}_{L}^{\prime}A\bm{\epsilon}_{S}+\bm{\epsilon}_{S}^{\prime}A^{\prime}A\bm{\epsilon}_{S})N^{2}\right)\\\
E\|\hat{A}(X+\bm{\epsilon})-A\|_{F}^{2}&\approx\mathrm{tr}\big{(}(E(\bm{\epsilon}_{L}^{\prime}\bm{\epsilon}_{L})-2E\left(\bm{\epsilon}_{S}^{\prime}A^{\prime}\bm{\epsilon}_{L}\right)+E\left(\bm{\epsilon}_{S}^{\prime}A^{\prime}A\bm{\epsilon}_{S}\right))N^{2}\big{)}\\\
&\leqslant\sigma^{2}W(X).\\\ \end{split}$ (35)

Here

$\begin{split}W(X|A,S,L)&:=\mathrm{tr}\big{(}(XSX^{\prime})^{-2}\big{(}XL^{2}X^{\prime}+d\cdot
XL^{\prime}LX^{\prime}+X(L^{\prime})^{2}X^{\prime}\\\
&\quad-2AXSLX^{\prime}-2\mathrm{tr}(A)\cdot
XS^{\prime}LX^{\prime}-2XS^{\prime}L^{\prime}X^{\prime}A\\\
&\quad+A^{\prime}AXS^{2}X^{\prime}+\mathrm{tr}(A^{\prime}A)\cdot
XS^{\prime}SX^{\prime}+X(S^{\prime})^{2}X^{\prime}A^{\prime}A)\\\
&\quad+\mathrm{tr}(L^{\prime}X^{\prime}XL-2S^{\prime}X^{\prime}A^{\prime}XL+S^{\prime}X^{\prime}A^{\prime}AXS)\cdot
I_{d}\big{)}\\\
&=\mathrm{vec}((XS^{\prime}X^{\prime})^{-2})^{\prime}\mathrm{vec}\big{(}XL^{2}X^{\prime}+d\cdot
XL^{\prime}LX^{\prime}+X(L^{\prime})^{2}X^{\prime}\\\
&\quad-2AXSLX^{\prime}-2\mathrm{tr}(A)\cdot
XS^{\prime}LX^{\prime}-2XS^{\prime}L^{\prime}X^{\prime}A\\\
&\quad+A^{\prime}AXS^{2}X^{\prime}+\mathrm{tr}(A^{\prime}A)\cdot
XS^{\prime}SX^{\prime}+X(S^{\prime})^{2}X^{\prime}A^{\prime}A\big{)}\\\
&\quad+\mathrm{tr}(L^{\prime}X^{\prime}XL-2S^{\prime}X^{\prime}A^{\prime}XL+S^{\prime}X^{\prime}A^{\prime}AXS)\mathrm{tr}(N^{2}).\end{split}$
(36)

By construction, $W(X)$ is a scalar that quantifies the MSE of $\hat{A}$ as a
function of $\sigma^{2}$, the maximum variance of $\epsilon_{ij}$ for all
$i=1,2,\dots,d$ and $j=1,2,\dots,n$. Smaller values of $W(X)$ imply better
practical identifiability in reconstructing $A$. Motivated by this fact, we
define the practical identifiability score (PIS) as the sample version of
$W(X)$. Specifically, PIS (denoted as $w^{*}$ in Equation (37)) is computed by
replacing $X$ and $A$ in Equation (36) with $Y$ and
$\hat{A}:=YLY^{\prime}(YSY^{\prime})^{-1}$, respectively:

$w^{*}(Y|S,L):=W(Y|\hat{A},S,L).$ (37)

Compared with Stanhope’s $\kappa$ and SCN, PIS depends not only on the
observed data ($Y$), but also the $S$ and $L$ matrix of the particular two-
stage method used in reconstructing $\hat{A}$, therefore it is a more accurate
indicator of practical identifiability. A simulation study was designed to
demonstrate this point in Section 4.2.

##  4 Simulation studies

###  4.1 ICIS is Inversely Correlated with the Relative Estimation Error (REE)

Figure 3: ICIS ($w_{0}^{*}$) is inversely correlated with the relative
estimation error (REE). In all four subfigures, the $x$-axis is $1/w_{0}^{*}$,
$y$-axis is REE. Each dot represents one of the 100 repetitions of SIM1. Top
two subfigures are generated from noisy data ($Y$); bottom two subfigures uses
noise-free data ($X$). The right subfigures are the zoomed-in version of the
left subfigures.

We design SIM1 to demonstrate that the ICIS is inversely correlated with
estimation error, as predicted in Section 3.2. In this simulation, the system
matrix and its two invariant subspaces are:

$A=\begin{pmatrix}-0.1&3&0\\\ -3&-0.1&0\\\ 0&0&-0.5\end{pmatrix},\qquad
L_{1}=\mathrm{span}\left(\mathbf{e}_{1},\mathbf{e}_{2}\right),\quad
L_{2}=\mathrm{span}\left(\mathbf{e}_{3}\right).$

Here $\mathbf{e}_{k}$ is the $k$th column vector of $I_{3}$, a.k.a. the $k$th
natural basis vector of $\mathbb{R}^{3}$. We generate $\mathbf{x}_{0}$ from
$N(0_{3},I_{3})$ first, then standardize it to have unit length to reduce the
variation in ICIS due to different $|\mathbf{x}_{0}|$. This is equivalent to
sampling $\mathbf{x}_{0}$ from a uniform distribution on $S^{2}$.

Once $\mathbf{x}_{0}$ is generated, we compute $X_{\cdot
j}:=e^{t_{j}A}\mathbf{x}_{0}$ at a time grid with range $[0,6]$ and step
$\Delta t=0.1$, i.e., $t_{j}=0,0.1,\dots,6$ for $j=1,\dots,61$. We also add a
small noise $\epsilon_{ij}\sim N(0,\sigma^{2})$, $\sigma=0.05$ to each
observation to create noisy data $Y=X+\bm{\epsilon}$. A functional two-stage
method based on cubic splines with roughness penalty $\lambda=0.001$ is used
to estimate $\hat{A}$ from both noisy ($Y$) and noise-free ($X$) data.

The accuracy of estimation is measured by relative estimation error (REE),
defined as follows

$\mathrm{REE}(\hat{A},A):=\frac{\|\hat{A}-A\|_{F}}{\|A\|_{F}}.$ (38)

We repeat SIM1 for 100 times, each with randomly generated $\mathbf{x}_{0}$
and measurement error. We find that ICIS ($w_{0}^{*}$) is strongly negatively
correlated with REE. The Spearman correlation between these two quantities is
$\rho_{1}=-0.803$ for the noisy data and $\rho_{2}=-0.843$ for the noise-free
data. This inverse correlation is visualized in Figure 3. Other than a few
outliers, $1/w_{0}^{*}$ has an almost perfect linear relationship with
$\mathrm{REE}$ in the noise-free case (the second row of Figure 3). The
correlation between $1/w_{0}^{*}$ and $\mathrm{REE}$ is weaker but still quite
apparent for the noisy data (the first row of Figure 3).

###  4.2 Using SCN and PIS to Classify Identifiable and Un-identifiable
Systems

We design SIM2 to demonstrate that, when data collected at all time points are
available, SCN and PIS have better performance in classifying identifiable and
un-identifiable Systems than ICIS and Stanhope’s $\kappa$. SIM2 contains one
identifiable case and two unidentifiable cases, which are described as
follows.

  1. 1.

In each one of 200 repetitions, we generate two $4\times 4$-dimensional system
matrices $A$ and $B$, and two initial conditions $\mathbf{x}_{0}^{(a)}$ and
$\mathbf{x}_{0}^{(b)}$ on $S^{3}$.

  2. 2.

Both $A$ and $B$ have one pair of complex eigenvalues and two real
eigenvalues. The eigenvalues of $A$,
$(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4})$, are generated in this way

$\begin{gathered}\lambda_{1},\lambda_{2}=-0.1\pm bi,\quad
b\sim\mathrm{Unif}([2,4]).\\\
\lambda_{3}\sim\mathrm{Unif}([-0.8,-0.4]),\quad\lambda_{4}\sim\mathrm{Unif}([-2,-1.2]).\end{gathered}$
(39)

The eigenvalues of $B$ are set to be
$(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{3})$, namely, $B$ has a pair of
_repeated eigenvalues_ by construction. Therefore it is not identifiable at
any initial condition.

  3. 3.

After generating the eigenvalues, we sample an orthogonal matrix $Q$ from the
standardized Haar measure (the uniform distribution) on the orthogonal group
$O(4)$, and create two system matrices $A$ and $B$ as follows:

$A=Q\begin{pmatrix}-0.1&b&0&0\\\ -b&-0.1&0&0\\\ 0&0&\lambda_{3}&0\\\
0&0&0&\lambda_{4}\\\ \end{pmatrix}Q^{\prime},\quad
B=Q\begin{pmatrix}-0.1&b&0&0\\\ -b&-0.1&0&0\\\ 0&0&\lambda_{3}&0\\\
0&0&0&\lambda_{3}\\\ \end{pmatrix}Q^{\prime}.$ (40)

  4. 4.

Like SIM1, we sample $\mathbf{x}_{0}^{(a)}$ from the uniform distribution on
$S^{d-1}$ ($d=4$ in this case). Furthermore, we only keep those
$\mathbf{x}_{0}^{(a)}$ with relatively large ICIS, namely
$w_{0}^{*}(A,\mathbf{x}_{0}^{(a)})>0.2$. This ensures the practical
identifiability of $A$ at $\mathbf{x}_{0}^{(a)}$.

  5. 5.

Once $\mathbf{x}_{0}^{(a)}$ is sampled, we define

$\mathbf{x}_{0}^{(b)}=\frac{(I_{4}-Q_{\cdot 4}Q_{\cdot
4}^{\prime})\mathbf{x}_{0}^{(a)}}{\big{|}(I_{4}-Q_{\cdot 4}Q_{\cdot
4}^{\prime})\mathbf{x}_{0}^{(a)}\big{|}}.$ (41)

By construction, $\mathbf{x}_{0}^{(b)}$ is a unit vector such that
$\mathbf{x}_{0}^{(b)}\perp Q_{\cdot 4}$, so that ICIS equals zero in this
case. According to our theoretical derivations, $A$ is not identifiable at
$\mathbf{x}_{0}^{(b)}$ due to ill-positioned initial conditions.

  6. 6.

We compute three sets of solution trajectories: case (A) corresponds with
$(A,\mathbf{x}_{0}^{(a)})$, case (B) with $(A,\mathbf{x}_{0}^{(b)})$, and case
(C) with $(B,\mathbf{x}_{0}^{(a)})$.

For each case, we compute ICIS ($w_{0}^{*}$), SCN ($\tau$), PIS ($w^{*}$), and
Stanhope’s $\kappa$, based on both noisy and noise-free data. The results are
illustrated in Figures 4 and 5. We find that for noise-free data, SCN, PIS,
and $\kappa$ perform very well, with almost perfect area under the curve (AUC)
in receiver operating characteristic (ROC) analyses. However, ICIS only has a
relatively small AUC=0.773. This is not a surprise at all because ICIS is
designed to detect un-identifiability issues associated with ill-positioned
_initial conditions_ (case B), not un-identifiable systems that have _repeated
eigenvalues_ (case C). This fact is also revealed in the corresponding boxplot
in Figure 4 (second column).

For data with noise, $\kappa$ is almost uninformative (AUC=0.503), but SCN,
which is a smoothed extension of $\kappa$, works very well (AUC=0.946). It
suggests that taking the smoothing effect into the consideration in SCN
improves its utility as a classifier of identifiable systems.

While SCN has significantly better performance than $\kappa$ and ICIS, it is
still an _ad hoc_ metric of practical identifiability that does not account
for the uncertainty in $\hat{A}$ due to measurement error. In contrast, PIS is
designed based on rigorous asymptotic analysis on the variance of $\hat{A}$,
therefore PIS has the best performance (AUC=0.962). That being said, we need
to point out that from the computational perspective, SCN is more efficient
and numerically robust, because SCN does not contain $(YSY^{\prime})^{-1}$ and
$(YSY^{\prime})^{-2}$ terms used in PIS, which could have numerical issues if
the dimension of the ODE system is large. In summary, SCN could be considered
as a simplified version of PIS that is less vulnerable to computational
issues.

Figure 4: Values of four identifiability measures (ICIS, SCN, PIS, and
Stanhope’s $\kappa$) computed from 200 repetitions of SIM2. For all four
scores, smaller values imply better identifiability. The first row of
subfigures are computed from noisy data ($Y$), the second row of subfigures
are computed from noise-free data ($X$). Each subfigure has three cases: (A)
is the identifiable case generated by $(A,\mathbf{x}_{0}^{(a)})$; (B) is the
unidentifiable case generated by $(A,\mathbf{x}_{0}^{(b)})$; (C) is the
unidentifiable case generated by $(B,\mathbf{x}_{0}^{(b)})$.  Figure 5: ROC
curves of four identifiability measures (ICIS, SCN, PIS, and Stanhope’s
$\kappa$) in classifying case A (identifiable systems) from cases B and C
(both are unidentifiable systems) in SIM2, with 200 repetitions. The left
panel use data with noise; the right panel use noise-free data.

Both noisy and noise-free data for all three cases were illustrated in Figure
6. Notably, visual examinations did not reveal apparent differences between
the three cases, suggesting that the identifiability of the ODE system does
not depend on obvious features in the solution trajectories.

Figure 6: An illustration of the fitted solution curves in SIM2. Dots
represent noisy data $y_{ij}$, four colors represent four dimensions. Solid
curves are the noise-free solutions ($x(t)$) and dotted curves are the
smoothed curves ($\hat{x}(t)$) fitted by roughness penalized splines. The same
roughness penalty ($\lambda=0.001$) was used in all three cases. Overall, the
fitted curves agree with the noise-free data well.

##  5 Conclusions

Classical identifiability analyses for ODE systems typically depend on the
availability of solution trajectories from arbitrarily many initial points.
However, in many real world problems, the system matrix must be estimated from
just one observed trajectory. In this case, identifiability depends not only
on the properties of $A$, but also the initial condition $\mathbf{x}_{0}$. In
this case, the $(A,\mathbf{x}_{0})$-identifiability used in our study is more
appropriate than classical identifiability measures.

We develop an explicit formula of all matrices that are unidentifiable with
$A$ at a given $\mathbf{x}_{0}$ in this study. It enables researchers to gain
better insight into identifiability analysis and help them design more
practical simulation studies.

Another notable finding of our study is that when $A$ is coupled (not
diagonal), an identifiability issue in just a one-dimensional invariant
subspace could cause issues in many other elements of $A$ (e.g., Example 1).
Consequently, identifiability analyses that only depend on the topology of the
network are insufficient in practice.

For high-dimensional cases, even if $\mathbf{x}_{0}$ is generated in a
“completely random” fashion (e.g., $\mathbf{x}_{0}\sim N(0_{d},I_{d})$), by
chance, one invariant subspace of $A$ may have very little information, which
in turn leads to practical identifiability issues. In fact, we are able to
prove that when $d\to\infty$, $\mathrm{ICIS}(A,\mathbf{x}_{0})\to 0$ for a
large class of random ODE systems, which suggests that the practical
identifiability properties of low-dimensional and high-dimensional systems are
fundamentally different. We believe it will be rewarding to derive more
accurate convergence rates for $\mathrm{ICIS}$ as a function of $d$ in a
future study. It will require combining advanced techniques in random matrix
theory, especially for ensembles of asymmetric matrices (e.g., Conjecture S2.1
in Supplementary Text) in which the $Q$ matrix are no longer orthogonal, with
the identifiability analysis of ODE systems.

In this study, we also developed two scores, SCN and PIS, that use the entire
dataset obtained at all time points, to quantify the practical identifiability
for real world applications. Both SCN and PIS are more accurate than
Stanhope’s $\kappa$ when noise is present in the data, as shown by extensive
simulation studies.

While our methods are developed for homogeneous systems, it should be
relatively easy to generalize them for the following inhomogeneous linear ODE
system

$\begin{cases}D\mathbf{x}(t)=A\mathbf{x}(t)+b,\quad
b\in\mathbb{R}^{d},&t\in(0,T],\\\
\mathbf{x}(0)=\mathbf{x}_{0}\in\mathbb{R}^{d}.\end{cases}$ (42)

This is because Equation (42) can be transformed into an _equivalent
homogeneous system_ with a simple mathematical technique. Let
$\mathbf{z}(t)=(\mathbf{x}(t),1)^{\prime}=(x_{1}(t),\dots,x_{d}(t),\,1)^{\prime}\in\mathbb{R}^{d+1}$.
It satisfies the following ODE

$\begin{cases}D\mathbf{z}(t)=\breve{A}\mathbf{z}(t),\qquad t\in(0,T],\\\
\mathbf{z}(0)=(\mathbf{x}_{0},1)^{\prime}\in\mathbb{R}^{d+1}.\end{cases}\qquad\breve{A}:=\begin{pmatrix}A&b\\\
0_{d}&0\end{pmatrix}.$ (43)

Therefore, the identifiability of Equation (42) is the same as the
identifiability of Equation (43), which is a homogeneous equation with the
constraint that the last row of $\breve{A}$ must be zeros. Let
$M_{0}=\left\\{B\in M_{(d+1)\times(d+1)}:\,M_{(d+1)\cdot}=0_{d+1}\right\\}$ be
the set of $(d+1)\times(d+1)$-dimensional matrices such that their last rows
equal $0_{d+1}$. The unidentifiability class associated with system
$(A,b,\mathbf{x}_{0})$, denoted by $[A,b]_{\mathbf{x}_{0}}$, is the following
subset of $M_{d\times d}$:

$[A,b]_{\mathbf{x}_{0}}=[\breve{A}]_{(\mathbf{x}_{0},1)^{\prime}}\cap M_{0}.$
(44)

More future work is required to extend ICIS, SCN, and PIS for constrained
systems, so that they can be used as practical guidance for applications with
a priori information.

In the near future, we plan to extend our work to the following family of
nonlinear ODE system:

$D\mathbf{x}(t)=Af\left(\mathbf{x}(t)\right),\qquad\mathbf{x}(0)=\mathbf{x}_{0}.$
(45)

Here $f(\cdot):\mathbb{R}^{d_{1}}\to\mathbb{R}^{d_{2}}$ is a known locally
Lipschitz function of $\mathbf{x}(t)$, $A\in M_{d_{1}\times d_{2}}$ is the
system matrix that needs to be estimated. This system has been studied by
Stanhope and colleagues, and their main conclusion (Theorem (5.3) in [20]) is
very similar to that for the linear ODE systems: $A$ is identifiable at
$\mathbf{x}_{0}$ if and only if the solution curve is not confine in a proper
linear subspace of $\mathbb{R}^{d_{2}}$. To extend the SCN and PIS we
developed in this study to Equation (45), we will need to study the
sensitivity of an extended two-stage method that works for Equation (45).

Using linearization techniques, we believe SCN and PIS can be further extended
to other types of nonlinear systems. To this end, we need: a) to approximate a
nonlinear ODE system by a linear ODE at $\mathbf{x}(t)$; b) to propose a local
version of the $(A,\mathbf{x}_{0})$-identifiability that works in a
neighborhood of $A$ at $\mathbf{x}(t)$; c) to study the sensitivity of a
reasonable parameter estimator for such system, and propose an identifiability
score based on the useful information aggregated from all time points.

## References

  * [1] J. Butcher, Ordinary differential equations, in: Walter Gautschi, Vol. 3, Springer, 2014, pp. 7–8.
  * [2] D. Commenges, D. Jolly, J. Drylewicz, H. Putter, R. Thiébaut, Inference in HIV dynamics models via hierarchical likelihood, Computational Statistics & Data Analysis 55 (1) (2011) 446–456.
  * [3] H. De Jong, Modeling and simulation of genetic regulatory systems: a literature review, Journal of Computational Biology 9 (1) (2002) 67–103.
  * [4] P. W. Hemker, Numerical methods for differential equations in system simulation and in parameter estimation, Analysis and Simulation of Biochemical Systems 28 (1972) 59–80.
  * [5] N. S. Holter, A. Maritan, M. Cieplak, N. V. Fedoroff, J. R. Banavar, Dynamic modeling of gene expression data, Proceedings of the National Academy of Sciences 98 (4) (2001) 1693–1698.
  * [6] Y. Huang, D. Liu, H. Wu, Hierarchical bayesian methods for estimation of parameters in a longitudinal HIV dynamic system, Biometrics 62 (2) (2006) 413–423.
  * [7] M. Lavielle, A. Samson, A. Karina Fermin, F. Mentré, Maximum likelihood estimation of long-term HIV dynamic models and antiviral response, Biometrics 67 (1) (2011) 250–259.
  * [8] Z. Li, P. Li, A. Krishnan, J. Liu, Large-scale dynamic gene regulatory network inference combining differential equation models with local dynamic bayesian network analysis, Bioinformatics 27 (19) (2011) 2686–2691.
  * [9] T. Lu, H. Liang, H. Li, H. Wu, High-dimensional ODEs coupled with mixed-effects modeling techniques for dynamic gene regulatory network identification, Journal of the American Statistical Association 106 (496) (2011) 1242–1258.
  * [10] J. O. Ramsay, G. Hooker, D. Campbell, J. Cao, Parameter estimation for differential equations: a generalized smoothing approach (with discussion), Journal of the Royal Statistical Society 69 (5) (2007) 741–796.
  * [11] C. Moler, C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM review 45 (1) (2003) 3–49.
  * [12] Y. Huang, H. Wu, A Bayesian approach for estimating antiviral efficacy in HIV dynamic models, Journal of Applied Statistics 33 (2) (2006) 155–174.
  * [13] Y. Huang, H. Wu, E. P. Acosta, Hierarchical Bayesian inference for HIV dynamic differential equation models incorporating multiple treatment factors, Biometrical Journal 52 (4) (2010) 470–486.
  * [14] Z. Li, M. R. Osborne, T. Prvan, Parameter estimation of ordinary differential equations, IMA Journal of Numerical Analysis 25 (2) (2005) 264–285.
  * [15] H. Putter, S. Heisterkamp, J. Lange, F. De Wolf, A Bayesian approach to parameter estimation in HIV dynamical models, Statistics in Medicine 21 (15) (2002) 2199–2214.
  * [16] L. Wu, X. Qiu, Y.-x. Yuan, H. Wu, Parameter estimation and variable selection for big systems of linear ordinary differential equations: A matrix-based approach, Journal of the American Statistical Association 114 (526) (2019) 657–667.
  * [17] H. Xue, H. Miao, H. Wu, Sieve estimation of constant and time-varying coefficients in nonlinear ordinary differential equation models by considering both numerical error and measurement error, Annals of Statistics 38 (4) (2010) 2351–2387.
  * [18] H. Miao, X. Xia, A. S. Perelson, H. Wu, On identifiability of nonlinear ode models and applications in viral dynamics, SIAM review 53 (1) (2011) 3–39.
  * [19] A. THOWSEN, Identifiability of dynamic systems, International Journal of Systems Science 9 (7) (1978) 813–825.
  * [20] S. Stanhope, J. Rubin, D. Swigon, Identifiability of linear and linear-in-parameters dynamical systems from a single trajectory, SIAM Journal on Applied Dynamical Systems 13 (4) (2014) 1792–1815.
  * [21] X. Qiu, S. Wu, S. P. Hilchey, J. Thakar, Z.-P. Liu, S. L. Welle, A. D. Henn, H. Wu, M. S. Zand, Diversity in compartmental dynamics of gene regulatory networks: the immune response in primary influenza a infection in mice, PloS one 10 (9) (2015).
  * [22] X. Sun, F. Hu, S. Wu, X. Qiu, P. Linel, H. Wu, Controllability and stability analysis of large transcriptomic dynamic systems for host response to influenza infection in human, Infectious Disease Modelling 1 (1) (2016) 52–70.
  * [23] S. Wu, Z.-P. Liu, X. Qiu, H. Wu, High-dimensional ordinary differential equation models for reconstructing genome-wide dynamic regulatory networks, in: Topics in applied statistics, Springer, New York, NY, 2013, pp. 173–190.
  * [24] S. Wu, Z.-P. Liu, X. Qiu, H. Wu, Modeling genome-wide dynamic regulatory network in mouse lungs with influenza infection using high-dimensional ordinary differential equations, PloS one 9 (5) (2014).
  * [25] J. A. McCullers, J. L. McAuley, S. Browall, A. R. Iverson, K. L. Boyd, B. Henriques Normark, Influenza enhances susceptibility to natural acquisition of and disease due to streptococcus pneumoniae in ferrets, The Journal of infectious diseases 202 (8) (2010) 1287–1295.
  * [26] E. Tunali, T.-J. Tarn, New results for identifiability of nonlinear systems, IEEE Transactions on Automatic Control 32 (2) (1987) 146–154.
  * [27] A. M. Jeffrey, X. Xia, I. Craig, Identifiability of hiv/aids models, Deterministic and Stochastic models of AIDS epidemics and HIV infections with intervention (2005) 255–286.
  * [28] X. Xia, C. H. Moog, Identifiability of nonlinear systems with application to hiv/aids models, IEEE transactions on automatic control 48 (2) (2003) 330–336.
  * [29] J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6 (1965) 440.
  * [30] I. Gohberg, P. Lancaster, L. Rodman, Invariant subspaces of matrices with applications, SIAM, 2006.
  * [31] N. Lehmann, H.-J. Sommers, Eigenvalue statistics of random real matrices, Physical Review Letters 67 (8) (1991) 941–944.
  * [32] T. Tao, Topics in random matrix theory, Vol. 132, American Mathematical Society Providence, RI, 2012.
  * [33] G. Livan, M. Novaes, P. Vivo, Introduction to random matrices: theory and practice, Vol. 26, Springer, 2018.
  * [34] H. Weyl, The classical groups: their invariants and representations, Vol. 45, Princeton university press, 1946.
  * [35] A. Edelman, The probability that a random real gaussian matrix haskreal eigenvalues, related distributions, and the circular law, Journal of Multivariate Analysis 60 (2) (1997) 203–232.
  * [36] T. S. Ferguson, A course in large sample theory, Vol. 49, Chapman & Hall London, 1996.

Supplementary Text: Identifiability Analysis of Linear Ordinary Differential
Equation Systems with a Single Trajectory

##  S1 Structural Identifiability is Unattainable for Linear ODE Systems

In this section, we prove the following statement:

######  Proposition S1.

If the dimension $d$ is odd, there is no open subset $\Omega\subseteq
M_{d\times d}$ on which ODE system (3) is unconditionally identifiable. If the
dimension $d$ is even, system (3) is not unconditionally identifiable for all
$\Omega\subseteq M_{d\times d}$ such that $\lambda\left(M_{d\times
d}\setminus\Omega\right)=0$, where $\lambda(\cdot)$ is the Lebesgue measure on
$M_{d\times d}$.

###### Proof.

First, we note that the structural identifiability defined in 1.3 is a special
case of the so-called _unconditional identifiability_ defined in Definition
2.4 [20] when $\Omega$ is set to be an open and dense subset of $M_{d\times
d}$.

According to Corollary 3.9 in [20], ODE system (3) is unconditionally
identifiable on an open set $\Omega\subset M_{d\times d}$ iff for every
$A\in\Omega$, there is no left-eigenvector of $A$ that is orthogonal to every
$\mathbf{x}_{0}\in\mathbb{R}^{d}\setminus\left\\{0_{d}\right\\}$. That
immediately excludes matrices that has at least one real eigenvalue and
eigenvector, which includes all cases when $d$ is odd.

Now let us focus on the even-dimensional cases. Edelman showed in [35] that
for a random matrix $A\in M_{d\times d}$ with $i.i.d.$ normally distributed
entries (the Ginibre ensemble), the probability of $A$ having a real
eigenvalue is strictly greater than zero. Because the probability measure of
the Ginibre ensemble and the Lebesgue measure on $M_{d\times d}$ are
absolutely continuous with respect to each other, we know that we cannot find
$\Omega$ such that: a) $A$ is unconditionally identifiable on $\Omega$, and b)
$\lambda\left(M_{d\times d}\setminus\Omega\right)=0$. ∎

##  S2 Practically Unidentifiable High-dimensional ODEs

In this section, we move Theorem 3.1, which states that if the dimension is
high and system matrix $A$ is generated from a large class of random matrices,
the ICIS converges to zero in $L^{2}$ (and in probability) when $d\to\infty$.

First, we need to prove the following technical lemma.

######  Lemma S1.

Let $U=(U_{1},U_{2},\dots,U_{d})^{\prime}$, $U\sim\mathrm{Unif}(S^{d-1})$ be a
unit vector in $\mathbb{R}^{d}$ generated from the uniform distribution on
$S^{d-1}$. Let $S_{\min}:=\min_{i}U_{i}^{2}$.

We have

$\frac{2d^{3}}{\pi}\cdot S_{\min}\stackrel{{\scriptstyle
w}}{{\longrightarrow}}\mathrm{Weibull}\Big{(}1,\frac{1}{2}\Big{)}.\\\ $ (S.1)

Here
$\mathrm{Weibull}\Big{(}1,\frac{1}{2}\Big{)}(x)=1_{\mathbb{R}^{+}}(x)\cdot\left(1-e^{-x^{1/2}}\right)$
is a Weibull distribution with scale parameter 1 and shape parameter 1/2.

###### Proof.

Since $U\sim\mathrm{Unif}(S^{d-1})$, there exists $Z\sim N(0_{d},I_{d})$, such
that
$U=\frac{Z}{|Z|}=\left(\frac{Z_{1}}{|Z|},\dots,\frac{Z_{d}}{|Z|}\right)^{\prime}$.
Let $W_{\min}:=\min_{i}Z_{i}^{2}$. Because $|Z|$ and $|Z|^{2}$ are a constant
for all $i=1,2,\dots,d$, it is easy to see that
$S_{\min}=\frac{W_{\min}}{|Z|^{2}}$.

Based on the Fisher-Tippett-Gnedenko theorem [36] and notice that $W_{\min}$
has a lower bound ($W_{\min}\geqslant 0$), we can show that

$\begin{gathered}\frac{2d^{2}W_{\min}}{\pi}\stackrel{{\scriptstyle
w}}{{\longrightarrow}}\mathrm{Weibull}\Big{(}1,\frac{1}{2}\Big{)}.\end{gathered}$
(S.2)

In addition, we know that $|Z|^{2}:=\sum_{i=1}^{d}Z_{i}^{2}\sim\chi^{2}_{d}$.
Simple calculations show that when $d\to\infty$,

$\begin{gathered}E\frac{|Z|^{2}}{d}=1,\qquad\mathrm{var}\left(\frac{|Z|^{2}}{d}\right)=\frac{2}{d}\longrightarrow
0,\qquad\frac{|Z|^{2}}{d}\stackrel{{\scriptstyle
L^{2}}}{{\longrightarrow}}1.\\\ \end{gathered}$ (S.3)

Using Slutsky’s theorem, we have

$\begin{split}\frac{2d^{3}S_{\min}}{\pi}=\frac{\frac{2d^{2}W_{\min}}{\pi}}{\frac{|Z|^{2}}{d}}\stackrel{{\scriptstyle
w}}{{\longrightarrow}}\mathrm{Weibull}\Big{(}1,\frac{1}{2}\Big{)}.\end{split}$
(S.4)

∎

######  Corollary S2.

Let $\tilde{\mathbf{x}}_{0}\sim\mathrm{Unif}\left(r\cdot S^{d-1}\right)$ be a
uniformly distributed random variable on a sphere in $\mathbb{R}^{d}$ with
radius $d$, and $T\in O(d)$ be an arbitrary orthogonal matrix. Let

$S_{\min}(T,r)=\min_{i}\left(T\tilde{\mathbf{x}}_{0}\right)_{i}^{2}$

be the smallest of the squared elements in vector $T\tilde{\mathbf{x}}_{0}$.

We have

$\frac{2d^{3}}{\pi r^{2}}\cdot S_{\min}(T,r)\stackrel{{\scriptstyle
w}}{{\longrightarrow}}\mathrm{Weibull}\Big{(}1,\frac{1}{2}\Big{)}.\\\ $ (S.5)

###### Proof.

It is easy to show that (a)
$\left(T\tilde{\mathbf{x}}_{0}\right)_{i}^{2}=r^{2}\left(T\frac{\tilde{\mathbf{x}}_{0}}{r}\right)_{i}^{2}$,
therefore $S_{\min}(T,r)=r^{2}S_{\min}(T,1)$, and (b)
$\tilde{\mathbf{x}}_{0}\sim\mathrm{Unif}\left(r\cdot S^{d-1}\right)$ is
invariant under orthogonal transformation, i.e.,
$T\tilde{\mathbf{x}}_{0}\stackrel{{\scriptstyle
d}}{{=}}\tilde{\mathbf{x}}_{0}$, therefore $S_{\min}(T,1)=S_{\min}$. In
summary, $S_{\min}(T,r)\stackrel{{\scriptstyle d}}{{=}}r^{2}S_{\min}$ for
every $T\in O(d)$ and $r\in\mathbb{R}^{+}$, in which implies Equation (S.5). ∎

We are now ready to prove Theorem 3.1.

###### Proof of Theorem 3.1.

As before, we write $A=Q\Lambda Q^{-1}$ as the Jordan canonical decomposition
of $A$. Recall that ICIS is defined as
$w_{0}^{*}(A,\mathbf{x}_{0}):=\min_{k}|w_{0,k}|$, where
$w_{0,k}:=Q^{-1}\mathbf{x}_{0}$. Here $Q$ is the basis system of all invariant
subspaces of $A$. The symmetry assumption implies that: (a) $Q$ is orthogonal,
so $Q^{-1}=Q^{\prime}$, and (b) $\Lambda$ is diagonal. The orthogonal
invariance assumption implies that the marginal distribution of $Q$ must be
the standardized Haar measure on $O(d)$.

Let $T\in O(d)$ be an arbitrary orthogonal matrix and let
$\tilde{A}=TAT^{\prime}$. Apparently, its orthogonal decomposition is
$\tilde{A}=\tilde{Q}\Lambda\tilde{Q}^{\prime}$, for $\tilde{Q}=TQ$. Based on
the orthogonal invariance of $A$, we know that
$p(\tilde{Q})=p(Q)=\mathrm{Haar}(O(d))$, and
$w_{0}^{*}(\tilde{A},\mathbf{x}_{0})\stackrel{{\scriptstyle
d}}{{=}}w_{0}^{*}(A,\mathbf{x}_{0})$.

Let us define $\varphi(\mathbf{x}):\mathbb{R}^{d}\to\mathbb{R}^{+}$ be
$\varphi(\mathbf{x}):=\max_{i}|Qx_{i}|$. Using this notation,
$w_{0}^{*}(A,\mathbf{x}_{0})=\varphi(Q^{-1}\mathbf{x}_{0})$, and

$\begin{gathered}w_{0}^{*}(\tilde{A},\mathbf{x}_{0}):=w_{0}^{*}(TAT^{\prime},\mathbf{x}_{0})=\varphi(Q^{-1}T^{\prime}\mathbf{x}_{0})=w_{0}^{*}(A,T^{\prime}\mathbf{x}_{0}).\\\
w_{0}^{*}(\tilde{A},\mathbf{x}_{0})\stackrel{{\scriptstyle
d}}{{=}}w_{0}^{*}(A,\mathbf{x}_{0})\Longrightarrow
w_{0}^{*}(A,T^{\prime}\mathbf{x}_{0})\stackrel{{\scriptstyle
d}}{{=}}w_{0}^{*}(A,\mathbf{x}_{0}).\end{gathered}$ (S.6)

The above equation implies that the distribution of
$w_{0}^{*}(A,\mathbf{x}_{0})$ is statistically invariant under an arbitrary
orthogonal transformation of $\mathbf{x}_{0}$, therefore its distribution
depends only on $|\mathbf{x}_{0}|$. Note that the orbit of a fixed
$\mathbf{x}_{0}$ under all orthogonal transformations is
$|\mathbf{x}_{0}|\cdot S^{d-1}$, a sphere in $\mathbb{R}^{d}$ with radius
$|\mathbf{x}_{0}|$. Furthermore, if $T\sim\mathrm{Haar}(O(d))$, the
distribution of $T^{\prime}\mathbf{x}_{0}$ must be the uniform distribution on
$|\mathbf{x}_{0}|\cdot S^{d-1}$. In particular, if we let $T=Q^{\prime}$, we
have

$w_{0}^{*}(A,\mathbf{x}_{0})\stackrel{{\scriptstyle
d}}{{=}}w_{0}^{*}(\Lambda,\tilde{\mathbf{x}}_{0}),\qquad\tilde{\mathbf{x}}_{0}:=Q\mathbf{x}_{0}\sim\mathrm{Unif}\left(|\mathbf{x}_{0}|\cdot
S^{d-1}\right).$ (S.7)

Conditionally on $|\mathbf{x}_{0}|=r$, $r\in\mathbb{R}^{+}$ and using
Corollary S2, we know that

$\begin{split}E\left(w_{0}^{*}(A,\mathbf{x}_{0})^{2}\big{|}|\mathbf{x}_{0}|=r\right)&=E\left(w_{0}^{*}(\Lambda,\tilde{\mathbf{x}}_{0})^{2}\big{|}|\mathbf{x}_{0}|=r\right)\\\
&=ES_{\min}(T,r)\approx C_{\mathrm{Weibull}(1,1/2)}\cdot\frac{\pi
r^{2}}{2d^{3}}.\end{split}$ (S.8)

Here $C_{\mathrm{Weibull}(1,1/2)}$ is the second order moment of
$\mathrm{Weibull}(1,1/2)$.

On the other hand, Assumption (b) states that as a function of $d$,
$E|\mathbf{x}_{0}|^{2}$ grows at rate $o(d^{3})$. Therefore

$\begin{split}E\left(w_{0}^{*}(A,\mathbf{x}_{0})^{2}\right)&=E_{r}\left(E\left(w_{0}^{*}(A,\mathbf{x}_{0})^{2}\big{|}|\mathbf{x}_{0}|=r\right)\right)\\\
&\approx C_{\mathrm{Weibull}(1,1/2)}\cdot\frac{\pi
E|\mathbf{x}_{0}|^{2}}{2d^{3}}=o(1)\to 0.\end{split}$

As a special case, if the second order moments of $x_{0,i}$ are bounded by a
constant (as mentioned in Remark 3), $E|\mathbf{x}_{0}|^{2}$ will grow at a
speed of $O(d^{1})$, therefore
$E\left(w_{0}^{*}(A,\mathbf{x}_{0})^{2}\right)=O(d^{-2})$.

∎

Below we propose a conjecture on asymmetric matrices.

######  Conjecture S2.1.

We conjecture that Equation (23) is true for system matrix $A$ sampled from
the Ginibre ensemble and initial condition $\mathbf{x}_{0}$ such that the
second order moments of $x_{0,i}$ are bounded by a constant.

Although we were not able to formally prove Conjecture S2.1 due to technical
difficulties ($Q$ is no longer an orthogonal matrix for an asymmetric matrix
sampled from GinOE), this conjecture was numerically verified by us in
numerical experiments; for example, the mini-simulation described in Section
3.2.

##  S3 Some Expected Values Related to Random Matrices

Below we derive some expected values related to random matrices that are
useful in deriving the PIS.

######  Proposition S1.

Let $\bm{\epsilon}\in M_{d\times n}$ be a random matrix with $i.i.d.$ entries
$\epsilon_{ij}\sim N(0,\sigma^{2})$. Let $A\in M_{n\times d}$ and $B\in
M_{n\times n}$ be two deterministic matrices. We have

$E\left(\bm{\epsilon}A\bm{\epsilon}\right)=\sigma^{2}A^{\prime},\qquad
E\left(\bm{\epsilon}B\bm{\epsilon}^{\prime}\right)=\sigma^{2}\mathrm{tr}(B)\cdot
I_{d}.$ (S.9)

###### Proof.

Let $U=E\left(\bm{\epsilon}A\bm{\epsilon}\right)$ and
$V=E\left(\bm{\epsilon}B\bm{\epsilon}^{\prime}\right)$, we know that

$\begin{split}U_{ij}&=\sum_{l=1}^{n}\sum_{k=1}^{d}A_{lk}E(\epsilon_{il}\epsilon_{kj})=A_{ji}E\epsilon_{ij}^{2}=\sigma^{2}A_{ji}.\\\
V_{ij}&=\sum_{l=1}^{n}\sum_{k=1}^{n}B_{lk}E(\epsilon_{il}\epsilon_{jk})=\begin{cases}0,&i\neq
j,\\\
\sum_{l=1}^{n}B_{ll}E\epsilon_{il}^{2}=\sigma^{2}\mathrm{tr}(B),&i=j.\end{cases}\end{split}$

∎

##  S4 Proof of Main Theorems

Before we prove Theorem 3.2, we need to develop the following lemma that
establishes the relationship between the linear independence of a set of
functions and the invertibility of their pairwise inner product matrix.

######  Lemma S1.

Let $\mathbf{x}(t)\in\mathbb{R}^{d}$ be a set of continuous functions defined
on $[0,T]$. $\mathbf{x}(t)$ is linearly independent iff its pairwise inner
product matrix $\Sigma_{\mathbf{x}\mathbf{x}}:=\Sigma_{\mathbf{x}\mathbf{x}}$
is invertible.

###### Proof.

  * •

The “$\Longleftarrow$” part. Let us assume that member functions in
$\mathbf{x}(t)$ are linearly dependent. There exist nonzero vector
$\mathbf{c}$ such that $\mathbf{c}^{\prime}\mathbf{x}(t)\equiv 0$ for all
$t\in[0,T]$. Due to the bilinearity of the pairwise inner product matrix, we
have

$\mathbf{c}^{\prime}\Sigma_{\mathbf{x}\mathbf{x}}\mathbf{c}=\langle\mathbf{c}^{\prime}\mathbf{x}(t),\;\mathbf{c}^{\prime}\mathbf{x}(t)\rangle=\langle
0,\;0\rangle=0.$

  * •

The “$\Longrightarrow$” part. As we mentioned earlier,
$\Sigma_{\mathbf{x}\mathbf{x}}$ is singular if and only if there exists
nonzero $\mathbf{v}_{0}\in\mathbb{R}^{d}$ such that

$\mathbf{v}_{0}^{\prime}\Sigma_{\mathbf{x}\mathbf{x}}\mathbf{v}_{0}=0.$

Due to bilinearity, we know
$\langle\mathbf{v}_{0}^{\prime}\mathbf{x}(t),\;\mathbf{v}_{0}^{\prime}\mathbf{x}(t)\rangle=\|\mathbf{v}_{0}^{\prime}\mathbf{x}(t)\|_{L^{2}}^{2}=0$.
This implies that $\mathbf{x}(t)$ is essentially zero on $[0,T]$. Together
with the continuity assumption of $\mathbf{x}(t)$, we know that
$\mathbf{x}(t)$ must be zero for all $t\in[0,T]$.

∎

######  Corollary S2.

Assume that $\bm{\xi}(t):=(\xi_{1}(t),\xi_{2}(t),\dots,\xi_{d}(t))^{\prime}$
is a set of linearly independent functions, $C\in M_{d\times d}$ is a matrix.
Set of functions $\mathbf{x}(t):=C\bm{\xi}(t)$ is linearly independent iff
matrix $C$ is invertible.

###### Proof.

Based on the bilinearity, we know that

$\Sigma_{\mathbf{x}\mathbf{x}}=C\Sigma_{\bm{\xi}}C^{\prime}.$

Lemma S1 says that $\Sigma_{\bm{\xi}}$ is invertible. So
$\Sigma_{\mathbf{x}\mathbf{x}}$ is invertible (hence linearly independent) iff
$C$ is invertible. ∎

The next lemma shows that when all eigenvalues are distinct (which is our
assumption), the basic “building blocks” of $\mathbf{x}(t)$ are linearly
independent. To this end, we first define
$\mathbf{z}(t)=(z_{1}(t),z_{2}(t),\dots,z_{d}(t))^{\prime}$ to be a vector of
_fundamental solutions_ of Equation (3) in which $A$ has Jordan canonical form
specified in Equation (6). Elements in $\mathbf{z}(t)$ can be represented as
follows

$z_{i}(t)=\begin{cases}e^{c_{i}t},&k=1,2,\dots,K_{1},\\\ e^{a_{k}}\cos
b_{k}t,&k=K_{1}+\frac{i-K_{1}+1}{2},\;i=K_{1}+1,K_{1}+3,\dots,d-1.\\\
e^{a_{k}}\sin
b_{k}t,&k=K_{1}+\frac{i-K_{1}}{2},\;i=K_{1}+2,K_{1}+4,\dots,d.\end{cases}$
(S.10)

As a reminder, we point out that the relationship between $i$ and $k$ is
provided in Equation (8). In particular, the $k$th pair of complex eigenvalues
$a_{k}\pm b_{k}i$, for $k=K_{1}+1,\dots,K$, corresponds with a pair of
fundamental solutions $(z_{i(k)}(t),z_{i(k)+1}(t))$, where $i(k):=2k-K_{1}-1$.

######  Lemma S3.

As a set of functions, $\mathbf{z}(t)$ is linearly independent iff all the
eigenvalues of $A$ are distinct.

###### Proof.

It is easy to see that for real eigenvalues ($k=1,2,\dots,K_{1}$), we have

$D^{m}z_{k}(t)=c_{k}^{m}z_{k}(t).$

For complex eigenvalues ($k=K_{1}+1,\dots,K$), let

$\begin{gathered}r_{k}:=|a_{k}+b_{k}i|=\sqrt{a_{k}^{2}+b_{k}^{2}},\qquad\theta_{k}:=\arg(a_{k}+b_{k}i).\\\
a_{k,m}:=r_{k}^{m}\cos(m\theta_{k}),\qquad
b_{k,m}:=r_{k}^{m}\sin(m\theta_{k}).\end{gathered}$

We have

$\begin{split}D^{m}\begin{pmatrix}z_{i(k)}(t)\\\
z_{i(k)+1}(t)\end{pmatrix}&=\begin{pmatrix}a_{k}&-b_{k}\\\
b_{k}&a_{k}\end{pmatrix}^{m}\begin{pmatrix}z_{i(k)}(t)\\\
z_{i(k)+1}(t)\end{pmatrix}\\\ &=\begin{pmatrix}a_{k,m}&-b_{k,m}\\\
b_{k,m}&a_{k,m}\end{pmatrix}\begin{pmatrix}z_{i(k)}(t)\\\
z_{i(k)+1}(t)\end{pmatrix}\\\
&=\begin{pmatrix}a_{k,m}z_{i(k)}(t)-b_{k,m}z_{i(k)+1}(t)\\\
b_{k,m}z_{i(k)}(t)+a_{k,m}z_{i(k)+1}(t)\end{pmatrix}\\\
&=\begin{pmatrix}z_{i(k)}(t)&-z_{i(k)+1}(t)\\\
z_{i(k)+1}(t)&z_{i(k)}(t)\end{pmatrix}\begin{pmatrix}a_{k,m}\\\
b_{k,m}\end{pmatrix}.\end{split}$

Therefore, we can represent the $m$th derivative of $\mathbf{z}(t)$ as
$D^{m}\mathbf{z}(t)=Z(t)C_{m}$, where

$Z(t):=\begin{pmatrix}z_{1}(t)&&&&&\\\ &\ddots&&&&\\\ &&z_{K_{1}}(t)&&&\\\
&&&\begin{matrix}z_{K_{1}+1}(t)&-z_{K_{1}+2}(t)\\\
z_{K_{1}+2}(t)&z_{K_{1}+1}(t)\end{matrix}&&\\\ &&&&\ddots&\\\
&&&&&\begin{matrix}z_{d-1}(t)&-z_{d}(t)\\\
z_{d}(t)&z_{d-1}(t)\end{matrix}\end{pmatrix},$

and $C_{m}$ is a vector with the following elements

$\quad C_{m}:=\begin{pmatrix}c_{1}^{m}\\\ \vdots\\\ c_{K_{1}}^{m}\\\ \hline\cr
a_{K_{1}+1,m}\\\ b_{K_{1}+1,m}\\\ \hline\cr\vdots\\\ \hline\cr a_{K,m}\\\
b_{K,m}\\\ \end{pmatrix}.$

As a special case, it is easy to see that $C_{0}=1_{d}$. Let us define

$C:=\begin{pmatrix}C_{0}&C_{1}&\dots&C_{d}\end{pmatrix}.$

The Wronskian of $\mathbf{z}(t)$ can be represented as

$\begin{split}W(t)&=\left|Z(t)C\right|=\left|Z(t)\right|\cdot\left|C\right|.\end{split}$
(S.11)

It is easy to show that

$\begin{split}\left|Z(t)\right|&=\prod_{k=1}^{K_{1}}z_{k}(t)\cdot\prod_{k=K_{1}+1}^{K}\left(z_{k}^{2}+z_{k+1}^{2}\right)\\\
&=\exp\left(\Big{(}\sum_{k=1}^{K_{1}}c_{k}+\sum_{k=K_{1}+1}^{K}a_{k}\Big{)}t\right)\neq
0.\end{split}$ (S.12)

Therefore, $W(t)\neq 0$, hence elements in $\mathbf{z}(t)$ are linearly
independent if and only if $|C|\neq 0$.

The determinant of $C$ can be computed by the following technique. First, we
notice that

$\begin{gathered}\begin{pmatrix}1&i\\\
1&-i\end{pmatrix}\begin{pmatrix}a_{k,m}\\\
b_{k,m}\end{pmatrix}=\begin{pmatrix}a_{k,m}+b_{k,m}i\\\
a_{k,m}-b_{k,m}i\end{pmatrix}=\begin{pmatrix}\lambda_{i(k)}^{m}\\\
\lambda_{i(k)+1}^{m}\end{pmatrix},\qquad k=K_{1}+1,\dots,K.\\\ \end{gathered}$

Here $(\lambda_{i(k)},\lambda_{i(k)+1})$ are the pair of complex eigenvalues
associated with the $k$th block. Because $\left(\begin{smallmatrix}1&i\\\
1&-i\end{smallmatrix}\right)^{-1}=\left(\begin{smallmatrix}\frac{1}{2}&\frac{1}{2}\\\
\frac{-i}{2}&\frac{i}{2}\end{smallmatrix}\right)$, we can rewrite the $C$
matrix as follows

$\begin{split}C&=\begin{pmatrix}I_{K_{1}}&&&\\\
&\begin{matrix}\frac{1}{2}&\frac{1}{2}\\\
\frac{-i}{2}&\frac{i}{2}\end{matrix}&&\\\ &&\ddots&\\\
&&&\begin{matrix}\frac{1}{2}&\frac{1}{2}\\\
\frac{-i}{2}&\frac{i}{2}\end{matrix}\\\
\end{pmatrix}\begin{pmatrix}1&\lambda_{1}&\lambda_{1}^{2}&\dots&\lambda_{1}^{d-1}\\\
1&\lambda_{2}&\lambda_{2}^{2}&\dots&\lambda_{2}^{d-1}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
1&\lambda_{d}&\lambda_{d}^{2}&\dots&\lambda_{d}^{d-1}\\\
\end{pmatrix}\end{split}$ (S.13)

Therefore

$\left|C\right|=\left(\frac{i}{2}\right)^{K_{2}}\prod_{1\leqslant i<j\leqslant
d}(\lambda_{j}-\lambda_{i}),\qquad\left|C\right|\neq 0\text{ iff all
$\lambda_{i}$ are distinct.}$ (S.14)

∎

###### Proof of Theorem 3.2.

Again, I first prove the special case in which all eigenvalues are real. In
this case, we know that

$\begin{split}\mathbf{x}(t)&=e^{At}\mathbf{x}_{0}=Q\cdot\mathrm{diag}(e^{c_{i}t})\cdot
Q^{-1}\mathbf{x}_{0}\\\ &=Q\begin{pmatrix}w_{0,1}e^{c_{1}t}\\\
w_{0,2}e^{c_{2}t}\\\ \vdots\\\ w_{0,d}e^{c_{d}t}\\\
\end{pmatrix}=Q\cdot\mathrm{diag}(w_{0,i})\begin{pmatrix}e^{c_{1}t}\\\
\vdots\\\ e^{c_{d}t}\end{pmatrix}.\end{split}$ (S.15)

By Lemma S3, $(e^{c_{1}t},\dots,e^{c_{d}t})^{\prime}$ is linearly independent,
so $\mathbf{x}(t)$ is linearly independent iff $Q\cdot\mathrm{diag}(w_{0,i})$
is invertible. Now $Q$ is invertible based on our assumption, according to
Corollary S2, whether or not $\mathbf{x}(t)$ is linearly independent only
depends on the invertibility of $\mathrm{diag}(w_{0,i})$, which in turn is
equivalent to $w_{0,i}\neq 0$ for $i=1,2,\dots,d$. Using Theorem 2.4, we know
that $A$ is identifiable at $\mathbf{x}_{0}$ iff all $|w_{0,k}|\neq 0$, which
is equivalent to the invertibility of matrices $\mathrm{diag}(w_{0,i})$ and
$Q\cdot\mathrm{diag}(w_{0,i})$ (because $Q$ is invertible), which in turn is
equivalent to the linear independence of
$(e^{c_{1}t},\dots,e^{c_{d}t})^{\prime}$ and the invertibility of
$\Sigma_{\mathbf{x}\mathbf{x}}$. ∎

###### Proof of Theorem 3.3.

The “$\Longrightarrow$” direction. By assumption,
$D\mathbf{x}(t)=A\mathbf{x}(t)$, we have

$\Sigma_{D\mathbf{x},\mathbf{x}}=\langle
A\mathbf{x}(t),\;\mathbf{x}(t)\rangle=A\Sigma_{\mathbf{x}\mathbf{x}}.$

Because $\Sigma_{\mathbf{x}\mathbf{x}}$ is invertible, we can multiply both
sides of Equation (30) by $\Sigma_{\mathbf{x}\mathbf{x}}^{-1}$, therefore $A$
must be $\Sigma_{D\mathbf{x},\mathbf{x}}\Sigma_{\mathbf{x}\mathbf{x}}^{-1}$.

The “$\Longleftarrow$” direction. Because $\Sigma_{\mathbf{x}\mathbf{x}}$ is
invertible, based on Theorem 3.2, $A$ is identifiable at $\mathbf{x}_{0}$,
which means that there is only one $A$ that can produce the solution curve
$\mathbf{x}(t)$, therefore
$\Sigma_{D\mathbf{x},\mathbf{x}}\Sigma_{\mathbf{x}\mathbf{x}}^{-1}$ must be
the system matrix. ∎

###### Proof of Theorem 3.4.

Due to Theorem 3.3, we know that
$A=\Sigma_{D\mathbf{x},\mathbf{x}}\cdot\Sigma_{\mathbf{x}\mathbf{x}}^{-1}$.
Therefore, it suffices to show that $\hat{A}$ is a continuous matrix-valued
functional of $\hat{\mathbf{x}}(t)$ and $D\hat{\mathbf{x}}(t)$. Clearly, we
have

$\begin{split}\langle\hat{\mathbf{x}}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}\rangle-\langle\mathbf{x}_{i},\;\mathbf{x}_{i^{\prime}}\rangle&=\langle\hat{\mathbf{x}}_{i}-\mathbf{x}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}\rangle+\langle\mathbf{x}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}-\mathbf{x}_{i^{\prime}}\rangle,\\\
\left(\langle\hat{\mathbf{x}}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}\rangle-\langle\mathbf{x}_{i},\;\mathbf{x}_{i^{\prime}}\rangle\right)^{2}&\leqslant(\langle\hat{\mathbf{x}}_{i}-\mathbf{x}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}\rangle)^{2}+(\langle\mathbf{x}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}-\mathbf{x}_{i^{\prime}}\rangle)^{2}\\\
&\quad+2|\langle\hat{\mathbf{x}}_{i}-\mathbf{x}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}\rangle|\cdot|\langle\mathbf{x}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}-\mathbf{x}_{i^{\prime}}\rangle|\\\
&\leqslant\delta_{1}^{2}\left(\|\mathbf{x}_{i}\|+\|\mathbf{x}_{i^{\prime}}\|\right)^{2}.\\\
\|\Sigma_{\hat{\mathbf{x}}\hat{\mathbf{x}}}-\Sigma_{\mathbf{x}\mathbf{x}}\|_{F}^{2}&=\sum_{i,i^{\prime}=1}^{d}\left(\langle\hat{\mathbf{x}}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}\rangle-\langle\mathbf{x}_{i},\;\mathbf{x}_{i^{\prime}}\rangle\right)^{2}dt\\\
&\leqslant 2\delta_{1}^{2}\|\mathbf{x}\|^{2}.\end{split}$ (S.16)

$\begin{split}\langle D\hat{\mathbf{x}},\;\hat{\mathbf{x}}\rangle-\langle
D\hat{\mathbf{x}},\;\hat{\mathbf{x}}\rangle&=\langle
D\hat{\mathbf{x}}_{i}-D\mathbf{x}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}\rangle+\langle
D\mathbf{x}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}-\mathbf{x}_{i^{\prime}}\rangle,\\\
\left(\langle D\hat{\mathbf{x}},\;\hat{\mathbf{x}}\rangle-\langle
D\hat{\mathbf{x}},\;\hat{\mathbf{x}}\rangle\right)^{2}&\leqslant\left(\delta_{2}\|\hat{\mathbf{x}}_{i^{\prime}}\|+\delta_{1}\|D\mathbf{x}_{i}\|\right)^{2}\\\
&\leqslant\left(\delta_{2}\|\mathbf{x}_{i^{\prime}}\|+\delta_{1}+\delta_{1}\|D\mathbf{x}_{i}\|\right)\leqslant\delta\left(\|\mathbf{x}_{i^{\prime}}\|+\|D\mathbf{x}_{i}\|+1\right).\\\
\|\Sigma_{D\hat{\mathbf{x}},\hat{\mathbf{x}}}-\Sigma_{D\mathbf{x},\mathbf{x}}\|_{F}^{2}&=\sum_{i,i^{\prime}=1}^{d}\left(\langle
D\hat{\mathbf{x}}_{i},\;\hat{\mathbf{x}}_{i^{\prime}}\rangle-\langle
D\mathbf{x}_{i},\;\mathbf{x}_{i^{\prime}}\rangle\right)^{2}dt\\\
&\leqslant\delta^{2}\left(\|\mathbf{x}\|^{2}+\|D\mathbf{x}\|^{2}+1\right).\end{split}$
(S.17)

In other words, as matrix-valued functionals of $D\mathbf{x}$ and
$\mathbf{x}$, $\Sigma_{\hat{\mathbf{x}},\hat{\mathbf{x}}}$ and
$\Sigma_{D\hat{\mathbf{x}},\hat{\mathbf{x}}}$ are continuous w.r.t. the
$L^{2}$ metric. Based on Theorem 3.2, the assumption that the ODE system is
identifiable at $\mathbf{x}_{0}$ implies that $\Sigma_{\mathbf{x}\mathbf{x}}$
is invertible. Consequently,
$\hat{A}:=\hat{\Sigma}_{D\mathbf{x},\mathbf{x}}\cdot\hat{\Sigma}_{\mathbf{x}\mathbf{x}}^{-1}$
must be a matrix-valued continuous functional of $D\mathbf{x}$ and
$\mathbf{x}$. ∎

##  S5 Identifiability Issues Induced by Repeated Eigenvalues

In this subsection, we argue that when there exist repeated eigenvalues (real
or complex), $A$ will not be identifiable for any $\mathbf{x}_{0}$. We then
provide a representation of $[A]_{\mathbf{x}_{0}}$ in this case. Although it
is well known that the set of matrices with repeated eigenvalues has Lebesgue
measure zero in $M_{d\times d}$, In reality, $A$ may still have eigenvalues
with similar numerical values due to randomness, hence work presented in this
section is relevant for real world applications with uncertainty.

First, we would like to provide an interpretation based on the Jordan
canonical decomposition. Suppose there exists a real eigenvalue of $A$ with
multiplicity greater than one. WLOG, we can always rearrange the JCP so that
these repeated eigenvalues are ordered as the first $m_{1}$ eigenvalues,
$\lambda_{1}=\dots=\lambda_{m_{1}}$. The combined Jordan block for these
repeated eigenvalues is

$\Lambda_{1}=\begin{pmatrix}\lambda_{1}&&\\\ &\ddots&\\\
&&\lambda_{1}\end{pmatrix}=I_{m_{1}}\otimes J_{1},\qquad
J_{1}:=(\lambda_{1})_{1\times 1}.$

Similarly, if $A$ has $m_{1}$ pairs of repeated complex eigenvalues, we can
arrange the JCP so that they are the first $2m_{1}$ eigenvalues (denoted as
$a_{1}\pm b_{1}i$) in the semi-diagonal matrix $\Lambda$, thus the combined
Jordan block of them is the following $2m_{k}\times 2m_{k}$-dimensional semi-
diagonal matrix

$\Lambda_{1}=I_{m_{1}}\otimes J_{1},\qquad
J_{1}:=\begin{pmatrix}a_{1}&-b_{1}\\\ b_{1}&a_{1}\end{pmatrix}.$

In either case, the combined Jordan block of the repeated eigenvalue,
$\Lambda_{1}=I_{m_{1}}\otimes J_{1}$, is the top-left block of the middle
matrix in the JCP of $A$:

$\begin{gathered}A=\begin{pmatrix}Q_{1}&Q_{-1}\end{pmatrix}\begin{pmatrix}\Lambda_{1}&\\\
&\Lambda_{-1}\end{pmatrix}\begin{pmatrix}R_{1}\\\
R_{-1}\end{pmatrix}=Q_{1}\Lambda_{1}R_{1}+Q_{-1}\Lambda_{-1}R_{-1}.\end{gathered}$
(S.18)

Here $\Lambda_{-1}$ represents eigenvalues not in $\Lambda_{1}$, $Q_{1}$ is
the collection of the first $m_{1}$ (or $2m_{1}$, if the repeated eigenvalues
are complex) columns of $Q$, which is a basis for the invariant subspace
associated with $\Lambda_{1}$. Likewise, $Q_{-1}$ is the collection of
remaining column vectors of $Q$ and a basis of the invariant subspace
associated with $\Lambda_{-1}$. $R_{1}$ and $R_{-1}$ are the corresponding
left- and right-submatrices of $Q^{-1}$.

When considered as a $m_{1}\times m_{1}$-dimensional (if the repeated
eigenvalue is real) or $2m_{1}\times 2m_{1}$-dimensional (if the repeated
eigenvalue is complex) _subsystem_ , $\Lambda_{1}$ is not identifiable at any
initial point $\mathbf{v}\in\mathbb{R}^{m_{1}}$ and the solution curves
satisfies the following equation

$e^{tA}\mathbf{v}=e^{t(A+UDU^{\prime})}\mathbf{v}.$ (S.19)

Here $D$ is an arbitrary matrix and $U$ is a semi-orthogonal matrix that
depends on the initial condition $\mathbf{v}$. Technical details of these
results are summarized in Lemmas S1 and S2 below.

######  Lemma S1.

Let $\lambda_{1}\in\mathbb{R}$ and $A=\lambda_{1}I_{m_{1}}$,
$\mathbf{v}\in\mathbb{R}^{m_{1}}$. We have

$e^{tA}\mathbf{v}=e^{t(A+UDU^{\prime})}\mathbf{v}.$

Here $U\in M_{m_{1}\times(m_{1}-1)}$ is a semi-orthogonal matrix such that
$U^{\prime}\mathbf{v}=\mathbf{0}$, and $D$ is an arbitrary
$(m_{1}-1)\times(m_{1}-1)$-dimensional matrix.

###### Proof.

WLOG, we may assume that $\|\mathbf{v}\|=1$. By construction, matrix
$W:=\begin{pmatrix}\mathbf{v}&U\end{pmatrix}$ is orthogonal. It is easy to see
that

$\begin{split}A+UDU^{\prime}&=W\left(\lambda_{1}I_{m_{1}}\right)W^{\prime}+\begin{pmatrix}\mathbf{v}&U\end{pmatrix}\begin{pmatrix}0&\mathbf{0}\\\
\mathbf{0}&D\end{pmatrix}\begin{pmatrix}\mathbf{v}^{\prime}\\\
U^{\prime}\end{pmatrix}\\\
&=W\tilde{D}W^{\prime},\qquad\tilde{D}:=\begin{pmatrix}\lambda_{1}&\\\
&\lambda_{1}I_{m_{1}-1}+D\end{pmatrix}.\end{split}$

$\begin{gathered}e^{tA}=e^{t\lambda_{1}}I_{m_{1}},\qquad
e^{tA}\mathbf{v}=e^{t\lambda_{1}}\cdot\mathbf{v}.\\\
\begin{split}e^{t(A+UDU^{\prime})}\mathbf{v}&=We^{t\tilde{D}}W^{\prime}\mathbf{v}=\begin{pmatrix}\mathbf{v}&U\end{pmatrix}\begin{pmatrix}e^{t\lambda_{1}}&\\\
&e^{t(\lambda_{1}I_{m_{1}-1}+D)}\end{pmatrix}\begin{pmatrix}1\\\
\mathbf{0}\end{pmatrix}\\\
&=e^{t\lambda_{1}}\cdot\mathbf{v}=e^{tA}\mathbf{v}.\end{split}\end{gathered}$

∎

Likewise, for repeated complex eigenvalues, we have the following lemma.

######  Lemma S2.

Let $a\pm bi\in\mathbb{C}$ and
$A=I_{m_{1}}\otimes\left(\begin{smallmatrix}a&-b\\\
b&a\end{smallmatrix}\right)$. For $\mathbf{v}\in\mathbb{R}^{2m_{1}}$, we have

$e^{tA}\mathbf{v}=e^{t(A+UDU^{\prime})}\mathbf{v}.$

Here $D\in M_{(2m_{1}-2)\times(2m_{1}-2)}$ is an arbitrary matrix and $U\in
M_{2m_{1}\times(2m_{1}-2)}$ is a semi-orthogonal matrix such that

$U^{\prime}\begin{pmatrix}\mathbf{v}&\bar{\mathbf{v}}\end{pmatrix}=\mathbf{0}_{2m_{1}},\qquad\bar{\mathbf{v}}:=\left(I_{m_{1}}\otimes\left(\begin{smallmatrix}0&-1\\\
1&0\end{smallmatrix}\right)\right)\mathbf{v}.$ (S.20)

###### Proof.

WLOG, assume that $\|\mathbf{v}\|=1$. Let
$W:=\begin{pmatrix}\mathbf{v}&\bar{\mathbf{v}}&U\end{pmatrix}$, it is easy to
check that $W$ is an orthogonal matrix. For simplicity, let us denote matrix
$\left(\begin{smallmatrix}a&-b\\\ b&a\end{smallmatrix}\right)$ by
$\Lambda_{2}$ and matrix
$\begin{pmatrix}\mathbf{v}&\bar{\mathbf{v}}\end{pmatrix}$ by $V_{2}$. Based on
matrix analysis, we can derive

$\begin{gathered}\left(I_{m_{1}}\otimes\Lambda_{2}\right)V_{2}=V_{2}\Lambda_{2},\\\
W^{\prime}\left(I_{m_{1}}\otimes\Lambda_{2}\right)W=I_{m_{1}}\otimes\Lambda_{2}.\end{gathered}$
(S.21)

Therefore

$\begin{split}A+UDU^{\prime}&=W\left(I_{m_{1}}\otimes\Lambda_{2}\right)W^{\prime}+\begin{pmatrix}V_{2}&U\end{pmatrix}\begin{pmatrix}\mathbf{0}_{2\times
2}&\mathbf{0}_{2\times(2m_{1}-2)}\\\ \mathbf{0}_{(2m_{1}-2)\times
2}&D\end{pmatrix}\begin{pmatrix}V_{2}^{\prime}\\\ U^{\prime}\end{pmatrix}\\\
&=W\tilde{D}W^{\prime},\qquad\tilde{D}:=\begin{pmatrix}\Lambda_{1}&\\\
&I_{m_{1}-1}\otimes\Lambda_{2}+D\end{pmatrix}.\end{split}$

$\begin{gathered}e^{tA}=I_{m_{1}}\otimes e^{t\Lambda_{1}},\qquad
e^{tA}V_{2}=V_{2}\cdot e^{t\Lambda_{1}}.\\\
\begin{split}e^{t(A+UDU^{\prime})}V_{2}&=We^{t\tilde{D}}W^{\prime}V_{2}=\begin{pmatrix}V_{2}&U\end{pmatrix}\begin{pmatrix}e^{t\Lambda_{1}}&\\\
&e^{t(I_{m_{1}-1}\otimes\Lambda_{2}+D)}\end{pmatrix}\begin{pmatrix}I_{2}\\\
\mathbf{0}\end{pmatrix}\\\ &=e^{t\Lambda_{1}}\cdot
V_{2}=e^{tA}V_{2}.\end{split}\end{gathered}$

The last equality implies that
$e^{t(A+UDU^{\prime})}\mathbf{v}=e^{tA}\mathbf{v}$ because $\mathbf{v}$ is the
first column of $V_{2}$. ∎

As a direct consequence of Lemmas S1 and S2, we have the following theorem.

######  Theorem S3.

Let $A$ be a matrix with repeated real or complex eigenvalues but no Jordan
block with nilpotent components. WLOG, we arrange the JCF such that $\lambda$,
the repeated real eigenvalue or pair of complex eigenvalues, corresponds with
the first Jordan block:

$P\begin{pmatrix}J_{1}&\\\ &J_{2}\end{pmatrix}P^{-1},\qquad
P=\begin{pmatrix}P_{1}&P_{2}\end{pmatrix}.$

Here $P_{1}$ is the set of (generalized)-eigenvalues associated with the
repeated eigenvalues, and $P_{2}$ are other generalized eigenvalues. Depending
on whether the repeated eigenvalue is real or complex, $J_{1}$ has the
following form ($m$ is the multiplicity of $\lambda$)

$J_{1}=\begin{cases}\lambda I_{m},&\lambda\in\mathbb{R},\\\
I_{m}\otimes\begin{pmatrix}a&-b\\\ b&a\end{pmatrix},&\lambda=a\pm
bi\in\mathbb{C},\;b\neq 0.\end{cases}$ (S.22)

We claim that $A$ cannot be identifiable for any
$\mathbf{x}_{0}\in\mathbb{R}^{d}$. Specifically, the
$(A,\mathbf{x}_{0})$-unidentifiable class has the following structure

$\forall B\in[A]_{\mathbf{x}_{0}},\qquad
B=P\begin{pmatrix}J_{1}+UDU^{\prime}&\\\ &J_{2}\end{pmatrix}P^{-1}.$ (S.23)

Here $D$ is an _arbitrary_ matrix in $M_{(m-1)\times(m-1)}$ (for real
$\lambda$) or $M_{(2m-2)\times(2m-2)}$ (for complex $\lambda$), $U$ is a semi-
orthogonal matrix such that for
$\mathbf{v}:=\begin{pmatrix}I&\mathbf{0}\end{pmatrix}P^{-1}\mathbf{x}_{0}$,

$\begin{gathered}\begin{cases}U\in M_{m\times(m-1)},\quad
U^{\prime}\mathbf{v}=\mathbf{0}_{m-1},&\lambda\in\mathbb{R}\\\ U\in
M_{m\times(2m-2)},\quad
U^{\prime}\begin{pmatrix}\mathbf{v}&\bar{\mathbf{v}}\end{pmatrix}=\mathbf{0},&\lambda=a\pm
bi\in\mathbb{C},\;b\neq 0.\end{cases}\end{gathered}$

######  Remarks 4.

The assumption that no Jordan block of $A$ contain a _nilpotent_ component in
Theorem S3 is to ensure that $A$ can be decomposed in the form of Equation
(6), even if there are repeated eigenvalues. Without this assumption, $A$ may
have Jordan blocks that look like

$J_{1}=\begin{pmatrix}c_{1}&1&0\\\ 0&c_{1}&1\\\ 0&0&c_{1}\end{pmatrix},\qquad
J_{2}=\begin{pmatrix}\begin{matrix}a_{2}&-b_{2}\\\
b_{2}&a_{2}\end{matrix}&\begin{matrix}1&0\\\ 0&1\end{matrix}\\\
\begin{matrix}0&0\\\ 0&0\end{matrix}&\begin{matrix}a_{2}&-b_{2}\\\
b_{2}&a_{2}\end{matrix}\end{pmatrix}.$

We believe it is theoretically possible to prove the equivalence of Lemmas S1
and S2, but it requires much more tedious computations that are not directly
related to the main focus of our manuscript, so we decide to focus on those
$A$ with two classical types of Jordan blocks in this study.

##  S6 Prior Information and Identifiability

As shown in Equation (17), we can write the explicit form of
$[A]_{\mathbf{x}_{0}}$ as an affine subspace in $M_{d\times d}$ based on
Theorem 2.5. This fact not only can serve as a guidance for us to use prior
information about $A$ to resolve the identifiability issue, but also implies
that we need to check the compatibility of prior information first, because
not all prior information on $A$ is compatible with $[A]_{\mathbf{x}_{0}}$. To
this end, we will use the following explicit definition of “prior information”
on $A$ throughout this study.

######  Definition S6.1 (Subset prior information on $A$).

A piece of subset prior information on $A$ is a subset $\mathcal{S}\subset
M_{d\times d}$ to which $A$ belongs. If $\mathcal{S}$ is a linear (affine)
subspace of $M_{d\times d}$, we call it linear (affine) prior information on
$A$.

As a remark, subset prior information on $A$ is “harder” than a typical
Bayesian prior information, in which we are given a prior distribution of $A$
so $A_{ij}$ can still take arbitrary numerical values.

In practice, most examples of subset prior information is affine prior
information in which a subset of $x_{ij}$ are assigned known values. If all
these given values equal zero, it is linear prior information; if some of
these values are nonzero, it is affine prior information.

######  Example S4.

Let us assume that the condition $A_{13}=0$, which is linear prior
information, is given to us in Example 1.

Due to the explicit form of $[A]_{\mathbf{x}_{0}}$ in Equation (17), we know
that $A_{13}=\frac{-3+b}{4}$, therefore $b=3$, and $\tilde{A}$ defined in
Equation (18) must be the _unique_ system matrix for the solution curve.

On the other hand, we cannot impose conditions $A_{11}=0$ and $A_{13}=0$
simultaneously, because $A_{11}=0$ implies $\frac{1+b}{4}=0$ or $b=-1$, which
contradicts with $b=3$ derived from $A_{13}=0$.

######  Definition S6.2.

A piece of subset prior information $\mathcal{S}$ is said to be compatible
with $(A,\mathbf{x}_{0})$ if $\mathcal{S}\cap[A]_{\mathbf{x}_{0}}$ is non-
empty. It is proper subset prior information for $(A,\mathbf{x}_{0})$ if
$\mathcal{S}\cap[A]_{\mathbf{x}_{0}}$ contains only one unique matrix.

If the prior information is _correct_ , that is, the true system matrix
$A\in\mathcal{S}$, $\mathcal{S}$ must be compatible with $(A,\mathbf{x}_{0})$
because $A\in[A]_{\mathbf{x}_{0}}$ so $\mathcal{S}\cap[A]_{\mathbf{x}_{0}}$ is
non-empty. However, in practice we most often only have an imperfect estimate
of $A$ and $\mathbf{x}_{0}$, so the prior information $\mathcal{S}$ _could_ be
incorrect, thus it is useful to check whether it is compatible with
$(\hat{A},\hat{\mathbf{x}}_{0})$ or not. It is also useful to see if
$\mathcal{S}$ is proper for $(A,\mathbf{x}_{0})$. These questions are answered
in part by the next Theorem.

######  Theorem S1.

Let $\mathcal{S}$ be affine prior information represented in this way:

$S\,\mathrm{vec}\left(A-A_{0}\right)=0_{L},\qquad S\in M_{L\times d^{2}}.$

$\mathcal{S}$ is proper prior information if and only if

$\begin{gathered}\mathrm{rank}\big{(}(\tilde{S}|\mathbf{b})\big{)}=\mathrm{rank}(\tilde{S}).\\\
\tilde{S}:=S\left((I_{0}Q^{-1})^{\prime}\otimes(QI_{0})\right),\qquad\mathbf{b}:=S\,\mathrm{vec}(A_{0}-A).\end{gathered}$
(S.24)

Here $\mathrm{rank}(\cdot)$ should be understood as the row rank.

###### Proof.

According to Equation (14), we have

$\begin{gathered}S\,\mathrm{vec}\left(A+QI_{0}DI_{0}Q^{-1}-A_{0}\right)=0\Longrightarrow\tilde{S}\,\mathrm{vec}(D)=\mathbf{b}.\\\
\end{gathered}$

A necessary and sufficient condition for equation
$\tilde{S}\,\mathrm{vec}(D)=\mathbf{b}$ have a unique solution is that the row
rank of matrix $(\tilde{S}|\mathbf{b})$ equals the row rank of
$\mathrm{rank}(\tilde{S})$. ∎

##  S7 Two Examples of Two-stage Methods

###  S7.1 Simple Two-stage Method

In this approach, no smoothing is required. We consider $Y$ as a good
approximation of $\mathbf{x}(t)$ evaluated at the time grid, and use simple
differences in the time direction to estimate $D\mathbf{x}(t)$. The two inner
products are estimated by

$\begin{gathered}\hat{\Sigma}_{\mathbf{x}\mathbf{x}}=YI_{n}Y^{\prime},\quad
K^{\text{simple}}=I_{n}.\\\
\hat{\Sigma}_{D\mathbf{x},\mathbf{x}}=YL^{\text{simple}}Y^{\prime},\\\
L^{\text{simple}}:=\frac{1}{\Delta t}\begin{pmatrix}-1&&&&\\\ 1&-1&&&\\\
&1&\ddots&&\\\ &&\ddots&&-1\\\
&&&&1\end{pmatrix}_{n\times(n-1)}\begin{pmatrix}I_{n-1}&0_{n-1}\end{pmatrix}\\\
=\frac{1}{\Delta t}\begin{pmatrix}[ccccc|c]-1&&&&&0\\\ 1&-1&&&&\vdots\\\
&1&\ddots&&&\\\ &&\ddots&&-1&\vdots\\\ &&&&1&0\end{pmatrix}_{n\times
n}.\end{gathered}$ (S.25)

###  S7.2 Functional Two-stage Method

Let $\bm{\phi}(t)=\left(\phi_{1}(t),\dots,\phi_{B}(t)\right)$ be a basis
system defined on $[0,T]$,
$J_{\phi}:=\langle\bm{\phi}(t),\;\bm{\phi}(t)\rangle$ be the inner product
matrix of those basis functions so that a candidate solution curve can be
represented as $C_{d\times B}\bm{\phi}(t)$, where $C$ is a matrix of linear
coefficients. Let $[D]_{B\times B}$ be the matrix representation of the
differential operator such that $D\bm{\phi}(t)=[D]\bm{\phi}(t)$. Let $\lambda$
be the roughness penalty parameter, and $H_{\lambda}\in M_{n\times B}$ be the
“hat matrix” that maps the discrete data to the linear coefficients of the
fitted curves, namely, $\hat{\mathbf{x}}(t)=\hat{C}\bm{\phi}(t)$,
$\hat{C}=YH_{\lambda}$.

In this case, we have

$\begin{gathered}\hat{\Sigma}_{\mathbf{x}\mathbf{x}}=\langle\hat{\mathbf{x}}(t),\;\hat{\mathbf{x}}(t)\rangle=\hat{C}J_{\phi}\hat{C}^{\prime}=YH_{\lambda}J_{\phi}H_{\lambda}^{\prime}Y^{\prime},\quad
K^{\text{smooth}}=H_{\lambda}J_{\phi}H_{\lambda}^{\prime}.\\\
\hat{\Sigma}_{D\mathbf{x},\mathbf{x}}=\langle\hat{C}[D]\bm{\phi}(t),\;\hat{C}\bm{\phi}(t)\rangle=YH_{\lambda}[D]J_{\phi}H_{\lambda}^{\prime}Y^{\prime},\quad
L^{\text{smooth}}=H_{\lambda}[D]J_{\phi}H_{\lambda}^{\prime}.\end{gathered}$
(S.26)

##  S8 Additional Examples

We provide a few additional examples in this section so that the readers can
have a better understanding of the $(A,x_{0})$-identifiability.

######  Example S5.

By definition, a one-dimensional invariant subspace is just the line generated
by an eigenvector of $A$.

######  Example S6.

Assume that

$\begin{gathered}A=Q\begin{pmatrix}J&\\\ &B\end{pmatrix}Q^{-1},\qquad
J=\begin{pmatrix}a&b\\\ -b&a\end{pmatrix},\qquad B\in M_{(d-2)\times(d-2)}.\\\
Q=\begin{pmatrix}[c|c|c]\mathbf{v}_{1}&\mathbf{v}_{2}&C\end{pmatrix},\quad\mathbf{v}_{1},\;\mathbf{v}_{2}\in\mathbb{R}^{d},\quad
C\in M_{d\times(d-2)}.\end{gathered}$

Then $\mathrm{span}(\mathbf{v}_{1},\mathbf{v}_{2})$ is a 2-dimensional
invariant subspace (associated with $J$). This is because for any
$w_{1},w_{2}\in\mathbb{R}$

$Q\begin{pmatrix}w_{1}\\\ w_{2}\\\
\mathbf{0}_{d-2}\end{pmatrix}=\begin{pmatrix}\mathbf{v}_{1}&\mathbf{v}_{2}&C\end{pmatrix}\begin{pmatrix}w_{1}\\\
w_{2}\\\
\mathbf{0}_{d-2}\end{pmatrix}=w_{1}\mathbf{v}_{1}+w_{2}\mathbf{v}_{2}.$
$\begin{split}A(w_{1}\mathbf{v}_{1}+w_{2}\mathbf{v}_{2})&=Q\begin{pmatrix}J&\\\
&B\end{pmatrix}Q^{-1}\cdot Q\begin{pmatrix}w_{1}\\\ w_{2}\\\
\hline\cr\mathbf{0}_{d-2}\end{pmatrix}\\\ &=Q\begin{pmatrix}aw_{1}+bw_{2}\\\
-bw_{1}+aw_{2}\\\ \hline\cr\mathbf{0}_{d-2}\end{pmatrix}\\\
&=(aw_{1}+bw_{2})\mathbf{v}_{1}+(-bw_{1}+aw_{2})\mathbf{v}_{2}\in\mathrm{span}(\mathbf{v}_{1},\mathbf{v}_{2}).\end{split}$

######  Example S7.

In this example, we will illustrate that $A$ with a repeated eigenvalue is not
identifiable for any $\mathbf{x}_{0}\in\mathbb{R}^{d}$, as declared by Theorem
S3.

Let $A=I_{2}$ and $\mathbf{x}_{0}=r(\cos\theta,\sin\theta)^{\prime}$ for
arbitrary $r\in\mathbb{R}^{+}$ and $\theta\in[0,2\pi)$. The trajectory is

$\mathbf{x}(t|I_{2},\mathbf{x}_{0})=re^{t}\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}.$

Now let

$B=\begin{pmatrix}1+\sin^{2}\theta&-\sin\theta\cos\theta\\\
-\sin\theta\cos\theta&1+\cos^{2}\theta\end{pmatrix}.$

It is easy to see that

$\begin{gathered}B=T_{B}\begin{pmatrix}1&0\\\
0&2\end{pmatrix}T_{B}^{\prime},\qquad
T_{B}=\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix}.\\\
\begin{split}\mathbf{x}(t|B,\mathbf{x}_{0})&=e^{tB}\mathbf{x}_{0}=T_{B}\begin{pmatrix}e^{t}&0\\\
0&e^{2t}\end{pmatrix}T_{B}^{\prime}\mathbf{x}_{0}\\\
&=r\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}e^{t}&0\\\
0&e^{2t}\end{pmatrix}\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\\\ &=r\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}e^{t}\\\
0\end{pmatrix}=re^{t}\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}=\mathbf{x}(t|I_{2},\mathbf{x}_{0}).\end{split}\end{gathered}$

Therefore, $A$ and $B$ must be in the same unidentifiable class.

In the next example, we show that certain network topology (the sparsity
structure of $A$) always imply unidentifiability.

######  Example S8.

Let $A$ be a matrix with two rows (or columns) with all zeros. Elementary
linear algebra shows that $\lambda=0$ must be an eigenvalue of $A$ with
multiplicity greater or equal to two. As a specific example, consider

$A=\begin{pmatrix}3&4&5\\\ 0&0&0\\\ 0&0&0\\\ \end{pmatrix},\quad
A=T_{A}\begin{pmatrix}3&0&0\\\ 0&0&0\\\ 0&0&0\\\ \end{pmatrix}T_{A}^{-1},\quad
T_{A}=\begin{pmatrix}1&-\frac{4}{5}&-\frac{5}{\sqrt{34}}\\\ 0&\frac{3}{5}&0\\\
0&0&\frac{3}{\sqrt{34}}\end{pmatrix}.$

It is easy to see that for an arbitrary initial condition $\mathbf{x}_{0}$
with the following decomposition

$\mathbf{x}_{0}=T_{A}\begin{pmatrix}c_{1}\\\ c_{2}\\\
c_{3}\end{pmatrix}=c_{1}\begin{pmatrix}1\\\ 0\\\
0\end{pmatrix}+c_{2}\begin{pmatrix}-\frac{4}{5}\\\ \frac{3}{5}\\\
0\end{pmatrix}+c_{3}\begin{pmatrix}-\frac{5}{\sqrt{34}}\\\ 0\\\
\frac{3}{\sqrt{34}}\end{pmatrix},$

the trajectory is

$\begin{split}\mathbf{x}(t|A,\mathbf{x}_{0})&=e^{tA}\mathbf{x}_{0}=\begin{pmatrix}c_{1}e^{3t}\\\
0\\\ 0\\\ \end{pmatrix}.\end{split}$

Such solution can be generated by the following alternative system

$B=T_{A}\begin{pmatrix}3&0&0\\\ 0&c_{3}^{2}&-c_{2}c_{3}\\\
0&-c_{2}c_{3}&c_{2}^{2}\\\ \end{pmatrix}T_{A}^{-1},$

because

$\begin{split}e^{tB}\mathbf{x}_{0}&=T_{A}\begin{pmatrix}e^{3t}&0_{2}\\\
0_{2}^{\prime}&e^{t\left(\begin{smallmatrix}-c_{3}\\\
c_{2}\end{smallmatrix}\right)\left(\begin{smallmatrix}-c_{3}&c_{2}\end{smallmatrix}\right)}\\\
\end{pmatrix}\begin{pmatrix}c_{1}\\\ \hline\cr c_{2}\\\ c_{3}\end{pmatrix}\\\
&=T_{A}\begin{pmatrix}c_{1}e^{3t}\\\ 0_{2}\\\
\end{pmatrix}=\mathbf{x}(t|A,\mathbf{x}_{0}).\end{split}$

######  Example S9 (An open dense set with arbitrarily small measure).

Let $\mathbb{Q}_{d\times d}\subset M_{d\times d}$ be the set of matrices with
rational entries. It is clear that $\mathbb{Q}_{d\times d}$ is a countable
set, so we can enumerate all elements in $\mathbb{Q}_{d\times d}$ as
$\left\\{Q_{1},Q_{2},\dots\right\\}$. Let

$\Omega_{a}:=\bigcup_{i=1}^{\infty}B(Q_{i},a\cdot 2^{-i}),\qquad
B(q,r):=\left\\{x\in M_{d\times d}:\,\|x-q\|_{F}<r\right\\}.$

Apparently, $\Omega_{a}$ is open and dense (because $Q_{d\times d}$ is dense)
in $M_{d\times d}$. Furthermore, the Lebesgue of $Q_{d\times d}$ can be made
arbitrarily small because

$\lambda_{d\times
d}\left(\Omega_{a}\right)\leqslant\sum_{i=1}^{\infty}\lambda_{d\times
d}\left(B(Q_{i},a\cdot
2^{-i})\right)=\sum_{i=1}^{\infty}c_{d}a^{d^{2}}2^{-d^{2}i}<c_{d}a^{d^{2}}.$

Here $c_{d}:=\frac{\pi^{d^{2}/2}}{\Gamma(\frac{d^{2}}{2}+1)}$ is a constant
that represents the volume of the Frobenius unit ball $B(0_{d\times
d},1)\subset M_{d\times d}$. Because $a$ is an arbitrary positive number,
$\lambda_{d\times d}\left(\Omega_{a}\right)$ can be made arbitrarily small.
Due to the absolute continuity between the probability of GinOE and Lebesgue
measure, we can also select $a$ so that
$P_{\mathrm{GinOE}}\left(\Omega_{a}\right)$ is smaller than any arbitrary
positive number.