On Identification of Sparse Multivariable ARX Model:
A Sparse Bayesian Learning Approach ?
J. Jin, Y. Yuan† , W. Pan, D. L.T. Pham, C. J. Tomlin, A. Webb, J. Goncalves

arXiv:1609.09660v1 [] 30 Sep 2016

Abstract
Although system identification of complex networks is widely studied, most of the work focuses on modelling the dynamics of
the ground truth without exploring the topology because in many applications, the network topology is known as a priori. For
instance, in industry since the system is artificial, its topology is fixed and the system dynamics is identified for the purpose
of control. Nevertheless, in other cases, especially natural sciences, finding the correct network topology sometimes is even
more important than modelling potential dynamics since it helps scientists understand the underlying mechanism of biological
systems. More importantly, one of the most crucial features of such networks is sparsity which causes a huge difference if
no effort is put on exploring the sparse structure. To fill the gap, this paper begins with considering the identification of
sparse linear time-invariant networks described by multivariable ARX models. Such models possess relatively simple structure
thus used as a benchmark to promote further research. With identifiability of the network guaranteed, this paper presents
an identification method that infers both the Boolean structure of the network and the internal dynamics between nodes.
Identification is performed directly from data without any prior knowledge of the system, including its order. The proposed
method solves the identification problem using Maximum a posteriori estimation (MAP) but with inseparable penalties for
complexity, both in terms of element (order of nonzero connections) and group sparsity (network topology). Such an approach
is widely applied in Compressive Sensing (CS) and known as Sparse Bayesian Learning (SBL). We then propose a novel
scheme that combines sparse Bayesian and group sparse Bayesian to efficiently solve the problem. The resulted algorithm has
a similar form of the standard Sparse Group Lasso (SGL) while with known noise variance, it simplifies to exact re-weighted
SGL. The method and the developed toolbox can be applied to infer networks from a wide range of fields, including systems
biology applications such as signaling and genetic regulatory networks.
Key words: System Identification; Sparse Bayesian Learning; Polynomial Model.

1

Introduction

In system identification, the prediction error method
(PEM) equivalent to Maximum-likelihood (ML) is well
developed. It is applicable to a big family of black-box
models such as ARX, ARMAX and Box-Jenkins. Given
an infinite number of data points and some mild assumptions, PEM attains consistent and asymptotically estimates [11]. However, the main disadvantage for PEM is
that for most model structures it has to deal with nonlinear numerical optimization for the parameter estima? Junyang Jin and Alex Webb are with Circadian Signal
Transduction Group, Department of Plant Sciences, University of Cambridge. Ye Yuan is with School of Automation,
Huazhong University of Science and Technology. Wei Pan is
with Cardwell Investment Tech and Imperial College London. Duong L.T. Pham is with the Luxembourg Centre for
Systems Biomedicine. 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,
yye@hust.edu.cn.

Preprint submitted to Automatica

tion, which normally only guarantees local convergence.
Some exceptions exist for the models in linear regression form such as ARX and FIR model whose optimization problem is solved analytically [28]. The second issue is that without any prior knowledge about the system order, over-fitting is essentially inevitable. Hence,
most proper model is selected afterwards using AIC or
BIC criterion. As such, the procedure demands more
computation effort than necessary. Finally, when dealing with MIMO system, assuming no prior knowledge of
the topology, PEM generates full transfer matrices even
if the ground truth is sparse. The advanced technique,
MAP (Type I method), thus is proposed to penalize the
parameters. While it effectively relieves over-fitting, the
weighting variable which controls the trade-off between
data fitting and model complexity is evaluated by crossvalidation which increases the computation burden and
causes information waste.
In the machine learning literature, a well developed technique called Sparse Bayesian Learning (Type II method)
actually is a general form of MAP and suitable to solve
linear regression problems [1,5,20,25]. Since it applies in-

29 September 2016

separable prior, Type I method is a special case of SBL.
Compared with conventional Type I method, SBL has
superior performance and is tuning-free [25]. In previous research, it is already applied to solve nonlinear system identification [14, 15]. In terms of the model structure of nonlinear systems, element SBL or group SBL
(GSBL) is utilized alone to cope with different data set.
For a general SISO system, approximation of an IIR system by a FIR one using truncation usually is necessary
and the kernel method which generalizes SBL is proposed [6]. To infer a MIMO system described by a dynamical structure function (DSF), only the topology of
the network is inferred by Lasso algorithm using steady
state time-series data while internal dynamics remains
unknown [10]. A similar work is done to infer the multivariable ARX model with known polynomial order using Block Orthogonal Matching Pursuit (BOMP) which
focuses more on identifying network topology [16]. A
non-parametric technique via Bayesian approach is developed to infer the Granger causality and intended to
get correct topology and good prediction [7].

2

MODEL FORMULATION

The network considered is described by a multivariable
ARX model: A(z −1 )Y (t) = B(z −1 )U (t) + E(t), where
A(z −1 ) = I + A1 z −1 + ... + Ana z −na ,
B(z −1 ) = B1 z −1 + ... + Bnb z −nb .

(1)

Y (t) ∈ Rp are the nodes of the network, U (t) ∈ Rm denotes input, and E(t) ∈ Rp is i.i.d white Gaussian noise.
A(z −1 ) is a polynomial matrix showing the connectivity
of each node to other nodes including self-loops. Similarly, B(z −1 ) is a polynomial matrix relating the input
to the nodes. The boolean structure of the network is reflected by the nonzero elements in A(z −1 ) and B(z −1 ),
whereas the dynamics of the network are given by the
elements in these transfer matrices. The multivariable
ARX model defines the network topology and results in
a unique Input-Output map:

In this paper, we consider the identification of a network
described by a multivariable ARX model which has the
simplest model structure in the MIMO family. We also
assume no prior knowledge of the system order and network topology. By the virtue of the model structure, SBL
is introduced to solve this identification problem. To promote sparse network topology and search for lowest order possible, whereas Sparse Group Lasso (SGL) can be
used to solve the same problem as a Type I method, we
propose a novel scheme to combine SBL and GSBL. The
resultant algorithm indicates such approach is closely
linked to SGL method.

where

Y (t) = G(z −1 )U (t) + H(z −1 )E(t),

(2)

G(z −1 ) = A−1 (z −1 )B(z −1 )
H(z −1 ) = A−1 (z −1 ).

(3)

Another useful way to describe a linear network is by
dynamical structure function (DSF) which is uniquely
defined by ARX model in this case.
Y (t) = Q(z −1 )Y (t) + P (z −1 )U (t) + H(z −1 )E(t), (4)
where

Q(z −1 ) = I − Ã−1 (z −1 )A(z −1 )
P (z −1 ) = Ã−1 (z −1 )B(z −1 )

The paper is organized as follows. Section 2 introduces
the polynomial model and discusses its identifiability.
Section 3 formulates the network reconstruction problem. Section 4 solves the problem of network identification using a combination of SBL and GSBL. Section 5
considers an extended nonlinear polynomial model. Section 6 applies the method to a randomly generated network and shows how its topology and internal dynamics
can be reconstructed accurately. Finally, Section 7 concludes and discusses further development in this field.

H(z −1 ) = Ã−1 (z −1 )

(5)

Ã(z −1 ) = diag{diag{A(z −1 )}}.
Remark 1 : In general, the map between DSF and
Input-Output transfer matrix is not bijective. A DSF
uniquely defines a transfer matrix whereas a transfer matrix corresponds to infinite many DSF. Given
a transfer matrix, its associated DSF is unique if and
only if partial information is available about the matrix
[ P (z −1 ) H(z −1 ) ] [9, 22, 26]. For an ARX model, since
each node is only influenced by its own intrinsic noise,
the transfer matrix H(z −1 ) is diagonal so that DSF is
unique.

Notation: The notation in this paper is standard. I denotes the identity matrix. If L ∈ Rn×n , diag{L} denotes
a vector which consists of diagonal elements of matrix
L. If l ∈ Rn , diag{l} denotes a diagonal matrix with its
diagonal elements to be the vector l. trace{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 v k means array
power
of k of v. A vector y(t
h
i 1 : t2 ) denotes a row vector
y(t1 ) y(t1 + 1) · · · y(t2 ) . (.∗ ) means array multiplication.

The theory of identifiability of a general black-box SISO
system is clear [11]. In addition to some necessary conditions to the model structure, not all but at least one polynomial order of rational transfer functions must coincide
with the ground truth. The identifiability of a general
MIMO system is more subtle. It is discussed with specific model classes instead of a general model structure.

2

Also, the conditions vary a lot under different parameterization schemes. Nevertheless, we find the identifiability of a multivariable ARX model is straightforward.

tors:



y=


Lemma 2 : Consider the multivariable ARX model in
eq. (2). The degrees of polynomials are na ≥ n∗a and
nb ≥ n∗b where θ∗ , n∗a and n∗b are ground truth. Then
θ∗ corresponds to a globally identifiable θ of the model
structure.





yi (t)
..
.
yi (k + 1)





w1
..
.



,w = 




wp+m







y1 (t − k : t − 1) ... −yi (t − k : t − 1)


..
..

Φ=
.
...
.


y1 (1 : k)

h
iT

ij
ij

a
.
.
.
a

1
k

h
iT
ij
ij
wj =
.
.
.
b
b
1
k


h
iT


 aii . . . aii
1
k

The proof is straightforward. Given H(z −1 ), A(z −1 ) is
fixed which in turn deduces B(z −1 ) given G(z −1 ). The
lemma indicates an important point that no prior knowledge of the polynomial order is needed aside from an upper bound to guarantee the identifiability of the model
structure.

...

−yi (1 : k)

(8)

if j ≤ p − 1, j 6= i

if p ≤ j ≤ p + m − 1
if j = p + m

2

λ = E{ei (t) }.

The likelihood function based on Baye’s rule is thus:

3

p(y w, λ) =

RECONSTRUCTION PROBLEM FORMULATION

l=0

Remark 3 : The likelihood function of multiple independent experiment is just the multiplication of likelihood
of each experiment. If the experiment condition does not
change too much so that noise variance is same, the likelihood function is in the same form of (9) whose y and Φ
are stack of counterparts of individual experiment. Even
if this assumption does not apply, only little modification
is needed to fit the framework in this paper.

yi (t) denotes ith node, u(t) input and ei (t) i.i.d Gaussian
noise and:
Aii (z

)=

Aij (z

−1

)=

Bij (z −1 ) =
k≥

(9)

where Yt2t1 means data points of all the nodes and input
from t1 to t2 .

yi (t) = Ai1 (z −1 )y1 (t) + . . . + [1 − Aii (z −1 )]yi (t)
(6)
+ ... + Bim (z −1 )um (t) + ei (t).

−k
−k+1
−1
aii
+ aii
+ . . . + aii
+
1z
2z
kz
ij −k
ij
ij −1
−2
a1 z + . . . + a(k−1) z + ak z
−k
−2
−1
bij
+ . . . + bij
+ bij
1 z
kz
(k−1) z
max{n∗a , n∗b },

t−l−k
p(yi (t − l) Yt−l−1
, w, λ)

1
1
exp{ ky − Φwk22 }
=
2λ
(2πλ)(t−k)/2

We parameterize each node of multivariable ARX model
structure in the same form and set a corresponding upper
bound to the polynomial orders of ground truth as:

−1

t−k−1
Y

Note that this likelihood distribution is not presented in
a standard form of Gaussian distribution, which however
shows its logarithm is a quadratic function of w. By maximizing its logarithm with respect to w, we end with the
PEM (ML) method. With infinite data points, PEM is
guaranteed to converge to the ground truth model [11].
In practice, given limited data, PEM may suffer from
over-fitting due to the upper bound of polynomial orders
and also generate a full connected network even if the
true one is sparse. Therefore, penalties for both network
topology and model complexity are essential. Referring
to the parameterizations in eq. (8), sparse network can
be interpreted as group sparse w whereas sparsity within
each group indicates reduced order of polynomials. This
consideration naturally leads us to the usage of Maximum a-posteriori method (MAP). A direct framework
to achieve those two levels of sparsity is Sparse Group
Lasso (SGL) as a Type I MAP method. In this paper, we
resort to Sparse Bayesian Learning as a Type II MAP

1
(7)

where a and b denote the parameters in polynomial matrices A(z −1 ) and B(z −1 ) respectively, superscript ij the
polynomial of ijth element of matrices and subscript i
the index of ith coefficient of polynomial.

Assume that time-series data collected from discrete
time indices 1 to t for each node and input are available.
For the ith node, we define following matrices and vec-

3

method. It employs inseparable prior distribution which
generalizes the penalty of Type I. It is also shown to be
superior than most Type I MAP schemes. In addition, we
show later that our approach is closely related to SGL.

Remark 4 : If the true model is ARX type, then the
last group of w which presents auto-regression can be excluded from the group sparsity to improve the estimation
accuracy. If not, the resultant model can be an FIR OE
model.

4

RECONSTRUCTION VIA SPARSE BAYESIAN To promote both element and group sparsity, we multiLEARNING
ply (10) by (11) and normalize it to get a proper distribution:
4.1 Sparsity inducing priors
p(w) = C ∗ max N (w|0, B)N (w|0, Γ)ϕ(β)ϕ(γ)
γ≥0,β≥0

(12)
where C is the normalization constant and can be absorbed by functions ϕ(β) or ϕ(γ). Hence, we get a lower
bound of this prior as:

Full Bayesian treatment requires to introduce a prior
distribution for w. We define ha distributioni p(w) in a
P
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 | [25]. Estimating w using its posteriori mean is intractable in this
case because the posterior distribution p(w|y) is nonGaussian and not analytical. To simplify the problem,
sparse Bayesian learning approximates p(w|y) with a
Gaussian distribution so that the solution, E(w|y) can
be easily calculated . To do that, first note that 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 [12, 21].

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

If either βij or γi approaches 0, the corresponding Gaussian distribution becomes a Dirac delta function imposing element or group sparsity to w.
Remark 5 : The conventional way to promote both element and group sparsity is to use hierarchical Bayesian
by introducing two hyperparameters where one is regulated by the other. The main problem is how to choose
corresponding hyper-distributions, which is non-trivial.
Also, the hyperparameter which is deeper in hierarchy has
less impact on the inference procedure [8]. That means
the resultant penalty is probably not strong enough to impose group sparsity to the weighting vector w. As such,
multiplying two priors makes sense since both hyperparameters influence w directly. Nevertheless, it is shown
later that this scheme can be transferred into an equivalent hierarchical Bayesian format.

To induce element sparsity to w, we introduce a prior
p(w) as [14, 25]:
p(w) =

p+m
Y

p(wi ) =

i=1
k
Y

j=1

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

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

(13)

4.2

(10)

βi ≥0

Type II Maximization

Although the implicit prior p̂(w) is improper, we can still
get a normalized posteriori distribution of w as:

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

p̂(w|y) = R

where subscript i denotes ith group in a vector and j
the jth element in that group. β is a vector of hyperparameters which controls element sparsity of the vector
w. B is the covariance matrix of Gaussian distribution
and parameterized by vector β as B = diag{β}.

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

(14)

Clearly, p̂(w|y) is Gaussian since log(p(w|y)) is a
quadratic function of w:
p̂(w|y) = N (w|µ, Σ),

For group sparsity, the corresponding prior is [15, 27]:

(15)

where
p(w) =

p+m
Y
i=1

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


−1
Σ = (Γ−1 + B −1 ) + λ−1 ΦT Φ

−1 T
µ = λ−1 (Γ−1 + B −1 ) + λ−1 ΦT Φ
Φ y.

(11)

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

where γ is a vector of hyper-parameters which controls
group sparsity of w.

(16)

Now, the core of the problem is to choose proper hyperparameters β and γ so that p̂(w|y) is close to the real

4

Solve L1 as a DCP problem

4.3.1

p(w|y) under the selected criterion. 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 [2, 14, 25]. It is equivalent to estimating hyperparameters using the maximum likelihood method:
Z
(γ ∗ , β ∗ , λ∗ ) = arg min
p(y|w)|p(w) − p̂(w)|dw
β,γ,λ≥0
Z
= arg min −2 log p(y|w)p̂(w)dw

Let
u(w, λ, β, γ) = λ−1 ky − Φwk22 + wT (Γ−1 + B −1 )w

v(λ, β, γ) = − log |B + Γ| − log|λI + Φ(Γ−1 + B −1 )−1 ΦT |,
(19)
then the optimization problem L1 becomes:
min

β,γ,λ≥0

γ≥0,β≥0,λ≥0

= arg min −2 log p̂(y|β, γ, λ)

u(w, λ, β, γ) − v(λ, β, γ)

(20)

β,γ,λ≥0

(17)

Proposition 8 : Functions u(w, λ, β, γ) and v(λ, β, γ)
are both jointly convex with respect to their own variables.

Proposition 6 : Following (17), w, hyperparameters
and noise covariance, λ can be estimated by solving the
optimization problem:

proof: The derivation can be found in Section B of Appendix.

L1 :

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 convex-concave procedure (CCCP) which belongs
to majorization-minimization (MM) algorithm using the
linear majorization function [6,18]. For minx f (x) where
f (x) = u(x) − v(x) and u(x), v(x) are convex, we can
solve it iteratively by:

min λ−1 ky − Φwk22 + wT (Γ−1 + B −1 )w

β,γ,λ,w

+ log |B + Γ| + log|λI + Φ(Γ−1 + B −1 )−1 ΦT |
Subject to:
β ≥ 0, γ ≥ 0, λ ≥ 0
Or:
L2 :

−1
min y T λI + Φ(Γ−1 + B −1 )−1 ΦT
y

(18)

β,γ,λ

xn+1 = arg min u(x)− < x, ∇v(xn ) >,

+ log |B + Γ| + log|λI + Φ(Γ−1 + B −1 )−1 ΦT |
Subject to:
β ≥ 0, γ ≥ 0, λ ≥ 0

x

where < ·, · > denotes inner product.

where we set log ϕ(·) = const [25].

Therefore, the optimization problem L1 can be decomposed into sequential convex optimization problems:

The derivation can be found in Section A of Appendix.
Remark 7 : The optimization problem (18) can also
be derived using hierarchical Bayesian format. Assume
3
p(w|L) = N (w|0, L), p(L|P ) = ( 23 )p+m |P | 2 |P − L|1/2
5
and p(P ) ∝ |P |− 2 where L andRP are counterparts of B
and Γ respectively, then p(w) ∝ p(w|L)p(L|P )p(P )dLdP .
In SBL, normalized p(w|y, L, P ) is used to approximate
p(w|y). The value of L and P is evaluated by solving
MAP of marginal distribution of hyperparameters as
max p(L, P |y). After simple manipulation, the resultant
optimization problem is same with (18). In this hierarchical format, L takes charge of element sparsity of w.
p(L|P ) brings an underlying constraint 0  L  Γ which
controls sparsity of L. p(P ) is similar to a Jeffrey prior
thus enforcing group sparsity to L. Therefore, hyperparameter P imposes group sparsity to w indirectly via L.
4.3

(21)

Algorithm to solve Type II maximization

 n+1 n+1 n+1 n+1 
γ
,β
,λ
,w

= arg min u(w, λ, β, γ) −
γ,β,λ,w

− ∇Tβ v

[λn ,β n ,γ n ]

subject to:

β − ∇Tγ v

∂v
∂λ

[λn ,β n ,γ n ]

[λn ,β n ,γ n ]

5

(22)

β ≥ 0, γ ≥ 0, λ ≥ 0
where
∂v
= trace{∆}
∂λ
−∇β v = (γ + β)−1 + diag{ΦT ∆Φ}.∗ [γ 2 .∗ (γ + β)−2 ]
−

−

k
2
X
βij
(ΦT ∆Φ)qq
∂v
1
=
+
∂γi
γ + βij
(βij + γi )2
j=1 i

∆ = [λI + Φ(Γ−1 + B −1 )−1 ΦT ]−1
q = (i − 1)k + j.

Since optimization problems L1 and L2 are equivalent,
we begin with the solution to L1 .

γ

λ

(23)

Algorithm 1 Solve L1 using CCCP
1: Initialize β 0 , γ 0 , λ0
β
γ
2: Calculate gij , gi and g λ using (23) and (25)
3: for n = 1 : M ax do
4:
Solve the reweighted convex problem:

If we optimize γ, β and λ first, we get analytical expressions of their optimal solutions as functions of w:
|wij |
opt
βij
= q , γiopt
β
gij

where

kwi k2
ky − Φwk2
p
= p γ , λopt =
.
gi
gλ

(24)

wn+1 = arg min
w

β
= −∇ij
gij
β v|[λn ,β n ,γ n ]

giγ = −∇iγ v|[λn ,β n ,γ n ]
∂v
g λ = − |[λn ,β n ,γ n ] .
∂λ

p+m
k
q
q
XX
p
1
β
g λ ky − Φwk2 +
giγ kwi k2 + gij
|wij |
k
i=1 j=1

(25)
5:
6:

Therefore, algorithm (22) can be further simplified by
substituting (24):

7:
8:
9:
10:
11:

wn+1 = arg min
w

p

g λ ky

− Φwk2 +

p+m
k
XX
i=1 j=1

1
k

q

giγ kwi k2 +

q
β
|wij |.
gij

(27)
Update β n+1 , γ n+1 and λn+1 according to (24)
Prune out small βij and γi and remove corresponding columns of dictionary matrix Φ
β
Update gij
, giγ and g λ using (23) and (25)
if Any stopping criteria is satisfied then
Break
end if
end for

can be treated as a sharing problem and can be solved
using ADMM algorithms [3, 13].

(26)

The optimization problem (26) can be solved as a Second Order Cone Program (SOCP). It has a very similar form with SGL which blends `1 and `2 -norm as the
penalty [17]. The differences are an extra weight imposed
to the first data-dependent term with 2-norm instead of
`2 -norm thus leading to a very different function property. Such variation is due to the estimation of λ incorporated into the problem. If λ is known, then the resultant
optimization problem is same with standard SGL.

The optimization problem (26) is firstly transferred into
a sharing problem:
win+1 = arg min
p
g λ ky −

w
p+m
X
i=1

zi k2 +

subject to

i=1 j=1

1
k

q

giγ kwi k2 +

q
β
gij
|wij |

φi wi − zi = 0.

A special attention should be paid to the estimation of
the noise variance λ, which is already mentioned in [20,
24]. Under some circumstance, λ can be compensated by
two hyperparameters thus unidentifiable. That is mainly
caused by limited data points. A simple interpretation is
law of large numbers. A large number of data is essential
to fully explore the property of a distribution. Therefore,
one should not expect to get an accurate estimation of
the noise variance but rather treat it as a tuning parameter of the optimization problem. Under this consideration, other reasonable criteria to explore the value of λ
are recommended such as cross validation.

where

(28)

h
i
Φ = φ1 . . . φp+m ,

(29)

φi ∈ RM m×k .

Such a problem can be solved using a scaled ADMM
method [3]:
win+1 := arg min
wi

To summarize, the above procedure can be described as
Algorithm 1.

k q
q
X
β
giγ kwi k2 +
gij
|wij |
j=1

n

z n+1
4.3.2

p+m
k
XX

Solve CCCP by Alternating Direction Method of
Multipliers (ADMM)

+ (ρ/2)kφi wi − φi win + Φw − z n + un k22
p
:= arg min g λ ky − (p + m)zk2
z

+ [(p + m)ρ/2] kz − un − Φw

un+1 :=uk + Φw

If the network is large leading to high dimensional variables, solving (26) will be very slow due to hardware limitation. In this case, we split the optimization problem
into a series of independent sub-problems to reduce computation burden. It turns out that the split cost function

n+1

− z n+1 ,

where
n

Φw = 1/(p + m)

(30)
p+m
X
i=1

6

n+1 2
k2

φi win .

(31)

Algorithm 2 Solve L2 using EM
1: Initialize β 0 , γ 0 , λ0
2: for n = 1 : M ax do
3:
E step: Formulate p(w|y, β n , γ n , λn ) according to
(15) and (16)
4:
M step:

w-update is a standard SGL problem which can be efficiently solved using accelerated generalized gradient descent algorithm [17]. z- update is a Group Lasso proby
, the
lem and can be solved analytically. Let ẑ = z − p+m
original z update becomes:
1
y
n+1 2
ẑ n+1 := arg min kẑ +
− un − Φw
k2
ẑ 2
p+m
p
(32)
gλ
kẑk2 .
+ √
ρ p+m
so that
ẑ n+1 =
where



c
kck2 − σ
kck2



[β n+1 , γ n+1 , λn+1 ]
= arg min Ew|β n ,γ n ,λn {ln p(y, w|β, γ, λ)} .

γin+1 =

+

y
n+1
+ un + Φw
p+m
p
(34)
gλ
σ= √
ρ p+m
Note that w-update can be solved in parallel. z and uupdate are then solved in sequence after collecting wupdate.

λ

6:

Solve L2 using EM method

7:
8:
9:
10:

Now, to solve L2 , first note that it has a similar form with
a traditional SBL problem [20, 23] except an extra term
log|B + Γ| which plays a key role in our formulation. If
we remove this term, then by replacing [(B −1 + Γ−1 )−1 ]
with a new variable the resultant problem is exactly same
with SBL. The most common way to solve SBL are EM
methods by treating w as the hidden variable [2,24]. Following the basic procedure of EM method, the algorithm
is described in Algorithm 2.

n+1

=

ky − Φµn k22 + λn

Pp+m Pk
i=1

N

j=1

1 − τij Σnqq

,

(36)
n −1
where τij = (βij
) + (γin )−1 ,q = (i − 1)k + j and
N = k(p + m)
Prune out small βij and γi and remove corresponding columns of dictionary matrix Φ
if Any stopping criteria is satisfied then
Break;
end if
end for

Although EM method and CCCP both belong to the
class of majorization-minimization (MM) method and
are special cases of DCA(Difference of Convex functions
Algorithm), a DC function has infinite many DC decompositions which can greatly influence the performance
of the algorithm [19]. To solve a standard linear regression problem with the dictionary matrix Φ ∈ RN ×M
and N  M , the cost of EM method is O(M N 2 ). However, if the dimension of w increases either because of the
scale of the network or the model complexity, ADMM
approach may be more advantageous.

L2 can also be solved directly by setting the gradient of
cost function 0 in order to find stationary points [20]:
−

k
1X n
Σ + (µnij )2
k j=1 qq

n+1
βij
= Σnqq + (µnij )2

c=−

4.3.3

Update solutions of M step as:

5:

(33)

(35)

1 2
1
=0
2 (µij + Σqq ) + β
βij
ij

k
k
1 X 2
(µij + Σqq ) +
=0
2
γi j=1
γi
Pp+m Pk
−1
−1
ky − Φµk22 + λ i=1
j=1 1 − (βij + γi )Σqq
λ=
N
(37)

−

5

Simple extension to nonlinear ARX model

Most networks in reality are nonlinear such as genetic
regulation networks. To cope with this fact, we introduce
nonlinear terms to the multivariable ARX model.

There is no analytical solution to equations (37). However, we can search for optima by cycling through variables in sequence. Applying simple manipulations to (37)
results in the same iteration in Algorithm 2.

For ith node:
y
yi (t) = fi (t) + Di1
(z −1 )y1 (t) + . . . + [1 − Ni (z −1 )]yi (t)
u
+ ... + Dim
(z −1 )um (t) + ei (t).
(38)

Remark 9 L2 can also be solved using CCCP. The first
term of the cost function is convex since it is partial minimum of u(w, β, γ). Clearly, it coincides with Algorithm 1.
The only difference comes from the sequence of optimizing variables.

We assume fi is the linear combination of nonlinear functions which form a dictionary or selected kernel func-

7

Table 1
Inference of the network topology of 100 experiment

tions [14].
fi [y(t), u(t)] =

q
X

αji fij [y(t), u(t)]

(39)

j=1

To select most relevant nonlinear terms, sparse priors
are introduced based on the property of those terms. For
instance, if a group of nonlinear terms describes the potential transcription activity of a transcriptional factor,
then group sparsity can be imposed. In addition, if it is
not sure wether a transcriptional factor is repressive or
active, then element sparsity is useful to select the right
hill functions.

FP(Max)

Correct

Our method

95%

0%

99%

SBL

95%

3%

83%

GSBL

91%

0%

38%

calculate the infinite norm of estimation error in Table 2:
kxerr k∞ = max |xi |
xerr = xest − xtrue .

(40)

Table 2
Infinite norm of estimation error of 100 experiment

Remark 10 : The concept of nonlinear ARX model can
be extended to a more general case. In biological system,
the ground truth model is always nonlinear. To identify
a nonlinear system directly is intractable in most cases
and approximating such systems using linear ones is only
valid within a small region of dynamics. We hence propose an intermediate model between complete nonlinear
and linear system so that it is more informative than linear models but still tractable at the same time. We propose
a state space model where all the hidden states (hidden
nodes) are described by a linear system while dynamics
of observable states (manifest nodes) also includes nonlinear dictionary functions whose variable only involves
observable states. Inference of a sparse network described
by such a state space model still needs further research.
However, one simplified version is nonlinear multivariable ARX model. All the observable states are nodes in
ARX model and nonlinear terms are only associated with
these nodes.
6

TP(Min)

Mean

Min

Max

Our method

2e-1

<1e-5

1.0

SBL

4e-1

<1e-5

4.9

GSBL

1.2

9.6e-1

2.2

Table 1 indicates that our method and GBL successfully exclude redundant correlation among nodes. Our
method is superior than the other two by ensuring 99%
correct inference of the network topology. GSBL is even
worse than SBL in lower rate of TP and total correct inference. Simulation shows that GSBL intends to search
for the correct network topology by setting the corresponding groups of parameters 0 but the resultant nonzero groups do not have sparse structure. Hence, the
polynomial order is not reduced if a big upper bound is
used. Also, GSBL is not able to explore the sparse profile within each group even if the ground truth polynomial coefficient is sparse and the lowest upper bound is
used. As a result, its estimation accuracy is lowest among
three approaches. In contrast to GSBL, SBL pays more
attention to finding accurate estimation of parameters
within each group. Its ability to infer the network topology and provide accurate parameter estimation seems
satisfactory. However, when applied to real data in practice such as biological time series, one does not expect
it to yield group sparse result, which in turn generates
full-connected networks. Our method has the best performance since it not only explores the network topology but also pursues the sparsity in each group to find
the lowest polynomial order possible. Although sometimes its parameter estimation is biased, the network
topology is inferred correctly with high probability. The
biased estimation may be due to two possibilities. The
major reason is the dictionary matrix correlates with the
noise in the model. That means some good properties
of SBL when dealing with linear regression problems no
long apply in this case [25]. Additional simulation using the same framework to solve standard linear regression problems confirms this concern since the algorithm
is capable to recover the weighting vector of a linear regression model. As such, deeper research in the performance of SBL when applied to system identification is

Simulation

We constructed ground truth ARX models with 10 nodes
and applied the technique to identify them. To compare
with other traditional methods, GSB and SBL were also
employed at the same time. The network topology of independent experiment was fixed while polynomial coefficients were generated randomly with maximum order
to be 3. We set the upper bound of polynomial order to
be 6. The exciting input was designed to be Gaussian
random noise. 20 time-series data points were collected
for each node and input in each experiment.
We first compare the inference of boolean structure of
networks using difference methods in Table 1. The first
column records minimum True Positive of multiple experiment, second column maximum False Positive and
third column rate of correct inference (100% TP with
0% FP).
To evaluate the estimation accuracy of parameters, we

8

required. The second reason is that both CCCP and EM
methods only guarantee convergence to a suboptimal solution. Simulation reveals that by changing the initial
point, some of the biased estimation can be fixed.
7

[10] D. Hayden, Y. Chang, J. Goncalves, and C.J. Tomlin. Sparse
network identifiability via compressed sensing. Automatica,
68:9–17, 2016.
[11] L. Ljung. System Identification: Theory for User. Prentice
Hall, 1998.
[12] J. Palmer, D.P. Wipf, K. Kreutz-Delgado, and B.D. Rao.
Variational em algorithms for non-gaussian latent variable
models. Adv. Neural Inf. Process. Syst., 18:1059, 2006.

CONCLUSION AND DISCUSSION

This paper proposes a method to identify the multivariable ARX model. The identification problem is formulated 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, including its topology or system order. The newly
proposed method penalises the complexity of the network by introducing both group and element sparsity.
These two kinds of sparse priors are combined and its
corresponding optimization problem is closely related to
Type I Sparse Group Lasso method.

[13] W. Pan, A. Sootla, and G.B. Stan. Distributed reconstruction
of nonlinear networks: An admm approach. IFAC Proc. Vol,
19:3208–3213, 2014.
[14] W. Pan, Y. Yuan, J. Goncalves, and G. Stan. A sparse
bayesian approach to the identification of nonlinear statespace systems. IEEE Trans. Automat. Contr., 61(1):182–187,
2016.
[15] W. Pan, Y. Yuan, L. Ljung, J. Goncalves, and
G. Stan. Identifying biochemical reaction networks using
heterogeneous datasets. IEEE Conference on Decision and
Control, 2015.
[16] B. Sanandaji, T. Vincent, and M. Wakin. Exact topology
identification of large-scale interconnected dynamical systems
from compressive observations. Proceedings of the 2011
American Control Conference, pages 649–656, 2011.

Further development includes two aspects. The first is
to find out the theoretical guarantees of the algorithm
performance. The second question is how to extend this
framework to infer a more general network model describes by a dynamical structural function. The main obstacle here is that SBL normally demands to formulate
the problem into a linear regression format so that the
likelihood is a quadratic function of parameters. Without this, variational inference or approximate Bayesian
computation may be more appropriate.

[17] Noah. Simon, J. Friedman, T. Hastie, and R. Tibshirani. A
sparse-group lasso. Journal of Computational and Graphical
Statistics, 22:231–245, 2013.
[18] B.K. Sriperumbudur. On the convergence of the concaveconvex procedure. Nips., 1:1–9, 2009.
[19] H. Le Thi, T. Pham Dinh nad H. Le, and X. Vo. Dc
approximation approaches for sparse optimization. European
Journal of Operational Research, pages 26–46, 2015.
[20] M. Tipping. Sparse bayesian learning and the relevance
vector mach. J. Mach. Learn. Res., 1:211–244, 2001.

References

[21] M. Wainwright and M. Jordan.
Graphical models,
exponential families, and variational inference. Foundations
and Trends in Machine Learning, pages 1–305, 2008.

[1] S. Babacan, S. Nakajima, and M. Do. Bayesian group-sparse
modeling and variational inference. Signal Process. IEEE
Trans., 62(11):2906–2921, 2014.
[2] C. Bishop. Pattern recognition and machine learning.
Springer New York, 2006.
[3] S. Boyd. Distributed optimization and statistical learning
via the alternating direction method of multipliers. Found.
Trends Mach. Learn., 3(1):1–122, 2010.
[4] S. Boyd and L. Vandenberghe.
Convex Optimization.
Cambridge University Press, 2004.
[5] R. Chen, C. Chu, S. Yuan, and Y. Wu. Bayesian sparse group
selection. Journal of Computational and Graphical Statistics,
pages 1–29, 2015.
[6] 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.,
59(11):2933–2945, 2014.
[7] A. Chiuso and G. Pillonetto. A bayesian approach to sparse
dynamic network identification. Automatica, pages 1553–
1565, 2012.
[8] R. Giri and B.D. Rao. Type i and type ii bayesian methods
for sparse signal recovery using scale mixtures. IEEE
Transactions on Signal Processing, pages 3418–3428, 2016.
[9] J. Goncalves and S. Warnick. Necessary and sufficient
conditions for dynamical structure reconstruction of lti
networks. IEEE Trans. Automat. Contr., 53(7):1670–1674,
2008.

[22] H. Weerts, A. Dankers, and P. Hof. Identifiability in dynamic
network identification. IFAC-PapersOnLine, 48(28):1409–
1414, 2015.
[23] D.P. Wipf and B.D. Rao. Sparse bayesian learning for basis
selection. Signal Processing, IEEE Transactions on Signal
Processing, 52(8):2153–2164, 2004.
[24] D.P. Wipf and B.D. Rao. An empirical bayesian strategy
for solving the simultaneous sparse approximation problem.
IEEE Transactions on Signal Processing, 55(7):3704–3716,
2007.
[25] D.P. Wipf, B.D. Rao, and S. Nagarajan. Latent variable
bayesian models for promoting sparsity. IEEE Trans. Inf.
Theory., 57(9):6236–6255, 2011.
[26] Y. Yuan, G. Stan, S. Warnick, and J. Goncalves. Robust
dynamical networkstructure reconstruction. Automatica,
47:1230–1235, 2011.
[27] 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., 61(8):2009–2015, 2013.
[28] Y. Zhu and H. Hjalmarsson. The box-jenkins steiglitzmcbride algorithm. Automatica, 65:170–182, 2016.

9

A

Derivation of Type II maximization problem

Using (A.3) and (A.5),the integral in (A.1) becomes:
− 2 log

Firstly, note that:

− 2 log
− 2 log
− 2 log

Z

= Ey + log I + (Γ−1 + B −1 )−1 λ−1 ΦT Φ
+ log |(Γ + B)| + log |λI|

N (y|Φw, λI)N (w|0, B)N (w|0, Γ)ϕ(β)ϕ(γ)dw

Z

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

= Ey + log (Γ−1 + B −1 ) + λ−1 ΦT Φ
+ log |λI| + log|B| + log |Γ|

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

Z

Z

= Ey + log I + λ−1 Φ(Γ−1 + B −1 )−1 ΦT
+ log |(Γ + B)| + log |λI|

exp(Ew )dw + log |λI| + log |B| + log |Γ|dw.

= Ey + log λI + Φ(Γ−1 + B −1 )−1 ΦT + log |(Γ + B)|

(A.1)

= min λ−1 k(y − Φw)k22 + wT (Γ−1 + B −1 )w

where

w

+ log λI + Φ(Γ−1 + B −1 )−1 ΦT + log |(Γ + B)|.
1
Ew = − [λ−1 k(y − Φw)k22 + wT (B −1 + Γ−1 )w],
2
(A.2)
ignoring all the constant terms.

(A.6)

As a result, we get:
L1 :

min λ−1 k(y − Φw)k22 + wT (Γ−1 + B −1 )w

β,γ,λ,w

+ log |B + Γ| + log |λI + Φ(Γ−1 + B −1 )−1 ΦT |
Subject to:
β ≥ 0, γ ≥ 0, λ ≥ 0.
Or:
(A.7)
L2 :

−1
min y T λI + Φ(Γ−1 + B −1 )−1 ΦT
y

By completing the square:

1
1
k(y − Φw)k22 − wT (Γ−1 + B −1 )w
2λ
2

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

−

(A.3)

β,γ,λ

+ log |B + Γ| + log |λI + Φ(Γ−1 + B −1 )−1 ΦT |
Subject to:
β ≥ 0, γ ≥ 0, λ ≥ 0.

where

−1
Σ = (Γ−1 + B −1 ) + λ−1 ΦT Φ

−1 T
µ = λ−1 (Γ−1 + B −1 ) + λ−1 ΦT Φ
Φ y

Ey = min λ−1 k(y − Φw)k22 + wT (Γ−1 + B −1 )w
w

−1
= y T λI + Φ(Γ−1 + B −1 )−1 ΦT
y

B

To see functions u(w, λ, β, γ) and v(λ, β, γ) are jointly
convex functions, we need to prove each term is jointly
convex. To check the convexity of u(w, λ, β, γ), we consider the epigraph of its two terms defined as epif =
{(x, t)|x ∈ domf, f (x) ≤ t} [4]. It is known that:

(A.4)

λI  0, λ−1 k(y − Φw)k22 < t
equivalent to
"
#
λI
y − Φw
0
(y − Φw)T
t

In addition, note that:

1
T −1
− 2 log exp − [(w − µ) Σ (w − µ)] dw
2
(A.5)
= − log |Σ|
Z

Convexity of DCP problem



(B.1)

so the term λ−1 k(y − Φw)k22 is jointly convex as is same
with wT (Γ−1 + B −1 )w since their epigraph is convex set
described by LMI [4].

= log (Γ−1 + B −1 ) + λ−1 ΦT Φ

For the function v(λ, β, γ), firstly note that − log | · | is a
convex function in S + . Since B + Γ is an affine function
of Γ and B, − log |B +Γ| is jointly convex with respect to

ignoring all the constant terms.

10

tion gij (γ, β) is concave. For g 1 , it is an affine function. According to the lemma above, log|λI + Φ(Γ−1 +
B −1 )−1 ΦT | is jointly concave with respect to λ, γ and β.

Γ and B. The second term of the function v seems more
complex. To prove its convexity, we first consider the
following lemma. A similar lemma with lower dimension
of function domain can be found in [4]:
Lemma 11 If a function f (x) = h(g1 (x), ..., gk (x))
where h(z1 , ..., zk ) : Rk → R and gi : Rn → R, then
f is concave if h is concave and nondecreasing in each
argument and gi are concave.
The proof of this lemma is straightforward by checking
the Hessian matrix:
∇2 f = JgT ∇2 hJg +

k
X
∂h
∇2 gl ,
∂zk g

(B.2)

l=1

where Jg denotes Jacobian matrix of function g =
[g1 , ..., gk ]T .
Let h(x1 , x) = log|x1 I + Φdiag(x)ΦT | and g 1 (λ, β, γ) =
γ β
λ, gij (λ, β, γ) = γi i+βijij with x ∈ R(p+m)k , i ∈ [1, p + m]
and j ∈ [1, k], then h(g 1 , g11 , g12 , ..., g(p+m)k ) = log |λI+
Φ(Γ−1 + B −1 )−1 ΦT |. Obviously, h(·) is jointly concave
with respect to x1 and x.
Since the Hessian matrix of g 1 is 0, we only need to check
the gradient of h(·) with respect to x:


∂h
= trace ΦT (x1 I + Φdiag(x)ΦT )−1 Φdiag(ei )
∂xi
= ∆ii ,
(B.3)
where ∆ = ΦT (x1 I +Φdiag(x)ΦT )−1 Φ and ei is a vector
with its ith element 1 and all the others 0.
Since matrix ∆ is at least semi-positive definite, its diagonal elements must be non-negative. As a result, h(x)
is nondecreasing in each argument. We finally calculate
the Hessian matrix, H of gij (β, γ). Note that matrix
Hqq =

∂ 2 gij
2 ,
∂βij

Hll =

∂ 2 gij
∂γi2

and Hql = Hlq =

∂ 2 gij
∂βij γi

with

all the other elements 0 where q = (i − 1)k + j + 1 and
l = (p + m)k + i + 1. It is always possible to find a permutation matrix such that:
"
#
Ĥ 0
T
P HP =
(B.4)
0 0
where


Ĥ = 

−

2γi2 (γi +βij ) 2γi βij (γi +βij )
(γi +βij )4
(γi +βij )4
2
2βij
(γi +βij )
?
− (γi +βij )4



.

(B.5)

Obviously, matrix Ĥ is semi-negative definite so that
H is also semi-negative definite, which indicates func-

11

