An empirical Bayes approach to identification of modules in
dynamic networks

arXiv:1705.01337v1 [] 3 May 2017

Niklas Everitt a, Giulio Bottegal a, Håkan Hjalmarsson a
a

ACCESS Linneaus Center, School of Electrical Engineering, KTH Royal Institute of Technology, Sweden

Abstract
We present a new method of identifying a specific module in a dynamic network, possibly with feedback loops. Assuming
known topology, we express the dynamics by an acyclic network composed of two blocks where the first block accounts for
the relation between the known reference signals and the input to the target module, while the second block contains the
target module. Using an empirical Bayes approach, we model the first block as a Gaussian vector with covariance matrix
(kernel) given by the recently introduced stable spline kernel. The parameters of the target module are estimated by solving
a marginal likelihood problem with a novel iterative scheme based on the Expectation-Maximization algorithm. Additionally,
we extend the method to include additional measurements downstream of the target module. Using Markov Chain Monte
Carlo techniques, it is shown that the same iterative scheme can solve also this formulation. Numerical experiments illustrate
the effectiveness of the proposed methods.
Key words: system identification, dynamic network, empirical Bayes, expectation-maximization.

1

Introduction

Networks of dynamical systems are everywhere, and
applications are in different branches of science, e.g.,
econometrics, systems biology, social science, and power
systems. Identification of these networks, usually referred to as dynamic networks, has been given increasing
attention in the system identification community, see
e.g., Materassi and Innocenti (2010), Van den Hof et al.
(2013), Hjalmarsson (2009).
In this paper, we use dynamic network to mean the interconnection of modules, where each module is a linear
time-invariant (LTI) system. The interconnecting signals are the outputs of these modules. In a graph interpretation, the interconnecting signals represent nodes
and the modules represent the edges of the graph. Moreover, we assume that exogenous measurable signals may
affect the dynamics of the network.
⋆ This work was supported by the Swedish Research Council
under contracts 2015-05285 and 2016-06079, and by the European Research Council under the advanced grant LEARN,
contract 267381.
Email addresses: neveritt@kth.se (Niklas Everitt),
bottegal@kth.se (Giulio Bottegal), hjalmars@kth.se
(Håkan Hjalmarsson).

Preprint submitted to Automatica

Two main problems arise in dynamic network identification. The first is unraveling the network topology (i.e.,
identify the edges of the graph), which can be seen as
a model structure selection problem. The second problems is the identification of one or more specific modules
in the network.
Some recent papers deal with both the aforementioned
problems (Materassi and Salapaka; 2012; Chiuso and Pillonetto;
2012; Materassi and Innocenti; 2010; Hayden et al.;
2014), whereas others are mainly focused on the
identification of a single module in the network
(Dankers et al.; 2013; Gunes et al.; 2014; Dankers et al.;
2015; Haber and Verhaegen; 2014; Torres et al.; 2014).
In particular, Dankers et al. (2013), and Van den Hof et al.
(2013) study the problem of understanding which of
the available output measurements should be used to
obtain a consistent estimate of a target module. In
Dankers et al. (2015) instead, errors-in-variables dynamic networks are considered, and methods that lead
to consistent module estimates are proposed. As observed in Van den Hof et al. (2013), dynamic networks
with known topology can be seen as a generalization
of simple compositions, such as systems in cascade,
series or feedback connection. Therefore, identification
techniques for dynamic networks may be derived by
extending methods already developed for simple structures. This is the idea underlying the method presented

4 May 2017

as the Expectation/Conditional-Maximization (ECM)
algorithm (Meng and Rubin; 1993), which alternates
the so called expectation step (E-step) with a series
of conditional-maximization steps (CM-steps). When
only the module input and output sensors are used, the
E-step admits an analytical expression, because joint
likelihood of the module output and the sensitivity function is Gaussian. As for the CM-steps, one has to solve
relatively simple optimization problems, which either
admit a closed form solution, or can be efficiently solved
using gradient descent strategies. Therefore, the overall
optimization scheme for solving the marginal likelihood
problem turns out computationally efficient.

in Van den Hof et al. (2013), which generalizes the twostage method, originally developed for closed-loop systems, to dynamic networks (Forssell and Ljung; 1999).
Instrumental variable methods for closed-loop systems
(Gilson and Van den Hof; 2005) are adapted to networks in Dankers et al. (2015). Similarly, the methodology proposed in Wahlberg et al. (2009) for the identification of cascaded systems is generalized to the context
of dynamic networks in Gunes et al. (2014). In that
work, the underlying idea is that a dynamic network
can be transformed into an acyclic structure, where any
reference signal of the network is the input to a cascaded
system consisting of two LTI blocks. In this alternative
system description, the first block captures the relation
between the reference and the noisy input of the target
module, the second block contains the target module.
The two LTI blocks are identified simultaneously using
the prediction error method (PEM) (Ljung; 1998). In
this setup, determining the model structure of the first
block of the cascaded structure may be complicated,
due to the possibly large number of interconnections in
the dynamic network. Furthermore, it requires knowledge of the model structure of essentially all modules
in the feedback loop. Therefore, in Gunes et al. (2014),
the first block is modeled by an unstructured finite impulse response (FIR) model of high order. The major
drawback of this approach is that, as is usually the case
with estimated models of high order, the variance of the
estimated FIR model is high. The uncertainty in the
estimate of the FIR model of the first block will in turn
decrease the accuracy of the estimated target module.

The second main contribution of the paper deals with
the case where more sensors spread in the network
are used in the identification of the target module.
Adding information through addition of measurements
used in the identification process has the potential to
further reduce the variance of the estimated module
(Everitt et al.; 2017). The downside is that an additional
measurement comes with another module to estimate,
also increasing the number of parameters to estimate.
To keep the number of additional parameters to estimate low, we propose a method that exploits regularization, modeling as a Gaussian process also the impulse
response of the path linking the target module to any
additional sensor. In this case, however, the measured
outputs and the unknown paths do not admit a joint
Gaussian description. As a consequence, the E-step of
the ECM method does not admit an analytical expression, as opposed to the one-sensor case described above.
To overcome this issue, we use Markov Chain Monte
Carlo (MCMC) techniques (Gilks et al.; 1995) to solve
the integral associated with the E-step. In particular, we
design an integration scheme based on the Gibbs sampler (Geman and Geman; 1984) that, in combination
with the ECM method, builds up a novel identification
method for the target module reminiscent of the so
called empirical Bayes Gibbs sampling (Casella; 2001).

The objective of this paper is to propose a method
for the identification of a module in dynamic networks
that circumvents the high variance that is due to the
high order model of the first block. The main contributions of this paper are two-fold. First, we discuss
the case where only the sensors directly measuring the
input and the output of the target module are used
in the identification process. Following a recent trend
in system identification, we use regularization to control the variance (Chen et al.; 2012). In particular, by
exploiting the equivalence between regularization and
Gaussian process regression (Pillonetto et al.; 2014),
we model the impulse response of the first block as a
zero-mean stochastic process. The covariance matrix
is given by the recently introduced first-order stable
spline kernel (Pillonetto and De Nicolao; 2010), whose
structure is parametrized by two hyperparameters. An
estimate of the target module is then obtained by empirical Bayes (EB) arguments, that is, by maximization
of the marginal likelihood of the available measurements (Pillonetto et al.; 2014). This likelihood depends
not only on the parameter of the target module, but
also on the kernel hyperparameters and the variance
of the measurement noise. Therefore, it is required
to estimate all these quantities. This is done by designing a novel iterative solution scheme based on an
EM-type algorithm (Dempster et al.; 1977), known

The effectiveness of the proposed methods is demonstrated through numerical experiments. The methods
proposed in this paper are close in spirit to some
recently proposed kernel-based techniques for blind
system identification (Bottegal et al.; 2015) and Hammerstein system identification (Risuleo et al.; 2015).
A part of this paper has previously been presented in
Everitt, Bottegal, Rojas and Hjalmarsson (2016). More
specifically, the case where only the sensors directly
measuring the input and the output of the target module are used in the identification process where partly
covered in Everitt, Bottegal, Rojas and Hjalmarsson
(2016), whereas, the method where more sensors spread
in the network are used in the identification of the target
module is completely novel.
The paper is organized as follows. In the next section,
we introduce the dynamic network model and we give

2

the problem statement. In Section 3 we present the identification strategy. In Section 4, we describe the solution scheme based on the ECM algorithm. Additional
measurements are added in Section 5, and we present
the MCMC based scheme to estimate the target module. Section 6 reports the results of Monte Carlo experiments. Some conclusions end the paper.

are unknown. To ensure stability and causality of the
network the following assumptions hold for all networks
considered in this paper.

1.1

Assumption 2.2 The sensitivity path S(q)

Assumption 2.1 The network is well posed in the sense
that all principal minors of limq→∞ (I − G(q)) are nonzero (Van den Hof et al.; 2013).

Notation

S(q) := (I − G(q))−1

Given a sequence of scalars {a(t)}m
t=1 , we denote by a
its vector representation a = [a(1) · · · a(m)]T ∈ ℜm .
Given a vector a ∈ ℜm , we define by Tn (a) the m ×
n lower triangular Toeplitz matrix whose elements are
the entries of a. Lower case letters indicate, in general,
column vectors and, when there is no confusion, capital
letters indicate their Toeplitz form, so given a ∈ ℜm , we
have that A = Tn (a), where the number n of columns is
consistent with the rest of the formula. The symbol “⊗”
denotes the standard Kronecker product of two matrices.

is stable.
Assumption 2.3 The reference variables {rl (t)} are
mutually uncorrelated and uncorrelated with the measurement noise e(t).
Thus, we can write
w̃(t) = S(q)r(t) + e(t) .

2

(2)

Problem Statement

2.1

We define a Nj as the set of indices of internal variables
that have a direct causal connection to wj , i.e., i ∈ Nj
if and only if Gji (q) 6= 0. Without loss of generality, we
assume that Nj = 1, 2, . . . , p, where p is the number of
direct causal connections to wj (we may always rename
the nodes so that this holds). The goal is to identify
module Gj1 (q) given N measurements of the reference
r(t), the “output” w̃j (t) and the set of p neighbor signals
in Nj . To this end, we express w̃j , the measured output
of module Gj1 (q) as

Dynamic networks

We consider dynamic networks that consist of L scalar
internal variables wj (t), j = 1, . . . , L and L scalar external reference signals rl (t), l = 1, . . . , L, that can
be manipulated by the user. Some of the reference
signal may not be present, i.e., they may be identically zero. Define R as the set of indices of reference
signals that are present. In the dynamic network, the
internal variables are considered nodes and transfer
functions are the edges. Introducing the vector notation
w(t) := [w1 (t) . . . wL (t)]T , r(t) := [r1 (t) . . . rL (t)]T ,
the dynamics of the network are defined by the equation
w(t) = G(q)w(t) + r(t) ,

w̃j (t) =

X

Gji (q)wi (t) + rj (t) + ej (t) .

(3)

i∈Nj

The above equation depends on the internal variables
wi (t), i ∈ Nj , which we we only have noisy measurement
of; these can be expressed as

(1)

where


0

G12 (q)



 G21 (q)
G(q) = 
 .
 ..

GL1 (q)

0
..

.

···

···
..
.
..
.

G1L (q)
..
.

w̃i (t) = wi (t) + ei (t) =



X

Sil (q)rl (t) + ei (t) .

(4)

l∈R




,

G(L−1)L (q)

GL(L−1) (q)
0

where Sil (q) is the transfer function path from reference
rl (t) to output w̃i (t). Together, (3) and (4) allow us to
express the relevant part of the network, possibly containing feedback loops, as a direct acyclic graph with
two blocks connected in cascade. Note that, in general,
the first block depends on all other blocks in the network. Therefore, accurate low order parametrization of
this block depends on global knowledge of the network.

where Gji (q) is a proper rational transfer function for
j = 1, . . . , L, i = 1, . . . , L. The internal variables w(t)
are measured with additive white noise, that is
w̃(t) = w(t) + e(t) ,

Example 2.1 As an example consider the network depicted in Figure 1, where, using (3) and (4), the acyclic
graph of Figure 2 can describe the relevant dynamics,
when wj = w3 is the output and we wish to identify
G31 (q).

where e(t) ∈ RL is a stationary zero-mean Gaussian
white-noise process
noise covariance ma with diagonal
2
trix Σe = diag σ12 , . . . , σL
. We assume that the σi2
3

r4
+

r2
w4

G14

+

w1

+

G21

r1

G31
w2

G12

+

G32

w3

Fig. 1. Network example of 4 internal variables and 2 reference signals.

S12 (q) S14 (q)
r4

S22 (q) S24 (q)

+
+

w1
w2

w1
+
G21 (q)

w2
+

Fig. 3. Basic network of 1 reference signal and 2 internal
variables.

(Gunes et al.; 2014). The main idea underlying SMPE
is that if, the two prediction errors (5) and (7) are simultaneously minimized, the variance will be decreased
(Wahlberg et al.; 2009). In the SMPE method, the prediction error of the measurement w̃j depends explicitly
on α and is given by
X
X
εj (t, θ, α) = w̃j (t) −
Gji (q, θ)
Sil (q, α)rl (t) . (8)

G23

G43

r2

S11 (q)

i∈Nj

G31 (q)
+ w3

l∈R

The method proceeds to minimize

G32 (q)



N
X ε2 (t, α)
1 X  ε2j (t, θ, α)
i
 . (9)
VN (θ, α) =
+
N t=1
σj2
σi2

Fig. 2. Direct acyclic graph of part of the network in Figure 1.

i∈Nj

In the following, we briefly review two standard methods for closed-loop identification that we will use as a
starting point to derive the methodology described in
the paper.
2.2

In (Gunes et al.; 2014), the noise variances are assumed
known, and how to estimate the noise variances is not
analyzed. As an initial estimate of the parameters θ and
α, the minimizers of the two-stage method can be taken.

A two stage method

The first stage of the two-stage method (Van den Hof et al.;
2013), proceeds by finding a consistent estimate ŵi (t)
of all nodes wi (t) in Nj . This is done by high-order
modeling of {Sil } and estimating it from (4) using
the prediction error method. The prediction errors are
constructed as
X
Sil (q, α)rl (t) ,
(5)
εi (t, α) = w̃i (t) −

The main drawback is that the least-squares estimation
of S may still induce high variance in the estimates.
Additionally, if each of the ns estimated transfer functions in S is estimated by the first n impulse response
coefficients, the number of estimated parameters in S
alone is ns · n. Already for relatively small dimensions
of S the SMPE method is prohibitively expensive. To
handle this, a frequency domain approach is taken in
Dankers and Van den Hof (2015). In this paper, we will
instead use regularization to reduce the variance and the
complexity.

where α is a parameter vector. The resulting estimate
Sil (q, α̂) is then used to obtain the node estimate as

3

l∈R

ŵi (t) =

X

Sil (q, α̂)rl (t) .

(6)

In this section we derive our approach to the identification of a specific module based on EB. For ease of exposition, we give a detailed derivation in the one-referenceone-module case. The extension to general dynamic networks follows along similar arguments.

l∈R

In a second stage, the module of interest Gj1 (q) (and
the other modules in Nj ) is parameterized by θ and estimated from (3), again using the prediction error method.
The prediction errors are now constructed as
X
εj (t, θ) = w̃j (t) − rj (t) −
Gji (q, θ)ŵi (t) . (7)

We consider a dynamic network with one non-zero reference signal r1 (t). Without loss of generality, we assume that the module of interest is G21 (q), and hence
G22 (q), . . . , G2L (q) are assumed zero (We can always rename the signals such that this holds). The setting we
consider has been illustrated in Figure 3. We parametrize
the target module by means of a parameter vector θ ∈
Rnθ . Using the vector notation introduced in the previous section, we denote by w̃1 the stacked measurements
w̃1 (t) before the module of interest G21 (q, θ), and by w̃2

i∈Nj

2.3

Empirical Bayes estimation of the module

Simultaneous minimization of prediction errors

It is useful to briefly introduce the simultaneous
minimization of prediction error method (SMPE)

4

the stacked output of this module w̃2 (t). We define the
impulse response coefficients of G21 (q, θ) by the inverse
discrete-time Fourier transform
Z π
1
gθ (t) :=
G21 (ejω , θ)ejωt dω .
(10)
2π −π

The parameter β regulates the decay velocity of the realizations from (15), whereas, λ tunes their amplitude.
In this context, Kβ is usually called a kernel (due to
the connection between Gaussian process regression and
the theory of reproducing kernel Hilbert space, see e.g.
Rasmussen and Williams (2006) for details) and determines the properties of the realizations of s. In particular, the stable spline kernel enforces smooth and BIBO
stable realizations (Pillonetto and De Nicolao; 2010).

Similarly we define s11 as the impulse response coefficients of S11 (q), where S11 (q) is, as before, the sensitivity path from r1 (t) to w1 (t), and e1 (t) and e2 (t) are
the measurement noise sources (which we have assumed
white and Gaussian). Their variance is denoted by σ12
and σ22 , respectively. We rewrite the dynamics as
w̃1 =
S11 r1 + e1 ,
w̃2 = Gθ S11 r1 + e2 ,

3.2

Since s11 is assumed stochastic, it admits a probabilistic description jointly with the vector of observations z,
parametrized by the vector

(11)

h
i
η = σ12 σ22 λ β θ .

where Gθ is the N ×N lower triangular Toeplitz matrix of
the N first impulse response samples gθ . The same notation holds for the impulse response s11 and its Toeplitzmatrix version S11 = TN (s11 ). We further rewrite (11)
as
w̃1 =
R1 s11 + e1 ,
w̃2 = Gθ R1 s11 + e2 .

z = Wθ s11 + e ,

Wθ =

R1
Gθ R1

#

(12)

e=

" #
e1
e2

Note that e is a random vector such that
#
"
 T
σ12 I 0
.
Σe := E ee =
0 σ22 I
3.1

p

"

z
s11

#

;η

!

∼N

#!
" # "
Σz Σzs
0
,
,
Σsz λKβ
0

(18)

T
where Σz = Wθ λKβ WθT + Σe , and Σzs = Σsz
=
Wθ λKβ . It is instrumental to derive the posterior distribution of s11 given the measurement vector z. It is
given by (Anderson and Moore; 1979)

p(s11 |z; η) ∼ N (P WθT Σe−1 z, P ) ,
P = (WθT Σe−1 Wθ + (λKβ )−1 )−1 ,

(13)

(19)
(20)

and it is also parametrized by the vector η.
The module identification strategy we propose in this
paper relies on an empirical Bayes approach. We introduce the marginal probability density function (pdf) of
the measurements
Z
p(z; η) = p(z, s11 ) ds11 ∼ N (0, Σz ) ,
(21)

(14)

Bayesian model of the sensitivity path

To reduce the variance in the sensitivity estimate (and
also reduce the number of estimated parameters), we
cast our problem in a Bayesian framework and model the
sensitivity function as a zero-mean Gaussian stochastic
vector (Rasmussen and Williams; 2006), i.e.,
p(s11 ; λ, Kβ ) ∼ N (0, λKβ ) .

(17)

In particular, having assumed a Gaussian distribution of
the noise, the joint description is also Gaussian, that is,

where R1 = TN (r1 ). For computational purposes, we
only consider the first n samples of s11 , where n is large
enough such that the truncation captures the dynamics
of the sensitivity S11 (q) well enough. Let z := [w̃1T w̃2T ]T ;
we rewrite (12) as
"

The marginal likelihood estimator

that is, the pdf of the measurements after having integrating out the dependence on the sensitivity path s11 .
Then, we can define the (log) marginal likelihood (ML)
criterion as the maximum of the marginal pdf defined
above

(15)

η̂ = arg max p(z; η)

(22)

η

The structure of the covariance matrix is given by the
first-order stable spline kernel (Pillonetto and De Nicolao;
2010):
{Kβ }i,j = β max(i,j) ,

β ∈ [0, 1) .

= arg min log det Σz + z T Σz−1 z ,
η



whose solution provides also an estimate of θ and thus
of the module of interest.

(16)

5

4

Proposition 4.1 Define

Computation of the solution of the marginal
likelihood criterion

o
n
(k)
Qβ (β) = log det Kβ + n log Tr Kβ−1 Ŝ11 .

Problem (22) is nonlinear and may involve a large number of decision variables, if nθ is large. In this section,
we derive an iterative solution scheme based on the
Expectation/Conditional-Maximization (ECM) algorithm (Meng and Rubin; 1993), which is a generalization of the standard Expectation-Maximization (EM)
algorithm. In order to employ EM-type algorithms, one
has to define a latent variable; in our problem, a natural
choice is s11 . Then, a (local) solution to (22) is achieved
by iterating over the following steps:

Then
β̂ (k+1) = arg min Qβ (β) ,
λ̂(k+1) =

o
n
1
(k)
.
Tr Kβ̂−1
Ŝ
(k+1) 11
n

(30)

Therefore, the update of the scaling hyperparameter is
available in closed-form, while the update of β requires
the solution of a scalar optimization problem in the domain [0, 1], an operation that requires little computational effort, see Bottegal et al. (2016) for details.

(23)

We are left with the maximization of the function
(k)
Q0 (σ12 , σ22 , θ). In order to simplify this step, we split
the optimization problem into constrained subproblems
that involve fewer decision variables. This operation
is justified by the ECM paradigm, which, under mild
conditions (Meng and Rubin; 1993), guarantees the
same convergence properties of the EM algorithm even
when the optimization of Q(k) (η) is split into a series of
constrained subproblems. In our case, we decouple the
update of the noise variances from the update of θ. By
means of the ECM paradigm, we split the maximiza(k)
tion of Q0 (σ12 , σ22 , θ) in a sequence of two constrained
optimization subproblems:

where the expectation is taken with respect to the
posterior of s11 when the estimate η (k) is used, i.e.,
p(s11 |z, η̂ (k) ) ;
(M-step) Solve the problem
η̂ (k+1) = arg max Q(k) (η) .

(29)

β∈[0,1)

(E-step) Given an estimate η̂ (k) (computed at the k-th
iteration of the algorithm), compute
Q(k) (η) := E [log p(z, s11 ; η)] ,

(28)

(24)

η

First, we turn our attention on the computation of the E(k)
step, i.e., the derivation of (23). Let ŝ11 and P̂ (k) be the
posterior mean and covariance matrix of s11 , computed
(k) (k)T
(k)
from (19) using η̂ (k) . Define Ŝ11 := P̂ (k) + ŝ11 ŝ11 .
The following proposition provides an expression for the
function Q(k) (η).

(k)

θ̂(k+1) = arg max Q0 (σ12 , σ22 , θ)

(31)

θ

2(k)

2(k)

s.t. σ12 = σ̂1

2(k)

Lemma 4.1 Let η̂ (k) = [σ̂1 σ̂2 λ̂(k) β̂ (k) θ̂(k) ] be an
estimate of η after the k-th iteration of the EM method.
Then
1 (k)
1
Q(k) (η) = − Q0 (σ12 , σ22 , θ) − Q(k)
(λ, β) ,
2
2 s

2(k+1)
σ̂1
,

2(k+1)
σ̂2

=

2(k)

, σ22 = σ̂2

(k)
arg max Q0 (σ12 ,
2
2
σ1 , σ2

σ12 , θ)

,
(32)

s.t. θ = θ̂(k+1) .

(25)

The following result provides the solution of the above
problems.

where

2


(k)
(k)
Q0 (σ12 , σ22 , θ) = log det Σe + z T Σe−1 z − 2z T Wθ ŝ11
o
n
(k)
,
(26)
+ Tr WθT Σe−1 Wθ Ŝ11
n
o
(k)
−1
Qs(k) (λ, β) = log det λKβ + Tr (λKβ ) Ŝ11 .

Proposition 4.2 Introduce the matrix D ∈ RN ×N
such that Da = vec(TN (a)), for any a ∈ RN . Define
(k)

Â(k) = DT (R1 Ŝ11 R1T ⊗ IN )D

(33)

(k)
TN (R1 ŝ11 )T w̃2

(34)

b̂

(27)

(k)

=

.

Then
Having computed the function Q(k) (η), we now focus on
its maximization. We first note that the decomposition
(25) shows that the kernel hyperparameters can be updated independently of the rest of the parameters:

θ̂(k+1) = arg min gθT Â(k) gθ − 2b̂(k)T gθ .

(35)

θ

The closed form updates of the noise variances are as

6

the module of interest, and modeled in time domain as
gθ1 , . . . , gθp . For any i = 1, . . . , p, we can write

follows
2(k+1)

σ̂1

2(k+1)
σ̂2

o
n
1
(k)
,
kw̃1 − R1 ŝ11 k22 + Tr R1 P̂ (k) R1T
N
1
(k)
kw̃2 − Gθ̂(k+1) R1 ŝ11 k22
=
N
o
n
. (36)
+ Tr Gθ̂(k+1) R1 P̂ (k) R1T GTθ̂(k+1)
=

w̃i = Rl1 sil1 + . . . + Rlm silm + eki = Rsi + eki , (38)
where R := [Rl1 . . . Rlm ] and si = [sTil1 . . . sTilm ]T . Using these definitions we can also write (cf. (3))
w̃j = rj + Gθ1 Rs1 + . . . + Gθp Rsp + ej .

Defining also w = [w̃1T . . . w̃pT ]T , s = [sT1 . . . sTp ]T , Gθ =
[Gθ1 . . . Gθp ], ew = [eT1 . . . eTp ]T , we obtain the following
expression for the network dynamics

Each variance is the result of the sum of one term that
measures the adherence of the identified systems to the
data and one term that compensates for the bias in the
estimates introduced by the Bayesian approach. The update of the parameter θ involves a (generally) nonlinear
least-squares problem, which can be solved using gradient descent strategies. Note that, in case the impulse response gθ is linearly parametrized (e.g., it is an FIR system or orthonormal basis functions are used (Wahlberg;
1991)), then the update of θ is also available in closedform.

w = (Ip ⊗ R)s + ew
w̃j − rj = Gθ (Ip ⊗ R)s + ej ,

4.1

z = Wθ s + e, Wθ =

"

(Ip ⊗ R)

#

Gθ (Ip ⊗ R)

, e=

" #
ew
ej

. (41)

Each sensitivity path sil is given a prior of the form
(15), with hyperparameters λil and βil , assuming mutual
independence between the sensitivity paths. Although
it may appear more sensible to incorporate some correlation among the sensitivity paths, at present, it is
not clear how this can be done using Gaussian priors.
Some recent work suggests to enrich the stable spline
kernel with a component enforcing low McMillan degree
(Prando et al.; 2014). Furthermore, as we will see, assuming independent priors allows the kernel hyperparameters to be updated independently. Introducing Λ
as the diagonal matrix with elements corresponding to
{λil }, and similarly, defining Kβ with diagonal elements
{Kβil }, we have

(37)

Identification algorithm

The proposed method for module identification can be
summarized in the following steps.
(1) Find an initial estimate of η̂ (0) , set k = 0.
(k)
(2) Compute ŝ11 and P̂ (k) from (19).
(3) Update the kernel hyperparameters using (30),
(29).
(4) Update the vector θ solving (35).
(5) Update the noise variances from (36).
(6) Check if the algorithm has converged. If not, set
k = k + 1 and go back to step 2.

p(s; Λ, Kβ ) ∼ N (0, (Λ ⊗ In )Kβ ) .

(42)

We collect all the parameters characterizing the model
into the vector η. It follows that
p(z; η) ∼ N (0, Σz ) ,

The method can be initialized in several ways. One option is to first estimate Ŝ11 (q) by an empirical Bayes
method using only r1 and w̃1 . Then, ŵ1 is constructed
from (6), using the obtained Ŝ11 (q). Finally, G is estimated using the prediction error method, using ŵ1 as
input and w̃2 as output.
4.2

(40)

or, with z = [wT (w̃j − rj )T ]T

Example 4.1 Assume that the linear parametrization
gθ = Lθ, L ∈ RN ×nθ , is used, then

−1
LT b̂(k) .
θ̂(k+1) = LT Â(k) L

(39)

(43)

where Σz = Wθ (Λ ⊗ In )Kβ WθT + Σe , and

Σe = diag σ12 , . . . , σp2 , σj2 ⊗ IN .

(44)

η̂ = arg max log p(z; η) .

(45)

Therefore, we can define the following ML criterion
η

Extension to general structures

In this section, we generalize the previous algorithm
to a general network structure with m ≤ L reference signals {rl1 (t), . . . , rlm (t)}, and p ≤ L modules
{Gj1 (q), . . . , Gjp (q)} sharing the same output w̃j (t) as

Having set the notation, we outline the ECM algorithm
for this general setting below. To this end, note that
p(s|z; η) = N (ŝ, P) ,

7

(46)

where

r1

ŝ = PWθT Σe−1 z ,
−1

P = WθT Σe−1 Wθ + ((Λ ⊗ In )Kβ )−1
.

S11 (q)

w1
+
G21 (q)

w2
+

F (q)

w3
+

(47)
Fig. 4. Basic network of 1 reference signal and 3 internal
variables.

(48)

We use again the notation ŝ(k) , P̂(k) , Ŝ(k) to mean the
estimates of the corresponding quantities at iteration k.

only a few extra hyperparameters need to be estimated
and the dimensions of (35) remain the same in the update of θ.

Proposition 4.3 Let η collect all the parameters characterizing (45), and let η (k) be its estimate after the k-th
iteration of the ECM method. Then the estimate η (k+1)
is obtained by means of the following updates.

As reference signals can be added with little effort, a
natural question is if also output measurements “downstream” of the module of interest can be added with little effort. In Example 2.1 the measurement w4 is such
a measurement that, with the same strategy as before,
can be expressed as

Hyperparameters: Define
n
o
(k)
Qβij (β) = log det Kβ +n log Tr Kβ−1 Ŝij
, (49)
(k)

where Ŝij is the n×n diagonal block of Ŝ (k) corresponding
to the path sij . Then λij and βij are updated as in
Proposition 4.1, for any i, j.
Module parameters: Define


Â(k) := (Ip ⊗ D)T R̄Ŝ(k) R̄T ⊗ IN (Ip ⊗ D) , (50)
h
i
(k)
b̂(k)T := (w̃j − rj )T TN (Rŝ1 ) . . . TN (Rŝ(k)
p ) ,(51)
where D is as in Proposition 4.2 and R̄ = Ip ⊗ R.
Then
θ̂(k+1) = arg min gθT Â(k) gθ − 2b̂(k)T gθ ,

(52)

θ

where gθ := [gθT1 . . . gθTp ]T .
Noise variances:
o
n
1 
2(k+1)
(k)
(k)
=
σ̂i
kw̃i − Rŝi k22 + Tr RP̂i RT
N 
1
2(k+1)
kw̃j − rj − Gθ̂(k+1) (I ⊗ R)ŝ(k) k22
=
σ̂j
N n
o
,
(53)
+ Tr Gθ̂(k+1) R̄P̂(k) R̄T GTθ̂(k+1)

w4 (t) = G43 (q)w3 (t) + r4 (t) .

(54)

Using this measurement for the purpose of identification
would require the identification of G43 (q) in addition
to the previously considered modules. The signal w4 (t)
contains information about w3 (t), and thus information
about the module of interest. The price we have to pay
for this information is the additional parameters to estimate and, as we will see, another layer of complexity.
To extend the previous framework to include additional
measurements after the module of interest, let us consider the case where we would like to include only one additional measurement, in this context denoted by w̃3 (t);
the generalization to more sensors is straightforward but
notationally heavy. Let the path linking the target module to the additional sensor be denoted by F (q), with
impulse response f . Furthermore, let us for simplicity
consider the one-reference-signal-one-input case again,
i.e., (11), (12). The setting we consider has been illustrated in Figure 4. We model also this module using a
Bayesian framework by interpreting f as a zero-mean
Gaussian stochastic vector, i.e.,

(k)

where Pi is the nm × nm diagonal block of P(k) ,
(k)
corresponding to the covariance matrix of ŝi .
5

p(f ; λf , Kβf ) ∼ N (0, λf Kβf ) ,

Including additional sensors

(55)

where again Kβf is the first-order stable spline kernel
(16). We introduce the following variables

By using the kernel-based approach adopted above, the
sensitivity paths could be modeled with only a few hyperparameters while still keeping the module of interest
parametric. One potential benefit with this approach is
that including another reference signal will not increase
the number of estimated parameters significantly. Although the complexity of the problem increases slightly,

h
i
σ = σ12 σ22 σ32 ,
h
iT
z = w̃1T w̃2T w̃3T
,

zf = w̃3 .

8

(56)
(57)
(58)

For given vales of θ, s and f , we construct


R



Ws = 
 Gθ R  ,
F Gθ R
Wf = TN (Gθ Rs11 ) ,
Σ = diag {σ} ⊗ IN .

(M-step) Solve the problem
η̂ (k+1) = arg max Q(k) (η) .
(59)

While the M-Step remains substantially unchanged, the
E-step requires more attention. Now, it requires the computation of the integral

(60)
(61)

E [log p(z, s11 , f ; η)] =
(67)
Z
log p(z, s11 , f ; η)p(s11 , f |z, η̂ (k) ) ds11 df,

Notice that the last internal variable w3 can be expressed
as
w3 = F Gθ S11 r
= Gθ F S11 r
= Gθ F Rs11
= Gθ RF s11
= Gθ Rv,

which does not admit an analytical solution, because the
posterior distribution p(s11 , f |z, η̂ (k) ) is non-Gaussian
(it does not have an analytical form, in fact). However,
using Markov Chain Monte Carlo (MCMC) techniques,
we can compute an approximation of the integral by
sampling from the joint posterior density (also called a
target distribution)

(62)

where commutation of the matrices follows from the fact
that they are lower-triangular Toeplitz matrices, and
v := F s11 . For ease of exposition, we will also use the
notation v = f ∗ s11 .

p(s11 , f |z; η).

where λs , βs are the hyperparameters of the prior of s11 ,
that is

• p(s11 |f, z; η). Using (59), we write the linear model
z = Ws s11 + e,

(63)

(69)

where e = [eT1 eT2 eT3 ]T . Then, given f , the vectors
s11 and z are jointly Gaussian, so that

(64)

η

p(s|f, z; η) ∼ N (ms , Ps ) ,

does not admit an analytical expression, since the integral (64) is intractable. To treat this problem, again we
resort to EM-type methods. In this case, the latent variables to add to the problem are both s11 and f , so that
the EM method has to alternate between the following
steps.

(70)

with
Ps = WsT Σ −1 Ws + (λs Kβs )−1
ms = Ps WsT Σ −1 z .

(E-step) Given an estimate η̂ (k) (computed at the k-th
iteration of the algorithm), compute
Q(k) (η) := E [log p(z, s11 , f ; η)] ,

(68)

As pointed out before, (68) does not admit a closedform expression and hence direct sampling is a hard
task. However, if it is easy to draw samples from the
conditional probability distributions, samples of (68)
can be easily drawn using the Gibbs sampler. In Gibbs
sampling, each conditional is considered the state of a
Markov chain; by iteratively drawing samples from the
conditionals, the Markov chain will converge to its stationary distribution, which corresponds to the target distribution. In our problem, the conditionals of (68) are
as follows

The key difficulty in this setup is that the description
of the measurements and the system description with
both s11 and f no longer admit a jointly Gaussian probabilistic model, because v in (62) is the result of the convolution of two Gaussian vectors. In fact, a closed-form
expression is not available. This fact has a detrimental effect in our empirical Bayes approach, because the
marginal likelihood estimator of
h
i
η = σ λs βs λf βf θ ,

η̂ = arg max p(z; η)
η
Z
= arg max p(z, s11 , f ; η) ds11 df,

(66)

η

−1

• p(f |s, z; η). Given s and r, all sensors but the last
becomes redundant. Using (60) we write the linear
model
zf = Wf f + e3 ,
(71)

(65)

which shows that

where the expectation is taken with respect to the
target distribution when the estimate η (k) is used,
i.e., p(s11 , f |z, η̂ (k) ) ;

p(f |s11 , z; η) ∼ N (mf , Pf ),

9

(72)

with

Define
Pf =

WfT Wf
+ (λf Kβf )−1
σ32

mf = Pf

!−1

Q̃s (λ, β, x, X) := log det λKβ

+ Tr (λKβ )−1 (xxT + X) ,
1
2
Q̃z (σ 2 , z, x, X) := N log σ 2 + 2 kz − Rxk2
σ

1
+ 2 Tr RXRT ,
σ
1
2
2
Q̃f (σ , z, θ, x, X) := N log σ 2 + 2 kz − Gθ Rxk2
σ

1
+ 2 Tr Gθ RXRT GTθ .
σ

WfT
zf .
σ32

The following algorithm summarizes the Gibbs sampler
used for dynamic network identification.
Algorithm 1 Gibbs sampler for a dynamic network.
Initialization: compute initial value of s0 and f 0 . For
k = 1 to M + M0 :

Then

M
−2Q(k) (η) = lim Q̃s (λs , βs , sM
s , Ps ) ,
M→∞

(1) Draw the sample sk from p(s|f k−1 , z; η);
(2) Draw the sample f k from p(f |sk , z; η);

+Q̃s (λf , βf , fsM , PfM ) ,
M
+Q̃z (σ12 , w̃1 , sM
s , Ps ) ,
2
M
+Q̃f (σ2 , w̃2 , θ, ss , PsM ) ,
+Q̃f (σ32 , w̃3 , θ, vsM , PvM ) .

In this algorithm, M0 is the number of initial samples
that are discarded, which is also known as the burn-in
(Meyn and Tweedie; 2009). These samples are discarded
since the Markov chain needs a certain number of samples to converge to its stationary distribution.
5.1

The CM-steps are now very similar to the previous
method and follows by similar reasoning as in the proof
of Proposition 4.2.

The ECM method with additional sensor

Proposition 5.2 Let η̂ (k) be the parameter estimate obM T
tained at the k:th iteration. Define SsM = sM
s (ss ) +
M
M
M M T
M
Ps , Sv = vs (vs ) + Pv ,

We now discuss the computation of the E-step and the
CM-steps using the Gibbs sampler scheme introduce
above.

Âs
Âv
b̂s
b̂v

Proposition 5.1 Introduce the mean and covariance
quantities
sM
s =
fsM =
vsM =
PsM =
PfM

1
M
1
M
1
M
1
M

1
=
M

PvM =

1
M

MX
0 +M

sk ,

(73)

fk ,

(74)

vk ,

(75)

k
M T
(sk − sM
s )(s − ss ) ,

(76)

k=M0 +1
MX
0 +M
k=M0 +1
MX
0 +M
k=M0 +1
MX
0 +M
k=M0 +1
MX
0 +M

(f −

k=M0 +1
MX
0 +M

fsM )(f k

−

fsM )T

(v k − vsM )(v k − vsM )T ,

,

= DT (RSsM RT ⊗ IN )D ,
= DT (RSvM RT ⊗ IN )D ,
T
= TN (RsM
s ) w̃2 ,
= TN (RvsM )T w̃3 .

Then the updated parameter vector η̂ (k+1) is obtained as
follows


1
1
(k+1)
T
θ̂
= arg min gθ
Âs + 2 Âv gθ
σ22
σ
θ
 3
1 T
1
−2
(80)
b̂ + 2 b̂Tv gθ .
σ22 s
σ3
The closed form updates of the noise variances are
2(k+1)

k

(79)

σ̂1

(77)

2(k+1)

σ̂2
(78)

k=M0 +1

2(k+1)

σ̂3

where sk , f k and v k = sk ∗ f k are samples drawn using
Algorithm 1.

10



1
2
M T
,
kw̃1 − RsM
s k2 + Tr RPs R
N 
1
2
=
kw̃2 − Gθ̂(k+1) RsM
s k2
N
n
o
+ Tr Gθ̂(k+1) RPsM RT GTθ̂(k+1)
,

1
kw̃3 − Gθ̂(k+1) RvsM k22
=
N
o
n
. (81)
+ Tr Gθ̂(k+1) RPvM RT GTθ̂(k+1)
=

SMPE: The method is initialized by the two-stage
method, exactly as NEB. Then, the cost function (9),
with a slight modification, is minimized. The modification of the cost function comes from that, as mentioned
before, the SMPE method assumes that the noise variances are known. To make the comparison fair, also
the noise variances need to be estimated. By maximum
likelihood arguments, the logarithm of the determinant
of the complete noise covariance matrix is added to the
cost function (9) and the noise variances are included in
θ, the vector of parameters to estimate. The tolerance
is set to θ̂(k+1) − θ̂(k) / θ̂(k) < 10−10 .

The kernel hyperparameters are updated through (29) and
(30) for both s11 and f .
5.2

Identification algorithm

The proposed method for module identification can be
summarized in the following steps.
(1) Find an initial estimate of η̂ (0) , set k = 0.
(2) Compute the quantities (73)-(78) using Algorithm
1.
(3) Update the kernel hyperparameters using (30),
(29).
(4) Update the vector θ solving (80).
(5) Update the noise variances from (81).
(6) Check if the algorithm has converged. If not, set
k = k + 1 and go back to step 2.

The simulations were run in Julia, a high-level, highperformance dynamic programming language for technical computing (Bezanson et al.; 2017).
6.1

As can be seen, the main difference with the one-inputone-sensor algorithm (see Section 5.2) is that Step 2 of
the algorithm requires a heavier computational burden,
because of the integration via Gibbs sampling. Nevertheless, as will be seen in the next section, this pays off in
terms of performance in identifying the target module.
6

Closed-loop identification

The first Monte Carlo simulation is from a system operating in closed loop with an unknown low order controller
with N = 200 data samples. This setting is slightly different to the standard closed-loop setting in that the
measurement noise of w̃2 is not fed back in the loop, and
that the signals w1 and w2 are treated completely symmetric. The noise to signal ratio are all set to 1. The true
plant and true controller are chosen such that the sensitivity function S(q −1 ) has an impulse response that can
be well approximated by n = 100 impulse response coefficients. The closed loop is depicted in Figure 5, where

Numerical experiments

In this section, we present results from two Monte Carlo
simulations to illustrate the performance of the proposed method, which we abbreviate as Network Empirical Bayes (NEB) and its extension NEBX outlined in
Section 5, and we compare with SMPE (see Section 2.3).
We consider the network case of Example 2.1 and a simple closed loop network. The reference signals used are
zero-mean unit-variance Gaussian white noise. The noise
signals ek are zero-mean Gaussian white noise with variances such that noise to signal ratios Var wk /Var ek are
constant. The setting of the compared methods are provided in some more details below, where the model order
of the plant G(q) is known for both the SMPE method
and the proposed NEB method.

r1
+

w1

G

+

w2

C
Fig. 5. Closed loop network of first Monte Carlo simulation.

G(q, θ) =

b1 q −1 + b2 q −2
,
1 + a1 q −1 + a2 q −2

(82)
NEB: The method is initialized by the two-stage
method. First, Ŝ(q) is estimated by least-squares. SecThe controller C is given by
ond, G is estimated using MORSM (Everitt, Galrinho and Hjalmarsson;
2016) from the simulated signal ŵ obtained from (6)
0.8 + 0.4q −1 − 0.5q −2
and w̃j . MORSM is an iterative method that is asymp,
(83)
C(q, θ) =
1 + 0.5q −1 + 0.2q −2
totically efficient for open loop data. Then, the iterative
method outlined in Section 4.1 is employed with the
with the parameter vector θ = [b1 , b2 , a1 , a2 ], and true
stopping criteria η̂ (k+1) − η̂ (k) / η̂ (k) < 10−10 .
parameters θ0 = [0.4, 0.5, −0.4, 0.3].
NEBX: The method is initialized by NEB. f 0 is obtained
The two methods are compared using the fit of the imby an empirical Bayes method using simulated input and
pulse response coefficients of g according to
measured output of f . Then, the iterative method outlined in Section 5 is employed with the stopping criteria
g 0 − ĝ 2
η̂ (k+1) − η̂ (k) / η̂ (k) < 10−10 , or a maximum of 50
F IT = 1 −
(84)
iterations.
kg 0 k2

11

for a set of parameters that were chosen such that all
modules are stable and {S12 (q), S24 (q), S22 (q), S24 (q)}
are stable and can be well approximated with 70
impulse response coefficients. Two reference signals,
r2 (t) and r4 (t) are available and N = 200 data
samples are used with the goal to estimate G31 (q)
and G32 . In total 6 transfer functions are estimated,
{S12 (q), S24 (q), S22 (q), S24 (q), G31 (q) and G32 (q)},
where {S12 (q), S24 (q), S22 (q), S24 (q)} are each parameterized by n = 75 impulse response coefficients in all
methods. For NEBX also G43 (q) is estimated by n = 75
impulse response coefficients. The noise to signal ratio
at each measurement is set to Var wk /Var ek = 0.1 and
the additional measurement used in NEBX has a lower
noise to signal ratio of Var w4 /Var e4 = 0.01.

For this example, the proposed NEB method achieves a
higher fit, on average, than the SMPE method, cf. the
box plot of Figure 6. Comparing the fits obtained at
each Monte Carlo run (see Figure 7), it can be seen that
NEB consistently performs at least as good as SMPE
for almost every Monte Carlo run and in some runs considerably better. From the sample means and variance
reported in Table 1, it can be seen that, in general, the
estimates produced by NEB have smaller variance than
SMPE while their mean values are similar.
1

0.95

The fits of the impulse responses of G31 and G32 for the
experiment are shown as a boxplot in Figure 8 and Figure 10 respectively. Comparing the fits obtained at each
Monte Carlo run (see Figure 11 and Figure 11), the proposed NEB and NEBX methods are competitive with
the SMPE method for this network. In many cases, the
SMPE method failed to produce a reasonable estimate
as 10 percent of the Monte Carlo runs gave a negative
fit and were removed before the impulse response fits,
boxplots and parameter sample means and variances
were computed. From the sample means and variance reported in Table 2 and Table 3, it can be seen that, in general, the estimates produced by NEB and NEBX have,
in general, significantly smaller variance than SMPE,
while the mean values are roughly the same. Recalling
that one of the motivations of the proposed methods was
to reduced the variance induced by the high order modeling of the sensitivity paths, both the closed-loop example and network example gives some support for this
motivation.

0.9

0.85

SMPE

NEB

Fig. 6. Box plot of the fit of the impulse response of G
obtained by the SMPE, and NEB methods respectively.

FITNEB

1

0.95

In almost all of the Monte Carlo runs, NEBX outperformed NEB in this simulation. However, NEBX is significantly more computationally expensive than NEB.

NEB
0.85

0.9
0.95
FITSMPE

1

1

Fig. 7. Each fit of the impulse response coefficients of G for
NEB compared with SMPE for 100 Monte Carlo simulations.
The black line represents y = x, i.e., when SMPE performs
equally good as NEB. Note the scaling of the x-axis of this
figure.

6.2

0.95
0.9

Dynamic network example

0.85

This Monte Carlo simulation compares the NEB method
and NEBX with the SMPE method on data from the
network of Example 2.1, illustrated in Figure 1, where
each of the modules are of second order, i.e.,
Gij (q) =

0.8

b1 q −1 + b2 q −2
,
1 + a1 q −1 + a2 q −2

SMPE

NEB

NEBX

Fig. 8. Box plot of the fit of the impulse response of G31 obtained by the methods SMPE, NEB and NEBX respectively.

12

Table 1
Sample mean and sample variance of the parameters estimates for Ĝ for compared methods.
b01 = 0.2
b02 = 0.3
a01 = 0.4

a02 = 0.5

Method

Es b̂1

N · Vars b̂1

Es b̂2

N · Vars b̂2

Es â1

N · Vars â1

Es â2

N · Vars â2

SMPE

0.21

0.43

0.31

0.93

0.50

3.4

0.16

2.8

NEB

0.20

0.22

0.31

0.26

0.68

2.9

0.23

2.0

Table 2
Sample mean and sample variance of the parameters estimates for Ĝ31 for the three compared methods.
b01 = 0.2
b02 = 0.3
a01 = 0.4
a02 = 0.5
Method

Es b̂1

N · Vars b̂1

Es b̂2

N · Vars b̂2

Es â1

N · Vars â1

Es â2

N · Vars â2

SMPE

0.20

0.088

0.28

0.075

0.36

1.6

0.53

0.85

NEB

0.21

0.049

0.29

0.070

0.36

0.94

0.52

0.62

NEBX

0.20

0.024

0.29

0.036

0.40

0.60

0.50

0.52

Table 3
Sample mean and sample variance of the parameters estimates for Ĝ32 for the three compared methods.
b01 = 0.4
b02 = 0.5
a01 = 0.5
a02 = 0.15
Method

Es b̂1

N · Vars b̂1

Es b̂2

N · Vars b̂2

Es â1

N · Vars â1

Es â2

N · Vars â2

SMPE

0.34

1.9

0.44

2.1

0.60

5.0

0.23

3.0

NEB

0.34

0.30

0.44

0.30

0.65

1.0

0.26

0.84

NEBX

0.36

0.11

0.45

0.16

0.63

0.68

0.25

0.55

FITNEB/NEBX

1

1
0.8

0.95

0.6
0.4

NEB
NEBX

0.9
0.9

0.95

0.2
1

FITSMPE

SMPE

NEB

NEBX

Fig. 10. Box plot of the fit of the impulse response of G32 obtained by the methods SMPE, NEB and NEBX respectively.

Fig. 9. Fit of impulse response coefficients of G31 for SMPE
compared with NEB and NEBX respectively for 100 Monte
Carlo simulations. The black line represents y = x, i.e., when
SMPE performs equally good as NEB and NEBX.

7

Conclusion

In this paper, we have addressed the identification of a
module in dynamic networks with known topology. The
problem is cast as the identification of a set of systems
in series connection. The second system corresponds to
the target module, while the first represents the dynamic
relation between exogenous signals and the input and
the target module. This system is modeled following a
Bayesian kernel-based approach, which enables the identification of the target module using empirical Bayes ar-

guments. In particular, the target module is estimated
using a marginal likelihood criterion, whose solution is
obtained by a novel iterative scheme designed through
the ECM algorithm. The method is extended to incorporate measurements downstream of the target module,
which numerical experiments suggest increases performance.

13

FITNEB/NEBX

1

Proof of Proposition 5.1
Using Bayes’ rule we can decompose the complete likelihood as
log p(z, s11 , f ; η) = log p(z|s11 , f ; η)
+ log p(s11 ; η) + log p(f ; η) ,

0.8

and we will analyze each term in turn. First, note that
NEB
NEBX
0.4

0.6
FITSMPE

0.8

−2 log p(s11 |η) = log det λs Kβs + sT11 (λs Kβs )−1 s11

= log det λs Kβs + Tr (λs Kβs )−1 s11 sT11

1

Replacing s11 sT11 with its sample estimate yields the first
term in (79). Similarly,

Fig. 11. Fit of impulse response coefficients of G32 for SMPE
compared with NEB and NEBX respectively for 100 Monte
Carlo simulations. The black line represents y = x, i.e., when
SMPE performs equally good as NEB and NEBX. Note the
scaling of the x-axis of this figure.

A


−2 log p(f |η) = log det λf Kβf + Tr (λf Kβf )−1 f f T .

Replacing f f T with its sample estimate yields the second term in (79). Finally,

Appendix
−2 log p(z|t, s11 ; η) = log det Σ + (z − ẑ)T Σ −1 (z − ẑ) ,
(A.4)

Proof of Lemma 4.1
From Bayes’ rule it follows that
log p(z, s11 ; η̂

(k)

) = log p(z|s11 , η̂

with
(k)

)+log p(s11 ; η̂

(k)

),



Rs






ẑ := 
 Gθ Rs .
Gθ Rv

with (neglecting constant terms)
1
1
log p(z|s11 , η) ∝ − log det Σe − kz − Wθ s11 k2Σ −1
e
2
2
1 T
1
−1
log p(s11 ; η) ∝ − log det λKβ − s11 (λKβ ) s11 .
2
2

The first term of (A.4) is N times the sum of the logarithms of the noise variances squared. The second term
of (A.4) decomposes into a sum of the (weighted) error
of each signal. Then, the first weighted error is given by

Now we have to take the expectation w.r.t. the posterior
p(s11 |w̃2 ; η̂ (k) ). Developing the second term in the first
equation above and recalling that
n
o
(k)
Ep(s11 |w̃2 ; η̂(k) ) [sT11 As11 ] = Tr AŜ11 ,


2
2
σ12 kw̃1 − Rsk2 = kw̃1 k2 − 2w̃1T Rs + Tr RssT RT

.

Replacing s and ssT with their respective estimates gives
the third term in (79), with the corresponding noise variance term of (A.4) added. Similar calculations on the remaining two weighted errors in (A.4) gives the last two
terms in (79). This concludes the proof.

the statement of the lemma readily follows.
Proof of Proposition 4.2
(k)

In (26), fix Σe to the value Σ̂e (computed inserting
2(k)
2(k)
σ1 and σ2 ). We obtain the θ-dependent terms (A.1)
and (A.2) (after multiplying by a factor −2), where k1
and k2 contain terms independent of θ. Recalling the
definitions of Â(k) and b̂(k) , (35) readily follows.

References
Anderson, B. and Moore, J. (1979). Optimal Filtering,
Prentice-Hall, Englewood Cliffs, N.J., USA.
Bezanson, J., Edelman, A., Karpinski, S. and Shah, V. B.
(2017). Julia: A fresh approach to numerical computing, SIAM Review 59(1): 65–98.
Bottegal, G., Aravkin, A. Y., Hjalmarsson, H. and Pillonetto, G. (2016). Robust EM kernel-based methods
for linear system identification, Automatica 67: 114–
126.

Now, let θ be fixed at the value θ̂(k+1) . The function (26)
can be rewritten as (A.3) (after multiplying by a factor
−2). The results (36) follow by minimizing (A.3) with
respect to σ12 and σ22 . Differentiating w.r.t. σ12 and σ22
and calculating the zeros.

14

(−1)

2
2
(k)
(k)
(k)
Wθ ŝ11 = − 2(k) y T Gθ R1 ŝ11 + k1
−2z T Σ̂e(k)
= − 2(k) y T TN (R1 ŝ11 )gθ + k1
(A.1)
σ2
σ2



−1
o
n
1
1
(k)
(k)
(k)
Tr WθT Σe(k)
Wθ Ŝ11 = 2(k) Tr Gθ R1 Ŝ11 R1T GTθ + k2 = 2(k) vec(Gθ )T (R1 Ŝ11 R1T ⊗ IN )vec(Gθ ) + k2
σ2
σ2
1
(k)
(A.2)
= 2(k) gθT DT (R1 Ŝ11 R1T ⊗ IN )Dgθ + k2 ,
σ2

kw̃1 k22
kw̃2 k22
2w̃1T
2w̃2T
(k)
(k)
(k)
Q0 (σ12 , σ22 , θ̂(k+1) ) = N (log σ12 + log σ22 ) +
+
−
R
ŝ
−
G (k+1) R1 ŝ11
1
11
σ12
σ22
σ12
σ22 θ̂
o
o
n
n
1
1
(k)
(k)
+ 2 Tr R1T R1 Ŝ11 + 2 Tr R1T GTθ̂(k+1) Gθ̂(k+1) R1 Ŝ11
σ1
σ2

(A.3)

steiglitz-mcbride method, submitted to Automatica
(arXiv:1610.08534) .
Forssell, U. and Ljung, L. (1999). Closed-loop identification revisited, Automatica 35: 1215–1241.
Geman, S. and Geman, D. (1984). Stochastic relaxation,
Gibbs distributions, and the Bayesian restoration of
images, IEEE Transactions on pattern analysis and
machine intelligence (6): 721–741.
Gilks, W., Richardson, S. and Spiegelhalter, D. (1995).
Markov Chain Monte Carlo in Practice, Chapman &
Hall/CRC Interdisciplinary Statistics, Taylor & Francis.
Gilson, M. and Van den Hof, P. M. J. (2005). Instrumental variable methods for closed-loop system identification, Automatica 41(2): 241–249.
Gunes, B., Dankers, A. and Van den Hof, P. M. J.
(2014). A variance reduction technique for identification in dynamic networks, Proceedings of the 19th
IFAC World Congress.
Haber, A. and Verhaegen, M. (2014). Subspace identification of large-scale interconnected systems, IEEE
Transactions on Automatic Control 59(10): 2754–
2759.
Hayden, D., Yuan, Y. and Gonçalves, J. (2014). Network
reconstruction from intrinsic noise: Minimum-phase
systems, Proceedings of the 2014 American Control
Conference, pp. 4391–4396.
Hjalmarsson, H. (2009). System identification of complex and structured systems, European J. of Control
15(3-4): 275–310.
Ljung, L. (1998). System identification, Springer.
Materassi, D. and Innocenti, G. (2010). Topological
identification in networks of dynamical systems, IEEE
Transactions on Automatic Control 55(8): 1860–1871.
Materassi, D. and Salapaka, M. V. (2012). On the problem of reconstructing an unknown topology via locality properties of the Wiener filter, IEEE Transactions
on Automatic Control 57(7): 1765–1777.
Meng, X. L. and Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A general

Bottegal, G., Risuleo, R. S. and Hjalmarsson, H. (2015).
Blind system identification using kernel-based methods, IFAC-PapersOnLine 48(28): 466–471.
Casella, G. (2001). Empirical Bayes Gibbs sampling,
Biostatistics 2(4): 485–500.
Chen, T., Ohlsson, H. and Ljung, L. (2012). On the estimation of transfer functions, regularizations and gaussian processes - revisited, Automatica 48(8): 1525–
1535.
Chiuso, A. and Pillonetto, G. (2012). A Bayesian approach to sparse dynamic network identification, Automatica 48(8): 1553–1565.
Dankers, A. and Van den Hof, P. M. J. (2015). Nonparametric identification in dynamic networks, Proceedings of the 54th IEEE Conference on Decision and
Control, pp. 3487–3492.
Dankers, A., Van den Hof, P. M. J., Bombois, X. and
Heuberger, P. S. (2015). Errors-in-variables identification in dynamic networks - Consistency results for an
instrumental variable approach, Automatica 62: 39–
50.
Dankers, A., Van den Hof, P. M. J. and Heuberger,
P. S. C. (2013). Predictor input selection for direct
identification in dynamic networks, Proceedings of the
52nd IEEE Annual Conference on Decision and Control, IEEE, pp. 4541–4546.
Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977).
Maximum likelihood from incomplete data via the EM
algorithm, J. of the royal statistical society. Series B
(methodological) pp. 1–38.
Everitt, N., Bottegal, G., Rojas, C. R. and Hjalmarsson, H. (2016). Identification of modules in dynamic
networks: An empirical bayes approach, Proceedings
of the 55th IEEE Annual Conference on Decision and
Control, IEEE, pp. 4612–4617.
Everitt, N., Bottegal, G., Rojas, C. R. and Hjalmarsson,
H. (2017). Variance analysis of linear simo models
with spatially correlated noise, Automatica 77: 68–81.
Everitt, N., Galrinho, M. and Hjalmarsson, H.
(2016). Optimal model order reduction with the

15

framework, Biometrika 80(2): 267–278.
Meyn, S. and Tweedie, R. L. (2009). Markov chains and
stochastic stability; 2nd ed., Cambridge Mathematical
Library, Cambridge Univ. Press, Leiden.
Pillonetto, G. and De Nicolao, G. (2010). A new kernelbased approach for linear system identification, Automatica 46(1): 81–93.
Pillonetto, G., Dinuzzo, F., Chen, T., De Nicolao, G. and
Ljung, L. (2014). Kernel methods in system identification, machine learning and function estimation: A
survey, Automatica 50(3): 657–682.
Prando, G., Chiuso, A. and Pillonetto, G. (2014).
Bayesian and regularization approaches to multivariable linear system identification: the role of rank
penalties, Proceedings of the 53rd IEEE Annual Conference on Decision and Control, pp. 1482–1487.
Rasmussen, C. and Williams, C. (2006). Gaussian Processes for Machine Learning, The MIT Press.
Risuleo, R. S., Bottegal, G. and Hjalmarsson, H. (2015).
A kernel-based approach to Hammerstein system
identification, IFAC-PapersOnLine 48(28): 1011–
1016.
Torres, P., van Wingerden, J. W. and Verhaegen, M.
(2014). Output-error identification of large scale
1D-spatially varying interconnected systems, IEEE
Transactions on Automatic Control 60(1): 130–142.
Van den Hof, P. M. J., Dankers, A., Heuberger, P. S. C.
and Bombois, X. (2013). Identification of dynamic
models in complex networks with prediction error
methods - basic methods for consistent module estimates, Automatica 49(10): 2994–3006.
Wahlberg, B. (1991). System identification using Laguerre models, IEEE Transactions on Automatic Control 36: 551–562.
Wahlberg, B., Hjalmarsson, H. and Mårtensson, J.
(2009). Variance results for identification of cascade
systems, Automatica 45(6): 1443–1448.

16

