An algorithm for computing estimates for parameters of an ARMA-model from noisy measurements of inputs and outputs

An algorithm for computing estimates for parameters of an

ARMA-model from noisy measurements of inputs and outputs

Citation for published version (APA):

Vregelaar, ten, J. M. (1987). An algorithm for computing estimates for parameters of an ARMA-model from noisy measurements of inputs and outputs. (Memorandum COSOR; Vol. 8713). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/1987

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

Department of Mathematics and Computing Science

cas

OR-memorandum 87-13

An algorithm for computing estimates for parameters of an ARMA-model from noisy measurements

of inputs and outputs

by

J.M. ten Vregelaar

Eindhoven University of Technology

Department of Mathematics and Computing Science P.O. Box 513

5600 MB Eindhoven The Netherlands

Eindhoven, May 1987 The Netherlands

2

-An algorithm for computing estimates for parameters of an ARMA·modeI from noisy measurements of inputs and outputs

ABSTRACT

By applying some iterative algorithm a nonlinear minimization problem is solved in order to obtain estimates for the parameters of an ARMA-model. The algorithm contains a special Q - R decomposition using the structure of the problem in an optimal way.

1. Introduction

When one introduces some least squares estimation method for the parameters of an ARMA-model based on a set of noisy measurements of inputs and outputs, the following minimization problem is involved (cf. Ten Vregelaar '87):

(1.1a)

where

(1.1b) Here z is the known (s + r) (N + m) (column) vector of measurements,

e

is the s x (ps + (q + 1)sr) matrix of unknown parameters (cf. (3.2», and

(1.2a) with

D _N +.-_.-D _N T (D _N D _N 1')-1 _{• (1.2b)} The sN x (s + r) (N + m) matrix DN(8) has been introduced in order to denote the model equa-tion as

(1.3) where ~ is the vector of (unknown) true inputs and outputs.

For the moment we do not specify DN and , any further, but DN will be rowregular for all

a

anyway. Later on, we will choose DN and , conveniently. The integers r ,s ,m and N are sup-posed to be known, actually r ,8 and m will be small in comparison to N.

The solutions of (1.1) represent least squares estimates for the unknown matrix of parameters

a.

Since (1.1) has no closed-form solution, an iterative algorithm has to be applied to obtain esti-mates, provided an initial estimate for

a.

In this paper we will describe an algorithm to compute the object function J N and its vector of derivatives J N' for given

a.

The iteration procedure is supposed to be based on both I_N and I_{N ',}

If iN denotes any component of I_N',

(. representing

a!j

for any element of the matrix

e),

iN(a)

=

2zT DN+ DN PNi Z holds, with

I denoting the (8 + r) (N + m) identity matrix.

(1.4a)

4

-2. Application of a Q - R decomposition

We are concerned with the computation of I_N(8) and I_N '(8) for 8 known. Let us omit the subscript N in DN and P_{N •}

Since D is rowregular, DT has full column rank. Therefore it admits a Q - R decomposition

(2.1)

for some orthogonal matrix Q and sN square upper triangular regular R (it will be more con-venient here to construct a lower triangular R, as one can see from the next section).

In this section we will deal with the use of this Q - R decomposition for the computation of

J N(9) and J N' (9). In the next section we will describe how to perform efficiently the Q - R

decomposition for a suitable choice of D •

If

then

Hence

using the orthogonality of Q. Furthermore, from (1.1b)

IN

=

IIQl

T zlll •

II • U denoting the Euclidean norm.

Defining the sN -vector

implies

and from (1.4a)

(2.2)

(2.3)

(2.4)

(2.5)

From (2.2) and (2.5) we obtain

Introducing the sN -vector u := Ql T Z • we summarize.

By adding the column vector z to D, it follows

with v := Q2 T Z (sm + r(N + m) vector) for completeness. Then I

_N

and an arbitrary component

iN

of its derivative vector satisfy

IN = U U 112 •

iN

=

2

AT

D

(z - DT A) where R 1..= u .

(2.7)

(2.8)

(2.9) (2.10) (2. 11 a) (2. 11 b) (2. 11 c)

From the Q - R decomposition we only need R and u to compute _{I N and IN'. The solution of} A from (2.11c) can be done easily, since R turns out to be lower triangular.

6 -3. Performing the Q - R decomposition

According to relation (2.1) in Ten Vregelaar '87 we

can

denote the model equation as

where D=

• •

D (9)

C

=

0 ,

•

•

•

•

•

and the corresponding column

c

=

[~]

~o ~l

.

••• ~'"

.

•

(3.1b)

~N+m-l1

with 1']

=

1']1 (N + m) s -vector of true outputs

1']0

;N+m-1

and;

=

;1

;0

(N + m) r-vector of true inputs.

.

• ~o ~1 ••• ~'"

(3.1 a)

(3.1c)

The unknown parameterrnatrices to be estimated are collected in the s x (ps + (q + 1)sr) matrix

a,

(3.2) (the (Xi and ~j are resp. s x s and s x r matrices).

In (3.1b) we have

no

:= -Is • (XK := 0 for K > P and ~K := 0 for K > q . (3.3)

It is obvious that the sN x (s + r) (N + m) matrix D is rowregular, so its transposed DT has full column rank.

The integers p, q. r ,s are known, m := max (P. q) and N + m represents the number of meas-urements of inputs and outputs.

Introducing ai := Clj T and bi := ~l T (i = 0,1, ... ,m), we obtain for DT

.

.

, (s + r) (N + m) x sN (3.4)

consisting of two block-Toeplitz band matrices.

We split up the process of transfoffiling D T into two steps:

QDT~ [~l

(3.5a)

and

p

[~

1

~ [~

1 '

(3.5b)

where Q and P are orthogonal,

S 1,1

.

s=

,s(N+m)xsN (3.6a)

8

-•

,sNxsN (3.6b)

•

(the SiJ and rjJ are S x S blocks),

Defining QT := P !l (orthogonal), (3.5) implies

QT DT =

[~]

Slei2.1:

n

DT

=

[~].

G

_•

G

•

•

• am am

•

• • • ao

•

• ao •

•

am a m-2 am am-l S ", S ", SN+m,N S S '" _°N+m

°1

~ (3,7) DT ~ _~ ,

..

~ 0 0 0

@,

_{o ·}

•

@

•

bm- 1 • • • bm- 1 0 0

•

b_o 0 b_o

•

bm- 1 bm- 1 bm- 2

Some explanation is needed.

In

any

step 0,. Householder transfonnations are used in order to transfonn [::

1

into

~o].

where the transformed pairs

~

1

are called

[b~~1

1

(j

=

0.1 •...• m) with b -I

=

0 and b.

=

0 if

i>1.

Now, the block Toeplitz structure is maintained by applying each Householder transformation on all pairs [::

1

along the diagonals. moreover on all corresponding pairs [:;

J.

j

=

0,1, ... ,m.

The result of each OJ is a O-diagonal in the lower matrix and a block-row of S.

By using and maintaining the block Toeplitz structure the effect of any OJ is known com-pletely by just computing the effect of S Householder matrices (all (r + 1) x (r + 1» on all pairs

~;].

j O.I •...• m.

-

10-rabO O

]

In order to keep the blocks ao and thus Sk,k (Ie = 1.2 •... ,N) lower triangular, column S of ~

is attacked first. The nonzero diagonal element of the ao block will make the column (size r) of bo zero. After applying this Householder transformation on all pairs

~;].

J

= O.I ... m.

we

proceed wilh column

s -

1 of

[~:

l

etc.

~:p[~]= ~]

There will be orthogonal matrices

Pl>P

_{z •...}

,P

N acting on S as follows (cf. (3.6»:

$1.1 • $",+1,1 • _• ·s

.

_$

...

~ [~l·

r .• r ...

*

*

o

*

*

o

0

*

*

o

*

*

o

0

(3.8)

Defining

I'

=

PN PN-1".PZP

1 and P = diag(P, l) (orthogonal).

we obtain P

[g)

=

[~]

.

Each

Pi

consists of s (ms + 1) x (ms + 1) Householder transfonnations. The result is a zero (block) column and a block-row of R. The square appearing in the matrices has size ms x ms.

The first step of

PI

is making zero the very last column (size ms) of S by using the House-holder matrix detennined by this column and the (nonzero) diagonal element of the SN,N -block.

Proceed with the last but one column. again using the diagnonal element of the $N,N -block.

After $ Householders we continue with

P

2• now using the $N-l,N-l-block etc.

As a consequence. the diagonal blocks of R. the rj,i are lower triangular and regular. Obvi-ously the same holds for R itself.

The detennination of R is the result of the two steps described above. In order to compute IN

and J N' we also need the sN -vector u.

From (2.10) it is clear that u is detennined by perfonning the same transfonnations on the measurement vector

(cf. (3.1c»

as we have applied on the columns of D T •

Since R is a lower triangular and band matrix the computation of i... from R i...

=

u is easy to carry out.

Determination of IN and IN' is now straightforward (cf. (2. lla,b».

From (3.lb) the N nonzero elements of

D

are equal to 1

(. =

a!.

,any i

=

1,2, ... ,pS2 + (q + l)sr)

12 -4. Discussion

An upper bound for the number of operations to obtain R in (3.8) is

sCm + 1) (r + 1)2 (N + m) + sCm + 1) (ms + 1)2 N ,

so with respect to N the algorithm is of order 0 (N).

A Q - R algorithm for a general m

x

n matrix (m ;;:: n), based on Householder matrices, is of order n2(m - ;) (cf. Golub and Van Loan '83, p. 148).

From (3.1b) it follows DDT

=

dm • m-l where dk =

L

(CJ.;Cll+k + ~'~i~)' i=Q (4.1)

SO DDT is symmetric, block-Toeplitz and band matrix. In literature some algorithms are known to invert such matrices or to solve linear systems, containing those matrices (cf. Jain '78, Trench '74, Meek '83).

For our purpose, the special Q - R decomposition introduced here, will be superior to the

existing algorithms if r ,s , m < < N.

For the special case r = s = 1, the Householder transformations of step 1 in section 3 reduce to Givens transformations.

The matrix of second derivatives of J N (9) is not used in the iterative algorithm to find the minimum of IN because of its form (cf. Ten Vregelaar '87)

(PJ_N

aejae.

= 2zT_{[PI(D i l(DD}T_)-lDi _{pI - D+(DiD+Di + DiD+Dj)Pl - D+Dipl(Dil(D+l] z , (4.2)}

.

a

withD' ' : : : - D

_{. ae}

The minimization problem is solved by using the Broyden-Fletcher-Goldfarb-Shanno formula (cf. Scales '85 p. 89-90).

If some relations exist between the parameters 0.1> ... ,

a.

_{p •} ~o ... ~q we have to solve a con-strained optimization problem. The described algorithm to compute the object function and its gradient can be applied for this case as well.

References

Golub and Van Loan

'83

Golub G.H., Van Loan, C.F.,

14

-Matrix Computations, North Oxford Academic, Oxford, 1983.

Jain '78 Jain, A.K .•

Fast inversion of banded Toeplitz matrices by circular decompositions, IEEE Trans. Acoust, Speech, Signal Processing, vol. ASSP-26, 2, pp. 121-126, 1978.

Meek, '83

Meek,D.S.,

The inverses of Toeplitz band matrices, Linear Algebra Appl. 49, pp. 117-129, 1983.

Scales '85

Scales, L.E.,

Introduction to Non-Linear Optimization, MacMillan,London,1985.

Ten

Vre~elaar '87 Ten Vregelaar, J.M.,

Asymptotic properties of an estimation method for an ARMA-model with noisy measurements of inputs and outputs, Memorandum-COSOR, Techn. Univ. Eindhoven, in preparation.

Trench '74

Trench, W.F.,

Inversion of Toeplitz band matrices,