RECURSIVE EXTENDED COMPENSATED LEAST SQUARES

BASED ALGORITHM FOR

ERRORS-IN-VARIABLES IDENTIFICATION

Tomasz Larkowski, Jens G. Linden and Keith J. Burnham

Control Theory and Applications Centre, Coventry University, Priory Street, Coventry, CV1 5FB, U.K.

Keywords:

Bias compensation, Errors-in-variables, Parameter estimation, Recursive algorithms, System identiﬁcation.

Abstract:

An algorithm for the recursive identiﬁcation of single-input single-output linear discrete-time time-invariant

errors-in-variables system models in the case of white input and coloured output noise is presented. The

approach is based on a bilinear parametrisation technique which allows the model parameters to be estimated

together with the auto-correlation elements of the input/output noise sequences. In order to compensate for the

bias in the recursively obtained least squares estimates, the extended bias compensated least squares method

is used. An alternative for the online update of the associated pseudo-inverse of the extended observation

covariance matrix is investigated, namely an approach based on the matrix pseudo-inverse lemma and an

approach based on the recursive extended instrumental variables technique. A Monte-Carlo simulation study

demonstrates the appropriateness and the robustness against noise of the proposed scheme.

1 INTRODUCTION

The errors-in-variables (EIV) approach forms an ex-

tension of the standard output error system setup

in which it is postulated that only the output mea-

surements are uncertain. In the EIV framework all

measured signals, hence, including the system in-

put, are assumed to be contaminated with noise, see

(S¨oderstr¨om, 2007) for the recent survey on this sub-

ject. The EIV framework can offer advantages over

the classical approach, mainly when the description

of the internal laws governing a system is of prime in-

terest, e.g. application areas in chemistry, image pro-

cessing, fault detection etc., see (S¨oderstr¨om, 2007;

Markovsky and Van Huffel, 2007) for further details.

One of the EIV techniques that has been shown

to be robust and to yield relatively precise estimates

is the extended compensated least squares (ECLS)

method. The approach is based on the extended bias

compensated least squares (EBCLS) and utilises sep-

arable nonlinear least squares to solve the resulting

overall identiﬁcation problem. The method was ﬁrst

proposed in (Ekman, 2005a), which considered the

case of white input and output noise sequences and

subsequently extended to handle the case of coloured

output noise in (Ekman et al., 2006). Further anal-

ysis, considering a generalised framework, has been

carried out in (Mahata, 2007).

Alternatively, by exploiting the property that the

overall optimisation problem is bilinear in the un-

knowns, see (Ljung, 1999), which in this case corre-

sponds to the model parameters and the input/output

noise auto-correlation elements, the principle of bi-

linear parametrisation can be utilised. The resulting

scheme, termed here the extended bilinear parametri-

sation method (EBPM) involves solving iteratively

two ordinary least squares problems, see (Larkowski

et al., 2008) for details. Although the quality of the

parameters obtained by the EBPM is comparable to

the quality of the estimates yielded by the ECLS, an

important distinction is that the EBPM is signiﬁcantly

less computationally demanding than the ECLS tech-

nique.

The bilinear parametrisation method was ﬁrst

utilised to solve the EIV identiﬁcation problem in

a recursive manner in (Ekman, 2005b) for the case

of white input and output noise. It has also been

exploited in (Ikenoue et al., 2008) for the case of

coloured input and output noise sequences and for the

purpose of ofﬂine as well as online estimation. How-

ever, in both cases the term ‘bilinear parametrisation’

has not been explicitly stated. In (Ekman, 2005b) the

constructed recursive algorithm is not computation-

ally attractive, since its complexity at each iteration is

actually greater than that of the corresponding batch

algorithm applied in an ofﬂine manner at each recur-

142

Larkowski T., G. Linden J. and J. Burnham K. (2009).

RECURSIVE EXTENDED COMPENSATED LEAST SQUARES BASED ALGORITHM FOR ERRORS-IN-VARIABLES IDENTIFICATION.

In Proceedings of the 6th International Conference on Informatics in Control, Automation and Robotics - Signal Processing, Systems Modeling and

Control, pages 142-147

DOI: 10.5220/0002211901420147

 SciTePress

sion. Whereas, in (Ikenoue et al., 2008) due to a spe-

cial choice of the instruments, the resulting algorithm

is not causal, in general, hence its recursive imple-

mentation yields delayed estimates.

In this paper a recursive realisation of the EBPM

is presented for a discrete-time linear time-invariant

(LTI) single-input single-output (SISO) system model

in the case of the white input and coloured output

noise and it is demonstrated that the above men-

tioned shortcomings may be avoided. The bias of

the recursively calculated least squares (LS) estima-

tor is removed at each recursion via the extended bias

compensated least squares (EBCLS) technique. The

online update of the pseudo-inverse of the overde-

termined observation matrix is realised by consider-

ing an alternative, namely an approach based on the

pseudo-inverse lemma, see (Feng et al., 2001) and an

approach based on the recursive extended instrumen-

tal variables technique, see (Friedlander, 1984). The

two resulting algorithms are analysed and compared

with their ofﬂine counterpart via a Monte-Carlo sim-

ulation study. It is shown that the instrumental vari-

ables approach is the more preferable due to its supe-

rior robustness and improved convergence properties,

in general.

2 NOTATION AND PROBLEM

STATEMENT

Consider a discrete-time LTI SISO system repre-

sented by the difference equation

A(q

−1

= B(q

−1

, (1)

where the polynomials A(q

−1

) and B(q

−1

) are given

A(q

−1

) , 1+ a

−1

+ . . . + a

−n

, (2a)

B(q

−1

) , b

−1

+ . . . + b

−n

(2b)

with q

−1

being the backward shift operator, deﬁned

by q

−1

, x

k−1

. The unknown noise-free input and

noise-free output signals denoted u

and y

, respec-

tively, are related to the available noisy variables, de-

noted u

and y

, such that

= u

+ ˜u

, y

= y

+ ˜y

, (3)

where ˜u

and ˜y

denote the input and output measure-

ment noise sequences, respectively. The following

standard assumptions, see e.g. (Ekman et al., 2006),

are introduced:

A1 The LTI system (1) is asymptotically stable, i.e.

A(q

−1

) has all zeros inside the unit circle.

A2 All system modes are observable and control-

lable, i.e. A(q

−1

) and B(q

−1

) share no common

factors.

A3 The system structure, i.e. n

and n

, is known a

priori and n

≥ n

A4 The true input u

is a zero mean, ergodic random

sequence persistently exciting and of sufﬁciently

high order, i.e. at least of order n

+ n

A5a The additive input noise sequence ˜u

of un-

known variance σ

˜u

is an ergodic zero mean white

process.

A5b The additive output noise sequence ˜y

an ergodic zero mean process characterised

by an unknown auto-covariance sequence

˜y

(0), r

˜y

(1), . . .}.

A6 The input/output noise sequences are mutually

uncorrelatedand uncorrelated with signals u

and

By postulating that the output noise sequence exhibits

an arbitrary degree of correlation allows for measure-

ment sensor uncertainties to be taken into account, as

well as potential disturbances in the process.

The system parameter vector is denoted

θ ,





∈ R

, (4a)

a ,



. . . a



∈ R

, (4b)

b ,



. . . b



∈ R

, (4c)

where n

= n

+ n

. The extended regressor vectors

for the k-th measured data are deﬁned as



−y



∈ R

, (5a)



−y



∈ R

, (5b)

where





∈ R

, (5c)



−y

k−1

. . . − y

k−n



∈ R

, (5d)



k−1

. . . u

k−n



∈ R

. (5e)

The noise contributions in the corresponding regres-

sor vectors are denoted by a tilde, i.e.

[·], whereas the

noise-free signals are denoted by a zero subscript, i.e.

[·]

. From (3) it follows that

. (6)

The notation Σ

is used as a general notion for the

covariance matrix of the vectors g

and d

, whereas

is utilised for a covariance vector with f

being

a scalar. The corresponding estimates are denoted by

a hat. In addition, 0

g×d

denotes the null matrix of

arbitrary dimension g × d and a single index is used

RECURSIVE EXTENDED COMPENSATED LEAST SQUARES BASED ALGORITHM FOR

ERRORS-IN-VARIABLES IDENTIFICATION

143

in the case of a column vector as well as in the case of

a square matrix, e.g. the identity matrix I

. The auto-

correlation elements, denoted r

˜y

(·) are deﬁned as

˜y

(τ) , E [ ˜y

˜y

k−τ

], (7)

where E[·] is the expected value operator. Introducing

ρ ,



˜u



∈ R

, (8a)



˜y

(0) . . . r

˜y

)



∈ R

, (8b)

the dynamic identiﬁcation problem in the EIV frame-

work considered here is formulated as:

Problem 1. (Dynamic EIV identiﬁcation problem)

Given N samples of the measured signals, i.e. {u

}

k=1

and {y

}

k=1

, determine the vector

Θ ,





∈ R

. (9)

3 REVIEW OF APPROACHES

This section brieﬂy reviewsthe EBCLS technique and

the ofﬂine EBPM algorithm.

3.1 Extended Bias Compensated Least

Squares

Denoting an estimate by

[·], a solution of the system

(1)-(3) in the LS sense is given by

†

xϕ

, (10)

where [·]

†

is the pseudo inverse operator deﬁned by

†

, (A

−1

, x

∈ R

denotes an arbitrary in-

strumental vector with n

≥ n

. Due to the mea-

surement noise, unless the elements of x

are uncor-

related with

, the solution obtained is biased. In

order to achieve an unbiased estimate of θ, a bias

compensation procedure is required to be carried out

(S¨oderstr¨om, 2007). This consideration yields the

EBCLS estimator deﬁned as

EBCLS



xϕ

− Σ

˜x



†



− ξ

˜x˜y



. (11)

Note that Σ

˜x

and ξ

˜x˜y

, in general, are functions of ρ,

which, in turn, will depend on the elements contained

in the instrument vector x

3.2 Extended Bilinear Parametrisation

Method

The bilinear parametrisation method is applicable

for problems that are bilinear in the parameters, see

(Ljung, 1999) for details, and it is presented here

in accordance with the development proposed in

(Larkowski et al., 2008).

Based on the EBCLS rule given by (11) a bilinear

(in the parameters) cost function can be formulated,

i.e.

Θ = argmin

V(Θ), (12)

where

V(Θ) ,



− ξ

˜x˜y

−



xϕ

− Σ

˜x





. (13)

Note that the instruments x

must be chosen such that

the resulting problem is soluble, i.e. the total number

of unknowns is less than or equal to the total num-

ber of equations, see (Larkowski et al., 2008) for a

detailed treatment. Alternatively, utilising the bilin-

earity property, (13) can be re-expressed as

V(Θ) =



−

xϕ

θ−Wρ



, (14)

where W , S

− S

(θ) ∈ R

×(n

+2)

such that S

ρ ,

˜x˜y

and S

(θ)ρ , Σ

˜x

θ.

It is observed that for ﬁxed ρ (i.e. the expres-

sions Σ

˜x

and ξ

˜x ˜y

) the cost function (13) is linear in

θ. Analogously, for ﬁxed θ (i.e. the matrix W) the

cost function (14) is linear in ρ. Consequently, a nat-

ural approach is to treat (13) and (14) as separate LS

problems, cf. (Ljung, 1999). This leads to a two-step

algorithm where the LS solutions of the sub-problems

deﬁned by (13) and (14) are obtained at each iteration.

Furthermore, local convergence of such algorithm is

guaranteed, see (Ljung, 1999).

4 RECURSIVE EXTENDED

BILINEAR

PARAMETRISATION METHOD

This section presents the proposed recursive realisa-

tion of the EBPM technique, denoted REBPM. First

the problem of an online update of the parameter vec-

tor is addressed. Subsequently, two approaches for

updating the pseudo-inverse of the extended obser-

vation matrix are considered. Finally, the problem

of calculating the input noise variance and the auto-

correlation elements of the output noise is discussed.

4.1 Recursive Update of Parameter

Vector

Considering (11) and by making use of (10) it follows

that

θ =

†

xϕ



˜x

θ− ξ

˜x˜y



. (15)

ICINCO 2009 - 6th International Conference on Informatics in Control, Automation and Robotics

144

It is remarked that the expression

†

xϕ



˜x

θ− ξ

˜x˜y



represents the bias of the LS estimator. Since the true

value of θ on the right hand side of (15) is unknown, a

natural approach is to utilise the most recent estimate,

i.e. the previous value. This leads to the following

recursive EBCLS scheme

EBCLS



xϕ



†



˜x

k−1

EBCLS

− ξ

˜x˜y



. (16)

Despite that an inevitable error is introduced by as-

suming

EBCLS

≈

k−1

EBCLS

, the above approach also

known as the stationary iterativeLS principle (Bj¨orck,

1996), has been successfully employed in several re-

cursive as well as iterative algorithms.

4.2 Recursive Update of Pseudo-inverse

Considering equation (16), it is observed that a

recursive update of the pseudo-inverse of

xϕ

as well

as of the LS estimate, i.e.

, is required. This

problem can be tackled by two approaches described

below.

Approach based on the Matrix Pseudo-inverse

Lemma - REBPM

. The ﬁrst, i.e. direct approach

is to utilise an extension of the matrix inverse lemma,

namely the matrix pseudo-inverse lemma, see (Feng

et al., 2001). This allows the recursive computation of

the expression

†

xϕ

as well as the corresponding

The algorithm can be summarised as follows:

k−1

+ L



− ϕ

k−1



, (17a)



k−1

xϕ



†

k− 1+ ϕ



k−1

xϕ



†

, (17b)



xϕ



†

k− 1





k−1

xϕ



†

− L



k−1

xϕ



†



, (17c)

xϕ

k−1

xϕ



−

k−1

xϕ



, (17d)

k−1



−

k−1



. (17e)

The main shortcoming of the pseudo-inverse ap-

proach when dealing with practical applications re-

sults from its relatively high sensitivity with respect

to the initialisation of the pseudo-inverse of the ma-

trix

xϕ

. This issue is not trivial and can lead to a

divergence of the overall algorithm. In order to appro-

priately initialise the expression



xϕ



†

, it is required

that the pseudo-inverseis computed ofﬂineafter an ar-

bitrary number, denoted α, of measurements is taken

and before the recursive algorithm commences oper-

ation.

Remark 1. It is noted that the uniqueness of



xϕ



†

in the case of recursive approaches is not always

guaranteed when utilising equations (17), see (Lin-

den, 2008) for further details. As a consequence, the

corresponding estimate of θ

may not represent the

optimal, in terms of the minimum variance, solution

to the overdetermined set of equations given by (10).

Approach based on Extended Instrumental Vari-

ables - REBPM

. An alternative to employing the

matrix pseudo-inverse lemma, an approach based on

the recursive extended instrumental variables tech-

nique, see (Friedlander, 1984), can be utilised in or-

der to obtain, albeit indirectly, a recursive update of



xϕ



†

. Deﬁne





xϕ



xϕ



−1

. (18)

In this approach the expression P

is updated recur-

sively, rather than the total pseudo-inverse



xϕ



†

The algorithm can be summarised as:

k−1

+ K



− φ

k−1



, (19a)

= P

k−1



+ φ

k−1



−1

, (19b)



−x

1 0



, (19c)





, (19d)

k− 1



k−1

xϕ



, (19e)



(k− 1)x

k−1



, (19f)

= P

k−1

− K

k−1

(19g)

with

xϕ

and

updated as in equations (17d) and

(17e), respectively. Since it is the expression P

which

is obtained recursively, hence, in order to calculate



xϕ



†

, for the recursive bias compensation equation

(16), an additional matrix producthas to be computed,

i.e.



xϕ



†

= P



xϕ



. (20)

Consequently, the pseudo-inverse of

xϕ

is obtained

in an indirect manner. Moreover, note that the recur-

sive algorithm (19a) requires an inverse of the matrix

of dimension 2 × 2 at each recursion. This, however,

does not signiﬁcantly increase the associated compu-

tational burden. On the other hand, the important ad-

vantages of this algorithm are that, ﬁrstly, it can be

easily initialised and, secondly, it is relatively insensi-

tive to the quality of the initial values. With reference

RECURSIVE EXTENDED COMPENSATED LEAST SQUARES BASED ALGORITHM FOR

ERRORS-IN-VARIABLES IDENTIFICATION

145

to (Friedlander, 1984), in the case of no a priori in-

formation the initialisation can be performed as

zφ

=µ



−n

)×n



, P

, θ

×1

. (21)

The scalar parameter µ allows the speed of conver-

gence to be adjusted, hence affects the ‘smoothness’

Θ (i.e. large value of µ corresponds to the slow

convergence and smooth parameters). Further algo-

rithmic details ensuring that the update of the matrix

, given by (19g), is (semi-) positive deﬁnite are ad-

dressed in (Friedlander, 1984).

4.3 Determination of Noise

Auto-correlation Elements

Since the matrix W

is sparse, in general, the compu-

tational effort involved in its pseudo-inverse is negli-

gible when compared to that of

xϕ

. Therefore, it is

the pseudo-inverse of

xϕ

which forms a crucial bot-

tleneck of the overall algorithm. Consequently, a re-

cursive computation of

is not considered here and

its estimate is determined ofﬂine at each recursion by

solving (14) in the LS sense, i.e.





†



−

xϕ

EBCLS



. (22)

5 SIMULATION STUDIES

This section addresses a numerical analysis of the

two proposed recursive realisations of the EBPM ap-

proach, namely REBPM

and REBPM

, when ap-

plied for the purpose of identifying a SISO discrete-

time LTI second order system within the EIV frame-

work. The system to be identiﬁed is described by

θ =



−1.5 0.7 1.0 0.5



(23)

with the input generated by

= 0.5u

k−1

+ β

, (24)

where β

is a white, zero mean sequence of unity vari-

ance. The input noise sequence is zero mean, white of

variance σ

˜u

and the coloured output noise sequence is

generated by

˜y

= 0.7˜y

k−1

+ γ

, (25)

where γ

is zero mean, white and of variance σ

. In

the case of both algorithms the instrumental vector is

based on the instruments proposed in (Ekman et al.,

2006), i.e. built from delayed inputs and delayed out-

puts, and utilised with n

= 10.

Table 1: Results of the estimation of model parameters and

auto-correlation elements of the noise sequences.

true EBPM REBPM

REBPM

SNR ≈ 11dB

−1.500 −1.501±0.041 −1.504±0.051 −1.494±0.023

0.700 0.701±0.045 0.705±0.056 0.694±0.024

1.000 0.998±0.039 0.996±0.045 1.001±0.038

0.500 0.500±0.072 0.495±0.083 0.508±0.051

˜u

0.100 0.100±0.054 0.095±0.065 0.124±0.052

˜y

(0) 3.922 3.273±2.349 2.647±3.902 3.834±1.376

˜y

(1) 2.745 2.168±1.938 1.631±3.275 2.618±1.174

˜y

(2) 1.922 1.540±0.949 1.250±1.612 1.721±0.715

− 0.001±0.001 0.004±0.005 0.001± 0.001

− 0.097±0.120 1.187±3.464 0.143± 0.197

Λ − 0 2 0

T − − 1.381±0.102 1.663±0.134

The robustness of the two algorithms is examined

via a Monte-Carlo simulation study comprising of

100 runs. The mean values of the estimates obtained

at the last recursion, i.e. for k = N are recorded and

compared with the corresponding results produced by

the ofﬂine EBPM. The the overall quality of the esti-

mators is assessed via the following two performance

criteria:



− θ



, e



− ρ



, (26)

where λ denotes the λ-th Monte-Carlo run. Prior to

the calculation of the performance indeces e

and e

the possible outliers are removed from the data. An

estimate is classiﬁed as an outlier if



> 10. The

number of outliers is denoted by Λ. Additionally, a

computation time in seconds, denoted T, is recorded.

The initial values of the parameters are set as fol-

lows: α = 50 for the REBPM

and µ = 100 for the

REBPM

. In order to provide a fair comparison, in

the case of the REBPM

, the bias compensation phase

is enabled from sample 50 onwards, although the ex-

pressions

and



xϕ



†

are recursively calculated

from the commencement of the algorithm. The values

of the noise parameters are chosen as σ

˜u

= 0.1 and

= 2.0. Consequently, the noise auto-correlation

vector is given by

ρ =



3.922 2.745 1.922 0.100



, (27)

which yields an approximately equal signal-to-noise

ratio (SNR) of around 11dB on both the input and

the output signals. The results expressed obtained

in terms of mean value ± standard deviation are pre-

sented in Table 1. It is observed that the mean values

of the model parameters, obtained by the algorithms,

see e

, are relatively accurate and close to the true val-

ues and are also characterised by acceptable standard

deviations. In the case of e

the estimates

ρ are rel-

atively less precise, especially those produced by the

REBPM

ICINCO 2009 - 6th International Conference on Informatics in Control, Automation and Robotics

146

In general, comparison of the two recursive real-

isations of the EBPM reveals that it is the REBPM

which produces the more accurate results overall.

Moreover, it is noted that in the case of the REBPM

the algorithm diverged twice, producing two outliers.

In terms of the computational burden, the time re-

quired by the REBPM

is slightly greater when com-

pared to that of REBPM

, i.e. the former technique is

faster by approximately 17% with respect to the latter

method.

In general, the experiments carried out seem to

suggest that the REBPM

is more advantageous than

the REBPM

due to a simpler initialisation, greater

robustness and an absence of convergence problems,

at least under the conditions considered here.

6 CONCLUSIONS

A recursive realisation of the extended bilinear

parametrisation method for the identiﬁcation of dy-

namical linear discrete-time time-invariant single-

input single-output errors-in-variables models has

been proposed. Two alternative approaches for the

online update of the pseudo-inverse of the extended

observation covariance matrix have been considered.

The ﬁrst approach is based on the pseudo-inverse ma-

trix lemma, whereas the second is constructed within

the framework of the extended instrumental variables

technique. For the cases considered, the two resulting

algorithms appear to be relatively robust and they are

also found to yield precise estimates of the model pa-

rameters. Results suggest that the instrumental vari-

ables based approach would appear to be the superior

of the two developed algorithms.

REFERENCES

Bj¨orck,

A. (1996). Numerical Methods for Least Squares

Problems. SIAM, Philadelphia.

Ekman, M. (2005a). Identiﬁcation of linear systems with

errors in variables using separable nonlinear least

squares. In Proc. of 16th IFAC World Congress,

Prague, Czech Republic.

Ekman, M. (2005b). Modeling and Control of Bilinear Sys-

tems: Applications to the Activated Sludge Process.

PhD thesis, Uppsala University, Sweden.

Ekman, M., Hong, M., and S¨oderstr¨om, T. (2006). A sep-

arable nonlinear least-squares approach for identiﬁca-

tion of linear systems with errors in variables. In 14th

IFAC Symp. on System Identiﬁcation, Newcastle, Aus-

tralia.

Feng, D., Zhang, H., Zhang, X., and Bao, Z. (2001). An

extended recursive least-squares algorithm. Signal

Proc., 81(5):1075–1081.

Friedlander, B. (1984). The overdetermined recursive in-

strumental variable method. IEEE Trans. on Auto-

matic Control, 29(4):353–356.

Ikenoue, M., Kanae, S., Yang, Z., and Wada, K.

(2008). Bias-compensation based method for errors-

in-variables model identiﬁcation. In Proc. of 17th

IFAC World Congress, pages 1360–1365, Seul, South

Korea.

Larkowski, T., Linden, J. G., Vinsonneau, B., and Burn-

ham, K. J. (2008). Identiﬁcation of errors-in-variables

systems via extended compensated least squares for

the case of coloured output noise. In The 19th Int.

Conf. on Systems Engineering, pages 71–76, Las Ve-

gas, USA.

Linden, J. G. (2008). Algorithms for recursive Frisch

scheme identiﬁcation and errors-in-variables ﬁlter-

ing. PhD thesis, Coventry University, UK.

Ljung, L. (1999). System Identiﬁcation - Theory for the

User. Prentice Hall PTR, New Jersey, USA, 2nd edi-

tion.

Mahata, K. (2007). An improved bias-compensation ap-

proach for errors-in-variables model identiﬁcation.

Automatica, 43(8):1339–1354.

Markovsky, I. and Van Huffel, S. (2007). Overview of to-

tal least-squares methods. Signal Proc., 87(10):2283–

2302.

S¨oderstr¨om, T. (2007). Errors-in-variables methods in sys-

tem identiﬁcation. Automatica, 43(6):939–958.

RECURSIVE EXTENDED COMPENSATED LEAST SQUARES BASED ALGORITHM FOR

ERRORS-IN-VARIABLES IDENTIFICATION

147