Online Identification of Time-Varying Systems: a Bayesian approach

arXiv:1609.07393v1 [cs.SY] 23 Sep 2016

G. Prando, D. Romeres and A. Chiuso

Abstract— We extend the recently introduced regularization/Bayesian System Identification procedures to the estimation
of time-varying systems. Specifically, we consider an online
setting, in which new data become available at given time steps.
The real-time estimation requirements imposed by this setting
are met by estimating the hyper-parameters through just one
gradient step in the marginal likelihood maximization and by
exploiting the closed-form availability of the impulse response
estimate (when Gaussian prior and Gaussian measurement
noise are postulated). By relying on the use of a forgetting factor,
we propose two methods to tackle the tracking of time-varying
systems. In one of them, the forgetting factor is estimated
by treating it as a hyper-parameter of the Bayesian inference
procedure.

I. I NTRODUCTION
The identification of time-varying systems plays a key role in
different applications, such as adaptive and model predictive
control, where a good real-time tracking of the system to
be controlled is necessary. In addition, the detection of
changes or drifts in plant parameters is crucial in terms
of process monitoring and fault detection. Online System
Identification (SysId) and the estimation of time-varying
systems are typically strictly connected problems: one would
like to exploit the new data that become available in order
to track possible changes in the system dynamics.
Recursive Prediction Error Method (RPEM), a variant of
the classical PEM [1], [2], represents nowadays a wellestablished technique, through which the current estimate can
be efficiently updated, as soon as new data are provided.
RPEM are parametric approaches, relying on Recursive
Least-Squares (RLS) routines, which compute the parameter
estimate by minimizing a functional of the prediction errors
[1]. An extension of these approaches to the identification of
time-varying systems involves the adoption of a forgetting
factor, through which old data become less relevant in the
estimation criterion. Convergence and stability properties of
Forgetting Factor RPEM have been well-studied within the
SysId community [3], [4].
Alternative approaches model the coefficients trajectories as
stochastic processes [5], thus exploiting Kalman filtering [6]
or Bayesian inference [7] for parameter estimation. Combinations of bases sequences (e.g. wavelet basis [8]) have also
been considered to model the parameters time evolution.
The above-mentioned parametric procedures share the criticality of the model selection complexity: this step is especially crucial when the model complexity has to be modified
This work has been partially supported by the FIRB project “Learning
meets time” (RBFR12M3AC) funded by MIUR.
† Dept. of Information Engineering, University of Padova (e-mail:
{romeresd,prandogi,chiuso}@dei.unipd.it)

†

in response to changes in the true system dynamics. In
addition, classical complexity selection rules (e.g. crossvalidation or information criteria) may not be applicable in
online settings, due to the excessive computational effort
they require. Model complexity issues have been partially
addressed in the SysId community through the recent introduction of non-parametric methods, relying on Gaussian
processes and Bayesian inference [9], [10]. In this framework
model complexity is tuned in a continuous manner by
estimating the hyper-parameters which describe the prior
distribution chosen by the user [11]. This property makes
these new techniques appealing for the online identification
of time-varying systems: indeed, model complexity can be
continuously adapted whenever new data become available.
In a previous work [12] we started exploring this research
direction by adapting the newly introduced Bayesian procedures to an online identification setting. The methodologies
proposed in [12] are extended in this new paper by dealing
with time-varying systems. Two approaches, relying on the
use of a forgetting factor, are proposed; in particular, following the approach in [13], we investigate the online estimation
of the forgetting factor by treating it as a hyper-parameter
of the Bayesian inference procedure. These techniques are
experimentally compared with the classical parametric counterparts: the results appear favourable and promising for the
methods we propose.
The paper is organized as follows. Sec. II presents the online
identification framework and the challenges we will try to
address. Sec. III provides a brief review of parametric realtime identification techniques, while Sec. IV illustrates the
Bayesian approach to linear SysId, both in the batch and
online scenarios. In particular, Sec. IV-C focuses on the
estimation of time-varying systems. Experimental results are
reported in Sec. V, while conclusions are drawn in Sec. VI.
II. P ROBLEM F ORMULATION
Consider a dynamical system described through an outputerror model, i.e.:
y(t) = [h ∗ u] (t) + e(t),

y(t), u(t) ∈ R

(1)

where h(t) denotes the model impulse response and e(t) is
assumed to be a zero-mean Gaussian noise with variance σ 2 .
SysId techniques aim at estimating the impulse response h of
N
the system, once a set D := {y(t), u(t)}t=1
of measurements
of its input and output signals is provided.
In this work we consider an online setting, in which a new
set of input-output measurements becomes available every T
time steps. Specifically, let us define the variable i := k/T by
assuming w.l.o.g. that k is a multiple of T , and the ith −dataset

iT
as Di = {u(t), y(t)}t=(i−1)T
+1 .
We suppose that at time k an impulse response estimate ĥ(i)
has been computed using the data coming from a collection
S
iT
of previous datasets il=1 Dl = {u(t), y(t)}t=1
; at time k + T
new data Di+1 become available and we would like to update
the previous estimate ĥ(i) by exploiting them. In addition
we assume that the underlying system undergoes certain
variations that we would like to track: this situation could
often arise in practice, due to e.g. variations of the internal
temperature, of the masses (e.g. after grasping an object).
Furthermore, online applications typically require that the
new estimate is available before the new dataset Di+2 is
provided, thus limiting the computational complexity and the
memory storage of the adopted estimation methods.
In this paper, the recently proposed Bayesian approach to
SysId [10] is adapted in order to cope with the outlined
online setting. Its performances are compared with the ones
achieved using classical parametric approaches.
Remark 1: We stress that in the remainder of the paper
we will use the indexes k and iT interchangeably.

III. PARAMETRIC A PPROACH
Standard parametric approaches to SysId rely on the a-priori
choice of a model class M (e.g. ARX, ARMAX, OE, etc.),
which is completely characterized by a parameter θ ∈ Rm .
A. Batch Approach
N
In the batch setting, when a dataset D = {y(t), u(t)}t=1
is
provided, the identification procedure reduces to estimate θ
by minimizing the sum of squared prediction errors:
1 N
θ̂ = arg minm VN (θ , D) = arg minm ∑ (y(t) − ŷ(t|θ ))2 (2)
θ ∈R
θ ∈R 2 t=1
where ŷ(t|θ ) denotes the one-step ahead predictor [1].
B. Online Approach
The extension of these procedures to an online setting relies
on RLS (or pseudo LS) methods.
For ease of notation, let us assume T = 1 in this section.
Suppose that at time k + 1 a new input-output data pair Di+1
is provided; then θ̂ (i) is updated as:
−1

θ̂ (i+1) = θ̂ (i) + µ (i+1) Q(i+1) ∇θ Vk+1 (θ̂ (i) ,

l=1 Dl )

Si+1

(3)

where
Dl ) denotes the gradient of the
loss function computed in the previous estimate and in the
new data; µ (i+1) ∈ R and Q(i+1) ∈ Rm×m are appropriate
scalings which assume different shapes according to the
specific algorithm which is adopted (see [2] and [1], Ch.
11, for further details). Notice that (3) is simply a scaled
S
gradient step w.r.t. the loss function Vk+1 (θ , i+1
l=1 Dl ).
S
∇θ Vk+1 (θ̂ (i) , i+1
l=1

C. Dealing with time-varying systems
In order to cope with time-varying systems, a possible
strategy involves the inclusion of a forgetting factor γ̄ in
the loss function Vk (θ , D):
Vk (θ , D) =
γ

1 k k−t
∑ γ̄ (y(t) − ŷ(t|θ ))2 ,
2 t=1

γ̄ ∈ (0, 1]

(4)

In this way old measurements become less relevant for
the computation of the estimate. A recursive update of the
estimate θ̂ (i) (as the one in (3)) can be derived ([1], Ch. 11).
As an alternative, a sliding window approach can be adopted:
at each time step only the last Nw data are used for computing
the current estimate (with Nw being the window length).
However, since this approach does not admit an update rule
as the one in (3), the computational complexity of the new
estimate will depend on the window length.
A crucial role in the application of parametric SysId techniques is played by the model order selection step: once a
model class M is fixed, its complexity has to be chosen
using the available data. This is typically accomplished
by estimating models with different complexities and by
applying tools such as cross-validation or information criteria
to select the most appropriate one. However, the estimation of
multiple models may be computationally expensive, making
this procedure not suited for the online identification of timevarying systems. Indeed, in this framework, it should ideally
be applied every time new data become available.
The recently proposed approach to SysId, relying on regularization/Bayesian techniques, overpasses the above-described
issue by jointly performing estimation and order selection. Next section will illustrate how the batch regularization/Bayesian method can be tailored to the online identification of time-varying systems.
IV. R EGULARIZATION /BAYESIAN A PPROACH
A. Batch Approach
We discuss how the regularization/Bayesian technique
works in the standard batch setting, i.e. when data D =
N
{y(t), u(t)}t=1
are given. For future use, let us define the
vector YN = [y(1) ... y(N)]> ∈ RN .
According to the Bayesian estimation, the impulse response
h is considered as a realization of a stochastic process with
a prior distribution pη (h), depending on some parameters
η ∈ Ω. The prior pη (h) is designed in order to account for
some desired properties of the estimated impulse response,
such as smoothness and stability [9], [10]. In the Bayesian
framework, the parameters η are known as hyper-parameters
and they need to be estimated from the data, e.g. by optimizing the so-called marginal likelihood (i.e. the likelihood
once the latent variable h has been integrated out) [11]:
Z

η̂ = arg max pη (YN ) = arg max
η∈Ω

η∈Ω

p(YN |h)pη (h)dh

(5)

Once the hyper-parameters η have been estimated, the
minimum variance estimate of h needs to be computed; it
coincides with the posterior mean given the observed data:
Z

ĥ := Eη̂ [h|YN ] =

h

p(YN |h)pη̂ (h)
dh
pη̂ (YN )

(6)

In the SysId context, h is typically modelled as a zeromean Gaussian process (independent of the noise e(t)) with
covariance E [h(t), h(s)] = K̄η (t, s) (aka kernel in the Machine
Learning literature) [9], [14]. Thanks to this asssumption, the
marginal likelihood pη (YN ) is Gaussian and the estimate (6)

is available in closed form.
Furthermore, for simplicity the IIR model in (1) can be accurately approximated by a FIR model of order n, whenever n is
chosen large enough to catch the relevant components of the
system dynamics. By collecting in h := [h(1) · · · h(n)]> ∈
Rn the first n impulse response coefficients, the following
Gaussian prior can be defined:
pη (h) ∼ N (0, Kη ),

η ∈ Ω ⊂ Rd , Kη ∈ Rn×n

(7)

as 1st or 2nd order optimization algorithms [15] or the
Expectation-Maximization (EM) algorithm [16], [17]. Since
these methods may require a large number of iterations
before reaching convergence, they may be unsuited for
online applications. We should recall that, when adopted
for marginal likelihood optimization, each iteration of these
algorithms has a computational complexity of O(n3 ), due to
the objective function evaluation. Specifically, f(i+1)T (η) can
be robustly evaluated as [18]

The hyper-parameters η can then be estimated by solving
η̂ = arg min − ln pη (YN ) = arg min fN (η)
η∈Ω

u(−1)
..
.

u(N) u(N − 1)

(8)

···
..
.
···


u(−n + 1)

..

.
u(N − n + 1)

(9)
(10)

(11)

In the batch setting we are considering the quantities
u(−n + 1), ..., u(0) can be either estimated or set to zero.
Here, we follow the latter option. Once η̂ has been computed,
the corresponding minimum variance estimate is given by
b : = Eη̂ [h|YN ] = arg min (YN − ΦN h)> (YN − ΦN h) + σ 2 h> K −1 h
h
η̂
∈ n
h R

2 −1 −1 >
= (Φ>
N ΦN + σ Kη̂ ) ΦN YN

+ σ −2 ( Y

η∈Ω

fN (η) = YN> Σ(η)−1YN + ln det Σ(η)
2
Σ(η) = ΦN Kη Φ>
N + σ IN
where ΦN ∈ RN×n :

u(0)

ΦN :=  ...

f(i+1)T (η) =((i + 1)T − n) ln σ 2 + 2 ln |S|
(i+1)

− kS−1 L>Ye (i+1) k22 )

(16)

where L and S are Cholesky factors: Kη =: LL> and σ 2 In +
L> R(i+1) L =: SS> (whose computation is O(n3 )).
To tackle the real-time constraints, the approach proposed
in [12] is adopted: η̂ (i+1) is computed by running just one
iteration of a Scaled Gradient Projection (SGP) algorithm
(a 1st order optimization method) applied to solve problem
(8) [15]. Algorithm 1 summarizes its implementation. Notice
that it is initialized
with the previous estimate η̂ (i) (obtained
Si
using the data l=1 Dl ) which is likely to be close to a local
optimum of the objective function fiT (η) ≡ fk (η). If the
number of new data T << k, it is reasonable to suppose
that arg minη∈Ω fiT (η) ≈ arg minη∈Ω f(i+1)T (η). Therefore,
by just performing one SGP iteration, η̂ (i+1) will be sufficiently close to a local optimum of f(i+1)T (η).

(12)

b in (12) can be computed once
Remark 2: The estimate h
a noise variance estimate σ̂ 2 is available. For this purpose,
σ 2 can be treated as a hyper-parameter and estimated by
solving (8) or it can be computed from a LS estimate of h.
In this work the latter option is adopted.
B. Online Approach
We now adapt the batch technique described in Sec. IVA to the online setting outlined in Sec. II. At time k + T ,
(i+1)T
when data Di+1 = {u(t), y(t)}t=iT +1 are provided, the current
b(i) is updated through formula
impulse response estimate h
(12), once the data matrices are enlarged with the new data
and a new hyper-parameter estimate η̂ (i+1) is computed. The
data matrices are updated through the following recursions


(i+1)T > (i+1)T
(i)
ΦiT +1 (13)
R(i+1) := Φ>
(i+1)T Φ(i+1)T = R + ΦiT +1

>
(i+1)T
(i+1)T
e (i)
Ye (i+1) := Φ>
YiT +1
(14)
(i+1)T Y(i+1)T = Y + ΦiT +1


(i+1)
(i)
(i+1)T > (i+1)T
>
Y
:= Y(i+1)T
Y(i+1)T =Y + YiT +1
YiT +1
(15)
(i+1)T

where Y(i+1)T = [y(1) · · · y(iT + T )]> ∈ R(i+1)T , YiT +1 =
[y(iT + 1) · · · y(iT + T )]; Φi is defined as in (11) with N
(i+1)T
replaced by (i + 1)T , while ΦiT +1 has the same structure
of matrix (11) but it contains the data from iT − n + 1 to
(i + 1)T . The computational cost of (13)-(15) is, O(n2 T ),
O(nT ) and O(T 2 ), respectively.
The minimization of f(i+1)T (η) in (9), needed to determine
η̂ (i+1) , is typically performed through iterative routines, such

Algorithm 1 1-step Scaled Gradient Projection (SGP)

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:

Inputs: previous estimates {η̂ (i) , η̂ (i−1) }, ∇ fiT (η̂ (i−1) ),
(i+1)
2
R(i+1) , Ye (i+1) , Y
, σ̂ (i+1)
Initialize: c = 10−4 , δ = 0.4
Compute ∇ f(i+1)T (η̂ (i) )
r(i−1) ← η̂ (i) − η̂ (i−1)
w(i−1) ← ∇ f(i+1)T (η̂ (i) ) − ∇ fiT (η̂ (i−1) )
Approximate the inverse Hessian of f(i+1)T (η̂ (i) ) as
B(i) = α (i) D(i) (using the procedure outlined in [15])
Project onto the feasible set:
z ← ΠΩ,D(i) ( η̂ (i) − B(i) ∇ f(i+1)T (η̂ (i) ) )
∆η̂ (i) ← z − η̂ (i)
ν ←1
if
f(i+1)T (η̂ (i) + ν∆η̂ (i) ) ≤ f(i+1)T (η̂ (i) ) +
cν∇ f(i+1)T (η̂ (i) )> ∆η̂ (i) then
Go to step 12
else
ν ← δν
η̂ (i+1) ← η̂ (i) + ν∆η̂ (i)
Output: η̂ (i+1)

The key step in Algorithm 1 is 4, where the inverse Hessian
is approximated as the product between the positive scalar
α (i) ∈ R+ and the diagonal matrix D(i) ∈ Rd×d . α (i) is chosen
by alternating the so-called Barzilai-Borwein (BB) rules [19]:
>

(i)

α1 :=

r(i−1) r(i−1)
,
>
r(i−1) w(i−1)

>

(i)

α2 :=

r(i−1) w(i−1)
>
w(i−1) w(i−1)

(17)

with r(i−1) and w(i−1) specified at steps 2 and 3 of Algorithm
1. The definition of D(i) depends on the constraints set and
on the objective function. The authors in [15] exploit the
following decomposition of ∇η f(i+1)T (η) (defined in (9)):
∇η f(i+1)T (η) = V (η) −U(η),

V (η) > 0, U(η) ≥ 0 (18)

to specify D(i) . We refer the interested reader to [15] for
further details.
The projection operator adopted at step 5 of Algorithm 1 is
ΠΩ,D(i) (z) = arg minx∈Ω

−1
(x − z)> D(i) (x − z)

In this section we deal with the identification of timevarying systems: specifically, estimators have to be equipped
with tools through which past data become less relevant for
the current estimation. In the following we propose two
routines which combine the “online Bayesian estimation”
above sketched with the ability to “forget” past data.
1) Fixed Forgetting Factor: Following a classical practice
in parametric SysId (see Sec. III), we introduce a forgetting
factor γ̄ ∈ (0, 1] into the data we are provided in order to base
the estimation mainly on the more recent data. Specifically,
we assume that the first k data are generated according to
the following linear model:
ḠkYk = Ḡk Φk h + E, E = [e(1)...e(k)]> ∼ N (0, σ 2 Ik ) (20)

where Ḡk Ḡk =: Γ̄k and Γ̄k := diag γ̄ k−1 , γ̄ k−2 , ..., γ̄ 0 . Therefore, when adopting the regularized regression criterion (12),
the estimate at time k is computed as:
(21)

= arg minn (Yk − Φk h)> Γ̄k (Yk − Φk h) + σ 2 h> Kη̂−1 h

h∈R
>
= (Φk Γ̄k Φk + σ 2 Kη̂−1 )−1 Φ>
k Γ̄kYk

(22)

Correspondingly, the hyper-parameters are estimated solving:
η̂ = arg min Yk> Ḡk Σγ̄ (η)−1 ḠkYk + ln det Σγ̄ (η)

(23)

η∈Ω

2
Σγ̄ (η) = Ḡk Φk Kη Φ>
k Ḡk + σ Ik

Algorithm 2
identification
In particular,
and η̂ (i) are

Inputs: forgetting factor γ̄, previous estimates
(i) (i) (i)
{η̂ (i) , η̂ (i−1) }, previous data matrices {Rγ̄ , Yeγ̄ ,Yγ̄ },
(i+1)T

C. Dealing with time-varying systems

h∈R t=1

Algorithm 2 Online Bayesian SysId: Fixed Forgetting Factor

(19)

Remark 3: Besides SGP, in [12] other inverse Hessian
approximations are investigated (e.g. the BFGS formula).
In this work we only consider the SGP approximation,
since it appears preferable to the others, according to the
experiments we performed (both in the time-invariant and
-variant domain). [12] also considers the EM algorithm as
an alternative to 1st order optimization methods to solve
problem (8). Even if the results reported for EM in [12]
are comparable to the ones achieved through SGP, the latter
approach appears superior to EM in the time-varying setting
we are considering.

k

bγ̄ := arg min ∑ γ̄ k−t y(t) − Φtt h 2 + σ 2 h> K −1 h
h
η̂
n

solving, respectively, (21) and (23); these estimates are
then online updated after the new data Di+1 are provided. Once γ̄ is chosen by the user, it is inserted in
(i+1)
e
the data matrices Rγ̄
:= Φ>
(i+1)T Γ̄(i+1)T Φ(i+1)T , Yγ̄ (i+1) :=
>
>
Φ(i+1)T Γ̄(i+1)T Y(i+1)T , Yγ̄ (i+1) := Y(i+1)T Γ̄(i+1)T Y(i+1)T , updated at steps 1-3 of the algorithm.

(24)

illustrates the online implementation of the
procedure based on equations (22) and (23).
b(i)
it assumes that at time k the estimates h
available and they have been computed by

new data Di+1 = {u(t), y(t)}t=iT +1


(i+1)
(i)
(i+1)T >
(i+1)T
1: Rγ̄
← γ̄ T Rγ̄ + ΦiT +1
Γ̄T ΦiT +1

>
e (i+1) ← γ T Yeγ(i) + Φ(i+1)T Γ̄T Y (i+1)T
2: Y
γ̄
iT +1
iT +1

>
(i+1)
(i)
(i+1)T
(i+1)T
← γ̄ TYγ̄ + YiT +1
Γ̄T YiT +1
3: Y γ̄
−1

4:

b(i+1) ← R(i+1) Ye (i+1)
h
γ̄
γ̄
LS
2

5: σ̂ (i+1) ←

1
(i+1)T −n



(i+1)
(i+1)> b(i+1) b(i+1)> (i+1) b(i+1)
Ȳγ̄
− 2Yeγ̄
hLS + hLS
Rγ̄
hLS

η̂ (i+1) ← arg minη∈Ω f(i+1)T (η) (use Algorithm 1)

−1
2
(i+1)
b(i+1) ← R(i+1) + σ̂ (i+1) K −1
7: h
Yeγ̄
(i+1)
γ̄
γ̄
η̂
b(i+1) , η̂ (i+1)
Output: h

6:

2) Treating the Forgetting Factor as a Hyper-parameter:
The Bayesian framework provides the user with the possibility to treat the forgetting factor as a hyper-parameter and
to estimate it by solving:
η̂, γ̂ = arg minη∈Ω,γ∈(0,1] fk (η, γ)
fk (η, γ) = Yk> Gk Σ(η, γ)−1 GkYk + ln det Σ(η, γ)
2
Σ(η, γ) = Gk Φk Kη Φ>
k Gk + σ Ik

(25)
(26)

(27)

where Gk Gk =: Γk and Γk := diag
.
Remark 4: Notice that the model (20) is equivalent to
>

Yk = Φk h + Eγ̄ , Eγ̄ = eγ̄ (1), ..., eγ̄ (k) ∼ N (0, σ 2 Γ̄−1
k )
γ k−1 , γ k−2 , ..., γ 0

Therefore, treating the forgetting factor as a hyper-parameter
is equivalent to modeling the noise with a non-constant
variance and to give to the diagonal entries of the covariance
matrix an exponential decaying structure.
The online implementation of this approach is detailed in
Algorithm 3, where

>
(i)
(i−1)
b(i) ΦiT
Rγ̂γ := γ̂ (i) Rγ̂γ
+ ΦiT
Γ
(28)
T
(i−1)T +1
(i−1)T +1
b(i) = diag((γ̂ (i) )T −1 , .., (γ̂ (i) )0 ). Ye (i) and Yγ̂γ(i) are analowith Γ
T
γ̂γ
gously defined.
We should stress that the objective function in (26) does not
admit the decomposition (18); we have
∂ fk (η, γ)
= V (η, γ) +U(η, γ),
∂γ

V (η, γ) > 0, U(η, γ) ≥ 0

Thus, when γ is treated as an hyper-parameter, Algorithm
1 is run setting D(i) = Id at step 4; α (i) is still determined
alternating the BB rules (17).

Algorithm 3 Online Bayesian SysId: Forgetting Factor as a
hyper-parameter
Inputs: previous estimates {η̂ (i) , η̂ (i−1) , γ̂ (i) , γ̂ (i−1) }, pre(i) (i) (i)
vious data matrices {Rγ̂γ , Yeγ̂γ ,Yγ̂γ }, new data Di+1 =
(i+1)T

{u(t), y(t)}t=iT +1



(i)
(i+1)T >
(i+1)T
← γ T Rγ̂γ + ΦiT +1
ΓT ΦiT +1

>
eγ(i+1) ← γ T Ye (i) + Φ(i+1)T ΓT Y (i+1)T
2: Y
iT +1
iT +1
γ̂γ


(i+1)
(i)
(i+1)T >
(i+1)T
T
← γ Yγ̂γ + YiT +1
ΓT YiT +1
3: Y γ
(i+1)
(i)
(i)
−1
b
4: h
← (R ) Ye

1:

(i+1)

Rγ

LS

5: σ̂ 2

(i+1)

γ̂γ

←

1
(i+1)T −n

 (i)γ̂γ

(i)
b(i+1)
b(i+1) + (h
b(i+1) )> R(i) h
Yγ̂γ − 2(Yeγ̂γ )> h
LS
LS
γ̂γ LS

η̂ (i+1) , γ̂ (i+1) ← arg minη∈Ω,γ∈(0,1] f(i+1)T (η, γ)
(use Algorithm
1)

−1
(i+1)
b(i+1) ← R(i+1) + σ̂ 2(i+1) (γ̂ (i+1) ) K −1
7: h
Yeγ̂γ
(i+1)
γ̂γ
η̂
b(i+1) , η̂ (i+1)
Output: h

6:

V. E XPERIMENTAL R ESULTS
In this section we test the online algorithms for parametric
and Bayesian SysId described in Sec. III and IV. Their
performance are compared through a Monte-Carlo study over
200 time-varying systems.
A. Data
200 datasets consisting of 3000 input-output measurement
pairs are generated. Each of them is created as follows:
the first 1000 data are produced by a system contained in
the data-bank D4 (used in [20]), while the remaining 2000
data are generated by perturbing the D4-system with two
additional poles and zeros. These are chosen such that the
order of the D4-system changes, thus creating a switch on
the data generating system at time k = 1001.
The data-bank D4 consists of 30th order random SISO
dicrete-time systems having all the poles inside a circle of
radius 0.95. These systems are simulated with a unit variance
band-limited Gaussian signal with normalized band [0, 0.8].
A zero mean white Gaussian noise, with variance adjusted
so that the Signal to Noise Ration (SNR) is always equal to
1, is then added to the output data.
B. Estimators
The parametric estimators are computed with the roe Matlab routine, using the BIC criterion for the model complexity
selection. In the following this estimator will be denoted as
PEM BIC. Furthermore, as a benchmark we introduce the
parametric oracle estimator, called PEM OR, which selects
the model complexity by choosing the model that gives the
best fit to the impulse response of the true system. The order
selection is performed every time a new dataset becomes
available: multiple models with orders ranging from 1 to 20
are estimated and the order selection is performed according
to the two criteria above-described.
For what regards the methods relying on Bayesian inference

we adopt a zero-mean Gaussian prior with a covariance
matrix (kernel) given by the so-called TC-kernel:
KηTC (k, j) = λ min(β k , β j ),

η = [λ , β ]

(29)

with Ω = {(λ , β ) : λ ≥ 0, 0 ≤ β ≤ 1} [14]. The length n
of the estimated impulse responses is set to 100. In the
following, we will use the acronym TC to denote these
methods. Furthermore, the notation OPT will refer to the
standard Bayesian procedure, in which the SGP algorithm
adopted to optimize the marginal likelihood fk (η) is run until
the relative change in fk (η) is less than 10−9 . From here on,
the online counterpart (illustrated in Sec. IV) will be referred
to as the 1-step ML. We will also use the acronyms TC FF
when a fixed forgetting factor is adopted, TC est FF when
the forgetting factor is estimated as a hyper-parameter.
For each Monte Carlo run, the identification algorithms are
300
initialized using the first batch of data Dinit = {u(t), y(t)}t=1
.
After this initial step, the estimators are updated every T = 10
(i+1)T
time steps, when new data Di+1 = {u(t), y(t)}t=iT are
provided. The forgetting factor in the TC FF and PEM
methods is set to 0.998, while its estimation in TC est FF
method is initialized with 0.995.
C. Performance
The purpose of the experiments is twofold. First, we will
compare the two routines we have proposed in Sec. IV-C
to explicitly deal with time-varying systems. Second, we
will compare the parametric and the Bayesian identification
approaches while dealing with time-varying systems.
As a first comparison, we evaluate the adherence of the
b to the true one h, measured as
estimated impulse response h

b 
b = 100 · 1 − kh − hk2
(30)
F (h)
khk2
b at four time instants.
Figure 1 reports the values of F (h)
It is interesting to note that immediately before the change
in the data generating system (k = 1000) the TC methods
slightly outperform the ideal parametric estimator PEM OR.
After the switch (occurring at k = 1001), among regularization/Bayesian routines TC est FF recovers the fit performance a bit faster than TC FF; even at regime it outperforms
the latter because it can choose forgetting factor values that
retain a larger amount of data.
We also observe how the 1-step ML procedures and the
corresponding OPT routines provide analogous performance
at each time step k, validating the method we propose for
online estimation and confirming the results in [12].
The unrealistic PEM OR represents the reference on the
achievable performance of the PEM estimators; it outperforms TC methods in the transient after the switch, while it
has comparable performance at regime. Instead, the recursive
PEM BIC estimator performs very poorly.
As a second comparison, Table I reports the computational
cumulative time of the proposed algorithms in terms of mean
and standard deviation after the estimators are fed with all
3000
the data D = {u(t), y(t)}t=1
. The 1-step ML methods are
one order of magnitude faster than the corresponding OPT

k = 1000

100

update of the Bayesian estimate, resembling the one which
is available for parametric techniques.

b
F (h)

80

R EFERENCES

60
40
TC OPT FF

TC FF

TC est FF

PEM OR

PEM BIC

PEM OR

PEM BIC

PEM OR

PEM BIC

PEM OR

PEM BIC

k = 1300

100

b
F (h)

80
60
40
TC OPT FF

TC FF

TC est FF

k = 1800

100

b
F (h)

80
60
40
TC OPT FF

TC FF

b
F (h)

TC est FF

k = 3000

100
80
60
40
TC OPT FF

TC FF

TC est FF

b achieved at four time instants k (correFig. 1: Fit F (h)
sponding to the number of data available for the estimation).

mean
std

OPT FF
6.70
1.28

TC
FF
0.44
0.03

est FF
5.45
0.67

OR
18.44
0.69

PEM
BIC
18.44
0.69

TABLE I: Computational cumulative time after data D =
3000
{u(t), y(t)}t=1
are used: mean and std over 200 datasets.
ones. The TC est FF estimator is slower than TC FF: this
should be a consequence of having set D(i) = Id in Algorithm
1. On the other hand the RPEM estimators are three times
slower than the OPT ones, thus appearing not particularly
appealing for online applications.
VI. C ONCLUSION AND F UTURE W ORK
We have adopted recently developed SysId techniques relying on the Gaussian processes and Bayesian inference to
the identification of time-varying systems. Specifically, we
have focused on an online setting by assuming that new data
become available at predefined time instants. To tackle the
real-time constraints we have modified the standard Bayesian
procedure: hyper-parameters are estimated by performing
only one gradient step in the corresponding marginal likelihood optimization problem. In order to cope with the timevarying nature of the systems to be identified, we propose
two approaches, based on the use of a forgetting factor.
One of them treats the forgetting factor as a constant, while
the other estimates it as a hyper-parameter of the Bayesian
inference procedure.
We believe that the preliminary investigation performed in
this work may pave the way for further research in this
topic. A future research direction could consider the recursive

[1] L. Ljung, System Identification - Theory for the User, 2nd ed. Upper
Saddle River, N.J.: Prentice-Hall, 1999.
[2] L. Ljung and T. Söderström, Theory and Practice of Recursive
Identificationn, ser. Signal Processing, Optimization, and Control. The
MIT Press, 1983.
[3] S. Bittanti, P. Bolzern, and M. Campi, “Convergence and exponential
convergence of identification algorithms with directional forgetting
factor,” Automatica, vol. 26, no. 5, pp. 929–932, 1990.
[4] L. Guo, L. Ljung, and P. Priouret, “Performance analysis of the
forgetting factor rls algorithm,” International journal of adaptive
control and signal processing, vol. 7, no. 6, pp. 525–538, 1993.
[5] G. C. Chow, “Random and changing coefficient models,” in Handbook
of Econometrics, 1st ed. Elsevier, 1984, vol. 2, ch. 21, pp. 1213–1245.
[6] L. Guo, “Estimating time-varying parameters by the kalman filter
based algorithm: stability and convergence,” Automatic Control, IEEE
Transactions on, vol. 35, no. 2, pp. 141–147, 1990.
[7] A. H. Sarris, A Bayesian approach to estimation of time-varying
regression coefficients. National Bureau of Economic Research, 1973.
[8] M. K. Tsatsanis and G. B. Giannakis, “Time-varying system identification and model validation using wavelets,” Signal Processing, IEEE
Transactions on, vol. 41, no. 12, pp. 3512–3523, 1993.
[9] G. Pillonetto, A. Chiuso, and G. De Nicolao, “Prediction error
identification of linear systems: a nonparametric Gaussian regression
approach,” Automatica, no. 47, pp. 291–305, 2011.
[10] G. Pillonetto, F. Dinuzzo, T. Chen, G. De Nicolao, and L. Ljung, “Kernel methods in system identification, machine learning and function
estimation: a survey,” Automatica, vol. 50, no. 3, pp. 657–682, 2014.
[11] G. Pillonetto and A. Chiuso, “Tuning complexity in kernel-based
linear system identification: the robustness of the marginal likelihood
estimator,” Automatica, vol. 58, pp. 106–117, 2015.
[12] D. Romeres, G. Prando, G. Pillonetto, and A. Chiuso, “On-line
bayesian system identification,” in Proc. of ECC 2016, 2016.
[13] G. Pillonetto and A. Aravkin, “A new kernel-based approach for identification of time-varying linear systems,” in 2014 IEEE International
Workshop on Machine Learning for Signal Processing (MLSP), Sept
2014, pp. 1–6.
[14] T. Chen, H. Ohlsson, and L. Ljung, “On the estimation of transfer
functions, regularizations and Gaussian processes - revisited,” Automatica, vol. 48, no. 8, pp. 1525–1535, 2012.
[15] S. Bonettini, A. Chiuso, and M. Prato, “A scaled gradient projection
methods for Bayesian learning in dynamical systems,” SIAM Journal
on Scientific Computing, p. in press, 2015.
[16] G. Bottegal, A. Y. Aravkin, H. Hjalmarsson, and G. Pillonetto, “Robust
em kernel-based methods for linear system identification,” Automatica,
vol. 67, pp. 114–126, 2016.
[17] G. Bottegal, R. S. Risuleo, and H. Hjalmarsson, “Blind system identification using kernel-based methods,” IFAC-PapersOnLine, vol. 48,
no. 28, pp. 466 – 471, 2015.
[18] T. Chen and L. Ljung, “Implementation of algorithms for tuning
parameters in regularized least squares problems in system identification,” Automatica, vol. 49, no. 7, pp. 2213–2220, 2013.
[19] J. Barzilai and J. M. Borwein, “Two-point step size gradient methods,”
IMA Journal of Numerical Analysis, vol. 8, no. 1, pp. 141–148, 1988.
[20] T. Chen, M. Andersen, L. Ljung, A. Chiuso, and G. Pillonetto, “System
identification via sparse multiple kernel-based regularization using
sequential convex optimization techniques,” IEEE Transactions on
Automatic Control, 2014.

