# Weight Least Squares Algorithm for Rational Models with Outliers

Authors: Lixin Lv; Donglei Lu
Journal: Complexity
 (2020)

---

## Abstract

A weight least squares algorithm is developed for rational models with outliers in this paper. Different weights are assigned for each cost function, and by calculating the derivatives of these cost functions, the parameter estimates can be estimated. Compared with the traditional least squares algorithm, the proposed algorithm can remove the bad effect caused by the outliers, thus has more accurate parameter estimates. A simulation example is proposed to validate the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

System identification can be roughly divided into two kinds: the structure identification and the parameter estimation [1–4]. Compared to the parameter estimation, structure identification is more difficult, because one should choose the available products from numerous product pool. Once the structure is determined, parameter estimation is involved [5–7]. Recently, a lot of parameter estimation algorithms are proposed, including the least squares (LS) algorithm [8, 9], the gradient algorithm [10, 11], the particle swarm optimization algorithm, and the expectation maximization algorithm [12, 13]. Among these algorithms, the LS algorithm is to find a vector that is a local minimizer to a function that is a sum of squares; thus, it is the simplest algorithm and is widely used.Although the LS algorithm is applied to many different kinds of nonlinear models, e.g., bilinear models [14, 15], polynomial nonlinear models [16, 17], hard nonlinear models, and rational nonlinear models [18, 19], it also brings a challenging problem: heavy computational effort caused by the inverse matrix calculation. In order to overcome this difficulty, some online algorithms are proposed. For example, Liu and Ding proposed an auxiliary model-based recursive generalized least squares algorithm for multivariate output-error autoregressive system [20], Zhu provided an implicit least squares algorithm for nonlinear rational models [21], Wang et al. developed a recursive LS (RLS) algorithm for Hammerstein models [22]. However, all the systems in the above literature have normal data, when some data of the systems are abnormal, those algorithms may be invalid.Processes in industry usually suffer from outliers in measurement data [23–25], which make the identification of the process a challenging problem. Thus, a reliable estimation of the systems with outliers is essential to efficient process identification [26, 27]. Zhao et al. provided a VB approach for ARX models with outliers and assumed that all the variables are available and the precise interval of time-delays is known a priori [28]. Jin et al. proposed a conventional approach to nonlinear process models with outliers [29]. Both these two algorithms are off-line algorithms, which update the parameters through all the collected data, and thus have heavy computational efforts especially for large-scale systems.In this paper, a weight LS algorithm is proposed for rational models with outliers. By introducing different weights for each cost function, the cost functions of the outliers can be neglected. Furthermore, to reduce the computational efforts, a weight recursive least squares (W-RLS) is also derived. Both these two algorithms can yield the optimal parameter estimates. Briefly, the paper is listed as follows. Section2 introduces the rational model and the traditional LS algorithm. Section 3 develops the weight LS algorithm and the weight RLS algorithm. Section 4 provides a simulation example. Finally, Section 5 gives a conclusion of this paper.


## 2. Rational Model with Outliers

The rational model with outliers is written by(1)y1t=αt+e1tβt,(2)y2t=αt+e2tβt,where y1t is the normal output, y2t is the abnormal output, e1t a normal stochastic white noise with zero mean and variance σ2, e2t the outliers of the rational model, whose noise to signal ratio is larger, and αt and βt are the structures of the rational model that can be expressed as(3)αt=φTtθa,βt=ψTtθb.The information vectorsφt=φ1t,…,φntT and ψt=ψ1t,…,ψmtT are the products of past inputs ut−1,ut−2,… and past outputs yt−1,yt−2,…, and the structures of φt and ψt are known in prior, θa and θb are the unknown parameters to be estimated and can be expressed as(4)θa=a1,a2,…,anT,θb=b1,b2,…,bmT.In application, the number of the abnormal data is far smaller than that of the normal data. Assume that we have collectedL input-output data, the number of the normal data is L1 and the abnormal data is L2, L1+L2=L,L1≫L2. Without loss of generality, assume that b1=1, then the rational model is transformed into(5)y¯t=ϕtυ+eit,i=1,2,y¯t=ytψ1t,ϕt=φ1t,…,φnt,−ytψ2t,…,−ytψmtT,υ=a1,…,an,b2,…,bmT.CollectL input and output data and define(6)Y¯L=y¯L,y¯L−1,…,y¯1T∈ℝL,ΦL=ϕL,…,ϕ1T∈ℝL×m+n−1,EL=eiL,…,ei1T∈ℝL,where i=1 or 2. We have(7)Y¯L=ΦLυ+EL.Define the cost function(8)Jυ=12Y¯L−ΦLυ2,then the following LS algorithm can be obtained(9)υ^=ΦTLΦL−1ΦTLY¯L.Remark 1.
In the LS algorithm, all the data are used to update the parameter estimates; thus, the outliers cause extremely bad influence to the estimates.


## 3. Weight Based Least Squares Algorithm

In order to get more accurate parameter estimates, one should weak/eliminate the bad influence caused by the outliers. Therefore, a weight least squares algorithm is proposed in this section.


### 3.1. Weight LS Algorithm

Rewrite the cost function as(10)Jυ=12Y¯L−ΦLυ2,which can be transformed into(11)Jυ=12∑t=1Ly¯t−ϕtυ2.Clearly, each cost function plays the same role in estimating the parameters when using the traditional LS algorithm. For this reason, we can introduce different weights for each cost function, e.g., small weight for the cost function of the outliers and large weight for the cost function of the normal data.It follows that the new cost function can be written by(12)Jυ=12∑t=1LλL−ty¯t−ϕtυ2,0<λ<1.Taking the derivative ofJυ with respect to υ and then equating it to zero give(13)υ^=∑t=1LλL−tϕtϕTt−1∑t=1LλL−tϕty¯t.Remark 2.
Based on equations (12) and (13), we can get that each cost function has different weights, and the latest data have larger weights.
Taking the conditional expectation on both sides of (13) yields(14)Eυ^=E∑t=1LλL−tϕtϕTt−1∑t=1LλL−tϕty¯t=E∑t=1LλL−tϕtϕTt−1∑t=1LλL−tϕtϕtυ+et=υ+∑t=1LλL−tϕtϕTt−1E∑t=1LλL−tϕtet.
Sinceet is independent of ϕt, which means that(15)∑t=1LλL−tϕtϕTt−1E∑t=1LλL−tϕtet=0.
Therefore, we have(16)Eυ^=υ.
Unfortunately, the W-LS algorithm needs to perform the inverse matrix calculation∑t=1LλL−tϕtϕTt−1, if the order m+n−1 is large, and the computational effort is heavy. To overcome this difficulty, we derive a weight recursive least squares algorithm in the following subsection.


### 3.2. Weight Recursive Least Squares Algorithm

Define(17)P−1L=∑t=1LλL−tϕtϕTt=λL−1ϕ1ϕT1+λL−2ϕ2ϕT2+⋯+λϕL−1ϕTL−1+ϕLϕTL,P−1L−1=∑t=1L−1λL−1−tϕtϕTt=λL−2ϕ1ϕT1+λL−3ϕ2ϕT2+⋯+ϕL−1ϕTL−1.The above two equations give rise to(18)P−1L=λP−1L−1+ϕLϕTL.Lemma 1.
[30]. Assume that A∈ℝn×n, B∈ℝn×m, and C∈ℝm×n, then the following equality holds:(19)A+BC−1=A−1−A−1BI+CA−1B−1CA−1.
According to Lemma1, PL can be written by(20)PL=P−1L−1=λP−1L−1+ϕLϕTL−1=1λPL−1−1/λ2PL−1ϕLϕTLPL−11+1/λϕTLPL−1ϕL=1λPL−1−PL−1ϕLϕTLPL−1λ2+λϕTLPL−1ϕL.The parameter estimates by using the W-LS algorithm are rewritten by(21)υ^L=∑t=1LλL−tϕtϕTt−1∑t=1LλL−tϕty¯t=PL∑t=1LλL−tϕty¯t=PL∑t=1L−1λL−tϕty¯t+ϕLy¯L=PLλ∑t=1L−1λL−1−tϕty¯t+ϕLy¯L=PLλ∑t=1L−1λL−1−tϕtϕTt∑t=1L−1λL−1−tϕtϕTt−1∑t=1L−1λL−1−tϕty¯t+ϕLy¯L=PLλ∑t=1L−1λL−1−tϕtϕTtυ^L−1+ϕLy¯L=PLλP−1L−1υ^L−1+ϕLy¯L.According to equation (18), we have(22)λP−1L−1=P−1L−ϕLϕTL.Therefore, equation (21) is simplified to(23)υ^L=PLP−1L−ϕLϕTLυ^L−1+ϕLy¯L=υ^L−1+PLϕLy¯L−ϕTLυ^L−1.The parameter estimatesυ^t by using the weight recursive least squares (W-RLS) algorithm are listed as follows:(24)υ^t=υ^t−1+Ptϕty¯t−ϕTtυ^t−1,(25)ϕt=−ytψ2t,…,−ytψmt,φ1t,…,φntT,(26)y¯t=ytψ1t,(27)Pt=1λPt−1−Pt−1ϕtϕTtPt−1λ2+λϕTtPt−1ϕt,0<λ<1.The W-RLS algorithm consists of the following iterations:(1)
Letut=0,yt=0, t≤0, and give a small positive number ε and a positive number 0<λ<1.(2)
Lett=1υ^0=1/p0 with 1 being a column vector whose entries are all unity and p0=106, P0=p0I.(3)
Collectyt and ut.(4)
Formy¯t by equation (26).(5)
Formϕt according to equation (25).(6)
ComputePt by equation (27).(7)
Update the parameter estimation vectorυ^t by equation (24).(8)
Compareυ^t and υ^t−1: if υ^t−υ^t−1/υ^t⩽ε, then terminate the procedure and obtain υ^t; otherwise, increase t by 1 and go to step 3.


### 3.3. Property of the W-RLS Algorithm

Lemma 2.
Assume thatA∈ℝn×n and B∈ℝn×n are both symmetric positive semidefinite matrices, then A+B is also a symmetric positive semidefinite matrix.Proof.
SinceA and B are symmetric positive semidefinite matrices, for a random nonzero vector X∈ℝn, the following inequalities hold:(28)XTAX≥0,XTBX≥0.
It gives rise to(29)XTA+BX≥0,which means that A+B is also a symmetric positive semidefinite matrix.Lemma 3.
Assume thatA∈ℝn×n and B∈ℝn×n are both positive definite symmetric matrices. The largest eigenvalue of the matrix A is a, and the largest eigenvalue of the matrix B is b. Then the largest eigenvalue of the matrix A+B is a+b.Proof.
SinceA and B are positive definite symmetric matrices, and the largest eigenvalues are a and b, respectively. Then we have(30)aI−A≥0,bI−B≥0.
Clearly,aI−A and bI−B are symmetric positive semidefinite matrices. We have(31)a+bI−A+B≥0.
Therefore, the largest eigenvalue of the matrixA+B is a+b.Theorem 2.
For the rational model in (2), the parameter estimate υ^t by using the W-RLS algorithm is expressed by equations (24)–(27). Then we have(32)limt⟶∞υ^t=υ.Proof.
Rewriteυ^t as follows:(33)υ^t=υ^t−1+Ptϕty¯t−ϕTtυ^t−1,and the following equation can be obtained(34)υ^t−υ=υ^t−1−υ+PtϕtϕTtυ+vt−ϕTtυ^t−1.
Taking the conditional expectation on both sides of (34) gets(35)Eυ^t−υ=EI−PtϕtϕTtυ^t−υ+EPtϕtet.
e
t is Gaussian noise and satisfies(36)Eet=0.
It follows that(37)Eυ^t−υ=EI−PtϕtϕTtυ^t−υ.
Since(38)Pt=P−1t−1=∑i=1tλt−iϕiϕTi−1=∑i=1t−1λt−iϕiϕTi+ϕtϕTt−1.
According to Lemmas2 and 3 we have(39)I−PtϕtϕTt<1.
Assume that the largest eigenvalue of the matrixI−PtϕtϕTt is λt, and let(40)λmax=λ1,λ2,…,.
Clearly,0<λmax<1.
It follows that(41)limt⟶∞υ^t−υ=limt⟶∞λmaxtυ^0−υ=0.
Then the parameter estimate by using the W-RLS algorithm is asymptotic convergent.Remark 3.
Unlike the W-LS algorithm, the W-RLS algorithm avoids the inverse matrix calculation, thus has less computational efforts.Remark 4.
At the sampling instantt, the W-RLS algorithm uses the data up to and including time t, while the stochastic gradient algorithm only uses the current data yt,ut; thus, the W-RLS algorithm has a quicker convergence rate than that of the stochastic gradient algorithm.


## 3.1. Weight LS Algorithm

Rewrite the cost function as(10)Jυ=12Y¯L−ΦLυ2,which can be transformed into(11)Jυ=12∑t=1Ly¯t−ϕtυ2.Clearly, each cost function plays the same role in estimating the parameters when using the traditional LS algorithm. For this reason, we can introduce different weights for each cost function, e.g., small weight for the cost function of the outliers and large weight for the cost function of the normal data.It follows that the new cost function can be written by(12)Jυ=12∑t=1LλL−ty¯t−ϕtυ2,0<λ<1.Taking the derivative ofJυ with respect to υ and then equating it to zero give(13)υ^=∑t=1LλL−tϕtϕTt−1∑t=1LλL−tϕty¯t.Remark 2.
Based on equations (12) and (13), we can get that each cost function has different weights, and the latest data have larger weights.
Taking the conditional expectation on both sides of (13) yields(14)Eυ^=E∑t=1LλL−tϕtϕTt−1∑t=1LλL−tϕty¯t=E∑t=1LλL−tϕtϕTt−1∑t=1LλL−tϕtϕtυ+et=υ+∑t=1LλL−tϕtϕTt−1E∑t=1LλL−tϕtet.
Sinceet is independent of ϕt, which means that(15)∑t=1LλL−tϕtϕTt−1E∑t=1LλL−tϕtet=0.
Therefore, we have(16)Eυ^=υ.
Unfortunately, the W-LS algorithm needs to perform the inverse matrix calculation∑t=1LλL−tϕtϕTt−1, if the order m+n−1 is large, and the computational effort is heavy. To overcome this difficulty, we derive a weight recursive least squares algorithm in the following subsection.


## 3.2. Weight Recursive Least Squares Algorithm

Define(17)P−1L=∑t=1LλL−tϕtϕTt=λL−1ϕ1ϕT1+λL−2ϕ2ϕT2+⋯+λϕL−1ϕTL−1+ϕLϕTL,P−1L−1=∑t=1L−1λL−1−tϕtϕTt=λL−2ϕ1ϕT1+λL−3ϕ2ϕT2+⋯+ϕL−1ϕTL−1.The above two equations give rise to(18)P−1L=λP−1L−1+ϕLϕTL.Lemma 1.
[30]. Assume that A∈ℝn×n, B∈ℝn×m, and C∈ℝm×n, then the following equality holds:(19)A+BC−1=A−1−A−1BI+CA−1B−1CA−1.
According to Lemma1, PL can be written by(20)PL=P−1L−1=λP−1L−1+ϕLϕTL−1=1λPL−1−1/λ2PL−1ϕLϕTLPL−11+1/λϕTLPL−1ϕL=1λPL−1−PL−1ϕLϕTLPL−1λ2+λϕTLPL−1ϕL.The parameter estimates by using the W-LS algorithm are rewritten by(21)υ^L=∑t=1LλL−tϕtϕTt−1∑t=1LλL−tϕty¯t=PL∑t=1LλL−tϕty¯t=PL∑t=1L−1λL−tϕty¯t+ϕLy¯L=PLλ∑t=1L−1λL−1−tϕty¯t+ϕLy¯L=PLλ∑t=1L−1λL−1−tϕtϕTt∑t=1L−1λL−1−tϕtϕTt−1∑t=1L−1λL−1−tϕty¯t+ϕLy¯L=PLλ∑t=1L−1λL−1−tϕtϕTtυ^L−1+ϕLy¯L=PLλP−1L−1υ^L−1+ϕLy¯L.According to equation (18), we have(22)λP−1L−1=P−1L−ϕLϕTL.Therefore, equation (21) is simplified to(23)υ^L=PLP−1L−ϕLϕTLυ^L−1+ϕLy¯L=υ^L−1+PLϕLy¯L−ϕTLυ^L−1.The parameter estimatesυ^t by using the weight recursive least squares (W-RLS) algorithm are listed as follows:(24)υ^t=υ^t−1+Ptϕty¯t−ϕTtυ^t−1,(25)ϕt=−ytψ2t,…,−ytψmt,φ1t,…,φntT,(26)y¯t=ytψ1t,(27)Pt=1λPt−1−Pt−1ϕtϕTtPt−1λ2+λϕTtPt−1ϕt,0<λ<1.The W-RLS algorithm consists of the following iterations:(1)
Letut=0,yt=0, t≤0, and give a small positive number ε and a positive number 0<λ<1.(2)
Lett=1υ^0=1/p0 with 1 being a column vector whose entries are all unity and p0=106, P0=p0I.(3)
Collectyt and ut.(4)
Formy¯t by equation (26).(5)
Formϕt according to equation (25).(6)
ComputePt by equation (27).(7)
Update the parameter estimation vectorυ^t by equation (24).(8)
Compareυ^t and υ^t−1: if υ^t−υ^t−1/υ^t⩽ε, then terminate the procedure and obtain υ^t; otherwise, increase t by 1 and go to step 3.


## 3.3. Property of the W-RLS Algorithm

Lemma 2.
Assume thatA∈ℝn×n and B∈ℝn×n are both symmetric positive semidefinite matrices, then A+B is also a symmetric positive semidefinite matrix.Proof.
SinceA and B are symmetric positive semidefinite matrices, for a random nonzero vector X∈ℝn, the following inequalities hold:(28)XTAX≥0,XTBX≥0.
It gives rise to(29)XTA+BX≥0,which means that A+B is also a symmetric positive semidefinite matrix.Lemma 3.
Assume thatA∈ℝn×n and B∈ℝn×n are both positive definite symmetric matrices. The largest eigenvalue of the matrix A is a, and the largest eigenvalue of the matrix B is b. Then the largest eigenvalue of the matrix A+B is a+b.Proof.
SinceA and B are positive definite symmetric matrices, and the largest eigenvalues are a and b, respectively. Then we have(30)aI−A≥0,bI−B≥0.
Clearly,aI−A and bI−B are symmetric positive semidefinite matrices. We have(31)a+bI−A+B≥0.
Therefore, the largest eigenvalue of the matrixA+B is a+b.Theorem 2.
For the rational model in (2), the parameter estimate υ^t by using the W-RLS algorithm is expressed by equations (24)–(27). Then we have(32)limt⟶∞υ^t=υ.Proof.
Rewriteυ^t as follows:(33)υ^t=υ^t−1+Ptϕty¯t−ϕTtυ^t−1,and the following equation can be obtained(34)υ^t−υ=υ^t−1−υ+PtϕtϕTtυ+vt−ϕTtυ^t−1.
Taking the conditional expectation on both sides of (34) gets(35)Eυ^t−υ=EI−PtϕtϕTtυ^t−υ+EPtϕtet.
e
t is Gaussian noise and satisfies(36)Eet=0.
It follows that(37)Eυ^t−υ=EI−PtϕtϕTtυ^t−υ.
Since(38)Pt=P−1t−1=∑i=1tλt−iϕiϕTi−1=∑i=1t−1λt−iϕiϕTi+ϕtϕTt−1.
According to Lemmas2 and 3 we have(39)I−PtϕtϕTt<1.
Assume that the largest eigenvalue of the matrixI−PtϕtϕTt is λt, and let(40)λmax=λ1,λ2,…,.
Clearly,0<λmax<1.
It follows that(41)limt⟶∞υ^t−υ=limt⟶∞λmaxtυ^0−υ=0.
Then the parameter estimate by using the W-RLS algorithm is asymptotic convergent.Remark 3.
Unlike the W-LS algorithm, the W-RLS algorithm avoids the inverse matrix calculation, thus has less computational efforts.Remark 4.
At the sampling instantt, the W-RLS algorithm uses the data up to and including time t, while the stochastic gradient algorithm only uses the current data yt,ut; thus, the W-RLS algorithm has a quicker convergence rate than that of the stochastic gradient algorithm.


## 4. Example

Consider the rational model proposed in [31],(42)yt=0.2yt−1+0.1yt−1ut−1+ut−1+et1+y2t−1+y2t−2.Then, one can get(43)yt=0.2yt−1+0.1yt−1ut−1+ut−1−yty2t−1−yty2t−2+et,where ut is a input sequence with zero mean and unit variance et is a white noise with zero mean and variance σ2. In simulation, we collect 1800 input-output data, and the data from 500−600 are outliers. The noise to signal ratio is 260.2% when t∈500,600, and 2.78% when t∈0,499∪601,1800. The simulation data are shown in Figure 1, t∈1,1000, it shows the noise and output data from 500 to 600 are outliers.Figure 1
The simulation data.
(a)
(b)
(c)Apply the LS and W-LS algorithms to the rational modelλ=0.9. The parameter estimates and their estimation errors τ:=υ^−υ/υ are shown in Table 1.Table 1
The parameter estimates and errors.
Algorithms
a
1
a
2
a
3
b
2
b
3
τ
%
LS
0.15181
0.07996
0.79217
0.51845
0.51375
41.06137
W-LS
0.12532
0.12418
1.00963
0.98697
0.98437
4.67610
True values
0.20000
0.10000
1.00000
1.00000
1.00000Next, use the RLS and the W-RLS algorithms for this rational modelsλ=0.9, the parameter estimates and their estimation errors τ:=υ^−υ/υ are shown in Table 2 and Figure 2.Table 2
The parameter estimates and errors.
Algorithms
t
a
1
a
2
a
3
b
2
b
3
τ
%
RLS
100
0.19112
0.12398
0.95403
0.89402
0.84710
11.07001
300
0.18911
0.11181
0.95851
0.89135
0.91121
8.42857
600
0.02324
−0.01496
0.46701
−0.37791
−0.29143
113.00607
900
0.04782
0.00113
0.49319
−0.31770
−0.24903
108.43468
1200
0.05971
0.01433
0.51901
−0.25419
−0.21230
104.03382
1500
0.07107
0.02377
0.54218
−0.19289
−0.17885
99.91320
1800
0.08303
0.02983
0.57281
−0.12658
−0.12806
94.83039
W-RLS
100
0.19100
0.12362
0.94507
0.88143
0.82101
12.77180
300
0.18310
0.09174
0.95825
0.88932
0.92104
8.21469
600
−0.07442
−0.12394
0.34539
−0.54102
−0.30873
123.36106
900
0.10473
0.04088
0.61333
−0.03655
0.05515
83.55360
1200
0.17230
0.08650
0.91166
0.78224
0.75087
19.68922
1500
0.17447
0.09915
0.92978
0.83961
0.77437
16.41848
1800
0.18480
0.09849
0.97901
0.95148
0.93699
4.79017
True values
0.20000
0.10000
1.00000
1.00000
1.00000Figure 2
The parameter estimation errorsτ versus.t.Then the following findings can be obtained: (1) Table1 shows that the W-LS algorithm is more effective than the LS algorithm; (2) Table 2 and Figure 2 demonstrate that the parameter estimates by using the W-RLS algorithm are more accurate than those by using the RLS algorithm.


## 5. Conclusions

This paper proposes a weight least squares algorithm for rational models with outliers. By assigning the weights for each data set, the proposed algorithm can obtain more accurate parameter estimates when compared with the traditional least squares algorithm. Furthermore, a weight recursive least squares algorithm is derived to decrease the computational efforts. The convergence analysis and simulation example show that the proposed algorithm is effective.



---

*Source: 8963691-2020-04-14.xml*