On Identification of Dynamical Structure Functions:
A Sparse Bayesian Learning Approach

arXiv:1605.09543v1 [cs.SY] 31 May 2016

J. Jin, W. Pan, D. L.T. Pham, Y. Yuan† , C. J. Tomlin, A. Webb, J. Gonçalves
Abstract— This paper considers the identification of linear
time-invariant networks, also known as dynamic structure
functions. Assuming identifiability of the network addressed
in previous work, this paper presents an identification method
that infers both the Boolean structure of the network and
the transfer functions between nodes. The identification is
performed directly from data and without any prior knowledge
of the system, including its order. The method is to formulate
the identification as a linear regression problem together with
penalties for complexity, both in terms of element (order of
nonzero connections) and group sparsity (network topology).
We then propose a novel scheme that combines sparse Bayesian
and sparse group Bayesian to efficiently solve the problem. The
method and the developed toolbox can now be used to infer
networks from a wide range of fields, including systems biology
applications such as signalling and genetic regulatory networks.

I. INTRODUCTION
Network identification of biochemical reaction networks
has been a central topic in systems biology. Its first objective consists of reconstructing Boolean structures or network topologies. The next step aims at identifying the
internal dynamics between nonzero nodes [2], [3], [13]–
[15]. Such dynamics can be complex, involving nonlinear
ODEs that typically obey biochemical kinetic laws such
as mass action functions, Michaelis-Menten functions and
Hill functions [10]. Building such networks demands prior
knowledge of the exact type of nonlinearities and full state
measurements, which is often unrealistic [16]. Sometimes,
however, the dynamics may be reasonably approximated by
linear systems. A major advantage is that linear systems no
longer require full state measurements. For a network, linear
dynamics between nodes can be described by dynamical
structure functions (DSF), which focuses on the casual relationships among measured states (nodes), while suppressing
the structure of hidden (not-measured) states [11], [12]. This
paper infers DSF from time-series data of the measured
nodes.
Since data usually contains noise, modelling DSF can
lead to model over-fitting. In particular, a simple linear
Junyang Jin and Alex Webb are with Circadian Signal Transduction
Group, Department of Plant Sciences, University of Cambridge. Wei Pan
is with Cardwell Investment Tech and Imperial College London. Duong
L.T. Pham is with the Luxembourg Centre for Systems Biomedicine.
Ye Yuan and Claire J. Tomlin are with the Department of Electrical Engineering and Computer Sciences, UC Berkeley. Jorge Gonçalves
is with the Department of Engineering, University of Cambridge and
the Luxembourg Centre for Systems Biomedicine. † For correspondence,

yy311@berkeley.edu

Due to page limit, the supporting information of this paper is
available from http://hybrid.eecs.berkeley.edu/˜yeyuan/
cdc16.pdf

regression is likely to return a fully connected network with
high dimensional subsystems between nodes. As in classical
system identification, to avoid these situations we need to
penalise for complexity, both in terms of network topology
and individual system dimension. Methods to recover such
DSF can be classified into two categories [18]. Type I
method, called sparse group lasso, solves the problem by
introducing a penalty that blends lasso and group lasso [20].
Type II uses instead the Bayesian perspective called sparse
Bayesian learning (SBL). While element SBL [16], [18],
[22], [25] and group SBL [2], [17], [19], [21] have already
been developed to find the sparse or group-sparse solution of
a linear regression problem, no similar approach is available
to search for a solution possessing both features at the
same time. This paper proposes a scheme to combine SBL
and group SBL, to generate sparse networks with internal
dynamics described by strictly proper transfer functions with
the lowest order possible.
The paper is organised as follows. Section II introduces the
DSF and discusses its identifiability. Section III formulates
the network reconstruction problem and casts it as a linear regression problem. Section IV solves the problem of network
identification using a combination of SBL and group SBL.
Section V extends the formulation to a more general case
and demonstrates the underlying mechanisms of our method.
Section VI applies the method to a randomly generated
network and shows how its topology and internal dynamics
can be reconstructed almost exactly. Finally, Section VII
concludes and discusses further development in this field.
Notation: The notation in this paper is standard. I denotes
the identity matrix. For L ∈ Rn×n , diag{L} denotes a vector
with its elements to be the diagonal elements of matrix L.
tr{L} denotes the trace of matrix L. A matrix L  0 means
L is positive semi-definite. A vector v ≥ 0 means each
element of the vector is non-negative. A vector y(t1 : t2 )
denotes a row vector y(t1 ) y(t1 + 1) · · · y(t2 ) .
II. MODEL FORMULATION

In a particular experiment, the choice of measurements and
what non-measured nodes (hidden nodes) represent influence
the network. This paper assumes the states of the system have
been defined so that the network formed by their interactions
is well defined. In biology, for example, the states typically
are concentrations (or number) of molecules. Consider a
linear model given by its state space dynamics (without
noise):
x(t + 1) = Ax(t) + Bu(t)
(1)
y(t) = Cx(t)

where x ∈ Rn are the states of the system, u ∈ Rm the
inputs, and y ∈ Rp the measurements. To simplify notation,
the paper assumes that the first p states correspond to the
measurements (outputs), followed by hidden states, so that
C can be written as C = [I 0]. For example, in systems
biology there are typically a relatively small number of
measurements and a large number of non-measured (hidden)
species. Hence, the system can be partitioned as

 

 

y(t + 1)
A11 A12
y(t)
B1
=
+
u(t)
h(t + 1)
A21 A22
h(t)
B2
(2)
where h ∈ Rn−p are the hidden states. To represent the system only at the measurements (nodes) level while encoding
the hidden nodes via transfer functions between nodes, we
use the notion of dynamical structure functions (DSF) [1],
[11], [12]. From equation (10), DSF can be written as:
Y
where

and

=

(3)

QY + P U

Q = (zI − D)−1 (W − D)
P = (zI − D)−1 V

where each element in E is i.i.d white Gaussian noise and
they are independent.






e1 (t)
y1 (t)
u1 (t)
 . 




..
Y =  ...  , U = 
 , E =  .. 
.
(7)
ep (t)
yp (t)
um (t)
Q = [Qij (z −1 )],
Qii = 0,
P = [Pij (z −1 )], H = diag(Hii (z −1 )).

We regard each row of the function as an ARMAX
process:
yi (t)
= Qi1 y1 (t) + . . . + Qip yp (t) + . . . + Pim um (t) + Hii Ei (t)
Du (z −1 )
De (z −1 )
Dy (z −1 )
= i1 −1 y1 (t) + . . . + im −1 um (t) + ii −1 ei (t)
Ni (z )
Ni (z )
Ni (z )
(8)
so that
Ni (z −1 )yi (t)

(4)

y
u
e
= Di1
(z −1 )y1 (t) + . . . + Dim
(z −1 )um (t) + Dii
(z −1 )ei
| {z }
MA process

−1

W = A11 + A12 (zI − A22 )

yi (t)

A21

V = A12 (zI − A22 )−1 B2 + B1

(5)

D = diag(W11 , W22 . . . Wpp )

Here, Q is a matrix of transfer functions showing the
connectivity of each node to other nodes. In addition, its
diagonal elements are zero, since we are interested in network connectivity, not self loops. Similarly, P is a matrix
of transfer functions relating the inputs to the nodes. Note
that each element of Q and P is strictly proper indicating
the causality of the network. The boolean structure of the
network is reflected by the nonzero elements in Q and
P , whereas the dynamics of the network are given by the
transfer functions of their elements.
Given time-series data of nodes, the dynamical structure
function is not necessarily identifiable unless partial structure
information is available [1], [11], [12]:
Theorem 1: [11]. (Reconstruction with partial structure):
Given a p × m transfer function G, the dynamical structure
function can be reconstructed from partial structure information
if and only if p − 1 elements in each column of

0
Q P
are known that uniquely
specify
the component


of (Q,P) in the null space of G0 I .
This paper assumes the network is identifiable. In particular,
a sufficient condition for identifiability is that P is diagonal,
i.e., there are p experiments with only one input set to be
non-zero for each experiment.
III. RECONSTRUCTION PROBLEM
FORMULATION
A. Problem Formulation
We add the process and measurement noise into the DSF
model as:
Y = QY + P U + HE
(6)

y
= Di1
(z −1 )y1 (t) + . . . + [1 − Ni (z −1 )]yi (t)

u
e
+ Dim
(z −1 )um (t) + Dii
(z −1 )ei .

(9)

where
Ni (z −1 ) = αi1 z −k + αi2 z −k+1 + . . . + αik z −1 + 1
y
Dij
(z −1 ) = λyij1 z −k + . . . + λyij(k−1) z −2 + λyijk z −1
u
Dij
(z −1 ) = λuij1 z −k + . . . + λyij(k−1) z −2 + λuijk z −1
e
Dii
(z −1 ) = λeii1 z −k + . . . + λyii(k−1) z −1 + λeiik .

(10)

Remark 1: All the terms which are known to guarantee
identifiability can be moved to the left hand side of eq. (9).
The expressions in the rest of the paper are changed accordingly.
Assume that m experiments are conducted with M timeseries data points sampled for each node and input using
Theorem 1. For the ith node, defining


yic (M )


..
Yic = 

.
yic (k + 1)


y1c (M − k : M − 1) ... −yic (M − k : M − 1)


..
..
Aci = 

.
...
.
...
−yic (1 : k)
y1c (1 : k)
  y
T
y

if j <= p − 1
  λij1 . . . λijk 
T
j
u
u
wi =
λij1 . . . λijk
if p <= j <= m + p − 1

T
 
αi1 . . . αik
if j = p + m
 c
T
c
e
−1 c
c
εi (t) = Dii (z )ei (t), Ξi = εi (M ) . . . εci (k + 1)
(11)





wi1
Yi1




..
Yi =  ...  , wi = 

.
Yim
wip+m
 1 
 1 
Ξi
Ai
 .. 
 .. 
Ai =  .  , Ξi =  .  ,
Am
Ξm
i
i


(12)

p(wij ) = max N (wij |0, βij )ϕ(βij ),

eq. (9) can be written as a linear regression model:
Yi = Ai wi + Ξi

virtually all sparse priors namely super Gaussian can be
presented in a variational form which yields a lower bound
for the sparse prior p(w). There are two types of variational
representations; while we apply the convex type which is
more general than the integral one [25].
To induce element sparsity to w, we introduce a prior p(w)
as [16], [18]:

(13)

To reconstruct the whole network, we need to solve p
independent linear regression problems in the same form due
to p nodes. For simplicity, we use a unified expression for
all the nodes:
y = Aw + Ξ
(14)
Now, we aim to solve this linear regression problem given
measured data recorded in Y and A. Note that if the real
network topology is sparse and the system order k is set
to be much bigger than the true one (since we don’t know
the true system order, instead, we set an upper bound for
it), the weighting vector w is both group-sparse (network
is sparse) and element-sparse in each group. Therefore, we
decide to solve the above optimization (14) utilising these
prior knowledge using sparse Bayesian learning.
IV. RECONSTRUCTION VIA GE SPARSE
BAYESIAN LEARNING
A. Sparsity inducing priors
From a Bayesian perspective, all the unknown variables
in the linear regression model (14) are treated as stochastic
random variables and assigned with a probability density
function (PDF) individually. According to eq. (12), each
block in the noise vector Ξ is a random variable following a
Multivariate Gaussian distribution. Since these blocks present
noise from different experiments, they are independent Gaussian. For eq. (9), as the noise in each experiment is generated
by a Moving-Average (MA) process, the covariance matrix
for each block possesses a symmetric Toeplitz structure. As a
result, the vector Ξ ∼ N (0, Π) itself is Multivariate Gaussian
whose covariance matrix has a block diagonal structure with
each block to be a symmetric Toeplitz matrix. In this case, the
likelihood of Y given the weighting vector w is p(Y |w) =
N (Y |Aw, Π). Full Bayesian treatment also requires to introduce a prior distribution for w. Weh define a distribution
i
P
p(w) in general form as: p(w) ∝ exp − 12 j g(wj ) . If w
has certain sparse properties, a prior inducing sparsity like
Generalized Gaussian, Student’s t and Logistic is assigned to
p(w). The function g(·) of such priors is usually a concave,
non-decreasing function of |wj | [18].
Estimating w using maximum a posteriori (MAP) is
intractable in this case because the posterior distribution
p(w|y) is non-Gaussian and not analytical. To simplify the
problem, sparse Bayesian learning approximates p(w|y) with
a Gaussian distribution so that the solution is just E(w|y)
and can be calculated easily. To do that, first note that

βij ≥0

p(wi ) =

k
Y

j=1

p(w) =

p+m
Y
i=1

where

p(wij ) = max N (wi |0, Bi )ϕ(βi ),
βi ≥0

(15)

p(wi ) = max N (w|0, B)ϕ(β).
β≥0




w1
wi1

 . 

..
w=
 , wi =  .. 
.
wik
wp+m




β1
βi1


 . 
..
β=
 , βi =  .. 
.




βik


βp+m


... 0
..  .
..
.
. 
Bp+m
. . . βik
(16)
For group sparsity, the corresponding prior is [8], [17]:
B1
 ..
B= .
0

...
..
.
...

0
..
.



βi1

 ..
 Bi =  .
0

p(wi ) = max N (w|0, γi I)ϕ(γi )
γi ≥0

p(w) =

p+m
Y
i=1

where




γ=




Γ=

γ1
..
.



p(wi ) = max N (w|0, Γ), ϕ(γ)
γ≥0



γi

 ..
 , Γi =  .
γp+m
0

Γ1 . . .
0

..
..
..
.
.
.
.
0

...

(17)

...
..
.
...


0
..  ∈ Rk×k ,
. 
γi

(18)

Γp+m

Remark 2: Note that (17) is for a normal group-sparse w.
Considering the dynamical structure function, the last group
of w contains parameters of the denominator of the transfer
function and hence cannot be all zero. So in fact, it should
be excluded from the group-sparse prior.
To promote both element and group sparsity, we multiply (15) by (17) and normalize it:
p(w) = C ∗ max N (w|0, B)N (w|0, Γ)ϕ(β)ϕ(γ) (19)
γ≥0,β≥0
where C is the normalization constant and can be absorbed
by functions ϕ(β) or ϕ(γ). Hence:

p(w) = max N (w|0, B)N (w|0, Γ)ϕ(β)ϕ(γ).
γ≥0,β≥0

(20)

Proposition 1: Following (25), the estimation w, hyperparameters and noise covariance matrix, Π can be estimated
by solving the optimization problem:

We now get a lower bound of the prior as:
p(w) ≥ N (w|0, B)N (w|0, Γ)ϕ(β)ϕ(γ)
= p̂(w).

(21)

If either βij or γi is 0, the corresponding Gaussian
distribution becomes a Dirac delta function imposing element
or group sparsity to w.
B. Type II Maximization
Although the implicit prior p̂(w) is improper, we can still
get a normalized posteriori distribution of w as:
p(y|w)p̂(w)
p(y|w)p̂(w)dw
N (y|Aw, Π)N (w|0, B)N (w|0, Γ)ϕ(β)ϕ(γ)
R
=
.
p(y|w)p̂(w)dw
(22)
Clearly, p(w|y) is Gaussian since log(p(w|y)) is a
quadratic function of w:
p̂(w|y) = R

p(w|y) = N (w|µw , Σw )

(23)

where


−1
Σw = (Γ−1 + B −1 ) + AT Π−1 A

−1 T −1
µw = (Γ−1 + B −1 ) + AT Π−1 A
A Π y.

(24)

Now, the core of the problem is to choose proper hyperparameters β and γ so that p̂(w|y) ≈ p(w|y) to some extend.
One way is to minimize the misaligned mass between p(w)
and p̂(w) weighted by the marginal likelihood p(y|w) which
is also called evidence maximization or Type II maximization
[9], [16], [18]. It is equivalent to estimating all the unknown
variables using maximum likelihood method:
(γ, β, Π) = arg
= arg

min

β,γ≥0,Π0

max

Z

Z

p(y|w)|p(w) − p̂(w)|dw

p(y|w)p̂(w)dw
β,γ≥0,Π0
Z
= arg min −2 log p(y|w)p̂(w)dw
min

β,γ≥0,Π0

−2 log p̂(y|β, γ, Π)

= min (y − Aw)T Π−1 (y − Aw) + wT (Γ−1 + B −1 )w
β,γ,Π,w

+ log |B + Γ| + log|Π + A(Γ−1 + B −1 )−1 AT |

Subject to:

β ≥ 0, γ ≥ 0, Π  0

(26)
where we set p(·) = −2 log ϕ(·) = const.
Proof: The derivation can be found in Section A of
Appendix.
Now, assume
L = (Γ−1 + B −1 )−1
(27)
so that
where

L ≤ Γ, B = (L−1 − Γ−1 )


li1
 ..
Li =  .
0

...
..
.
...



0
..  , L = 

. 
lik

L1
..
.
0


... 0
.. 
..
.
. 
. . . Lp

(28)

(29)

Assuming l = diag(L) and substituting (27) into (26) lead
to:
L(γ, l, Π, w)
= min (y − Aw)T Π−1 (y − Aw) + wT L−1 w
Γ,L,Π,w

+ log|Π + ALAT | + 2log|Γ| − log|Γ − L|

(30)

subject to:

0  L  Γ, γ ≥ 0, Π  0.

C. Algorithm to solve Type II maximization
Let:
u(w, Π, l, γ)

= (y − Aw)T Π−1 (y − Aw) + wT L−1 w − log|Γ − L|

v(Π, l, γ)


= − log|Π + ALAT | + 2log|Γ| .

β,γ≥0,Π0

= arg

L(γ, β, Π, w)

The optimization problem (30) becomes:
(25)

Remark 3: Π matrix has a block-diagonal symmetric
Toeplitz structure. It is parameterized by mk(2M − k + 1)/2
variables. Estimation of Π is very important because the
recovery performance can be poor if estimated Π is suboptimal [8]. To avoid over-fitting, we assume 1) white
Gaussian noise e(t) in each experiment has same variance
(number of variables reduces to k(2M − k + 1)/2). 2)
Each diagonal block of Π is well approximated by its main
diagonal and sub-diagonal elements. Consequently, Π is only
parameterized by two variables.

L(γ, l, Π, w) =

min

Γ0,0LΓ,Π0

(31)

u(w, Π, l, γ) − v(Π, l, γ).

(32)
Proposition 2: Functions u(w, Π, l, γ) and v(Π, l, γ) are
both jointly convex with respect to their own variables [5].
Proof: The derivation can be found in Section B of
Appendix.
As a result, the cost function actually is a difference of
two convex functions and thus is a difference of convex
programming (DCP) problem. It can be solved using sequential convex optimization techniques. Here, we use convexconcave procedure (CCCP) which belongs to majorizationminimization (MM) algorithm using the linear majorization

function [4], [24]. For minx f (x) where f (x) = u(x)−v(x),
we can solve it iteratively by:
x

n+1

n

= arg min u(x)− < x, ∇f (x ) >
x

−

Γ,L,Π,w
T
∇l v(Πn , ln , γ n )l

subject to:

− ∇Tγ v(Πn , ln , γ n )γ

(34)

T

T −1

(35)

To alleviate the computation burden at each iteration and
generate reasonable estimation of Π, we separate the update
of variables in (32) into two steps using Nonlinear Gaussian
Seidel method [23].
Step 1:


wn+1
 ln+1  = arg
min
u(w, Πn , l, γ)
Γ0,0LΓ
n+1
γ
Step 2:

− ∇Tl v(Πn , ln , γ n )l − ∇Tγ v(Πn , ln , γ n )γ.
(36)

Πn+1 = arg min u(wn+1 , Π, ln+1 , γ n+1 )
Π0


− tr ∇Π v(Πn , ln+1 , γ n+1 )Π .

M m×k

Ap+m

.



,

(39)




w
 l  = arg min u(w, Π, l, γ) − ∇Tl v(Πn , ln , γ n )l
Γ,L
γ
(40)
X
T
n n
n
− ∇γ v (Π , l , γ )γ +
g(γi , li )
where

γp−1



...

Next, we transform (36) into a unconstrained problem by
introducing a indicator function:

g(γi , li ) =

−∇l v = diag{A (Π + ALA ) A}
 −1 
γ1
 .. 
−1
−1
−∇γ v = 2kγ , γ =  .  .

A1

Ai ∈ R

T −1

−∇Π v = (Π + ALA )



i

0  L  Γ, γ ≥ 0, Π  0,

where

A=

(33)

where < ·, · > denotes inner product.
Therefore, the optimization problem (30) can be decomposed into sequential convex optimization problems:

γ n+1 , ln+1 , Πn+1 , wn+1
= arg min u(w, Π, l, γ) − tr [∇Π v(Πn , ln , γ n )Π]

and

(37)

D. Solve CCCP by Alternating Direction Method of Multipliers (ADMM)
If the network is large, solving (36) using toolbox like
CVX may encounter memory difficulty because the dimension of w is huge. In this case, we consider to split
the optimization problem into a series of independent subproblems with lower variable dimension. It turns out that the
split cost function can be treated as a sharing problem and
as a result can be solved using ADMM algorithm [6], [7].
First of all, we partition and rearrange variables in (36):


w1
 l1 




 γ1 
w



..
 l ⇒
(38)

 ∈ R(p+m)(2k+1)×1 ,
.


γ
 wp+m 


 lp+m 
γp+m



0 if li ≥ 0, γi ≥ 0, li ≤ γi
+∞ else

(41)

Now, we can split the cost function as:



wi
 li  = arg min
Γ,L,w
γi
"p+m
#
X
T
T
wiT L−1
i wi − log|Γi − Li | − ∇li li − ∇γi γi + g(γi , li )
i=1

+ (y −

p+m
X

Ai wi )T Π−1 (y −

p+m
X

zi )T Π−1 (y −

i=1

p+m
X

Ai wi ).

i=1

(42)
The optimization problem (42) can be then expressed as
a sharing problem [6]:


wi
 li  = arg min
Γ,L,w
γi
"p+m
#
X
T −1
T
wi Li wi − log|Γi − Li | − ∇li li − ∇γi γi + g(γi , li )
i

+ (y −

i=1

subject to

p+m
X

zi )

i=1

Ai wi − zi = 0.

(43)
T
Let xi = wiT liT γiT
, such problem can be solved
using scaled ADMM method as [6]:


xn+1
:= arg
i
−

wiT L−1
i wi − log|Γi − Li |

min

γi ≥0,0≤li ≤γi ,wi
∇Tli li − ∇γi γi +

z n+1 := arg min(y −
z

+ (ρ/2)

n
X
i

X
i

(ρ/2)kAi wi − zin + uni k22
X
zi )T Π−1 (y −
zi )
i

kzi − uni − Ai win+1 k22

un+1
:=uni + Ai win+1 − zin+1 .
i

(44)

Clearly, the z-update still requires to solve a problem in
M m(p + m) variables. However, we can further simplify
the formulation so that the number of variables is reduced
to M m [6]:

It is easy to see that prior (48) is still super-Gaussian which
means it is a sparsity prior. If we rewrite the lower bound
of prior as p̂(w) = p(w|L)p(L|Γ)P (Γ) where p(L|Γ) =
ϕ(l; γ) and p(Γ) = ϕ(γ), the prior actually employs hierarchical Bayesian format. Obviously, l takes charge of
n+1
T −1
T
xi
:= arg min wi Li wi − log|Γi − Li | − ∇li li − ∇γi γi element sparsity of w. p(L|Γ) brings an underlying constraint
γi ,li ,wi
0  L  Γ which controls sparsity of L. ϕ(γ) = |Γ|−1
n
n
n 2
n
is similar to a Jeffrey prior inducing element sparsity to γ
+ (ρ/2)kAi wi − Ai wi + Aw − z + u k2
thus enforcing group sparsity to l via p(L|Γ). Therefore, two
subject to:
hyper-parameters in (20) controls element and group sparsity
γi ≥ 0, 0 ≤ li ≤ γi
of w directly while in (48), γ imposes group sparsity to w
T
indirectly via l.
z n+1 := arg min [y − (p + m)z] Π−1 [y − (p + m)z]
z
The discussion above indicates that it is possible to conn+1 2
+ [(p + m)ρ/2] kz − un − Aw
k2
sider more features of w at the same time by using hierarn+1
chical Bayesian prior. We can assume n hyper-parameters
un+1 :=uk + Aw
− z n+1 ,
(45) each of which stands for one feature of w and impose such
feature to the hyper-parameter ranked ahead. Consequently,
where
p+m
X
the prior of w induces all these features to the estimation.
n
(46)
Ai win .
Aw = 1/(p + m)
VI. S IMULATION
i=1

To guarantee the identifiability of the network, we assign a
special structure to the matrix P . We assume there are same
number of input as nodes and set P matrix to be diagonal.
T
xn+1
:= arg min wiT L−1
i
i wi − log|Γi − Li | − ∇li li − ∇γi γi Such a network is identifiable following immediately from
γi ,li ,wi
n
the Theorem 1 because p − 1 elements in each column of P
+ (ρ/2)kAi wi − Ai win + Aw − z n + un k22
are known 0.We consider here an LTI system with dynamic
subject to:
structure functions Q and P as in eq. (10) with 10 nodes.
γi ≥ 0, 0 ≤ li ≤ γi
We set the maximal order of denominators in Q, P to


n+1
−1
n+1
n
k
=
5. Time series data from 10 experiment with 30 meaz
:= [2(p + m)I + ρΠ]
2y + ρΠu + ρΠAw
surements for each is collected when the input to the system
n+1
un+1 :=un + Aw
− z n+1 .
is Gaussian random noise. Since there are 10 measured states
(47) and each measured state is controlled by one input, the
Note that x-update can be solved in parallel independently. dictionary matrix has (5 ∗ 11) = 55 columns and 30 rows.
z and u-update are then solved in sequence after collecting We can correctly identify the topology of Q by the proposed
x-update.
method.
For comparison purpose, the Type I sparse group lasso
V. E XTEND TO GENERAL FORMULATION
[20] is also applied. Figure 1 shows that Type I algorithm
Inspired by the final optimization problem (30), we can has a poor performance, as it picks many incorrect links
deduce a general formulation of sparsity prior which leads and misses a correct link. Numerical values of true and
to the same Type-II maximization problem.
identified dynamical structure functions are given in the
Let:
Supporting Information 1 . The estimated LTI model indicates
the performance of the proposed algorithm can be further
p(w) = max N (w|0, B)ϕ(β; γ)ϕ(γ)
β≥0,γ≥0
(48) improved if a proper rule to prune out small elements of w
is adopted.
z-update can be solved analytically so the algorithm ends
with:

so that
p̂(w) = N (w|0, B)ϕ(β; γ)ϕ(γ).

(49)

Following the same Bayesian perspective, we end with the
optimization problem:
(γ, l, Π, w)
= arg

min

(y − Aw)T Π−1 (y − Aw) + wT L−1 w

Γ0,0LΓ,Π0
T

+ log |Π + ALA | − 2 log ϕ(l; γ) − 2 log ϕ(γ).

(50)
Clearly, if we set ϕ(l; γ) = |Γ − L|1/2 and ϕ(γ) = |Γ|−1 ,
(50) is same with (30).

VII. CONCLUSION AND DISCUSSION
This paper proposes a method to identify the dynamical
structure function assuming its identifiability. The problem
is formulated as a linear regression model and solved using
sparse Bayesian learning. Given measured time series data
from the nodes of the network, the method is able to infer
the topology and internal dynamics of the network without
any prior knowledge of the system, including its topology
or system order. The newly proposed method penalises
the complexity of the network by introducing both group
1 the supporting information of this paper is available from http://
hybrid.eecs.berkeley.edu/˜yeyuan/cdc16.pdf

R EFERENCES

Fig. 1: Topology of Q matrices. Red x: identified nodes. Blue
diamond: True.
and element sparsity. These two kinds of sparse priors are
combined and further extended to a general form which
reveals their rationality.
To guarantee identifiability, we assumed the data originated from p (same number as measurements) experiments
and that only one input is non-zero for each experiment.
However, in reality it may not be feasible to run so many
experiments, or having such specific inputs. In that case, the
network would probably not be identifiable. Nevertheless,
future work will investigate which classes of systems can still
reliably be inferred with this method. We expect there will be
networks, especially those that are sparse, that can achieve
good estimation with a much smaller number of inputs.
VIII. ACKNOWLEDGEMENT
Duong L.T. Pham and Jorge Goncalves acknowledge
the support from the Luxembourgish National Research
Fund (FNR 14/BM/8231540). Ye Yuan and Claire J. Tomlin are supported in part by NSF under CPS:ActionWebs
(CNS-0931843) and CPS:FORCES (CNS1239166), by
NASA under grants NNX12AR18A and UCSCMCA14-022 (UARC), by ONR under grants N00014-12-10609, N000141310341 (Embedded Humans MURI), and
MIT 5710002646 (SMARTS MURI), and by AFOSR under
grants UPenn-FA9550-10-1-0567 (CHASE MURI) and the
SURE project.

[1] H. Weerts, A. Dankers, P. Hof Identifiability in dynamic network
identification, IFAC-PapersOnLine, 48(28), 2015, pp. 1409-1414.
[2] B. Sanandaji, T. Vincent, M. Wakin A review of sufficient conditions
for structure identification in interconnected systems, Proceedings of
the 16th IFAC Symposium on System Identification, 2012, pp. 16231628.
[3] N. Everitt, G. Bottegal, C. Rojas, H. Hjalmarsson On the Effect of
Noise Correlation in Parameter Identification of SIMO Systems, IFACPapersOnLine, vol. 48(28), 2015, pp. 326-331.
[4] T. Chen, M.S. Andersen, L. Ljung, A. Chiuso, and G. Pillonetto,
System identification via sparse multiple kernel-based regularization
using sequential convex optimization techniques, IEEE Trans. Automat. Contr., vol. 59(11), 2014, pp. 2933-2945.
[5] S. Boyd, and L. Vandenberghe Convex Optimization, Cambridge
University Press; 2004.
[6] S. Boyd, Distributed Optimization and Statistical Learning via the
Alternating Direction Method of Multipliers, Found. Trends Mach.
Learn., vol. 3(1), 2010, pp. 1-122.
[7] W. Pan, A. Sootla, and G. B. Stan, Distributed reconstruction of
nonlinear networks: An ADMM approach, IFAC Proc. Vol, vol. 19,
2014, pp. 3208-3213.
[8] Z. Zhang and B. D. Rao, Extension of SBL algorithms for the recovery
of block sparse signals with intra-block correlation, IEEE Trans.
Signal Process., vol. 61(8), 2013, pp. 2009-2015.
[9] C. Bishop. Pattern recognition and machine learning, Springer New
York; 2006.
[10] D. J. Wilkinson. Stochastic Modelling for Systems Biology, CRC
Press; 2012.
[11] J. Goncalves and S. Warnick Necessary and sufficient conditions
for dynamical structure reconstruction of LTI networks, IEEE Trans.
Automat. Contr., vol. 53(7), 2008, pp. 1670-1674.
[12] Y. Yuan, G. Stan, S. Warnick, and J. Goncalves, Robust dynamical
network structure reconstruction, Automatica, vol. 47, 2011, pp. 1230–
1235.
[13] 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, 2013, pp. 2994-3006.
[14] A. Dankers, P. M. J. Van den Hof, P. S. C. Heuberger and X.
Bombois, Dynamic network identification using the direct predictionerror method, Proceedings of 51st IEEE Conference on Decision and
Control, 2012.
[15] D. Materassi and M. V. Salapaka, On the problem of reconstructing
an unknown topology via locality properties of the wiener filter,IEEE
Trans. Automat. Contr., vol. 57(7), 2012, pp. 1765–1777.
[16] W. Pan, Y. Yuan, J. Goncalves, and G. Stan A Sparse Bayesian
Approach to the Identification of Nonlinear State-Space Systems,
IEEE Trans. Automat. Contr., vol. 61(1), 2016, pp. 182-187.
[17] W. Pan, Y. Yuan, L. Ljung, J. Goncalves, and G. Stan Identifying
biochemical reaction networks using heterogeneous datasets, CDC,
2015.
[18] D. P. Wipf, B. D. Rao, and S. Nagarajan Latent Variable Bayesian
Models for Promoting Sparsity, IEEE Trans. Inf. Theory., vol. 57(9),
2011, pp. 6236-6255.
[19] S. Babacan, S. Nakajima, and M. Do Bayesian group-sparse modeling
and variational inference, Signal Process. IEEE Trans., vol. 62(11),
2014, pp. 2906-2921.
[20] N. Simon, J. Friedman, T. Hastie, and R. Tibshirani A Sparse-Group
Lasso, J. Comput. Graph. Stat., vol. 22(2), 2013, pp. 231-245.
[21] K. Lee and R. Chen BSGS: Bayesian Sparse Group Selection, The R
Journal., vol. 7(2), 2015, pp. 122-133.
[22] M. Tipping Sparse Bayesian Learning and the Relevance Vector Mach,
J. Mach. Learn. Res., vol. 1, 2001, pp. 211-244.
[23] C. Nazionale and D. Ricerche On the convergence of the block
nonlinear Gauss-Seidel method under convex constraints, 1998
[24] B. K. Sriperumbudur On the Convergence of the Concave-Convex
Procedure, Nips., vol. 1, 2009, pp. 1-9.
[25] J. Palmer, D. Wipf, K. Kreutz-Delgado, and B. Rao Variational EM
algorithms for non-Gaussian latent variable models, Adv. Neural Inf.
Process. Syst., vol. 18, 2006, pp. 1059.

A PPENDIX

Inserting (56) back to (25), we get:

A. Derivation of Type II maximization problem
Firstly, note that:
Z
− 2 log p(y|w)p̂(w)dw
Z
− 2 log N (y|Aw, Π)N (w|0, B)N (w|0, Γ)ϕ(β)ϕ(γ)dw
Z
− 2 log exp(Ew )dw + log |Π| + log |B| + log |Γ|dw
where

(51)

1
Ew = − [(y − Aw)T Π−1 (y − AW ) + (wT (B −1 + Γ−1 )w)]
2
(52)
ignoring all the constant terms.
By completing the square:
1
1
− (y − Aw)T Π−1 (y − Aw) − wT (Γ−1 + B −1 )w
2
2

1
T −1
= − (w − µw ) Σw (w − µw ) + Ey
2
(53)
where

−1
Σw = (Γ−1 + B −1 ) + AT Π−1 A

−1 T −1
µw = (Γ−1 + B −1 ) + AT Π−1 A
A Π y

Ey = min(y − Aw)T Π−1 (y − Aw) + wT (Γ−1 + B −1 )w
w
(54)
In addition, note that:


Z
1
− 2 log exp − [(w − µw )T Σ−1
(w
−
µ
)]
dw
w
w
2
(55)
= − log |Σw |
= log (Γ−1 + B −1 ) + AT Π−1 A

ignoring all the constant terms.
Using (53) and (55),the integral in (51) becomes:
Z
− 2 log p(y|w)p̂(w)dw
= Ey + log (Γ−1 + B −1 ) + AT Π−1 A

+ log |Π| + log|B| + log |Γ|

= Ey + log (Γ−1 + B −1 ) + AT Π−1 A − log|(Γ−1 + B −1 )|

+ log |(Γ−1 + B −1 )| + log |Π| + log |B| + log |Γ|
= Ey + log I + (Γ−1 + B −1 )−1 AT Π−1 A

+ log |(Γ + B)| + log |Π|

= Ey + log I + Π−1 A(Γ−1 + B −1 )−1 AT

+ log |(Γ + B)| + log |Π|

= Ey + log Π + A(Γ−1 + B −1 )−1 AT + log |(Γ + B)|
= min(y − Aw)T Π−1 (y − Aw) + wT (Γ−1 + B −1 )w
w

+ log Π + A(Γ−1 + B −1 )−1 AT + log |(Γ + B)|.

(56)

L(γ, β, Π, w)

= min (y − Aw)T Π−1 (y − Aw) + wT (Γ−1 + B −1 )w
β,γ,Π,w

+ log |B + Γ| + log |Π + A(Γ−1 + B −1 )−1 AT |

Subject to:

β ≥ 0, γ ≥ 0, Π  0.

(57)

B. Convexity of DCP problem
To see functions u(w, Π, l, γ) and v(Π, l, γ) are jointly
convex functions, we need to prove each term is jointly
convex.
For u(w, Π, l, γ), we see that
Π  0,

(y − Aw)T Π−1 (y − Aw) < t

equivalent to

Π
(y − Aw)T

y − Aw
t



(58)
0

so the term (y − Aw)T Π−1 (y − Aw) is jointly convex as is
same with wT L−1 w.
For the term − log |Γ − L|, firstly − log | · | is a convex
function in S + . Since Γ − L is an affine function of Γ and
L, − log |Γ − L| is jointly convex. The convexity of function
v(Π, l, γ) can be seen in the same way as function − log |Γ−
L|.

Supporting Information
C. Dynamical structure function
True dynamical structure functions that were used for data generation
y1 (t) =

1 + 0.01z −1
0.1427 + 0.275z −1 + 0.42z −2
u1 (t) +
e1 (t)
−1
−2
−3
−1
1 + 0.4z − 0.27z + 0.04z
1 + 0.4z − 0.27z −2 + 0.04z −3
−0.0153 + 0.3862z −1
1 + 0.01z −1 + 0.02z −2
u
(t)
+
e2 (t)
2
1 + 0.5z −1 + 0.06z −2
1 + 0.5z −1 − 0.06z −2

y2 (t) =
y3 (t) =

−0.3z −1 − 0.5z −2
0.831 − 0.707z −1 + 0.61z −2
1
y2 +
u3 (t) +
e3 (t)
−1
−2
−3
1 + 1.55z + 0.8z + 0.15z
1 + 1.55z −1 + 0.8z −2 + 0.15z −3
1 + 1.55z −1 + 0.8z −2 + 0.15z −3
y4 (t) =

0.8195 − 0.68z −1 + 0.59z −2
1 + 0.89z −1
−0.2z −1 − 0.1z −2
y2 +
u4 (t) +
e4 (t)
−1
−2
−1
−2
1 − 0.7z + 0.1z
1 − 0.7z + 0.1z
1 − 0.7z −1 + 0.1z −2
y5 (t) =

−1.502 + 1.627z −1
1 + 0.075z −1
u
(t)
+
e5 (t)
5
1 + 1.4z −1 + 0.98z −2
1 + 1.4z −1 + 0.98z −2
y6 (t) =

y9 (t) =

1.214
1 + 0.02z −1
u6 (t) +
e6 (t)
−1
1 + 0.35z
1 + 0.35z −1

y7 (t) =

−0.35
−0.7642
1 + 0.03z −1
y
(t)
+
u
(t)
+
e7 (t)
6
7
1 + 0.55z −1
1 + 0.55z −1
1 + 0.55z −1

y8 (t) =

−0.55
0.038
1 + 0.049z −1
y6 (t) +
u8 (t) +
e8 (t)
−1
−1
1 + 0.49z
1 + 0.49z
1 + 0.49z −1

−0.35
−0.49
0.9347
1 + 0.01z −1
y6 (t) +
y8 (t) +
u9 (t) +
e9 (t)
−1
−1
−1
1 + 0.6z
1 + 0.6z
1 + 0.6z
1 + 0.6z −1

0.8
1 + 0.04z −1
−0.25
y
(t)
+
u
(t)
+
e10 (t)
6
10
1 + 0.8z −1
1 + 0.8z −1
1 + 0.8z −1
The identified dynamic structure functions by the proposed method are
y10 (t) =

y1 (t) =

0.15 + 0.28z −1 + 0.402z −2
u1 (t)
1 + 0.366z −1 − 0.26z − 2 + 0.035z −3

y2 (t) =
y3 (t) =

−0.02 − 0.3805z −1
u2 (t)
1 + 0.469z −1 − 0.05383z −2

0.8195 − 0.68z −1 + 0.59z −2
−0.3005z −1 − 0.48z −2
y
+
u3 (t)
2
1 + 1.56z −1 + 0.7913z −2 + 0.105z −3
1 + 1.56z −1 + 0.7913z −2 + 0.105z −3
y4 (t) =

−0.1909z −1 − 0.67z −2
0.8195 − 0.68z −1 + 0.59z −2
y2 +
u4 (t)
−1
−2
1 − 0.751z + 0.1089z
1 − 0.751z −1 + 0.1089z −2
y5 (t) =

−1.532 + 1.67z −1
u5 (t)
1 + 1.4z −1 + 0.98z −2

y6 (t) =
y7 (t) =
y8 (t) =

1.19
u6 (t)
1 + 0.36z −1

−0.2958
−0.7625
y6 (t) +
u7 (t)
1 + 0.54z −1
1 + 0.54z −1

−0.55
0.0341
y6 (t) +
u8 (t)
1 + 0.475z −1 + 0.026z −2
1 + 0.475z −1 + 0.026z −2

y9 (t) =

−0.4829 + 0.0085z −1
0.9147
−0.3413
y6 (t) +
y8 (t) +
u9 (t)
−1
1 + 0.602z
1 + 0.602z −1
1 + 0.602z −1
y10 (t) =

−0.21
0.76 + 0.04z −1
y
(t)
+
u10 (t)
6
1 + 0.73z −1
1 + 0.73z −1

From the identified dynamic structure functions, it can be seen that the causal relationships between (y3 , y2 ); (y4 , y2 );
(y7 , y6 );(y8 , y6 );(y9 , y8 , y6 ) are almost correctly identified.

