arXiv:1412.4056v2 [] 19 May 2016

Blind system identification
using kernel-based methods?
Giulio Bottegal, Riccardo S. Risuleo, and Håkan Hjalmarsson

∗

March 15, 2018

Abstract
We propose a new method for blind system identification (BSI). Resorting to a Gaussian regression framework, we model the impulse response of the unknown linear system as a realization of a Gaussian process. The structure of the covariance matrix (or kernel) of such a process
is given by the stable spline kernel, which has been recently introduced
for system identification purposes and depends on an unknown hyperparameter. We assume that the input can be linearly described by few
parameters. We estimate these parameters, together with the kernel hyperparameter and the noise variance, using an empirical Bayes approach.
The related optimization problem is efficiently solved with a novel iterative scheme based on the Expectation-Maximization (EM) method. In
particular, we show that each iteration consists of a set of simple update
rules. We show, through some numerical experiments, very promising
performance of the proposed method.

1

Introduction

In many engineering problems where data-driven modeling of dynamical systems
is required, the experimenter may not have access to the input data. In these
cases, standard system identification tools such as PEM [1] cannot be applied
and specific methods, namely blind system identification (BSI) methods (or blind
deconvolution, if one is mainly interested in the input), need be employed [2].
BSI finds applications in a wide range of engineering areas, such as image
reconstruction [3], biomedical sciences [4] and in particular communications [5,
6], for which literally hundreds of methods have been developed. It would be
impossible to give a thorough literature review here.
∗ G.

Bottegal, R. S. Risuleo and H. Hjalmarsson are with the ACCESS Linnaeus
Center, School of Electrical Engineering, KTH Royal Institute of Technology, Sweden
(risuleo;bottegal;hjalmars@kth.se). This work was supported by the European Research Council under the advanced grant LEARN, contract 267381 and by the Swedish Research Council under contract 621–2009–4017.

1

Clearly, the unavailability of the input signal makes BSI problems generally
ill-posed. Without further information on the input sequence or the structure
of the system, it is impossible to retrieve a unique description of the system [7].
To circumvent (at least partially) this intrinsic non-uniqueness issue, we shall
assume some prior knowledge on the input. Following the framework of [8]
and [9], we describe the input sequence using a number of parameters considerably smaller than the length of the input sequence; see Section 2 for details
and applications.
The main contribution of this paper is to propose a new BSI method. Our
system modeling approach relies upon the kernel-based methods for linear system identification recently introduced in a series of papers [10, 11, 12, 13].
The main advantage of these methods, compared to standard parametric methods [1], is that the user is not required to select the model structure and order
of the system, an operation that might be difficult if little is known about the
dynamics of the system. Thus, we model the impulse response of the unknown
system as a realization of a Gaussian random process, whose covariance matrix,
or kernel, is given by the so called stable spline kernel [10, 14], which encodes
prior information on BIBO stability and smoothness. Such a kernel depends
on a hyperparameter which regulates the exponential decay of the generated
impulse responses.
In the kernel-based framework, the estimate of the impulse response can be
obtained as its Bayes estimate given the output data. However, when applied
to BSI problems, such an estimator is a function of the kernel hyperparameter,
the parameters characterizing the input and the noise variance. All these parameters need to be estimated from data. In this paper, using empirical Bayes
arguments, we estimate such parameters by maximizing the marginal likelihood
of the output measurements, obtained by integrating out the dependence on
the system. In order to solve the related optimization problem, which is highly
non-convex, involving a large number of variables, we propose a novel iterative
solution scheme based on the Expectation-Maximization method [15]. We show
that each iteration of such a scheme consists of a sequence of simple updates
which can be performed using little computational efforts. Notably, our method
is completely automatic, since the user is not required to tune any kind parameter. This in contrast with the BSI methods recently proposed in [8, 9], where,
although the system is retrieved via a convex optimization problem, the user
is required to select some regularization parameters and the model order. The
method derived in this paper follows the same approach used in [16], where a
novel method for Hammerstein system identification is described.
The paper is organized as follows. In the next section, we introduce the BSI
problem and we state our working assumptions. In Section 3, we give a background on kernel-based methods, while in Section 4 we describe our approach
to BSI. Section 5 presents some numerical experiments, and some conclusions
end the paper.

2

2

Blind system identification

We consider a SISO linear time-invariant discrete-time dynamic system (see
Figure 1)
+∞
X
yt =
gi ut−i + vt ,
(1)
i=0
+∞
{gt }t=0

where
is a strictly causal transfer function (i.e., g0 = 0) representing
the dynamics of the system, driven by the input ut . The measurements of the
output yt are corrupted by the process vt , which is zero-mean white Gaussian
noise with unknown variance σ 2 . For the sake of simplicity, we will also hereby
assume that the system is at rest until t = 0.
vt
ut

g

+

yt

Figure 1: Block scheme of the system identification scenario.
We assume that N samples of the output measurements are collected, and
N
denote them by {yt }t=1 . The input u(t) is not directly measurable and only
some information about it is available. More specifically, we assume we know
N −1
that the input, restricted to the N time instants {ut }t=0 , belongs to a certain
subspace of RN and thus can be written as
u = Hx ,

(2)


T
where u = u0 · · · uN −1 . In the above equation, H ∈ RN ×p is a known
matrix with full column rank and x ∈ Rp , p ≤ N , is an unknown vector characterizing the evolution of u(t). Below we report two examples of inputs generated
in this way.
Piecewise constant inputs with known switching instants
Consider a piecewise constant input signal u(t) with known switching instants
T1 , T2 . . . Tp , with Tp = N . The levels the input takes in between the switching
instants are unknown and collected in the vector x. Then, the input signal can

3

be expressed as
u0
u1
..
.


1
 1

 

  ..

 .

 

uT1 −1  1
 

 uT1  
 

 ..  


u= . 
=
uT2 −1  
 

 .  
 ..  

 
 uT  
 p−1  
 .  
 ..  
uTp






1
..
.
1
..

.






 

 x1
 
  x2 
  .  = Hx .
 . 
 .

 xp


1

.. 
.

(3)

1

with
H = diag {1T1 , 1T2 −T1 , . . . , 1Tp −Tp−1 } , H ∈ RN ×p

(4)

where 1m denotes a column vector of length m with all entries equal to 1:
The vector x needs to be estimated from output data. Applications of BSI
with piecewise constant inputs are found in room occupancy estimation [17] and
nonintrusive appliance load monitoring (NIALM) [18, 19].
Combination of known sinusoids
Assume that u is composed by the sum of p sinusoids with unknown amplitude
and known frequencies ω1 , . . . , ωp . Then in this case we have


sin(ω1 ) · · ·
sin(ωp )


..
..
H=
(5)
,
.
.
sin(N ω1 ) · · ·

sin(N ωp )

with ω1 , . . . , ωp such that H is full column rank. The vector x represents
the amplitude of the sinusoids. Applications of this setting are found in blind
channel estimation [9].

2.1

Problem statement

We state our BSI problem as the problem of obtaining an estimate of the impulse
n
N
response gt for n time instants, namely {gt }t=1 , given {yt }t=1 and H. Recall
that, by choosing n sufficiently large, these samples can be used to approximate
gt with arbitrary accuracy [1]. To achieve our goal we will need to estimate the

T
input u = u0 · · · uN −1 ; hence, we might also see our problem as a blind
deconvolution problem.

4

Remark 1 The identification method we propose in this paper can be derived
also in the continuous-time setting, using the same arguments as in [10]. However, for ease of exposition, here we focus only on the discrete-time case.
In a condition of complete information, this problem can be solved by least
squares [1] or using regularized kernel-based approaches [10], [11], [12].

2.2

Identifiability issues

It is well-known that BSI problems are not completely solvable (see e.g. [2,
5]). This because the system and the input can be determined up to a scaling
factor, in the sense that every pair (αu, α1 g), α ∈ R, can describe the output
dynamics equally well. Hence, we shall consider our BSI problem as the problem
of determining the system and the input up to a scaling factor. Another possible
way out for this issue is to assume that kgk2 or g1 are known [20].

3

Kernel-based system identification

In this section we briefly review the kernel-based approach introduced in [10,
11] and show how to readapt it to BSI problem.
Let us first introduce the following vector notation


 
 
 
u0
y1
g1
v1
 .. 
 .. 
 .. 
 .. 
u :=  .  , y :=  .  , g :=  .  , v :=  . 
uN −1

yN

and the operator Tn (·) that, given
Toeplitz matrix, e.g.

u0
0
 u1
u
0


..
Tn (u) =  ...
.

uN −2 uN −3
uN −1 uN −2

gn

vN

a vector of length N , maps it to an N × n

0

···
···
..
.

···
···

uN −n+1
···

0
0
..
.






 ∈ RN ×n .

0 
uN −n

We shall reserve the symbol U for Tn (u). Then, the input-output relation for
the available samples can be written
y = Ug + v .

(6)

In this paper we adopt a Bayesian approach to the BSI problem. Following
a Gaussian process regression approach [21], we model the impulse response as
follows
g ∼ N (0, λKβ )
(7)
where Kβ is a covariance matrix whose structure depends on a shaping parameter β, and λ ≥ 0 is a scaling factor, which regulates the amplitude of the
5

realizations from (7). Given the identifiability issue described in Section 2.2, λ
can be arbitrarily set to 1. In the context of Gaussian regression, Kβ is usually
called a kernel and its structure is crucial in imposing properties on the realizations drawn from (7). An effective choice of kernel for system identification
purposes is given by the so-called stable spline kernels [10, 11]. In particular, in
this paper we adopt the so-called first-order stable spline kernel (or TC kernel
in [12]), which is defined as
{Kβ }i,j := β max(i,j) ,

(8)

where β is a scalar in the interval [0, 1). Such a parameter regulates the decaying
velocity of the generated impulse responses.
Recall the assumption introduced in Section 2 on the Gaussianity of noise.
Due to this assumption, the joint distribution of the vectors y and g is Gaussian,
provided that the vector x (and hence the input u), the noise variance σ 2 and
the parameter β are given. Let us introduce the vector


θ := xT σ 2 β
∈ Rp+2 ,
(9)
which we shall call hyperparameter vector. Then we can write
  

  
0
Σy Σyg
y
θ ∼N
,
,
p
0
Σgy Kβ
g

(10)

where Σyg = ΣTgy = U Kβ and Σy = U Kβ U T + σ 2 I. It follows that the posterior
distribution of g given y (and θ) is Gaussian, namely
p(g|y, θ) = N (Cy, P ) ,
where

P =

UT U
+ Kβ−1
σ2

−1
,

C=P

(11)
UT
.
σ2

(12)

From (11), the impulse response estimator can be derived as the Bayesian estimator [22]
ĝ = E[g|y, θ] = Cy .
(13)
Clearly, such an estimator is a function of θ, which needs to be determined
from the available data y before performing the estimation of g. Thus, the BSI
algorithm we propose in this paper consists of the following steps.
1. Estimate the hyperparameter vector θ.
2. Obtain ĝ by means of (13).
In the next section, we discuss how to efficiently compute the first step of
the algorithm.

6

4

Estimation of the hyperparameter vector

An effective approach to choose the hyperparameter vector characterizing the
impulse response estimator (13) relies on Empirical Bayes arguments [23]. More
precisely, since y and g are jointly Gaussian, an efficient method to choose θ
is given by maximization of the marginal likelihood [24], which is obtained by
integrating out g from the joint probability density of (y, g). Hence, an estimate
of θ can be computed as follows
θ̂ = arg max log p(y|θ) .
θ

(14)

Solving (14) in that form can be hard, because it is a nonlinear and nonconvex problem involving a large number (p + 2) of decision variables. For
this reason, we propose an iterative solution scheme which resorts to the EM
method. To this end, we define the complete likelihood
L(y, g|θ) := log p(y, g|θ) ,

(15)

which depends also on the missing data g. Then, the EM method provides θ̂ by
iterating the following steps:
(E-step) Given an estimate θ̂k after the k-th iteration of the scheme, compute
Q(θ, θ̂k ) := Ep(g|y, θ̂k ) [L(y, g|θ)] ;

(16)

θ̂k+1 = arg max Q(θ, θ̂k ) .

(17)

(M-step) Compute
θ

The iteration of these steps is guaranteed to converge to a (local or global)
maximum of (14) [25], and the iterations can be stopped if kθ̂k+1 − θ̂k k2 is
below a given threshold.
Assume that, at iteration k + 1 of the EM scheme, the estimate θ̂k of θ is
available. Using the current estimate of the hyperparameter vector, we construct the matrices Ĉ k and P̂ k using (12) and, accordingly, we denote by ĝ k the
estimate of g computed using (13), i.e. ĝ k = Ĉ k y and the linear prediction of y
as ŷ k = U ĝ k . Furthermore, let us define


Âk = −H T RT (P̂ k + ĝ k ĝ kT ) ⊗ IN RH
T

b̂k = H T TN (ĝ k ) y ,

(18)

where R ∈ RN n×N is a matrix such that, for any u ∈ RN :
Ru = vec(Tn (u)) .

(19)

Having introduced this notation, we can state the following theorem, which
provides a set of upgrade rules to obtain the hyperparameter vector estimate
θ̂k+1 .
7

Theorem 1 Let θ̂k be the estimate of the hyperparameter vector after the k-th
iteration of the EM scheme. Then


θ̂k+1 = x̂k+1,T σ̂ 2,k+1 β̂ k+1
(20)
can be obtained performing the following operations:
• The input estimate is updated computing
−1 k

x̂k+1 = −(Ak )

b ;

• The noise variance is updated computing
h
i
1 
ky − ŷ k k22 + Tr Û k+1 P̂ k Û k+1,T
,
σ̂ 2,k+1 =
N

(21)

(22)

where Û k+1 denotes the Toeplitz matrix of the sequence ûk+1 = H x̂k+1 ;
• The kernel shaping parameter is updated solving
β̂ k+1 = arg min Q(β, θ̂k ) ,

(23)

h
i
Q(β, θ̂k ) = log det Kβ + Tr Kβ−1 (P̂ k + ĝ k ĝ kT ) .

(24)

β∈[0, 1)

where

Hence, the maximization problem (14) reduces to a sequence of very simple
optimization problems. In fact, at each iteration of the EM algorithm, the
input can be estimated by computing a simple update rule available in closedform. The same holds for the noise variance, whereas the update of the kernel
hyperparameter β does not admit any closed-form expression. However, it can
be retrieved by solving a very simple scalar optimization problem, which can be
solved efficiently by grid search, since the domain of β is the interval [0, 1).
It remains to establish a way to set up the initial estimate θ̂0 for the EM
method. This can be done by just randomly choosing the entries of such a
vector, keeping the constraints β̂ 0 ∈ [0, 1), σ̂ 2,0 > 0.
Below, we provide our BSI algorithm.

5

Numerical experiments

We test the proposed BSI method by means of Monte Carlo experiments. Specifically, we perform 6 groups of simulations where, for each group, 100 random
systems and input/output trajectories are generated. The random systems are
generated by picking 20 zeros and 20 poles with random magnitude and phase.
The zero magnitude is less than or equal to 0.95 and the pole magnitude is no
larger than 0.92. The inputs are piecewise constant signals and the number of
switching instants is p = 10, 20, 30, 40, 50, 60, depending on the group of experiments. We generate 200 input/output samples per experiment. The output
8

Algorithm: Bayesian kernel-based EM Blind System Identification
N
Input: {yt }t=1 , H
n

N

Output: {ĝ}t=1 , {ût }t=1

1. Initialization: randomly set θ̂0 = x̂0T

σ̂ 2,0

β̂ 0



2. Repeat until convergence:
(a) E-step: update P̂ k , Ĉ k from (12) and ĝ k from (13);
(b) M-step: update the parameters:
• x̂k+1 from (21);
• σ̂ k+1 from (22),
• β̂ k+1 from (23)
n

N

3. Compute {ĝ}t=1 from (13) and {ût }t=1 = H x̂;

is corrupted by random noise whose variance is such that σ 2 = var(U g)/10, i.e.
the noise variance is ten times smaller than the variance of the noiseless output.
The goal of the experiments is to estimate n = 50 samples of the generated
impulse responses. We compare the following three estimators.
• B-KB. This is the proposed Bayesian kernel-based BSI method that estimates the input by marginal likelihood maximization, using an EM-based
scheme. The convergence criterion of the EM method is kθ̂k+1 − θ̂k k2 <
10−3 .
• NB-LS. This is an impulse response estimator based on the least squares
criterion. Here, we assume that the input is known, so the only quantity
to be estimated is the system. Hence, this corresponds to an unbiased
FIR estimator [1].
• NB-KB. This is the Bayesian kernel-based system identification method
introduced in [10] and revisited in [12]. Like the estimator NB-LS, this
method has knowledge of the input and in fact corresponds to the estimator B-KB when x is known and not estimated.
The performance of the estimators are evaluated by means of the output fitting
score
kÛi ĝi − Ui gi k2
,
(25)
F IT = 1 −
kUi gi − Ui gi k2
where, at the i-th Monte Carlo run, Ui and Ûi are the Toeplitz matrices of the
true and estimated inputs (Ûi = Ui if the method needs not estimate x), gi and
ĝi are the true and estimated systems and Ui gi is the mean of Ui gi .

9

Figure 2 shows the results of the six group of experiments. As expected,
the estimator NB-KB, which has access to the true input, gives the best performance for all the values of p and in fact is independent of such a value.
Surprisingly, for p = 10 amd p = 20, the proposed BSI method outperforms
the least squares estimator that knows the true input. An example of one such
Monte Carlo experiments in reported in Figure 3. As one might expect, there
p = 20
1

0.95

0.95

0.9

0.9
Fit

Fit

p = 10
1

0.85

0.85

0.8

0.8

0.75

0.75

0.7

B−KB

NB−LS

0.7

NB−KB

B−KB

1

1

0.95

0.95

0.9

0.9

0.85

0.8

0.75

0.75
B−KB

NB−LS

0.7

NB−KB

B−KB

p = 50
1

0.95

0.95

NB−KB

0.9
Fit

0.9
Fit

NB−LS
p = 60

1

0.85

0.85

0.8

0.8

0.75

0.75

0.7

NB−KB

0.85

0.8

0.7

NB−LS
p = 40

Fit

Fit

p = 30

B−KB

NB−LS

0.7

NB−KB

B−KB

NB−LS

NB−KB

Figure 2: Results of the Monte Carlo simulations for different values of p, namely
the number of different levels in the input signals.
is a performance degradation as p increases, since the blind estimator has to
estimate more parameters. Figure 4 shows the median of the fitting score of
each group of experiments as function of p. It appears that, approximately,
there is a linear trend in the performance degradation.

10

Input Signal

0.12

true
B−KB

0.1
0.08
0.06
0.04
0.02
0

0

20

40

60

80

100

120

140

160

180

200

Impulse Response
0.6

true
B−KB

0.4
0.2
0
−0.2
−0.4
5

10

15

20

25

30

35

40

45

50

Output

15

True
B−KB

10
5
0

−5
0

20

40

60

80

100

120

140

160

180

200

Figure 3: Example of one Monte Carlo realization with p = 20. Top panel: true
vs estimated (normalized) input. Middle panel: true vs estimated (normalized)
impulse response. Bottom panel: true vs predicted output.

0.9
0.88

Fit

0.86
0.84
0.82
0.8
0.78
10

20

30

40

50

60

p

Figure 4: Median of the fitting score for each group of Monte Carlo simulations
as function of p.

11

6

Conclusions

In this paper we have proposed a novel blind system identification algorithm.
Under a Gaussian regression framework, we have modeled the impulse response
of the unknown system as the realization of a Gaussian process. The kernel
chosen to model the system is the stable spline kernel. We have assumed that the
unknown input belongs to a known subspace of the input space. The estimation
of the input, together with the kernel hyperparameter and the noise variance,
has been performed using an empirical Bayes approach. We have solved the
related maximization problem resorting to the EM method, obtaining a set of
update rules for the parameters which is simple and elegant, and permits a fast
computation of the estimates of the system and the input. We have shown,
through some numerical experiments, very promising results.
We plan to extend the current method in two ways. First, a wider class of
models of the system, such as the Box-Jenkins model, will be considered. We
shall also attempt to remove the assumption on the input belonging to a known
subspace by adopting suitable Bayesian models.

A

Proof of Theorem 1

First note that p(y, g|θ) = p(y|g, θ)p(g|θ). Hence we can write the complete
likelihood as
L(y, g|θ) = log p(y|g, θ) + log p(g|θ)
(26)
and so
1
N
2
log σ 2 − 2 ky − U gk
2
2σ
1
1
− log det Kβ − g T Kβ−1 g
2
2

N
1
= − log σ 2 − 2 y T y + g T U T U g − 2y T U g
2
2σ
1
1
− log det Kβ − g T Kβ−1 g .
2
2

L(y, g|θ) = −

We now proceed by taking the expectation of this expression with respect to
the random variable g|y, θ̂k . We obtain the following components


N
N
2
(a) : E − log σ = − log σ 2
2
2


1 T
1 T
(b) : E − 2 y y = − 2 y y
2σ
2σ




1 T T
1 T
k
k kT
(c) : E − 2 g U U g = Tr − 2 U U (P̂ +ĝ ĝ )
2σ
2σ


1 T
1 T k
(d) : E 2 y U g = 2 y U ĝ
σ
σ
12



1
1
(e) : E − log det Kβ = − log det Kβ
2
2


i
1
1 h
(f ) : E − g T Kβ−1 g = − Tr Kβ−1 (P̂ k + ĝ k ĝ kT )
2
2
It follows that Q(θ, θ̂k ) is the summation of the elements obtained above. By
inspecting the structure of Q(θ, θ̂k ), it can be seen that such a function splits
in two independent terms, namely
Q(θ, θ̂k ) = Q1 (x, σ 2 , θ̂k ) + Qβ (β, θ̂k ) ,

(27)

Q1 (x, σ 2 , θ̂k ) = (a) + (b) + (c) + (d)

(28)

where
2

is function of x and σ , while
Qβ (β, θ̂k ) = (e) + (f )

(29)

depends only on β and corresponds to (24). We now address the optimization
of (28). To this end we write
1
Qx (x, θ̂k ) + Qσ2 (σ 2 , θ̂k )
(30)
σ2  


1
1
= 2 Tr − U T U (P̂ k +ĝ k ĝ kT ) + y T U ĝ k
σ
2
N
1
−
log σ 2 − 2 y T y .
2
2σ
This means that the optimization of Q1 can be carried out first with respect
to x, optimizing only the term Qx , which is independent of σ 2 and can be
written in a quadratic form
1
Qx (x, θ̂k ) = xT Âk x + b̂kT x .
(31)
2
To this end, first note that, for all v1 ∈ Rn , v2 ∈ Rm ,
Q1 (x, σ 2 , θ̂k ) =

Tm (v1 )v2 = Tn (v2 )v1 .
Recalling (19), we can write
h
i
T
Tr U T U (P̂ k+ĝ k ĝ kT ) =vec(U ) ((P̂ k+ĝ k ĝ kT ) ⊗ IN )vec(U )


1
= − uT RT (P̂ k + ĝ k ĝ kT ) ⊗ IN Ru
2


1
= − xT H T RT (P̂ k + ĝ k ĝ kT ) ⊗ IN RHx ,
2

(32)

(33)

where the matrix in the middle corresponds to Âk defined in (18). For the linear
term we find
y T U ĝ k = y T TN (ĝ k )u = y T TN (ĝ k )Hx ,
(34)
so that the term b̂kT in (18) is retrieved and the maximizer x̂k+1 is as in (21).
Plugging back x̂k+1 into (28) and maximizing with respect to σ 2 we easily find
σ̂ 2,k+1 corresponding to (22). This concludes the proof.
13

References
[1]

L. Ljung. System Identification, Theory for the User. Prentice Hall, 1999.

[2]

K. Abed-Meraim, W. Qiu, and Y. Hua. “Blind system identification”.
Proc. IEEE 85.8 (1997), pp. 1310–1322. doi: 10.1109/5.622507.

[3]

G. R. Ayers and J. C. Dainty. “Iterative blind deconvolution method and
its applications”. Opt. Lett. 13.7 (1988), p. 547. doi: 10.1364/ol.13.
000547.

[4]

D. B. McCombie, A. T. Reisner, and H. H. Asada. “Laguerre-Model Blind
System Identification: Cardiovascular Dynamics Estimated From Multiple
Peripheral Circulatory Signals”. IEEE Trans. Biomed. Eng. 52.11 (2005),
pp. 1889–1901. doi: 10.1109/tbme.2005.856260.

[5]

F. Gustafsson and B. Wahlberg. “Blind equalization by direct examination
of the input sequences”. IEEE Trans. Commun. 43.7 (1995), pp. 2213–
2222. doi: 10.1109/26.392964.

[6]

E Moulines, P Duhamel, J.-F. Cardoso, and S Mayrargue. “Subspace
methods for the blind identification of multichannel FIR filters”. IEEE
Trans. Signal Process. 43.2 (1995), pp. 516–525. doi: 10.1109/78.348133.

[7]

L Tong, R.-W. Liu, V. C. Soon, and Y.-F. Huang. “Indeterminacy and
identifiability of blind identification”. IEEE Trans. Circuits Syst. 38.5
(1991), pp. 499–509. doi: 10.1109/31.76486.

[8]

H. Ohlsson, L. J. Ratliff, R. Dong, and S. S. Sastry. “Blind Identification
Via Lifting”. Proc. IFAC World Cong. 2014. doi: 10.3182/20140824-6za-1003.02567.

[9]

A. Ahmed, B. Recht, and J. Romberg. “Blind Deconvolution Using Convex
Programming”. IEEE Trans. Inform. Theory 60.3 (2014), pp. 1711–1732.
doi: 10.1109/tit.2013.2294644.

[10]

G. Pillonetto and G. De Nicolao. “A new kernel-based approach for linear
system identification”. Automatica 46.1 (2010), pp. 81–93. doi: doi:10.
1016/j.automatica.2009.10.031.

[11]

G. Pillonetto, A. Chiuso, and G. De Nicolao. “Prediction error identification of linear systems: a nonparametric Gaussian regression approach”.
Automatica 47.2 (2011), pp. 291–305. doi: 10.1016/j.automatica.2010.
11.004.

[12]

T. Chen, H. Ohlsson, and L. Ljung. “On the estimation of transfer functions, regularizations and Gaussian processes—Revisited”. Automatica 48.8
(2012), pp. 1525–1535. doi: 10.1016/j.automatica.2012.05.026.

[13]

G. Pillonetto, F. Dinuzzo, T. Chen, G. D. Nicolao, and L. Ljung. “Kernel
methods in system identification, machine learning and function estimation: A survey”. Automatica 50.3 (2014), pp. 657–682. doi: 10.1016/j.
automatica.2014.01.001.

14

[14]

G. Bottegal and G. Pillonetto. “Regularized spectrum estimation using
stable spline kernels”. Automatica 49.11 (2013), pp. 3199–3209. doi: 10.
1016/j.automatica.2013.08.010.

[15]

A. P. Dempster, N. M. Laird, and D. B. Rubin. “Maximum likelihood
from incomplete data via the EM algorithm”. J. R. Stat. Soc. Ser. B.
Stat. Methodol. (1977), pp. 1–38.

[16]

R. S. Risuleo, G. Bottegal, and H. Hjalmarsson. “A kernel-based approach
to Hammerstein system identication”. Proc. IFAC Symp. System Identification (SYSID). Vol. 48. 28. 2015, pp. 1011–1016. doi: doi:10.1016/j.
ifacol.2015.12.263.

[17]

A. Ebadat, G. Bottegal, D. Varagnolo, B. Wahlberg, and K. H. Johansson.
“Estimation of building occupancy levels through environmental signals
deconvolution”. Proc. ACM Workshop Embedded Systems For EnergyEfficient Buildings (BuildSys). Association for Computing Machinery (ACM),
2013. doi: 10.1145/2528282.2528290.

[18]

G. Hart. “Nonintrusive appliance load monitoring”. Proceedings of the
IEEE 80.12 (1992), pp. 1870–1891. doi: 10.1109/5.192069.

[19]

R. Dong, L. Ratliff, H. Ohlsson, and S. S. Sastry. “A dynamical systems
approach to energy disaggregation”. 52nd IEEE Conference on Decision
and Control. Institute of Electrical & Electronics Engineers (IEEE), 2013.
doi: 10.1109/cdc.2013.6760891.

[20]

E. W. Bai and D. Li. “Convergence of the iterative Hammerstein system identification algorithm”. IEEE Trans. Autom. Control 49.11 (2004),
pp. 1929–1940.

[21]

C. Williams and C. Rasmussen. Gaussian processes for machine learning.
2006.

[22]

B. D. O. Anderson and J. B. Moore. Optimal filtering. Courier Corporation, 2012.

[23]

J. Maritz and T. Lwin. Empirical bayes methods. Chapman and Hall London, 1989. isbn: 9780412277603.

[24]

G. Pillonetto and A. Chiuso. “Tuning complexity in kernel-based linear
system identification: The robustness of the marginal likelihood estimator”. Proc. European Control Conf. (ECC). 2014, pp. 2386–2391. doi:
10.1109/ECC.2014.6862629.

[25]

G. McLachlan and T. Krishnan. The EM algorithm and extensions. Vol. 382.
John Wiley and Sons, 2007. isbn: 9780471201700.

15

