Open AccessFeature PaperArticle

A Recursive Least-Squares Algorithm for the Identification of Trilinear Forms

Camelia Elisei-Iliescu

¹,

Laura-Maria Dogariu

¹,

Constantin Paleologu

^1,*

Jacob Benesty

Andrei-Alexandru Enescu

¹ and

Silviu Ciochină

Department of Telecommunications, University Politehnica of Bucharest, 1-3, Iuliu Maniu Blvd., 061071 Bucharest, Romania

INRS-EMT, University of Quebec, Montreal, QC H5A 1K6, Canada

Author to whom correspondence should be addressed.

Algorithms 2020, 13(6), 135; https://doi.org/10.3390/a13060135

Submission received: 23 March 2020 / Revised: 26 May 2020 / Accepted: 27 May 2020 / Published: 1 June 2020

Download

Browse Figures

Versions Notes

Abstract

High-dimensional system identification problems can be efficiently addressed based on tensor decompositions and modelling. In this paper, we design a recursive least-squares (RLS) algorithm tailored for the identification of trilinear forms, namely RLS-TF. In our framework, the trilinear form is related to the decomposition of a third-order tensor (of rank one). The proposed RLS-TF algorithm acts on the individual components of the global impulse response, thus being efficient in terms of both performance and complexity. Simulation results indicate that the proposed solution outperforms the conventional RLS algorithm (which handles only the global impulse response), but also the previously developed trilinear counterparts based on the least-mean- squares algorithm.

Keywords:

adaptive filters; recursive least-squares (RLS) algorithm; system identification; tensor decomposition; trilinear forms

1. Introduction

Currently, there is an increasing interest in developing methods and algorithms that exploit tensor decompositions and modelling [1,2]. These techniques become of significant importance in many real-world scenarios, e.g., when dealing with large amounts of data, processing multidimensional signals, or solving high-dimensional system identification problems. Many important applications rely on such tensor-based techniques, which can be successfully used in the fields of big data [3], source separation [4], machine learning [5], multiple-input multiple-output (MIMO) communication systems [6], and beamforming [7].

Tensor decompositions and their related applications are frequently addressed based on multilinear signal processing techniques [8,9]. For example, in the context of system identification scenarios, the problems can be formulated in terms of identifying multilinear forms. As particular cases, we can mention the bilinear and trilinear forms, where the decomposition is performed using two and three components, respectively. Since the Wiener filter and adaptive algorithms represent popular methods to address system identification problems, their applicability was also extended to the multilinear framework. Among the recent related works, we can mention the iterative Wiener filter for bilinear forms [10] and the subsequent adaptive filtering methods [11,12,13], together with their extensions to trilinear forms [14,15,16,17].

In this context, the work in [17] provided a system identification framework based on tensor decomposition, which was suitable for the trilinear approach. This work presents the iterative Wiener filter and least-mean-squares (LMS) adaptive algorithms tailored for trilinear forms. Among those, the normalized LMS (NLMS) algorithm for trilinear forms, namely NLMS-TF, represents a practical solution in terms of both performance and complexity. Nevertheless, it is known that the recursive least-squares (RLS) algorithm [18] could provide improved performance as compared to the LMS-based algorithms, especially in terms of the convergence rate. This represents the motivation behind the solution proposed in this paper, which targets the development of the RLS algorithm for trilinear forms, namely RLS-TF. Therefore, the goal of this work is mainly twofold. First, we aim to design a faster converging algorithm as compared to the recently developed NLMS-TF. Second, we intend to outline the performance features of the RLS-TF algorithm as compared to the conventional RLS benchmark.

Besides the conventional adaptive algorithms (which usually act as supervised methods), we should also note that unsupervised algorithms can be used in conjunction with tensor-based approaches, e.g., [19,20]. In this context, we can mention the single-channel blind source separation problem, which targets the identification of individual source signals from a single mixture recording. Such an approach can be found in [19], where the parameters of a so-called “imitated-stereo” mixture model (i.e., one real and the other virtual microphones, which results in an artificial mixing system of dual channels) were found by applying tensor estimation and exploiting sparsity features. The method proposed in [20] also exploits sparsity (using the Gibbs distribution) for multichannel source separation, by solving the underdetermined convolutive mixture separation, while considering the reverberations of the surrounding environment. Another appealing field where such tensor-based techniques can be applied is video processing. For example, unsupervised algorithms can be used for detecting anomalies in a video sequence, e.g., detecting micro defects while employing a thermography imaging system [21].

Currently, in many real-world applications related to MIMO systems, speech processing, and image/video processing, the received signals are tensors or can be grouped in a tensorial form. Consequently, as outlined before, utilizing estimation techniques from tensor algebra is beneficial. Moreover, in most of these applications, the underlying parameters to be estimated are sparse, so that specific features could be exploited, thus bringing additional advantages, especially for systems of large dimensions. The tensor-based RLS algorithm proposed in this paper represents an improved solution (in terms of convergence rate) as compared to the existing solution based on the NLMS algorithm [17], but also an efficient version (in terms of computational complexity) as compared to the conventional RLS algorithm [18].

The rest of the paper is organized as follows. Section 2 provides a background on third-order tensors, which represents the framework of trilinear forms. The proposed RLS-TF algorithm is developed in Section 3. Simulation results are provided in Section 4, outlining the main performance features of the designed solution. Finally, Section 5 concludes this paper and discusses several perspectives for future works.

2. Third-Order Tensors

In this section, we provide a brief summary on tensors, outlining the main definitions and operations, while also establishing the notation. A tensor can be defined as a multidimensional array of data [22,23]. For example, a matrix and a vector can be referred to as second- and first-order tensors, respectively. Tensors of order three or higher are called higher order tensors. In the following, the notation used for a tensor, a matrix, a vector, and a scalar is

A

A

a

, and a, respectively. In this work, we are focusing only on real-valued third-order tensors, i.e.,

A \in R^{L_{1} \times L_{2} \times L_{3}}

, so that the array dimension is

L_{1} \times L_{2} \times L_{3}

The entries of a tensor can be referred to by using multiple indices. In our case, for a third-order tensor, the first and second indices

l_{1}

(with

l_{1} = 1, 2, \dots, L_{1}

) and

l_{2}

(with

l_{2} = 1, 2, \dots, L_{2}

) correspond to the row and column, respectively (like in a matrix); in addition, the third index

l_{3}

(with

l_{3} = 1, 2, \dots, L_{3}

) corresponds to the tube and describes its depth. Consequently, these three indices describe the three different modes. In terms of notation, the entries of the different order tensors are denoted by

{(A)}_{l_{1} l_{2} l_{3}} = a_{l_{1} l_{2} l_{3}}

{(A)}_{l_{1} l_{2}} = a_{l_{1} l_{2}}

, and

{(a)}_{l_{1}} = a_{l_{1}}

In the case of a matrix, the vectorization operation leads to a vector that contains the concatenated columns of the original matrix. In the case of a third-order tensor, the so-called matricization operation concatenates the “slices” of the tensor and produces a large matrix. Of course, the result depends on which index, all the elements of which are considered first. Thus, the matricization can be performed along three different modes [8,9]. For example, considering the mode row and then varying the columns and the tubes, we obtain:

\begin{matrix} A_{[1]} & = A_{:, 1 : L_{2}, 1 : L_{3}} \\ = [\begin{matrix} A_{: : 1} & \dots & A_{: : L_{3}} \end{matrix}], \end{matrix}

where

A_{[1]} \in R^{L_{1} \times L_{2} L_{3}}

and the matrices

A_{: : l_{3}} \in R^{L_{1} \times L_{2}}

with

l_{3} = 1, 2 \dots, L_{3}

denote the frontal slices. Similarly, we can take the mode column and then vary the lines and the tubes, which results in

A_{[2]} = A_{1 : L_{1}, :, 1 : L_{3}}

, with

A_{[2]} \in R^{L_{2} \times L_{1} L_{3}}

. Finally, we can take the mode tube and then vary the rows and the columns, in order to obtain

A_{[3]} = A_{1 : L_{1}, 1 : L_{2}, :}

, where

A_{[3]} \in R^{L_{3} \times L_{1} L_{2}}

. The ranks of

A_{[1]}

A_{[2]}

, and

A_{[3]}

represent the mode-1, mode-2, and mode-3 ranks of the tensor, respectively. Furthermore, the vectorization of a tensor (e.g., following mode-1) is

\begin{matrix} vec (A) & = vec (A_{[1]}) \\ = [\begin{matrix} vec (A_{: : 1}) \\ ⋮ \\ vec (A_{: : L_{3}}) \end{matrix}] \in R^{L_{1} L_{2} L_{3}} . \end{matrix}

Nevertheless, there are some fundamental differences between the rank of a matrix

A \in R^{L_{1} \times L_{2}}

and the rank of a tensor

A \in R^{L_{1} \times L_{2} \times L_{3}}

. For example, the rank of

A

can never be larger than

\min {L_{1}, L_{2}}

, while the rank of

A

can be greater than

\min {L_{1}, L_{2}, L_{3}}

. A rank-1 tensor (of dimension

L_{1} \times L_{2} \times L_{3}

) is defined as:

\begin{matrix} B & = b_{1} \circ b_{2} \circ b_{3}, \end{matrix}

(1)

where

b_{1}

b_{2}

, and

b_{3}

are vectors of lengths

L_{1}

L_{2}

, and

L_{3}

, respectively, ∘ is the vector outer product, and the elements of

B

are given by

{(B)}_{l_{1} l_{2} l_{3}} = b_{1, l_{1}} b_{2, l_{2}} b_{3, l_{3}}

, where

b_{i, l_{i}}

are the elements of the vector

b_{i}

, with

i = 1, 2, 3

and

l_{i} = 1, 2, \dots, L_{i}

. In this case, it can be easily verified that:

\begin{matrix} vec (B) & = b_{3} \otimes b_{2} \otimes b_{1}, \end{matrix}

(2)

where ⊗ denotes the Kronecker product [24]. In this context, the rank of a tensor

A

, denoted

rank (A)

, is defined as the minimum number of rank-1 tensors that generate

A

as their sum. For example, if:

\begin{matrix} A & = \sum_{r = 1}^{R} a_{1 r} \circ a_{2 r} \circ a_{3 r}, \end{matrix}

(3)

where

a_{1 r}

a_{2 r}

, and

a_{3 r}

are vectors of lengths

L_{1}

L_{2}

, and

L_{3}

, respectively, then

rank (A) = R

when R is minimal and (3) is called the canonical polyadic decomposition (CPD) of

A

Another important operation is the multiplication of a tensor with a matrix [1,2], which can also be defined in different ways. For example, the mode-1 product between the tensor

A \in R^{L_{1} \times L_{2} \times L_{3}}

and the matrix

M_{1} \in R^{M_{1} \times L_{1}}

gives the tensor:

\begin{matrix} U & = A \times_{1} M_{1}, U \in R^{M_{1} \times L_{2} \times L_{3}}, \end{matrix}

(4)

whose entries are

u_{m_{1} l_{2} l_{3}} = \sum_{l_{1} = 1}^{L_{1}} a_{l_{1} l_{2} l_{3}} m_{m_{1} l_{1}}

(for

m_{1} = 1, 2, \dots, M_{1}

) and

U_{[1]} = M_{1} A_{[1]}

. Similarly, the mode-2 product between the same tensor

A

and the matrix

M_{2} \in R^{M_{2} \times L_{2}}

leads to the tensor:

\begin{matrix} \underset{̲}{U} & = A \times_{2} M_{2}, \underset{̲}{U} \in R^{L_{1} \times M_{2} \times L_{3}}, \end{matrix}

(5)

with the entries

{\underset{̲}{u}}_{l_{1} m_{2} l_{3}} = \sum_{l_{2} = 1}^{L_{2}} a_{l_{1} l_{2} l_{3}} m_{m_{2} l_{2}}

(for

m_{2} = 1, 2, \dots, M_{2}

) and

{\underset{̲}{U}}_{[2]} = M_{2} A_{[2]}

. Finally, the mode-3 product between the tensor

A

and the matrix

M_{3} \in R^{M_{3} \times L_{3}}

results in the tensor:

\begin{matrix} \bar{U} & = A \times_{3} M_{3}, \bar{U} \in R^{L_{1} \times L_{2} \times M_{3}}, \end{matrix}

(6)

having the entries

{\bar{u}}_{l_{1} l_{2} m_{3}} = \sum_{l_{3} = 1}^{L_{3}} a_{l_{1} l_{2} l_{3}} m_{m_{3} l_{3}}

(for

m_{3} = 1, 2, \dots, M_{3}

), while

{\bar{U}}_{[3]} = M_{3} A_{[3]}

. Clearly, we can multiply the tensor

A

with the three previously defined matrices

M_{1}

M_{2}

, and

M_{3}

. In this case, we get the tensor:

\begin{matrix} T & = A \times_{1} M_{1} \times_{2} M_{2} \times_{3} M_{3} \\ = A \times_{1} M_{1} \times_{3} M_{3} \times_{2} M_{2} \\ = A \times_{2} M_{2} \times_{1} M_{1} \times_{3} M_{3} \\ = A \times_{2} M_{2} \times_{3} M_{3} \times_{1} M_{1} \\ = A \times_{3} M_{3} \times_{1} M_{1} \times_{2} M_{2} \\ = A \times_{3} M_{3} \times_{2} M_{2} \times_{1} M_{1}, \end{matrix}

(7)

where

T \in R^{M_{1} \times M_{2} \times M_{3}}

. As a consequence, considering

b_{1}

b_{2}

, and

b_{3}

three vectors of lengths

L_{1}

L_{2}

, and

L_{3}

, respectively, the multiplication of the tensor

A

with the transposed vectors gives the scalar:

\begin{matrix} c & = A \times_{1} b_{1}^{T} \times_{2} b_{2}^{T} \times_{3} b_{3}^{T} \\ = \sum_{l_{1} = 1}^{L_{1}} \sum_{l_{2} = 1}^{L_{2}} \sum_{l_{3} = 1}^{L_{3}} a_{l_{1} l_{2} l_{3}} b_{1 l_{1}} b_{2 l_{2}} b_{3 l_{3}}, \end{matrix}

(8)

where

^{T}

is the transpose operator. It is easy to check that (8) is trilinear with respect to

b_{1}

b_{2}

, and

b_{3}

At this point, we can also define the inner product between two tensors

A

and

B

of the same dimension (

L_{1} \times L_{2} \times L_{3}

), which is:

\begin{matrix} 〈 A, B 〉 & = \sum_{l_{1} = 1}^{L_{1}} \sum_{l_{2} = 1}^{L_{2}} \sum_{l_{3} = 1}^{L_{3}} a_{l_{1} l_{2} l_{3}} b_{l_{1} l_{2} l_{3}} \\ = {vec}^{T} (B) vec (A) . \end{matrix}

(9)

Therefore, Expression (8) can also be written in a more convenient way, i.e.,

\begin{matrix} c & = 〈 A, B 〉 \\ = {vec}^{T} (B) vec (A) \\ = {vec}^{T} (b_{1} \circ b_{2} \circ b_{3}) vec (A) \\ = {(b_{3} \otimes b_{2} \otimes b_{1})}^{T} vec (A), \end{matrix}

(10)

where

B = b_{1} \circ b_{2} \circ b_{3}

(see (1)). Moreover, if

A

is also a rank-1 tensor, i.e.,

A = a_{1} \circ a_{2} \circ a_{3}

, where the vectors

a_{i}

have the lengths

L_{i}

(with

i = 1, 2, 3

), then:

\begin{matrix} 〈 A, B 〉 & = b_{1}^{T} a_{1} \times b_{2}^{T} a_{2} \times b_{3}^{T} a_{3} . \end{matrix}

(11)

Furthermore, it is easy to check that:

\begin{matrix} B \times_{1} M_{1} & = M_{1} b_{1} \circ b_{2} \circ b_{3}, \\ B \times_{2} M_{2} & = b_{1} \circ M_{2} b_{2} \circ b_{3}, \\ B \times_{3} M_{3} & = b_{1} \circ b_{2} \circ M_{3} b_{3}, \end{matrix}

where the matrices

M_{1}

M_{2}

, and

M_{3}

were previously defined (related to (4)–(6)).

The short background on tensors provided before and the main related operations (e.g., matricization, vectorization, rank, and different types of product) aim to facilitate the development that follows in Section 3. It is also important to outline that the trilinear forms result in the context of the decomposition of third-order tensors. Extension to higher order tensors and multilinear forms could be straightforward when dealing with rank-1 tensors. Otherwise, in the general case, decomposing higher rank higher order tensors (see (3)) raises additional difficulties, as will be briefly pointed out at the end of Section 5.

3. RLS Algorithm for Trilinear Forms

In the following, for the sake of consistency with the development of the NLMS-TF algorithm, we will keep the framework and notation from [17]. Therefore, let us consider the output of a multiple-input single-output (MISO) system (with real-valued data) at the discrete-time index n defined as:

\begin{matrix} y (n) & = X (n) \times_{1} h_{1}^{T} \times_{2} h_{2}^{T} \times_{3} h_{3}^{T} \\ = \sum_{l_{1} = 1}^{L_{1}} \sum_{l_{2} = 1}^{L_{2}} \sum_{l_{3} = 1}^{L_{3}} x_{l_{1} l_{2} l_{3}} (n) h_{1, l_{1}} h_{2, l_{2}} h_{3, l_{3}}, \end{matrix}

(12)

where the tensorial form

X (n) \in R^{L_{1} \times L_{2} \times L_{3}}

groups the input signals, with:

\begin{matrix} {[X (n)]}_{l_{1} l_{2} l_{3}} & = x_{l_{1} l_{2} l_{3}} (n), l_{i} = 1, 2, \dots, L_{i}, i = 1, 2, 3, \end{matrix}

and the vectors

h_{i}, i = 1, 2, 3

(of lengths

L_{1}

L_{2}

, and

L_{3}

, respectively) define the three impulse responses, i.e.,

\begin{matrix} h_{i} = {[\begin{matrix} h_{i, 1} & h_{i, 2} & \dots & h_{i, L_{i}} \end{matrix}]}^{T}, i = 1, 2, 3 . \end{matrix}

As we can notice,

y (n)

represents a trilinear form (see (12) as compared to (8)), because it is a linear function of each of the vectors

h_{i}, i = 1, 2, 3

, if the other two are fixed.

Next, we can also introduce a rank-1 tensor of dimension

L_{1} \times L_{2} \times L_{3}

, using the three impulse responses of the MISO system:

\begin{matrix} H & = h_{1} \circ h_{2} \circ h_{3}, \end{matrix}

(13)

whose elements are:

\begin{matrix} {(H)}_{l_{1} l_{2} l_{3}} & = h_{1, l_{1}} h_{2, l_{2}} h_{3, l_{3}}, l_{i} = 1, 2, \dots, L_{i}, i = 1, 2, 3 . \end{matrix}

Consequently, the output signal results in:

\begin{matrix} y (n) = {vec}^{T} (H) vec [X (n)], \end{matrix}

(14)

where:

\begin{matrix} vec (H) & = [\begin{matrix} vec (H_{: : 1}) \\ ⋮ \\ vec (H_{: : L_{3}}) \end{matrix}], \end{matrix}

(15)

\begin{matrix} vec [X (n)] & = [\begin{matrix} vec [X_{: : 1} (n)] \\ ⋮ \\ vec [X_{: : L_{3}} (n)] \end{matrix}], \end{matrix}

(16)

with

H_{: : l_{3}}

and

X_{: : l_{3}} (n)

(

l_{3} = 1, 2 \dots, L_{3}

) being the frontal slices of

H

and

X (n)

, respectively. At this point, we can introduce the notation:

\begin{matrix} h & ≜ vec (H) = h_{3} \otimes h_{2} \otimes h_{1}, \end{matrix}

(17)

\begin{matrix} x (n) & ≜ vec [X (n)], \end{matrix}

(18)

where

h

and

x (n)

are two long vectors, each of them having

L_{1} L_{2} L_{3}

elements. Thus, the output signal can also be expressed as:

\begin{matrix} y (n) = h^{T} x (n) . \end{matrix}

(19)

The main goal is to estimate the output of the MISO system, which is usually corrupted by an additive noise. Hence, the reference signal results in:

\begin{matrix} d (n) & = y (n) + w (n) \\ = h^{T} x (n) + w (n), \end{matrix}

(20)

where

w (n)

is a zero-mean additive noise, which is uncorrelated with the input signals. Alternatively, we could estimate the global impulse response

h

, using an adaptive filter

\hat{h} (n)

(of length

L_{1} L_{2} L_{3}

). At this point, we may also define the error signal:

\begin{matrix} e (n) & = d (n) - \hat{y} (n) \\ = d (n) - {\hat{h}}^{T} (n - 1) x (n), \end{matrix}

(21)

which represents the difference between the reference signal and the estimated signal,

\hat{y} (n)

Based on (17), we can notice that the global impulse response

h

(of length

L_{1} L_{2} L_{3}

) results based on a combination of the shorter impulse responses

h_{i}, i = 1, 2, 3

, of lengths

L_{1}

L_{2}

, and

L_{3}

, respectively. Consequently, the estimated impulse response can also be decomposed as:

\begin{matrix} \hat{h} (n) = {\hat{h}}_{3} (n) \otimes {\hat{h}}_{2} (n) \otimes {\hat{h}}_{1} (n), \end{matrix}

(22)

where

{\hat{h}}_{i} (n), i = 1, 2, 3

are three adaptive filters (with

L_{1}

L_{2}

, and

L_{3}

coefficients, respectively), which aim to model the individual impulse responses

h_{i}, i = 1, 2, 3

. Nevertheless, we should notice that there is no unique solution related to the decomposition in (22). It is obvious that, for any constants

η_{1}, η_{2}

, and

η_{3}

, with

η_{1} η_{2} η_{3} = 1

, we have:

\begin{matrix} h & = h_{3} \otimes h_{2} \otimes h_{1} \\ = η_{3} h_{3} \otimes η_{2} h_{2} \otimes η_{1} h_{1} . \end{matrix}

(23)

Hence,

η_{i} h_{i}, i = 1, 2, 3

also represent a set of solutions. However, the global impulse response

h

is identified with no scaling ambiguity.

Following the decomposition from (22), we can easily verify that:

\begin{matrix} \hat{h} (n) & = {\hat{H}}_{32} (n) {\hat{h}}_{1} (n) \end{matrix}

(24)

\begin{matrix} = {\hat{H}}_{31} (n) {\hat{h}}_{2} (n) \end{matrix}

(25)

\begin{matrix} = {\hat{H}}_{21} (n) {\hat{h}}_{3} (n), \end{matrix}

(26)

where:

\begin{matrix} {\hat{H}}_{32} (n) & = {\hat{h}}_{3} (n) \otimes {\hat{h}}_{2} (n) \otimes I_{L_{1}}, \end{matrix}

(27)

\begin{matrix} {\hat{H}}_{31} (n) & = {\hat{h}}_{3} (n) \otimes I_{L_{2}} \otimes {\hat{h}}_{1} (n), \end{matrix}

(28)

\begin{matrix} {\hat{H}}_{21} (n) & = I_{L_{3}} \otimes {\hat{h}}_{2} (n) \otimes {\hat{h}}_{1} (n), \end{matrix}

(29)

and

I_{L_{i}}, i = 1, 2, 3

are the identity matrices of sizes

L_{1} \times L_{1}

L_{2} \times L_{2}

, and

L_{3} \times L_{3}

, respectively. Moreover, introducing the notation:

\begin{matrix} x_{32} (n) & = {\hat{H}}_{32}^{T} (n - 1) x (n), \end{matrix}

(30)

\begin{matrix} x_{31} (n) & = {\hat{H}}_{31}^{T} (n - 1) x (n), \end{matrix}

(31)

\begin{matrix} x_{21} (n) & = {\hat{H}}_{21}^{T} (n - 1) x (n), \end{matrix}

(32)

the error signal from (21) can be expressed in three equivalent ways as:

\begin{matrix} e (n) & = d (n) - {\hat{h}}_{1}^{T} (n - 1) x_{32} (n), \end{matrix}

(33)

\begin{matrix} = d (n) - {\hat{h}}_{2}^{T} (n - 1) x_{31} (n), \end{matrix}

(34)

\begin{matrix} = d (n) - {\hat{h}}_{3}^{T} (n - 1) x_{21} (n) . \end{matrix}

(35)

Based on the least-squares error criterion [18] applied in the context of (20) and (21), the conventional RLS algorithm is derived from:

\begin{matrix} J [\hat{h} (n)] & = \sum_{i = 1}^{n} λ^{n - i} {[d (i) - {\hat{h}}^{T} (n) x (i)]}^{2}, \end{matrix}

(36)

where

λ \leq 1

is a positive constant known as the forgetting factor. On the other hand, based on (24)–(26), the cost function from (36) can be expressed in three different ways, targeting the optimization of the individual components, i.e.,

\begin{matrix} J_{{\hat{h}}_{3}, {\hat{h}}_{2}} [{\hat{h}}_{1} (n)] & = \sum_{i = 1}^{n} λ_{1}^{n - i} {[d (i) - {\hat{h}}_{1}^{T} (n) x_{32} (i)]}^{2}, \end{matrix}

(37)

\begin{matrix} J_{{\hat{h}}_{3}, {\hat{h}}_{1}} [{\hat{h}}_{2} (n)] & = \sum_{i = 1}^{n} λ_{2}^{n - i} {[d (i) - {\hat{h}}_{2}^{T} (n) x_{31} (i)]}^{2}, \end{matrix}

(38)

\begin{matrix} J_{{\hat{h}}_{2}, {\hat{h}}_{1}} [{\hat{h}}_{3} (n)] & = \sum_{i = 1}^{n} λ_{3}^{n - i} {[d (i) - {\hat{h}}_{3}^{T} (n) x_{21} (i)]}^{2}, \end{matrix}

(39)

where

λ_{1}

λ_{2}

, and

λ_{3}

are the individual forgetting factors. The previous cost functions suggest a “trilinear” optimization strategy [25], where we assume that two components are fixed during the optimization of the third one. Consequently, based on the minimization of (37)–(39) with respect to

{\hat{h}}_{1} (n)

{\hat{h}}_{2} (n)

, and

{\hat{h}}_{3} (n)

, respectively, the following set of normal equations are obtained:

\begin{matrix} R_{32} (n) {\hat{h}}_{1} (n) & = p_{32} (n), \end{matrix}

(40)

\begin{matrix} R_{31} (n) {\hat{h}}_{2} (n) & = p_{31} (n), \end{matrix}

(41)

\begin{matrix} R_{21} (n) {\hat{h}}_{3} (n) & = p_{21} (n), \end{matrix}

(42)

where:

\begin{matrix} R_{32} (n) & = \sum_{i = 1}^{n} λ_{1}^{n - i} x_{32} (i) x_{32}^{T} (i) \\ = λ_{1} R_{32} (n - 1) + x_{32} (n) x_{32}^{T} (n), \\ p_{32} (n) & = \sum_{i = 1}^{n} λ_{1}^{n - i} x_{32} (i) d (i) \\ = λ_{1} p_{32} (n - 1) + x_{32} (n) d (n), \\ R_{31} (n) & = \sum_{i = 1}^{n} λ_{2}^{n - i} x_{31} (i) x_{31}^{T} (i) \\ = λ_{2} R_{31} (n - 1) + x_{31} (n) x_{31}^{T} (n), \\ p_{31} (n) & = \sum_{i = 1}^{n} λ_{2}^{n - i} x_{31} (i) d (i) \\ = λ_{2} p_{31} (n - 1) + x_{31} (n) d (n), \\ R_{21} (n) & = \sum_{i = 1}^{n} λ_{3}^{n - i} x_{21} (i) x_{21}^{T} (i) \\ = λ_{2} R_{21} (n - 1) + x_{21} (n) x_{21}^{T} (n), \\ p_{21} (n) & = \sum_{i = 1}^{n} λ_{3}^{n - i} x_{21} (i) d (i) \\ = λ_{3} p_{21} (n - 1) + x_{21} (n) d (n) . \end{matrix}

Solving (40)–(42), the updates of the individual filters result in:

\begin{matrix} {\hat{h}}_{1} (n) & = {\hat{h}}_{1} (n - 1) + R_{32}^{- 1} (n) x_{32} (n) e (n) \\ = {\hat{h}}_{1} (n - 1) + k_{32} (n) e (n), \end{matrix}

(43)

\begin{matrix} {\hat{h}}_{2} (n) & = {\hat{h}}_{2} (n - 1) + R_{31}^{- 1} (n) x_{31} (n) e (n) \\ = {\hat{h}}_{2} (n - 1) + k_{31} (n) e (n), \end{matrix}

(44)

\begin{matrix} {\hat{h}}_{3} (n) & = {\hat{h}}_{3} (n - 1) + R_{21}^{- 1} (n) x_{21} (n) e (n) \\ = {\hat{h}}_{3} (n - 1) + k_{21} (n) e (n), \end{matrix}

(45)

where

k_{32} (n) = R_{32}^{- 1} (n) x_{32} (n)

k_{31} (n) = R_{31}^{- 1} (n) x_{31} (n)

, and

k_{21} (n) = R_{21}^{- 1} (n) x_{21} (n)

are the Kalman gain vectors, while the error signal can be computed based on (33). At this point, the main task is to update the inverse of the matrices

R_{32} (n)

R_{31} (n)

, and

R_{21} (n)

efficiently. The solution relies on the matrix inversion lemma [18], which leads to the following updates:

\begin{matrix} R_{32}^{- 1} (n) & = \frac{1}{λ_{1}} [I_{L_{1}} - k_{32} (n) x_{32}^{T} (n)] R_{32}^{- 1} (n - 1), \end{matrix}

(46)

\begin{matrix} R_{31}^{- 1} (n) & = \frac{1}{λ_{2}} [I_{L_{2}} - k_{31} (n) x_{31}^{T} (n)] R_{31}^{- 1} (n - 1), \end{matrix}

(47)

\begin{matrix} R_{21}^{- 1} (n) & = \frac{1}{λ_{3}} [I_{L_{3}} - k_{21} (n) x_{21}^{T} (n)] R_{21}^{- 1} (n - 1) . \end{matrix}

(48)

Therefore, the Kalman gain vectors are evaluated as:

\begin{matrix} k_{32} (n) & = \frac{R_{32}^{- 1} (n - 1) x_{32} (n)}{λ_{1} + x_{32}^{T} (n) R_{32}^{- 1} (n - 1) x_{32} (n)}, \end{matrix}

(49)

\begin{matrix} k_{31} (n) & = \frac{R_{31}^{- 1} (n - 1) x_{31} (n)}{λ_{2} + x_{31}^{T} (n) R_{31}^{- 1} (n - 1) x_{31} (n)}, \end{matrix}

(50)

\begin{matrix} k_{21} (n) & = \frac{R_{21}^{- 1} (n - 1) x_{21} (n)}{λ_{3} + x_{21}^{T} (n) R_{21}^{- 1} (n - 1) x_{21} (n)} . \end{matrix}

(51)

For initialization, we can choose:

\begin{matrix} {\hat{h}}_{1} (0) & = [\begin{matrix} 1 \\ 0_{L_{1} - 1} \end{matrix}], \end{matrix}

(52)

\begin{matrix} {\hat{h}}_{2} (0) & = \frac{1}{L_{2}} 1_{L_{2}}, \end{matrix}

(53)

\begin{matrix} {\hat{h}}_{3} (0) & = \frac{1}{L_{3}} 1_{L_{3}}, \end{matrix}

(54)

where

0_{N}

and

1_{N}

are two vectors of length N, all elements of which are equal to zero and one, respectively.

The proposed RLS algorithm for trilinear forms, namely RLS-TF, is summarized in Table 1. It could also be interpreted as an extension of the RLS-based algorithm tailored for bilinear forms, which was presented in [12]. However, if the MISO system identification problem results based on (12), it is natural to use the RLS-TF algorithm, which is designed for the identification of third-order (rank-1) tensors.

In terms of computational complexity, it can be noticed that the proposed RLS-TF algorithm combines the solutions provided by three RLS-based filters, i.e.,

{\hat{h}}_{1} (n)

(of length

L_{1}

{\hat{h}}_{2} (n)

(of length

L_{2}

), and

{\hat{h}}_{3} (n)

(of length

L_{3}

). Since the complexity of a regular RLS-based algorithm is proportional to the square of the filter length, the overall complexity of the RLS-TF algorithm roughly results in

O (L_{1}^{2} + L_{2}^{2} + L_{3}^{2})

. On the other hand, the system identification problem can also be handled based on the conventional RLS algorithm, which results following (20) and (21), together with the cost function from (36). However, in this case, there is a single adaptive filter

\hat{h} (n)

, of length

L = L_{1} L_{2} L_{3}

, so that the computational complexity is of the order of

O (L^{2})

. This could be much more computationally expensive as compared to the proposed RLS-TF algorithm.

Basically, the RLS-TF algorithm “transforms” a large system identification problem of length

L = L_{1} L_{2} L_{3}

into three “smaller” problems of lengths

L_{1}

L_{2}

, and

L_{3}

, respectively, with advantages in terms of both performance (as will be shown in simulations) and complexity. As outlined before, the proposed RLS-TF algorithm combines the solutions provided by three adaptive filters of lengths

L_{1}

L_{2}

, and

L_{3}

, respectively, while the conventional RLS algorithm deals with a single filter of length

L = L_{1} L_{2} L_{3}

, which is usually much longer. Since the length of the filter highly influences the main performance criteria, i.e., convergence rate and misadjustment [18], the proposed algorithm is able to outperform the conventional one in terms of both criteria. In other words, the shorter the length, the faster the convergence and the lower the misadjustment. This expected behaviour will be supported by the simulation results provided in the next section.

Finally, we should observe that there are some extra operations specific to the RLS-TF algorithm. For example, the “input” signals

x_{32} (n)

x_{31} (n)

, and

x_{21} (n)

result based on (30)–(32), which rely on (27)–(29). Furthermore, the global impulse response (if required within the application) can be evaluated based on (22). These operations require Kronecker products, but the related computational complexity is moderate, i.e., of the order of

O (L_{1} L_{2} L_{3}) = O (L)

The detailed computational complexities of the proposed RLS-TF algorithm and other benchmark algorithms (i.e., the conventional RLS and NLMS algorithms) are summarized in Table 2. For a better visualization, the computational complexities are also illustrated in Figure 1, in terms of the number of multiplications and additions (per iteration), for different values of

L_{1}

; the other lengths are fixed to

L_{2} = 8

and

L_{3} = 4

(similar to the experimental setup from Section 4). As we can notice, the computational complexity of the conventional RLS algorithm was significantly greater, while the computational amount of the proposed RLS-TF algorithm was closer to the conventional NLMS algorithm, especially for higher lengths.

4. Simulation Results

Simulations were performed in the framework of a tensor-based system identification problem, which resulted following the MISO model defined by (12) and (20) and was similar to the setup used in [17]. The input signals that form the third-order tensor

X (n)

are AR(1) processes; they are generated by filtering white Gaussian noises through a first-order system

1 / (1 - 0.9 z^{- 1})

. The additive noise

w (n)

is white and Gaussian; its variance was set to

σ_{w}^{2} = 0.01

The third-order system used in the simulations and its components (

h_{1}

h_{2}

, and

h_{3}

) are depicted in Figure 2. First, the component

h_{1}

is an impulse response from the G168 Recommendation [26], of length

L_{1} = 64

; it is provided in Figure 2a. Second, in Figure 2b, the component

h_{2}

is a random impulse response (with Gaussian distribution) of length

L_{2} = 8

. Third, the coefficients of the last component, i.e., the impulse response

h_{3}

, are depicted in Figure 2c; those were evaluated as

h_{3, l_{3}} = {0.5}^{l_{3} - 1}, l_{3} = 1, 2, \dots, L_{3}

, using the length

L_{3} = 4

. Therefore, the global impulse response from Figure 2d resulted as

h = h_{3} \otimes h_{2} \otimes h_{1}

, and its length was

L = L_{1} L_{2} L_{3} = 2048

. This global impulse response resembled a channel with echoes, e.g., like an acoustic echo path [27]. Finally, the third-order (rank-1) tensor

H

of dimension

L_{1} \times L_{2} \times L_{3}

could be formed according to (13). In order to test the tracking capabilities of the algorithms, an abrupt change of the system was introduced in the middle of each experiment, by changing the sign of the coefficients of each impulse response.

As shown in Section 3, the proposed RLS-TF algorithm was designed to estimate the individual components of the global system. However, we could identify

h_{1}

h_{2}

, and

h_{3}

up to some scaling factors, as explained using (23). Therefore, to evaluate the identification of these individual impulse responses, a proper performance measure is the normalized projection misalignment (NPM) [28]:

\begin{matrix} NPM [h_{1}, {\hat{h}}_{1} (n)] & = 1 - {[\frac{h_{1}^{T} {\hat{h}}_{1} (n)}{{∥h_{1}∥}_{2} {∥{\hat{h}}_{1} (n)∥}_{2}}]}^{2}, \end{matrix}

(55)

\begin{matrix} NPM [h_{2}, {\hat{h}}_{2} (n)] & = 1 - {[\frac{h_{2}^{T} {\hat{h}}_{2} (n)}{{∥h_{2}∥}_{2} {∥{\hat{h}}_{2} (n)∥}_{2}}]}^{2}, \end{matrix}

(56)

\begin{matrix} NPM [h_{3}, {\hat{h}}_{3} (n)] & = 1 - {[\frac{h_{3}^{T} {\hat{h}}_{3} (n)}{{∥h_{3}∥}_{2} {∥{\hat{h}}_{3} (n)∥}_{2}}]}^{2}, \end{matrix}

(57)

where

{∥\cdot∥}_{2}

denotes the Euclidean norm. On the other hand, the global impulse response

h

results without any scaling ambiguity. Consequently, we can use a performance measure based on the regular normalized misalignment (NM):

\begin{matrix} NM [h, \hat{h} (n)] & = \frac{{∥h - \hat{h} (n)∥}_{2}^{2}}{{∥h∥}_{2}^{2}}, \end{matrix}

(58)

which is also equivalent to

{∥H - \hat{H} (n)∥}_{F}^{2} / {∥H∥}_{F}^{2}

, where

\hat{H} (n) = {\hat{h}}_{1} (n) \circ {\hat{h}}_{2} (n) \circ {\hat{h}}_{3} (n)

and

{∥\cdot∥}_{F}

denotes the Frobenius norm (the Frobenius norm of a third-order tensor

A

is defined as

{∥A∥}_{F} = \sqrt{〈 A, A 〉} = {∥A_{[1]}∥}_{F} = {∥A_{[2]}∥}_{F} = {∥A_{[3]}∥}_{F}

The simulation results should provide answers to several important questions, as follows. (i) What is the influence of the forgetting factors on the performance of the proposed RLS-TF algorithm? (ii) What are the advantages of the RLS-TF algorithm over the previously developed NLMS-TF counterpart [17]? (iii) What are the advantages of the RLS-TF algorithm over the conventional RLS benchmark? The following three experiments are designed to address these issues.

In the first experiment, the performance of the proposed RLS-TF algorithm was analysed with respect to its main parameters, i.e., the forgetting factors

λ_{1}

λ_{2}

, and

λ_{3}

. In the case of an RLS-based algorithm, the value of the forgetting factor is usually related to the filter length, following a well-known rule of thumb, as shown in Table 1 (see “Initialization”). In our case, the forgetting factors of the RLS-TF algorithm were set to

λ_{1} = 1 - 1 / (K L_{1})

λ_{2} = 1 - 1 / (K L_{2})

, and

λ_{3} = 1 - 1 / (K L_{3})

. As we can notice, the value of each forgetting factor depended on the length of its related filter (i.e.,

L_{1}

L_{2}

, or

L_{3}

), but also on the constant K. This tuning parameter could be used to adjust the values of the forgetting factors, as indicated in Figure 3 and Figure 4. Clearly, a higher value of K would result in a higher value of the forgetting factor (i.e., closer to one). We could expect that a higher value of the forgetting factor would reduce the misalignment, but slowing down the convergence/tracking [29]. On the other hand, reducing the forgetting factor improves the convergence/tracking, but increasing the misalignment. This behaviour was supported by the results depicted in Figure 3 and Figure 4, in terms of NPM and NM, respectively. As we can notice, the value

K = 20

lead to a good compromise between the performance criteria, so that it would be used in the following experiments.

Next, we compare the performance of the proposed RLS-TF algorithm with its previously developed counterpart based on the NLMS algorithm, i.e., NLMS-TF [17]. The overall performance of this algorithm is mainly controlled by its normalized step-sizes, which are positive constants smaller than one. Using notation similar to that involved in [17], we set equal values for these parameters, i.e.,

α_{1} = α_{2} = α_{3} = α

. In the case of the NLMS-TF algorithm, the fastest convergence mode was obtained when

α_{1} + α_{2} + α_{3} \approx 1

, so that we could use

α = 0.33

. Smaller values of the normalized step-sizes (e.g.,

α = 0.1

) reduced the convergence/tracking, but led to a lower misalignment. As shown in Figure 5 and Figure 6 (in terms of NPM and NM, respectively), the RLS-TF algorithm clearly outperformed the NLMS-TF counterpart, achieving a faster convergence rate and tracking, together with low misalignment.

Finally, in the last experiment, we investigated the comparison between the RLS-TF solution and the conventional RLS algorithm. As explained in the last part of Section 3 (related to the computational complexity), the conventional RLS algorithm could also be used for the identification of the global impulse response of length L, based on the cost function from (36). This algorithm uses a single forgetting factor, which can also be set as

λ = 1 - 1 / (K L)

, where K is the same tuning constant. The influence of the value of

λ

on the performance of the algorithm is also related to the well-known compromise between low misalignment and fast tracking. In the experiment reported in Figure 7, the conventional RLS algorithm uses two values of the forgetting factor, which were set by varying the tuning constant to

K = 1

and

K = 20

. As we can notice, even if the largest value of the forgetting factor (obtained for

K = 20

) led to a lower misalignment, the tracking capability of the conventional RLS algorithm was significantly reduced. Clearly, the tracking was improved when using a smaller forgetting factor (corresponding to

K = 1

), but the misalignment of the conventional RLS algorithm was much higher in this case. On the other hand, the RLS-TF algorithm outperformed by far its conventional counterpart, in terms of both performance criteria. Moreover, the complexity of the conventional RLS algorithm, i.e.,

O (L^{2})

, was prohibitive for practical implementations, due to the long length of the global filter (

L = 2048

). On the other hand, the RLS-TF algorithm worked on the individual components and combined the solutions of three shorter filters, of lengths

L_{1}

L_{2}

, and

L_{3}

(with

L_{1} L_{2} L_{3} = L

); thus, it was much more computationally efficient. As a trivial example related to the last experiment given in Figure 7, we could mention that the simulation time (using MATLAB) of the RLS-TF algorithm was less than one minute, while the conventional RLS algorithm took hours to reach the final result.

5. Conclusions and Future Works

In this paper, we explored the identification of trilinear forms using the RLS algorithm. The approach was developed in the framework of an MISO system identification problem, based on the decomposition and modelling of third-order tensors. The resulting RLS-TF algorithm was tailored for the identification of such trilinear forms in a more efficient way as compared to the conventional RLS algorithm. Moreover, the proposed RLS-TF algorithm outperformed its previously developed NLMS-TF counterpart in terms of the main performance criteria, providing a faster convergence rate and tracking, together with low misalignment.

The ideas presented in this paper could be further exploited in an extended framework, aiming to identify more general forms of global impulse responses, which cannot be decomposed as rank-1 tensors. Several preliminary results can be found in [30,31,32], but they are applicable only in the bilinear context (i.e., second-order tensors). The extension to the trilinear case represents a very challenging problem, since finding (and approximating) the rank of a higher order tensor is a much more sensitive task. Furthermore, it would be interesting to extend other versions of the NLMS and RLS algorithms (e.g., based on variable step-sizes [33] and variable forgetting factors [34], respectively) to the trilinear forms.

Finally, it would be useful to evaluate how the proposed algorithm could benefit (in terms of implementation) from the current technology of the tensor processing unit (TPU) and the TensorFlow software [35]. This aspect could bring additional advantages, especially in the framework of specific applications related to machine learning/vision, neural networks, and artificial intelligence.

Author Contributions

Conceptualization, C.E.-I.; methodology, L.-M.D.; investigation, C.P.; validation, J.B.; software, A.-A.E.; formal analysis, S.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by UEFISCDI grant number PN-III-P1-1.1-PD-2019-0340.

Conflicts of Interest

The authors declare no conflict of interest.

References

Comon, P. Tensors: A brief introduction. IEEE Signal Process. Mag. 2014, 31, 44–53. [Google Scholar] [CrossRef] [Green Version]
Cichocki, A.; Mandic, D.; De Lathauwer, L.; Zhou, G.; Zhao, Q.; Caiafa, C.; Phan, H.A. Tensor decompositions for signal processing applications: From two-way to multiway component analysis. IEEE Signal Process. Mag. 2015, 32, 145–163. [Google Scholar] [CrossRef] [Green Version]
Vervliet, N.; Debals, O.; Sorber, L.; De Lathauwer, L. Breaking the curse of dimensionality using decompositions of incomplete tensors: Tensor-based scientific computing in big data analysis. IEEE Signal Process. Mag. 2014, 31, 71–79. [Google Scholar] [CrossRef]
Boussé, M.; Debals, O.; De Lathauwer, L. A tensor-based method for large-scale blind source separation using segmentation. IEEE Trans. Signal Process. 2017, 65, 346–358. [Google Scholar] [CrossRef]
Sidiropoulos, N.; De Lathauwer, L.; Fu, X.; Huang, K.; Papalexakis, E.; Faloutsos, C. Tensor decomposition for signal processing and machine learning. IEEE Trans. Signal Process. 2017, 65, 3551–3582. [Google Scholar] [CrossRef]
Da Costa, M.N.; Favier, G.; Romano, J.M.T. Tensor modelling of MIMO communication systems with performance analysis and Kronecker receivers. Signal Process. 2018, 145, 304–316. [Google Scholar] [CrossRef] [Green Version]
Ribeiro, L.N.; de Almeida, A.L.; Mota, J.C.M. Separable linearly constrained minimum variance beamformers. Signal Process. 2019, 158, 15–25. [Google Scholar] [CrossRef]
De Lathauwer, L. Signal Processing Based on Multilinear Algebra. Ph.D. Thesis, Katholieke Universiteit Leuven, Leuven, Belgium, 1997. [Google Scholar]
Kolda, T.G.; Bader, B.W. Tensor decompositions and applications. SIAM Rev. 2009, 51, 455–500. [Google Scholar] [CrossRef]
Benesty, J.; Paleologu, C.; Ciochină, S. On the identification of bilinear forms with the Wiener filter. IEEE Signal Process. Lett. 2017, 24, 653–657. [Google Scholar] [CrossRef]
Paleologu, C.; Benesty, J.; Ciochină, S. Adaptive filtering for the identification of bilinear forms. Digit. Signal Process. 2018, 75, 153–167. [Google Scholar] [CrossRef]
Elisei-Iliescu, C.; Stanciu, C.; Paleologu, C.; Benesty, J.; Anghel, C.; Ciochină, S. Efficient recursive least-squares algorithms for the identification of bilinear forms. Digit. Signal Process. 2018, 83, 280–296. [Google Scholar] [CrossRef]
Dogariu, L.-M.; Ciochină, S.; Paleologu, C.; Benesty, J. A connection between the Kalman filter and an optimized LMS algorithm for bilinear forms. Algorithms 2018, 11, 211. [Google Scholar] [CrossRef] [Green Version]
Ribeiro, L.N.; de Almeida, A.L.F.; Mota, J.C.M. Identification of separable systems using trilinear filtering. In Proceedings of the 2015 IEEE 6th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), Cancun, Mexico, 13–16 December 2015; pp. 189–192. [Google Scholar]
Rupp, M.; Schwarz, S. A tensor LMS algorithm. In Proceedings of the 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Brisbane, QLD, Australia, 19–24 April 2015; pp. 3347–3351. [Google Scholar]
Dogariu, L.-M.; Ciochină, S.; Benesty, J.; Paleologu, C. An iterative Wiener filter for the identification of trilinear forms. In Proceedings of the 2019 42nd International Conference on Telecommunications and Signal Processing (TSP), Budapest, Hungary, 1–3 July 2019; pp. 88–93. [Google Scholar]
Dogariu, L.-M.; Ciochină, S.; Benesty, J.; Paleologu, C. System identification based on tensor decompositions: A trilinear approach. Symmetry 2019, 11, 556. [Google Scholar] [CrossRef] [Green Version]
Haykin, S. Adaptive Filter Theory, 4th ed.; Prentice-Hall: Upper Saddle River, NJ, USA, 2002. [Google Scholar]
Parathai, P.; Tengtrairat, N.; Woo, W.L.; Gao, B. Single-channel signal separation using spectral basis correlation with sparse nonnegative tensor factorization. Circuits Syst. Signal Process. 2019, 38, 5786–5816. [Google Scholar] [CrossRef]
Woo, W.L.; Dlay, S.S.; Al-Tmeme, A.; Gao, B. Reverberant signal separation using optimized complex sparse nonnegative tensor deconvolution on spectral covariance matrix. Digit. Signal Process. 2018, 83, 9–23. [Google Scholar] [CrossRef] [Green Version]
Gao, B.; Lu, P.; Woo, W.L.; Tian, G.Y.; Zhu, Y.; Johnston, M. Variational Bayes sub-group adaptive sparse component extraction for diagnostic imaging system. IEEE Trans. Ind. Electron. 2018, 65, 8142–8152. [Google Scholar] [CrossRef]
Kiers, H.A.L. Towards a standardized notation and terminology in multiway analysis. J. Chemom. 2000, 14, 105–122. [Google Scholar] [CrossRef]
Kroonenberg, P. Applied Multiway Data Analysis; Wiley: Hoboken, NJ, USA, 2008. [Google Scholar]
Van Loan, C.F. The ubiquitous Kronecker product. J. Comput. Appl. Math. 2000, 123, 85–100. [Google Scholar] [CrossRef] [Green Version]
Bertsekas, D.P. Nonlinear Programming, 2nd ed.; Athena Scientific: Belmont, MA, USA, 1999. [Google Scholar]
Digital Network Echo Cancellers; ITU-T Recommendations G.168; ITU: Geneva, Switzerland, 2002.
Gay, S.L.; Benesty, J. (Eds.) Acoustic Signal Processing for Telecommunication; Kluwer Academic Publisher: Boston, MA, USA, 2000. [Google Scholar]
Morgan, D.R.; Benesty, J.; Sondhi, M.M. On the evaluation of estimated impulse responses. IEEE Signal Process. Lett. 1998, 5, 174–176. [Google Scholar] [CrossRef]
Ciochină, S.; Paleologu, C.; Benesty, J.; Enescu, A.A. On the influence of the forgetting factor of the RLS adaptive filter in system identification. In Proceedings of the 2009 International Symposium on Signals, Circuits and Systems, Iasi, Romania, 9–10 July 2009; pp. 205–208. [Google Scholar]
Paleologu, C.; Benesty, J.; Ciochină, S. Linear system identification based on a Kronecker product decomposition. IEEE/ACM Trans. Audio Speech Lang. Process. 2018, 26, 1793–1808. [Google Scholar] [CrossRef]
Elisei-Iliescu, C.; Paleologu, C.; Benesty, J.; Ciochină, S. A recursive least-squares algorithm based on the nearest Kronecker product decomposition. In Proceedings of the ICASSP 2019—2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Brighton, UK, 12–17 May 2019; pp. 4843–4847. [Google Scholar]
Elisei-Iliescu, C.; Paleologu, C.; Benesty, J.; Stanciu, C.; Anghel, C.; Ciochină, S. Recursive least-squares algorithms for the identification of low-rank systems. IEEE/ACM Trans. Audio Speech Lang. Process. 2019, 27, 903–918. [Google Scholar] [CrossRef]
Benesty, J.; Rey, H.; Rey Vega, L.; Tressens, S. A non-parametric VSS NLMS algorithm. IEEE Signal Process. Lett. 2006, 13, 581–584. [Google Scholar] [CrossRef]
Paleologu, C.; Benesty, J.; Ciochină, S. A robust variable forgetting factor recursive least-squares algorithm for system identification. IEEE Signal Process. Lett. 2008, 15, 597–600. [Google Scholar] [CrossRef]
Jouppi, N.P.; Young, C.; Patil, N.; Patterson, D.; Agrawal, G.; Bajwa, R.; Bates, S.; Bhatia, S.; Boden, N.; Borchers, A.; et al. In-datacenter performance analysis of a tensor processing unit. In Proceedings of the 44th Annual International Symposium on Computer Architecture, Toronto, ON, Canada, 24–28 June 2017; pp. 1–12. [Google Scholar]

Figure 1. Computational complexity of the proposed RLS-TF algorithm, as compared to the conventional RLS and NLMS algorithms, as a function of

L_{1}

; the other dimensions are set to

L_{2} = 8

and

L_{3} = 4

: (a) number of multiplications per iteration and (b) number of additions per iteration.

Figure 1. Computational complexity of the proposed RLS-TF algorithm, as compared to the conventional RLS and NLMS algorithms, as a function of

L_{1}

; the other dimensions are set to

L_{2} = 8

and

L_{3} = 4

: (a) number of multiplications per iteration and (b) number of additions per iteration.

Figure 2. The components of the third-order system used in simulations: (a)

h_{1}

is the first impulse response (of length

L_{1} = 64

) from the G168 Recommendation [26]; (b)

h_{2}

is a randomly generated impulse response (of length

L_{2} = 8

), with Gaussian distribution; (c) the impulse response

h_{3}

(of length

L_{3} = 4

), with the coefficients computed as

h_{3, l_{3}} = {0.5}^{l_{3} - 1}

, with

l_{3} = 1, 2, \dots, L_{3}

; and (d) the global impulse response (of length

L = L_{1} L_{2} L_{3} = 2048

) results based on (17),

h = h_{3} \otimes h_{2} \otimes h_{1}

Figure 2. The components of the third-order system used in simulations: (a)

h_{1}

is the first impulse response (of length

L_{1} = 64

) from the G168 Recommendation [26]; (b)

h_{2}

is a randomly generated impulse response (of length

L_{2} = 8

), with Gaussian distribution; (c) the impulse response

h_{3}

(of length

L_{3} = 4

), with the coefficients computed as

h_{3, l_{3}} = {0.5}^{l_{3} - 1}

, with

l_{3} = 1, 2, \dots, L_{3}

; and (d) the global impulse response (of length

L = L_{1} L_{2} L_{3} = 2048

) results based on (17),

h = h_{3} \otimes h_{2} \otimes h_{1}

Figure 3. Normalized projection misalignment (NPM) evaluated based on (55)–(57), in dB, for the identification of the individual impulse responses from Figure 2a–c, using the RLS-TF algorithm with different values of the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

(varying the value of K): (a)

NPM [h_{1}, {\hat{h}}_{1} (n)]

, (b)

NPM [h_{2}, {\hat{h}}_{2} (n)]

, and (c)

NPM [h_{3}, {\hat{h}}_{3} (n)]

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

(varying the value of K): (a)

NPM [h_{1}, {\hat{h}}_{1} (n)]

, (b)

NPM [h_{2}, {\hat{h}}_{2} (n)]

, and (c)

NPM [h_{3}, {\hat{h}}_{3} (n)]

Figure 4. Normalized misalignment (NM) evaluated based on (58), in dB, for the identification of the global impulse response

h

(of length

L = 2048

) from Figure 2d, using the RLS-TF algorithm with different values of the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

(varying the value of K).

Figure 4. Normalized misalignment (NM) evaluated based on (58), in dB, for the identification of the global impulse response

h

(of length

L = 2048

) from Figure 2d, using the RLS-TF algorithm with different values of the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

(varying the value of K).

Figure 5. Normalized projection misalignment (NPM) evaluated based on (55)–(57), in dB, for the identification of the individual impulse responses from Figure 2a–c, using the NLMS-TF algorithm [17] (with different normalized step-sizes

α_{1} = α_{2} = α_{3} = α

) and the RLS-TF algorithm (with the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

, where

K = 20

): (a)

NPM [h_{1}, {\hat{h}}_{1} (n)]

, (b)

NPM [h_{2}, {\hat{h}}_{2} (n)]

, and (c)

NPM [h_{3}, {\hat{h}}_{3} (n)]

α_{1} = α_{2} = α_{3} = α

) and the RLS-TF algorithm (with the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

, where

K = 20

): (a)

NPM [h_{1}, {\hat{h}}_{1} (n)]

, (b)

NPM [h_{2}, {\hat{h}}_{2} (n)]

, and (c)

NPM [h_{3}, {\hat{h}}_{3} (n)]

Figure 6. Normalized misalignment (NM) evaluated based on (58), in dB, for the identification of the global impulse response

h

(of length

L = 2048

) from Figure 2d, using the NLMS-TF algorithm [17] (with different normalized step-sizes

α_{1} = α_{2} = α_{3} = α

) and the RLS-TF algorithm (with the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

, where

K = 20

Figure 6. Normalized misalignment (NM) evaluated based on (58), in dB, for the identification of the global impulse response

h

(of length

L = 2048

) from Figure 2d, using the NLMS-TF algorithm [17] (with different normalized step-sizes

α_{1} = α_{2} = α_{3} = α

) and the RLS-TF algorithm (with the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

, where

K = 20

Figure 7. Normalized misalignment (NM) evaluated based on (58), in dB, for the identification of the global impulse response

h

(of length

L = 2048

) from Figure 2d, using the conventional RLS algorithm (with different values of the forgetting factor

λ = 1 - 1 / (K L)

, varying the value of K) and the RLS-TF algorithm (with the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

, where

K = 20

Figure 7. Normalized misalignment (NM) evaluated based on (58), in dB, for the identification of the global impulse response

h

(of length

L = 2048

) from Figure 2d, using the conventional RLS algorithm (with different values of the forgetting factor

λ = 1 - 1 / (K L)

, varying the value of K) and the RLS-TF algorithm (with the forgetting factors

λ_{i} = 1 - 1 / (K L_{i}), i = 1, 2, 3

, where

K = 20

Table 1. RLS algorithm for trilinear forms (RLS-TF).

Initialization:

Set

{\hat{h}}_{1} (0)

{\hat{h}}_{2} (0)

, and

{\hat{h}}_{3} (0)

based on (52)–(54)

R_{32}^{- 1} (0) = \frac{1}{δ_{1}} I_{L_{1}}, R_{31}^{- 1} (0) = \frac{1}{δ_{2}} I_{L_{2}}, R_{21}^{- 1} (0) = \frac{1}{δ_{3}} I_{L_{3}}

(Regularization parameters : δ_{1} > 0, δ_{2} > 0, δ_{3} > 0)

λ_{1} = 1 - \frac{1}{K L_{1}}, λ_{2} = 1 - \frac{1}{K L_{2}}, λ_{3} = 1 - \frac{1}{K L_{3}}

(Tuning constant : K \geq 1)

For n = 1, 2, \dots, number of iterations :

Compute

{\hat{H}}_{32} (n - 1)

{\hat{H}}_{31} (n - 1)

, and

{\hat{H}}_{21} (n - 1)

based on (27)–(29)

Compute

x_{32} (n)

x_{31} (n)

, and

x_{21} (n)

based on (30)–(32)

Evaluate the error signal

e (n)

based on (33)

Compute

k_{32} (n)

k_{31} (n)

, and

k_{21} (n)

based on (49)–(51)

Update

{\hat{h}}_{1} (n)

{\hat{h}}_{2} (n)

, and

{\hat{h}}_{3} (n)

based on (43)–(45)

Update

R_{32}^{- 1} (n)

R_{31}^{- 1} (n)

, and

R_{21}^{- 1} (n)

based on (46)–(48)

Evaluate

\hat{h} (n)

based on (22)

Table 2. Computational complexity of the RLS-TF algorithm, as compared to the conventional RLS and NLMS algorithms.

Algorithms	×	+	÷
RLS	$2 L^{2} + 2 L$	$2 L^{2} + L + 1$	1
RLS-TF	$4 L + 3 (L_{1}^{2} + L_{2}^{2} + L_{3}^{2}) + 3 (L_{1} + L_{2} + L_{3}) + \min (L_{1}, L_{2}, L_{3})$	$3 L + 2 (L_{1}^{2} + L_{2}^{2} + L_{3}^{2}) + L_{1} + L_{2} + L_{3} + \min (L_{1}, L_{2}, L_{3})$	3
NLMS	$2 L + 2$	$2 L + 3$	1

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Elisei-Iliescu, C.; Dogariu, L.-M.; Paleologu, C.; Benesty, J.; Enescu, A.-A.; Ciochină, S. A Recursive Least-Squares Algorithm for the Identification of Trilinear Forms. Algorithms 2020, 13, 135. https://doi.org/10.3390/a13060135

AMA Style

Elisei-Iliescu C, Dogariu L-M, Paleologu C, Benesty J, Enescu A-A, Ciochină S. A Recursive Least-Squares Algorithm for the Identification of Trilinear Forms. Algorithms. 2020; 13(6):135. https://doi.org/10.3390/a13060135

Chicago/Turabian Style

Elisei-Iliescu, Camelia, Laura-Maria Dogariu, Constantin Paleologu, Jacob Benesty, Andrei-Alexandru Enescu, and Silviu Ciochină. 2020. "A Recursive Least-Squares Algorithm for the Identification of Trilinear Forms" Algorithms 13, no. 6: 135. https://doi.org/10.3390/a13060135

APA Style

Elisei-Iliescu, C., Dogariu, L. -M., Paleologu, C., Benesty, J., Enescu, A. -A., & Ciochină, S. (2020). A Recursive Least-Squares Algorithm for the Identification of Trilinear Forms. Algorithms, 13(6), 135. https://doi.org/10.3390/a13060135

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

A Recursive Least-Squares Algorithm for the Identification of Trilinear Forms

Abstract

1. Introduction

2. Third-Order Tensors

3. RLS Algorithm for Trilinear Forms

4. Simulation Results

5. Conclusions and Future Works

Author Contributions

Funding

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI