IEEE TRANSACTIONS ON SIGNAL PROCESSING

1

Asymptotic Analysis of Semi-Parametric
Weighted Null-Space Fitting Identification

arXiv:1708.03947v2 [cs.SY] 31 Aug 2017

Miguel Galrinho, Cristian R. Rojas, Member, IEEE, and Håkan Hjalmarsson, Fellow, IEEE

the number of parameters become very large), guaranteeing
that the noise spectrum is captured. If the global minimum
of the PEM cost function is found, in open loop this will
asymptotically not affect the dynamic-model estimates; in
closed loop, consistency is attained but not efficiency. The
problem is that, because the noise model might require many
parameters, the PEM cost function will potentially have more
local minima, thus aggravating the numerical search for the
global minimum.
Some methods use a semi-parametric approach, estimating a non-parametric noise model estimated in a first step,
which is then used in a second step to estimate the dynamic
model. Possible approaches have been delineated both in the
frequency [9]–[11] and time domains [12]–[14].
The weighted null-space fitting (WNSF) method [15] also
first estimates a non-parametric model, and then reduces it
to a full parametric model (i.e., dynamic and noise models).
Moreover, WNSF does not apply non-linear optimization
techniques, but uses weighted least squares iteratively. In this
Index Terms—System identification, least squares
sense, the method can be seen as belonging to the family
of iterative least-squares methods. These methods date back
to [16], and have later been denoted as iterative quadratic maxI. I NTRODUCTION
The prediction error method (PEM) is a benchmark for imum likelihood (IQML) methods, with applications to filter
estimation of linear parametric models. If the model orders design [17,18] and identification of dynamical systems [19]–
are correct and the noise is Gaussian, PEM with a quadratic [21]. The Steiglitz-McBride method [22] belongs also to this
cost function is asymptotically efficient [1]: the asymptotic class of methods, being equivalent to IQML for an impulsecovariance of the estimates is the Cramér-Rao (CR) bound— input case [23]. However, unlike in the aforementioned identification works, WNSF is asymptotically efficient in one iteration
the lowest covariance attainable by a consistent estimator.
Two models can typically be distinguished in a parametric in open and closed loop for Box-Jenkins systems [24].
Instead of simultaneously estimating the dynamic and noise
model structure: the dynamic model and the noise model.
models,
WNSF also allows to disregard the parametric noise
Because the noise sequence is often the result of different
model,
reducing
the non-parametric model estimate to obtain
noise contributions aggregated in a complex or even intractable
a
parametric
dynamic
model only. This guarantees asymptotimanner, the concept of a “correct order” for the noise model
cally
efficient
estimates
in open loop, and consistent estimates
is often inappropriate in practice. If the noise-model order is
in
closed
loop
with
optimal
asymptotic covariance when an
chosen too small with PEM, the dynamic-model estimate will
infinite-order
noise
model
is
used. This distinguishes this
be consistent in open loop, but biased in closed loop.
method
from
[13,14],
which
are
developed for open loop. The
The bias issue with closed-loop data is not exclusive of
asymptotic
properties
of
the
proposed
method correspond to
PEM. Instrumental variable methods [2] require the reference
the
asymptotic
properties
of
PEM
with
an infinite-order noise
signal to construct the instruments in closed loop [3,4]. For
model
[25]
in
both
open
and
closed
loop,
but performed with
classical subspace methods [5,6], the bias issue in closed loop
a
robust
numerical
procedure.
has been overcome by more recent algorithms [7,8].
In [24], an extensive study of the fully parametric WNSF
With PEM, the bias issue can in theory be solved by
method
has been conducted, including theoretical analysis and
letting the noise model structure be more flexible (i.e., letting
simulations. The case where WNSF is used with no parametric
Automatic Control Lab and ACCESS Linnaeus Center, School of Electrical
noise-model estimate—denoted semi-parametric WNSF—has
Engineering, KTH Royal Institute of Technology, SE-100 44 Stockholm, been addressed and illustrated in [24], but a formal proof
Sweden. (e-mail: {galrinho, crro, hjalmars}@kth.se.)
of the asymptotic properties has not been provided. In this
This work was supported by the Swedish Research Council under contracts
2015-05285 and 2016-06079.
paper, we perform the analysis of this setting, which has

Abstract—Standard system identification methods often provide biased estimates with closed-loop data. With the prediction
error method (PEM), the bias issue is solved by using a noise
model that is flexible enough to capture the noise spectrum.
However, a too flexible noise model (i.e., too many parameters)
can cause additional numerical problems for PEM. In this
paper, we perform a theoretical analysis of the weighted nullspace fitting (WNSF) method when a parametric noise model
is not estimated. With this method, the system is first captured
using a non-parametric ARX model, which is then reduced to a
parametric model of interest. In the reduction step, a noise model
does not need to be estimated if it is not of interest. In open loop,
this still provides asymptotically efficient estimates of the dynamic
model. In closed loop, the estimates are consistent, and their
covariance is optimal for a non-parametric noise model. In this
paper, we prove these results, which require additional technical
details compared with the case with a full parametric model
structure. In particular, we use a geometric approach for variance
analysis, deriving a new result that will be instrumental to our
end. Finally, we use a simulation study to illustrate the benefits
of the method when the noise model cannot be parametrized by
a low-order model.

IEEE TRANSACTIONS ON SIGNAL PROCESSING

2

additional difficulties because of the non-squared dimension The transfer functions Go , Ho , and Ho−1 are assumed to be
of some matrices. For this purpose, we use the geometric stable. The polynomials Lo and Fo —as well as Co and Do —do
approach of [26], extending the results therein to the setting not share common factors.
We parametrize Go (q) as
of our problem. Thus, our contributions are: first, we derive
a result for variance analysis in system identification; second,
l1 q −1 + · · · + lml q −ml
L(q, θ)
=
,
G(q, θ) =
we use this result to perform a theoretical analysis of the semiF (q, θ)
1 + f1 q −1 + · · · + fmf q −mf
parametric WNSF method, proving consistency and deriving
the asymptotic covariance of the estimates in open and closed where

⊤
θ = f 1 · · · f mf l 1 · · · l ml
loop. This analysis is complemented with a simulation study
to illustrate the benefits of the method.
is the parameter vector to estimate, with known orders mf
and ml . The noise model orders mc and md are not known.
II. P RELIMINARIES
Assumption 2 (Noise). The noise sequence {et } is a stochasEssentially, the same notation, definitions, and assumptions tic process that satisfies
used in [24] apply to this paper. For convenience, we reproduce
E[et |Ft−1 ] = 0, E[e2t |Ft−1 ] = σo2 , E[|et |10 ] ≤ C, ∀t.,
the notation and assumptions here; for the definitions of fN stability and fN -quasi-stationarity, see [24,27].
Assumption 3 (Input). The input sequence {ut } is defined by
ut = −K(q)yt + rt under the following conditions.
1) The sequence {rt } is assumed to q
be independent of {et },
A. Notation
qP
log N
fN -quasi-stationary with fN =
n
2
N , and uniformly
• kxk =
k=1 |xk | , with x an n × 1 vector.
bounded: |rt | ≤ C, ∀t.
• kAk = supx6=0 kAxk / kxk, with A a matrix and x a
2) Let the spectral factorization of the power spectral density
vector of p
appropriate dimensions.
of {rt } be Φr (z) = Fr (z)Fr (z −1 ). Then, Fr (q) is
Trace(AA∗ ) (i.e., the Frobenius norm), with
• ||A||F =
assumed to be fN -stable, with fN = 1.
√
A a matrix.
q R
3)
The closed loop system is fN -stable with fN = N .
π
1
iω )Ḡ∗ (eiω )dω, with
Trace
Ḡ(e
• Ḡ(q)
:=
2π −π
H2
4) The feedback transfer function K(z) is bounded on the
Ḡ(q) a transfer matrix.
unit circle.
• Ḡ(q)
:= supω ||Ḡ(eiω )||.
5) The spectral density of the process {[rt et ]⊤ } is bounded
H∞
• C denotes any constant, which need not be the same in
from below by the matrix δI, δ > 0 (this implies an
different expressions.
informative experiment).
−1
• Γn (q) = [q
· · · q −n ]⊤ , where q −1 is the backward
Operation in open loop is obtained by taking K(q) = 0.
time-shift operator.
∗
Alternatively to (1), the true system can be written as
• A is the complex conjugate transpose of the matrix A.
• Tn,m (X(q)) is the lower-triangular Toeplitz matrix of size
Ao (q)yt = Bo (q)ut + et ,
n × m (mP
≤ n) with first column [x0 · · · xn−1 ]⊤ , where
k
X(q) = ∞
k=0 xk q . The dimension n may be infinity, where
∞
X
denoted T∞,m (X(q)).
1
A
(q)
:=
=:
1+
aok q k ,
o
• Ex denotes expectation of the random vector x.
H
(q)
o
P
N
k=1
1
• Ēxt := lim N
(3)
t=1 Ext .
∞
N →∞
X
Go (q)
o
k
• xN = O(fN ) means that the function xN tends to zero
Bo (q) :=
=:
bk q
Ho (q)
at a rate not slower than
f
,
as
N
→
∞,
w.p.1.
N
k=1
Rπ
1
iω
∗ iω
• hX(q), Y (q)i := 2π −π X(e )Y (e )dω, where X(q)
are stable (Assumption 1). In an intermediate step, WNSF
and Y (q) are transfer matrices of appropriate sizes.
estimates truncated A (q) and B (q), using the ARX model
o

o

n

A(q, η )yt = B(q, η n )ut + et ,

B. Assumptions
Assumption 1 (True system and parametric model). The
system has scalar input {ut }, scalar output {yt }, and is
subject to scalar noise {et }. These signals are related by
yt = Go (q)ut + Ho (q)et ,

(1)

where Go (q) and Ho (q) are rational functions according to
Go (q) =

o
l1o q −1 + · · · + lm
q −ml
Lo (q)
l
=
,
o
−1
o
Fo (q)
1 + f1 q + · · · + fmf q −mf

1 + co1 q −1 + · · · + comc q −mc
Co (q)
=
.
Ho (q) =
Do (q)
1 + do1 q −1 + · · · + domd q −md

(2)

(4)

where
⊤

(5)
η n = a1 · · · an b 1 · · · b n ,
n
n
X
X
A(q, η n ) = 1 +
ak q −k , B(q, η n ) =
bk q −k .
k=1

k=1

Because the order needs to be infinite for the system to be
in the ARX model set, we make the model order n depend
on the sample size N —denoted n = n(N )—according to the
following assumption.
Assumption 4 (ARX-model order). The ARX model order is
selected according to:

IEEE TRANSACTIONS ON SIGNAL PROCESSING

3

D1. n(N ) → ∞, as N → ∞;
D2. n4+δ (N )/N → 0, for some δ > 0, as N → ∞.

III. S EMI -PARAMETRIC W EIGHTED N ULL -S PACE F ITTING

C. Prediction Error Method
The prediction error method minimizes a cost function of
the prediction errors


L(q, θ)
−1
(6)
ut ,
εt (θ, ζ) = H (q, ζ) yt −
F (q, θ)
where H(q, ζ) is some noise model parametrization, function
of a parameter vector ζ. Using a quadratic cost function, which
is optimal when the noise sequence is Gaussian, the PEM
estimate of the parameters is obtained by minimizing
N
1 X1 2
J(θ, ζ) =
ε (θ, ζ),
N t=1 2 t

(7)

where N is the sample size.
Let H(q, ζ) be such that there exists ζ = ζo for which
PEM
H(q, ζo ) = Ho (q). Denoting by θ̂N
the parameter vector θ
PEM
that (together with some ζ) minimizes (7), the estimate θ̂N
is asymptotically distributed as [1]
√
−1
PEM
N (θ̂N
− θo ) ∼ AsN (0, σo2 MPEM
),
where N stands for the Gaussian distribution. The covariance
matrix satisfies
Z π
1
Ω̄(eiω )Φru (eiω )Ω̄∗ (eiω )dω =: M,
(8)
MPEM ≥
2π −π
with (arguments often omitted for notational simplicity)

 Go
− Ho Fo Γmf
Ω̄ =
1
Ho Fo Γml
= So (q)rt ,

yt = (ϕnt )⊤ η n + et ,

ϕnt = −yt−1

· · · −yt−n

ut−1

· · · ut−n

(9)

⊤

.

Then, the least-squares estimate of η n is obtained by
n
n −1 n
η̂N
= [RN
] rN ,

(10)

where
n
RN
=

N
1 X n n ⊤
ϕ (ϕ ) ,
N t=n+1 t t

n
rN
=

N
1 X n
ϕ yt .
N t=n+1 t

(11)

The asymptotic distribution of the estimates is [27]

√
 −1 
n
,
− η̄ n ) ∼ AsN 0, σo2 R̄n
N (η̂N

where



n
RN
→ R̄n := Ē ϕnt (ϕnt )⊤ , as N → ∞, w.p.1,
n
rN
→ r̄n := Ē [ϕnt yt ] ,
 −1 n
n
r̄ ,
η̂N
→ η̄ n := R̄n

as N → ∞, w.p.1,

as N → ∞, w.p.1.

For the second step, we obtain an estimate of G(q, θ), from
the non-parametric ARX-model estimate. For this purpose, we
use (2) and (3) to write
Fo (q)Bo (q) − Lo (q)Ao (q) = 0.

and Φru the spectrum of
urt

The WNSF method consists of three steps [24]. First, we
estimate a non-parametric ARX model, with least squares.
Second, we reduce this estimate to a parametric model, with
least squares. Third, we re-estimate the parametric model, with
weighted least squares. We now detail the procedure for each
step, without estimating a parametric noise model.
For the first step, consider (4) in the regression form

(12)

Because we are not interested in estimating a parametric noise
model, we do not consider the part of (3) from where such a
model could be obtained. By convolution, (12) can be written
in matrix form as

where So (q) = [1 + K(q)Go (q)]−1 is the sensitivity function.
In open loop (in which case Φru is simply the input specbno − Qn (ηon )θo = 0,
(13)
trum), MPEM = M , and it corresponds to the CR bound. In
closed loop, MPEM = M when the number of parameters in ζ where ηon is given by (5) evaluated at the true coefficients
tends to infinity. In this case, M does not correspond to the CR of (3), bno consists of the last n coefficients of ηon , and
bound, but to the best possible covariance (from a prediction


Qn (η n ) = −Qfn (η n ) Qln (η n ) ,
error perspective) with a non-parametric noise model [25].
Qln (η n ) = Tn,ml (A(q, η n )), Qfn (η n ) = Tn,mf (B(q, η n )).
The interest of estimating a non-parametric noise model in
(14)
closed loop is that even if the noise spectrum needs to be
captured by a high-order model, it will still be possible to
n
Motivated by (13), we replace ηon by its estimate η̂N
(and the
obtain an estimate of the dynamic model G(q, θ). However,
same for bno , which is a part of ηon ) and obtain an estimate of
estimating a non-parametric noise model simultaneously with
θ with least squares:
a parametric dynamic model with PEM is not realistic. The rea−1 ⊤ n n
LS
n
n
son is that, as the number of parameters in H(q, ζ) increases,
θ̂N
= Q⊤
Qn (η̂N )b̂N .
(15)
n (η̂N )Qn (η̂N )
the prediction error (6) becomes a more complicated function
For the third step, we re-estimate θ taking into account the
of ζ, and finding the global minimum of the non-convex
n
n
errors
in η̂N
. As ηon is replaced by η̂N
in (13), the residuals
cost function (7) is a more difficult task. Consequently, the
n
n
n
n
are
given
by
b̂
−
Q
(η̂
)θ
=
T
(θ
)(η̂
n N o
n o
result that PEM with a non-parametric noise model provides
N
N − ηo ), where


estimates with covariance corresponding to M may not always
Tn (θ) = −Tnl (θ) Tnf (θ) ,
be useful in practice. With WNSF, this setting can be handled
Tnl (θ) = Tn,n (L(q, θ)), Tnf (θ) = Tn,n (F (q, θ)). (16)
without increasing the difficulty of the problem.

IEEE TRANSACTIONS ON SIGNAL PROCESSING

4

The variance of the estimate is then minimized by solving
a weighted least-squares problem, where the weighting is
the inverse of the covariance of the residuals. Although this
covariance is dependent on the true parameters, a consistent
estimate is available by
LS
LS
LS
) = Tn (θ̂N
)[Rn ]−1 Tn⊤ (θ̂N
),
Wn−1 (θ̂N

(17)

and the estimate obtained in this step is thus given by
−1

n
LS
n
WLS
n
LS n
)W
(
θ̂
)Q
(η̂
)
θ̂N
= Q⊤
(η̂
Q⊤
n N
n N
n N
n (η̂N )Wn (θ̂N )b̂N .
(18)
The algorithm may be summarized as follows.
Algorithm 1. The semi-parametric WNSF method consists of
the following steps:
1) compute a non-parametric ARX-model estimate with (10);
2) compute a parametric dynamic-model estimate with (15);
3) re-compute a parametric dynamic-model estimate
with (18).

n(N )

n
n
,
where M (η n , θ) := Q⊤
n (η )Wn (θ)Qn (η ) and η̂N := η̂N
recalling that n is a function of N according to Assumption 4 (for notational simplicity, we use only n instead
of n(N ) in matrix subscripts). To study consistency and
asymptotic distribution,
limit value of (19) and the asymp√ the WLS
− θo ) will be analyzed. The
totic distribution of N (θ̂N
challenge in this analysis (compared to [24]) is in how to
LS
treat the matrix Wn (θ̂N
), and consequently also the matrix
LS
LS
LS
M (η̂N , θ̂N ), which contains Wn (θ̂N
). From (17), Wn (θ̂N
)=
LS
n ⊤ LS −1
[Tn (θ̂N )R̄ Tn (θ̂N )] , which in [24] may be written as
LS
LS
LS
LS
Wn (θ̂N
) = Tn−⊤ (θ̂N
)R̄n Tn−1 (θ̂N
), as the matrix Tn (θ̂N
)
is square therein. However, in the semi-parametric version
LS
treated in this paper, Tn (θ̂N
) is not square, and we cannot
LS
analyze Wn (θ̂N ) by taking inverses of the individual matrices
it consists of.
To deal with this issue, we use the approach in [26], writing
the aforementioned matrices as projections of the rows of
some matrix onto the subspace spanned by the rows of another
matrix. This will be applied to the limit value of the matrix
LS
M (η̂N , θ̂N
), defined by

Optionally, we may continue to iterate, potentially improving the estimation quality for finite sample size. However,
M̄ (ηo , θo ) := lim Q⊤
(η n )[Tn (θo )R̄n Tn⊤ (θo )]−1 Qn (ηon ).
n→∞ n o
we show in the next section that Algorithm 1 has the same
(20)
asymptotic properties as PEM with an infinite-order noise Writing Qn (η n ), Tn (θo ), and R̄n (defined in (14), (16),
o
model. However, WNSF estimates the non-parametric noise and (11), respectively) in the frequency domain as
model in a separate step, as part of an ARX model, which


Z π 
−Bo ∗
Γm 0
is linear in the model parameters. Thus, it does not make Q⊤ (η n ) = 1
Γn dω,
n o
Ao
2π −π 0 Γm
the problem computationally more difficult, unlike if (7) is
 ∗
 ∗

Z π
 −Lo 0

1
minimized with an arbitrary large number of noise model
Γn 0
Γn Γn
dω,
Tn (θo ) =
parameters.
0
Fo∗
0 Γ∗n
2π −π


Z π 
1
−Go So Fr −Ho σo
Γn 0
R̄n =
IV. T HEORETICAL A NALYSIS
So Fr
0
2π −π 0 Γn
∗  ∗

Γn 0
−G
S
F
−H
σ
o o r
o o
In this section, we perform a theoretical analysis of the semi·
dω,
So Fr
0
0 Γ∗n
parametric WNSF method. In particular, we will show that
Step 3 in Algorithm 1 provides a consistent and asymptotically we may express M̄ (ηo , θo ) as
efficient estimate. For that, we will need some auxiliary results.
M̄ (ηo , θo ) =
A. Auxiliary Results
To show the aforementioned results, we will need that the
estimate obtained in Step 2 of Algorithm 1 is consistent. For
that, we have the following result.
LS
Theorem 1. Let Assumptions 1, 2, 3, and 4 hold, and θ̂N
be
given by (15). Then,
LS
θ̂N
→ θo , as N → ∞, w.p.1.

Proof. The result is a specific case of [24, Theorem 1].
Although consistency of Step 2 with semi-parametric WNSF
is a specific case of the fully-parametric results [24], consistency and asymptotic distribution of Step 3 require additional
results. To get further insight into why, we begin by writing,
for the estimate from Step 3,
WLS
θ̂N

− θo

LS
LS
n(N )
= M −1 (η̂N , θ̂N
)Q⊤
),
n (η̂N )Wn (θ̂N )Tn (θo )(η̂N − ηo
(19)

lim hγ, Ψn i[hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψn i]−1 hΨn , γi,

n→∞

(21)

where



−Go (q)So (q)Γn (q) −σo Ho (q)Γn (q)
,
So (q)Fr (q)Γn (q)
0

F ∗ (q) 
Γn (q) 0 ,
Ψn = ∗ o ∗
So (q)Fr (q)


So (q)Fr (q) −Bo (q)Γmf (q) 0
γ=
.
Ao (q)Γml (q) 0
Fo (q)
Ωn =

(22)

Following the approach in [26], we recognize that the term
hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψn i in (21) can be written as
hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψn i = hProjSΩn Ψn , ProjSΩn Ψn i,
where ProjSΩn Ψn denotes the projection of the rows of Ψn
onto the subspace spanned by the rows of Ωn . As n → ∞,
the dimensions of the matrix Ωn increase, and the subspace
spanned by its rows approaches H2 . Then, the limit value of
the projection will be the causal part of the projected matrix.

IEEE TRANSACTIONS ON SIGNAL PROCESSING

5

For a simplified case, suppose that Ψn were causal and
V. S IMULATIONS
that its dimension did not depend on n (i.e., Ψn = Ψ). In
In this section, we perform two simulation studies. In the
this case, we would have limn→∞ hProjSΩn Ψ, ProjSΩn Ψi = first, we illustrate the asymptotic properties of the method.
hΨ, Ψi. In turn, we would then have that M̄ (ηo , θo ) = In the second, we illustrate how estimating a non-parametric
hγ, ΨihΨ, ΨihΨ, γi = hProjSΨ γ, ProjSΨ γi. If we now reintro- noise model with WNSF may be useful in scenarios where a
duce that the dimension of Ψ depends on n (Ψ = Ψn ), and we low-order parametrization for the noise model does not capture
assume that the rows of Ψn span H2 as n → ∞, we have that the noise spectrum accurately enough.
M̄ (ηo , θo ) = limn→∞ hProjSΨn γ, ProjSΨn γi = hγ, γi. These
arguments follow from results in [26]. However, handling the
A. Illustration of asymptotic properties
dimensional increase of Ψn with n requires additional techniAccording to the results in Section IV, semi-parametric
cal developments. One of the key results for the asymptotic
WNSF
is asymptotically efficient in open loop, with the
analysis in this paper is that the aforementioned result (i.e., that
asymptotic
covariance of the dynamic-model estimates given
M̄ (ηo , θo ) = hγ, γi) still holds when the dimensions of Ψn
2
−1
by
σ
M
.
In closed loop, the asymptotic covariance is still
o
increase with n. This is considered in the following theorem.
given by σo2 M −1 , but in this case it does not correspond to
Theorem 2. Let
the CR bound, but to the optimal asymptotic covariance when


the noise-model order tends to infinity.
F1 (q)Γn (q) F2 (q)Γn (q)
Ωn =
,
To illustrate this, we perform open- and closed-loop simuF3 (q)Γn (q)
0


lations
such that the closed-loop data are generated by
 ∗

F (q)Γmf (q) 0
Ψn = F4 (q)Γn (q) 0 , γ = 5
,
F6 (q)Γml (q) 0
K(q)Ho (q)
1
rt −
et ,
ut =
P∞ (j) −k
1 + K(q)Go (q)
1 + K(q)Go (q)
where Fj (q) =
f
q
(j
=
{1,
...,
6})
are
expok=0 k
Go (q)
Ho (q)
(j)
(4)
yt =
rt +
et ,
nentially stable (i.e., |fk | < Cλk ∀j, λ < 1), f0 6= 0,
1 + K(q)Go (q)
1 + K(q)Go (q)
−1
−1
and F3 (q) and F4 (q) are exponentially stable. Then, if
and the open-loop data by
||hΩn , Ωn i−1 || < C ∀n,
1
rt ,
1 + K(q)Go (q)
yt = Go (q)ut + Ho (q)et ,

ut =

lim hγ, Ψn i[hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψn i]−1 hΨn , γi

n→∞

= hγ, γi. (23)
Proof. See Appendix B.
B. Consistency and Asymptotic Efficiency
Using the auxiliary results derived above, we now show
that semi-parametric WNSF is consistent and asymptotically
efficient. Regarding consistency of Step 3 in Algorithm 1, we
have the following result.
Theorem 3. Let Assumptions 1, 2, 3, and 4 hold. Then,
WLS
→ θo , as N → ∞, w.p.1.
θ̂N

Proof. See Appendix C.
Theorem 3 implies that semi-parametric WNSF provides consistent estimates of θo .
Regarding the asymptotic distribution and covariance of
Step 3 in Algorithm 1, we have the following result.
Theorem 4. Let Assumptions 1, 2, 3, and 4 hold. Then,
√
WLS
N (θ̂N
− θo ) ∼ AsN (0, σo2 M −1 ),
where M is given by (8).
Proof. See Appendix D.
As consequence of Theorem 4, the semi-parametric WNSF
method summarized in Algorithm 1 is asymptotically efficient
in open loop, with M corresponding to the CR bound. In
closed loop, M corresponds to the best covariance attainable
by PEM with an infinite-order noise model [25].

where {rt } and {et } are independent Gaussian white sequences with unit variance, K(q) = 1, and
Go (q) =

q −1 + 0.1q −2
,
1 − 0.5q −1 + 0.75q −2

Ho (q) =

1 + 0.7q −1
.
1 − 0.9q −1

We perform 1000 Monte Carlo runs, with sample sizes N ∈
{300, 600, 1000, 3000, 6000, 10000}. We apply WNSF with an
ARX model of order 50 with the open- and closed-loop data.
Performance is evaluated by the mean-squared error of the
estimated parameter vector of the dynamic model,
WLS
MSE = ||θ̂N
− θ̄o ||2 ,

As this simulation has the purpose of illustrating asymptotic
properties, initial conditions are assumed known—that is, the
sums in (11) start at t = 1 instead of t = n + 1.
The results are presented in Fig. 1, with the average MSE
plotted as function of the sample size (closed loop in solid
line, open loop in dashed line). We plot also σo2 Trace[M −1 ]/N
(dotted line), which the average MSE attains both in open and
closed loop: because the data were generated such that Φru , the
spectrum of (9), is the same for both data sets, both scenarios
have the same asymptotic covariance, in accordance to our
theoretical results.
B. Random noise model
When a low-order parametrization of the noise model is
not enough to capture the noise spectrum, the noise model
may require many parameters. With PEM, a simultaneous

IEEE TRANSACTIONS ON SIGNAL PROCESSING

6

The alternative is to use a high-order model for the noise
model. For example,

10−2

MSE

H(q, ζ n ) = 1 +

n
X

ζk q −k ,

(26)

k=1

10

−3

or
H(q, ζ n ) =

10

−4

103

N

104

Fig. 1. Illustration of asymptotic properties: theoretical asymptotic MSE
(dotted) and average MSE for the parameter estimates as function of sample
size obtained with semi-parametric WNSF in closed loop (solid) and open
loop (dashed).

1+

where

ζ n = ζ1

1
Pn

k=1 ζk q

· · · ζn

⊤

−k

,

(27)

.

In particular, the choice (27) has the same structure as the
noise model estimated by semi-parametric WNSF.
Motivated by these alternatives, we compare the following
methods:
• semi-parametric WNSF, with non-parametric model order
estimate of the dynamic model and a long noise model is
n = 200;
not numerically robust due to the non-convexity of the cost
• PEM, with default MATLAB initialization, and noise
function. The semi-parametric WNSF is appropriate to deal
model (25) with mh ∈ {1, 2, ...30}, where the order mh
with this scenario, because the noise spectrum is captured
is chosen using the AIC or BIC criterion (denoted PEMaic
beforehand with a non-parametric ARX model.
and PEMbic , respectively);
Modeling the correct noise spectrum is particularly impor• PEM, with default MATLAB initialization, and noise
tant in closed loop, where the estimates of the plant will be
model (27) with n = 200 (denoted PEMnp , where ‘np’
biased if the noise model is not flexible enough to capture
stands
for non-parametric).
the noise spectrum. For this reason, we consider a closed-loop
PEM uses the implementation in MATLAB2016b System
setting, where data are generated by
Identification Toolbox. All the methods use a maximum of 100
1
K(q)Ho (q)
iterations, but stop early upon convergence (default settings for
ut =
rt −
et
1 + K(q)Go (q)
1 + K(q)Go (q)
PEM, 10−4 as tolerance for the normalized relative change in
Ho (q)
Go (q)
the parameter estimates for WNSF). The search algorithm used
rt +
et .
yt =
1 + K(q)Go (q)
1 + K(q)Go (q)
by PEM is chosen automatically. The noise model (26) was not
used with PEM for computational reasons: the optimization
The signals {rt } and {et } are Gaussian white noise sequences
becomes extremely slow as stability of the inverse of the noise
with variances 1 and 4, respectively. The system is given by
model when computing the prediction errors (6) is difficult
1.0q −1 − 0.80q −2
to fulfill with so many parameters, while the inverse of any
Go (q) =
,
estimate of (27) is always stable.
1 − 0.95q −1 + 0.90q −2
We use sample sizes N ∈ {1000, 5000, 10000} and perform
the controller by K(q) = 0.2, and the true noise model by
100 Monte Carlo runs. Performance is evaluated by the FIT
N
−1
of the impulse response of the dynamic model, given by
X


Ho (q) = 1 +
λk q −k
(24)
kgo − ĝk
k=1
FIT = 100 1 −
,
(28)
kgo − mean(go )k
−0.2k
with λk = wk e
, where wk is drawn from a Gaussian
distribution with zero mean and unit variance. Here, differently in percent, where go is a vector with the impulse response
than Assumption 1, stability of Ho (q) is not ensured. However, parameters of Go (q) (mean(go ) is its average), and similarly
this is not an issue if the noise is Gaussian, as there must exists for ĝ but for the estimated model. In (28), sufficiently long
an inversely stable Ho (q) for which the noise sequence has the impulse responses are taken to make sure that the truncation
of their tails does not affect the FIT.
same spectrum.
The FITs for the different sample sizes are shown in Fig. 2.
It may be difficult to find an appropriate low-order
For
N = 1000, WNSF has the most robust performance, with
parametrization for the noise model (24) in the form
smaller whiskers than the remaining algorithms and no lowc1 q −1 + · · · + cmh q −mh
H(q, ζ) =
,
(25) performance occurrences, which occur when PEM is used with
d1 q −1 + · · · + dmh q −mh
AIC/BIC or with a a non-parametric noise model. Among the
PEM alternatives, an AIC/BIC criterion with a Box-Jenkins
where

⊤
model with orders up to 30 performed better than using a nonζ = c1 · · · cmh d1 · · · dmh .
parametric noise model. For N = 5000, WNSF and PEM with
In this case, a Box-Jenkins model is estimated. The most non-parametric noise model have similar median performance,
appropriate noise model order may be chosen by using some put the algorithm for PEM failed more often. Here, PEM with
information criterion, such as the Akaike Information Criterion AIC/BIC had no low-performance outliers, but the median
performance was poorer than for WNSF and PEMnp . Similar
(AIC) or the Bayesian Information Criterion (BIC) [1].

IEEE TRANSACTIONS ON SIGNAL PROCESSING

N = 1000

7

N = 5000

N = 10000

100

100

50

80

80

60

60

FIT

100

0

↓1

↓3 ↓4 ↓8

↓4

PEMaic

WNSF

↓4 ↓5 ↓4

PEMbic

PEMnp

Fig. 2. FITs for given methods and sample sizes with 100 Monte Carlo runs.
TABLE I
AVERAGE COMPUTATIONAL TIMES IN SECONDS FOR WNSF, THE SEARCH
AMONG ALL ORDERS REQUIRED FOR PEM AIC AND PEM BIC , AND PEM NP .
N

1000

5000

10000

WNSF
PEMaic,bic
PEMnp

1.03
30.3
641

0.907
26.8
133

1.29
38.1
236

conclusions can be drawn for N = 10000, where we observe
that PEM does not necessarily provide better results with more
data samples, potentially due to numerical problems.
Overall, WNSF showed more robust performance among the
sample sizes used. However, an even more evident advantage
is the computational time. Table I shows the average times,
in seconds, required for the identification using WNSF, PEM
with all the orders computed for AIC and BIC, and PEMnp ,
for the different sample sizes (all the computations were
performed in the same computer). Here, we observe how
WNSF requires much lower computational time than the
alternatives. This is a consequence of PEM estimating the
noise model in the non-linear optimization procedure, while
in WNSF the high-order model is estimated in a previous
least-squares step, which is numerically robust. Moreover, the
time required for WNSF and PEM with AIC/BIC does not
change significantly with N , unlike with PEMnp . In this case,
the time does not necessarily decrease for smaller N . The
problem arising when using smaller N is that the cost function
more likely becomes ill-conditioned at some parameter values
during the optimization.
VI. C ONCLUSION

order noise model is required, both in terms of performance
and computational time. With WNSF, this separation always
occurs, as the method first estimates a high-order ARX model.
Then, a low-order noise model does not need to be obtained,
as the noise spectrum has been captured in the first step.
We also provide a theoretical analysis of the asymptotic
properties. To this end, we extend the geometric approach
in [26], deriving a more general result that may also be
useful for variance analysis of other methods. We show that
semi-parametric WNSF provides consistent estimates of the
dynamic model with closed-loop data. With open-loop data,
the estimates are also asymptotically efficient; with closedloop data, the asymptotic covariance of the estimates corresponds to the best possible covariance with a non-parametric
noise model. This gives WNSF attractive features in terms of
flexibility of noise-model structure and asymptotic properties:
if a correct parametric noise model is estimated, the dynamicmodel estimates are asymptotically efficient; if not, they are
consistent and optimal for a non-parametric noise model. We
used a simulation study to illustrate how semi-parametric
WNSF is an appropriate method for scenarios where the
noise model cannot be accurately modeled with a low-order
parametrization.
A PPENDIX A
AUXILIARY R ESULTS FOR THE P ROOF OF T HEOREM 2
The following results will be used to prove Theorem 2.
Extension of Cauchy-Schwarz inequality for transfer-matrix
inner products: Let X(q) and Y (q) be transfer matrices and
x and y be vectors of appropriate dimensions. Then, we have
Z π
1
x⊤ X(eiω )Y ∗ (eiω )ydω
||hX, Y i||2 =
sup
||x||=1,||y||=1 2π −π
q
x⊤ hX, Xix y ⊤ hY, Y iy
≤
sup
||x||=1,||y||=1

= ||hX, Xi|| ||hY, Y i||.

(29)
Bound for spectral norm of inner product of transfer matrices: Let X(q) be a transfer matrix. Then, we have
p
||hX, Xi|| ≤ ||hX, Xi||F = Trace(hX, Xi2 )
(30)
≤ TracehX, Xi = ||X||2H2 ,

where the second inequality follows from Trace(A2 ) ≤
[Trace(A)]2 for a positive semi-definite matrix A.
for Toeplitz operators of stable filters: Let X(q) :=
PBound
∞
k
x
q
−∞ k . From [28, Chapter 3], we have
Z π
1
Γ∞ Γ∗∞ X(eiω )dω = ||X(q)||H∞
2π −π
Z π
1
≥
Γn Γ∗n X(eiω ) , (31)
2π −π

Many standard system identification methods provide biased
estimates with closed-loop data. In the particular case of PEM,
the bias issue is avoided by choosing a noise-model order that
is large enough to capture the noise spectrum. An appropriate
order is often difficult to choose, and making it arbitrarily where the inequality is due to the matrix on the right-hand
large increases the numerical problems of PEM. In [24], where side of (31) being a block of the infinite-dimensional matrix
WNSF is analyzed, it is stated that the method can be used on the left-hand side.
without a parametric noise-model estimate, named the semiA PPENDIX B
parametric WNSF. In this paper, we deepen this discussion.
P
ROOF
OF T HEOREM 2
A simulation study illustrates the importance of separating
the dynamic- and noise-model identification when a highIn this appendix, we prove Theorem 2.

IEEE TRANSACTIONS ON SIGNAL PROCESSING

8

Inner Projection: hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψn i

Replacing (37) in (23), we obtain

Let

Ψn = Ψcn + Ψan ,

(32)

lim hγ, Ψn i[hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψn i]−1 hΨn , γi

n→∞

= lim hγ, Ψn i[hΨcn , Ψcn i + ∆n ]−1 hΨn , γi, (38)

where

n→∞

(4)
f0 e−iω
(4)
(4) −2iω
+ f1 e−iω
f0 e





:= 
..


.
Pn−1 (4) −kiω
(4) −niω
f0 e
+ k=1 fn−k e

 P∞ (4)
iωk
0
k=0 fk+1 e

 P∞ (4) iωk
0
 k=0 fk+2 e
a

Ψn := 
..
.. 
.

.
.
P∞ (4) iωk
0
k=0 fk+n e

Alternatively, Ψcn can be written as Ψcn = Pn Γn
Pn = Tn×n (F4 (q)). Using (32), we can write
Ψcn

and using the Sherman-Morrison-Woodbury formula, we rewrite (38) as


0

0
.. 
,
.
0

−1

lim hγ, Ψn i[hΨcn , Ψcn i + ∆n ]−1 hΨn , γi

n→∞

= lim hγ, Ψn ihΨcn , Ψcn i−1 hΨn , γi
n→∞

+ lim hγ, Ψn ihΨcn , Ψcn i−1 ∆n
n→∞

· [I + hΨcn , Ψcn i−1 ∆n ]−1 hΨcn , Ψcn i−1 hΨn , γi.

We want to show that the second term on the right-hand
side of (39), for which we can write

0 , where

hΩn , Ψn i
c
= hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψcn i
+ hΨan , Ωn ihΩn , Ωn i−1 hΩn , Ψcn i
+ hΨcn , Ωn ihΩn , Ωn i−1 hΩn , Ψan i
+ hΨan , Ωn ihΩn , Ωn i−1 hΩn , Ψan i
= hΨcn , Ωn ihΩn , Ωn i−1 hΩn , Ψcn i

||hγ, Ψn ihΨcn , Ψcn i−1 ∆n [I + hΨcn , Ψcn i∆n ]−1
· hΨcn , Ψcn i−1 hΨn , γi||

≤ ||hγ, Ψn i||2 ||hΨcn , Ψcn i−1 ||2
· ||[I + hΨcn , Ψcn i−1 ∆n ]−1 || ||∆n ||,

hΨn , Ωn ihΩn , Ωn i

hΨan , Ωn i

(33)

because
=0=
∀n, as Ωn is causal and
Ψan is anti-causal.
Now, we will construct an approximation of Ψcn whose
rows can be written using a linear combination of the rows
−1 c
−1 c
hofPΩn . Let row i iof F3 Ψn be given by F3 Ψn (i) =:
(i)
(i) −k
∞
0 , where |βk | ≤ Cλk with λ < 1, and
k=1 βk hq
i
Pn
(i)
−k
Ψ̂cn (i) :=
0 be the row i of a matrix
k=1 βk F3 q
Ψ̂cn (i). We re-write the right side of (33) as

hΨcn , Ωn ihΩn , Ωn i−1 hΩn , Ψcn i = Λn + ∆(1)
n ,

(34)

where
c
c
−1
∆(1)
hΩn , Ψcn i
n =hΨn − Ψ̂n , Ωn ihΩn , Ωn i

+ hΨcn , Ωn ihΩn , Ωn i−1 hΩn , Ψcn − Ψ̂cn i

+ hΨcn − Ψ̂cn , Ωn ihΩn , Ωn i−1 hΩn , Ψcn − Ψ̂cn i

By construction, Ψ̂cn (i) is a linear combination of the last n
rows of Ωn ; therefore, Ψ̂cn ∈ SΩn and
(35)

Using (33), (34), and (35), we write
hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψn i = hΨ̂cn , Ψ̂cn i + ∆(1)
n .

(36)

Moreover, we can re-write (36) as
hΨn , Ωn ihΩn , Ωn i
where ∆n =

(1)
∆n

+

hΩn , Ψn i =

(2)
∆n

hΨcn , Ψcn i

+ ∆n ,

tends to zero as n tends to infinity. We start by considering
the term ∆n , for which we will need that
Z π
1
hΨcn , Ψcn i =
Pn Γn Γ∗n Pn⊤ dω = Pn Pn⊤ .
(41)
2π −π
||∆n || ≤ 2||hΩn , Ωn i|| ||hΩn , Ωn i−1 || ||Pn || ||Ψcn − Ψ̂cn ||H2
+ ||hΩn , Ωn i|| ||hΩn , Ωn i−1 || ||Ψcn − Ψ̂cn ||2H2

+ 2||Pn || ||Ψcn − Ψ̂cn ||H2 + ||Ψcn − Ψ̂cn ||2H2 . (42)

For row i of Ψcn − Ψ̂cn , we have that
P∞ (i) −k Pn
(i)
c,i
− k=1 βk F3 q −k ||
||Ψc,i
n − Ψ̂n || = ||F3
k=1 βk q
P∞
(i)
≤ ||F3 || || k=n+1 βk q −k || ≤ Cλn

Then,

√
||Ψcn − Ψ̂cn ||H2 ≤ C nλn → 0 as n → ∞.

Moreover, using (29), we have
Z π
1
||Pn || =
Γn Γ∗n F4∗ dω ≤ ||F4∗ ||H∞ = C,
2π −π

Λn =hΨ̂cn , Ωn ihΩn , Ωn i−1 hΩn , Ψ̂cn i

Λn = hΨ̂cn , Ψ̂cn i.

(40)

Using also (29) and (30), we can write

hΩn , Ψan i

−1

(39)

(37)

and

c
c
c
c
c
c
c
c
c
c
∆(2)
n = hΨn −Ψ̂n , Ψn i+hΨn , Ψn −Ψ̂n i+hΨn −Ψ̂n , Ψn −Ψ̂n i.

(43)

(44)

||hΩn , Ωn i||

Z π 
1
Γn Γ∗n (|F1 |2 + |F2 |2 ) Γn Γ∗n F1 F3∗
dω
=
Γn Γ∗n F3 F1∗
Γn Γ∗n |F3 |2
2π −π
≤ || |F1 |2 + |F2 |2 ||H∞ + 2||F1 F3 ||H∞ + || |F3 | ||H∞
= C,

and (by assumption) ||hΩn , Ωn i−1 || < C. Together with (42)
and (43),
||∆n || → 0 as n → ∞.
(45)
Then, if the remaining matrix norms in the right-hand side
of (40) are bounded, this term will tend to zero as n → ∞.
For hγ, Ψn i, we have

Z π 
1
Γn Γ∗n F4 F5
dω
||hγ, Ψn i|| =
2π −π Γn Γ∗n F4 F6
(46)
≤ ||F4 F5 ||H∞ + ||F4 F6 ||H∞ = C.

IEEE TRANSACTIONS ON SIGNAL PROCESSING

9

Also, the inverse of hΨcn , Ψcn i exists because Pn is lower
triangular with non-zero diagonal entries, and
||hΨcn , Ψcn i−1 || = Pn−⊤ Pn−1
Z π
Z π
1
1
∗ −1
=
Γn Γn F4 dω
Γn Γ∗n F4−∗ dω
2π −π
2π −π
≤ ||F4 ||2H∞ = C
(47)
Using (39), (40), (45), (46), and (47), we have
−1

lim hγ, Ψn i[hΨn , Ωn ihΩn , Ωn i

n→∞

hΩn , Ψn i]

−1

hΨn , γi

= lim hγ, Ψn ihΨcn , Ψcn i−1 hΨn , γi. (48)
n→∞

Outer projection: hγ, Ψn ihΨcn , Ψcn i−1 hΨn , γi: Recalling
that Ψn can be written as (32), we use that γ is causal and
Ψan anti-causal to conclude that, analogously to (33),
lim hγ, Ψn ihΨcn , Ψcn i−1 hΨn , γi

n→∞

=

lim hγ, Ψcn ihΨcn , Ψcn i−1 hΨcn , γi.
n→∞

(49)

Now, row i of γ can be written as γ(i)
=:
P∞ (i) −k
(i)
k
s
q
0
],
where
|s
|
≤
Cλ
,
λ
<
1
due
to
[
k=1 k
k
P
(i)
exponential stability. Let also γ n (i) := [ nk=1 sk q −k 0 ]
be row i of a matrix γ n . We re-write the right side of (49) as
lim hγ, Ψcn ihΨcn , Ψcn i−1 hΨcn , γi

n→∞

= lim hγ
n→∞

n

, Ψcn ihΨcn , Ψcn i−1 hΨcn , γ n i

+

lim ∆(3) ,
n→∞ n

(50)

where
n
c
c
c −1
∆(3)
hΨcn , γi
n = hγ − γ , Ψn ihΨn , Ψn i

+ hγ, Ψcn ihΨcn , Ψcn i−1 hΨcn , γ − γ n i

+ hγ − γ n , Ψcn ihΨcn , Ψcn i−1 hΨcn , γ − γ n i.
(3)

Using a similar approach for ∆n as we did for ∆n , it can
be shown that
||∆(3)
n ||

n

≤ C(||γ − γ ||H2 + ||γ −

γ n ||2H2 )

→ 0 as n → ∞.
(51)

Thus, (50) reduces to

= lim hγ n , Ψcn ihΨcn , Ψcn i−1 hΨcn , γ n i. (52)
n→∞

Finally, we have that
hγ

, Ψcn ihΨcn , Ψcn i−1 hΨcn , γ n i

n

n

= hProjSΨc γ , ProjSΨc γ i
n

n

= hγ n , γ n i,

(53)
where the last equality follows from γ n ∈ SΨcn , as consequence of Pn being invertible. Thus, replacing (53) in (52),
we have
lim hγ, Ψcn ihΨcn , Ψcn i−1 hΨcn , γi = lim hγ n , γ n i = hγ, γi,

n→∞

WLS
||θ̂N
− θo || → 0, as N → ∞, w.p.1.

(54)

To do this, we use (19) to write
LS
WLS
) kQn (η̂N )k
θ̂N
− θo ≤ M −1 (η̂N , θ̂N
LS
· Wn (θ̂N
) kTn (θo )k η̂N − ηon(N ) . (55)

From [24], we have that
||η̂N − ηon(N ) || → 0 as N → ∞, w.p.1,
||Tn (θo )|| ≤ C

∀n,

(56)
(57)

and that, w.p.1, there exists N̄ such that
||Qn (η̂N )|| ≤ C

∀N > N̄ .

(58)

LS
Wn (θ̂N
)

Then, we have left to show that
is bounded and
LS
M (η̂N , θ̂N
) is invertible for sufficiently large N .
LS
We begin by considering the inverse of Wn (θ̂N
), for which
we have
LS
n −1 ⊤ LS
||Tn (θ̂N
)[RN
] Tn (θ̂N )|| ≤ ||Tn (θo )[R̄n ]−1 Tn⊤ (θo )||
LS
n −1 ⊤ LS
+ ||Tn (θ̂N
)[RN
] Tn (θ̂N ) − Tn (θo )[R̄n ]−1 Tn⊤ (θo )||.
(59)
In turn, it can be shown that

Tn (θo )[R̄n ]−1 Tn⊤ (θo ) = hΨn , Ωn ihΩn , Ωn i−1 hΩn , Ψn i,
where Ψn and Ωn are given by (22), which satisfy the
conditions of Theorem 2 with F1 (q) = −Go (q)So (q),
F2 (q) = −σo Ho (q), F3 (q) = So (q)Fr (q), and F4 (q) =
Fo (q)/[So (q)Fr (q)] From (37), (41), (44), and (45), we have
that there is n̄ such that
||Tn (θo )[R̄n ]−1 Tn⊤ (θo )|| = ||hΨcn , Ψcn i + ∆n || ≤ C ∀n > n̄.
(60)
Then, and using also (47), we have
||W̄ (θo )|| := ||[Tn (θo )[R̄n ]−1 Tn⊤ (θo )]−1 ||

= ||(hΨcn , Ψcn i + ∆n )−1 || ≤ C

∀n > n̄. (61)

Concerning the second term on the right-hand side of (59), we
can write
LS
n −1 ⊤ LS
||Tn (θ̂N
)[RN
] Tn (θ̂N ) − Tn (θo )[R̄n ]−1 Tn⊤ (θo )||

lim hγ, Ψcn ihΨcn , Ψcn i−1 hΨcn , γi

n→∞

n

A PPENDIX C
P ROOF OF T HEOREM 3
We will show that

n→∞

which, together with (48), (49), (50), and (51) implies (23), as
we wanted to show.

LS
n −1
LS
≤ ||Tn (θ̂N
) − Tn (θo )|| ||[RN
] || ||Tn (θ̂N
)||

LS
n −1
+ ||Tn (θ̂N
) − Tn (θo )|| ||[RN
] || ||Tn (θo )||

(62)

n −1
+ ||Tn (θo )||2 ||[RN
] − [R̄n ]−1 ||

Now, the results in [24] apply to (62). In particular, there is
N̄ such that
n −1
n
||[RN
] || ≤ C ∀N > N̄ ,
||[RN
]|| ≤ C ∀N > N̄ ,
LS
||Tn (θ̂N
)|| ≤ C

∀N > N̄

(63)

w.p.1, and that
LS
||Tn (θ̂N
) − Tn (θo )|| → 0 as N → ∞, w.p.1,

n −1
n
n
n −1
||[RN
] − [R̄n ]−1 || ≤ ||R̄N
|| ||R̄n − RN
|| ||[RN
] ||
→ 0 as N → ∞, w.p.1.

IEEE TRANSACTIONS ON SIGNAL PROCESSING

10

Together with (57), we conclude that
LS
n −1 ⊤ LS
||Tn (θ̂N
)[RN
] Tn (θ̂N ) − Tn (θo )[R̄n ]−1 Tn⊤ (θo )||
→ 0 as N → ∞, w.p.1. (64)

Using (64), (60), and (59), there is N̄ such that
LS
n −1 ⊤ LS
||Tn (θ̂N
)[RN
] Tn (θ̂N )|| ≤ C

∀N > N̄
n −1

w.p.1.

LS
||W (θ̂N
) − W̄ (θo )|| → 0, as N → ∞, w.p.1.

by
=

∀N > N̄ .

LS
LS
M (η̂N , θ̂N
) = Q⊤
n (η̂N )Wn (θ̂N )Qn (η̂N ),
√
LS
LS
n(N )
x(θo ; η̂N , θ̂N
) = N Q⊤
).
n (η̂N )Wn (θ̂N )Tn (θo )(η̂N − ηo

If we assume that
(65)

Moreover, (65) and (61) imply that, w.p.1,
LS
||W (θ̂N
)|| ≤ C

In this appendix, we prove Theorem 4. We begin by
reformulating (19) as
√
WLS
LS
LS
N (θ̂N
− θo ) = M −1 (η̂N , θ̂N
)x(θo ; η̂N , θ̂N
),
where

Tn⊤ (θo ),
LS
W (θ̂N
)

Because of (64) and invertibility of Tn (θo )[R̄ ]
continuity of eigenvalues there is N̄ such that
LS
n −1 ⊤ LS −1
[Tn (θ̂N
)[RN
] Tn (θ̂N )] exists for all N > N̄ , and

A PPENDIX D
P ROOF OF T HEOREM 4

(66)

LS
x(θo ; η̂N , θ̂N
) ∼ AsN (0, X),

we have that, using (71) and [29, Lemma B.4],
√

WLS
N (θ̂N
− θo ) ∼ AsN 0, M −1 XM −1 .

LS
Having shown (66), we have left to show that M (η̂N , θ̂N
)
is invertible for sufficiently large N , in order to show (54). We will proceed to show that (72) is verified with
Recall the definition (20), which can alternatively be written
X = σo2 lim Q⊤
(η n )W̄n (θo )Qn (ηon ) = σo2 M,
as (21), where γ is given by (22). Then, from Theorem 2 with
n→∞ n o

Bo (q)So (q)Fr (q)
,
Fo (q)

(72)

(73)

Ao (q)So (q)Fr (q)
,
Fo (q)

where the second equality follows directly from (20) and (67).
We now proceed to show the first equality.
The idea is, as in [24], to use [27, Theorem 7.3]. With
we have that
LS
this purpose, we first show that x(θo ; η̂N√, θ̂N
) has the same
M̄ (ηo , θo ) = M.
(67)
asymptotic distribution and covariance as N Υn (η̂N − η̄ n(N ) ),
n
where M is given by (8). Because the inverse of M corre- where Υ is a deterministic matrix. From [24], it follows
LS
sponds to the CR bound for an open-loop problem with input that x(θo ; η̂N , θ̂N ) has the same asymptotic covariance and
ut = So (q)rt , and the CR bound exists for an informative distribution as
√
experiment [1], we conclude that M̄ (ηo , θo ) is invertible. Then,
n(N )
LS
N Q⊤
)Wn (θ̂N
)Tn (θo )(η̂N − η̄ n(N ) ).
(74)
n (ηo
we analyze the difference
LS
The next step is to show that Wn (θ̂N
) can be replaced by
LS
n
n
||M (η̂N , θ̂N
) − Q⊤
n (ηo )W̄n (θo )Qn (ηo )||
W̄n (θo ) in (74) without affecting the asymptotic distribution
n
LS
≤ ||Qn (η̂N ) − Qn (ηN
)|| ||W (θ̂N
)|| ||Qn (η̂N )||
and covariance. However, this step does not follow directly
n
LS
n
from [24], as inverses of the matrices that compose Wn cannot
+ ||Qn (η̂N ) − Qn (ηN )|| ||W (θ̂N )|| ||Qn (ηo )||
be taken individually here, because Tn is not square. In the
n
LS
)||2 ||W (θ̂N
) − W̄ (θo )||
+ ||Qn (ηN
following, we do this for a non-square Tn matrix.
From [24], we have that
First, re-write (74) as
n
√
||Qn (η̂N ) − Qn (ηN
)|| → 0, as N → ∞, w.p.1,
n(N )
N Q⊤
)W̄n (θo )Tn (θo )(η̂N − η̄ n(N ) )
n (ηo
n
√
||Qn (ηo )|| ≤ C ∀n.
(68)
n(N )
LS
+ N Q⊤
)[Wn (θ̂N
) − W̄n (θo )]Tn (θo )(η̂N − η̄ n(N ) ).
n (ηo
(75)
Together with (65) and (58), we conclude that
F5 (q) = −

F6 (q) =

Then, it follows from [29, Lemma B.4] that (74) and the
first term in (75) have the same asymptotic distribution and
→ 0, as N → ∞, w.p.1. (69) covariance if the second term in (75) tends to 0 as N → ∞
Using (69), invertibility of M̄ (ηo , θo ), and continuity of eigen- w.p.1. We proceed to show this. Consider
√
LS
) is
values, we have that there is N̄ such that M (η̂N , θ̂N
LS
(ηon(N ) )[Wn (θ̂N
) − W̄n (θo )]Tn (θo )(η̂N − η̄ n(N ) )||
|| N Q⊤
√ n ⊤
invertible for all N > N̄ ,
LS
≤ N ||Qn (ηon(N ) )|| ||Wn (θ̂N
) − W̄n (θo )||
−1
LS
||M (η̂N , θ̂N )|| ≤ C ∀N > N̄ ,
(70)
· ||Tn (θo )|| ||η̂N − η̄ n(N ) ||
√
LS
) − W̄n (θo )|| ||η̂N − η̄ n(N ) ||,
≤ C N ||Wn (θ̂N
and, using also (67) and (20),
(76)
LS
M −1 (η̂N , θ̂N
) → M −1 as N → ∞, w.p.1.
(71) where the last inequality follows from (57) and (68). Writing
LS
n
n
||M (η̂N , θ̂N
) − Q⊤
n (ηo )W̄n (θo )Qn (ηo )||

Finally, using (70), (66), (58), (57), (56), and (55), we
conclude that (54) is satisfied, as we wanted to show.

LS
LS
LS
Wn (θ̂N
)−W̄n (θo ) = W̄n (θo )[W̄n−1 (θo )−Wn−1 (θ̂N
)]Wn (θ̂N
),

IEEE TRANSACTIONS ON SIGNAL PROCESSING

11

LS
and because (66) and (61) guarantee that W̄n (θo ) and Wn (θ̂N
)
are bounded (in the latter, for sufficiently large N ), it follows
from (76) that
√
n(N )
LS
|| N Q⊤
)[W
) − W̄n (θo )]Tn (θo )(η̂N − η̄ n(N ) )||
n (ηo
√ n (θ̂N−1
LS
≤ C N ||Wn (θ̂N
) − W̄n−1 (θo )|| ||η̂N − η̄ n(N ) ||.

LS
Now, the term ||Wn−1 (θ̂N
) − W̄n−1 (θo )|| corresponds to (62);
so, using [24]
!
r
log N
LS
2
(1 + d(N )) ,
||Tn (θ̂N ) − Tn (θo )|| = O
n (N )
N
(77)
(57) and (63), the first two terms on the right-hand side of (62)
decay with (77). For the third term, we first write
n −1
n
n
n −1
||[RN
] − [R̄n ]−1 || ≤ ||R̄N
|| ||R̄n − RN
|| ||[RN
] ||. (78)

In [27], it is shown that
r

log N
n2 (N )
n
||R̄n − RN
|| = O 2 n2 (N )
+C
N
N

!

.

(79)

Then, using also (78), (63), and (57), we have that the third
term on the right-hand side of (62) decays according to (79).
Thus, we have that
!
r
log
N
LS
,
(80)
||Wn (θ̂N ) − W̄n (θo )|| = O 2 n2 (N )
N
considering only the slowest-decaying term. Then, from (76),
(80), and [27]
!
r
n(N ) log N
n(N )
[1 + d(N )] ,
||η̂N − η̄
|| = O
N
it follows that
√
LS
) − W̄n (θo )|| ||η̂N − η̄ n(N ) || → 0
C N ||Wn (θ̂N

as N → ∞, w.p.1.

Finally, using (76) and (75), it follows from [29, Lemma B.4]
that
√
n(N )
N Q⊤
)W̄n (θo )Tn (θo )(η̂N − η̄ n(N ) )
(81)
n (ηo
LS
and (74)—and, in turn, x(θo ; η̂N , θ̂N
)—have the same asymptotic distribution and covariance. Thus, we will analyze (81)
LS
instead of x(θo ; η̂N , θ̂N
).
Applying [27, Theorem 7.3] to (81) with Υn =
n(N )
Q⊤
)W̄n (θo )Tn (θo )—and recalling that it has the same
n (ηo
LS
asymptotic distribution and covariance as x(θo ; η̂N , θ̂N
)—we
LS
have that x(θo ; η̂N , θ̂N ) is distributed according to (72) with
n(N )
X = lim Q⊤
)W̄n (θo )Tn (θo )σo2
n (ηo
n→∞

n ⊤
· R̄N
Tn (θo )W̄n (θo )Qn (ηon(N ) )

(82)

n(N )
= σo2 lim Q⊤
)W̄n (θo )Qn (ηon(N ) ) = σo2 M,
n (ηo
n→∞

where the last equality follows from (20). Then, replacing (82)
in (73), we obtain
√

WLS
N (θ̂N
− θo ) ∼ AsN 0, M −1 .

R EFERENCES
[1] L. Ljung. System Identification. Theory for the User. Prentice-Hall,
1999.
[2] T. Söderström and P. Stoica. Instrumental Variable Methods for System
Identificaiton. Springer, 1983.
[3] M. Gilson and P. van den Hof. Instrumental variable methods for closedloop system identification. Automatica, 41(2):241–249, 2005.
[4] M. Gilson, H. Garnier, P. Young, and P. van den Hof. Refined
instrumental variable methods for closed-loop system identification. In
15th IFAC Symp. on System Identification, Saint-Malo, France, 2009.
[5] M. Verhaegen and P. DeWilde. Subspace model identification, part I:
The output-error state-space model identification class of algorithms.
International Journal of Control, 56:1187–1210, 1992.
[6] P. van Overschee and B. de Moor. N4SID: Subspace algorithms for
identification of combined deterministic-stochastic systems. Automatica,
30:75–93, 1994.
[7] M. Verhaegen. Subspace model identification, part I: The outputerror state-space model identification class of algorithms. Automatica,
29(4):1027–1040, 1993.
[8] M. Jansson. Subspace identification and ARX modeling. In IFAC
Symposium on System Identification, Rotterdam, The Netherlands, 2003.
[9] J. Schoukens, Y. Rolain, G. Vandersteen, and R. Pintelon. User friendly
Box-Jenkins identification using nonparametric noise models. In 50th
IEEE Conference on Decision and Control, Orlando, Florida, USA,
2011.
[10] R. Pintelon, J. Schoukens, G. Vandersteen, and K. Barb. Estimation of
nonparametric noise and FRF models for multivariable systems—Part I:
theory. Mech. Systems and Signal Processing, 24(3):573–595, 2010.
[11] R. Pintelon, J. Schoukens, G. Vandersteen, and K. Barb. Estimation
of nonparametric noise and FRF models for multivariable systems—
Part II: extensions, applications. Mech. Systems and Signal Processing,
24(3):596–616, 2010.
[12] Y. Zhu. Multivariable System Identification for Process Control. Pergamon, 2001.
[13] Y. Zhu and H. Hjalmarsson. The Box-Jenkins Steiglitz-McBride algorithm. Automatica, 65:170–182, 2016.
[14] N. Everitt, M. Galrinho, and H. Hjalmarsson. Open-loop asymptotically
efficient model reduction with the Steiglitz-McBride method. under
review at Automatica (arXiv:1610.08534).
[15] M. Galrinho, C. R. Rojas, and H. Hjalmarsson. A weighted least-squares
method for parameter estimation in structured models. In 53rd IEEE
Conference on Decision and Control, pages 3322–3327, Los Angeles,
USA, 2014.
[16] C. Sanathanan and J. Koerner. Transfer function synthesis as a ratio
of two complex polynomials. IEEE Transactions on Automatic Control,
8(1):56–58, 1963.
[17] A. G. Evans and R. Fischl. Optimal least squares time-domain synthesis
of recursive digital filters. IEEE Transactions on Audio and Electroacoustics, 21(1):61–65, 1973.
[18] Y. Bresler and A. Macovski. Exact maximum likelihood parameter estimation of superimposed exponential signals in noise. IEEE Transactions
on Signal Processing, 34(5):1081–1089, 1986.
[19] A. K. Shaw. Optimal identification of discrete-time systems from
impulse response data. IEEE Trans. on Signal Processing, 42(1):113–
120, 1994.
[20] A. K. Shaw, P. Misra, and R. Kumaresan. Identification of a class of
multivariable systems from impulse response data: Theory and computational algorithm. Circuits, Systems and Signal Processing, 13(6):759–
782, 1994.
[21] P. Lemmerling, L. Vanhamme, S. van Huffel, and B. de Moor. IQMLlike algorithms for solving structured total least squares problems: a
unified view. Signal Processing, 81:1935–1945, 2001.
[22] K. Steiglitz and L. E. McBride. A technique for the identification of
linear systems. IEEE Trans. on Automatic Control, 10:461–464, 1965.
[23] J. H. McClellan and D. Lee. Exact equivalence of the SteiglitzMcBride iteration and IQML. IEEE Transactions on Signal Processing,
39(2):509–5012, 1991.
[24] M. Galrinho, C. R. Rojas, and H. Hjalmarsson. Parametric identification
using weighted null-space fitting. submitted to IEEE Transactions on
Automatic Control (arXiv 1708.03946).
[25] U. Forssell and L. Ljung. Closed-loop identification revisited. Automatica, 35:1215–1241, 1999.
[26] H. Hjalmarsson and J. Mårtensson. A geometric approach to variance
analysis in system identification. IEEE Transactions on Automatic
Control, 56(5):983–997, 2011.

IEEE TRANSACTIONS ON SIGNAL PROCESSING

[27] L. Ljung and B. Wahlberg. Asymptotic properties of the least-squares
method for estimating transfer functions and disturbance spectra. Advances in Applied Probabilities, 24:412–440, 1992.
[28] V. Peller. Hankel operators and their applications. Springer, 2003.
[29] T. Söderström and P. Stoica. System Identification. Prentice Hall, 1989.

PLACE
PHOTO
HERE

Miguel Galrinho was born in 1988. He received
his M.S. degree in aerospace engineering in 2013
from Delft University of Technology, The Netherlands, and the Licentiate degree in automatic control
in 2016 from KTH Royal Institute of Technology,
Stockholm, Sweden.
He is currently a PhD student at KTH, with the
Department of Automatic Control, School of Electrical Engineering, under supervision of Professor
Håkan Hjalmarsson. His research is on least-squares
methods for identification of structured models.

PLACE
PHOTO
HERE

Cristian R. Rojas (M’13) was born in 1980. He
received the M.S. degree in electronics engineering
from the Universidad Técnica Federico Santa Marı́a,
Valparaı́so, Chile, in 2004, and the Ph.D. degree in
electrical engineering at The University of Newcastle, NSW, Australia, in 2008.
Since October 2008, he has been with the Royal
Institute of Technology, Stockholm, Sweden, where
he is currently Associate Professor at the Department
of Automatic Control, School of Electrical Engineering. His research interests lie in system identification

and signal processing.

Håkan Hjalmarsson (M’98–SM’11–F’13) was
born in 1962. He received the M.S. degree in electrical engineering in 1988, and the Licentiate and
Ph.D. degrees in automatic control in 1990 and
PLACE
1993, respectively, all from Linköping University,
PHOTO
Linköping, Sweden.
HERE
He has held visiting research positions at California Institute of Technology, Louvain University,
and at the University of Newcastle, Australia. His research interests include system identification, signal
processing, control and estimation in communication
networks and automated tuning of controllers.
Dr. Hjalmarsson has served as an Associate Editor for Automatica (1996–
2001) and for the IEEE Transactions on Automatic Control (2005–2007),
and has been Guest Editor for the European Journal of Control and Control
Engineering Practice. He is a Professor at the School of Electrical Engineering,
KTH, Stockholm, Sweden. He is a Chair of the IFAC Coordinating Committee
CC1 Systems and Signals. In 2001, he received the KTH award for outstanding
contribution to undergraduate education. He is co-recipient of the European
Research Council advanced grant.

12

