Low-Complexity Recursive Least-Squares Adaptive Algorithm Based on Tensorial Forms

Abstract

Modern solutions for system identification problems employ multilinear forms, which are based on multiple-order tensor decomposition (of rank one). Recently, such a solution was introduced based on the recursive least-squares (RLS) algorithm. Despite their potential for adaptive systems, the classical RLS methods require a prohibitive amount of arithmetic resources and are sometimes prone to numerical stability issues. This paper proposes a new algorithm for multiple-input/single-output (MISO) system identification based on the combination between the exponentially weighted RLS algorithm and the dichotomous descent iterations in order to implement a low-complexity stable solution with performance similar to the classical RLS methods.

1. Introduction

Tensor decomposition is an important topic for the development of modern technologies such as big data, machine learning, communication systems with multiple inputs and/or multiple outputs, respectively source separation/beamforming [1,2,3]. When implementing system identification tasks, the problems can be reformulated to split the arithmetic workload into multiple smaller system identification problems. Considerable recent research was dedicated to identify multilinear forms [4,5,6] with two, three, or even more components using the Wiener filter [7,8,9] and adaptive algorithms based on the least-mean-square (LMS) method [10,11,12,13] or on the recursive least-squares (RLS) systems [10,14,15]. Despite the attractive performance of the latter, its classical implementations are considered too costly for most practical applications and the preferred workhorse remains the LMS or the normalized-LMS (NLMS) method, which is not sensitive to the scaling of its input signal.

The traditional least-squares based adaptive systems require complexities proportional to the second or third power of the filter’s length in order to compute the corresponding matrix inverse and the associated coefficients. The filter update employs a full spectrum of mathematical operations, including divisions, and it is often plagued by numerical instability. Despite its prohibitive nature in terms of required chip area, the decorrelation properties of the RLS (translated in practice into superior convergence speeds and overall performance) pushed researchers to develop new corresponding versions with attractive compromise between complexity and performance. However, for some algorithms, such as Fast-RLS (FRLS) [16], the negative aspects still tend to overtake the value brought by the reduction in terms of complexity and deny the practical applications a reliable least-squares algorithm.

For the system identification implementations, which can be decomposed using multiple-input/single-output (MISO) setups, several tensor based models were recently proposed [11]. One such decomposition uses the RLS method based on Woodburry’s identity [15] to split the unknown system determination into multiple smaller adaptive systems. Despite the fact that the overall complexity is reduced, the algorithm previously introduced as tensor-based RLS (RLS-T) is still dependent on the square of each filter’s length, and it is also prone to inherit the problems of classical least-squares solutions [14].

Recently, amongst the newly proposed RLS methods, the combination with the dichotomous coordinate descent (DCD) iterations has been established as a possible stable alternative with lower complexity traits, proportional to the adaptive filter’s length. This paper proposes a tensorial decomposition for multilinear forms based on the RLS-DCD method, which was extensively studied and widely accepted as an efficient solution for least-squares adaptive systems [17,18]. The DCD portion of the algorithm exploits the nature of the corresponding correlation matrices and estimates with minimal computational effort, in an iterative manner, an incremental vector for each of the adaptive filter to be updated. The operations are performed using only additions and bit-shifts, and any overall requirement for divisions is eliminated [19,20]. For the identification of an unknown MISO system based on its tensorial form, multiple RLS-DCD smaller adaptive filters will be employed, which inherit the performance of the classical RLS versions and require lower arithmetic efforts.

The paper is organized as follows: Section 2 introduces the framework of the multilinear forms. In Section 3, the RLS-T is described with its corresponding notations and the newly RLS-DCD-T adaptive systems are introduced as an alternative for the former. The new adaptive model is validated in Section 4 through several simulation results and conclusions are presented in Section 5.

2. System Model

Considering N individual channels described by the following vectors:

$$\begin{array}{c}\hfill {\mathbf{h}}_i={\left[h_{i,1}{h}_{i,2}\cdots h_{i,L_i}\right]}^T,\end{array}$$

(1)

each one of length $L_i,i=1,2,\dots ,N$, where the superscript ${}^T$ denotes the transpose operator, they compose a multiple-input/single-output (MISO) system. For efficiency reasons, in the rest of the paper, the index i will repeatedly have the same significance as the one mentioned for Equation (1) (i.e., to identify elements associated with the i-th channel). Its corresponding real valued output is represented at the discrete-time index n by the expression:

$$\begin{array}{c}\hfill y\left(n\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\dots l_N}\left(n\right)h_{1,l_1}{h}_{2,l_2}\cdots h_{N,l_N},\end{array}$$

(2)

where the input signals from (2) can be aggregated using the tensorial form $\mathcal{X}\left(n\right)\in {\mathbb{R}}^{L_1\times L_2\times \cdots \times L_N}$, with the constituent elements ${\left(\mathcal{X}\right)}_{l_1{l}_2\dots l_N}\left(n\right)=x_{l_1{l}_2\dots l_N}\left(n\right)$. Such linearly separable systems are subject of research for applications like nonlinear acoustic echo cancellation [21,22], beamforming [3,23], or source separation techniques [1,2].

Consequently, the input–output equation can be written as

$$\begin{array}{c}\hfill y\left(n\right)=\mathcal{X}\left(n\right){\times }_1{\mathbf{h}}_1^T{\times }_2{\mathbf{h}}_2^T{\times }_3\cdots {\times }_N{\mathbf{h}}_N^T,\end{array}$$

(3)

where ${\times }_i$ is the mode-i product explained in [1], and $y\left(n\right)$ is a linear function of each of the vectors ${\mathbf{h}}_i$, when the other $N-1$ vectors are fixed (i.e., $y\left(n\right)$ connotes a multilinear form). Due to this hypothesis, the structure presented in (2) and (3) can be interpreted as an extension of the bilinear form [24].

Moreover, we denote the tensor $\mathcal{H}\in {\mathbb{R}}^{L_1\times L_2\times \cdots \times L_N}$, with the elements ${\left(\mathcal{H}\right)}_{l_1,l_2,\dots ,l_N}=h_{1,l_1}{h}_{2,l_2}\cdots h_{N,l_N}$, as the result of:

$$\begin{array}{c}\hfill \mathcal{H}={\mathbf{h}}_1\circ {\mathbf{h}}_2\circ \cdots \circ {\mathbf{h}}_N,\end{array}$$

(4)

where ∘ symbolizes the outer product, defined as ${\mathbf{h}}_1\circ {\mathbf{h}}_2={\mathbf{h}}_1{\mathbf{h}}_2^T$, with the corresponding elements ${\left({\mathbf{h}}_1\circ {\mathbf{h}}_2\right)}_{i,j}=h_{1,i}{h}_{2,j}$. At this stage, it is noteworthy that the Kronecker product between two vectors (denoted using ⊗) can be expressed as ${\mathbf{h}}_2\otimes {\mathbf{h}}_1=\mathrm{vec}\left({\mathbf{h}}_1\circ {\mathbf{h}}_2\right)$, where $\mathrm{vec}(·)$ is the vectorization operation.

Therefore, we write

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

(5)

and we reformulate the output signal as

$$\begin{array}{c}\hfill y\left(n\right)={\mathrm{vec}}^T\left(\mathcal{H}\right)\mathrm{vec}\left[\mathcal{X}\left(n\right)\right],\end{array}$$

(6)

where

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

(7)

with ${\mathbf{X}}_{::\dots :l_i}\left(n\right)\in {\mathbb{R}}^{L_1\times L_2\times \cdots \times L_{N-1}}$ representing the frontal slices of $\mathcal{X}\left(n\right)$ and $\mathbf{x}\left(n\right)$ denoting the input vector of length L = $L_1{L}_2\cdots L_N$. In addition,

$$\begin{array}{c}\hfill \mathbf{g}=\mathrm{vec}\left(\mathcal{H}\right)={\mathbf{h}}_N\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{h}}_1\end{array}$$

(8)

symbolizes the global impulse response of length $L_1{L}_2\cdots L_N$. The decomposition from (8) is not given by a unique set of estimators corresponding to the individual members. Exemplifying, let us consider N real-valued constants $c_i$, having the property that ${\prod }_{i=1}^N{c}_i=1$. Then,

$$\begin{array}{c}\hfill c_N{\mathbf{h}}_N\otimes c_{N-1}{\mathbf{h}}_{N-1}\otimes \cdots \otimes c_1{\mathbf{h}}_1={\mathbf{h}}_N\otimes {\mathbf{h}}_{N-1}\otimes \cdots \otimes {\mathbf{h}}_1.\end{array}$$

(9)

Nevertheless, the global impulse response $\mathbf{g}$ can be identified with no scaling ambiguity.

Furthermore, considering (6)–(8), $y\left(n\right)$ can be rewritten as

$$\begin{array}{c}\hfill y\left(n\right)={\mathbf{g}}^T\mathbf{x}\left(n\right).\end{array}$$

(10)

The target is to identify the unknown global system $\mathbf{g}$ given by its individual components ${\mathbf{h}}_i$ by using the reference signal:

$$\begin{array}{c}\hfill d\left(n\right)={\mathbf{g}}^T\mathbf{x}\left(n\right)+w\left(n\right),\end{array}$$

(11)

where the measurement noise, $w\left(n\right)$, is uncorrelated with the incoming signal. The reference signal from (11) is described by the variance:

$$\begin{array}{c}\hfill {\sigma }_d^2={\mathbf{g}}^{T}E\left[\mathbf{x}\left(n\right){\mathbf{x}}^T\left(n\right)\right]\mathbf{g}+{\sigma }_w^2={\mathbf{g}}^T\mathbf{R}\mathbf{g}+{\sigma }_w^2.\end{array}$$

(12)

$E[·]$ denotes mathematical expectation, $\mathbf{R}=E\left[\mathbf{x}\left(n\right){\mathbf{x}}^T\left(n\right)\right]$, and ${\sigma }_w^2=E\left[w^2\left(n\right)\right]$.

The error signal has the expression:

$$\begin{array}{c}\hfill e\left(n\right)=d\left(n\right)-{\widehat{\mathbf{g}}}^T\mathbf{x}\left(n\right),\end{array}$$

(13)

where $\widehat{\mathbf{g}}$ describes an estimator of the global impulse response.

The Wiener approach for the estimate $\widehat{\mathbf{g}}$ [14] is not suitable in nonstationary conditions. Moreover, the requirement of real-time processing makes it only a starting point in providing more adequate solutions based on adaptive algorithms. The majority of adaptive systems rely on the LMS-based methods, which offer a satisfactory balance between performance and arithmetic workload. However, state-of-the-art literature provides several versions of the RLS adaptive filter as appealing alternatives with better performance indicators (i.e., tracking speed, convergence rate). Despite of the advantages which the RLS-based versions have over their LMS-based counterparts, most of the RLS algorithms cannot provide stable functionality combined with proper arithmetic complexity. In recent years, the integration between the dichotomous coordinate descend (DCD) iterative methods [20] with the RLS adaptive filters [14] successfully contested the prohibitive tag associated with the RLS family implementations. The RLS-DCD delivered promising results with acceptable processing requirements and performance similar to the classical RLS versions.

Now, having (13), we can define the mean-squared error (MSE) criterion:

$$\begin{array}{c}\hfill J\left(\widehat{\mathbf{g}}\right)=E\left[e^2\left(n\right)\right]={\sigma }_d^2-2{\widehat{\mathbf{g}}}^T\mathbf{p}+{\widehat{\mathbf{g}}}^T\mathbf{R}\widehat{\mathbf{g}},\end{array}$$

(14)

with $\mathbf{p}=E\left[\mathbf{x}\left(n\right)d\left(n\right)\right]$ representing the cross-correlation vector between $\mathbf{x}\left(n\right)$ and $d\left(n\right)$. Our focus is the minimization of $J\left(\widehat{\mathbf{g}}\right)$, and it represents the optimization criterion. In order to do this, we may use the well-known solution [14] given by the Wiener filter:

$$\begin{array}{c}\hfill {\widehat{\mathbf{g}}}_{\mathrm{W}}={\mathbf{R}}^{-1}\mathbf{p}.\end{array}$$

(15)

With the iterative Wiener filter [8], we are able to speculate the decomposition of the global impulse response, considering the topic of the multilinear forms. The optimization criterion consists on a block coordinate descent method applied on the individual components [25]. Comparing to the standard Wiener filter given by Equation (15), the iterative one boosts the performance. It is suitable, especially when a small amount of data are used in order to estimate the statistics of the signal (i.e., $\mathbf{R}$ and $\mathbf{p}$).

3. Tensor-Based RLS Algorithms

The faster convergence rate is one of the main advantages of the RLS methods, with respect to the performances obtained by the LMS-based family of algorithms. However, this convergence aspect comes together with a very high computational complexity. It was previously demonstrated that tensor-based algorithms could produce better results than the classical RLS approach, by splitting the long filter associated with the unknown system identification problem into multiple smaller filters, i.e., into multiple smaller system identification processes. The main advantages of this category of adaptive algorithms are two-fold: a significant decrease in number of mathematical operations and improved convergence rates.

3.1. Tensor-Based Recursive Least Squares Algorithm (RLS-T)

Starting with the desired (11) and error (13) signals, we can apply the least-squares error criterion [14] for the tensor-based decomposition. The RLS algorithm can be applied by following the minimization of the cost function:

$$\begin{array}{c}\hfill J\left[\widehat{\mathbf{g}}\left(n\right)\right]=\sum _{k=1}^n{\lambda }^{n-k}{\left[d\left(k\right)-{\widehat{\mathbf{g}}}^T\left(n\right)\mathbf{x}\left(k\right)\right]}^2,\end{array}$$

(16)

with the forgetting factor of the algorithm, $0<\lambda \le 1$.

First, let us consider the estimated impulse responses of the channels, ${\widehat{\mathbf{h}}}_i\left(n\right)$. Therefore, we may define the corresponding a-priori error signals:

$$\begin{array}{cc}\hfill e_{{\widehat{\mathbf{h}}}_i}\left(n\right)& =d\left(n\right)-y_{{\widehat{\mathbf{h}}}_i}\left(n\right)\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & =d\left(n\right)-{\widehat{\mathbf{h}}}_i^T(n-1){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right),\hfill \end{array}$$

(18)

where $y_{{\widehat{\mathbf{h}}}_i}\left(n\right)$ are the outputs of the individual channels, and we used the form $\mathbf{x}\left(n\right)$ defined in (7) to express

$$\begin{array}{c}\hfill {\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)={\left[{\widehat{\mathbf{h}}}_N(n-1)\otimes {\widehat{\mathbf{h}}}_{N-1}(n-1)\otimes \cdots \otimes {\widehat{\mathbf{h}}}_{i+1}(n-1)\otimes {\mathbf{I}}_{L_i}\otimes {\widehat{\mathbf{h}}}_{i-1}(n-1)\otimes \cdots \otimes {\widehat{\mathbf{h}}}_2(n-1)\otimes {\widehat{\mathbf{h}}}_1(n-1)\right]}^T\mathbf{x}\left(n\right),\end{array}$$

(19)

with ${\mathbf{I}}_{L_i}$ denoting the identity matrices of sizes $L_i\times L_i$. It can be verified that $e_{{\widehat{\mathbf{h}}}_1}\left(n\right)=e_{{\widehat{\mathbf{h}}}_2}\left(n\right)=\cdots =e_{{\widehat{\mathbf{h}}}_N}\left(n\right)$.

Considering (19), the cost function in (16) can be reworked in N alternative ways by performing the optimization procedure of the individual impulse responses, i.e.,

$$\begin{array}{c}\hfill J_{{\widehat{\mathbf{h}}}_N,{\widehat{\mathbf{h}}}_{N-1},\dots ,{\widehat{\mathbf{h}}}_{i+1},{\widehat{\mathbf{h}}}_{i-1},\dots ,{\widehat{\mathbf{h}}}_2,{\widehat{\mathbf{h}}}_1}\left[{\widehat{\mathbf{h}}}_i\left(n\right)\right]=\sum _{k=1}^n{\lambda }_i^{n-k}{\left[d\left(k\right)-{\widehat{\mathbf{h}}}_i^T\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(k\right)\right]}^2,\end{array}$$

(20)

where the ${\lambda }_i$’s are the individual forgetting factors. These impulse responses can be processed by following a multilinear optimization strategy [25], considering the fact that $N-1$ components are fixed, and the optimization method is applied to the remaining one. Therefore, the focus is the minimization of (20) with respect to ${\widehat{\mathbf{h}}}_i\left(n\right)$, which leads to a set of normal equations:

$$\begin{array}{c}\hfill {\mathbf{R}}_i\left(n\right){\widehat{\mathbf{h}}}_i\left(n\right)={\mathbf{p}}_i\left(n\right),\end{array}$$

(21)

where

$$\begin{array}{cc}\hfill {\mathbf{R}}_i\left(n\right)& =\sum _{k=1}^n{\lambda }_i^{n-k}{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(k\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^T\left(k\right)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & ={\lambda }_i{\mathbf{R}}_i(n-1)+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^T\left(n\right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\mathbf{p}}_i\left(n\right)& =\sum _{k=1}^n{\lambda }_i^{n-k}{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(k\right)d\left(k\right)\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & ={\lambda }_i{\mathbf{p}}_i(n-1)+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)d\left(n\right).\hfill \end{array}$$

(25)

We can achieve N update relations linked to the individual filters by solving Equation (21). The obvious requirement for its corresponding solutions is the computation of the inverse matrices ${\mathbf{R}}_i\left(n\right)$. The straightforward approach is considered too costly for real-time practical applications. A well-known and accepted compromise, which is often used as a reference, is the method provided by the matrix inversion lemma (also known as Woodbury’s matrix identity) [14]. We can express the inverse matrices as

$$\begin{array}{c}\hfill {\mathbf{R}}_i^{-1}\left(n\right)=\frac{1}{{\lambda }_i}\left[{\mathbf{I}}_{L_i}-{\mathbf{k}}_i\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^T\left(n\right)\right]{\mathbf{R}}_i^{-1}(n-1),\end{array}$$

(26)

where

$$\begin{array}{c}\hfill {\mathbf{k}}_i\left(n\right)=\frac{{\mathbf{R}}_i^{-1}(n-1){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)}{{\lambda }_i+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^T\left(n\right){\mathbf{R}}_i^{-1}(n-1){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)}\end{array}$$

(27)

denotes the Kalman gain vectors.

Consequently, the filter updates can be performed using

$$\begin{array}{cc}\hfill {\widehat{\mathbf{h}}}_i\left(n\right)& ={\widehat{\mathbf{h}}}_i(n-1)+{\mathbf{R}}_i^{-1}\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)e_{{\widehat{\mathbf{h}}}_i}\left(n\right)\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & ={\widehat{\mathbf{h}}}_i(n-1)+{\mathbf{k}}_i\left(n\right)e_{{\widehat{\mathbf{h}}}_i}\left(n\right),\hfill \end{array}$$

(29)

where ${\mathbf{k}}_i\left(n\right)={\mathbf{R}}_i^{-1}\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)$ are the Kalman gain vectors. The error signal can be easily computed by using (18).

The arithmetic complexity associated with the computation of each set of filter coefficients is proportional to the square of the corresponding filter length. The usage of a single adaptive filter is similarly prohibitive for hardware implementations.

In the end, the estimate of the resulting filter can be obtained using

$$\begin{array}{c}\hfill \widehat{\mathbf{g}}\left(n\right)={\widehat{\mathbf{h}}}_N\left(n\right)\otimes {\widehat{\mathbf{h}}}_{N-1}\left(n\right)\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\left(n\right).\end{array}$$

(30)

The tensor-based RLS (RLS-T) algorithm is available in Algorithm 1. This approach has been deeply analyzed in the paper [26] by using $N=3$ (i.e., third-order tensors). The current sub-section described a generalization for $N\ge 3$, thus providing even more efficient implementation solutions based on RLS adaptive filters. However, the RLS-T form requires for each of the N filters a computational workload proportional to $L_i^2$. When considering applications where one of the filters still has a considerable amount of weights, the cost of practical implementations remains quite high.

Algorithm 1: RLS-T algorithm
Step	Actions
	$\mathrm{Initialization}:$
	Set ${\widehat{\mathbf{h}}}_i\left(0\right)$
0	${\mathbf{R}}_i^{-1}\left(0\right)={\delta }_i^{-1}{\mathbf{I}}_{L_i},{\delta }_i>0$
	${\lambda }_i=1-\frac{1}{K_i{L}_i},K_i\ge 1$
	$\mathrm{For}n=1,2,\dots ,\mathrm{number}\mathrm{of}\mathrm{iterations}:$
1	Compute ${\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)$, based on (19)
2	$e_{{\widehat{\mathbf{h}}}_i}\left(n\right)=d\left(n\right)-{\widehat{\mathbf{h}}}_i^T(n-1){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)$
3	${\mathbf{k}}_i\left(n\right)=\frac{{\mathbf{R}}_i^{-1}(n-1){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)}{{\lambda }_i+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^T\left(n\right){\mathbf{R}}_i^{-1}(n-1){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)}$
4	${\widehat{\mathbf{h}}}_i\left(n\right)={\widehat{\mathbf{h}}}_i(n-1)+{\mathbf{k}}_i\left(n\right)e_{{\widehat{\mathbf{h}}}_i}\left(n\right)$
5	${\mathbf{R}}_i^{-1}\left(n\right)=\frac{1}{{\lambda }_i}\left[{\mathbf{I}}_{L_i}-{\mathbf{k}}_i\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^T\left(n\right)\right]{\mathbf{R}}_i^{-1}(n-1)$
	$\widehat{\mathbf{g}}\left(n\right)={\widehat{\mathbf{h}}}_N\left(n\right)\otimes {\widehat{\mathbf{h}}}_{N-1}\left(n\right)\otimes \cdots \otimes {\widehat{\mathbf{h}}}_1\left(n\right)$

3.2. Tensor-Based Recursive Least-Squares Dichotomous Coordinate Descent Algorithm (RLS-DCD-T)

We propose to use the combination between the dichotomous coordinate descent (DCD) iterations and the RLS method as an alternative solution for the systems of equations in (21). The RLS-DCD algorithm [20] was employed in the past due to its low complexity arithmetic workloads and improved numerical stability. The advantages obtained by using the generalized tensorial model applied with low-complexity RLS algorithms can lead to high convergence/tracking speeds and acceptable computational requirements, an overall design which is suitable for efficient hardware applications.

As was similarly considered in [20,27], we assume that each of the N filters corresponding to the tensor-based decomposition has generated a result at time index $n-1$ for its associated systems in (21). A residual vector for each solution can be expressed as

$$\begin{array}{c}\hfill {\mathbf{r}}_i(n-1)={\mathbf{p}}_i(n-1)-{\mathbf{R}}_i(n-1){\widehat{\mathbf{h}}}_i(n-1)\end{array}$$

(31)

and several other notations can be employed to track the changes of the entities in (21) with respect to n:

$$\begin{array}{cc}\hfill \Delta {\mathbf{R}}_i\left(n\right)& ={\mathbf{R}}_i\left(n\right)-{\mathbf{R}}_i(n-1),\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \Delta {\mathbf{p}}_i\left(n\right)& ={\mathbf{p}}_i\left(n\right)-{\mathbf{p}}_i(n-1),\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \Delta {\mathbf{h}}_i\left(n\right)& ={\mathbf{h}}_i\left(n\right)-{\widehat{\mathbf{h}}}_i(n-1),\hfill \end{array}$$

(34)

where $\Delta {\mathbf{R}}_i\left(n\right)$ is a $L_i\times L_i$ matrix, and $\Delta {\mathbf{p}}_i\left(n\right)$, respectively $\Delta {\mathbf{h}}_i\left(n\right)$, are $L_i$-valued column vectors.

The classical methods for directly solving the systems in (21) can be shifted towards determining the estimates $\Delta {\widehat{\mathbf{h}}}_i\left(n\right)$ and adding the vectors to the estimates generated by the previous filter update process ${\widehat{\mathbf{h}}}_i(n-1)$, in order to obtain new filter coefficients in ${\widehat{\mathbf{h}}}_i\left(n\right)$. Consequently, for time index n, (21) can be rewritten as:

$$\begin{array}{c}\hfill {\mathbf{R}}_i\left(n\right)\left[{\widehat{\mathbf{h}}}_i(n-1)+\Delta {\mathbf{h}}_i\left(n\right)\right]={\mathbf{p}}_i\left(n\right).\end{array}$$

(35)

From (35), we can start formulating new systems of equations, similar to (21), with the purpose of determining the estimates $\Delta {\widehat{\mathbf{h}}}_i\left(n\right)$. On the left side of (35), we extract the term

$$\begin{array}{c}\hfill {\mathbf{R}}_i\left(n\right)\Delta {\mathbf{h}}_i\left(n\right)={\mathbf{p}}_i\left(n\right)-{\mathbf{R}}_i\left(n\right){\widehat{\mathbf{h}}}_i(n-1)={\mathbf{p}}_{0,i}\left(n\right).\end{array}$$

(36)

Furthermore, (36) can be developed by using (32) and (33), then (31), in order to eliminate ${\mathbf{p}}_i\left(n\right)$, respectively ${\mathbf{R}}_i\left(n\right)$, and we obtain:

$$\begin{array}{c}\hfill {\mathbf{p}}_{0,i}\left(n\right)={\mathbf{r}}_i(n-1)+\Delta {\mathbf{p}}_i\left(n\right)-\Delta {\mathbf{R}}_i\left(n\right){\widehat{\mathbf{h}}}_i(n-1).\end{array}$$

(37)

The reinterpretation of (21), which led to (36), generated N new systems of equations with the form

$$\begin{array}{c}\hfill {\mathbf{R}}_i\left(n\right)\Delta {\mathbf{h}}_i\left(n\right)={\mathbf{p}}_{0,i}\left(n\right),\end{array}$$

(38)

where $\Delta {\widehat{\mathbf{h}}}_i\left(n\right)$ are the new solutions that will be computed. The approach is equivalent to (21) because the new filter coefficients are determined with

$$\begin{array}{c}\hfill {\widehat{\mathbf{h}}}_i\left(n\right)={\widehat{\mathbf{h}}}_i(n-1)+\Delta {\widehat{\mathbf{h}}}_i\left(n\right).\end{array}$$

(39)

Moreover, the solutions obtained from (38) are based on ${\widehat{\mathbf{h}}}_i(n-1)$ estimated at the previous time index, and the computation of the N increment vectors $\Delta {\widehat{\mathbf{h}}}_i\left(n\right)$ is expected to require considerably less computational resources than the direct approaches. This makes the form in (38) suitable for iterative methods, such as the DCD iterations, which can also update the values of the residual vectors ${\mathbf{r}}_i\left(n\right)$ using a minimum computational effort. They are required at the next time index, as can be noticed in (37). The solution for this apparent issue is to use (33) and (34), then ${\mathbf{p}}_{0,i}$, in

$$\begin{array}{c}\hfill {\mathbf{r}}_i\left(n\right)={\mathbf{p}}_i\left(n\right)-{\mathbf{R}}_i\left(n\right){\widehat{\mathbf{h}}}_i\left(n\right),\end{array}$$

(40)

which leads to the conclusion that

$$\begin{array}{c}\hfill {\mathbf{r}}_i\left(n\right)={\mathbf{p}}_{0,i}\left(n\right)-{\mathbf{R}}_i\left(n\right)\Delta {\widehat{\mathbf{h}}}_i\left(n\right).\end{array}$$

(41)

The form in (41) fits well with the solutions that must be determined for (38) because most of the DCD versions presented in corresponding literature [19,20] can compute with a reduced arithmetic workload (and using the same entities) both the target estimates $\Delta {\widehat{\mathbf{h}}}_i\left(n\right)$, respectively the residual vectors ${\mathbf{r}}_i\left(n\right)$.

Finally, in order to replace the classical RLS methods employed in the tensorial structure with the previously described approach, a convenient way of computing ${\mathbf{p}}_{0,i}\left(n\right)$ is also necessary. When the cost functions in (20) are considered, the corresponding exponentially weighted version of the RLS algorithm uses (23), respectively (25), to update the elements in (21) needed for the direct computation of ${\widehat{\mathbf{h}}}_i\left(n\right)$. We will employ (32) in (23), respectively (33) in (25), to obtain

$$\begin{array}{cc}\hfill \Delta {\mathbf{R}}_i\left(n\right)& =({\lambda }_i-1){\mathbf{R}}_i(n-1)+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^T\left(n\right),\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill \Delta {\mathbf{p}}_i\left(n\right)& =({\lambda }_i-1){\mathbf{p}}_i(n-1)+d\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right).\hfill \end{array}$$

(43)

Therefore, by replacing $\Delta {\mathbf{R}}_i\left(n\right)$ and $\Delta {\mathbf{p}}_i\left(n\right)$ in (37), then using (31) and (18), we can write

$$\begin{array}{cc}\hfill {\mathbf{p}}_{0,i}\left(n\right)& ={\mathbf{r}}_i(n-1)+({\lambda }_i-1){\mathbf{p}}_i(n-1)+d\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & -({\lambda }_i-1){\mathbf{R}}_i(n-1){\widehat{\mathbf{h}}}_i(n-1)-{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^T\left(n\right){\widehat{\mathbf{h}}}_i(n-1)\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & ={\lambda }_i{\mathbf{r}}_i(n-1)+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)e_{{\widehat{\mathbf{h}}}_i}\left(n\right).\hfill \end{array}$$

(45)

We can summarize the newly proposed tensor-based RLS-DCD-T algorithm in Algorithm 2 using (23), (18), (45), and (39). For each of the corresponding channels, the RLS-DCD-T is designed with an overall complexity proportional to the length of the associated adaptive filter, in terms of additions and multiplications. No divisions are needed to perform the filter update process or to generate the output information.

Algorithm 2: Exponential weighted RLS-T algorithm for one channel
Step	Actions	Complexity ‘×’ & ‘+’
	$\mathrm{Initialization}:$
	Set ${\widehat{\mathbf{h}}}_i\left(0\right)={\mathbf{0}}_{L_i\times 1},{\mathbf{r}}_i\left(0\right)={\mathbf{0}}_{L_i\times 1}$
0	${\mathbf{R}}_i^{-1}\left(0\right)={\xi }_i{\mathbf{I}}_{L_i},{\xi }_i>0$
	${\lambda }_i=1-\frac{1}{K_i{L}_i},K_i\ge 1$
	$\mathrm{For}n=1,2,\dots ,\mathrm{number}\mathrm{of}\mathrm{iterations}:$
1	Compute ${\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)$, based on (19)	L + ${\displaystyle \sum _{i=1}^{N-1}}{\displaystyle \prod _{j=i}^N}{L}_j$ & $L-1$
2	${\mathbf{R}}_i^{\left(0\right)}\left(n\right)={\lambda }_i{\mathbf{R}}_i^{\left(0\right)}(n-1)+{\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}^{\left(0\right)}\left(n\right)$	$2L_i$ & $L_i$
3	$y_{{\widehat{\mathbf{h}}}_i}\left(n\right)={\widehat{\mathbf{h}}}_i^T(n-1){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)$	$L_i$ & $L_i-1$
4	$e_{{\widehat{\mathbf{h}}}_i}\left(n\right)=d\left(n\right)-y_{{\widehat{\mathbf{h}}}_i}\left(n\right)$	0 & 1
5	${\mathbf{p}}_{0,i}\left(n\right)={\lambda }_i{\mathbf{r}}_i(n-1)+e_{{\widehat{\mathbf{h}}}_i}\left(n\right){\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)$	$2L_i$ & $L_i$
6	${\mathbf{R}}_i\left(n\right)\Delta {\mathbf{h}}_i\left(n\right)={\mathbf{p}}_{0,i}\left(n\right)\stackrel{\mathrm{DCD}}{\to }\Delta {\widehat{\mathbf{h}}}_i\left(n\right),{\mathbf{r}}_i\left(n\right)$	0 & $N_{u,i}$(2$L_i+1$)+$M_b$
7	${\widehat{\mathbf{h}}}_i\left(n\right)={\widehat{\mathbf{h}}}_i(n-1)+\Delta {\widehat{\mathbf{h}}}_i\left(n\right)$	0 & $L_i$

3.3. DCD Method and Arithmetic Complexity

In Algorithm 2, the first step of each adaptive filter corresponding to the N channels is the computation the input vectors ${\mathbf{x}}_{{\widehat{\mathbf{h}}}_i}\left(n\right)$. It requires L + ${\displaystyle \sum _{i=1}^{N-1}}{\displaystyle \prod _{j=i}^N}{L}_j$ multiplications, respectively $L-1$ additions. When considering all the channels, the overall number of additions is $N(L-1)$, and the dominant term for the multiplications is $2NL$. The rest of the steps are expected to require arithmetic efforts proportional to the individual filter lengths $L_i$.

Considering that the input signals have the time-shift property [19,27], the matrix updates in step 2 can be performed by moving the upper-left $L_i\times L_i$ sub-matrix of ${\mathbf{R}}_i\left(n\right)$ to the lower-right $L_i\times L_i$ sub-matrix, then computing the first column and transposing the result in order to also insert the values in the first line of ${\mathbf{R}}_i\left(n\right)$. This procedure ensures a complexity (in terms of multiplications and additions) proportional to the filter lengths $L_i$. Similar magnitudes apply when quantifying arithmetic efforts corresponding to the computation of each filter output, respectively each of the ${\mathbf{p}}_{0,i}\left(n\right)$ vectors, in steps 3 and 5. Step 4 can be considered negligible in terms of arithmetic workload, as it requires a number of N simple subtractions (for all channels).

At step 6, the DCD iterations are used for each channel independently to determine the increment vectors $\Delta {\widehat{\mathbf{h}}}_i\left(n\right)$ and to update the values comprising the corresponding residual vectors ${\mathbf{r}}_i\left(n\right)$. The method exploits the positive-definite nature of ${\mathbf{R}}_i\left(n\right)$ and performs a maximum number of allowed $N_{u,i}$ numerical updates on $\Delta {\widehat{\mathbf{h}}}_i\left(n\right)$, which is initialized at every time index with the zero valued column vector ${\mathbf{0}}_{L_i\times 1}$.

The DCD variant we propose to use (mostly known as the DCD with a leading element) is presented in Algorithm 3. For $N_{u,i}$ times (or less—see the conditions for the while loop in Algorithm 3), it applies a greedy approach by selecting the maximum absolute value of ${\mathbf{p}}_{0,i}\left(n\right)$, and it decides on performing (or not) an update on the corresponding position (or coordinate) in $\Delta {\mathbf{h}}_i\left(n\right)$. A comparison is made between the selected value in ${\mathbf{p}}_{0,i}\left(n\right)$ and the value situated on the same position from the main diagonal of ${\mathbf{R}}_i\left(n\right)$ [19]. The coefficient update and its corresponding decision are also influenced by the so-called step-size (denoted by $\alpha$), which is changed by halving it and is initialized considering the values of $\Delta {\mathbf{h}}_i\left(n\right)$—they are expected to be situated in the interval $\left[-H;H\right]$. By choosing for H a power of 2, all the multiplications performed with the step size $\alpha$ can be replaced by bit-shifts. A short analysis of step 5 in Algorithm 3 reveals that the updates made for the values of $\Delta {\mathbf{h}}_i\left(n\right)$ actually represent bits set to 1 in the corresponding binary representations. A natural choice of $H=1$ determines that the values of the step size are always negative powers of 2 and implicitly influence which bit is to be updated on the chosen position of $\Delta {\mathbf{h}}_i\left(n\right)$ (in case of a successful iteration). The values of $m_i$ are upper limited by the number of bits chosen to represent the magnitudes in the solution vectors, denoted by $M_b$ (the same value applies for all coefficients from all filters). Consequently, the DCD stops when the maximum number of allowed updates $N_{u,i}$ is reached or any more changes to the solution values would be made for $m_i>M_b$. In other words, the latter situation is equivalent to saying that the step-size is too small with respect to the number of bits chosen to represent the values comprising $\Delta {\mathbf{h}}_i\left(n\right)$, and implicitly ${\widehat{\mathbf{h}}}_i\left(n\right)$. Finally, it can be noticed in Algorithm 3 that a successful iteration does not only perform a change on the solution vector—it additionally updates the associated residual vector ${\mathbf{r}}_i\left(n\right)$, which is initialized with zeros at time index n, and it will be used at the next filter iteration.

Algorithm 3: The DCD iterations with a leading element and overall complexity
Step	Action	Complexity ‘+’
	$\mathrm{Initialization}:{\alpha }_i=H_i/2,m_i=1$
	$\Delta {\widehat{\mathbf{h}}}_i={\mathbf{0}}_{L_i\times 1},{\mathbf{r}}_i={\mathbf{0}}_{L_i\times 1}$
	$\mathrm{For}(q=1...N_{u,i}):$
1	$p=arg{max}_{k=1,2,\dots ,L_i}\left\{|r_{k,i}|\right\}$	$N_{u,i}$× ($L_i-1$)
2	$\mathrm{While}((m_i\le M_{b,i})$ && $\left(\right|r_{p,i}|\ge ({\alpha }_i/2)R_{p,p,i})):$	(≤ $M_b$) + (≤ $N_{u,i}$)
	$m_i++,{\alpha }_i={\alpha }_i/2$	0
3	$\mathrm{if}{m}_i>M_b\stackrel{}{\to }\mathrm{the}\mathrm{algorithm}\mathrm{stops}$	0
4	$\Delta {\widehat{h}}_{p,i}=\Delta {\widehat{h}}_{p,i}+\mathrm{sign}\left(r_{p,i}\right){\alpha }_i$	$N_{u,i}$
5	${\mathbf{r}}_i={\mathbf{r}}_i-\mathrm{sign}\left(r_{p,i}\right){\alpha }_i{\mathbf{R}}_i^{\left(p\right)}$	$N_{u,i}$×$L_i$
		Overall(worse case):
		$N_{u,i}$(2$L_i+1$) + $M_b$

The arithmetic processes needed for the hardware implementation of the leading DCD do not require multiplications. In addition, considering the fact that the quantity of additions is variable due to the greedy behavior, we approach this estimate in terms of worse case (or the upper limits). It can be noticed in Algorithm 3 that the highest computational workload is necessary when all $N_{u,i}$ allowed updates are performed and the decision process passes through all possible values for $m_i$ (the while loop). The status of $m_i$ and the value of the step size $\alpha$ can be tracked using shift registries. Moreover, any apparent multiplication operations involving $\alpha$ are actually replaceable with bit-shifts.

In step 1 from Algorithm 3, the determination of p requires $L_i-1$ comparisons, which are equivalent to additions. This process can be necessary for a maximum amount of $N_{u,i}$ times. In a similar manner, the comparison in step 2 can be employed no more than $N_{u,i}+M_b$ total times (equivalent to the maximum numbers of successful and unsuccessful DCD iterations). Finally, steps 4 and 5 are performed when a successful iteration is triggered. Consequently, for maximum $N_{u,i}$ times, one addition is spent on updating a bit for some value in $\Delta {\mathbf{h}}_i\left(n\right)$ (step 4) and $L_i$ additions are used for updating ${\mathbf{r}}_i\left(n\right)$ (step 5).

Rounding up, from a hardware implementation point of view, the leading DCD iterations require for $N_{u,i}$ successful updates $N_{u,i}$(2$L_i+1$) additions (corresponding to steps 1, 4, 5 and the true results of the second condition in the while loop), plus a maximum of $M_b$ additions for the false results of the condition in step 2. The rest of the apparent arithmetic workload can be performed using bit-shift operations.

Consequently, the complexity of the Exponential Weighted RLS-T algorithm presented in Algorithm 2 reflects its split functionality design. Excepting step 1, the overall computational effort is a sum of values proportional to the individual filter lengths, in terms of multiplications, respectively additions. Considering the fact that decompositions can be performed such that $L_i\ll L$, the proposed reduction in complexity for steps 2–7 represents a migration from a setup difficult to implement on hardware platforms to an attractive solution for multiple adaptive systems configurations.

Using the conventional RLS family of algorithms [14] (i.e., the direct estimation of the global impulse response) could be very costly for large values of L, since the computational complexity order would be $\mathcal{O}\left(L^2\right)$. On the other hand, the computational complexity of the RLS-T algorithms is proportional to ${\sum }_{i=1}^N\mathcal{O}\left(L_i^2\right)$, which could be much more advantageous when $L_i\ll L$. In addition, since the RLS-T operates with shorter filters, improved performance is expected, as compared to the conventional RLS algorithm. The RLS-DCD-T brings an extra layer of efficiency by performing the same tasks with workloads of order ${\sum }_{i=1}^N\mathcal{O}\left(L_i\right)$. These observations are also supported by the simulation results provided in the next section.

4. Simulations and Practical Considerations

Several sets of experiments were performed in the context of the system identification problem with the MISO configuration. The input signals are white Gaussian noises with sub-unitary values. For some scenarios, before being fed to the input of the system, the signals are also filtered through an auto-regressive AR(1) system with the pole 0.85, in order to create highly-correlated input sequences [11]. For $d\left(n\right)$, the signal-to-noise (SNR) ratio is experimentally set to 15 dB.

In the context of adaptive filtering, the convergence behavior and performance analysis of any adaptive algorithm consider the Gaussian input as a benchmark. When the input is white Gaussian noise, the conditional number of the correlation matrix of the input signal is equal to one, since this matrix is a diagonal one and all its eigenvalues are equal. On the other hand, in case of more correlated inputs, the conditional number of this matrix increases, resulting in the challenges in terms of convergence rate. However, the recursive least-squares (RLS) algorithm is less influenced by the character of the input signal (in terms of the correlation degree), as compared to the family of least-mean-square (LMS) algorithms [14]. Consequently, it is important to show that the proposed RLS-DCD-T algorithm also inherits this feature. This is the reason for including both cases in our experiments, i.e., white Gaussian noise and AR(1) process as inputs.

The criterion employed to analyze the performance is the normalized misalignment, corresponding to the global system identification process. It is defined for each time index, using the global filter estimate $\widehat{\mathbf{g}}\left(n\right)$ and the true impulse response $\mathbf{g}$, as $20{\log }_{10}\left[{∥\mathbf{g}-\widehat{\mathbf{g}}\left(n\right)∥}_2/{∥\mathbf{g}∥}_2\right]$, where ${∥·∥}_2$ denotes the ${\ell }_2$ norm and $\mathbf{g}$ is computed based on (8).

The experiments use global impulse responses which can be decomposed using (8) into $N=3$, respectively $N=4$, smaller systems of lengths $L_i,i=1,2,\dots ,N$ (depending on the case). The first and the fourth small impulse responses are taken as the first $L_1$, respectively $L_4$, samples of the the B4 and B1 impulse responses from the G.168 ITU-T Recommendation [28]. The second and third systems are generated using $(-0.5^{j_1})/2$, respectively $(0.8^{j_2})/2$, where $j_1=0...L_2-1$ and $j_2=0...L_3-1$. In order to evaluate the tracking capabilities of the presented algorithms, the global system is suddenly modified after the initial convergence is achieved by changing the sign of the coefficients corresponding to one of the unknown channels (the longest one, when the case).

The experimental results will mainly compare the performances of the RLS-DCD-T, respectively the RLS-T, adaptive algorithms. For both tensor-based methods, the adaptive filter lengths are chosen as $L_i=2^{q_i}$ and the forgetting factors are set as ${\lambda }_i=1-1/\left(2^{k_i}{L}_i\right)=1-1/\left(2^{k_i+q_i}\right)$. From a hardware implementation point of view, the multiplications with lambda can be performed more efficiently with only bit-shifts and subtractions [18,29]. Furthermore, for the RLS-DCD-T, its DCD specific parameters are chosen considering several aspects. Firstly, the tap values for the estimated systems are sub-unitary ($H=1$), and they are represented using $M_b=16$ bits—one for the sign and the rest, up to 16, hold the magnitude information. Secondly, the maximum number of allowed updates was successfully used in the past with values lower than 10 for adaptive filters with hundreds and thousands of coefficients [18,27]. Taking into account that the selected scenarios have filters with the corresponding lengths $L_i$ of maximum tens of coefficients, we chose not to display simulations with $N_{u,i}>1$ because they showed no improvement in performance, even for the longer adaptive filters in the tensorial decompositions. For each of the tensor-based update processes (i.e., for each simulated channel), this setting leads to arithmetic complexities proportional to $L_i$.

Although the tensor-based adaptive filters are initialized using the method recommended in [11], the initial convergence is not shown, as it is not considered relevant. For comparison reasons, the classical RLS method based on Woodburry’s identity for the direct identification of the unknown global impulse response is also simulated as a reference, as it is an already validated/established solution in scientific literature. The corresponding forgetting factors were chosen in a similar manner as the equivalent parameters of the tensor based methods (i.e., $\lambda =1-1/\left(2^{k}L\right)$).

The first experiments were performed with a configuration of $N=3$ filters for the identification of a global system with the length $L=1024$ ($L_1=64$, $L_2=L_3=4$), where the input signal is Gaussian noise. It can be noticed in Figure 1 that the RLS-DCD-T with $N_{u,i}=1$ is enough to generate practically identical performances in comparison to the RLS-T algorithm. Moreover, the RLS adaptive method with different forgetting factors was simulated. It is obvious from the illustration that decreasing the forgetting factor for the RLS in order to increase the tracking speed deteriorates the system identification performance at a convergence state. It is also noteworthy that the individual adaptive filters are connected, i.e., the solution of one filter at time index n depends on the estimates of the other $N-1$ filters at time index $n-1$, due to Equation (19). Consequently, the estimation errors also cumulate in this manner, which could explain the steady-state variation of the tensor-based algorithms (RLS-T and RLS-DCD-T), as compared to the conventional RLS algorithm.

The experiment was repeated for an AR(1) process in Figure 2. Despite the fact that the RLS-DCD-T slightly loses performance in comparison to the RLS-T, similar conclusions can be applied when performing the comparison with the RLS. The tensor based methods outperform the RLS adaptive filter, which cannot simultaneously match them in terms of tracking speed, respectively performance when convergence is achieved.

For the third and fourth scenarios, the MISO configuration was changed to a decomposition based on $N=4$ filters. The global system has $L=2048$ coefficients and for the tensor-based systems $L_1=16$, $L_2=8$, respectively $L_3=L_4=4$. In Figure 3, the input signal is Gaussian noise. The performances corresponding to the tensor based algorithms are again very similar. In comparison to the approach which directly identifies the global impulse response (i.e., the RLS algorithm), the RLS-DCD-T and RLS-T have much better tracking capabilities, when targeting for all the methods similar performance after convergence is reached. When repeating the experiment for an AR(1) sequence (Figure 4), it can be considered clear that, for both situations, any attempt of improving the tracking capabilities of the RLS adaptive filter by reducing its memory (i.e., by using a smaller value of $\lambda$) leads to a compromise with the quality the unknown global system identification at the convergence state.

The final scenario analyzes the robustness of the RLS-DCD-T adaptive algorithm in low SNR conditions. The employed MISO configuration is for a global impulse response with $L=4096$, where the tensor based decomposition is performed for $N=4$ channels with $L_1=L_2=L_3=L_4=8$. For an interval of 5000 iterations, the SNR changes from the value of 15 dB to −15 dB. The performances corresponding to the RLS-DCD-T, respectively the RLS-T, are illustrated in Figure 5. The RLS-DCD-T demonstrates the same robustness as the RLS-T when affected by low SNR conditions, with a similar number of iterations required to update its filter coefficients and to achieve again the previous misalignment for the global system estimate $\widehat{\mathbf{g}}\left(n\right)$.

5. Conclusions

This paper proposed a low-complexity adaptive algorithm for the identification of unknown systems based on tensorial decompositions. We used as a starting point the RLS-T adaptive method, and we developed a new algorithm based on the combination between the exponentially weighted RLS and the DCD iterations. The newly RLS-DCD-T method benefits from the low computational requirements of the DCD and provides performance comparable with other established versions of tensorial based RLS methods. The reduction in complexity for the adaptive filter update process is twofold. Firstly, the RLS-DCD-T benefits from the tensorial approach and solves multiple smaller system identification problems instead of a single problem, which is prohibitive for practical applications. Secondly, the usage of the DCD iterations allows for the coefficient updates to be performed using only bit-shifts and additions.

In future works, we aim to further develop the combination between the least-squares methods and the DCD working within the tensorial decomposition framework. The DCD will be paired with sliding windows versions of the RLS. Moreover, we will develop various methods for increasing the robustness of the RLS-DCD-T in very low SNR conditions, with minimal additions to the arithmetic complexity.