# Nonlinear Blind Identification with Three-Dimensional Tensor Analysis

Authors: I. Cherif; S. Abid; F. Fnaiech
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper deals with the analysis of a third-order tensor composed of a fourth-order output cumulants used for blind identification of a second-order Volterra-Hammerstein series. It is demonstrated that this nonlinear identification problem can be converted in a multivariable system with multiequations having the form ofAx+By=c. The system may be solved using several methods. Simulation results with the Iterative Alternating Least Squares (IALS) algorithm provide good performances for different signal-to-noise ratio (SNR) levels. Convergence issues using the reversibility analysis of matrices A and B are addressed. Comparison results with other existing algorithms are carried out to show the efficiency of the proposed algorithm.


---

## Body

## 1. Introduction

Nonlinear system modeling based on real-world input/output measurements is so far used in many applications. The appropriate model and the determination of corresponding parameters using the input/output data are owned to apply a suitable and efficient identification method [1–7].Hammerstein models are special classes of second-order Volterra systems where the second-order homogenous Volterra kernel is diagonal [8]. These systems have been successfully used to model nonlinear systems in a number of practical applications in several areas such as chemical process, biological process, signal processing, and communications [9–12], where, for example, in digital communication systems, the communication channels are usually impaired by a nonlinear intersymbol interference (ISI). Channel identification allows compensating the ISI effects at the receivers.In [13], a penalty transformation method is developed. Indeed a penalty function is formed by equations relating the unknown parameters of the model with the autocorrelations of the signal. This function is then included in the cost function yielding to an augmented Lagrangian function. It has been demonstrated that this approach gives good identification results for a nonlinear systems. However, this approach is still sensitive to additive Gaussian noise because the 2nd-order moment is used as a constraint. Authors, in [7], overcame this sensitivity by using 4th-order cumulants as a constraint instead of 2nd-order moments in order to smooth out the additive Gaussian noise. But the proposed approach which is based on a simplex-genetic algorithm becomes so long and computationally complex.The main drawback of identification with Volterra series lies on the parametric complexity and the need to estimate a very big number of parameters. In many cases, Volterra series identification problem may be well simplified using the tensor formulation [10–12, 14].Authors, in [10], used a parallel factor (PARAFAC) decomposition of the kernels to derive Volterra-PARAFAC models yielding an important parametric complexity reduction for Volterra kernels of order higher than two. They proved that these models are equivalent to a set of parallel Wiener models. Consequently, they proposed three adaptive algorithms for identifying these proposed Volterra-PARAFAC models for complex-valued input/output signals, namely, the extended complex Kalman filter, the complex least mean square (CLMS) algorithm, and the normalized CLMS algorithm.In this paper, the algorithm derived in [14] is extended to be applied to blind identification of a general second-order Volterra-Hammerstein system. The main idea is to develop a general expression for each direction slices of a cubic tensor and then express the tensor slices in an unfolded representation. The three-dimensional tensor elements are formed by the fourth-order output cumulants. This yields to an Iterative Alternating Least Square (IALS) algorithm which has the benefit over the original Volterra filters in terms of implementation and complexity reduction. A convergence analysis based on matrices reversibility study is given showing that the proposed IALS algorithm converges to optimal solutions in the least mean squares sense. Furthermore, some simulation results and comparisons with different existing algorithms are provided.The present work is organized as follows; in Section2, a brief study of the three-dimensional tensor is presented. In Section 3, the model under study and the related output cumulants are then proposed, whereas, in Section 4 the decomposition analysis of the cumulant tensor is developed. In Sections 5 to 8, we give, respectively, the proposed blind identification algorithm, the convergence study, some simulation results, and at the end some main conclusions are drawn.


## 2. Three-Dimensional Tensor and Different Slice Expressions

A three-dimensional tensorC∈ℂM×M×M can be expressed by(2.1)C=∑i=1M∑j=1M∑k=1MCijkei(M)∘ej(M)∘ek(M),
where Cijk is the tensor value in the position (i,j,k) of the cube with dimension M, ep(M) denotes the pth canonical basis vector with dimension M, and the symbol ∘ stands for the outer product (Figure 1).Figure 1
Cubic tensor representation.A cubic tensor can be always sliced along three possible directions (horizontal, vertical, and frontal) as depicted in Figure2. This yields, in each case, to M matrices of M×M dimensions.Different direction slices of a cubic tensor.
(a)
horizontal slices(b)
vertical slices(c)
frontal slicesThe expression of theith slice in the horizontal direction is given by(2.2)Ci••=∑j=1M∑k=1MCijkej(M)∘ek(M)=∑j=1M∑k=1MCijkej(M)ek(M)T.In the same manner, the other matrix expressions along with the vertical and frontal directions are expressed, respectively, by(2.3)C•j•=∑i=1M∑k=1MCijkei(M)ek(M)T;C••k=∑i=1M∑j=1MCijkei(M)ej(M)T.It is important to express the tensor slices in an unfolded representation, obtained by stacking up the 2D matrices. Hence, three unfolded representations ofC are obtained. For the horizontal, the vertical, and the frontal directions, we get, respectively,(2.4)C[1]=(C1..C2..⋮CM..);C[2]=(C.1.C.2.⋮C.M.);C[3]=(C..1C..2⋮C..M).We note that each matrixC[p]:p=1,2,3 is an (M×M,M) one.


## 3. Nonlinear System Model and Output Cumulants Analysis

We focus on the identification of a second-order Volterra-Hammerstein model with finite memory as it is given in [14]:(3.1)y(n)=∑k=0Mh1(k)u(n-k)+∑k=0Mh2(k)u2(n-k);h1(0)=1;h2(0)=1,hi(M)≠0,
where u(n) is the input of the system, assumed to be a stationary zero mean Gaussian white random process with E(u2(n))=γ2. M stands for the model order.The Hammerstein coefficients vectorsh1 and h2 are defined by(3.2)hp=[hp(0),hp(1),…,hp(M)]T;p=1;2.As we evoked in the Introduction, identification algorithms based on the computation of 2nd-order output cumulants are sensitive to additive Gaussian noise because 2nd-order cumulants of this latter are in general different to zero. Since the 4th-order cumulants of additive Gaussian noise is null, it will be interesting to use the 4th-order output cumulants to derive identification algorithms. But this will introduce another problem which is the computation complexity. In this paper, we will overcome this shortcoming by using a tensor analysis.To determine the kernels of this model, we will generate the fourth-order output cumulants. For this purpose, we need to use the standard properties of cumulants and the Leonov-Shiryaev formula for manipulating products of random variables.The fourth-order output cumulant is given by [15]:(3.3)c4y(i1,i2,i3)=cum⁡[y(n+i1),y(n+i2),y(n+i3),y(n)]=8γ23{ϕ0(i1,i2,i3,0)+ϕ0(i1,i3,i2,0)+ϕ0(i2,i3,i1,0)+ϕ0(i1,0,i2,i3)+ϕ0(i2,0,i1,i3)+ϕ0(i3,0,i1,i2)}+48γ24∑l=0Mh2(l+i1)h2(l+i2)h2(l+i3)h2(l),
where(3.4)ϕ0(i1,i2,i3,i4)=∑l=0Mh1(l+i1)h1(l+i2)h2(l+i3)h2(l+i4).
It is easy to verify that c4y(i1,i2,i3)=0 for all |i1|,|i2|,|i3|>M.All the nonzero terms ofc4y(i1,i2,i3) are obtained for (i1,i2,i3)∈[-M,M]3. Such a choice allows us to construct a maximal redundant information, in which the fourth-order cumulants are taken for time lags i1,i2, and i3 within the range [-M,M].In the sequel we shall present an analysis of a 3rd-order tensor composed of the 4th-order output cumulants.


## 4. Formulation and Analysis of a Cumulant Cubic Tensor

Let us define the three-dimensional tensorC(4,y)∈ℂ(2M+1)×(2M+1)×(2M+1), in which the element in position (i,j,k) corresponds to c4y(i1,i2,i3), with i=i1+M+1; j=i2+M+1; k=i3+M+1. As i1,i2,i3∈[-M,M], we get i,j,k∈[1,2M+1]. Thus,
(4.1)Cijk=c4y(i-M-1,j-M-1,k-M-1)=8γ23{ϕ0(i-M-1,j-M-1,k-M-1,0)+ϕ0(i-M-1,k-M-1,j-M-1,0)+ϕ0(j-M-1,k-M-1,i-M-1,0)+ϕ0(i-M-1,0,j-M-1,k-M-1)+ϕ0(j-M-1,0,i-M-1,k-M-1)+ϕ0(k-M-1,0,i-M-1,j-M-1)}+48γ24∑l=0Mh2(l+i-M-1)h2(l+j-M-1)×h2(l+k-M-1)h2(l),ϕ0(·,·,·,·) is given by (3.4). It follows that(4.2)Cijk=8γ23{∑l=0Mh1(l+i-M-1)h1(l+j-M-1)h2(l+k-M-1)h2(l)+∑l=0Mh1(l+i-M-1)h1(l+k-M-1)h2(l+j-M-1)h2(l)+∑l=0Mh1(l+j-M-1)h1(l+k-M-1)h2(l+i-M-1)h2(l)+∑l=0Mh1(l+i-M-1)h1(l)h2(l+j-M-1)h2(l+k-M-1)+∑l=0Mh1(l+j-M-1)h1(l)h2(l+i-M-1)h2(l+k-M-1)+∑l=0Mh1(l+k-M-1)h1(l)h2(l+i-M-1)h2(l+j-M-1)}+48γ24∑l=0Mh2(l+i-M-1)h2(l+j-M-1)h2(l+k-M-1)h2(l).Then, expression of the tensorC will be given by
(4.3)C=8γ23∑i=12M+1∑j=12M+1∑k=12M+1{∑l=0Mh1(l+i-M-1)h1(l+j-M-1)h2(l+k-M-1)h2(l)+∑l=0Mh1(l+i-M-1)h1(l+k-M-1)h2(l+j-M-1)h2(l)+∑l=0Mh1(l+j-M-1)h1(l+k-M-1)h2(l+i-M-1)h2(l)+∑l=0Mh1(l+i-M-1)h1(l)h2(l+j-M-1)h2(l+k-M-1)+∑l=0Mh1(l+j-M-1)h1(l)h2(l+i-M-1)h2(l+k-M-1)+∑l=0Mh1(l+k-M-1)h1(l)h2(l+i-M-1)h2(l+j-M-1)}×ei(2M+1)∘ej(2M+1)∘ek(2M+1)+48γ24×∑i=12M+1∑j=12M+1∑k=12M+1∑l=0Mh2(l+i-M-1)⋅h2(l+j-M-1)h2(l+k-M-1)h2(l)×ei(2M+1)∘ej(2M+1)∘ek(2M+1).The mathematical development of the expression (4.3) yields to(4.4)C=8γ23{∑l=0Mh2(l)h1•l∘h1•l∘h2•l+∑l=0Mh2(l)h1•l∘h2•l∘h1•l+∑l=0Mh2(l)h2•l∘h1•l∘h1•l+∑l=0Mh1(l)h1•l∘h2•l∘h2•l+∑l=0Mh1(l)h2•l∘h1•l∘h2•l+∑l=0Mh1(l)h2•l∘h2•l∘h1•l}+48γ24∑l=0Mh2(l)h2•l∘h2•l∘h2•l,
where(4.5)hp•l=∑m=12M+1hp(l+m-M-1)em(2M+1),p=1;2.This notation leads to define two channel matricesH1;H2∈ℂ(2M+1)×(M+1) as follows:(4.6)Hp≜H(hp)=[hp•0,hp•1,…,hp•M]=(00⋯hp(0)⋮⋮⋱⋮0hp(0)⋯hp(M-1)hp(0)hp(1)⋯hp(M)⋮⋮⋱⋮hp(M-1)hp(M)⋯0hp(M)0⋯0),
with p=1;2, and ℋ(·) is the operator that builds a special Hankel matrix from the vector argument as shown above.Let us compute now the different slices of the proposed tensor.


### 4.1. Horizontal Slices Expressions

From (2.2) and (4.3), we get(4.7)Ci••=8γ23{∑l=0Mh2(l)h1(l+i-M-1)h1•lh2•lT+∑l=0Mh2(l)h1(l+i-M-1)h2•lh1•lT+∑l=0Mh2(l)h2(l+i-M-1)h1•lh1•lT+∑l=0Mh1(l)h1(l+i-M-1)h2•lh2•lT+∑l=0Mh1(l)h2(l+i-M-1)h1•lh2•lT+∑l=0Mh1(l)h2(l+i-M-1)h2•lh1•lT}+48γ24∑l=0Mh2(l)h2(l+i-M-1)h2•lh2•lT,which can be written as(4.8)Ci••=8γ23{H1diag⁡i⁡(H1)Σ2T+Σ2diag⁡i⁡(H1)H1+H1diag⁡i(Σ2)H1T+H2diag⁡i(Σ1)H2T+Σ1diag⁡i(H2)H2T+H2diag⁡i(H2)Σ1T}+48γ24H2diag⁡i(Σ2)H2T,
where Σ1=H1diag(h1); Σ2=H2diag(h2);diagn(·) is the diagonal matrix formed by the nth line of its argument.It can easily be demonstrated that(4.9)Cm••=C•m•=C••m=Cm.It follows that(4.10)Cm=8γ23{H1diag⁡m(H1)Σ2T+Σ2diag⁡m(H1)H1T+H1diag⁡m(Σ2)H1T+H2diag⁡m(Σ1)H2T+Σ1diag⁡m(H2)H2T+H2diag⁡m(H2)Σ1T}+48γ24H2diag⁡m(Σ2)H2T,
with m=1,…,2M+1.The expression of the unfolded tensor representation is given by(4.11)C[U]=(C1C2⋮C2M+1)=8γ23{(H1diag⁡1⁡(H1)Σ2TH1diag⁡2(H1)Σ2T⋮H1diag⁡2M+1(H1)Σ2T)+⋯}︸the 6th terms with respect toCmin(4.10)+48γ24{(H2diag⁡1(Σ2)H2TH2diag⁡2(Σ2)H2T⋮H2diag⁡2M+1(Σ2)H2T)}.To develop this expression, we need the following property.Property 1.
LetA be the matrix with dimensions (M,N) and B the matrix with dimensions (M’,N), then
(4.12)(Adiag⁡1(B)Adiag⁡2(B)⋮Adiag⁡M′(B))=B⋄A,
where ⋄ stands for the Khatri-Rao product.It becomes that(4.13)C[U]=8γ23{(H1⋄H1)diag⁡(h2)H2T+(H1⋄H2)diag⁡(h2)H1T+(H2⋄H1)diag⁡(h2)H1T+(H1⋄H2)diag⁡(h1)H2T+(H2⋄H1)diag⁡(h1)H2T+(H2⋄H2)diag⁡(h1)H1T}+48γ24((H2⋄H2)diag⁡(h2)H2T).


## 4.1. Horizontal Slices Expressions

From (2.2) and (4.3), we get(4.7)Ci••=8γ23{∑l=0Mh2(l)h1(l+i-M-1)h1•lh2•lT+∑l=0Mh2(l)h1(l+i-M-1)h2•lh1•lT+∑l=0Mh2(l)h2(l+i-M-1)h1•lh1•lT+∑l=0Mh1(l)h1(l+i-M-1)h2•lh2•lT+∑l=0Mh1(l)h2(l+i-M-1)h1•lh2•lT+∑l=0Mh1(l)h2(l+i-M-1)h2•lh1•lT}+48γ24∑l=0Mh2(l)h2(l+i-M-1)h2•lh2•lT,which can be written as(4.8)Ci••=8γ23{H1diag⁡i⁡(H1)Σ2T+Σ2diag⁡i⁡(H1)H1+H1diag⁡i(Σ2)H1T+H2diag⁡i(Σ1)H2T+Σ1diag⁡i(H2)H2T+H2diag⁡i(H2)Σ1T}+48γ24H2diag⁡i(Σ2)H2T,
where Σ1=H1diag(h1); Σ2=H2diag(h2);diagn(·) is the diagonal matrix formed by the nth line of its argument.It can easily be demonstrated that(4.9)Cm••=C•m•=C••m=Cm.It follows that(4.10)Cm=8γ23{H1diag⁡m(H1)Σ2T+Σ2diag⁡m(H1)H1T+H1diag⁡m(Σ2)H1T+H2diag⁡m(Σ1)H2T+Σ1diag⁡m(H2)H2T+H2diag⁡m(H2)Σ1T}+48γ24H2diag⁡m(Σ2)H2T,
with m=1,…,2M+1.The expression of the unfolded tensor representation is given by(4.11)C[U]=(C1C2⋮C2M+1)=8γ23{(H1diag⁡1⁡(H1)Σ2TH1diag⁡2(H1)Σ2T⋮H1diag⁡2M+1(H1)Σ2T)+⋯}︸the 6th terms with respect toCmin(4.10)+48γ24{(H2diag⁡1(Σ2)H2TH2diag⁡2(Σ2)H2T⋮H2diag⁡2M+1(Σ2)H2T)}.To develop this expression, we need the following property.Property 1.
LetA be the matrix with dimensions (M,N) and B the matrix with dimensions (M’,N), then
(4.12)(Adiag⁡1(B)Adiag⁡2(B)⋮Adiag⁡M′(B))=B⋄A,
where ⋄ stands for the Khatri-Rao product.It becomes that(4.13)C[U]=8γ23{(H1⋄H1)diag⁡(h2)H2T+(H1⋄H2)diag⁡(h2)H1T+(H2⋄H1)diag⁡(h2)H1T+(H1⋄H2)diag⁡(h1)H2T+(H2⋄H1)diag⁡(h1)H2T+(H2⋄H2)diag⁡(h1)H1T}+48γ24((H2⋄H2)diag⁡(h2)H2T).


## 5. Blind Identification System with Cumulant Tensor

To estimate the Volterra-Hammerstein kernels and to avoid the computation ofHp:p=1,2, we will use the following Khatri-Rao property to propose an Iterative Alternating Least Square (IALS) procedure.Property 2.
If matricesA∈ℂm×n and B∈ℂn×m and vector d∈ℂn are such that X=Adiag(d)B, then it holds that vec(X)=(BT⋄A)d, where vec(·) stands for the vectorizing operator.Applying this property to (4.13), it is straightforward to write(5.1)vec⁡(C[U])=8γ23{(H2⋄H1⋄H1)h2+(H1⋄H1⋄H2)h2+(H1⋄H2⋄H1)h2+(H2⋄H1⋄H2)h1+(H2⋄H2⋄H1)h1+(H1⋄H2⋄H2)h1}+48γ24((H2⋄H2⋄H2)h2)=8γ23{(H2⋄H1⋄H2)+(H2⋄H2⋄H1)+(H1⋄H2⋄H2)}h1+{8γ23((H2⋄H1⋄H1)+(H1⋄H1⋄H2)+(H1⋄H2⋄H1))+48γ24((H2⋄H2⋄H2))}h2.LetA and B be(5.2)A=8γ23{(H2⋄H1⋄H2)+(H2⋄H2⋄H1)+(H1⋄H2⋄H2)},B=8γ23((H2⋄H1⋄H1)+(H1⋄H1⋄H2)+(H1⋄H2⋄H1))+48γ24((H2⋄H2⋄H2)).The problem of the blind nonlinear identification will be expressed as(5.3)Ah1+Bh2=vec⁡(CU).This system can be solved using several methods. We propose to resolve it using the Iterative Alternating Least Square algorithm (IALS).


## 6. Cost Functions and Iterative Alternating Least Square Algorithm

To apply the IALS algorithm, we suppose alternatively thatAh1 or Bh2 is a constant vector. Then, we get two cost functions to be minimized. Assuming that the vector Bh2 is constant, the first cost function will be expressed by(6.1)E1(h1)=‖(vec⁡(CU)-Bh2)-Ah1‖2.For the second cost function, we assume thatAh1 is constant; thus(6.2)E2(h2)=‖(vec⁡(CU)-Ah1)-Bh2‖2.The application of the least mean squares algorithm to these two functions leads to the following solutions:(6.3)h1=A#(vec⁡(CU)-Bh2),h2=B#(vec⁡(CU)-Ah1),
where the subscript # denotes the Moore-Penrose pseudoinverse of the corresponding matrix.Finally, the different steps of the proposed IALS algorithm are summarized in Algorithm1.Algorithm 1:Different steps of the new blind identification algorithm-based cumulant tensor analysis.
Initializeh1 and h2 as random variables (estimates ĥ1(0);ĥ2(0)).Forn≥1,(i) build Hankel matrices using (4.5), for p=1;2,Ĥ1(n)=ℋ(ĥ1(n)),Ĥ2(n)=ℋ(ĥ2(n)),(ii) compute matrices estimateÂ and B̂ asÂ(n)=8γ23{(Ĥ2⋄Ĥ1⋄Ĥ2)+(Ĥ2⋄Ĥ2⋄Ĥ1)+(Ĥ1⋄Ĥ2⋄Ĥ2)}(n),B̂(n)={8γ23((Ĥ2⋄Ĥ1⋄Ĥ1)+(Ĥ1⋄Ĥ1⋄Ĥ2)+(Ĥ1⋄Ĥ2⋄Ĥ1))+48γ24((Ĥ2⋄Ĥ2⋄Ĥ2))}(n),(iii) minimize cost functions (6.1) and (6.2) so thatĥ1(n+1)=Â(n)#(vec(CU)-B̂(n)ĥ2(n)),ĥ2(n+1)=B̂(n)#(vec(CU)-Â(n)ĥ1(n)),(iv) reiterate until parametric error convergence∥vec((ĥ1(n+1)ĥ2(n+1))-(ĥ1(n)ĥ2(n)))∥∥vec(ĥ1(n+1)ĥ2(n+1))∥≤ε.The notationx̂ stands for the estimates of the parameter x.


## 7. Convergence Analysis

Equation (6.3) shows that the ALS algorithm converges to optimal solutions if and only if the Moore-Penroze pseudoinverse matrices A# and B# exist, which implies that matrices A and B must be full rank [14, 16]. To do this, we start by affirming that, due to the Hankel structure and the assumption that hi(M)≠0 (3.1), each of the matrices H1 and H2 is full rank. Then(7.1)rank⁡(H1)=rank⁡(H2)=M+1.Let us now find out the rank of matricesHi⋄Hj⋄Hk;i,j,k∈{1,2} obtained from Khatri-Rao product (5.1). We will make use of the following definition and property defining the k-rank of a matrix and the rank of a Khatri-Rao product of two matrices [17].Definition 7.1.
The rank of a matrixA∈ℂE×F (denoted by kA) is equal to k if and only if every k columns of A are linearly independent. Note that kA≤min(E,F), for all A.This means that the rank of the matrixA is the largest integer k for which every set containing k columns of A is independent.Property 3.
Consider the Khatri-Rao productA⋄B, where A is E×F and B is G×F. If neither A nor B contains a zero column (and hence kA≥1,kB≥1), then kA⋄B≥min(kA+kB-1,F).It follows that(7.2)kHi⋄Hj≥min⁡⁡(2M+1,M+1),∀i,j∈{1,2},
which is equivalent to(7.3)kHi⋄Hj≥M+1.Due to the definition of the Khatri-Rao product and the structure of the Hankel matricesHi;i∈{1,2}, we conclude that
(7.4)kHi⋄Hj=M+1=rank⁡(Hi⋄Hj).Consequently,(7.5)k(Hi⋄Hj)⋄Hk=M+1=rank⁡(Hi⋄Hj⋄Hk),∀i,j,k∈{1,2},
which means that each matrix Hi⋄Hj⋄Hk is full rank whatever the values taken by i and j in the set {1,2}.Let us now find out the rank of matricesA and B. For this purpose, we will study the structure of the matrix Hi⋄Hj⋄Hk. Recall that Hi is a (2M+1,M+1)matrix. Let Θi=[00⋯0]T be the zero column vector of dimension (2M+1); i=1,…,M.Then, for the matrixHi⋄Hj⋄Hk, we will have the following form:(7.6)Hi⋄Hj⋄Hk=(Θ1Θ1Θ1XMΘ2⋮XM-1Θ1⋮ΘM-1⋯Θ2Θ2ΘMX1⋮⋮X0ΘMΘM-1ΘM),whereXi stands for the column vector of dimension (2M+1)2(M+1) which is constituted by products of the kernels model arising from computation of the Khatri-Rao matrix product. We have seen that Hi⋄Hj⋄Hk is full rank. The sum of different matrices Hi⋄Hj⋄Hk has the same form of Hi⋄Hj⋄Hk whatever the system order and the values taken by i,j and k. Consequently, matrices A and B are full rank and then their pseudoinverse exist. We conclude that the IALS converges to an optimal solution in least mean squares sense.


## 8. Simulation Results

In this section, simulation results will be given to illustrate the performance of the proposed algorithm. Two identification Volterra-Hammerstein systems are considered:(8.1)System1:{z(n)=u(n)-0.25u(n-1)+0.9u(n-2)+u2(n)+0.5u2(n-1)-0.35u2(n-2),y(n)=z(n)+e(n),System2:{z(n)=u(n)-0.4u(n-1)+0.5u(n-2)+0.95u(n-3)+u2(n)+0.2u2(n-1)-0.7u2(n-2)+0.6u2(n-3),y(n)=z(n)+e(n).The input sequenceu(n) is assumed to be stationary, zero mean, white Gaussian noise with variance γ2=1. The noise signal e(n) is also assumed to be white Gaussian sequence and independent of the input. The parameter estimation was performed for two different signal-to-noise ratio (SNR) levels: 20 dB and 3 dB.The SNR is computed with the following expression:(8.2)SNR=E(z2(n))E(e2(n)).Fourth-order cumulants were estimated from different lengths of output sequences (N=4096 and N=16384) assuming perfect knowledge of the system model. To reduce the realization dependency, parameters were averaged over 500 Monte-Carlo runs. For each simulation, we give the curves representing the variation of the estimates along with the Monte-Carlo runs, and we resume exclusive results in different tables.System 1
Figures3 and 4 show the estimates of the different kernels of the proposed model, with the IALS algorithm, for N=4096 and for different SNR levels (3 dB and 20 dB).
The mean and the standard deviation of the estimated kernels against the true ones are shown in Table1.
Likewise, Figures5 and 6 show the estimates of the different kernels of System 1 for N=16384 and for SNR levels equal to 3 dB and 20 dB, while, in Table 2, the mean and the standard deviation of the estimated kernels against the true ones are shown.
From these results, we observe that the proposed IALS algorithm performs well generating estimates for a large variation of the SNR (from 20 dB to 3 dB). We also note that the standard deviation is relatively large and decreases with the number of the system observations.Table 1
True and estimated values of the kernels of System 1 forN=4096 (500 Monte-Carlo runs).
True parameters3 dB20 dBMeanSt. Dev.MeanSt. Dev.−0.25−0.24010.3242−0.24560.17780.90.99850.37730.88370.21430.50.51280.09450.48150.0589−0.35−0.34110.1046−0.35210.0655Table 2
True and estimated values of the kernels of System 2 forN=16384 (500 Monte-Carlo runs).
True parameters3 dB20 dBMeanSt. Dev.MeanSt. Dev.−0.25−0.24640.2520−0.25890.15760.90.92380.29470.88840.19760.50.49000.06790.50280.0482−0.35−0.36650.07530.34780.0508Estimates of the parameters of System 1 with the IALS algorithm forN=4096 and SNR=3dB.
(a)(b)(c)(d)Estimates of the parameters of System 1 with the IALS algorithm forN=4096 and SNR=20dB.
(a)(b)(c)(d)Estimates of the parameters of System 1 with the IALS algorithm forN=16384 and SNR=3dB.
(a)(b)(c)(d)Estimates of the parameters of System 1 with the IALS algorithm forN=16384 and SNR=20dB.
(a)(b)(c)(d)System 2
Figures7 and 8 show the estimates of the different kernels of the second proposed model for N=4096 and for different SNR (20 dB and 3 dB).
The mean and the standard deviation of the estimated kernels against the true ones are shown in Table3.
Figures9 and 10 show the estimates of the different kernels of System 2 for N=16384 and for different SNR (20 dB and 3 dB), while, in Table 4, the mean and the standard deviation of the estimated kernels against the true ones are shown. The mean and the standard deviation of the estimated kernels against the true ones, for the second system, are shown in Table 4.
From these results, we note also that the proposed algorithm provides good estimates for the proposed system. The number of observationsN affects the range of variation of the standard deviation values. Indeed, for important values of N, this range becomes so small. The method provides good estimates even for low levels of SNR. Furthermore, we note that the larger the Monte-Carlo runs number, the smaller the standard deviations are.Table 3
True and estimated values of the kernels of System 2 forN=4096 (500 Monte-Carlo runs).
True parameters3 dB20 dBMeanSt. Dev.MeanSt. Dev.−0.4−0.45120.6550−0.42930.15940.50.50740.76320.51730.19680.951.15560.53290.93900.15450.20.23580.25250.21570.0958−0.7−0.67940.2855−0.68140.12600.60.58720.22680.59620.0986Table 4
True and estimated values of the kernels of System 2 forN=16384(500 monte-carlo runs).
True parameters3 dB20 dBMeanSt. Dev.MeanSt. Dev.−0.4−0.43780.3824−0.41210.12580.50.51160.53460.51220.14320.950.94500.46140.95070.10730.20.18470.11200.21430.0750−0.7−0.68770.1796−0.69410.07720.60.60140.13230.59830.0604Estimates of the parameters of System 2 with the IALS algorithm forN=4096 and SNR=3dB.
(a)(b)(c)(d)(e)(f)Estimates of the parameters of System 2 with the IALS algorithm forN=4096 and SNR=20dB.
(a)(b)(c)(d)(e)(f)Estimates of the parameters of System 2 with the IALS algorithm forN=16384 and SNR=3dB.
(a)(b)(c)(d)(e)(f)Estimates of the parameters of System 2 with the IALS algorithm forN=16384 and SNR=20dB.
(a)(b)(c)(d)(e)(f)


### 8.1. Comparison with Existing Methods

The performance of the previous algorithm was compared with two works: the algorithm proposed in [9] (will be noted as BIL to blind identification with linearization) and the Lagrange Programming Neural Network (LPNN) proposed in [13].(i)
In [14], the problem of blind identification was converted into a linear multivariable form using Kronecker product of the output cumulants. This can be described by the following equations:
(8.3)C¯yk(q,τ)=b(τ)(Γ¯kw)p×pbT(q),
where Cyk(τ1,τ2,…,τk-1) denotes the output cumulants sequence of order k, Γkw is the intensity (zero lag cumulant) of order k of the vector W which is formed in its turn by the different powers of input, and b is the kernel vector.For τ=0, this becomes
(8.4)C¯yk(q,0)=b(0)(Γ¯kw)p×pbT(q)=(Γ¯kw)p×pbT(q).
Different important scenarios were discussed and successfully resolved. Here, we are interested in the case of Gaussian input when the input statistics are known. Despite the efficiency of the proposed method, the resulting algorithms are in general cumbersome especially for the high series order (As confirmed by authors). For more details, see [14].(ii)
In their work [13], authors tried to determine the different Volterra kernels and the variance of the input from the autocorrelation estimates ρ[k] and the third-order moments estimates μ[k,l] of the system output, using the Lagrange Programming Neural Network (LPNN). As the LPNN is essentially designed for general nonlinear programming, they expressed the identification problem as follows:
(8.5)Minimize:L(f)=∑i∑j(μ[i,j]-M[i,j,f])2,Subjectto:ρ[i]=R[i,j],
where R[i,j] is the autocorrelation function of the real process y[n] and M[i,j] is its third order moment sequence. f is the vector formed by the unknown parameters of the Volterra model and the unknown variance of the driving noise.So the Lagrangian function will be written as
(8.6)L(f,λ)=L(f)+∑iλi(ρ[i]-R[i,f]).To improve the convergence and the precision of the algorithm, authors extended the preceding function by defining the Augmented Lagrangian Function such as(8.7)L(f,λ)=L(f)+∑iλi(ρ[i]-R[i,f])+β∑i(ρ[i]-R[i,f])2,
where {βk} is a penalty parameter sequence satisfying 0<βk<βk+1forallk,βk→∞. So the back-propagation algorithm can be established using the Lagrange multiplier.The performance of the new proposed algorithm was compared with these two algorithms. Each of these algorithms was used to identify the two models presented above (8.1), for the case of Gaussian excitation, N=16384 samples, and for the tow proposed SNR levels 3 dB and 20 dB.Figures11 and 12 show a comparison of the standard deviations given by each algorithm. We note that these results may vary considerably depending on the number of the output observations. These results show that the new proposed algorithm performs well. For a small number of unknown parameters, we note that all algorithms give in general the same STD values and these values decrease by increasing the SNR values. We note furthermore that BIL algorithm is so complex for programming in comparison with the LPNN and the new one. For big number of unknown parameters, the BIL algorithm becomes very computational complex and even the LPNN, while the new algorithm keeps its simplicity and provides good parameters with very competitive STD values.Comparison of the standard deviations (STDs) of the new algorithm against those of the LPNN and the BIL algorithms: System 1.
(a)(b)Comparison of the standard deviations (STDs) of the new algorithm against those of the LPNN and the BIL algorithms: System 2.
(a)(b)


## 8.1. Comparison with Existing Methods

The performance of the previous algorithm was compared with two works: the algorithm proposed in [9] (will be noted as BIL to blind identification with linearization) and the Lagrange Programming Neural Network (LPNN) proposed in [13].(i)
In [14], the problem of blind identification was converted into a linear multivariable form using Kronecker product of the output cumulants. This can be described by the following equations:
(8.3)C¯yk(q,τ)=b(τ)(Γ¯kw)p×pbT(q),
where Cyk(τ1,τ2,…,τk-1) denotes the output cumulants sequence of order k, Γkw is the intensity (zero lag cumulant) of order k of the vector W which is formed in its turn by the different powers of input, and b is the kernel vector.For τ=0, this becomes
(8.4)C¯yk(q,0)=b(0)(Γ¯kw)p×pbT(q)=(Γ¯kw)p×pbT(q).
Different important scenarios were discussed and successfully resolved. Here, we are interested in the case of Gaussian input when the input statistics are known. Despite the efficiency of the proposed method, the resulting algorithms are in general cumbersome especially for the high series order (As confirmed by authors). For more details, see [14].(ii)
In their work [13], authors tried to determine the different Volterra kernels and the variance of the input from the autocorrelation estimates ρ[k] and the third-order moments estimates μ[k,l] of the system output, using the Lagrange Programming Neural Network (LPNN). As the LPNN is essentially designed for general nonlinear programming, they expressed the identification problem as follows:
(8.5)Minimize:L(f)=∑i∑j(μ[i,j]-M[i,j,f])2,Subjectto:ρ[i]=R[i,j],
where R[i,j] is the autocorrelation function of the real process y[n] and M[i,j] is its third order moment sequence. f is the vector formed by the unknown parameters of the Volterra model and the unknown variance of the driving noise.So the Lagrangian function will be written as
(8.6)L(f,λ)=L(f)+∑iλi(ρ[i]-R[i,f]).To improve the convergence and the precision of the algorithm, authors extended the preceding function by defining the Augmented Lagrangian Function such as(8.7)L(f,λ)=L(f)+∑iλi(ρ[i]-R[i,f])+β∑i(ρ[i]-R[i,f])2,
where {βk} is a penalty parameter sequence satisfying 0<βk<βk+1forallk,βk→∞. So the back-propagation algorithm can be established using the Lagrange multiplier.The performance of the new proposed algorithm was compared with these two algorithms. Each of these algorithms was used to identify the two models presented above (8.1), for the case of Gaussian excitation, N=16384 samples, and for the tow proposed SNR levels 3 dB and 20 dB.Figures11 and 12 show a comparison of the standard deviations given by each algorithm. We note that these results may vary considerably depending on the number of the output observations. These results show that the new proposed algorithm performs well. For a small number of unknown parameters, we note that all algorithms give in general the same STD values and these values decrease by increasing the SNR values. We note furthermore that BIL algorithm is so complex for programming in comparison with the LPNN and the new one. For big number of unknown parameters, the BIL algorithm becomes very computational complex and even the LPNN, while the new algorithm keeps its simplicity and provides good parameters with very competitive STD values.Comparison of the standard deviations (STDs) of the new algorithm against those of the LPNN and the BIL algorithms: System 1.
(a)(b)Comparison of the standard deviations (STDs) of the new algorithm against those of the LPNN and the BIL algorithms: System 2.
(a)(b)


## 9. Conclusion

In this paper, a new approach for blind nonlinear identification problem of a second-order Hammerstein-Volterra system is developed. Thanks to a matrix analysis of a cubic tensor composed of the fourth-order output cumulants, the nonlinear identification problem is reduced to a system having the following general form:Ax+By=c. This system is solved using the Iterative Alternating Least Square. A convergence analysis shows that matrices A and B are full rank which means that the IALS algorithm converges to optimal solutions in the least mean squares sense. Simulation results on two different systems show good performance of the proposed algorithm. It is noted also that the different values of the estimates improve with the number of the system observation even for small values of SNR. Comparison results with two algorithms show that the new proposed algorithm performs well and especially in the case of great number of unknown parameters. Extending the proposed algorithm for more input classes and for more general Volterra-Hammerstein systems remains an open problem, and it is now the subject matter of current works.



---

*Source: 284815-2012-06-13.xml*