Identification of Multilinear Systems: A Brief Overview

Abstract

Nonlinear systems have been studied for a long time and have applications in numerous research fields. However, there is currently no global solution for nonlinear system identification, and different used approaches depend on the type of nonlinearity. An interesting class of nonlinear systems, with a wide range of popular applications, is represented by multilinear (or multidimensional) systems. These systems exhibit a particular property that may be exploited, namely that they can be regarded as linearly separable systems and can be modeled accordingly, using tensors. Examples of well-known applications of multilinear forms are multiple-input/single-output (MISO) systems and acoustic echo cancellers, used in multi-party voice communications, such as videoconferencing. Many important fields (e.g., big data, machine learning, and source separation) can benefit from the methods employed in multidimensional system identification. In this context, this chapter aims to briefly present the recent approaches in the identification of multilinear systems. Methods relying on tensor decomposition and modeling are used to address the large parameter space of such systems.

1. Introduction

System identification is an important topic nowadays since it can be used in solving numerous problems [1]. The aim of system identification is to estimate an unknown model using the available and observed data, namely the input and output of the system. In this context, the well-known Wiener filter is a popular solution, along with the adaptive filters which can be derived starting from this approach.

In multilinear system identification, dealing with a large parameter space represents an important challenge [2, 3]. The huge length of the filter (hundreds or thousands of coefficients) is also a serious problem [4, 5]. The methods used for addressing these issues usually rely on tensor decomposition and modeling [2, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17], meaning that a high-dimension problem is rewritten as a combination of lower-dimension structures, using the Kronecker product decomposition [18].

In the context of multilinear forms identification, a few approaches were proposed recently, addressing the cases when the large system is decomposed into two or three smaller components (i.e., bilinear and trilinear forms, respectively) [10, 17, 18, 19, 20, 21, 22, 23]. The aforementioned solutions outperform their conventional counterparts, offering at the same time a lower computational complexity.

Motivated by the appealing performance of these previous developments, we extended the tensor decomposition technique to higher-order systems, and in this framework, this chapter presents a part of the work and results obtained recently by the authors in the context of multilinear system identification. An iterative Wiener filter and a family of LMS-based algorithms tailored for multilinear forms are presented. For more details on the results summarized here, the works [24, 25, 26] can be consulted.

Related to the work presented here, several other tensor-based solutions relying on the recursive least-squares (RLS) algorithm were also developed recently [27, 28]. Possible applications of such system identification frameworks can be encountered in topics such as big data [29], machine learning [14], but they may be also useful in nonlinear acoustic echo cancelation [30, 31], source separation [13, 32, 33], channel equalization [12, 34], array beamforming [16, 35], blind identification [36], object recognition [37, 38], and cardiac applications [39].

The rest of this chapter is organized in the following way. In Section 2, we introduce the system model for the multiple-input/single-output (MISO) system identification problem. In this context, Section 3 presents an iterative Wiener filter tailored for the identification of multilinear systems. Next, in Section 4, an LMS-based algorithm is presented, together with its normalized version, and then, in Section 5, the performance of these algorithms is proved through simulations. Finally, conclusions are drawn in Section 6.

2. System model in the multilinear framework

Let us consider a MISO system, whose output signal at the time index $t$ can be written as

$$y\left(t\right)=\sum _{l_1=1}^{L_1}\sum _{l_2=1}^{L_2}\cdots \sum _{l_N=1}^{L_N}{x}_{l_1{l}_2{l}_3\dots l_N}\left(t\right)h_{1,l_1}{h}_{2,l_2}\cdots h_{N,l_N},$$E1

where the individual channels are modeled by the vectors:

$${\mathbf{h}}_i={\left[h_{i,1}{h}_{i,2}\cdots h_{i,L_i}\right]}^T,i=1,2,\dots ,N,$$E2

the superscript ^T denotes the transpose operator, and the input signals may be expressed in a tensorial form as $\mathcal{X}\left(t\right)\in {\mathbb{R}}^{L_1\times L_2\times L_3\times \cdots \times L_N}$, with the elements ${\left(\mathcal{X}\right)}_{l_1{l}_2\dots l_N}\left(t\right)=x_{l_1{l}_2\dots l_N}\left(t\right)$. Consequently, the output signal becomes

$$y\left(t\right)=\mathcal{X}\left(t\right){\times }_1{\mathbf{h}}_1^T{\times }_2{\mathbf{h}}_2^T{\times }_3\cdots {\times }_N{\mathbf{h}}_N^T,$$E3

where ${\times }_i$ (for $i=1,2,\dots ,N$) denotes the mode-$i$ product [7]. It can be said that $y\left(t\right)$ is a multilinear form because it is a linear function of each of the vectors ${\mathbf{h}}_i,i=1,2,\dots ,N$, when the other $N-1$ vectors are fixed. In this context, $y\left(t\right)$ may be regarded as an extension of the bilinear form [19]. Next, let us define

$$\mathcal{H}={\mathbf{h}}_1\circ {\mathbf{h}}_2\circ \cdots \circ {\mathbf{h}}_N,$$E4

where $\circ$ is the vector outer product, i.e., ${\mathbf{h}}_1\circ {\mathbf{h}}_2={\mathbf{h}}_1{\mathbf{h}}_2^T$, ${\left({\mathbf{h}}_1\circ {\mathbf{h}}_2\right)}_{i,j}=h_{1,i}{h}_{2,j}$, $\mathrm{vec}\left({\mathbf{h}}_1\circ {\mathbf{h}}_2\right)={\mathbf{h}}_2\otimes {\mathbf{h}}_1$, and

$${\left(\mathcal{H}\right)}_{l_1,l_2,\dots ,l_N}=h_{1,l_1}{h}_{2,l_2}\cdots h_{N,l_N},$$E5

$$\mathrm{vec}\left(\mathcal{H}\right)={\mathbf{h}}_N\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{h}}_1,$$E6

where $\otimes$ denotes the Kronecker product and $\mathrm{vec}\left(\cdot \right)$ is the vectorization operation:

$$\mathrm{vec}\left(\mathcal{H}\right)=\left[\begin{array}{c}\mathrm{vec}\left({\mathbf{H}}_{::\dots :1}\right)\\ ⋮\\ \mathrm{vec}\left({\mathbf{H}}_{::\dots :L_N}\right)\end{array}\right],$$E7

$$\mathrm{vec}\left({\mathbf{H}}_{::\dots :l_i}\right)=\left[\begin{array}{c}\mathrm{vec}\left({\mathbf{H}}_{::\dots :1,i}\right)\\ ⋮\\ \mathrm{vec}\left({\mathbf{H}}_{::\dots :L_{N-1,i}}\right)\end{array}\right],$$E8

and so on, where ${\mathbf{H}}_{::\dots :l_i}\in {\mathbb{R}}^{L_1\times L_2\times L_3\times \cdots \times L_{N-1}}$ represent the frontal slices of the tensor $\mathcal{H}$. Therefore, the output signal can be expressed as

$$y\left(t\right)={\mathrm{vec}}^T\left(\mathcal{H}\right)\mathrm{vec}\left[\mathcal{X}\left(t\right)\right],$$E9

where

$$\mathrm{vec}\left[\mathcal{X}\left(t\right)\right]=\left[\begin{array}{c}\mathrm{vec}\left[{\mathbf{X}}_{::\dots :1}\left(t\right)\right]\\ ⋮\\ \mathrm{vec}\left[{\mathbf{X}}_{::\dots :L_N}\left(t\right)\right]\end{array}\right]=\mathbf{x}\left(t\right),$$E10

with ${\mathbf{X}}_{::\dots :l_i}\left(t\right)\in {\mathbb{R}}^{L_1\times L_2\times L_3\times \cdots \times L_{N-1}}$ being the frontal slices of the tensor $\mathcal{X}\left(t\right)$. Let us denote the global impulse response of length $L_1{L}_2\cdots L_N$ as

$$\mathbf{g}=\mathrm{vec}\left(\mathcal{H}\right)={\mathbf{h}}_N\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{h}}_1.$$E11

Here, an observation can be made: the solution of the decomposition in Eq. (11) is not unique [17, 24]. Despite this, no scaling ambiguity occurs in the identification of the global impulse response, $\mathbf{g}$.

Using Eqs. (9)–(11), we may rewrite $y\left(t\right)$ as

$$y\left(t\right)={\mathbf{g}}^T\mathbf{x}\left(t\right).$$E12

We aim to identify the global impulse response, $\mathbf{g}$. We can define the reference (or desired) signal as

$$d\left(t\right)={\mathbf{g}}^T\mathbf{x}\left(t\right)+w\left(t\right),$$E13

where $w\left(t\right)$ denotes the additive noise, which is uncorrelated with the input signals, and whose variance is

$${\sigma }_d^2={\mathbf{g}}^{T}E\left[\mathbf{x}\left(t\right){\mathbf{x}}^T\left(t\right)\right]\mathbf{g}+{\sigma }_w^2={\mathbf{g}}^T\mathbf{Rg}+{\sigma }_w^2,$$E14

with $E\left[\cdot \right]$ denoting mathematical expectation, $\mathbf{R}=E\left[\mathbf{x}\left(t\right){\mathbf{x}}^T\left(t\right)\right]$, and ${\sigma }_w^2=E\left[w^2\left(t\right)\right]$. Next, the error signal can be defined as

$$e\left(t\right)=d\left(t\right)-{\widehat{\mathbf{g}}}^T\mathbf{x}\left(t\right),$$E15

where $\widehat{\mathbf{g}}$ denotes an estimate of the global impulse response.

The optimization criterion is the minimization of the mean-squared error (MSE), which can be defined using Eq. (15):

$$J\left(\widehat{\mathbf{g}}\right)=E\left[e^2\left(t\right)\right]={\sigma }_d^2-2{\widehat{\mathbf{g}}}^T\mathbf{p}+{\widehat{\mathbf{g}}}^T\mathbf{R}\widehat{\mathbf{g}},$$E16

where $\mathbf{p}=E\left[d\left(t\right)\mathbf{x}\left(t\right)\right]$ denotes the cross-correlation vector between $d\left(t\right)$ and $\mathbf{x}\left(t\right)$. The solution to this minimization problem is given by the popular Wiener filter [40]:

$${\widehat{\mathbf{g}}}_{\mathrm{W}}={\mathbf{R}}^{-1}\mathbf{p}.$$E17

Relation (17) provides the global impulse response. In order to obtain the $N$ coefficient vectors ${\mathbf{h}}_i,i=1,2,\dots ,N$, the nonlinear equation set containing $L_1{L}_2\cdots L_N$ equations with $L_1+L_2+\cdots +L_N$ scalar variables needs to be solved:

$${\widehat{\mathbf{g}}}_{\mathrm{W}}={\widehat{\mathbf{h}}}_{\mathrm{W},N}\otimes {\widehat{\mathbf{h}}}_{\mathrm{W},N-1}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_{\mathrm{W},1}.$$E18

3. Multilinear iterative Wiener filter

It can be easily checked that

$$\begin{array}{c}\mathbf{g}={\mathbf{h}}_N\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{h}}_1\\ =\left({\mathbf{h}}_N\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{I}}_{L_1}\right){\mathbf{h}}_1\\ =\left({\mathbf{h}}_N\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{h}}_3\otimes {\mathbf{I}}_{L_2}\otimes {\mathbf{h}}_1\right){\mathbf{h}}_2\\ ⋮\\ =\left({\mathbf{h}}_N\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{I}}_{L_i}\otimes {\mathbf{h}}_{L_i-1}\otimes \cdots \otimes {\mathbf{h}}_1\right){\mathbf{h}}_i\\ ⋮\\ =\left({\mathbf{I}}_{L_N}\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{h}}_1\right){\mathbf{h}}_N,\end{array}$$E19

where ${\mathbf{I}}_{L_i}$ denotes the identity matrix of size $L_i\times L_i$.

Hence, the cost function given by Eq. (16) may be expressed in $N$ equivalent forms:

$$\begin{array}{l}J\left({\widehat{\mathbf{h}}}_1,{\widehat{\mathbf{h}}}_2,\dots ,{\widehat{\mathbf{h}}}_N\right)={\sigma }_d^2-2{\widehat{\mathbf{g}}}^T\mathbf{p}+{\widehat{\mathbf{g}}}^T\mathbf{R}\widehat{\mathbf{g}}\\ ={\sigma }_d^2-2{\widehat{\mathbf{h}}}_i^T{\left({\widehat{\mathbf{h}}}_N\otimes {\widehat{\mathbf{h}}}_{N-1}\otimes \cdots \otimes {\mathbf{I}}_{L_i}\otimes {\widehat{\mathbf{h}}}_{L_i-1}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\right)}^T\mathbf{p}\\ +{\widehat{\mathbf{h}}}_i^T{\left({\widehat{\mathbf{h}}}_N\otimes {\widehat{\mathbf{h}}}_{N-1}\otimes \cdots \otimes {\mathbf{I}}_{L_i}\otimes {\widehat{\mathbf{h}}}_{L_i-1}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\right)}^T\\ \times \mathbf{R}\left({\widehat{\mathbf{h}}}_N\otimes {\widehat{\mathbf{h}}}_{N-1}\otimes \cdots \otimes {\mathbf{I}}_{L_i}\otimes {\widehat{\mathbf{h}}}_{L_i-1}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\right){\widehat{\mathbf{h}}}_i\\ ={\sigma }_d^2-2{\widehat{\mathbf{h}}}_i^T{\mathbf{p}}_i+{\widehat{\mathbf{h}}}_i^T{\mathbf{R}}_i{\widehat{\mathbf{h}}}_i,i=1,2,\dots ,N,\end{array}$$E20

where

$${\mathbf{p}}_i={\left({\widehat{\mathbf{h}}}_N\otimes {\widehat{\mathbf{h}}}_{N-1}\otimes \cdots \otimes {\mathbf{I}}_{L_i}\otimes {\widehat{\mathbf{h}}}_{L_i-1}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\right)}^T\mathbf{p},$$E21

$$\begin{array}{l}{\mathbf{R}}_i={\left({\widehat{\mathbf{h}}}_N\otimes {\widehat{\mathbf{h}}}_{N-1}\otimes \cdots \otimes {\mathbf{I}}_{L_i}\otimes {\widehat{\mathbf{h}}}_{L_i-1}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\right)}^T\mathbf{R}\\ \times \left({\widehat{\mathbf{h}}}_N\otimes {\widehat{\mathbf{h}}}_{N-1}\otimes \cdots \otimes {\mathbf{I}}_{L_i}\otimes {\widehat{\mathbf{h}}}_{L_i-1}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\right).\end{array}$$E22

If all coefficients except ${\widehat{\mathbf{h}}}_i$ are kept fixed, we may define

$$\begin{array}{c}{J}_{{\widehat{\mathbf{h}}}_1,{\widehat{\mathbf{h}}}_2,\dots ,{\widehat{\mathbf{h}}}_{i-1},{\widehat{\mathbf{h}}}_{i+1},\dots ,{\widehat{\mathbf{h}}}_N}\left({\widehat{\mathbf{h}}}_i\right)={\sigma }_d^2-2{\widehat{\mathbf{h}}}_i^T{\mathbf{p}}_i+{\widehat{\mathbf{h}}}_i^T{\mathbf{R}}_i{\widehat{\mathbf{h}}}_i,\\ i=1,2,\dots ,N.\end{array}$$E23

The minimization of this convex cost function with respect to ${\widehat{\mathbf{h}}}_i$ yields

$${\widehat{\mathbf{h}}}_i={\mathbf{R}}_i^{-1}{\mathbf{p}}_i,i=1,2,\dots ,N.$$E24

Using this result, an iterative approach can be derived. A set of initial values ${\widehat{\mathbf{h}}}_i^{\left(0\right)},i=1,2,\dots ,N$ are chosen for starting the algorithm, and then we can compute

$${\mathbf{p}}_1^{\left(0\right)}={\left({\widehat{\mathbf{h}}}_N^{\left(0\right)}\otimes {\widehat{\mathbf{h}}}_{N-1}^{\left(0\right)}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_2^{\left(0\right)}\otimes {\widehat{\mathbf{I}}}_{L_1}\right)}^T\mathbf{p},$$E25

$${\mathbf{R}}_1^{\left(0\right)}={\left({\widehat{\mathbf{h}}}_N^{\left(0\right)}\otimes {\widehat{\mathbf{h}}}_{N-1}^{\left(0\right)}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_2^{\left(0\right)}\otimes {\widehat{\mathbf{I}}}_{L_1}\right)}^T\mathbf{R}\left({\widehat{\mathbf{h}}}_N^{\left(0\right)}\otimes {\widehat{\mathbf{h}}}_{N-1}^{\left(0\right)}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_2^{\left(0\right)}\otimes {\widehat{\mathbf{I}}}_{L_1}\right),$$E26

$$J_{{\widehat{\mathbf{h}}}_2,{\widehat{\mathbf{h}}}_3,\dots ,{\widehat{\mathbf{h}}}_N}\left({\widehat{\mathbf{h}}}_1^{\left(1\right)}\right)={\sigma }_d^2-2{\left({\widehat{\mathbf{h}}}_1^{\left(1\right)}\right)}^T{\mathbf{p}}_1^{\left(0\right)}+{\left({\widehat{\mathbf{h}}}_1^{\left(1\right)}\right)}^T{\mathbf{R}}_1^{\left(0\right)}\left({\widehat{\mathbf{h}}}_1^{\left(1\right)}\right).$$E27

The minimization of the cost function yields

$${\widehat{\mathbf{h}}}_1^{\left(1\right)}={\left({\mathbf{R}}_1^{\left(0\right)}\right)}^{-1}{\mathbf{p}}_1^{\left(0\right)}.$$E28

Using ${\widehat{\mathbf{h}}}_1^{\left(1\right)}$ and ${\widehat{\mathbf{h}}}_i^{\left(0\right)}$, $i=3,\dots ,N$, we can now compute ${\widehat{\mathbf{h}}}_2^{\left(1\right)}$. Then, the cost function becomes

$$J_{{\widehat{\mathbf{h}}}_1,{\widehat{\mathbf{h}}}_3,\dots ,{\widehat{\mathbf{h}}}_N}\left({\widehat{\mathbf{h}}}_2^{\left(1\right)}\right)={\sigma }_d^2-2{\left({\widehat{\mathbf{h}}}_2^{\left(1\right)}\right)}^T{\mathbf{p}}_2^{\left(1\right)}+{\left({\widehat{\mathbf{h}}}_2^{\left(1\right)}\right)}^T{\mathbf{R}}_2^{\left(1\right)}\left({\widehat{\mathbf{h}}}_2^{\left(1\right)}\right),$$E29

where

$${\mathbf{p}}_2^{\left(1\right)}={\left({\widehat{\mathbf{h}}}_N^{\left(0\right)}\otimes {\widehat{\mathbf{h}}}_{N-1}^{\left(0\right)}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_3^{\left(0\right)}\otimes {\widehat{\mathbf{I}}}_{L_2}\otimes {\widehat{\mathbf{h}}}_1^{\left(1\right)}\right)}^T\mathbf{p},$$E30

$${\mathbf{R}}_2^{\left(1\right)}={\left({\widehat{\mathbf{h}}}_N^{\left(0\right)}\otimes {\widehat{\mathbf{h}}}_{N-1}^{\left(0\right)}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_3^{\left(0\right)}\otimes {\widehat{\mathbf{I}}}_{L_2}{\widehat{\mathbf{h}}}_1^{\left(1\right)}\right)}^T\mathbf{R}\left({\widehat{\mathbf{h}}}_N^{\left(0\right)}\otimes {\widehat{\mathbf{h}}}_{N-1}^{\left(0\right)}\otimes \cdots \otimes {\widehat{\mathbf{h}}}_3^{\left(0\right)}\otimes {\widehat{\mathbf{I}}}_{L_2}{\widehat{\mathbf{h}}}_1^{\left(1\right)}\right).$$E31

The minimization of the cost function yields

$${\widehat{\mathbf{h}}}_2^{\left(1\right)}={\left({\mathbf{R}}_2^{\left(1\right)}\right)}^{-1}{\mathbf{p}}_2^{\left(1\right)}.$$E32

All the other estimates ${\widehat{\mathbf{h}}}_i^{\left(1\right)},i=3,4,\dots ,N$ can be computed in a similar manner. By further iterating up to iteration $n$, the estimates of the $N$ vectors are obtained. This minimization technique is called “block coordinate descent” [41].

4. LMS and NLMS algorithms for multilinear forms

The limitations of the Wiener filter (e.g., matrix inversion, statistics estimation) can restrict the applicability of the previously presented approach in real-world situations (for example, in nonstationary conditions, or when real-time processing is needed). Therefore, a better approach may be represented by adaptive filters. In this context, the well-known least-mean-square (LMS) algorithm is among the most popular solutions, due to its simplicity. In the following, a family of LMS-based algorithms for multilinear forms identification is presented.

By using the estimated impulse responses ${\widehat{\mathbf{h}}}_i\left(t\right),i=1,2,\dots ,N$, the corresponding a priori error signals can be defined

$$e_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)=d\left(t\right)-{\widehat{\mathbf{h}}}_1^T\left(t-1\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right),$$E33

$$e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)=d\left(t\right)-{\widehat{\mathbf{h}}}_2^T\left(t-1\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right),$$E34

$$⋮$$

$$e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)=d\left(t\right)-{\widehat{\mathbf{h}}}_N^T\left(t-1\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right),$$E35

where

$${\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)=\left[\right.{\widehat{\mathbf{h}}}_N\left(t-1\right)\otimes {\widehat{\mathbf{h}}}_{N-1}\left(t-1\right)\otimes \cdots \otimes {\widehat{\mathbf{h}}}_2\left(t-1\right)\otimes {\mathbf{I}}_{L_1}\left.\right]\mathbf{x}\left(t\right),$$E36

$${\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)=\left[\right.{\widehat{\mathbf{h}}}_N\left(t-1\right)\otimes {\widehat{\mathbf{h}}}_{N-1}\left(t-1\right)\otimes \cdots \otimes {\widehat{\mathbf{h}}}_3\left(t-1\right)\otimes {\mathbf{I}}_{L_2}\otimes {\widehat{\mathbf{h}}}_1\left(t-1\right)\left.\right]\mathbf{x}\left(t\right),$$E37

$$⋮$$

$${\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)=\left[\right.{\mathbf{I}}_{L_N}\otimes {\widehat{\mathbf{h}}}_{N-1}\left(t-1\right)\otimes \cdots \otimes {\widehat{\mathbf{h}}}_2\left(t-1\right){\widehat{\mathbf{h}}}_1\left(t-1\right)\left.\right]\mathbf{x}\left(t\right).$$E38

We can easily check that $e_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)=e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)=\dots =e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)$. Hence, the LMS updates of the individual filters will be

$$\begin{array}{c}{\widehat{\mathbf{h}}}_1\left(t\right)={\widehat{\mathbf{h}}}_1\left(t-1\right)-\frac{{\mu }_{{\widehat{\mathbf{h}}}_1}}{2}\times \frac{\partial e_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}^2\left(t\right)}{\partial {\widehat{\mathbf{h}}}_1\left(t-1\right)}\\ ={\widehat{\mathbf{h}}}_1\left(t-1\right)+{\mu }_{{\widehat{\mathbf{h}}}_1}{\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)e_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right),\end{array}$$E39

$$\begin{array}{c}{\widehat{\mathbf{h}}}_2\left(t\right)={\widehat{\mathbf{h}}}_2\left(t-1\right)-\frac{{\mu }_{{\widehat{\mathbf{h}}}_2}}{2}\times \frac{\partial e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}^2\left(t\right)}{\partial {\widehat{\mathbf{h}}}_2\left(t-1\right)}\\ ={\widehat{\mathbf{h}}}_2\left(t-1\right)+{\mu }_{{\widehat{\mathbf{h}}}_2}{\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right),\\ ⋮\end{array}$$E40

$$\begin{array}{c}{\widehat{\mathbf{h}}}_N\left(t\right)={\widehat{\mathbf{h}}}_N\left(t-1\right)-\frac{{\mu }_{{\widehat{\mathbf{h}}}_N}}{2}\times \frac{\partial e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}^2\left(t\right)}{\partial {\widehat{\mathbf{h}}}_N\left(t-1\right)}\\ ={\widehat{\mathbf{h}}}_N\left(t-1\right)+{\mu }_{{\widehat{\mathbf{h}}}_N}{\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right),\end{array}$$E41

where ${\mu }_{{\widehat{\mathbf{h}}}_i}>0,i=1,2,\dots ,N$ represent the step-size parameters. Relations (39–41) describe the LMS algorithm for multilinear forms (LMS-MF). The initialization of the estimated impulse responses could be

$${\widehat{\mathbf{h}}}_1\left(0\right)={\left[10\cdots 0\right]}^T,$$E42

$${\widehat{\mathbf{h}}}_k\left(0\right)=\frac{1}{L_k}{\left[11\cdots 1\right]}^T,k=2,3,\dots ,N.$$E43

The global filter estimate is obtained as

$$\widehat{\mathbf{g}}\left(t\right)={\widehat{\mathbf{h}}}_N\left(t\right)\otimes {\widehat{\mathbf{h}}}_{N-1}\left(t\right)\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\left(t\right).$$E44

We may also identify the global impulse response using the classical LMS algorithm:

$$\widehat{\mathbf{g}}\left(t\right)=\widehat{\mathbf{g}}\left(t-1\right)+{\mu }_{\widehat{\mathbf{g}}}\mathbf{x}\left(t\right)e\left(t\right),$$E45

$$e\left(t\right)=d\left(t\right)-\widehat{\mathbf{g}}\left(t-1\right)\mathbf{x}\left(t\right),$$E46

where ${\mu }_{\widehat{\mathbf{g}}}$ denotes the global step-size parameter, but this would involve a single long-length adaptive filter.

When choosing the constant values of the step-size parameters from Eqs. (39)–(41), we need to consider the compromise between convergence rate and steady-state misadjustment. In certain cases, it can be more useful to have variable step-size parameters. Hence, the update equations become

$${\widehat{\mathbf{h}}}_1\left(t\right)={\widehat{\mathbf{h}}}_1\left(t-1\right)+{\mu }_{{\widehat{\mathbf{h}}}_1}\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)e_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right),$$E47

$${\widehat{\mathbf{h}}}_2\left(t\right)={\widehat{\mathbf{h}}}_2\left(t-1\right)+{\mu }_{{\widehat{\mathbf{h}}}_2}\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right),$$E48

$$⋮$$

$${\widehat{\mathbf{h}}}_N\left(t\right)={\widehat{\mathbf{h}}}_N\left(t-1\right)+{\mu }_{{\widehat{\mathbf{h}}}_N}\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right).$$E49

Then, the a posteriori error signals can be defined as

$${\epsilon }_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)=d\left(t\right)-{\widehat{\mathbf{h}}}_1^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right),$$E50

$${\epsilon }_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)=d\left(t\right)-{\widehat{\mathbf{h}}}_2^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right),$$E51

$$⋮$$

$${\epsilon }_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)=d\left(t\right)-{\widehat{\mathbf{h}}}_N^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right).$$E52

After replacing Eq. (39) in Eq. (50), Eq. (40) in Eq. (51), and Eq. (41) in Eq. (52), and then canceling the a posteriori error signals, we get

$$e_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)\left[1-{\mu }_{{\widehat{\mathbf{h}}}_1}\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)\right]=0,$$E53

$$e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)\left[1-{\mu }_{{\widehat{\mathbf{h}}}_2}\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)\right]=0,$$E54

$$⋮$$

$$e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)\left[1-{\mu }_{{\widehat{\mathbf{h}}}_N}\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)\right]=0.$$E55

We assume that $e_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)\ne 0$, $e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)\ne 0,\dots ,$$e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)\ne 0$. Hence, the step-sizes become

$${\mu }_{{\widehat{\mathbf{h}}}_1}\left(t\right)=\frac{1}{{\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)},$$E56

$${\mu }_{{\widehat{\mathbf{h}}}_2}\left(t\right)=\frac{1}{{\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)},$$E57

$$⋮$$

$${\mu }_{{\widehat{\mathbf{h}}}_N}\left(t\right)=\frac{1}{{\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)}.$$E58

In the numerators of Eqs. (56)–(58), the normalized step-size parameters $0<{\alpha }_{{\widehat{\mathbf{h}}}_i}<1,i=1,2,\dots ,N$ can be used to achieve a good compromise between convergence rate and misadjustment. Some regularization constants denoted by ${\delta }_{{\widehat{\mathbf{h}}}_i}>0,i=1,2,\dots ,N$ are also introduced in the denominators of the step-sizes in order to ensure robust adaptation. Consequently, we obtain the update equations of the normalized LMS (NLMS) algorithm for multilinear forms (NLMS-MF):

$${\widehat{\mathbf{h}}}_1\left(t\right)={\widehat{\mathbf{h}}}_1\left(t-1\right)+\frac{{\alpha }_{{\widehat{\mathbf{h}}}_1}{\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)e_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)}{{\delta }_{{\widehat{\mathbf{h}}}_1}+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_2{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)},$$E59

$${\widehat{\mathbf{h}}}_2\left(t\right)={\widehat{\mathbf{h}}}_2\left(t-1\right)+\frac{{\alpha }_{{\widehat{\mathbf{h}}}_2}{\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)}{{\delta }_{{\widehat{\mathbf{h}}}_2}+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_3\dots {\widehat{\mathbf{h}}}_N}\left(t\right)},$$E60

$$⋮$$

$${\widehat{\mathbf{h}}}_N\left(t\right)={\widehat{\mathbf{h}}}_N\left(t-1\right)+\frac{{\alpha }_{{\widehat{\mathbf{h}}}_N}{\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)e_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}\left(t\right)}}{{\delta }_{{\widehat{\mathbf{h}}}_N}+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}^T\left(t\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_1{\widehat{\mathbf{h}}}_2\dots {\widehat{\mathbf{h}}}_{N-1}}\left(t\right)}.$$E61

The initialization of the individual impulse responses can be done using Eqs. (42, 43). We may also identify the global impulse response using the regular NLMS algorithm:

$$\widehat{\mathbf{g}}\left(t\right)=\widehat{\mathbf{g}}\left(t-1\right)+\frac{{\alpha }_{\widehat{\mathbf{g}}}\mathbf{x}\left(t\right)e\left(t\right)}{{\mathbf{x}}^T\left(t\right)\mathbf{x}\left(t\right)+{\delta }_{\widehat{\mathbf{g}}}},$$E62

where $0<{\alpha }_{\widehat{\mathbf{g}}}\le 1$ denotes the normalized step-size parameter, ${\delta }_{\widehat{\mathbf{g}}}>0$ represents the regularization constant, and $e\left(t\right)$ is defined in Eq. (46). However, this would involve a single (long length) adaptive filter, with $L=L_1{L}_2\cdots L_N$.

5. Experimental results

The purpose of this section is to illustrate through simulations the improved performance of the proposed solutions for the identification of multilinear forms. We performed experiments involving MISO system identification. As input signals, we used white Gaussian noises and AR(1) processes, obtained by filtering white Gaussian noises through a first-order system with the transfer function $1/\left(1-0.99z^{-1}\right)$. The noise is white and Gaussian, having the variance equal to ${\sigma }_w^2=0.01$. We use two different orders of the system ($N=4$ and $N=5$). The first individual impulse response, ${\mathbf{h}}_1$, of length $L_1=16$, is formed using the first 16 coefficients of the first network echo path from the ITU-T G168 Recommendation [42]. The second impulse response, ${\mathbf{h}}_2$, of length $L_2=8$, is randomly generated from a Gaussian distribution, whereas the other three impulse responses, ${\mathbf{h}}_3$, ${\mathbf{h}}_4$, and ${\mathbf{h}}_5$, of lengths $L_3=L_4=4$ and $L_5=2$, respectively, are obtained using an exponential decay based on the rule $h_{j,l_j}=a_j^{l_j-1}$, with $j=\mathrm{3,4,5}$, where $a_j$ takes the values 0.9, 0.5, and 0.1, respectively. Hence, the global impulse responses $\mathbf{g}$, computed using Eq. (11), have lengths 2048, when $N=4$, and 4096, when $N=5$. Figure 1 illustrates the individual impulse responses ${\mathbf{h}}_1$, ${\mathbf{h}}_2$, ${\mathbf{h}}_3$, and ${\mathbf{h}}_4$, as well as the resulting global impulse response $\mathbf{g}={\mathbf{h}}_4\otimes {\mathbf{h}}_3\otimes {\mathbf{h}}_2\otimes {\mathbf{h}}_1$, in the case when $N=4$.

Figure 1.

Impulse responses used in simulations of the multiple-input/single-output (MISO) system identification scenario (for $N=4$): (a) ${\mathbf{h}}_1$ (of length $L_1=16$) contains the first 16 coefficients of the first impulse response from G168 recommendation [42], (b) ${\mathbf{h}}_2$ (of length $L_2=8$) is a randomly generated impulse response, (c) ${\mathbf{h}}_3$ (of length $L_3=4$) has the coefficients computed as $h_{3,l_3}={0.9}^{l_3-1}$, with $l_3=1,2,\dots ,L_3$, (d) ${\mathbf{h}}_4$ (of length $L_4=4$) has the coefficients computed as $h_{4,l_4}={0.5}^{l_4-1}$, with $l_4=1,2,\dots ,L_4$, and (e) $\mathbf{g}$ (of length $L=L_1{L}_2{L}_3{L}_4=2048$) is the global impulse response, which results based on Eq. (11).

The measure of performance is the normalized misalignment (in dB) for the identification of the global impulse response, computed as $20{log}_{10}\left[{‖\mathbf{g}-\widehat{\mathbf{g}}\left(n\right)‖}_2/{‖\mathbf{g}‖}_2\right]$. A sudden change in the sign of the coefficients of ${\mathbf{h}}_1$ is introduced in the middle of each experiment, with the goal of observing the tracking capabilities of the proposed algorithms.

First, we aim to show comparatively the performances of the LMS-MF and LMS algorithms. When choosing the step-size parameter values, we need to take into account the theoretical upper bound, which for the conventional LMS is $2/\left(L{\sigma }_x^2\right)$, where ${\sigma }_x^2$ is the input signal’s variance [40]. In practical scenarios, this limit may not be usable, due to stability issues. The step-size parameters for the LMS-MF and LMS algorithms were chosen in our simulations in such a manner that similar misalignment values are obtained, in order to compare their convergence rate and tracking. The largest value of the step-size parameter for the conventional LMS algorithm is chosen close to its stability limit, such that the fastest convergence possible is obtained.

Figure 2 shows the case when $N=4$ and the input signals are white Gaussian noises. The resulting global filter length is $L=2048$. It can be seen that the LMS-MF achieves a higher convergence rate and faster tracking as compared to the conventional LMS algorithm. Of course, when the value of the step-size parameter decreases, the misalignment decreases, but at the same time, the convergence becomes slower.

Figure 2.

Performance of the LMS-MF and LMS algorithms (using different step-size parameters), for the identification of the global impulse response $\mathbf{g}$. The input signals are white Gaussian noises, $N=4$, and $L=2048$.

Next, in Figure 3, $N=4$ and the input signals are highly-correlated AR(1) processes. The performance gain of the LMS-MF with respect to the conventional LMS algorithm in terms of both convergence rate/tracking and steady-state misalignment is even higher in this scenario.

Figure 3.

Performance of the LMS-MF and LMS algorithms (using different step-size parameters), for the identification of the global impulse response $\mathbf{g}$. The input signals are AR(1) processes, $N=4$, and $L=2048$.

When the system order increases, the improvement in performance brought by the LMS-MF is even more apparent. This can be seen in Figure 4, where $N=5$ and the input signals are AR(1) processes. The resulting global impulse response length is $L=4096$. It is easily observed that the proposed algorithm is superior to the conventional LMS in terms of both convergence rate/tracking and misalignment.

Figure 4.

Performance of the LMS-MF and LMS algorithms (using different step-size parameters), for the identification of the global impulse response $\mathbf{g}$. The input signals are AR(1) processes, $N=5$, and $L=4096$.

In the following, we aim to illustrate the performance of the NLMS-MF and NLMS algorithms in the identification of the global system. Since the step-size parameter does no longer have a constant value, the normalized algorithms can work better in nonstationary environments. The fastest-convergence bound for the value of the normalized step-size parameter of the conventional NLMS algorithm is 1 [40].

Figure 5 illustrates the case when the inputs are white Gaussian noises and $N=4$, leading to a length of the global system of $L=2048$. Similar to the case of the LMS-MF and LMS algorithms from Figure 2, it can be concluded that the performance of the NLMS-MF algorithm is significantly better than the one of its conventional counterparts. This is even more apparent in the case of smaller normalized step-size values.

Figure 5.

Performance of the NLMS-MF and NLMS algorithms (using different normalized step-size parameters), for the identification of the global impulse response $\mathbf{g}$. The input signals are white Gaussian noises, $N=4$, and $L=2048$.

The improvement offered by the proposed approach is even more significant for correlated inputs. In Figure 6, the input signals are AR(1) processes. It is noticed that even when the NLMS-MF algorithm uses lower values for the normalized step-sizes, it can still outperform the NLMS algorithm working in the fastest convergence mode.

Figure 6.

Performance of the NLMS-MF and NLMS algorithms (using different normalized step-size parameters), for the identification of the global impulse response $\mathbf{g}$. The input signals are AR(1) processes, $N=4$, and $L=2048$.

The same conclusion applies when the order $N$ is increased. Figure 7 shows the case when $N=5$ and, hence, the length of the global impulse response is $L=4096$. Again, the NLMS-MF algorithm achieves a significantly better convergence rate and tracking with respect to the conventional NLMS algorithm.

Figure 7.

Performance of the NLMS-MF and NLMS algorithms (using different normalized step-size parameters), for the identification of the global impulse response $\mathbf{g}$. The input signals are AR(1) processes, $N=5$, and $L=4096$.

Next, we aim to show the influence of the normalized step-size values on the performance of the proposed algorithm. In Figure 8, the order of the system is $N=4$ and the input signals are AR(1) processes. As it was expected, lower values of the normalized step-sizes improve the misalignment, but at the cost of a slower convergence and tracking.

Figure 8.

Performance of the NLMS-MF algorithm using different normalized step-size parameters (with equal values of ${\alpha }_i,i=1,2,\dots ,N$), for the identification of the global impulse response $\mathbf{g}$. The input signals are AR(1) processes, $N=4$, and $L=2048$.

The last experiment involving the NLMS-MF algorithm aims to show the performance in the case when the normalized step-size parameters ${\alpha }_i$ ($i=1,2,\dots ,N$) take different values for each individual filter. In Figure 9, the value of the normalized step-size parameter for the first filter, of highest length, $L_1$, is ${\alpha }_1=0.25$, whereas the other normalized step-sizes ${\alpha }_j,j=\mathrm{2,3,4}$ are varied. The system order is $N=4$ and the input signals are AR(1) processes. The compromise between convergence rate and misalignment can be seen again. Since the convergence is mostly influenced by the individual filter with the highest length [20], different values of the normalized step-sizes may be used in different scenarios.

Figure 9.

Performance of the NLMS-MF algorithm using different normalized step-size parameters (with different values of ${\alpha }_i,i=1,2,\dots ,N$), for the identification of the global impulse response $\mathbf{g}$. The input signals are AR(1) processes, $N=4$, and $L=2048$.

Due to the important improvement in performance brought by the adaptive tensor-based LMS algorithms, observed through experiments, these algorithms may represent appealing solutions for the identification of long-length separable system impulse responses.

6. Conclusions

In this chapter, we have presented a decomposition-based approach for dealing with the identification of high-dimension MISO systems. Unlike the conventional method, which is based on the identification of the global system impulse response, our solution focuses on regarding the system as an $N$-order tensor and thus estimating $N$ shorter filters. At the end, the solutions are combined (“tensorized” together) into the original high-dimension tensor. Based on the tensor decomposition technique, an iterative Wiener filter was proposed, along with a family of LMS-based algorithms suitable for the identification of such systems, also called multilinear systems. In addition to the lower computational complexity, the proposed solutions achieve superior performance as compared to their conventional counterparts from the point of view of convergence rate, tracking capability, and steady-state misadjustment. Experiments have also shown the performance improvement of the proposed adaptive algorithms for multilinear forms in long-length system identification in different scenarios, even for highly correlated input signals.

Acknowledgments

This work was supported by a grant from the Romanian Ministry of Education and Research, CNCS–UEFISCDI, Project Number PN-III-P1-1.1-PD-2019-0340, within PNCDI III.