# Exhaustive Search and Power-Based Gradient Descent Algorithms for Time-Delayed FIR Models

Authors: Hua Chen; Yuejiang Ji
Journal: Complexity
 (2022)

---

## Abstract

In this study, two modified gradient descent (GD) algorithms are proposed for time-delayed models. To estimate the parameters and time-delay simultaneously, a redundant rule method is introduced, which turns the time-delayed model into an augmented model. Then, two GD algorithms can be used to identify the time-delayed model. Compared with the traditional GD algorithms, these two modified GD algorithms have the following advantages: (1) avoid a high-order matrix eigenvalue calculation, thus, are more efficient for large-scale systems; (2) have faster convergence rates, therefore, are more practical in engineering practices. The convergence properties and simulation examples are presented to illustrate the efficiency of the two algorithms.


---

## Body

## 1. Introduction

System identification plays an important role in control theory and application [1–3]. When the model of a dynamic system is established, one can design robust controllers for such a model to predict its dynamics in the future. There exist many identification algorithms, for example, the least squares (LS) algorithm [4, 5], the gradient descent (GD) algorithm [6, 7], and the particle swarm optimization (PSO) algorithm [8, 9]. When the considered model has a high order, the LS algorithm and the PSO algorithm are inefficient for their heavy computational efforts [10–12]. The GD algorithm has few computational efforts, but with slow convergence rates [13, 14]. To increase the convergence rate of the GD algorithm, two ways are usually performed: (1) design a more suitable direction [15–17]; (2) calculate a better step size [18, 19]. In [20], the best step size of a GD algorithm is given, which involves the eigenvalue calculation. For a high-order matrix, computing its eigenvalues is challenging. To deal with this problem, plenty of suboptimal step size calculating methods are developed, for example, the stochastic GD algorithm [21, 22], the forgetting factor GD algorithm [18, 19], the projection algorithm, and the steepest GD algorithm [23, 24]. Although these algorithms can increase the convergence rates, they are all sensitive to the considered model. That is, one should design different step sizes for different kinds of models.Time delay is normal in engineering practices. The data of a dynamic system are usually collected by a sensor and then transmitted via a communication channel; they may encounter time-delay due to network congestion [25, 26]. For the time-delayed model identification, Chen proposed a redundant rule-based off-line algorithm that can estimate the parameters and time delay simultaneously [27]. Since the off-line algorithm cannot update the parameters with newly arrived data, Zhang et al. developed a redundant rule-based recursive LS (RLS) algorithm for bilinear time-delayed systems, the RLS algorithm is an online algorithm [28]. This paper focuses on time-delayed model identification and aims to develop some novel identification algorithms which have fast convergence rates and less computational efforts.Inspired by the PSO algorithm and the power method [29], we propose two modified GD algorithms for time-delayed models: one is an exhaustive search method which chooses the step size based on the PSO algorithm, and the other is the power-based GD algorithm which computes the step size using power method. These two algorithms can get a better step size in each iteration without eigenvalue calculation. Therefore, the proposed algorithms have faster convergence rates when compared with the traditional GD algorithm.This paper is organized as follows: Section2 describes the time-delayed model and the traditional LS and GD algorithms. In Section 3, two modified GD algorithms are proposed. Section 4 proves the convergence properties of the two algorithms. Section 5 gives two simulation examples. Finally, conclusions are presented in Section 6.


## 2. Problem Statement

First, some notations are denoted as follows:I denotes an identity matrix of the appropriate sizes; X means the norm of the matrix X and is written as X:=ρXTX; ρXTX stands for the spectral radius of the matrix XTX; the superscript T is defined as the matrix transpose; λmaxM and λminM denote the maximum and minimum eigenvalues of a matrix M, respectively.


### 2.1. Time-Delayed Model

Consider the following time-delayed model,(1)yt=∑i=1Ngiut−i−τ+vt,where yt and ut are the output and input, respectively; vt is a Gaussian white noise, and satisfies vt∼N0,δ2; gi,i=1,…,N are the unknown parameters need to be estimated; τ is an unknown time delay.Since the time-delayτ is unknown, the corresponding information vector ut−1−τ,…,ut−N−τ of yt is unavailable which leads to the traditional GD algorithm being impossible. To deal with this dilemma, we use the redundant rule method. Assume that the upper bound of the time delay is M , and this assumption is rational and feasible. For example, when using the RIP protocol in the network, the maximum flop is 16.Rewrite the time-delayed model as follows:(2)yt=g1ut−1+g2ut−2+⋯+gτut−τ+gτ+1ut−τ−1+⋯+gτ+Nut−τ−N+gτ+N+1ut−τ−N−1+⋯+gτ+N+Mut−N−M+vt.Define the parameter vectorG and the information vector ψ as(3)G=g1,…,gτ,gτ+1,…,gτ+N,gτ+N+1,…,gN+MT∈ℝM+N,ψt=ut−1,…,ut−τ,ut−τ−1,…,ut−τ−N,ut−τ−N−1,…,ut−N−MT∈ℝM+N.The augmented parameter vector is decomposed into the following three parts:(4)1.redundant part:G1=g1,…,gτT,2.true part:G2=gτ+1,…,gτ+NT,3.redundant part:G3=gτ+N+1,…,gN+MT.Remark 1.
Since the two corresponding information vectors of the parameter vectors play a less role in the output, the redundant partsG1 and G2 are both zero vectors. If the parameter estimates of G converge to the true values, the two redundant parts equal zero vectors, and then we can obtain the time-delay estimates based on this special structure.


### 2.2. LS and GD Algorithms

Rewrite the augmented model of the time-delayed model as(5)yt=ψTtG+vt.CollectL sets of input and output data and define(6)YL=yL,yL−1,…,y1T∈ℝL,ΦTL=ψL,ψL−1,…,ψ1T∈ℝL×P,P=M+N,VL=vL,vL−1,…,v1T∈ℝL.It gives rise to(7)YL=ΦTLG+VL.Define the cost function as follows:(8)JG=12YL−ΦTLG2.Using the LS algorithm to estimate the parameters, it follows that(9)G∧=ΦLΦTL−1ΦLYL.The LS algorithm should perform a matrix inverse calculation which may lead to heavy computational efforts, especially for large-scale systems, e.g.,P is large.To avoid the matrix inverse calculation, the traditional GD (T-GD) algorithm is introduced [20],(10)G∧k=G∧k−1+γkΦLYL−ΦTLGk−1,0<γk〈2λmaxΦLΦTL.The T-GD algorithm does not need to compute the inverse of the information matrixΦLΦTL, but it requires calculating the eigenvalues of the information matrix to choose a suitable step size to keep the T-GD algorithm convergent. When ΦLΦTL has a high order, computing its eigenvalues is also a challenging problem.


## 2.1. Time-Delayed Model

Consider the following time-delayed model,(1)yt=∑i=1Ngiut−i−τ+vt,where yt and ut are the output and input, respectively; vt is a Gaussian white noise, and satisfies vt∼N0,δ2; gi,i=1,…,N are the unknown parameters need to be estimated; τ is an unknown time delay.Since the time-delayτ is unknown, the corresponding information vector ut−1−τ,…,ut−N−τ of yt is unavailable which leads to the traditional GD algorithm being impossible. To deal with this dilemma, we use the redundant rule method. Assume that the upper bound of the time delay is M , and this assumption is rational and feasible. For example, when using the RIP protocol in the network, the maximum flop is 16.Rewrite the time-delayed model as follows:(2)yt=g1ut−1+g2ut−2+⋯+gτut−τ+gτ+1ut−τ−1+⋯+gτ+Nut−τ−N+gτ+N+1ut−τ−N−1+⋯+gτ+N+Mut−N−M+vt.Define the parameter vectorG and the information vector ψ as(3)G=g1,…,gτ,gτ+1,…,gτ+N,gτ+N+1,…,gN+MT∈ℝM+N,ψt=ut−1,…,ut−τ,ut−τ−1,…,ut−τ−N,ut−τ−N−1,…,ut−N−MT∈ℝM+N.The augmented parameter vector is decomposed into the following three parts:(4)1.redundant part:G1=g1,…,gτT,2.true part:G2=gτ+1,…,gτ+NT,3.redundant part:G3=gτ+N+1,…,gN+MT.Remark 1.
Since the two corresponding information vectors of the parameter vectors play a less role in the output, the redundant partsG1 and G2 are both zero vectors. If the parameter estimates of G converge to the true values, the two redundant parts equal zero vectors, and then we can obtain the time-delay estimates based on this special structure.


## 2.2. LS and GD Algorithms

Rewrite the augmented model of the time-delayed model as(5)yt=ψTtG+vt.CollectL sets of input and output data and define(6)YL=yL,yL−1,…,y1T∈ℝL,ΦTL=ψL,ψL−1,…,ψ1T∈ℝL×P,P=M+N,VL=vL,vL−1,…,v1T∈ℝL.It gives rise to(7)YL=ΦTLG+VL.Define the cost function as follows:(8)JG=12YL−ΦTLG2.Using the LS algorithm to estimate the parameters, it follows that(9)G∧=ΦLΦTL−1ΦLYL.The LS algorithm should perform a matrix inverse calculation which may lead to heavy computational efforts, especially for large-scale systems, e.g.,P is large.To avoid the matrix inverse calculation, the traditional GD (T-GD) algorithm is introduced [20],(10)G∧k=G∧k−1+γkΦLYL−ΦTLGk−1,0<γk〈2λmaxΦLΦTL.The T-GD algorithm does not need to compute the inverse of the information matrixΦLΦTL, but it requires calculating the eigenvalues of the information matrix to choose a suitable step size to keep the T-GD algorithm convergent. When ΦLΦTL has a high order, computing its eigenvalues is also a challenging problem.


## 3. Two Modified GD Algorithms

In this section, two modified GD algorithms are developed which aim to avoid eigenvalue calculation and to increase the convergence rate.


### 3.1. Exhaustive Search-Based GD Algorithm

The PSO algorithm is an intelligent search algorithm, which assigns plenty of particles (initial parameter estimates) first, and then computes the personal best estimates and the global best estimates in each iteration [30, 31]. If the number of the particles is larger, the estimates can easily achieve the true values. Inspired by the PSO algorithm, an exhaustive search-based GD algorithm is developed in this subsection. Its basic idea is to assign several step sizes for a negative direction in each iteration, and the smallest cost function has the best step size.Assume that the parameter estimate in the iterationk−1 is G∧k−1, the parameter estimate in the iteration k is computed by(11)G∧k=G∧k−1+γΦLYL−ΦTLGk−1.If we assign a random step size for the above GD algorithm, we can find that (1) a small step size will have a slow convergence rate; (2) a large step size may lead to divergence of the GD algorithm. To choose a suitable step size and to avoid the eigenvalue calculation, we assignS step sizes for the GD algorithm in each iteration.Define an interval in an iterationk as(12)Tk=0,Sk.ChooseS uniformly distributed terms between 0,Sk, that is(13)γk1=SkS,⋯,γki=iSkS,⋯,γkS=Sk.Based on theS step sizes, we have S parameter estimates, that is(14)G∧k1=G∧k−1+γk1ΦLYL−ΦTLGk−1⋮G∧ks=G∧k−1+γkSΦLYL−ΦTLGk−1.Among theS parameter estimates, we will choose the best one. Once the S parameter estimates in iteration k have been obtained, the S corresponding cost functions are computed by(15)JG∧k1=12YL−ΦTLG∧k12⋮JG∧kS=12YL−ΦTLG^kS2.Let(16)G∧kbest=argminG∧kiJG∧k1,…,JG∧kS.That is, the smallest cost function has the best parameter estimate in iterationk.Then, the steps of the exhaustive search-based GD (ES-GD) algorithm are listed as follows:Algorithm 1: ES-GD algorithm.
InitialiseG∧0=1/p0, p0=106,1 is a vector whose entries all equal to 1Collect measurable datau1,…,uL and y1,…,yLAssign the value forSrepeatfork=1,2,⋯, doAssign0,Sk,Chooseγk1,…,γkSUpdateG∧ki,i=1,…,SComputeJG∧ki,i=1,…,SCompareJG∧ki,i=1,…,S and choose G∧kbestLetG∧k=G∧kbestenduntil convergenceRemark 2.
The same as the PSO algorithm, a largerS can lead to a more accurate parameter estimate G∧k. However, two problems exist for a poor Sk: (1) if Sk is small, all the step sizes can make JG∧ki〈JG∧k−1,i=1,2,…,S, in this case, the step.size is quite small, we then assign Sknew=2Skold ; (2) if Sk is too large, all the step sizes lead to JG∧ki〉JG∧k−1,i=1,2,…,S , in this case, we should assign Sknew=1/2Skold to keep the ES-GD algorithm convergent.Remark 3.
The ES-GD algorithm uses the exhaustive search method to choose the step size; the “best” step size in each iteration is better than the step size which is randomly chosen. However, we have no confidence in the “best” step size because it is not the best one. In addition, a largerS can make the “best” step size closer to the true one, but a larger S also leads to heavier computational efforts.


### 3.2. Power-Based GD Algorithm

In [20], the authors have given the best step size for the cost function(17)JG=12YL−ΦTLG2,is(18)γbest=2λmaxΦLΦTL+λminΦLΦTL.Therefore, to get the actual best step sizeγbest, one should compute both the maximum and minimum eigenvalues of the information matrix ΦLΦTL.Since the eigenvalues of a high-order matrix are difficult to compute, next, we introduce the power method. The power method can get the maximum eigenvalue of a matrix using an iterative method.For simplicity, let(19)Φ=ΦLΦTL.Assign an initial non-zero vectorx0, and use the following iterative function to get a sequence xk,(20)xk=Φxk−1.Let(21)xk=xk1,…,xkMT.The following lemma is obtained.Lemma 1.
For a symmetric positive definite (SPD) matrixΦ, the sequence xk is computed by (4). Then, the maximum eigenvalue of Φ is computed by(22)λmaxΦ=xkixk−1i,i=1,or,2,3,…,M.Proof.
SinceΦ is SPD, it has M eigenvalues λ1,…,λM, and their corresponding eigenvectors are d1,…,d2. Let(23)di=1,i=1,2,…,M,and d1,d2,…,dM are linearly independent. There exist M constants α1,α2,…,αM which are not all equal to zero, and the initial vector x0 can be written by(24)x0=∑i=1Mαidi.
Without loss of generality, assume that the eigenvalues ofΦ satisfy(25)λ1〉λ2〉⋯〉λM〉0.
Based on (20), it gives rise to(26)xk=Φxk−1=Φkx0=∑i=1MαiΦkdi=∑i=1Mαiλikdi.
It follows that(27)xk=λ1kα1d1+∑i=2Mλiλ1kαidi.
Sinceλ1〉λi,i=2,3,…,M, when k⟶∞, we have(28)limk⟶∞λiλ1k=0.
Then, (27) can be rewritten by(29)xk=λ1kα1d1.
Therefore, we can get that(30)λ1=xkixki−1,k⟶∞.
The proof is completed.
The power method can only get the maximum eigenvalue ofΦ. However, to get the best step size, one also should compute the minimum eigenvalue of Φ. Next, we introduce an effective method to compute the minimum eigenvalue of Φ.
Once the maximum eigenvalue ofΦ is obtained, we then assign a new term λ¯ as follows:(31)λ¯=λ1+δ,δ〉0,where δ is a positive constant which is chosen on a case by case basis.
Define a new matrixΨ(32)Ψ=λ¯I−Φ.
Then, the following lemma can be obtained.Lemma 2.
For a symmetric positive definite matrixΦ , its eigenvalues are λ1〉λ2〉⋯〉λM〉0. A matrix Ψ is defined by (6). Then, the eigenvalues of the matrix Ψ are(33)λiΨ=λ¯−λi,i=1,2,…,M,0〈λ1Ψ〈λ2Ψ〈⋯〈λMΨ.Proof.
For an SPD matrixΦ, there exists a nonsingular matrix Q which can guarantee(34)Φ=Q−1diagλ1,λ2,…,λMQ.
Then, the matrixΨ is written by(35)Ψ=Q−1diagλ¯,λ¯,…,λ¯Q−Q−1diagλ1,λ2,…,λMQ=Q−1diagλ¯−λ1,λ¯−λ2,…,λ¯−λMQ=Q−1diagλ1Ψ,λ2Ψ,…,λMΨQ.
Sinceλ¯=λ1+δ and δ〉0, we have(36)0〈λ1Ψ〈λ2Ψ〈⋯〈λMΨ.
For the SPD matrixΨ, using the power method can obtain its maximum eigenvalue λMΨ, and then the minimum eigenvalue λM of the matrix Φ can be computed by(37)λM=λ¯−λMΨ.
When the maximum and minimum eigenvalues of the matrixΦ are obtained, we can get the best step size.
The steps of the power-based GD (P-GD) algorithm are listed as follows:Algorithm 2: P-GD algorithm.
InitialiseG∧0=1/p0 , p0=106,1 is a vector whose entries all equal to 1Collect measurable datau1,…,uL and y1,…,yLUse the power method to computeλ1Assign a positive constantδ based on λ1Construct an SPD matrixΨUse the power method to compute the maximum eigenvalueλMΨ of ΨCalculateλM based on λMΨCompute the best step sizeγbestrepeatfork=1,2,⋯, doUpdateG∧kenduntil convergenceRemark 4.
If the maximum eigenvalueλ1 is not much bigger than λ2, to compute the maximum eigenvalue λ1 is time-consuming. Because the value of λ2/λ1k will take more iterations to converge to zero.Remark 5.
The choice of the positive constantδ is very important; to get the maximum λMΨ quickly, we would do better to choose a small δ. On the other hand, due to the estimation error, a small δ may lead Ψ not be an SPD matrix.Remark 6.
Recently, a novel GD algorithm, termed as fractional stochastic GD algorithm, has been proposed for parameter estimation. This algorithm is a well complement to the traditional GD algorithm, which can be widely used for different kinds of models [32–34].


## 3.1. Exhaustive Search-Based GD Algorithm

The PSO algorithm is an intelligent search algorithm, which assigns plenty of particles (initial parameter estimates) first, and then computes the personal best estimates and the global best estimates in each iteration [30, 31]. If the number of the particles is larger, the estimates can easily achieve the true values. Inspired by the PSO algorithm, an exhaustive search-based GD algorithm is developed in this subsection. Its basic idea is to assign several step sizes for a negative direction in each iteration, and the smallest cost function has the best step size.Assume that the parameter estimate in the iterationk−1 is G∧k−1, the parameter estimate in the iteration k is computed by(11)G∧k=G∧k−1+γΦLYL−ΦTLGk−1.If we assign a random step size for the above GD algorithm, we can find that (1) a small step size will have a slow convergence rate; (2) a large step size may lead to divergence of the GD algorithm. To choose a suitable step size and to avoid the eigenvalue calculation, we assignS step sizes for the GD algorithm in each iteration.Define an interval in an iterationk as(12)Tk=0,Sk.ChooseS uniformly distributed terms between 0,Sk, that is(13)γk1=SkS,⋯,γki=iSkS,⋯,γkS=Sk.Based on theS step sizes, we have S parameter estimates, that is(14)G∧k1=G∧k−1+γk1ΦLYL−ΦTLGk−1⋮G∧ks=G∧k−1+γkSΦLYL−ΦTLGk−1.Among theS parameter estimates, we will choose the best one. Once the S parameter estimates in iteration k have been obtained, the S corresponding cost functions are computed by(15)JG∧k1=12YL−ΦTLG∧k12⋮JG∧kS=12YL−ΦTLG^kS2.Let(16)G∧kbest=argminG∧kiJG∧k1,…,JG∧kS.That is, the smallest cost function has the best parameter estimate in iterationk.Then, the steps of the exhaustive search-based GD (ES-GD) algorithm are listed as follows:Algorithm 1: ES-GD algorithm.
InitialiseG∧0=1/p0, p0=106,1 is a vector whose entries all equal to 1Collect measurable datau1,…,uL and y1,…,yLAssign the value forSrepeatfork=1,2,⋯, doAssign0,Sk,Chooseγk1,…,γkSUpdateG∧ki,i=1,…,SComputeJG∧ki,i=1,…,SCompareJG∧ki,i=1,…,S and choose G∧kbestLetG∧k=G∧kbestenduntil convergenceRemark 2.
The same as the PSO algorithm, a largerS can lead to a more accurate parameter estimate G∧k. However, two problems exist for a poor Sk: (1) if Sk is small, all the step sizes can make JG∧ki〈JG∧k−1,i=1,2,…,S, in this case, the step.size is quite small, we then assign Sknew=2Skold ; (2) if Sk is too large, all the step sizes lead to JG∧ki〉JG∧k−1,i=1,2,…,S , in this case, we should assign Sknew=1/2Skold to keep the ES-GD algorithm convergent.Remark 3.
The ES-GD algorithm uses the exhaustive search method to choose the step size; the “best” step size in each iteration is better than the step size which is randomly chosen. However, we have no confidence in the “best” step size because it is not the best one. In addition, a largerS can make the “best” step size closer to the true one, but a larger S also leads to heavier computational efforts.


## 3.2. Power-Based GD Algorithm

In [20], the authors have given the best step size for the cost function(17)JG=12YL−ΦTLG2,is(18)γbest=2λmaxΦLΦTL+λminΦLΦTL.Therefore, to get the actual best step sizeγbest, one should compute both the maximum and minimum eigenvalues of the information matrix ΦLΦTL.Since the eigenvalues of a high-order matrix are difficult to compute, next, we introduce the power method. The power method can get the maximum eigenvalue of a matrix using an iterative method.For simplicity, let(19)Φ=ΦLΦTL.Assign an initial non-zero vectorx0, and use the following iterative function to get a sequence xk,(20)xk=Φxk−1.Let(21)xk=xk1,…,xkMT.The following lemma is obtained.Lemma 1.
For a symmetric positive definite (SPD) matrixΦ, the sequence xk is computed by (4). Then, the maximum eigenvalue of Φ is computed by(22)λmaxΦ=xkixk−1i,i=1,or,2,3,…,M.Proof.
SinceΦ is SPD, it has M eigenvalues λ1,…,λM, and their corresponding eigenvectors are d1,…,d2. Let(23)di=1,i=1,2,…,M,and d1,d2,…,dM are linearly independent. There exist M constants α1,α2,…,αM which are not all equal to zero, and the initial vector x0 can be written by(24)x0=∑i=1Mαidi.
Without loss of generality, assume that the eigenvalues ofΦ satisfy(25)λ1〉λ2〉⋯〉λM〉0.
Based on (20), it gives rise to(26)xk=Φxk−1=Φkx0=∑i=1MαiΦkdi=∑i=1Mαiλikdi.
It follows that(27)xk=λ1kα1d1+∑i=2Mλiλ1kαidi.
Sinceλ1〉λi,i=2,3,…,M, when k⟶∞, we have(28)limk⟶∞λiλ1k=0.
Then, (27) can be rewritten by(29)xk=λ1kα1d1.
Therefore, we can get that(30)λ1=xkixki−1,k⟶∞.
The proof is completed.
The power method can only get the maximum eigenvalue ofΦ. However, to get the best step size, one also should compute the minimum eigenvalue of Φ. Next, we introduce an effective method to compute the minimum eigenvalue of Φ.
Once the maximum eigenvalue ofΦ is obtained, we then assign a new term λ¯ as follows:(31)λ¯=λ1+δ,δ〉0,where δ is a positive constant which is chosen on a case by case basis.
Define a new matrixΨ(32)Ψ=λ¯I−Φ.
Then, the following lemma can be obtained.Lemma 2.
For a symmetric positive definite matrixΦ , its eigenvalues are λ1〉λ2〉⋯〉λM〉0. A matrix Ψ is defined by (6). Then, the eigenvalues of the matrix Ψ are(33)λiΨ=λ¯−λi,i=1,2,…,M,0〈λ1Ψ〈λ2Ψ〈⋯〈λMΨ.Proof.
For an SPD matrixΦ, there exists a nonsingular matrix Q which can guarantee(34)Φ=Q−1diagλ1,λ2,…,λMQ.
Then, the matrixΨ is written by(35)Ψ=Q−1diagλ¯,λ¯,…,λ¯Q−Q−1diagλ1,λ2,…,λMQ=Q−1diagλ¯−λ1,λ¯−λ2,…,λ¯−λMQ=Q−1diagλ1Ψ,λ2Ψ,…,λMΨQ.
Sinceλ¯=λ1+δ and δ〉0, we have(36)0〈λ1Ψ〈λ2Ψ〈⋯〈λMΨ.
For the SPD matrixΨ, using the power method can obtain its maximum eigenvalue λMΨ, and then the minimum eigenvalue λM of the matrix Φ can be computed by(37)λM=λ¯−λMΨ.
When the maximum and minimum eigenvalues of the matrixΦ are obtained, we can get the best step size.
The steps of the power-based GD (P-GD) algorithm are listed as follows:Algorithm 2: P-GD algorithm.
InitialiseG∧0=1/p0 , p0=106,1 is a vector whose entries all equal to 1Collect measurable datau1,…,uL and y1,…,yLUse the power method to computeλ1Assign a positive constantδ based on λ1Construct an SPD matrixΨUse the power method to compute the maximum eigenvalueλMΨ of ΨCalculateλM based on λMΨCompute the best step sizeγbestrepeatfork=1,2,⋯, doUpdateG∧kenduntil convergenceRemark 4.
If the maximum eigenvalueλ1 is not much bigger than λ2, to compute the maximum eigenvalue λ1 is time-consuming. Because the value of λ2/λ1k will take more iterations to converge to zero.Remark 5.
The choice of the positive constantδ is very important; to get the maximum λMΨ quickly, we would do better to choose a small δ. On the other hand, due to the estimation error, a small δ may lead Ψ not be an SPD matrix.Remark 6.
Recently, a novel GD algorithm, termed as fractional stochastic GD algorithm, has been proposed for parameter estimation. This algorithm is a well complement to the traditional GD algorithm, which can be widely used for different kinds of models [32–34].


## 4. Convergence Properties of the Two Modified GD Algorithms

The convergence properties of the two modified GD algorithms are given in the following which offer theory guidance for researchers.


### 4.1. Convergence Analysis of the ES-GD Algorithm

Rewrite the ES-GD algorithm as follows:(38)G∧k=G∧k−1+γΦLYL−ΦTLGk−1.SubtractingG on both sides of the above equation yields(39)E∧k=E∧k−1−γΦLΦTLEk−1+γΦLVL=I−γΦLΦTLEk−1+γΦLVL,where E∧k=G∧k−G. Since VL is a Gaussian white noise and is independent on ΦL, the above equation is simplified as(40)E∧k=I−γΦLΦTLEk−1.Based on the exhaustive search-based, in each iteration, we will find an optimalγ which guarantees(41)I−γΦLΦTL≤1.It gives rise to(42)E∧k≤Ek−1.Therefore, the ES-GD algorithm is convergent.


### 4.2. Convergence Analysis of the P-GD Algorithm

The P-GD algorithm is written by(43)G∧k=G∧k−1+2λmaxΦLΦTL+λminΦLΦTL⋅ΦLYL−ΦTLGk−1.Subtracting the true valueG on both sides of the above equation yields(44)E∧k=E∧k−1−2λmaxΦLΦTL+λminΦLΦTL⋅ΦLΦTLEk−1=I−2ΦLΦTLλmaxΦLΦTL+λminΦLΦTLE∧k−1.For simplicity, let(45)λ1=λmaxΦLΦTL,λM=λminΦLΦTL.Equation (44) is simplified as(46)E∧k=I−2ΦLΦTLλ1+λME∧k−1.For an SPD matrixΦLΦTL, there exists a matrix Q which can ensure(47)ΦLΦTL=Q−1diagλ1,λ2,…,λMQ.It follows that equation (46) can be transformed into(48)E∧k=Q−1Q−2λ1+λMQ−1diagλ1,λ2,…,λMQE∧k−1=Q−1diag1−λ1λ1+λM,1−λ2λ1+λM,…,1−λMλ1+λMQE∧k−1=Q−1diagλM−λ1λ1+λM,…,λ1−λMλ1+λMQE∧k−1.Clearly, all the absolute values in the diagonal matrix are smaller than 1, then we have(49)E∧k≤E∧k−1.Therefore, the P-GD algorithm is convergent.Remark 7.
In the P-GD algorithm, the maximum absolute value in the diagonal matrix isλ1−λM/λ1+λM, that is,(50)E∧k=λ1−λMλ1+λME∧k−1=τ−1τ+1E∧k−1,where τ=λ1/λM is the conditioned number of the matrix ΦLΦTL. If the matrix ΦLΦTL is ill-conditioned, no matter what the step size is, the convergence rates are always very slow. In this case, we can try to reconstruct a new information matrix ΦLΦTL or use the fractional stochastic GD algorithm proposed in [32–34] to increase the convergence rates.


## 4.1. Convergence Analysis of the ES-GD Algorithm

Rewrite the ES-GD algorithm as follows:(38)G∧k=G∧k−1+γΦLYL−ΦTLGk−1.SubtractingG on both sides of the above equation yields(39)E∧k=E∧k−1−γΦLΦTLEk−1+γΦLVL=I−γΦLΦTLEk−1+γΦLVL,where E∧k=G∧k−G. Since VL is a Gaussian white noise and is independent on ΦL, the above equation is simplified as(40)E∧k=I−γΦLΦTLEk−1.Based on the exhaustive search-based, in each iteration, we will find an optimalγ which guarantees(41)I−γΦLΦTL≤1.It gives rise to(42)E∧k≤Ek−1.Therefore, the ES-GD algorithm is convergent.


## 4.2. Convergence Analysis of the P-GD Algorithm

The P-GD algorithm is written by(43)G∧k=G∧k−1+2λmaxΦLΦTL+λminΦLΦTL⋅ΦLYL−ΦTLGk−1.Subtracting the true valueG on both sides of the above equation yields(44)E∧k=E∧k−1−2λmaxΦLΦTL+λminΦLΦTL⋅ΦLΦTLEk−1=I−2ΦLΦTLλmaxΦLΦTL+λminΦLΦTLE∧k−1.For simplicity, let(45)λ1=λmaxΦLΦTL,λM=λminΦLΦTL.Equation (44) is simplified as(46)E∧k=I−2ΦLΦTLλ1+λME∧k−1.For an SPD matrixΦLΦTL, there exists a matrix Q which can ensure(47)ΦLΦTL=Q−1diagλ1,λ2,…,λMQ.It follows that equation (46) can be transformed into(48)E∧k=Q−1Q−2λ1+λMQ−1diagλ1,λ2,…,λMQE∧k−1=Q−1diag1−λ1λ1+λM,1−λ2λ1+λM,…,1−λMλ1+λMQE∧k−1=Q−1diagλM−λ1λ1+λM,…,λ1−λMλ1+λMQE∧k−1.Clearly, all the absolute values in the diagonal matrix are smaller than 1, then we have(49)E∧k≤E∧k−1.Therefore, the P-GD algorithm is convergent.Remark 7.
In the P-GD algorithm, the maximum absolute value in the diagonal matrix isλ1−λM/λ1+λM, that is,(50)E∧k=λ1−λMλ1+λME∧k−1=τ−1τ+1E∧k−1,where τ=λ1/λM is the conditioned number of the matrix ΦLΦTL. If the matrix ΦLΦTL is ill-conditioned, no matter what the step size is, the convergence rates are always very slow. In this case, we can try to reconstruct a new information matrix ΦLΦTL or use the fractional stochastic GD algorithm proposed in [32–34] to increase the convergence rates.


## 5. Examples

Example 1.
Consider the following time-delayed model,(51)yt=g1ut−τ−1+g2ut−τ−2+g3ut−τ−3+g4ut−τ−4+g5ut−τ−5+vt=0.5ut−τ−1+0.67ut−τ−2−0.34ut−τ−3+0.23ut−τ−4+0.76ut−τ−5+vt.
Assume that the time delay isτ=2 and assign M=5, we have(52)yt=g0,1ut−1+g0,2ut−2+g1ut−3+g2ut−4+g3ut−5+g4ut−6+g5ut−7+g1,1ut−8+g1,2ut−9+g1,3ut−10+vt,G=0,0,0.5,0.67,−0.34,0.23,0.76,0,0,0T.
In simulation, we collect 500 sets of input and output data, whereut∼N0,1 and vt∼N0,0.12. Use the T-GD, ES-GD, and P-GD algorithms for the time-delayed model. The parameter estimates and their estimation errors ς=G−Gk/G are shown in Figure 1 and Table 1. The elapsed times of these three algorithms are illustrated in Table 2: the second row means that all the three algorithms run the same iteration, and the third row shows that the three algorithms have almost the same estimation error.
Assign a thresholdϱ=0.005. Compare the estimates (20-th iteration) with the threshold, if the absolute value of the estimate is smaller than the threshold, then it will be assigned as zero. We can get that the time-delay is 2.
In addition, we use the power method to compute the maximum and minimum eigenvalues of the information matrixΦLΦTL, and the estimates are shown in Figure 2.
From this simulation, we can obtain the following conclusions:(1)
The P-GD algorithm has the fastest convergence rates, then is the ES-GD algorithm, and the T-GD algorithm has the slowest convergence rates; this can be shown in Figure1;(2)
All the three algorithms can obtain the parameter estimates and the time-delay estimates simultaneously, as shown in Table1;(3)
The P-GD algorithm is the most effective algorithm, follows the ES-GD algorithm, and the last one is the T-GD algorithm; this can be shown in Table2;(4)
The power method can obtain the maximum eigenvalue and the minimum eigenvalue of the information matrix, as shown in Figure2.Figure 1
The parameter estimation errorsς versus k.Table 1
The parameter estimates and their estimation errors.
Algorithmskg0,1g0,2g1g2g3g4g5g1,1g1,2g1,3ς%T-GD1−0.000760.000140.049230.06527−0.033050.022490.074350.00028−0.000610.0003990.224592−0.001390.000270.093590.12417−0.062880.042800.141420.00050−0.001110.0007081.405465−0.002680.000620.202090.26859−0.136120.092590.305800.00092−0.002050.0013459.7945110−0.003510.001040.322270.42938−0.217830.148060.488580.00109−0.002480.0017235.7631815−0.003540.001320.393730.52564−0.266880.181290.597840.00096−0.002270.0017121.3989120−0.003280.001500.436220.58328−0.296320.201210.663160.00073−0.001860.0015712.81303ES-GD10.01755−0.029120.316940.47762−0.266690.158670.52597−0.043810.011150.0290831.1041220.01217−0.018850.432940.61302−0.326750.207610.68377−0.026560.008530.0185910.4764650.00153−0.003510.496230.66886−0.344690.229320.75506−0.001260.00332−0.000350.78800100.00023−0.002300.499480.67068−0.345180.230050.758120.001010.00266−0.002510.59350150.00021−0.002290.499510.67069−0.345190.230050.758150.001030.00265−0.002540.59388200.00021−0.002290.499510.67069−0.345190.230050.758150.001030.00265−0.002540.59389P-GD10.02557−0.042430.461840.69598−0.388620.231210.76643−0.063840.016250.042379.5312220.00246−0.001230.496990.66530−0.338470.229430.751100.001980.003270.000740.9687550.00022−0.002290.499500.67070−0.345200.230050.758140.001020.00266−0.002530.59415100.00021−0.002290.499510.67069−0.345190.230050.758150.001030.00265−0.002540.59389150.00021−0.002290.499510.67069−0.345190.230050.758150.001030.00265−0.002540.59389200.00021−0.002290.499510.67069−0.345190.230050.758150.001030.00265−0.002540.59389True values0.000000.000000.500000.67000−0.340000.230000.760000.000000.000000.000000.00000Table 2
The elapsed times of the three algorithms.
AlgorithmT-GDES-GDP-GDElapsed time (second)0.2482 (k = 20)0.3125 (k = 20)0.2876 (k = 20)Elapsed time (second)0.3516 (k = 40)0.1325 (k = 6)0.1008 (k = 4)Figure 2
The maximum and minimum eigenvalue estimates versusk.Example 2.
A water tank system with a communication channel is proposed for simulation, see Figure3, where ut is the position of the inlet water valve, and yt is the water level of Tank 2 and sampled by a pressure sensor. There exists a time-delay τ=3. The water tank system is modeled by the following model [35]:(53)yt=0.3yt−1−3−0.1yt−2−3+0.15yt−3−3−0.3ut−1+0.2ut−2+0.13ut−3+vt.
Using the T-GD, ES-GD, and P-GD algorithms for this model, the parameter estimates and their estimation errors are shown in Figure4.
This example also shows that the two modified GD algorithms have faster convergence rates than those of the T-GD algorithm.Figure 3
A water tank system.Figure 4
The parameter estimation errorsς versus k.


## 6. Conclusions

Two modified GD algorithms are proposed for systems with time-delay in this paper. The first is the ES-GD algorithm that does not require the eigenvalue calculation. The second is the P-GD algorithm, that can get the best step size by using the power method. These two modified GD algorithms have faster convergence rates than those of the T-GD algorithm, and can obtain the parameter estimates and time-delay estimates simultaneously. Thus, they can be widely used in engineering practices.In this paper, we only use the modified GD algorithms for the time-delayed systems. If the systems have other kinds of hidden variables, e.g., missing outputs and model identities, can these two algorithms be also effective? this topic will remain as an open issue in future.



---

*Source: 9244890-2022-09-08.xml*