Inverse Problems for Matrix Exponential in System Identification:
System Aliasing

arXiv:1605.06973v1 [cs.SY] 23 May 2016

Zuogon Yue, Johan Thunberg∗ and Jorge Gonçalves
Abstract— This note addresses identification of the A-matrix
in continuous time linear dynamical systems on state-space
form. If this matrix is partially known or known to have a
sparse structure, such knowledge can be used to simplify the
identification. We begin by introducing some general conditions
for solvability of the inverse problems for matrix exponential.
Next, we introduce “system aliasing” as an issue in the identification of slow sampled systems. Such aliasing give rise to nonunique matrix logarithms. As we show, by imposing additional
conditions on and prior knowledge about the A-matrix, the
issue of system aliasing can, at least partially, be overcome.
Under conditions on the sparsity and the norm of the A-matrix,
it is identifiable up to a finite equivalence class.

Then we consider the issue of system aliasing. When know
that the A-matrix is sparse, we observe that allowing for
“aliased” representations might lead to A-matrices that are
more sparse among the aliased solutions, thus exposing
more structures in the state generation. This note gives a
mathematical definition of “system aliases”, and study how
to select the sparest one among the aliases of the underlying
systems. Refer to [4] for more preliminaries and all proofs.
Let

I. I NTRODUCTION
Time-series models in engineering, economics and biology
can often be represented by state-space models. In the example of gene regulatory networks (evolving close to equilibria),
the structure of the right-hand side defines pathways, from
which conclusions can be drawn about possible diseases. For
continuous-time linear dynamics, the state-space form has an
A-matrix, which reveals the direct connections between the
states. In many dynamical systems, as the aforementioned
ones, the structure of this matrix is either partially known
or known to be sparse. This paper investigates what such
criteria can be used to identify the A-matrix.
Estimating continuous-time systems from discrete-time
measured data is an important part of the field of systems
identification, see e.g. [1]–[3]. However, with low sampling
rates, the identification of continuous-time systems becomes
particularly challenging, manifested in the lack of comprehensive studies. In the presence of “system aliasing”, the
discrete-time signals do not contain certain information about
the continuous-time signals. As a result, even though the
discrete-time system can be identified, the selection of the
continuous-time model may be ambiguous. This note sheds
light on how this issue of system aliasing complicates the
identification.
Before addressing the issues related to system aliasing,
we first recall some results on matrix logarithms and exponentials – the questions of existence and uniqueness are
addressed. The results are given as algebraic conditions for
obtaining a unique A-matrix within a set. In contrast to these
results, the note proceeds by first providing the minimal
sampling frequency such that system aliasing is avoided.

where Ā ∈ Rn×n and E ∈ S ⊆ Rn×n . This note addresses
properties that must hold for (A, D, S ) or (Ā, E, D, S )
in order to guarantee that E can be determined from
Dvec(exp(A)). In the following cited definitions and theorems, we adopt the notations in [5].

This work was supported by Fonds National de la Recherche Luxembourg
(9247977 & 8864515).
All authors are with Luxembourg Centre for Systems Biomedicine
(LCSB), University of Luxembourg, 6, avenue du Swing, L-4367 Belvaux,
Luxembourg.
∗ For correspondence, johan.thunberg@uni.lu

2

A = Ā + E and D ∈ Rp×n ,

A. Principal logarithm
Theorem 1 (principal logarithm [5, Thm. 1.31]). Let P ∈
Cn×n have no eigenvalues on R− . There is a unique logarithm X of P all of whose eigenvalues lie in the strip
{z : −π < im(z) < π}. We refer to X as the principal
logarithm of P and write X = Log(P ). If P is real then its
principal logarithm is real.
Let G(h) = {z ∈ C : −π/h < im(z) < π/h, h ∈ R}. We
denote the set of real matrices in Rn×n whose eigenvalues
lie in the strip G(1) by A (n). By restricting the set for which
Ā and E belong to A (n), it follows that
E = Log(exp(A)) − Ā,
is one-to-one. Throughout the text, the notations exp(·) and
e(·) are used interchangeably. We use log(·) for general primary matrix logarithms and Log(·) for principal logarithms.
Theorem 2 (Gantmacher [5, Thm. 1.27]). Let P ∈ Cn×n
be nonsingular with the Jordan canonical form
Z −1 P Z = J = diag(J1 , J2 , ..., Jp )

λk 1

.

λk . .
Jk = Jk (λk ) = 

..

.



(1a)



 ∈ Cmk ×mk .

1
λk
(1b)

Then all solutions to eA = P are given by
A = ZU diag(Lj11 , Lj22 , ..., Ljpp )U −1 Z −1 ,

(2)

where

2

Ljkk = log(Jk (λk )) + 2jk πiImk ;

log(Jk (λk )) denotes

f (λk ) f 0 (λk ) · · ·


..

.
f (λk )
f (Jk ) := 

.
..


f (mk −1) (λk )
(mk −1)!

..
.

f 0 (λk )
f (λk )

(3)

fĀ (E) = K(Ā, E)vec(E) and gĀ (E) = K(Ā, 0)vec(E).



II. P ROBLEM FORMULATIONS








A. General notation
Definition 2.

with f the principal branch of the logarithm, defined by
Im(log(z)) ∈ (−π, π]; jk is an arbitrary integer; and U
is an arbitrary nonsingular matrix that commutes with J.
Theorem 3 (classification of logarithms [5, Thm. 1.28]). Let
the nonsingular matrix P ∈ Cn×n have the Jordan canonical
form (1) with p Jordan blocks, and let s ≤ p be the number
of distinct eigenvalues of A. Then eA = P has a countable
infinity of solutions that are primary functions of P , given
by
(j )
−1
p)
,
(4)
Aj = Z diag(Lj11 , L2 2 , ..., L(j
p )Z
where Ljkk is defined in (3), corresponding to all possible
choices of the integers j1 , ..., jp , subject to the constraint
that ji = jk whenever λi = λk .
B. Fréchet Derivatives
Definition 1 ( [5]). The Fréchet derivative of the matrix
function f : Cn×n → Cn×n at a point X ∈ Cn×n is a linear
map
L
Cn×n −→ Cn×n
E 7−→ L(X, E)
such that for all E ∈ Cn×n

n
2
2
E (A, D, h, S ) = A∗ ∈ Rn×n : D ∈ Rn ×n , h ∈ R,
o
A∗ = arg min kDvec(exp(hA)) − Dvec(exp(hÃ))k2 ,
Ã∈S

where S ⊆ Rn×n contains A.

This general notation will be used to define an important
concept of system aliasing in continuous-time linear system
identification.
Consider the special case of the set S whose elements
A∗ are all in A (n). Moreover, to easily apply Fréchet
derivatives, with A expressed as A = Ā + E at Ā, we have
the following definitions.
Definition 3.

n
2
2
EL (Ā, E, D, S ) = E ∗ : Ā + E ∗ ∈ A (n), D ∈ Rn ×n ,

o
E ∗ = arg min kDvec(exp(Ā + E)) − Dvec(exp(Ā + Ẽ))k2 ,
Ẽ∈S

where S ⊆ Rn×n contains E.

If E is “sufficiently” small in norm, there is an approximated but easier problem that we could investigate, where
we can use Fréchet derivatives to approximate exp(Ā + E)
at Ā in the direction of E.
Definition 4.

f (X + E) − f (X) − L(X, E) = o(kEk).
The Fréchet derivative exists for matrix functions exp and
Log (principal logarithm) and it is unique. It holds that [5,
p. 238]
Z 1
Lexp (X, E) =
eX(1−s) EeXs ds,
0
Z 1
LLog (X, E) =
(t(X − I) + I)−1 E(t(X − I) + I)−1 )dt.
0

C. Fréchet derivatives for the vector representation

The vector representation of the Fréchet derivatives
Lexp (X, E) and LLog (X, E) have the structure of being
given as a matrix multiplied by vec(E).
vec(Lexp (X, E)) = K(X, 0)vec(E),

(5)

vec(LLog (X, E)) = K(X, 0)−1 vec(E),


K(X, E) = I ⊗ exp(X) ψ (X + E)T ⊕ (−X)) .

(6)
(7)

Here K(X, E) can be seen as an extension of the object
K(X) defined in [5, Thm. 10.13]. The vector representation
of eA = eĀ+E can be written as
vec(eĀ+E ) = vec(eĀ ) + K(Ā, E)E.

Moreover, let us define the maps fĀ : Rn×n → Rn and gĀ :
2
Rn×n → Rn by

(8)

n
2
2
ES (Ā, E, D, S ) = E ∗ : Ā + E ∗ ∈ A (n), D ∈ Rn ×n ,
 o
E ∗ = arg min kDK(Ā, 0) vec(E) − vec(Ẽ) k2 ,
Ẽ∈S

where S ⊆ Rn×n is a linear subspace containing E.

Here we use the first order approximation of the exponential matrix. Then we assume that S is a linear subspace
and we want to find out in what directions in this space the
first order approximation is good. To be more precise, the
question is for what (Ā, E, D, S) it holds that
!
sup

E 0 ∈ES (Ā,tE,D,S)

|t|

kE 0 − tEk

→0

ast → 0.

(9)

B. Continuous-time linear system identification
Consider the linear dynamical system
(
ẋ(t) = Ax(t),
y(t) = Cx(t),

(10)

where A = Ā + E; C ∈ Rp×n has full rank, x(t) ∈ Rn and
y(t) ∈ Rp .

With the general notation given in Section II-A, we can
give a definition on system aliasing only using the A matrix
and the sampling period h, which is not necessary to depend
on specific identification methods.

Im( i (A))

Im( i (A))h = ⇡

Definition 5 (System aliasing). Given A ∈ S and h ∈ R+ ,
if there exists Â 6= A ∈ E (A, I, h, S ) and Â is called system
alias of A with respect to S . By default, choose S (A) :=

Ã ∈ Rn×n : max{im(eig(Ã))} ≤ max{im(eig(A))} .

We are particularly interested in E (A, I, h, S ) = {A}, i.e.
there is no problem of system aliasing. Note that the concept
of system aliasing does not depend on specific data. It only
depends on system dynamics (e.g. the A-matrix in (10))
and sampling frequencies. If the D matrix is specifically
constructed by data instead of I, E (A, D, h, S ) = {A},
where A denotes the ground truth, tells that the underlying
system is identifiable from the given data (see Section III-B).
Obviously if we have system aliasing for the system with a
specific sampling frequency, without extra prior information
on A (see Section IV), the system is not identifiable.
III. N O SYSTEM ALIASING

1
2
3

1
8

.. 0
..
..

1 h
max
4

1
2

1 h

p

Im( i (A))h =

⇡

Fig. 1: The imaginary parts of all eigenvalues of A must
lie into (−π/h, π/h). λi (·) denotes the i-th eigenvalue of
A in Theorem 2. The symbols “×” denote the locations of
im(λi (A)). hmax is the maximal sampling period that allows
taking principal logarithms to estimate A, without facing
troubles from system aliasing.

A. The minimal sampling frequency
Provided with the definition of system aliasing, a
question comes first: for what (A, h) , it holds that
E (A, I, h, S (A)) = {A}.
To make principal matrix logarithm Log(·) well-defined,
assume that exp(hA) has no negative real eigenvalues. By
Theorem 1 and 2, it always holds that Log(exp(hA))/h ∈
E (A, I, h, S (A)). To avoid system aliasing, we have to
force Log(exp(hA))/h = A to be satisfied. It is equivalent
to eig(hA) ∈ G(1).
Given no other information on the system, consider the
identification problem of A using full-state measurement.
The only way to find the unique estimation is to decrease the
sampling period h until the ground truth falls into the strip
of G(h), and then use the principal logarithm, as illustrated
in Figure 1. Otherwise, we would be bothered by system
aliases of A and unable to make a decision, unless we know
extra prior information on A. For full-state measurement,
identifiability is guaranteed by selecting appropriate h such
that there is no system aliases. For the general case of
identification using output measurement, the issue is studied
in Section III-B.
Theorem 4 (Nyquist-Shannon-like sampling theorem). To
uniquely obtain A from Ad by taking the principal matrix
logarithm, where Ad is identified from sampled data, the
sampling frequency ω (rad/s) must satisfy
ω ≥ 2 max {| im (λi (A)) |, i = 1, . . . , n} .
Equivalently, the sampling period h should satisfy
h ≤ min {π/| im (λi (A)) |, i = 1, . . . , n} .
Proof. The proof immediately follows from the above explanation on principal matrix logarithms.

B. Partial information
Suppose all A’s in A (n), which implies there is no system
aliasing, i.e. the case in Definition 3. Now consider the
identifiability problem of (10) from data with low sampling
frequency 2π/h. To be precise, it is to find out for what
(Ā, E, D, S ), it holds that EL (Ā, E, D, S ) = {E}.
Lemma 5. For Ā, E, D, S 3 E, if there is linear subspace
L such that fĀ (S ) ⊆ L , then EL (Ā, E, D, S ) = {E} if
L ∩ ker(D) = {0}.

Lemma 6. For Ā, E, D and S 3 E, (9) holds if and only
if
−1
gĀ
(ker(D)) ∩ S = {0}.

Proposition 7. If

S ⊆ {Ẽ : im(Ẽ) ∪ im(ĀẼ) ∪ im(Ā2 Ẽ) ∪ . . .
and

∪ im(Ān−1 Ẽ) ⊆ im(C T )}.
D = (X T ⊗ C),

where X ∈ Rn×n is a non-singular matrix, then
EL (Ā, E, D, S ) = {E} and
!
 −1  A 
Ce X
C
E = log
X −1 − Ā,
Z
ZeĀ X
where Z is any matrix in R(n−p)×n such that im(Z T ) =
ker(C).
Proof.
Z 1
im(
eĀ(1−s) E 0 U (s)ds) ⊆ span{E 0 , ĀE 0 , Ā2 E 0 , . . . , Ān−1 E 0 }.
0

Now, if
S ⊆ {Ẽ : im(Ẽ) ∪ im(ĀẼ) ∪ im(Ā2 Ẽ) ∪ . . .
∪ im(Ān−1 Ẽ) ⊆ im(C T )},

for any E 0 ∈ S , it holds that K(Ā, E 0 )vec(E 0 ) ∈ im(I ⊗
C T ). But im(I ⊗ C T ) = im(DT ). Now, according to
Lemma 5, EL (Ā, E, D, S )R = {E}.
1
Furthermore, since im( 0 eĀ(1−s) Ee(Ā+E)s ds) ⊆ C T ,
ZeA = ZeĀ . It follows that
 
 A 
Ce X
C A
e X=
Z
ZeĀ X
Corollary 8. If
2

S ⊆ {Ẽ : im(Ẽ) ∪ im(ĀẼ) ∪ im(Ā Ẽ) ∪ . . .
∪ im(Ān−1 Ẽ) ⊆ im(X)}.

and
D = (C ⊗ X T ),
where p = n, C is a non-singular matrix and X ∈ Rn×k
where k ≤ n, then EL (ĀT , E T , D, S ) = {E} and
!T
 T −1 " T AT T #
X e C
X
C −T
,
A = log
T
Z
ZeĀ C T
where Z is any matrix in R(n−k)×n such that im(Z T ) =
ker(C).
Proof. Identify Ā with ĀT , E with E T , C with X T , and
X T with C in Proposition 7.
Remark 1. If the columns of X are linearly independent
and each column is an initial point for the system (10),
Corollary 8 can be used to provide an upper bound on the
number of measurements needed in order to guarantee that
EL (Ā, E, D, S ) = {E}.

Proposition 9. Suppose the system (10) is initialized at k
different initial points x0,1 , x0,2 , . . . , x0,k ∈ Rn , C ∈ Rp×n
is a full rank matrix and S is an l-dimensional linear
subspace. At the time t = 1, for initial point i, yi (1) =
CeA x0,i , where i ∈ {1, 2, . . . , k}.
For
almost
all
linearly
independent
vectors
x0,1 , x0,2 , . . . x0,k in Rn and almost all matrices Ā in
Rn×n , (9) holds where

T
D = x0,1 x0,2 . . . x0,k ⊗ C,

if and only if

l ≤ kp.
Proof. For linearly independent x0,1 , x0,2 , . . . x0,k ,

T
D = x0,1 x0,2 . . . x0,k ⊗ C

has rank equal to kp. Let {v1 , v2 , . . . , vl } be a basis for
S . Each vi can be written as vi = viker + viim where
K(Ā, 0)viker ∈ ker(D) and K(Ā, 0)viim ∈ im(DT ). For

almost all linearly independent x0,1 , x0,2 , . . . x0,k and Ā it
holds that the viim are linearly independent if and only if
l ≤ kp.
−1
gĀ
(ker(D)) ∩ S = {0}
if and only if the vi are linearly independent. Thus, according
to Proposition 6, (9) holds.
Remark 2. In Proposition 9, instead of having k different
initial points one can have one intial point and sample y(t)
at the times t = 1, t = 2 etc. for almost all Ā, E and x0 such
that [x0 , eA x0 , e2A x0 , . . . , e(k−1)A ] has full rank (9) holds
when


D = ( xT0 (eA x0 )T . . . (e(k−1)A x0 )T ⊗ C,
l ≤ kp and S is an l-dimensional linear subset.

IV. S YSTEM ALIASING AND BOUNDED CONSTRAINTS
In the previous section we hinted that the conditions on
no system aliasing follows as a consequence of bounded
eigenvalues. In this section we follow this path and explicitly
formulate an optimization problem to deal with identification
in the presence of system aliases.
Consider the case of full-state measurement, i.e. C = I
in (10), and the sampling period h is NOT chosen small
enough such that E (A, I, h, S (A)) = {A}. Then finding
out A among its aliases need extra information, for instance,
properties of A known a priori. Here assume that the ground
truth A is the sparest solution in E (A, I, h, S (κ)) and κ ∈ R
as an upper bound can be roughly estimated, where S (κ)
will be defined after giving Definition 6. It implies that A is
chosen by solving the following optimization problem
minimize
Â∈E (A,I,h,S (κ))

kÂk0 .

(11)

We need to calculate E (A, I, h, S (κ)) from data. Given
the measurement X1 = [x(h),
 x(2h), . . . , x(N h)], X2 =
[x(0), x(h), . . . , x (N − 1)h ], let Âd be an estimation of
the A-matrix in the corresponding discrete-time state space
representation. In the deterministic case1 as (10), Âd =
X1 X2T (X2 X2T )−1 and Âd = exp(hA) . Hence,
n
o
E (A, I, h, S (κ)) = Ã ∈ S (κ) : exp(hÃ) = Âd , (12)
and define

n
o
S := Ã ∈ Rn×n : exp(hÃ) = Âd .

To formulate S (κ), we need to introduce a special norm
of A, which is equivalent to the Frobenius norm up to a
change of coordinates.
Definition 6 (Z-weighted norm). Let hZ (A) = Z −1 AZ,
where Z is the matrix defined in Theorem 3. Then the norm
is defined as khZ (·)kF = k · kF ◦ hZ .
1 For stochastic cases, Â
d is consistently estimated by Prediction Error Minimization or Maximum Likelihood methods [1], and
limN →∞ E(Âd (N )) = exp(hA). If only finite samples are available,
we cannot obtain the exactly equivalent E (A, I, h, S (κ)) from data.

Since we assume that Âd is fixed, i.e., the data X is not
used in the optimization problems defined here, the matrix
Z is constant. One can observe that
khZ (Â)kF = vec(Â)T (Z T ⊗ Z −1 )T (Z T ⊗ Z −1 )vec(Â)
is a proper (Z T ⊗ Z −1 )T (Z T ⊗ Z −1 )-weighted vector norm
in terms of vec(Â). Using khZ (·)kF is on the one hand
simplifying the analysis we conduct throughout this section,
and on the other explicitly penalizes the imaginary part
of the eigenvalues without “distorting” them through the
transformation by Z.
Now we define S (κ) using the norm khZ (·)kF . The basic
idea is that one should exclude such A’s whose imaginary
parts of eigenvalues are too large, which implies their system
response will show wild fluctuation. To make our assumption
and the problem (11) practically meaningful, instead of
Rn×n , we restrict S to be a norm bounded subset
n
o
S (κ) = Ã ∈ Rn×n : khZ (Ã)kF ≤ κ .
(13)

In the following we will show that the feasible set of
(11) has only finite elements, which implies it can be solved
at least by brutal force methods. Recall that the set S is
countable according to Theorem 3.
Let M := diag(m1 , m2 , . . . , mp ), j := [j1 , j2 , . . . , jp ] and
β := [β1 , β2 , . . . , βp ], where log(λk ) , αk + iπβk , k =
1, . . . , p, and jk , λk are defined in Theorem 2. A function
I is defined as
I (j, δ) := δ T M δ + (2j + β)T M δ,

(14)

where j, δ ∈ Zp . Moreover, it satisfies I (j, δ) = I (0, j +
δ) − I (0, j), which follows by noticing
T

I (j, δ) = (δ + j + β/2) M (δ + j + β/2)
T

− (j + β/2) M (j + β/2) .

(15)

Moreover, let A0 denote a special matrix logarithm for which
all jk (k = 1, . . . , p) in (3) are equal to 0.
Definition 7 (equivalence relations). An equivalence relation
“∼” is defined on S as a binary relation: for any A1 , A2 ∈ S,
j(1) and j(2) are defined for A1 , A2 , respectively, we say
A1 ∼ A2 if I (j(1) , j(2) − j(1) ) = 0.
Lemma 10. Let S be the set defined in (12) and
parametrized by (4) in Theorem 3. For any A1 , A2 ∈ S,
khZ (A1 )kF = khZ (A2 )kF if and only if A1 ∼ A2 .
Proof. Let
(i)

j
Lkk

j (i)
Ai := Z diag(L1 , · · · , Lpp )Z −1 , where i =
:= log(Jk (λk )) + 2jk(i) πiImk , and all other
(i)
j1

1, 2,
notations are given in (4). By using (16) (see next page)
for A1 , A2 , we obtain
khZ (A1 )kF = khZ (A2 )kF ⇔ khZ (A1 )k2F −
= khZ (A2 )k2F − khZ (A0 )k2F ⇔ I (j(1) , j(2)

khZ (A0 )k2F
(1)
−j

) = 0,

which implies that A1 ∼ A2 by definition. The first equality
in (16) is due to the linear transformation hZ (Â).

Remark 3. It is not necessary that A0 is the principal matrix
logarithm (consider the case when the principal logarithm
does not exist), nor does it have to be the logarithm with the
smallest (weighted) Frobenius norm.
Lemma 11. Given any Ā ∈ S, there exist finite Ai ∈ S that
satisfies Ai ∼ Ā.

Proof. Let j denote [j1 , . . . , jp ] of Ā in (3), and j(i) denotes
(i)
(i)
[j1 , . . . , jp ] of Ai ∈ S. δ , j(i) − j, therefore δ ∈ Zp ,
where  denote the element-wise larger-or-equal relation. By
Definition 7, it is equivalent to show that I (j, δ) = 0 has
finite solutions, given j. We require δ to satisfy the following
condition:
s
(j + β/2)T M (j + β/2)
(17)
|δi + ji + βi /2| ≤
mi
for all i = 1, . . . , p. Otherwise, supposing that there exists
i ∈ {1, . . . , p} such that δi does not satisfy (17), we will
have
X
I (j, δ) = mi (δi + ji + βi /2)2 +
mk (δk + jk βk /2)2
−

X

k6=i

mk (jk + βk /2)2 >

k

X
k6=i

mk (δk + jk βk /2)2 ≥ 0.

(i)
Let S 0 := {Ai ∈ S : jk = δk + jk , δk satisfies (17)}. We
have {Ai ∈ S : Ai ∼ Ā} ⊆ S 0 and S 0 is a finite set.

Lemma 12. There exists finite Ai
khZ (Ai )kF ≤ κ.

∈

S such that

Proof. Let κ0 , khZ (A0 )kF . Then we need to show there
exists a finite number of Ai ∈ S such that khZ (Ai )k2F −
khZ (A0 )k2F ≤ κ2 − κ20 , which is equivalent to show that
there exists a finite number of solutions δ ∈ Z to I (0, δ) ≤
(κ2 − κ20 )/4π. δ must satisfy the following condition:
s
(β/2)T M (β/2) + (κ2 − κ20 )
|δi + βi /2| ≤
(18)
mi
for all i = 1, . . . , p. Otherwise, by supposing that there exists
i ∈ 1, . . . , p such that δi does not satisfy (18)leads to
I (0, δ) = mi (δi + βi /2)2 +
X
mk (δk + βk /2)2 − (β/2)T M (β/2)
k6=i

>

X
k6=i

mk (δk + βk /2)2 + (κ2 − κ20 ) ≥ κ2 − κ20 .

Note that the set of all δ ∈ Z that satisfies (18) is finite,
which finalizes the proof.
Remark 4. We couldPhave a more precise bound of δi . Let
2
µ(i) := minimize
k6=i mk (δk + βk /2) . Then δi (i =
δ

1, . . . , p) should satisfy
s
(β/2)T M (β/2) + (κ2 − κ20 ) − µ(i)
|δi + βi /2| ≤
. (19)
mi



(i)
(i)
j (i)
j (i)
j
j
khZ (Ai )k2F − khZ (A0 )k2F = tr diag∗ (L11 , · · · , Lpp ) diag(L11 , · · · , Lpp )


(0)
(0)
(0)
− tr diag∗ (L1 , · · · , L(0)
)
diag(L
,
·
·
·
,
L
)
p
p
1

 (i)
p
(i)
X
j ∗ j
(0)∗ (0)
=
tr Lkk Lkk − Lk Lk
=
=

k=1
p
X

k=1
p
X

(16)



(i)
(i)2
tr 2jk πi (log(Jk )∗ − log(Jk )) + 4π 2 jk Imk
(i)

(i)

4πjk mk (βk + jk ) = 4π I (0, j),

(i)

j , [j1 , . . . , jp(i) ].

k=1

Moreover, we have the solution to µ(i), i = 1, . . . , p:
µ(i) = min{I (0, δ/i ) :(δ/i )k = d−βk /2e or b−βk /2c,
k 6= i; (δ/i )i = 0}.

(20)

Proposition 13 (lower boundness of logarithms). Let S be
the set defined in (12). Given any Ā ∈ S, there exists
M (Ā) > 0, such that for any A ∈ {A ∈ S : A  Ā},
it holds that
khZ (A)kF − khZ (Ā)kF ≥ M.
Proof. Let j denote [j1 , . . . , jp ] of Ā in (3), Neqiv be the
number of A’s that satisfy A ∼ Ā. Note that |khZ (A)k2F −
khZ (Ā)k2F | = |(khZ (A)k2F − khZ (A0 )k2F ) − (khZ (Ā)k2F −
khZ (A0 )k2F )| = |I (j, δ)|, δ ∈ Z, which implies it is
equivalent to show that |I (j, δ)|, δ ∈ Z has a non-zero lower
bound if not considering the δ’s that result in I (j, δ) = 0.
We will prove it by contradiction. Assume this is not true,
i.e. ∀ > 0 there exists δ such that 0 < |I (j, δ)| < . It
implies that, arbitrarily given  > 0, there exists an infinite
number of δ such that I (j, δ) < , which is impossible since
I (0, j + δ) < I (0, j) +  (using the fact that I (j, δ) =
I (0, j + δ) − I (0, j)) has a finite number of solutions
provided by Lemma 12.
Proposition 14. Let S be the set defined in (12). For any
Ā ∈ S, there exist κl , κu ∈ R in S (κl , κu ) = {Ã ∈ Rn×n :
κl ≤ khZ (Ã)kF ≤ κu } such that (11) has a unique optimal
point in the sense of the equivalence relation in Definition 7.
Proof. It immediately follows by choosing
κl > max{0, khZ (Ā)kF − M (Ā)},

κu < khZ (Ā)kF + M (Ā),

where M (Ā) is the lower bound on the gap between Ā and
any A  Ā ∈ S, defined in Theorem 13.
V. C ONCLUSIONS
This paper addresses identification of continuous-time dynamical systems with sparse topologies. The key assumption
is that the sampling frequency is low. Under this assumption
a realization/identification problem comes to surface, which
has largely been overlooked in the community. First we

propose the minimal sampling frequency that guarantees
no system aliasing. Allowing system aliasing, one needs
to search over a collection of matrix logarithms to find
the sparsest one. We provide theoretical results for when a
unique solution exists up to a finite equivalence class.
A PPENDIX
P ROOF OF LEMMAS
A. Proof of Lemma 5
Proof. Suppose E 0 6= E and E 0 + Ā ∈ A (n), E 0 ∈ S .
0
Using (8) and the fact that eĀ+E 6= eĀ+E when Ā + E 0 ∈
A (n), we know that fĀ (E) 6= fĀ (E 0 ). Now, since
L ∩ ker(D) = {0}
it holds that DfĀ (E) 6= DfĀ (E 0 ), otherwise DfĀ (E) =
DfĀ (E 0 ) and Prker(D) (fĀ (E)) 6= Prker(D) (fĀ (E 0 )), which
means that
L ⊇ span{fĀ (E 0 ) − fĀ (E)} ⊆ ker(D).
Now,
kDvec(exp(Ā + E)) − Dvec(exp(Ā + E 0 ))k2

=kD(fĀ (E) − fĀ (E 0 ))k2 6= 0.
B. Proof of Lemma 6
Proof.
Only if: Suppose there is

E 0 ∈ (fĀ−1 (ker(D)) − {0}) ∩ S .
By using the fact that S is a linear subspace, one can show
that it holds that, E ∗ + sE 0 ∈ ES (Ā, tE, D, S ) for all s if
E ∗ ∈ ES (Ā, tE, D, S ). Now it is easy to show that (9) does
not hold.
If: Let
b(t) = DK(Ā, 0)vec(tE),
q(t) = D(vec(eA )) − D(vec(eĀ )) − b(t),
Q = DK(Ā, 0),
x = vec(Ẽ).

S is a linear subspace, let us assume its dimension is k.
There is a full rank matrix P ∈ Rn×k such that
The problem

S = {y : y = P v, v ∈ Rk }.

min kDvec(eA ) − Dvec(eĀ ) − DK(Ā, 0)vec(Ẽ)k2 },

Ẽ∈S

can equivalently be written as
(
min xT QT Qx − 2(b(t) + q(t))T Qx
s.t.
x = P v.
The solution to this problem is
x = P (P T QT QP )−1 P T QT (b(t) + q(t)).
Note that P T QT QP = P T K(Ā, 0)T DT DK(Ā, 0)P is
invertible since im(K(Ā, 0)P ) ∩ ker(D) = 0, K(Ā, 0) is
invertible and P has full rank. These facts combined with
the fact that P has no right-nullspace can be be used to
show that ES (Ā, tE, D, S ) contains only one element. It
holds that vec(tE) = P (P T QT QP )−1 P T QT (b(t)). Also,
q(t) is O(t2 ); (9) holds.
R EFERENCES
[1] L. Ljung, System Identification: Theory for the User, ser. Prentice-Hall
information and system sciences series. Prentice Hall PTR, 1999.
[2] H. Garnier and L. Wang, Identification of Continuous-time Models from
Sampled Data. Springer London, 2008.
[3] N. K. Sinha and G. P. Rao, Identification of continuous-time systems:
Methodology and computer implementation.
Springer Science &
Business Media, 2012, vol. 7.
[4] Y. Zuogong, T. Johan, and G. Jorge, “Inverse Problems for Matrix
Exponential in System Identification: System Aliasing,” arXiv preprint,
2016.
[5] N. Higham, Functions of Matrices. Society for Industrial and Applied
Mathematics, 2008.

