Systems Aliasing in Network Reconstruction:
Issues in low sampling frequencies

arXiv:1605.08590v3 [cs.SY] 28 Sep 2017

Zuogong Yue1 , Johan Thunberg1∗ , Lennart Ljung2 and Jorge Gonçalves1

Abstract—Dynamic network reconstruction is motivated by
applications from many fields. Due to the limits in experiments, especially in biology, it is required to use low-samplingfrequency data to perform network inference. The paper raises
the concept of “system aliasing” which plays a central role
when the sampling frequency is not high enough. It focuses
on the full-state measurement, which is supposed to the basic
case but has fairly shown complications. A sampling theorem
on system aliasing is presented, together with a test criterion.
In the cases of no system aliasing, an algorithm is provided to
reconstruct dynamic network from data measured with large
sampling periods. Moreover, it provides certain results in theory
for the cases when system aliases are presented. As the first
tentative study, the paper dedicates to manifest the intrinsic
complexity in this problem.

I. I NTRODUCTION
Many complex systems can be modeled as networks in
applications to reveal interactions between variables, and the
sparse structure is often presented in nature. An example of
the latter is the interactions between species such as genes
and proteins in human cells; such interactions can be modeled
by stochastic/ordinary differential equations, e.g. [1]. Such
network models in biology help to understand, for instance,
metabolic pathways, interactions between DNAs/proteins, and
furthermore contribute to pathology of, disease detection on
or even clinical treatment to complicate diseases, e.g. [2]. In
order to perform such an analysis, reconstruction of sparse
(Boolean) networks turns to be critical as more techniques
have been available to acquire time-series data.
There have been many papers on the study of dynamic
networks in different aspects: network identifiability [3],
[4], network module identifiability [5], network inference in
discrete-time approaches [6], [7], etc. With regard to network
inference, the factor that distinguishes itself from traditional
system identification [8] is the particular request on sparse
structures. To enhance sparsity, there are multiple methods
are available: LASSO [9], iterative reweighted l1 /l2 algorithms
[10], [11], Sparse Bayesian Learning [12], [13], etc.
It deserves to be emphasized that the discrete-time approach
for network inference is valid only if the sampling frequency is
high enough, where the discrete-time model shares the same
network structure as the continuous one that is the physical
process (here we assume the dynamical systems evolve in
This work was supported by Fonds National de la Recherche Luxembourg
(Ref. 9247977).
1 Zuogong Yue, Johan Thunberg and Jorge Gonçalves are with Luxembourg
Centre for Systems Biomedicine, 6, avenue du Swing, L-4367 Belvaux,
Luxembourg.
2 Lennart Ljung is with Department of Electrical Engineering, Linköping
University, Linköping, SE-58183, Sweden
∗ For correspondence, johan.thunberg@uni.lu

continuous time in nature). To use discrete-time methods, one
practical rule to choose sampling frequencies is taking ten
times the bandwidth of the underlying linear systems [8].
However, in biological systems, most time-series data are
sampled considerably slower than this empirical frequency,
e.g. “high time-resolution” time series in [14], which often
cannot be solved by increasing sampling rates due to strict
limits in biological experiments.
There have been papers studying the identification of
continuous-time systems, e.g. [15], [16]. However, most methods request a high sampling frequency to guarantee certain
simplifications on theoretical deduction or numerical calculations. Moreover, choosing too low sampling frequencies may
trigger the problem of “system aliasing”, that is, multiple
continuous-time systems produce the exactly same output
samples; meanwhile they show different network structures.
What’s more, it becomes particularly challenging in theory
and numerical computation to enhance sparsity on the network structures of continuous-time models. However, this is
inevitable since we have to resort to continuous-time models
to explore network structures due to the low sampling frequencies.
In this paper, we first reveal the challenges due to the
low sampling frequency by examples in Section III and then
present a definition of system aliasing. A Nyquist-Shannonlike sampling theorem was presented in Section IV to determine the minimal sampling frequency that avoids the effect of
system aliasing. Section V presents an algorithm to reconstruct
sparse networks in the case of no system aliasing using
low-sampling-frequency data. The case of system aliasing is
discussed in Section VI. The paper is closed by numerical
examples in Section VII.
II. P ROBLEM F ORMULATION
Consider a filtered probability space (Ω, F, (Ft )t∈[0,∞) , P),
where the filtration is always assumed to be complete. Let
{w(t) : t ≥ 0} be the n-dimensional standard Ft -Brownian
motion. The physical plant/process in our study, as a dynamical system in continuous time, is modeled by the following
stochastic differential equation
dx = Ax dt + Bu dt + Rdw,

(1)

where A ∈ Rn×n is stable, R ∈ Rn×n is symmetric
and semi-positive definitive, the initial x(t0 ) is a Gaussian
random variable with mean m0 and variance R0 , t0 ≥ 0, and
w(t) is interpreted as disturbance on the state variables (or
called process noise). The solution to (1) is an Ft -adapted n-

dimensional stochastic process x(t) = (x1 (t), · · · , xn (t))t≥t0
such that
Z t
Z t
R1 dw(s),
(Ax(s) + Bu(s)) ds +
x(t) = x(t0 ) +
t0

t0

where x(t0 ) ∈ Ft0 and see [17] for the definition of stochastic
integral. It is assumed that x(t0 ), w(t) are independent. The
solution x(t) is strong, that is, x(t) is adapted to Ftw := σ(wu :
u ≤ t) ∨ σ(N ) (i.e. the complete σ-field generated by {wu :
u ≤ t}; see [17] for details). An input signal {u(t), t ≥ t0 } has
been applied to the system, and the output y of the system is
observed at the discrete times t0 , t1 , ..., tN ,
y(tk ) = Cx(tk )

(2)

where C = [I 0] ∈ Rp×p , p ≤ n, tk , t0 + kh and h > 0
is the sampling period. Here the measurement noise is not
included mainly due to that we have not yet given a definition
of network models (see [3]) from state-space representations
with measurement noises. The stochastic difference equation
that relates the values of the state variable x in (1) at the
sampling instants [18, p. 82-85] [19, chap. 2] is given by
x(tk+1 ) = Ad x(tk ) + Bd u(tk ) + v(tk ),
where
Ad = exp(hA),

Bd =

Z

(3)

h

exp(sA)B ds,

(4)

0

exp(·) is the matrix exponential, and the Gaussian i.i.d. v(t)
has mean zero and covariance matrix
Z h
Rd =
exp(sA)R exp(sAT ) ds.
(5)
0

The linear dynamic network model of (1) is given as
y(t) = Q(q)y(t) + P (q)u(t) + H(q)e(t),

(6)

where Q(q), P (q) and H(q) are p×p, p×m and p×p matrices
of strictly-proper real-rational transfer functions respectively
in terms of q, q is the differential operator qx(t) = dx/dt,
and e(t) is the Gaussian white noise with zero mean and
E[e(t)T e(s)] = Iδ(t−s) (e.g. see [4]). The model (6) is called
Dynamical Structure Function (DSF), firstly proposed in [3].
The network model defines path diagrams which show the
interconnections between the elements of the output variable.
Definition 1 ( [7]). Let G = (V, E) be a digraph, where the
vertex set V = {y1 , . . . , yp , u1 , . . . , um } and the arc (directed
edge) set E is defined by
i) (yj , yi ) ∈ E ⇔ Qij (q) 6= 0,
ii) (uk , yi ) ∈ E ⇔ Pik (q) 6= 0,
iii) (yi , uk ) ∈
/ E, ∀i, k.
Let f be a map defined as
f:

E →
(yj , yi ) 7→

STF
Qij (q)

or

(uk , yi ) 7→ Pik (q),

where STF is a subset of single-input-single-output (SISO)
(strictly) proper real rational transfer functions. We call the
tuple N := (G, f ) a (linear) dynamic network, f the (linear)
dynamic capacity function of N , and G the underlying digraph
of N , which is also called (linear) Boolean dynamic network.

This article focuses on the full-state measurement case, i.e.
C = I, where A in (1) coincides with Q(q) in network
models (6). Concerning the network identifiability [3], we
assume B to be diagonal1 (DSF, [3]) or particularly B = I
(the model used in [5]). Let the measurement be denoted by
Y N , [y(t0 ), y(t1 ), . . . , y(tN )]. We summarize the main
problem in our study as follows:
Main Problem: Given the finite signal Y N in full-state
measurement (i.e. C = I), with probably large h (the sampling
period) and small N (the length of time series), infer the
dynamic network N (or the Boolean G), assuming that the
ground truth A is sparse and and B is diagonal.
Remark 1. The problem is challenging due to the following
two major reasons:
• Since h could be large, i.e. the sampling frequency is low,
we have to estimate A in order to determine G or N (see
Figure 1).
• Since N does not approach infinity, the estimation of
Ad from PEM (Prediction Error Minimization) or ML
(Maximum Likelihood) fails to identify G correctly by
taking matrix logarithms, even though PEM/ML gives
consistent estimation in theory.
Throughout the text, by default, we always deal with primary matrix functions, including exp (matrix exponential), log
(matrix logarithm) and Log (the principal matrix logarithm in
Theorem 4). Primary matrix functions refer to the ones defined
via Jordan Canonical Form or equivalently via Polynomial
Interpolation, Cauchy Integral Theorem [20, chap. 1]. The
primary notion of matrix functions is of particular interest and
the most useful in applications [20], [21].
III. S YSTEM A LIASING IN I DENTIFICATION
A. Observations on matrix logarithm
Supposing that Ad has been perfectly estimated from samples, the estimate of the A matrix for the continuous-time
system is straightforwardly calculated by solving
exp(hA) = Ad .

(7)

via matrix logarithm. However, referring to Theorem 18 [20],
the equation (7) has several (in fact infititely many) solutions.
Let us review the following observations on (7) to see the
troubles from low sampling frequencies (i.e. 1/h).
Observation 1: With the increase of h, the Boolean structures
of A (i.e. G determined from A) and Ad − I become more
and more different, as illustrated by Figure 1. The sampling
frequency (1/h) deserves to be emphasized as a core factor in
the categorization of different cases in our study:
• Case I: when h is “very small” such that Ad shares the
same Boolean structure as A + I. Indeed, one can see
2
it by exp(hA) = I + hA + h2! A2 + · · · . Hence we can
determine G by identifying discrete-time models;
• Case II: when h is “large” but the ground truth A is still
the principle matrix logarithm of Ad ;
1 This is not required in system identification, in which only the input-output
behavior is concerned. Any state-space realization could be feasible solutions.

Case III: when h is “even larger” such that the ground
truth A is no longer the principle logarithm of Ad .
The general network model (6) of Case I has been solved by
discrete-time approaches , e.g. see [6], [7]. Case II is what we
mainly studied in this paper. We call both Case I and II no
system aliasing, as defined and studied in later sections.
percentage of correct arcs (%)

•

A2

2
3
6660.0001 0.0001 0.1188 0.0001 777
666
7
6660.4982 0.0002 0.9595 0.3403 77777
6660.0000 0.0000 0.0002 0.0001 777
664
77
0.0001 0.5851 0.2239 0.00005

Fig. 4: h = 1, A1 = Log Ad1 /h and A2 = Log Ad2 /h.

90

2
666
666
666
666
664

80
70
60

Case I

Case II

Case III

50
40
30
10-2

10

0

10

2

10

4

10

6

Fig. 1: This example randomly chooses a sparse A and Ad =
exp(hA). The Boolean networks are determined from A and
Ad −I by Definition 1 (treating A and Ad −I as Q’s), denoted
by G(A), G(Ad − I). We use G(A) as the ground truth. Then
G(Ad − I) is compared with G(A) and the same arcs are
labeled as the correct ones.
Observation 2: Provided with a sparse A, the corresponding
Ad could be dense, as the example in Figure 2.

A

3
0.0149 0.0790
0 777
7
0.1932 0.1741 0.03137777
7
0.0325
0
0 77777
0.0020
0
0 5

2
6661.0000
666
6660.0004
6660.0000
664
0.0231

Ad
0.0178
1.2164
0.0359
0.0024

0.0804
0.1923
1.0030
0.0011

3
0.0003777
7
0.03467777
7
0.000577777
1.00005

Fig. 2: h = 1, Ad = exp (hA), A = Log Ad /h.

Observation 3: Provided with a sparse A, the corresponding
Ad can be also sparse. However, they have different Boolean
structures (i.e. zeros at different positions), as shown in
Figure 3.
2
666 0
666
6660.8611
6666 0
64
0.6714

2
3
666 0
0
0.1190
0 777
666
7
0
0.9597 0.34047777
6660.4984
777
6666 0
0
0
0
777
64
0
0.5853 0.2238
0 5

100

sampling periods (h)

2
666 0
666
666 0
666 0
664
0.0231

kAd2 Ad1k2
= 0.0138%
kAd1k2

A1

A

3
0
0
0.5303777
7
0 0.4849
0 7777
7
0 0.3935
0 77777
0
0
0 5

Ad

2
3
6661.1834
0
0
0.5624777
666
7
6660.9132 1.0000 0.5941 0.235277777
666 0
0
1.4821
0 77777
664
0.7120
0
0
1.18345

Fig. 3: h = 1, Ad = exp (hA), A = Log Ad /h.

Observation 4: Even though the norm difference of Ad1 , Ad2
has been very small, their matrix logarithms, e.g. A1 and A2
in Figure 4, have significantly different Boolean structures.
Observation 5: There exists more than one solution, e.g. A1
and A2 in Figure 5, and they have different Boolean structures.
Considering the problem formulation in Section II, one may
have already noticed the troubles on network reconstruction,
which originate from the matrix logarithms, due to the low
sampling frequency. The examples in Observation 1 clearly
show that why we have to resort to the continuous-time system

A1
3
3.7147 0.1778 0.1625 0.2690 777
7
0.8983 5.4242 0.5625 0.9199 7777
7
1.7815 0.2501 4.7786 1.122177777
0.8726 0.0043 0.6277 3.91625
A2

2
666 9.1872 1.7824
666
0
5.3854
666
666 16.4473 13.2820
664
13.7279 8.6466

3
2.9906 3.1019 777
7
8.4002 4.01127777
7
0
7.482877777
2.7852 3.26125

2
666
666
666
666
664

Ad
3
0.0300 0.0025 0.0023 0.0038 777
7
0.0127 0.0059 0.0079 0.0130 7777
7
0.0252 0.0035 0.0150 0.015877777
0.0123 0.0001 0.0089 0.0272 5

Fig. 5: h = 1, exp (hA1 ) = exp (hA2 ), A1 = Log Ad /h, and
A2 is the ground truth A.
identification to infer network structures. Observation 2 and 3
tell that there is no consistent relation between the sparsity of
A and Ad . Observation 4 points out that the Boolean structures
of the principle logarithms of two Ad1 and Ad2 close in
matrix norms could be significantly different. The example on
Observation 5 shows an even worse case: the sample period
is so large that the principle matrix logarithm is no longer A,
which appears as other branches of matrix logarithm of Ad ,
in which no robust algorithm has yet been available.
Remark 2. In a sum of the above observations, it tells us that,
in the identification of G or N ,
•
•

•

G should be determined from A instead of Ad , when h
is large (w.r.t. A); (see Observation 1, 4)
the sparsity penalty has to be imposed on A directly
instead of Ad , when h is large (w.r.t. A); (see Observation 2, 3, 4)
the A matrix should be estimated directly instead of via
taking matrix logarithm of Ad , in the presence of noise
and a limited length of signals. (see Observation 4)

B. Definitions
As shown in Section III-A, the A-matrix has to be identified
in network reconstruction when the sampling frequency is
low. In this scenario, a “good” case is that the ground truth
A is the principle matrix logarithm of Ad (i.e. Case II);
otherwise, it becomes particularly challenging (i.e. Case III),
e.g. Figure 5. To clarify this classification, we consequently
present an important concept in network reconstruction with
low sampling frequencies, “system aliasing”.
Let vec(X) denote the vectorization of the matrix X formed
by stacking the columns of X into a single column vector; and
ivec(·) is defined by ivec(vec(X)) = X. Im(x) denotes the
imaginary part of the complex number or vector x.

Definition 2.

n
2
2
E (A, M, h, S ) = A∗ ∈ Rn×n : M ∈ Rn ×n , h ∈ R,

o
A∗ = arg min kM vec(exp(hA)) − M vec(exp(hÃ))k2 ,
Ã∈S

be unable to make a decision, unless we know extra prior
information on A.
Im( i (A))

where S ⊆ Rn×n contains A.

With this general notation, we present a definition of system
aliasing only in terms of the A matrix in state-space representations and the sampling period h, which does not depend
on specific identification methods or data.

Definition 3 (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 = SA :=

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
(1)) and sampling frequencies. If the M matrix is specifically
constructed by data instead of I, E (A, M, h, S ) = {A},
where A denotes the ground truth, tells that the underlying
system is identifiable from the given data (see [22, Sec. IIIB]). Obviously if we have system aliasing for the system with
a specific sampling frequency, without extra prior information
on A, the system is always not identifiable.
IV. N O SYSTEM ALIASING : THE MINIMAL SAMPLING
FREQUENCY

Provided with the definition of system aliasing, a question
comes first: what (A, h) satisfies E (A, I, h, SA ) = {A}. To
answer this question, we need to introduce a theorem on matrix
logarithm.
Theorem 4 (principal logarithm [20, Thm. 1.31]). Let P ∈
Cn×n have no eigenvalues on R− . There is a unique logarithm
A of P all of whose eigenvalues lie in the strip {z : −π <
Im(z) < π}. We refer to A as the principal logarithm of P of
write A = Log(P ). If P is real then its principal logarithm
is real.
To make the principal matrix logarithm Log(·) be welldefined, we always assume that exp(hA) has no negative
real eigenvalues. Let G (h) = {z ∈ C : −π/h < Im(z) <
π/h, h ∈ R}. By Theorem 4 and 18, it always holds that
Log(exp(hA))/h ∈ E (A, I, h, SA ). To avoid system aliasing,
it implies that Log(exp(hA))/h = A , i.e. eig(A) ∈ G (h). It
is summarized as the following lemma.
Lemma 5. Let Ad = exp(hA), h ∈ R+ , which has no
negative real eigenvalues, and Â = Log(Ad )/h. Then A = Â
(i.e. E (A, I, h, SA ) = {A}) if and only if eig(A) ∈ G(h).

Given no other information on the system, consider the
identification problem of A using full-state measurement. It is
necessary to decrease the sampling period h until the ground
truth falls into the strip of G (h), and then the principal logarithm refers to the ground truth A, as illustrated in Figure 6.
Otherwise, we would be bothered by system aliases of A and

Im( i (A))h = ⇡

1
2
3

.. 0
..
..

1
8

1 h
max
4

1
2

1 h

p

Im( i (A))h =

⇡

Fig. 6: The imaginary parts of all eigenvalues of A must
lie into (−π/h, π/h). λi (·) denotes the i-th eigenvalue of
A in Theorem 18. 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.
Theorem 6 (Nyquist-Shannon-like sampling theorem). Consider the equidistant sampling. To uniquely reconstruct the
continuous-time system A from the corresponding discretetime system Ad by taking the principal matrix logarithm, the
sampling frequency ω (rad/s) must satisfy
ω ≥ 2 max {| Im (λi (A)) |, i = 1, . . . , n} .
Equivalently, the sampling period h (i.e. 2π/ω) should satisfy
h ≤ min {π/| Im (λi (A)) |, i = 1, . . . , n} .
Proof. The result immediately follows by verifying the condition eig(A) ∈ G(h) in Lemma 5.

The meaning of Theorem 6 in continuous-time system
identification can be understood by analogy with the NyquistShannon sampling theorem in signal processing. The NyquistShannon sampling theorem gives conditions on sampling frequencies, by looking at spectral information of signals, under
which the continuous signal can be uniquely reconstructed
from its discrete-time signal. As an analogy, Theorem 6
addresses that the continuous-time LTI system can be uniquely
reconstructed from the discrete-time system under a condition
that is built based on the spectrum of the A matrix.
Now we would like to show a property of matrix exponential
and logarithm, which further leads to a test criteria on system
aliasing. See Appendix B for the proofs of Lemma 7 and
Proposition 8.
Lemma 7. Considering h1 , h2 ∈ R+ and A ∈ Rn×n , let Â
be defined by Â = Log(exp(h1 A))/h1 . Then exp(h2 A) =
exp(h2 Â) if and only if h2 /h1 ∈ N or A = Â.

Proposition 8. Consider the dynamical system (1) without
inputs (i.e. B = 0), and two sampling periods h1 , h2 ∈ R+
such that h2 /h1 ∈
/ N. Let Â = Log(exp(h1 A))/h1 and Â 6=
A. The one-step prediction errors w.r.t. h2 are defined as
(tk ) = x(tk ) − exp(h2 A)x(tk−1 ),

where (tk , ϑ) := y(tk ) − ŷ(tk |tk−1 , ϑ), ŷ(tk |tk−1 , ϑ) denotes
the conditional mean of y(tk ), and Λ(tk , ϑ) the corresponding covariance matrix. The optimal prediction of y(tk ) (i.e.
ŷ(tk ), Λ(tk )) is obtained using Kalman filters (e.g. [15], [23]),

ˆ(tk ) = x(tk ) − exp(h2 Â)x(tk−1 ).
Assuming that E(x(tk−1 )) 6= 0, it yields
E((tk )) = 0,

E(ˆ
(tk )) 6= 0.

We have similar results for the case with inputs, as stated
in Proposition 22 in Appendix B, where we no longer require
E(x(tk )) 6= 0 due to the benefits from inputs. Meanwhile,
according to the condition (39), it is possible that a carefully
designed input signal invalidates the test criterion that is
built by evaluating E(ˆ
(tk )), which in practice may not be a
problem. The results in Proposition 8 and Proposition 22 can
be understood by Figure 7, where the output prediction of Â
(that is estimated from samples in h1 ) presents different values
from that of A in another sampling period h2 and it results in
that the expectation of one-step prediction errors is no longer
zero. The test criterion on system aliasing is summarized as
follows:

output

2
cont sys A^
cont sys A
dist sys Ad

1
0

h1

h2

-1
0

0.5

1

1.5

2

w(t), x(t0 ) are jointly Gaussian, the negative logarithmic
likelihood function is
PN
L(ϑ) = −2 log p(Y N |ϑ) = k=1 log det Λ(tk , ϑ)
PN
(8)
+ k=1 T (tk , ϑ)Λ−1 (tk , ϑ)(tk , ϑ) + const,

2.5

3

time (second)

Fig. 7: System responses of the dynamical system (1) with
A, Â, where Ad = exp(h1 A), Â = Log(Ad )/h1 , and the input
signal is the square wave with period 2h1 . The dots in deep
blue are samples with the sampling period h2 .
Test Criterion (system aliasing). Identify A, B by PEM or
ML (denoting the estimates by Â, B̂) assuming no system
aliasing under the sampling period h1 , i.e. Â asymptotically
converges to Log(exp(h1 A)). Choose another sample period
h2 such that h2 /h1 ∈
/ N+ , and sample by h2 the system
responses with non-zero initial conditions or non-zero inputs
(assuming (39) is satisfied). Use Â, B̂ to calculate the onestep prediction errors {(tk )}. Perform t-test to obtain the pvalue to make decisions, where the null hypothesis is that
{(tk )} comes from a normal distribution with mean zero and
unknown variance. Rejecting the null hypothesis implies the
existence of system aliasing.
V. N O SYSTEM ALIASING : SPARSE NETWORK
RECONSTRUCTION

Considering the systems given by (1) and (2), the likelihood
function is determined by the multiplication rule for conditional probability [23] p(Y N |ϑ) , p(y(tN ), . . . , y(t1 )|ϑ) =
p(y(tN )|y(tN −1 ), ϑ) p(y(tN −1 )|y(tN −2 ), ϑ) · · · p(y(t1 )|ϑ),
where ϑ denotes the parameters under estimation, which
parameterizes A, B, R, m0 , R0 . With the assumption that

ŷ(tk |tk−1 ) = C x̂(tk |tk−1 )
x̂(tk |tk ) = x̂(tk |tk−1 ) + K(tk )(tk )
d
x̂(t|tk ) = Ax̂(t|tk ) + Bu(t), tk ≤ t ≤ tk+1
dt
K(tk ) = P (tk |tk−1 )C T Λ−1 (tk )
P (tk |tk ) = P (tk |tk−1 ) − K(tk )CP (tk |tk−1 )
d
P (t|tk ) = AP (t|tk ) + P (t|tk )AT + R, tk ≤ t ≤ tk+1
dt
Λ(tk ) = CP (tk |tk−1 )C T ,
(9)
where the initial condition is x̂(t1 |t0 ) = m0 , P (t1 |t0 ) = R0 .
Considering the equidistant sampling and assuming the input
is constant over the sampling periods, the matrix Λ and K
appears in (9) can be treated as constant matrices by using
steady-state Kalman filtering [23, Sec. 3.6].
Now consider the full-state measurement case (i.e. C =
I) and restrict the noise to process noise. The calculation of
prediction ŷ(tk |tk−1 ) becomes particularly simple since K in
(9) always equals the identity, which yields
(tk , ϑ) = y(tk ) − Ad y(tk−1 ) − Bd u(tk−1 ),
Λ(tk , ϑ) = Rd ,

(10)

where Ad , Bd , Rd are defined via A, B, R in (4), (5). Here we
resolve p(y(t1 )|ϑ) by using p(y(t1 )|y(t0 ), ϑ), where y(t0 ) is
treated to be the deterministic and hence is removed from the
conditional variables, and takes the first sample as its value.
This simplification is due to the fact that K ≡ I, P (tk |tk ) ≡ 0
and the measurement of x0 is available (using y(t0 )), which
also leads to that the best estimate of the distribution of
x0 is nothing better than a delta function (even if including
the probability assumption of x0 in maximum likelihood,
i.e. p(y(t1 )|ϑ) = p(y(t1 )|y(t0 ), ϑ)p(y(t0 )|ϑ) and p(y(t0 )|ϑ)
takes the Gaussian density with mean m0 and covariance
matrix R0 ). Alternatively, the likelihood function (8) can be
obtained directly by considering (3) without using Kalman
filtering. However, the above standard procedure values when
we deal with general cases (i.e. C 6= I). Noticing the particular
parameterization [8, p. 92, 206], maximizing likelihood can be
performed as follows2 :
PN
θ̂ = argminθ k=1 (tk , θ)T (tk , θ),
(11a)
PN
1
T
(11b)
R̂d = N k=1 (tk , θ̂) (tk , θ̂),
where θ is composed of A, B. Note that this is not generally
true for maximizing likelihood, which should be solved by
[8, p. 219] in general for multivariable cases. To estimate
2 Similar to [8, p. 219], one instead firstly minimizes the cost function
analytically with respect to ϑ for every fixed R. Due to the particular
parameterization, the resultant optimization no longer depends on R.

θ, instead of minimizing the prediction error as (11a), we
impose the l1 -penalty to favor the sparse solution in network
reconstruction. This is due to the observations in Section III-A:
the consistency of ML may fail to present us with a correct
network structure unless appropriate thresholds of zero for
each row of A, B are selected, which is hardly implemented
in practice.
Remark 3. If we include measurement noise3 or consider
the output measurement case C 6= I, the prediction includes
the Kalman filter gain K, which depends on A, R and the
covariance of measurement noise. It deserves to be emphasized that, due to the possible large sampling periods h,
the numerical tricks used in [15], [23] are no longer valid
to compute the gradient of the prediction error. We have
to analytically calculate the gradient as far as possible until
the numerical computation is no longer restricted by h. This
problem becomes fairly complicated.
A. The cost function in matrix forms and the gradients
The reconstruction algorithm is supposed to infer a sparse
network, i.e. A is sparse. Due to the nonlinear least-square
cost function (11a), it no longer satisfies the setup of Sparse
Bayesian Learning proposed in [12]. Here we enhance sparsity
by heuristically imposing the l1 -norm of A as the penalty to
the PEM cost function as the first tentative treatment.
Considering the measurement signal Y N , let
X+ , [y(t1 ), y(t2 ), . . . , y(tN )],
X− , [y(t0 ), y(t1 ), . . . , y(tN −1 )],
U− , [u(t0 ), u(t1 ), . . . , u(tN −1 )],
where X+ , X− ∈ Rn×N . The matrix form of the l1 regularised PEM problem is formulated as
Z h
2
exp(sA) ds BU− F
minimize X+ − exp(hA)X− −
A

0

+λ kAk1 ,

(12)
where λ ∈ R+ and h ∈ R+ is the fixed
and
known
sampling
Pn
period, and the l1 norm kAk1 := i,j=1 |Aij | (Aij denotes
the (i, j)-th element of A). To avoid dealing with tensors, we
use the vectorized form of (12) as follows:
minimize
A

T
vec(X+ ) − (X−
⊗ In ) vec(exp(hA))
Z h
2
T T
−(U−
B ⊗ In )
vec(exp(sA)) ds 2

+λ kvec(A)k1 ,

(13)

0

2

T
where vec(X+ ) ∈ RnN , vec(exp(Ah)) ∈ Rn , (X−
⊗ In ) ∈
nN ×n2
, and ⊗ is the Kronecker product.
R
The problem (12) is challenging in optimization by noticing
that it is: non-convex due to matrix exponential; not globally
Lipschitz; and non-differentiable. The intuitive idea here is to
use the the Gauss-Newton framework, in which each iteration
3 Assume that it is reasonable to determine the network by A, B similarly,
even though we don’t have network models well-defined from the state-space
representations with measurement noise.

is to solve a constrained l1 -regularized linear least square
problem. Let
T
r(A, B) , vec(X+ ) − (X−
⊗ I ) vec(exp(hA))
Z hn
T T
vec(exp(sA)) ds,
−(U−
B ⊗ In )

(14)

0

φ(A, B) := r(A, B) r(A, B), and f (A, B) , φ(A, B) +
λk vec(A))k1 , which is the objective function of (13).
Then minA φ(A, B) denotes the problem (13) without l1 penalisation. The gradient of φ(A, B) w.r.t. A is ∇A φ =
T
2JA
(A, B)r(A, B), where
Z h
T
T T
sK(sA) ds,
JA = −h(X−
⊗ In )K(hA) − (U−
B ⊗ In )
T

0

(15)
and K(A) is defined in Theorem 21. The function sK(sA), s ∈
[0, h] is integrable by noticing ksK(sA)k ≤ h exp(hkAk)
(k · k denotes any matrix norm). The matrix function K(hA)
and the integration can be calculated numerically given A. To
compute the gradient of φ(A, B) w.r.t. B, using the matrix
identity vec(AXB) = (B T ⊗ A) vec(X)
 R h again for the term
T T
of U in (12), it yields U− B ⊗ In 0 vec(exp(sA)) ds =

Rh
T
U−
⊗ 0 vec(exp(sA)) ds vec(B). Then it follows that
T
∇B φ = 2JB
(A, B)r(A, B), where
Z h
T
JB = U−
⊗
vec(exp(sA)) ds.
(16)
0

If we assume B is diagonal and ∇B φ is calculated w.r.t. each
T
diagonal element of B, then ∇B φ = 2 (I ⊗ 1)JB r(A, B),
2
where I is an n×n identity matrix and 1 is an n -dimensional
row vector of 1’s. In a sum, the gradient of φ(A, B) is


 T
∇A φ
J
T
∇φ =
=2 A
T r(A, B) , 2J (A, B)r(A, B).
∇B φ
JB
(17)
B. A special case: update A with fixed B
The subspace method in system identification presents us
with nice initial estimation of A0 , B0 (e.g. see [24]). Concerning the task of network reconstruction, we would like to infer
a sparse A from data. As a special case, we only update A by
solving (12) with B fixed to be B0 . For simplicity, in this subsection, let r(A, B0 ) , r(A), JA (A, B0 ) , J(A), φ(A, B0 ) ,
φ(A) and f (A, B0 ) , f (A).
A linear approximation of r(A) in a neighbourhood of a
given point Ac is r(Ac ) + J(Ac ) vec(A − Ac ). One may then
use this approximation and formulate a l1 -regularized linear
least squares problem
minimize kr(Ac ) + J(Ac ) vec(A − Ac )k22 + λk vec(A)k1 ,
A
(18)
which can be solved to obtain an approximate solution to (13).
Resolving it in an iterative way amounts to a Gauss-Newton
method. However, vec(A−Ac ) is not necessary to be a descent
direction of (13).
In the k-th iteration, to guarantee the step pk = vec(A−Ak )
being a descent direction of (13), the search direction pk is

instead computed from the following constrained optimization Solving (12) is to search a sparse A in the neighborhood of A0 .
problem
Moreover, we have the following propositions to guarantee fair
(
properties of this algorithm. See Appendix C for the proofs.
2
minimize
kr(Ak ) + J(Ak )pk k2 + λk vec(Ak ) + pk k1 ,
n2
p
∈R
k
(P1 )
Proposition 9. Let fˆ(Ak , pk ) denote the objective function
subject to supg∈∂f (Ak ) g T pk ≤ 0,
of (P1 ) and p∗k be its optimal point. If p∗k 6= 0 and
where ∂f (Ak ) denotes the subdifferential of f (A) at Ak , supg∈∂f (Ak ) g T p∗k = 0, then 0 ∈ ∂f (Ak ).
defined as
Proposition 10. Let fˆ(Ak , pk ) and p∗k be defined in Propo(19) sition 9. If p∗k = 0 and 0 ∈ argminp fˆ(Ak , pk ), then
∂f (Ak ) := {∇φ(Ak ) + λz : z ∈ J1 × · · · Jn2 },
k

0 ∈ ∂f (Ak ).

[−1, 1] if vec(Ak )i = 0
Proposition 9 guarantees that the step p∗k from solvJi := {1}
(20)
if vec(Ak )i > 0 ,


ing (P1 ) will always be a descent direction of (13) (i.e.
{−1}
if vec(Ak )i < 0
supg∈∂f (Ak ) g T p∗k < 0) until either it reaches the stationary
in which vec(Ak )i denotes the i-th element of vec(Ak ). One point or {p∗k } converges to zero. When {p∗k } approaches
may have noticed that the constraint in the problem (P1 ) to zero, there are two cases: one is Proposition 10 which
is the definition of descent direction for f at Ak , except guarantees that it reaches the stationary point; the other is
replacing < 0 with ≤ − to guarantee the existence of described as follows:
minimum. The problem (P1 ) is a convex optimization problem
p∗k = 0 is the unique optimal point of (P1 ) and
by noticing that supg∈∂f (Ak ) g T pk is a convex function, which
0∈
/ argminpk fˆ(Ak , pk ).
is a pointwise supremum over an infinite set of a linear function
[25, chap. 3]. To solve the problem (P1 ), we need to explore Regarding the second case, indeed, if there exists other optimal
∗∗
∗∗
the constraint and derive an equivalent form (see Appendix C point pk 6= 0 of (P1 ), we instead consider pk using
∗
Proposition 9. This is why we restrict pk = 0 to be the unique
for details), given as follows.
optimum
of (P1 ). In the second case, {p∗k } converges to zero
(
2
minimize
kr(A
)
+
J(A
)p
k
+
λk
vec(A
)
+
p
k
,
k
k
k
k
k
1
2
and
the
objective
value f (Ak ) also converges. However, in
pk ∈Rn2
(P10 )
T
theory,
we
fail
to
prove
that the limit point of {Ak +ivec(p∗k )}
subject to ḡ(Ak ) pk + λkW (Ak )pk k1 ≤ 0,
is a stationary point. In the sense of applications, it has been
where
fairly good since we are looking up a sparse solution in the
ḡ(Ak ) = ∇φ(Ak ) + λ sgn(Ak ),
(21) neighborhood of a fairly good estimate.
W (Ak ) = I − diag(| sgn(Ak )|),

the identify matrix I is of a compatible dimension, | ·| denotes
the element-wise absolute value, diag(v) denotes the diagonal
matrix built from vector v, and the sgn function for vectors
and matrices is extended from the standard signum function
for real numbers, defined as follows: when x ∈ Rn , sgn(x)
denotes a n-dimensional vector whose i-th element equals
sgn(xi ); and when X ∈ Rm×n , sgn(X) := sgn(vec(X)).
Now the problem (P 0 ) can be easily modeled using CVX in
MATLAB and solved by standard optimization solvers [26].
The iterate is updated via
vec(Ak+1 ) = vec(Ak ) + sk pk ,

(22)

where the step length sk is determined by backtracking line
search. Let f 0 (Ak ; pk ) denote the directional derivative of f
at Ak in the direction pk , by subcalculus,
f 0 (Ak ; pk ) := supg∈∂f (Ak ) g T pk
= ḡ(Ak )T pk + λkW (Ak )pk k1 .

(23)

Given α ∈ (0, 0.5), β ∈ (0, 1) and an initial value sk = 1, the
line search is to perform sk ← βsk until

f Ak + ivec(sk pk ) ≤ f (Ak ) + αsk f 0 (Ak ; pk ).
(24)

The whole iterative method for (13) is summarized in Algorithm 1. One has to note that this algorithm cannot guarantee
that the iterate will converges to the stationary point. It is lucky
that we have good initial values of A0 , B0 to start with that
is provided by the subspace method in system identification.

Algorithm 1 Modified Gauss-Newton for l1 -regularized nonlinear
least square problems
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:

given  > 0, tolerance δ > 0; α ∈ (0, 0.5), β ∈ (0, 1)
initialize A0 , B0 by the subspace method.
repeat
Calculate r(Ak ), J(Ak ), ∇φ(Ak ) using (14), (15).
Calculate ḡ(Ak ), W (Ak ) using (21).
Compute the search direction pk by solving problem (P 0 ).
Compute the directional derivative f 0 (Ak ; pk ) using (23).
Determine the step size sk by backtracking line search (24).
Update the iterate Ak by (22).
until kAk − Ak+1 k2 < δ

Remark 4. The proposed method can be considered as a
variant of the damped Gauss-Newton method. If the Jacobian
matrix J(Ak ) does not have full column rank, one could adopt
the Levenberg-Marquardt method and solve
minimize
kr(Ak ) + J(Ak )pk k22 + µk kpk k22 + λk vec(Ak ) + pk k1 ,
2
pk ∈Rn

subject to supg∈∂f (Ak ) g T pk ≤ −.

C. General cases: update both A and B
Considering the vectorized form (13), let θ denote the
optimal variables, i.e. θ , [vec(A)T vec(diag(B))T ]T , where
diag(B) denotes the diagonal elements of B. The other notations follow that r(A, B) , r(θ), J(A, B) , J(θ), φ(A, B) ,
φ(θ) and f (A, B) , f (θ). In the same way as Section V-B,

output

the approximated l1 -regularized linear least square problem is a proper (Z T ⊗Z −1 )T (Z T ⊗Z −1 )-weighted vector norm in
with constraints is written as
terms of vec(Â). Using khZ (·)kF is on the one hand simplify(
2
minimize
kr(θk ) + J(θk )pk k2 + λk vec(Ak ) + Λpk k1 , ing the analysis we conduct throughout this section, and on the
pk ∈R2n2
other explicitly penalizes the imaginary part of the eigenvalues
(P2 )
subject to supg∈∂f (θk ) g T pk ≤ 0,
without “distorting” them through the transformation by Z.
Now we define Sκ using the norm khZ (·)kF . The basic
2
2
where Λ = [I 0] is of dimension n ×2n , ∂f (θk ) denotes
idea is that one should exclude such A’s whose imaginary
the subdifferential of f (θ) at θk , defined as
parts of eigenvalues are too large, which implies their system
response will show wild fluctuation, as illustrated in Figure. 8.
∂f (θk ) := {∇φ(θk ) + λΛT z : z ∈ J1 × · · · Jn2 },
That’s why we need to consider a reasonable set Sκ rather
and Ji , i = 1, . . . , n2 are defined in (20). Equivalently, we
solve the following convex-constrained convex problem to
6
4
cont sys A_1
update A, B by θk+1 = θk + pk .
cont sys A_2
2
(
cont sys A_3
2
minimize
kr(θ
)
+
J(θ
)p
k
+
λk
vec(A
)
+
Λp
k
,
dist sys A_d
0
k
k k 2
k
k 1
0
2n2
p
∈R
k
(P2 )
-2
0
0.5
1
1.5
2
2.5
3
subject to ḡ T pk + λkW (Ak )Λpk k1 ≤ 0,
time (second), h = 1s
where ḡ , ∇φ(θ) + λΛT sgn(Ak ) and W (Ak ) = I −
diag(| sgn(Ak )|). The backtracking line search is equipped
in the same way as Section V-B and the algorithm trivially
follows by modifying Algorithm 1.

Fig. 8: An example of system responses of multiple system
aliases, i.e. Ai (i = 1, 2, 3) that satisfies exp(hAi ) = Ad , where
A2 is the ground truth but A1 is the principle logarithm, and
the input signal is a square wave of period 2h.

VI. S YSTEM ALIASING AND BOUNDED CONSTRAINTS

than Rn×n in (26). In practice, even if we know the sampling
frequency is not high enough to guarantee no system aliasing,
we should still believe that the measurements do not miss
too many fluctuation between samples. To make the constraint
in (25) practically meaningful, we restrict S to be a norm
bounded subset

Sκ = Ã ∈ Rn×n : khZ (Ã)kF ≤ κ .
(27)

In the previous section we hinted that the conditions for
no system aliasing follow as a consequence of bounded
eigenvalues. In this section we follow this path and study the
problem in the presence of system aliases.
Consider the case of system aliasing, i.e. h is NOT chosen
small enough such that E (A, I, h, SA ) = {A}. In order to
find out A among the aliases we need extra information,
for instance, the properties of A known a priori. Here we
assume that the ground truth A is the sparest solution in
E (A, I, h, Sκ ) and κ ∈ R as an upper bound that has
been prescribed. The set Sκ will be defined after giving
Definition 11. A can be searched by the criterion
(25)

minimize kÂk0 .

Â∈E (A,I,h,Sκ )

Here is a niche that is the calculation of E (A, I, h, Sκ ) from
data. By definition,

E (A, I, h, Sκ ) = Ã ∈ Sκ : exp(hÃ) = Ad ,
(26)

where Ad = exp(hA). Even if we know Ad has consistent
estimation via PEM or ML, considering the observations in
Section III-A, we know the workflow, that is estimating Ad
and then obtaining E by matrix logarithms, is not robust in
the presence of noise. In this section, we focus on studying
the possibility of searching A in the set of system aliases
E (A, I, h, S (κ)) using the prior information.

Definition 11 (Z-weighted norm). Let hZ (A) = Z AZ,
where Z is the matrix defined in Theorem 19. Then the norm
is defined as khZ (·)kF = k · kF ◦ hZ .
−1

To formulate Sκ , we introduce this special norm of A,
which is equivalent to the Frobenius norm up to a change of
coordinates. The matrix Z is constant, which can be obtained
by Jordan decomposition of Ad . One can observe that
khZ (Â)kF = vec(Â)T (Z T ⊗ Z −1 )T (Z T ⊗ Z −1 )vec(Â)

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

(28)

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) .

(29)

Moreover, let A0 denote a special matrix logarithm for which
all jk (k = 1, . . . , p) in (33) are equal to 0.
Definition 12 (equivalence relations). Let S denote the set of
all primary matrix logarithms

(30)
S := Ã ∈ Rn×n : exp(hÃ) = Ad .

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.
4 We use a different font type from j to avoid misunderstanding j as a scalar
variable for indexes.

Lemma 13. Let S be the set defined in (26) and parametrized
by (34) in Theorem 19. For any A1 , A2 ∈ S, khZ (A1 )kF =
khZ (A2 )kF if and only if A1 ∼ A2 .

Lemma 14. Given any Ā ∈ S, there exists a finite number of
Ai ∈ S that satisfies Ai ∼ Ā.

Lemma 15. There exists a finite number of Ai ∈ S such that
khZ (Ai )kF ≤ κ.

Proposition 16 (lower boundness of logarithms). Let S be the
set defined in (26). Given any Ā ∈ S, there exists M (Ā) > 0,
such that for any A ∈ {A ∈ S : A  Ā}, it holds that

contaminated by process noise and unable to give reasonable
sparse structure of A, clearly shown in Figure 10a. This is
mainly due to the effect of process noise on the estimation
through matrix logarithms. However, the direct logarithm of
Âd might also work well when the dimension is small (e.g.
dim(A) ≤ 6).
2

5

5

10

10

10

15

15

15

20

20

20

-1
-2

Proof. It immediately follows by choosing

VII. N UMERICAL EXAMPLES

1.5

0

0
-0.5

-1

-1
0

20

40

60

0.2

1

20

40

60

1

20

40

60

0

0
0

-1

-0.5

-0.5

-1
0

20

40

60

20

40

60
2

-0.5

0

1

-1

-2

0

-1.5

-4

-1

-2

-6

-2

0

20

40

60

0

20

40

60

2

0.1

1

1

0

40

60

20

40

60

0.5

1

0

0.5

-0.5

20

40

60

0

60

0

20

40

1.5

20

40

60
1

0

0.5

-1

0

-2

-0.5

-3

-1
0

20

40

60

0.5

1

0

-1

-0.1

-1
0

20

40

60

20

40

60

20

40

60

40

60

20

A^d
10

5

5

5

10

10

10

15

15

15

20

20

20

5
0

10

15

20

5

10

15

20

5

10

15

20

Fig. 10: Matrix A and Ad visualized in scaled colours.
A and Ad are the ground truth; Âalg is estimated by the
proposed method using λ = 0.04, and Âd,alg is calculated by
exp(hÂalg ); Âlogm is the principal logarithm of Âd divided
T
T −1
by h, and Âd = X+ X−
(X− X−
) .
The convergence is illustrated in Figure 11, with respect to
λ = 0.04 using Algorithm 1. Here the λ is chosen by performing network reconstruction using λ logarithmically ranging
from 10−4 to 100 and checking the sparsity of Â’s (users’ prior
knowledge) and the resultant prediction errors (whiteness,
mean, standard deviations). A better λ could be automatically
calculated by running the cross-validation technique, when the
amount of data allows. As performance visualization widely
used in bioinformatics, the Receiver Operating Characteristic
(ROC) curve and the Precision-Recall (P-R) curve are provided
in Figure 12.

0

20

40

60

0

20

40

60

10 0
pred. error
pred. error of groundtruth
pred. error of logm

10 -5

10 -10

-1

-1
20

15

10 5

0

0

10

60

1

-0.5

-0.5
0

40

2

0

0

-0.2
0

20

-0.5
0

1

60

5

0

-1
0

20

1.5

0.5
0

0

1

-1
40

40

0.5

1

20

20

1

60

15

-0.1
0

1.5

2

0

2

20

-0.5
0

2

-0.5
0

0

-2
0

0

0.1
0

-0.5
0

1

0.5

0.5

0

0

-2
0

0.3

0.5

1
0.5

10

A^d;alg

Ad (ground truth)

^ 2 kF
kX1 ! exp(hA)X

0

5

(b) The corresponding Ad ’s via exp(hA)

In this section, we use an numerical example to show
the effectiveness of the proposed algorithm. The data set for
identification was sampled with h = 2 from simulation of
stochastic differential equations, in which A is a stable sparse
matrix generated randomly. The initial values of states were
randomly generated, and the process noise is normal i.i.d. Here
we chose a low sampling frequency, large noise and limited
samples to generate a challenging time series for identification,
which gives a typical profile of time series data in biomedicine
(e.g. microarray data [14]), as illustrated in Figure 9.
1

20

-10

where M (Ā) is the lower bound on the gap between Ā and
any A  Ā ∈ S, defined in Theorem 16.

1

15

(a) The ground truth and the estimates of A; λ = 0.04

5

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

2

10

-5

κl > max{0, khZ (Ā)kF − M (Ā)},

1

1
0

5

Proposition 17. Let S be the set defined in (26). For any
Ā ∈ S, there exist κl , κu ∈ R in S (κl , κu ) = {Ã ∈ Rn×n :
κl ≤ khZ (Ã)kF ≤ κu } such that (25) has a unique optimal
point in the sense of the equivalence relation in Definition 12.

0.5

A^logm

5

khZ (A)kF − khZ (Ā)kF ≥ M.

1.5

A^alg

A (ground truth)

-2
0

20

40

60

0

20

40

60

Fig. 9: System response for sparse identification (24 outputs,
24 samples, random initial states).
The results are shown in Figure 10, together with their corresponding Ad ’s computed via matrix exponential. Note that
in order to estimate A, we have to use the matrix exponential
and matrix logarithm due to the large sampling period. The
simplest way to estimate A is taking the principal matrix
logarithm of least square solution Âd , which is, however,

10 -15

0

5

10

15

iteration

Fig. 11: Convergence with respect to λ = 0.01.
The numerical example was programmed and computed in
MATLAB, and the codes will be released in public on the
github.com/oracleyue. Considering computational efficiency,
we directly used vector/matrix norms instead of quadratic
forms in implementation (cf. [26, chap. 11.1]).

ROC Curve for DNR (AUC: 0.82977)

P-R Curve for DNR (AUC: 0.58541)

1

0.8

0.8

Precision

True positive rate

1

0.6
0.4
0.2
0
0.2

0.4

0.6

regularized
principal logarithm

0.4
0.2

regularized
principal logarithm

0

0.6

0.8

1

0
0

0.2

0.4

False positive rate

0.6

0.8

1

Recall

(a) ROC curve

(b) P-R curve

Fig. 12: The ROC curve and Precision-Recall curve of the
proposed algorithm.

VIII. C ONCLUSIONS
Continuous-time system identification is challenging when
only low-sampling-frequency data with limited lengths are
available. Unfortunately, this is a typical profile of time
series in biomedicine applications. To reconstruct the correct
dynamic networks from such time series, we have to identify
the continuous-time models with sparse network structures.
This paper studies the full-state measurement case, which
is supposed to be the basic case while it shows particular
complications. We first clarify the concept of system aliasing,
which is raised by low sampling frequencies. A theorem on
how to choose the sampling frequency to guarantee no system
aliasing is provided, together with a test criterion. In regard
to the “easy” case, i.e. no system aliasing, we present an
algorithm to reconstruct sparse dynamic network from fullstate measurements. In the case with system aliasing, the
possibility on searching among system aliases is manifested
in theory relying on the prior information on network sparsity.
The paper dedicates to show the challenges and attract more
people contributing to this study.

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



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

..
.

f 0 (λk )
f (λk )







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 19 (classification of logarithms [20, Thm. 1.28]).
Let the nonsingular matrix P ∈ Cn×n have the Jordan
canonical form (31) 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)
Aj = Z diag(Lj11 , L2 2 , ..., L(j
,
p )Z

(34)

Ljkk

where
is defined in (33), corresponding to all possible
choices of the integers j1 , ..., jp , subject to the constraint that
ji = jk whenever λi = λk .
If s < p then eA = P has nonprimary solutions. They form
parametrized families
(j )

−1 −1
p)
Aj (U ) = ZU diag(Lj11 , L2 2 , ..., L(j
Z ,
p )U

(35)

where jk is an arbitrary integer, U is an arbitrary nonsingular
matrix that commutes with J, and for each j there exist i and
k, depending on j, such that λi = λk while ji 6= jk .

Definition 20 (Fréchet Derivatives [20]). The Fréchet derivative of the matrix function f : Cn×n → Cn×n at a point
X ∈ C(n×n) is a linear mapping
L

Cn×n −→ Cn×n

E 7−→ L(A, E)

such that for all E ∈ C

A PPENDIX A

M ATRIX E XPONENTIAL AND L OGARITHM
Theorem 18 (Gantmacher [20, Thm. 1.27]). Let P ∈ C
be nonsingular with the Jordan canonical form
Z −1 P Z = J = diag(J1 , J2 , ..., Jp )

λk 1

.

λk . .
Jk = Jk (λk ) = 

..

.



n×n

(31a)

The Fréchet derivative is unique if it exists, and for matrix
functions exp (matrix exponential) and Log (principal matrix
logarithm) it exists. The Fréchet derivative of the function exp
[20] is
Z
eX(1−s) EeXs ds,

Lexp (X, E) =

(36)

0

which can be efficiently calculated by the Scaling-PadeSquaring method in [27]. It gives a linear approximation of
exp at a given point Ac in the direction E
ehA = eh(Ac +E) = ehAc + L(hAc , hE) + O(khEk2 ). (37)

Then all solutions to e = P are given by
A

A = ZU

f (A + E) − f (A) − L(A, E) = o(kEk).

1

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

diag(Lj11 , Lj22 , ..., Ljpp )U −1 Z −1 ,

n×n

(32)

The Fréchet derivative of the function Log [20] is
Z 1
LLog (X, E) =
(t(X − I) + I)−1 E(t(X − I) + I)−1 )dt,
0

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

(33)

and its efficient computation algorithm is provided in [28].

Theorem 21 (Kronecker representation [20, Thm. 10.13]). For
A ∈ Cn×n , vec(L(A, E)) = K(A) vec(E), where K(A) ∈
2
2
Cn ×n has the representations


(I ⊗ eA )ψ AT ⊕ (−A)







1 T
 AT /2
[A ⊕ (−A)]
(e
⊗ eA/2 ) sinch
K(A) =
2





T
1
1

A
A
T

[e ⊕ (−A)]
 (e ⊕ e )τ
2
2

where ψ(x) = (e − 1)/x and τ (x) = tanh(x)/x. The
third expression is valid if 21 kAT ⊕ (−A)k < π/2 for some
consistent matrix norm.
x

Here the operator ⊕ is the Kronecker sum, defined as
(A ⊕ B) = A ⊗ In + Im ⊗ B,

for A ∈ Rm×m and B ∈ Rn×n .

A PPENDIX B
T EST CRITERIA ON SYSTEM ALIASING
A. Test criteria for the cases with inputs
Proposition 22. Consider the dynamical system (1), and
two sampling periods h1 , h2 ∈ R+ such that h2 /h1 ∈
/
N.
Let
f
(A,
B,
h,
t
)
,
x(t
)
−
exp(hA)x(t
)
−
k
k
k−1
Rh
exp(sA)B ds u(tk−1 ), Â = Log(exp(h1 A))/h1 and B̂ be
0
a value such that E(f (Â, B̂, h1 , tk )|x(tk−1 ), u(tk−1 )) = 0 for
all tk . The one-step prediction errors w.r.t. h2 are defined as
(tk ) = f (A, B, h2 , tk ), ˆ(tk ) = f (Â, B̂, h1 , tk ). If


E exp(h2 A) − exp(h2 Â) x(tk ) +
Z h2

(39)

exp(sA)B − exp(sÂ)B̂ ds u(tk ) 6= 0,
0

we have

E((tk )) = 0,

E(ˆ
(tk )) 6= 0.

The proof follows trivially by evaluating the expectation
E(ˆ
(tk )). We no longer require E(x(tk )) 6= 0 since we can
take advantages of u(tk ) to satisfy (39). In the cases with
E(x(tk )) = 0 and non-zero inputs in expectation E(u(tk )) 6=
0, assuming A, Â are non-singular, the condition (39) can be
further simplified as

−1

exp(h2 Â) − I exp(h1 Â) − I
exp(h1 A) − I
(40)
6= exp(h2 A) − I.
Rh
The simplification follows from (39) by noticing 0 exs ds =
(ehx − 1)/x and A commutes with its matrix functions.
B. Proof Lemma 7
Proof. Let Ad := exp(h1 A), which has the Jordan canonical
form (31) (i.e. let P in Theorem 18) be Ad ). By Theorem 18
and 19, we have
h1 A = Z
h1 Â = Z

j
diag(Lj11 , . . . , Lpp )Z −1 ,
0
diag(L1 , . . . , L0p )Z −1 .

To compare exp(h2 Â) with exp(h2 A), we need to find their
Jordan canonical form by the definition of matrix exponential.

To calculate the eigenvalues of h2 Â, consider the determinant
|µ̂I − h2 Â| = 0 ⇔ |µ̂0 I − h1 Â| = 0, where µ̂0 , h1 /h2 µ̂ and
|µ̂0 I − h1 Â| = |µ̂0 I − diag(L01 , . . . , L0p )|. It is equivalent to
solve p equations |µ̂0k Ik − L0k | = 0 (k = 1, . . . , p), where Ik
denotes the identity matrix of the dimension compatible with
L0k . It yields that
|µ̂0k Ik − L0k
| = |µ̂0k Ik − log(Jk (λk ))|
log(λk ) log0 (λk ) · · ·

..

.
log(λk )
= µ̂0k Ik − 

.
..

= (µ̂0k − log(λk ))

mk

= 0,

∗
..
.







0
log (λk )
log(λk )

where Jk , λk , jk , mk , Imk are given in (31), (33); and hence
µk = h2 /h1 log(λk ) with geometric multiplicity mk . Similarly for h2 A, consider |µ0k Ik − Ljkk | = |µ0k Ik − log(Jk (λk )) −
2jk πiImk | = 0, where µ0k , h1 /h2 µk , µ is the eigenvalues of
h2 A, the integer jk is given in (33), and i is the imaginary unit.
It yields mk = h2 /h1 (log(λk ) + 2jk πi) with multiplicity mk .
Considering the special forms of Ljkk , we have the following
Jordan decomposition
diag(L01 , . . . , L0p ) = U diag(J1 (µ̂1 ), . . . , Jp (µ̂p ))U −1 ,
j

diag(Lj11 , . . . , Lpp ) = U diag(J1 (µ1 ), . . . , Jp (µp ))U −1 ,
where Jk (µ̂k ) and Jk (µk ) denote the corresponding Jordan
blocks. Therefore, exp(h2 A) = exp(h2 Â) is equivalent to
exp(Jk (µk )) = exp(Jk (µ̂k )) for any k = 1, . . . , p, which
implies exp(h2 /h1 2jk πi) = 0. It leads to the conditions:
h2 /h1 ∈ N+ , or jk ≡ 0, ∀k (i.e. Â = A).
C. Proof of Proposition 8
Proof. Considering the dynamical system (1), it is obvious that
E((tk )|x(tk−1 )) = 0. Now we evaluate the other expectation


E ˆ(t
k)
k ) = E x(tk+1 ) − exp(h2 Â)x(t

= E x(tk+1 ) − exp(h2 A)x(tk ) + 

exp(h2 A) − exp(h2 Â) x(tk )
= 0 + (exp(h2 A) − exp(h2 Â))E(x(tk )) 6= 0,

by Lemma 7.

A PPENDIX C
M ORE DETAILS ON THE OPTIMIZATION PROBLEMS
A. Equivalent forms of (P1 ), (P2 )
The following shows how to derive (P10 ) from (P1 ). Consider g ∈ ∂f (Ak ) in (P1 ), which implies g = ∇φ(Ak )+λz =
∇φ(Ak )+λ sgn(Ak )+λ(z −sgn(Ak )) , ḡ +λ(z −sgn(Ak )).
Recall that maxs∈[−1,1]n sT x = kxk1 where x is an ndimensional vector. Then we have
sup

g T pk = ḡ T pk + λ

g∈∂f (Ak )

= ḡ T pk + λ

T

sup
z∈J1 ×···×Jn2
T

sup

(z − sgn(Ak )) pk

z̄ (W (Ak )pk )

z̄∈[−1,1]l

= ḡ T pk + λkW (Ak )pk k1 ,

where W (Ak ) = I − diag(| sgn(Ak )|) and l = n2 −
card(sgn(Ak )).

Consider the optimization (P2 ), we go through the same
procedure and obtain

where



[−1, 1]
Jˆi = {1}


{−1}

sup g T pk = ḡ T pk + λkW (Ak )Λpk k1 ,

g∈∂f (θ)

where ḡ , ∇φ(θ) + λΛT sgn(Ak ).
B. Proof of Proposition 9
Proof. Without loss of generality, suppose that there exists p∗k 6= 0 and 0 ∈
/ arg minpk fˆ(Ak , pk ) such that
supg∈∂f (Ak ) g T p∗k = 0 and 0 ∈
/ ∂f (Ak ). (Indeed, if 0 ∈
ˆ
arg minpk f (Ak , pk ), we apply Proposition 10 and obtain 0 ∈
∂f (Ak ).) It implies that supg∈∂f (Ak ) g T αp∗k = 0, ∀α ∈ (0, 1].
Hence, f (Ak ) ≤ f (Ak + αp∗k ) for small enough α. Let
φ(Ak ), φ̂(Ak , pk ) be defined by

if vec(Ak )i + (pk )i = 0
if vec(Ak )i + (pk )i > 0 ,
if vec(Ak )i + (pk )i < 0

(pk )i denotes the i-th element of pk , and ∇φ̂(Ak , pk ) =
2J(Ak )T (r(Ak ) + J(Ak )pk ). Hence ∂p fˆ(Ak , 0) = ∂f (Ak ).
Therefore, 0 ∈ ∂p fˆ(Ak , p∗k ) with p∗k = 0 implies 0 ∈
∂f (Ak ).
A PPENDIX D
P ROOFS FOR BOUNDNESS OF SYSTEM ALIASES
A. Proof of Lemma 13
j

j (i)

(i)

Proof. Let Ai := Z diag(L11 , · · · , Lpp )Z −1 , where i =
(i)
jk

f (Ak ) = φ(Ak ) + λk vec(Ak )k1 ,
fˆ(Ak , pk ) = φ̂(Ak , pk ) + λk vec(Ak ) + pk k1 ,

(i)
1, 2, Lk := log(Jk (λk ))+2jk πiImk , and all other notations
are given in (34). By using (42) for A1 , A2 , we obtain

and hence φ(Ak + ivec(pk )) = φ̂(Ak , pk ) + o(kpk k) (k · k denotes any vector norm). For simplicity, without any ambiguity,
we use f (Ak + pk ) to represent f (Ak + ivec(pk )). Then we
have f (Ak + pk ) = fˆ(Ak , pk ) + o(kpk k), which yields

khZ (A1 )kF = khZ (A2 )kF ⇔ khZ (A1 )k2F − khZ (A0 )k2F

lim
α↓0

|f (Ak +

− fˆ(Ak , αp∗k )|
= 0.
αkp∗k k

αp∗k )

(41)

Now let us calculate this limit in a different way. Noting that
fˆ(Ak , p∗k ) ≤ fˆ(Ak , αp∗k ) (since p∗k ∈ arg minpk fˆ(Ak , pk ))
and fˆ(Ak , p∗k ) < fˆ(Ak , 0) = f (Ak ) (since 0 ∈
/
arg minpk fˆ(Ak , pk )), we have
fˆ(Ak , αp∗k ) − f (Ak ) ≥ fˆ(Ak , p∗K ) − f (Ak ) ≥ δ > 0.
Moreover, since f (Ak + αp∗k ) − f (Ak ) = φ(Ak + αp∗k ) −
φ(Ak ) + λ(kAk + αp∗k k1 − kAk k1 ) = o(αkp∗k k) + O(αkp∗k k),
there exists M ∈ R such that
lim
α↓0

|f (Ak + αp∗k ) − f (Ak )|
≤ M.
αkp∗k k

Now we recalculate the limit in (41) as follows
lim
α↓0

|f (Ak + αp∗k ) − fˆ(Ak , αp∗k )|
αkp∗k k

|f (Ak + αp∗k ) − f (Ak )|
≥ lim
−
α↓0
αkp∗k k
δ
≥ M − lim
= +∞,
α↓0 αkp∗
kk

f (Ak ) − fˆ(Ak , αp∗k )
αkp∗k k

which contradicts with (41).
C. Proof of Proposition 10

Proof. Since p∗k ∈ argminpk fˆ(Ak , pk ), it yields 0 ∈
∂p fˆ(Ak , p∗k ). Note that, in our discussion, Ak is always fixed,
and thus ∂p fˆ(Ak , pk ) denotes the subgradient of fˆ(Ak , ·) at
pk . Now let us write ∂p fˆ(Ak , pk ) explicitly
n
o
∂p fˆ(Ak , pk ) = ∇φ̂(Ak , pk ) + λẑ : ẑ ∈ Jˆ1 × · · · Jˆn2 ,

= khZ (A2 )k2F − khZ (A0 )k2F ⇔ I (j(1) , j(2) − j(1) ) = 0,

which implies that A1 ∼ A2 by definition. The first equality
in (42) is due to the linear transformation hZ (Â).
B. Proof of Lemma 14
Proof. Let j denote [j1 , . . . , jp ] of Ā in (33), 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 12, 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)
|δi + ji + βi /2| ≤
(43)
mi
for all i = 1, . . . , p. Otherwise, supposing that there exists
i ∈ {1, . . . , p} such that δi does not satisfy (43), we will have
X
I (j, δ) = mi (δi + ji + βi /2)2 +
mk (δk + jk βk /2)2
−

X

k6=i

2

mk (jk + βk /2) >

k

X
k6=i

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

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

C. Proof of Lemma 15
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| ≤
(44)
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)
(i)
Pp
jk ∗ jk
(0)∗ (0)
= k=1 tr Lk Lk − Lk Lk


Pp
(i)
(i)2
= k=1 tr 2jk πi (log(Jk )∗ − log(Jk )) + 4π 2 jk Imk
Pp
(i)
(i)
(i)
(i)
= k=1 4πjk mk (βk + jk ) = 4π I (0, j), j , [j1 , . . . , jp ].
for all i = 1, . . . , p. Otherwise, by supposing that there exists
i ∈ 1, . . . , p such that δi does not satisfy (44)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 (44) is finite, which
finalizes the proof.
D. Proof of Proposition 16
Proof. Let j denote [j1 , . . . , jp ] of Ā in (33), 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 15.
R EFERENCES
[1] B. Ø. Palsson, Systems biology: simulation of dynamic network states.
Cambridge University Press, 2011.
[2] Z. Bar-Joseph, A. Gitter, and I. Simon, “Studying and modelling
dynamic biological processes using time-series gene expression data,”
Nature Reviews Genetics, vol. 13, no. 8, pp. 552–564, 2012.
[3] J. Goncalves and S. Warnick, “Necessary and Sufficient Conditions
for Dynamical Structure Reconstruction of LTI Networks,” Automatic
Control, IEEE Transactions on, vol. 53, no. 7, pp. 1670–1674, 2008.
[4] D. Hayden, Y. Yuan, and J. Goncalves, “Network Identifiability from
Intrinsic Noise,” IEEE Transactions on Automatic Control, vol. PP,
no. 99, p. 1, 2016.
[5] P. M. J. Van den Hof, A. Dankers, P. S. C. Heuberger, and X. Bombois,
“Identification of dynamic models in complex networks with prediction
error methods—Basic methods for consistent module estimates,” Automatica, vol. 49, no. 10, pp. 2994–3006, 2013.
[6] A. Chiuso and G. Pillonetto, “A Bayesian approach to sparse dynamic
network identification,” Automatica, vol. 48, no. 8, pp. 1553–1565, 2012.
[7] Z. Yue, W. Pan, J. Thunberg, L. Ljung, and J. Goncalves, “Linear
Dynamic Network Reconstruction from Heterogeneous Datasets,” arXiv
preprint, 2016.
[8] L. Ljung, System Identification: Theory for the User, ser. Prentice-Hall
information and system sciences series. Prentice Hall PTR, 1999.
[9] M. Yuan and Y. Lin, “Model selection and estimation in regression with
grouped variables,” Journal of the Royal Statistical Society: Series B
(Statistical Methodology), vol. 68, no. 1, pp. 49–67, 2006.

(42)

[10] E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by
reweighted l1 minimization,” Journal of Fourier analysis and applications, vol. 14, no. 5-6, pp. 877–905, 2008.
[11] R. Chartrand and W. Yin, “Iteratively reweighted algorithms for compressive sensing,” in Acoustics, speech and signal processing, 2008.
ICASSP 2008. IEEE international conference on. IEEE, 2008, pp.
3869–3872.
[12] M. E. Tipping, “Sparse Bayesian learning and the relevance vector
machine,” The journal of machine learning research, vol. 1, pp. 211–
244, 2001.
[13] D. P. Wipf and B. D. Rao, “An empirical Bayesian strategy for solving
the simultaneous sparse approximation problem,” Signal Processing,
IEEE Transactions on, vol. 55, no. 7, pp. 3704–3716, 2007.
[14] F. He, H. Chen, M. Probst-Kepper, R. Geffers, S. Eifes, A. del Sol,
K. Schughart, A. Zeng, and R. Balling, “PLAU inferred from a correlation network is critical for suppressor function of regulatory T cells,”
Molecular Systems Biology, vol. 8, no. 1, nov 2012.
[15] L. Ljung and A. Wills, “Issues in sampling and estimating continuoustime models with stochastic disturbances,” Automatica, vol. 46, no. 5,
pp. 925–931, may 2010.
[16] H. Garnier, M. Mensler, and A. Richard, “Continuous-time model identification from sampled data: implementation issues and performance
evaluation,” International Journal of Control, vol. 76, no. 13, pp. 1337–
1357, 2003.
[17] J.-F. L. Gall, Brownian Motion, Martingales, and Stochastic Calculus
(Graduate Texts in Mathematics), 1st ed. Springer, 2016.
[18] K. J. Åström, Introduction to stochastic control theory.
Courier
Corporation, 2012.
[19] H. Garnier and L. Wang, Identification of Continuous-time Models from
Sampled Data. Springer London, 2008.
[20] N. Higham, Functions of Matrices. Society for Industrial and Applied
Mathematics, 2008.
[21] R. A. Horn and G. G. Piepmeyer, “Two applications of the theory of
primary matrix functions,” Linear Algebra and its Applications, vol. 361,
pp. 99–106, mar 2003.
[22] Z. Yue, J. Thunberg, and J. Goncalves, “Inverse Problems for Matrix
Exponential in System Identification: System Aliasing,” in 2016 22nd
Proceedings of International Symposium on Machematical Theory of
Networks and Systems, 2016.
[23] K. J. Åström, “Maximum likelihood and prediction error methods,”
Automatica, vol. 16, no. 5, pp. 551–574, 1980.
[24] M. Viberg, “Subspace-based state-space system identification,” Circuits,
Systems and Signal Processing, vol. 21, no. 1, pp. 23–37, 2002.
[25] S. Boyd and L. Vandenberghe, Convex optimization.
Cambridge
university press, 2004.
[26] M. C. Grant and S. P. Boyd, The CVX Users’ Guide, Release v2.1,
2015. [Online]. Available: http://web.cvxr.com/cvx/doc/CVX.pdf
[27] A. H. Al-Mohy and N. J. Higham, “Computing the Fréchet derivative
of the matrix exponential, with an application to condition number
estimation,” SIAM Journal on Matrix Analysis and Applications, vol. 30,
no. 4, pp. 1639–1657, 2009.
[28] A. H. Al-Mohy, N. J. Higham, and S. D. Relton, “Computing the Fréchet
derivative of the matrix logarithm and estimating the condition number,”
SIAM Journal on Scientific Computing, vol. 35, no. 4, pp. C394–C410,
2013.

