Regularized parametric system identification:
a decision-theoretic formulation
Johan Wågberg ∗1 , Dave Zachariah†1 , and Thomas B. Schön‡1

arXiv:1710.04009v1 [] 11 Oct 2017

1

Department of Information Technology, Uppsala University
October 12, 2017
Abstract

Parametric prediction error methods constitute a classical approach to the identification of linear dynamic systems with
excellent large-sample properties. A more recent regularized approach, inspired by machine learning and Bayesian
methods, has also gained attention. Methods based on this approach estimate the system impulse response with excellent small-sample properties. In several applications, however, it is desirable to obtain a compact representation
of the system in the form of a parametric model. By viewing the identification of such models as a decision, we
develop a decision-theoretic formulation of the parametric system identification problem that bridges the gap between
the classical and regularized approaches above. Using the output-error model class as an illustration, we show that
this decision-theoretic approach leads to a regularized method that is robust to small sample-sizes as well as overparameterization.

1

Introduction

Traditionally, methods that identify a parametric model of a dynamical system are formulated using the frequentist
paradigm. A classical parametric approach is the prediction error methods (PEM), which has excellent large-sample
properties [11, 3]. Computational advances have more recently enabled identification methods formulated in the
Bayesian paradigm [6, 5]. A recent nonparametric approach is the kernel regression method, which estimates system
impulse responses in a regularized manner with excellent small-sample properties [7, 8]. The choice of regularization
corresponds to the modelling of prior information about the system under consideration [9, 4].
The classical and the Bayesian approach differ radically in their underlying philosophies, practical interpretations,
and target models. In this paper we propose a way to bridge the gap by developing a decision-theoretic approach for
parametric system identification, cf. [2, 1, 10]. Specifically, by viewing the identification of a specific parametric
system model as a decision, we can unify classical frequentist and modern Bayesian identification approaches. The
framework enables comparisons and fruitful cross-fertilization of ideas from both approaches. Specifically, it dispels
• the need to specify a prior distribution of the model parameters, and

• the assumption that the system belongs to the considered model class.

Moreover, it achieves a natural regularization of the identification problem, which mitigates overfitting
• in the small-sample regime

• when using overparameterized models.
p
Notation: kgkW = g > Wg is a weighted `2 -norm. The operation of extracting the diagonal elements of a square
matrix as a vector is denoted diag (·).
∗ johan.wagberg@it.uu.se
† dave.zachariah@it.uu.se
‡ thomas.schon@it.uu.se

1

2

Model classes

2.1

Nonparametric model specification

Consider a discrete linear time-invariant and causal dynamic system with a single input and a single output (SISO).
The relationship between the input u(t) and the output y(t) of such a system can be written as
y(t) = G? (q)u(t) + e(t),

(1)

t = 1, . . . , N,

where G? (q) is the system transfer function, q is the shift operator and e(t) is a random zero-mean disturbance assumed
to be white with variance λ [3]. We assume that the input u(t) is known for all t and that u(t) ≡ 0 for t ≤ 0. Let

D = (u(1), y(1)), . . . , (u(N ), y(N ))
be input-output data collected at time instants t = 1, . . . , N .
The transfer function can be written as
G? (q) =

∞
X

g? (k)q −k .

(2)

k=0

The collection g? = {g? (k)}∞
k=0 is the impulse response of the system and a linear system is uniquely defined by its
impulse response, which constitutes a nonparametric model of the system. Let G denote the class of all causal SISO
systems. Thus for (1) we have that g? ∈ G, see Figure 1 for an illustration.
Since the input u(t) is zero for t ≤ 0, for finite time N , (1) can be written as a finite sum
y(t) =

t−1
X

k=0

where

(3)

g? (k)u(t − k) + e(t) ⇔ y = Hg ? + e,



u(1)
y(1)
 u(2)
 y(2) 



y =  .  and H =  .
 ..
 .. 


y(N )

u(1)
..
.

u(N ) u(N − 1)

..

.
···





,

u(1)

and Hi,j = 0 for j > i, contain D. The N first coefficients of the impulse response are


g? (0)
 g? (1) 


g? = 

..


.
g? (N − 1)

>
and
 e> = [e(1) e(2) . . . e(N )] contains the disturbance samples, with mean value E [e] = 0 and covariance
E ee = λI.

2.2

Parametric model specification

Estimating the transfer function or, equivalently, the impulse response g? can give useful information about the properties of the system. However, in many applications, e.g. control, a more compact representation of a system is more
useful. We consider a family of impulse responses parameterized by a vector θ and denote it as Gθ . One possible
family is to let the transfer function be a rational function in q −1 , i.e. the ratio of two polynomials,
Gθ (q) =

Bθ (q)
b0 q −nk + b1 q −nk −1 + · · · + bnb q −nk −nb
=
,
Fθ (q)
1 + f1 q −1 + · · · + fnf q −nf
2

(4)

G
g?

gθ
Gθ

Figure 1: The nonparametric class G contains all linear time-invariant causal SISO systems. The
parametrized class Gθ is a subset of G. Note that g? does not have to belong to Gθ .

g?

gθ

1.5
1
0.5
0
−0.5
−1

0

2

4

6

8

10

12

14

16

18

t
Figure 2: The first 20 coefficients of impulse responses. For sake of illustration we chose θ =

>
0.13 −2.65 2.92 −1.64 0.41 .



and θ = b0 , b1 , . . . , bnb , f1 , f2 , . . . , fnf . The orders of the two polynomials are determined by nb and nf , respectively, and an initial time delay is controlled by nk . Here we consider them fixed a priori, but they can be determined as
well. The transfer function (4) has a corresponding impulse response denoted gθ . The relation between the parametric
model class Gθ and the family of all causal SISO systems G is illustrated in Figure 1.
Example 1 (Two members in the model classes). To illustrate the members of the model classes, suppose g? ∈ G is a
second order system with no zeros, two poles in z = 0.64 ± 0.48i and a static gain of 2. Each member of a given class
Gθ is specified by the parameters

θ = b0

f1

f2

f3

f4

>

Note that this model class is overparameterized with respect to the unknown system. A specific example of gθ ∈ Gθ is
illustrated in Figure 2 along with g? .

3

Decision-theoretic formulation

Given the dataset D, we want to identify a parametric model in Gθ that is as close to the unknown system as possible.
The choice of a specific model, specified by θ, is viewed as a decision with an associated loss
L(θ) =

1
g − gθ
2 ?
3

2
,
W

(5)

where g ? is given in (3), g θ is the corresponding impulse response of the parametric model and W  0 is a weight
matrix. In real applications, however, the system impulse response is unknown. In lieu of g ? , we consider an impulse
response g as a random variable drawn from G and average the loss (5) of a decision θ over all possible values of g,
a.k.a. risk:
1
1
R(θ) = Eg|D [L(θ)] = tr (WΣ) + kg − g θ k2W ,
(6)
2
2
where g and Σ are the mean and covariance matrix of g, respectively, given D. The proof of the equality is given in
Appendix A.1. The expected value g represents our best guess of g ? and therefore it is reasonable to choose W based
on the precision of this guess.
b that minimizes the risk, i.e.,
The optimal decision is the parameter θ

b = arg min R(θ).
θ

(7)

θ

This decision rule generates different identification methods based on model choice for g, which determines the mean
g and precision matrix W. The problem (7) can be solved using numerical search algorithms, cf. [11, ch. 7]. Note that
the loss (5) and resulting decision-rule can equivalently be formulated in the frequency domain.
Result 1. Suppose we let g be the unbiased least squares estimate of g ? and take W to be its precision matrix, i.e.,
g = H† y

and

W = λ−1 H> H.

(8)

Then (7) corresponds to PEM. Note that the variance λ does not affect the decision and can therefore be set to unity.
Proof. See Appendix A.2.
Example 2 (Making a decision). To identify the system, we generate in input u(t) as a Gaussian white process with
unit variance. The unknown system g? is the same as in the previous example and the unknown disturbance variance
is λ = 1. We also consider the same overparameterized model class Gθ .
Using the resulting input-output data D, it is straightforward to compute g and W in (8). Here we used a gradient
b in (7). The impulse response of the generating system g and the mean g are plotted in
search algorithm to find θ
?
Figure 3 along with the risk-minimizing decision g θb given by (7). The weighting W is visualized here as uncertainty
p
bands around g. Specifically, when the the matrix is full rank, we plot a band ±2 diag (W−1 ), where the square root
is computed element wise. The band then corresponds a dispersion of two standard deviations. Narrower uncertainty
bands correspond to higher weights and therefore forces the optimal decision g θb to be close to g at the corresponding
coefficients.
Result 2. Suppose g and the disturbances are modeled as Gaussian distributions, such that g ∼ N (0, K). Then
g = (λ−1 H> H + K−1 )−1 H> λ−1 y

and

W = K−1 + λ−1 H> H,

(9)

corresponds to the precision matrix of g when conditioned on D. Note that λ need not be set directly but can be
absorbed into the parametrization of K.
Proof. Let N (x|µ, Σ) denote the probability density function (pdf) of a multivariate Gaussian distribution with expected value µ and covariance matrix Σ. Since p(y | g) = N (y|Hg, λI) and p(g) = N (g|0, K), the distribution
over g is conjugate to the data distribution and the result follows in a straightforward manner [9].
Example 2. (cont’d) A priori we may expect the unknown system in G to be stable. To model this prior information,
suppose we set the entries of the covariance matrix K as
Ki,j = 100 · 0.8(i+j)/2 · 0.7|i−j|

(10)

In the subsequent section, we discuss means to automatically select K. The specific choice above places high prior
weight on systems in G that exhibit exponentially decaying impulse responses and weak positive correlation between
adjacent coefficients. Figure 4 shows the result of this choice. Compared to Figure 3, the uncertainty bands are now
tighter around the mean g, which also follows g ? more closely. Consequently, the risk-minimizing decision g θb is a
better approximation of g ? . The underlying reason is the regularization achieved by g and W.
4

p
g ± 2 diag (W−1 )

4

g?

g

g θb

2

0

−2
−4

0

2

4

6

8

10

12

14

16

18

t
Figure 3: Visualization of g (blue), the weighting W (gray) and the final decision gθb (red) together with
the data-generating impulse response g ? (black). Here g falls outside of the range of the figure.

p
g ± 2 diag (W−1 )

4

g?

g

g θb

2

0

−2
−4

0

2

4

6

8

10

12

14

16

18

t
Figure 4: Visualization of g (blue), the weighting matrix W (gray) and the final decision gθ̂ (red) together
with the data-generating impulse response g ? (black).

5

The example illustrates a bridge between classical frequentist and modern Bayesian approaches to system identification
using a decision-theoretic approach outlined above. Note also that unlike a standard Bayesian approach, there is no
need to directly specify a prior distribution over the parameters θ nor does the data generating g ? have to reside in the
considered model class Gθ . Via the covariance matrix K, the framework incorporates our prior information that the
system is probably stable. This, however, begs the question of how to model such information by K? This question is
discussed in the subsequent section.

4

Incorporating prior information

Prior information about the system is encoded in the prior distribution over the impulse response, e.g.,

g ∼ N 0, K(η) ,

where the hyperparameter η specifies the covariance matrix. Let K(η)i,j denote the ijth element of the matrix, which
corresponds to a covariance function or ‘kernel’ of the impulse response coefficients. A popular choice is

>
(11)
K(η)i,j = cα(i+j)/2 ρ|i−j| , η = c α ρ ,

which represents exponentially decaying or increasing impulse responses and is thus suitable for modelling e.g. a priori
information about stable systems. For other model choices, see [8].

The specification of hyperparameters η can be achieved in different ways. The Bayesian approach is to set them a
priori, reflecting personalized beliefs about the system prior to seeing any data. A more pragmatic approach is to tune
η using the data D. Such methods include cross-validation, the SURE-estimator, and maximum likelihood [1, 9, 4, 8].
Here we consider the maximum likelihood approach
b = arg max `(η),
η
η

(12)

where `(η) is obtained by marginalizing out g from the nonparametric data model. For Gaussian processes, we have
that
1
1
(13)
`(η) = − ln |λI + HKH> | − kyk2(λI+HKH> )−1
2
2
which only depends on K(η)i,j evaluated at the observed time instants [9].
Example 3 (Setting hyperparameters). Consider again the dataset D. In (10) we used the so-called diagonal correlated

>
(DC) kernel with hyperparameter η 0 = 100 0.8 0.7 . Using this as initial guess and maximizing the marginal

>
b = 0.45 0.71 0.74 . We can visualize the prior by studying the marginal variance of each
likelihood we obtain η
coefficient in the impulse response and plotting this as an uncertainty band together with the impulse response of the
generating system g ? . While both priors encompass g ? , the prior with tuned hyperparameters is much tighter as seen
in Figure 5. The resulting resulting decision is illustrated in Figure 6, which is to be compared with Figure 4. We
see that the tuned hyperparameters lead to a more accurate g and weights W, which consequently improves the final
decision. Thus we can avoid the manual tuning of hyperparameters after choosing the form of K(η) that model the
prior information about the system.

5

Numerical illustrations

In this section, we evaluate the proposed decision-theoretic approach, which we call Bayesian risk minimization (BRM)
when incorporating prior information in the form of a distribution on the impulse response g. For the prior covariance
matrix K, we consider a DC covariance function (11) and tune the hyperparameters using (13), cf. Section 4. The
method is summarized in Algorithm 1. The identified models are compared with those of PEM.
We consider the following data generating system
G? (q) =

1−

1.82q −1

0.41
.
+ 2.04q −2 − 1.27q −3 + 0.46q −4
6

p
±2 diag (K(η 0 ))

p
±2 diag (K(b
η ))

g?

4
2
0
−2
−4
0

2

4

6

8

10

12

14

16

18

t
Figure 5: A comparison of the prior distribution of g using the arbitrarily chosen η in the example (light
b (13) (dark shaded).
shaded) and using η

p
g ± 2 diag (W−1 )

1.5

g?

g

g θb

1
0.5
0
−0.5
−1

0

2

4

6

8

10

12

14

16

18

t
Figure 6: Same scenario as in Figure 4 but with a covariance matrix K(b
η ) that is automatically tuned
using (13).

Algorithm 1 Method for model class Gθ
1:
2:

3:
4:
5:
6:

Input: D = {(u(1), y(1)), . . . (u(N ), y(N ))}
Form matrix H.
Tune hyperparameters η.
Form g and W using (9).
Solve problem in (7).
b
Output: θ
7

log

kg ? − g θb k2
kg ? k2

2

0

−2

−4
PEM
BRM
N = 30

PEM
BRM
N = 60

PEM
BRM
N = 120

Figure 7: Box-plots of errors for PEM and BRM. Numbers indicate the sample size N .

This fourth-order system has two resonance peaks in the system frequency response. Similar to Example 2, we generate
the input as a white process with variance 1 and consider disturbances with variance 2. For each experiment we obtain
a dataset D of length N and evaluate the decision by the normalized squared error
log

kg ? − g θb k2
kg ? k2

in logarithmic units for clear comparisons. We repeat the experiments 100 times and study the distribution of errors.
In the first evaluation, we vary the size of the dataset using small, medium and large N . The model class Gθ under
consideration matches the form of the unknown system. The results are show in Figure 7. Note that for small sample
sizes, N = 30, PEM-based decisions frequently produces errors greater 0, which is higher than simply assuming an
impulse response with zero coefficients. The effect of regularization in BRM is significant and the performance by
PEM is only matched by doubling the sample size to N = 60. When N = 120, the median errors of both methods are
comparable, but BRM has a narrower dispersion than PEM.
In the second evaluation, we vary the order nf of the model class Gθ when N = 60. Specifically, we consider nf to
be lower, equal to, and greater than the unknown system model order. The results are illustrated in Figure 8, where it
is seen that for nf = 2, PEM and BRM perform similarly. As the overparameterization increases, PEM overfits and
produces many decisions with normalized errors above 0. By contrast, the robustness of BRM to overparameterization
is notable.

6

Conclusions

We developed a decision-theoretic approach that can be used to identify parametric linear time-invariant models. This
approach bridges the gap between the existing classical and regularized identification approaches. We showed that it
leads to a parametric identification method with a simple data-adaptive regularization by assigning a prior distribution
over the system impulse response. This incorporates prior information of the system and dispels both the need to
assume a prior distribution over the model parameters as well as the assumption that the unknown system belongs to
the model class. The numerical experiments illustrated that the regularization yields a method that is robust to both
small sample-sizes and overparameterization.

8

2

kg ? − g θb k2
kg ? k2

1
0

log

−1
−2
−3
−4
PEM BRM
nf = 2

PEM BRM
nf = 4

PEM BRM
nf = 6

PEM BRM
nf = 8

Figure 8: Box-plots of errors for PEM and BRM. Numbers indicate the model order nf . Here N = 60
samples are used. The order of the data generating system is 4.

In future work, we will study how the approach can be extended to a broader model class. We also consider the
equivalent formulation of the risk function in the frequency domain.

A
A.1

Appendix
Risk expression

Let us start by noting that
kg − g θ k2W = kg − g + g − g θ k2W = kg − gk2W + kg − g θ k2W + 2(g − g)> (g − g θ ),
for any vector g of appropriate size. Now, if we choose g = Eg|D [g], we obtain
"
#


2
2
2
>
Eg|D kg − g θ kW = Eg|D kg − gkW + kg − g θ kW + 2(g − g) (g − g θ ) = tr (WΣ) + kg − g θ k2W .

A.2

PEM as a decision

The objective function used in the prediction error method can be written as
b = arg min 1 y − y
eθ
θ
2
θ

2

,

e θ is the predicted output of the parametric system given the input. The predicted output is given by
where y
e θ = Hg θ .
y

9

(14)

(15)

Let g = H− y, where H− is a pseudoinverse of H. We get
b = arg min 1 ky − y
e θ k2
θ
2
θ
1
e θ k2
= arg min ky − Hg + Hg − y
2
θ
i
1h
= arg min ky − Hgk2 + 2(y − Hg)> H(g − g θ ) + kH(g − g θ )k2
2
θ
,
,
g = H− y =⇒ (y − Hg)> H = (y − HH− y)> H = 0

1
= arg min kg − g θ k2H> H .
2
θ
Assuming that the generating system is an output-error (OE) system, that is y = Hg + e, we have that [11]
E [g] = H− Hg

and

Cov(g) = λ(H> H)− .

Note that λ−1 H> H is a precision matrix and therefore a straight-forward choice for W. Thus the risk
R(θ) =

1
kg − g θ k2λ−1 H> H .
2

(16)

leads to PEM since the risk-minimizing decision is invariant with respect to λ.

Acknowledgement
The authors would like to thank Prof. Petre Stoica and Prof. Lennart Ljung for valuable comments.

References
[1] J.O. Berger. Statistical decision theory and Bayesian analysis. Springer Series in Statistics. Springer, 1985.
[2] Herman Chernoff and Lincoln E Moses. Elementary decision theory. Dover, 1986.
[3] L. Ljung. System Identification: Theory for the User. Pearson Education, 1998.
[4] Kevin P Murphy. Machine learning: a probabilistic perspective. MIT press, 2012.
[5] Brett Ninness and Soren Henriksen. Bayesian system identification via Markov chain Monte Carlo techniques.
Automatica, 46(1):40–51, 2010.
[6] Václav Peterka. Bayesian system identification. Automatica, 17(1):41–53, 1981.
[7] Gianluigi Pillonetto and Giuseppe De Nicolao. A new kernel-based approach for linear system identification.
Automatica, 46(1):81–93, 2010.
[8] Gianluigi Pillonetto, Francesco Dinuzzo, Tianshi Chen, Giuseppe De Nicolao, and Lennart Ljung. Kernel methods in system identification, machine learning and function estimation: A survey. Automatica, 50(3):657–682,
2014.
[9] Carl Edward Rasmussen and Christopher KI Williams. Gaussian processes for machine learning, volume 1. MIT
press Cambridge, 2006.
[10] Christian Robert. The Bayesian choice: from decision-theoretic foundations to computational implementation.
Springer Science & Business Media, 2007.
[11] T. Söderström and P. Stoica. System identification. Prentice-Hall, Inc., 1988.
10

