A Modified Recursive Regularization Factor Calculation for Sparse RLS Algorithm with l₁-Norm

Abstract

In this paper, we propose a new calculation method for the regularization factor in sparse recursive least squares (SRLS) with $l_1$-norm penalty. The proposed regularization factor requires no prior knowledge of the actual system impulse response, and it also reduces computational complexity by about half. In the simulation, we use Mean Square Deviation (MSD) to evaluate the performance of SRLS, using the proposed regularization factor. The simulation results demonstrate that SRLS using the proposed regularization factor calculation shows a difference of less than 2 dB in MSD from SRLS, using the conventional regularization factor with a true system impulse response. Therefore, it is confirmed that the performance of the proposed method is very similar to that of the existing method, even with half the computational complexity.

1. Introduction

The sparse channel system refers to a model characterized by only a few numbers of non-zero taps when modeling the impulse response of the channel system, using the tapped-delay-line model. We can experience such channels in TV channels [1], wide-band radio communication [2], underwater sound channels [3], etc. [4]. There have been many recent works on the use of adaptive estimation for sparse channel estimation; subsequently, various LMS-type algorithms [5,6] such as sparse LMF algorithms [7,8] and sparse LMS/F algorithms [9,10] have been proposed for sparse channel estimation. The RLS algorithm is better suited to fast convergence rates than the LMS-type algorithm. Hence, the RLS-based sparse adaptive filtering algorithm is considered one of the most promising fast algorithms in many system estimation applications, such as channel estimation. There are some algorithms based on sparse RLS [11,12,13,14] as well as the Total Least Squares (TLS)-based algorithm [15,16]. Many sparse RLS algorithms are algorithmically comparable to the plain RLS algorithm; however, the updated equations in most sparse RLS (SRLS) algorithms are not intrinsically recursive, as in the plain RLS algorithm. Eksioglu and Tanc [13] proposed a fully recursive SRLS algorithm that is comparable to the plain RLS algorithm. The sparse RLS algorithms handle the sparsity with $l_1$-norm. Therefore, it is essential to properly select the regularization factor for $l_1$-norm. Many researchers have developed as well as proposed a proper selection for the regularization factor. The authors in [13] also proposed a regularization factor calculation method for the SRLS algorithm. Similar recursive regularization factor selection methods were used in [15,16,17,18,19]. However, the algorithm from [15,16,17,18] is not practical because the regularization factor in [15,16,17] assumed that the true system impulse response is known in advance and set the true system response in part to an arbitrary constant in [18]. The regularization factor selection method in [19] needs no true system impulse response. However, the regularization factor is recursively updated. Therefore, errors in updating the regularization factor are likely to propagate and accumulate. In [20], Lim proposed the regularization factor by estimating the sparsity of the estimated system. Although the regularization factor required no true system information, it added error to the system. Then Lim in [20] utilized it only to TLS-based system modeling. In [21], $l_1$-IWF (iterative Wiener filter) was proposed, which was a kind of a steepest descent method. $l_1$-IWF had a regularization factor without requiring a priori knowledge of the true system response. However, it did not show whether the regularization factor converges to the optimal regularization factor.

In this paper, the major contribution is that we propose a new regularization factor for the SRLS algorithm in [13] and show that it converges to the scaled optimal regularization factor. The proposed regularization factor does not require a priori knowledge of the true system response. The minor contribution is that the proposed regularization factor requires less computation complexity than the regularization in [13].

The remainder of the paper is organized as follows. In Section 2, we reformulate the sparse RLS. In Section 3, a new regularization factor is proposed. In this section, we show that the regularization factor in [13] requires a priori knowledge of the true system impulse response. We also show the calculation complexity of the proposed regularization factor. The simulation condition and results are described and illustrated in Section 4, and a discussion of the results follows in Section 5. We conclude in Section 6.

2. Problem Formulation

In this section, we reformulate the SRLS in [13]. Consider a sparse weight vector ${\mathbf{w}}_o\in {\mathrm{R}}^N$, which represents the channel by a delay line with N taps. By sparsity, the number of significant factors in ${\mathbf{w}}_o$, S, is much lower than its total dimensionality (that is, S ≪ N). The goal is to derive the sparse vector ${\mathbf{w}}_o$ based on the input signal vector $\mathbf{x}\left(n\right)\in {\mathrm{R}}^N$ and the received signal $y\left(n\right)$, which is assumed to be generated by a linear system, such as (1).

$$y\left(n\right)={\mathbf{w}}_o^T\mathbf{x}\left(n\right)+\eta \left(n\right),$$

(1)

where $\eta \left(n\right)$ is the additive noise. Consider the following standard RLS optimization problem, subject to a sparsity constraint.

$$\mathrm{minimize}\left(\widehat{\mathbf{w}}\left(n\right)\right)=\frac{1}{2}{\displaystyle \sum _{m=0}^n}{\lambda }^{n-m}e{\left(m\right)}^2,s.t.{∥\widehat{\mathbf{w}}\left(n\right)∥}_1\le c,$$

(2)

where $e\left(m\right)=y\left(m\right)-\widehat{\mathbf{w}}{\left(n\right)}^T\mathbf{x}\left(m\right)$, $\widehat{\mathbf{w}}\left(n\right)={\left[{\widehat{w}}_1\left(n\right),\cdots ,{\widehat{w}}_N\left(n\right)\right]}^T$, $\mathbf{x}\left(m\right)=\left[x\right(m),x(m-1),$ ${\cdots ,x(m-N+1)]}^T$, $\lambda$ is a forgetting factor, and ${∥\widehat{\mathbf{w}}\left(n\right)∥}_1={\sum }_{k=0}^{N-1}\left|{\widehat{w}}_k\left(n\right)\right|$. As the reformulation in [22], (2) becomes the following:

$$J\left(\widehat{\mathbf{w}}\left(n\right),\gamma \left(n\right)\right)=\frac{1}{2}\zeta \left(\widehat{\mathbf{w}}\left(n\right)\right)+\gamma \left(n\right)\left({∥\widehat{\mathbf{w}}\left(n\right)∥}_1-c\right),$$

(3)

where $\gamma \left(n\right)$ is a real-valued Lagrangian multiplier. We can find the optimal vector by minimizing the regularized cost function (3). The regularized cost function is convex and nondifferentiable; therefore, subgradient analysis is used instead of the normal gradient. When representing a subgradient of f at $\widehat{\mathbf{w}}$ with ${\nabla }^{s}f\left(\widehat{\mathbf{w}}\right)$, the subgradient of $J(\widehat{\mathbf{w}}\left(n\right),\gamma \left(n\right))$ with respect to $\widehat{\mathbf{w}}\left(n\right)$ is as follows.

$${\nabla }^{s}J\left(\widehat{\mathbf{w}}\left(n\right),\gamma \left(n\right)\right)=\frac{1}{2}\nabla \zeta \left(\widehat{\mathbf{w}}\left(n\right)\right)+\gamma \left(n\right){\nabla }^s\left({∥\widehat{\mathbf{w}}\left(n\right)∥}_1\right),$$

(4)

where ${\nabla }^s\left({∥\widehat{\mathbf{w}}\left(n\right)∥}_1\right)=\mathrm{sgn}\left(\widehat{\mathbf{w}}\left(n\right)\right)$ [13]. Here, $\mathrm{sgn}(·)$ means the element-wise sign function. Hence, for the optimal $\widehat{\mathbf{w}}\left(n\right)$ for $J(\widehat{\mathbf{w}}\left(n\right),\gamma \left(n\right))$, we set (4) to 0 at the optimal point. After evaluating the gradient $\nabla \zeta \left(\widehat{\mathbf{w}}\left(n\right)\right)$, (4) can be written as normal equations in matrix form as follows [13].

$$\mathsf{\Phi }\left(n\right)\widehat{\mathbf{w}}\left(n\right)=\mathbf{r}\left(n\right)-\gamma \left(n\right){\nabla }^s\left({∥\widehat{\mathbf{w}}\left(n\right)∥}_1\right),$$

(5)

where $\mathsf{\Phi }\left(n\right)={\displaystyle \sum _{m=0}^n}{\lambda }^{n-m}\mathbf{x}\left(m\right)\mathbf{x}{\left(m\right)}^T=\lambda \mathsf{\Phi }(n-1)+\mathbf{x}\left(n\right)\mathbf{x}{\left(n\right)}^T$ and $\mathbf{r}\left(n\right)={\displaystyle \sum _{m=0}^n}{\lambda }^{n-m}y\left(m\right)$ $\mathbf{x}\left(m\right)=\lambda \mathbf{r}(n-1)+y\left(n\right)\mathbf{x}\left(n\right)$. With $\mathbf{P}\left(n\right)={\mathsf{\Phi }}^{-1}\left(n\right)={\left(\lambda \mathsf{\Phi }(n-1)+\mathbf{x}\left(n\right)\mathbf{x}{\left(n\right)}^T\right)}^{-1}$ and using the matrix inversion lemma, we obtain the following:

$$\mathbf{P}\left(n\right)={\lambda }^{-1}\left(\mathbf{P}(n-1)-\mathbf{k}\left(n\right)\mathbf{x}{\left(n\right)}^T\mathbf{P}(n-1)\right),$$

(6)

where $\mathbf{k}\left(n\right)=\mathbf{P}(n-1)\mathbf{x}\left(n\right)/\left(\lambda +\mathbf{x}{\left(n\right)}^T\mathbf{P}(n-1)\mathbf{x}\left(n\right)\right)$ is the gain vector. Assuming that $\gamma (n-1)$ and ${\nabla }^s{∥\widehat{\mathbf{w}}\left(n-1\right)∥}_1$ do not change considerably over a single time step, the updated equation of $\widehat{\mathbf{w}}\left(n\right)$ can be approximately written as follows [13]:

$$\begin{array}{c}\widehat{\mathbf{w}}\left(n\right)\cong \widehat{\mathbf{w}}(n-1)+\mathbf{k}\left(n\right)\xi \left(n\right)-(1-\lambda )\gamma (n-1)\mathbf{P}\left(n\right){\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1\hfill \\ =\widehat{\mathbf{w}}(n-1)+\mathbf{k}\left(n\right)\xi \left(n\right)+(1-\lambda )\widehat{\gamma }(n-1)\mathbf{P}\left(n\right){\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1,\hfill \end{array}$$

(7)

where $\widehat{\gamma }\left(n\right)=-\gamma \left(n\right)$ and $\xi \left(n\right)=y\left(n\right)-\widehat{\mathbf{w}}{\left(n-1\right)}^T\mathbf{x}\left(n\right)$.

Algorithm 1 summarizes the conventional RLS algorithm and Algorithm 2 summarizes the $l_1$-norm sparse RLS algorithm.

Algorithm 1 conventional RLS algorithm

1: Initialize: $\lambda$, $\mathbf{x}\left(0\right)$, $y\left(0\right)$, ${\mathbf{w}}_{RLS}\left(0\right)=\mathbf{0}$, $\eta$, $\mathbf{P}\left(0\right)={\eta }^{-1}\mathbf{I}$ 2: for$n=1,2,\dots$do 3:     $\mathbf{k}\left(n\right)=\mathbf{P}(n-1)\mathbf{x}\left(n\right)/\left(\lambda +\mathbf{x}{\left(n\right)}^T\mathbf{P}(n-1)\mathbf{x}\left(n\right)\right)$ 4:     $e\left(n\right)=y\left(n\right)-{\mathbf{w}}_{RLS}{\left(n-1\right)}^T\mathbf{x}\left(n\right)$ 5:     $\mathbf{P}\left(n\right)={\lambda }^{-1}\left(\mathbf{P}(n-1)-\mathbf{k}\left(n\right)\mathbf{x}{\left(n\right)}^T\mathbf{P}(n-1)\right)$ 6:     ${\mathbf{w}}_{RLS}\left(n\right)={\mathbf{w}}_{RLS}(n-1)+\mathbf{k}\left(n\right)e\left(n\right)$ 7: end for

Algorithm 2$l_1$-norm sparse RLS algorithm

1: Initialize: $\lambda$, $\mathbf{x}\left(0\right)$, $y\left(0\right)$, $\widehat{\mathbf{w}}\left(0\right)=\mathbf{0}$, $\eta$, $\mathbf{P}\left(0\right)={\eta }^{-1}\mathbf{I}$ 2: for$n=1,2,\dots$do 3:     $\mathbf{k}\left(n\right)=\mathbf{P}(n-1)\mathbf{x}\left(n\right)/\left(\lambda +\mathbf{x}{\left(n\right)}^T\mathbf{P}(n-1)\mathbf{x}\left(n\right)\right)$ 4:     $\xi \left(n\right)=y\left(n\right)-\widehat{\mathbf{w}}{\left(n-1\right)}^T\mathbf{x}\left(n\right)$ 5:     $\mathbf{P}\left(n\right)={\lambda }^{-1}\left(\mathbf{P}(n-1)-\mathbf{k}\left(n\right)\mathbf{x}{\left(n\right)}^T\mathbf{P}(n-1)\right)$ 6:     $\mathrm{calculate}\widehat{\gamma }(n-1)$ 7:     $\widehat{\mathbf{w}}\left(n\right)=\widehat{\mathbf{w}}(n-1)+\mathbf{k}\left(n\right)\xi \left(n\right)+(1-\lambda )\widehat{\gamma }(n-1)\mathbf{P}\left(n\right){\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1$ 8: end for

3. Proposed Recursive Regularization Factor for Sparse RLS Algorithm

In this section, we derive the regularization factor, $\widehat{\gamma }(n-1)$, such that ${∥\widehat{\mathbf{w}}(n-1)∥}_1=c$, which means that the $l_1$-norm of $\widehat{\mathbf{w}}\left(n\right)$ is preserved for all the time steps in n. This property yields the time invariance in $\widehat{\mathbf{w}}\left(n\right)$ [23].

$$\begin{array}{c}\frac{d{∥\widehat{\mathbf{w}}\left(n\right)∥}_1}{dn}={\left(\frac{\partial {∥\widehat{\mathbf{w}}\left(n\right)∥}_1}{\partial \widehat{\mathbf{w}}\left(n\right)}\right)}^T\frac{d\widehat{\mathbf{w}}\left(n\right)}{dn}={\left({\nabla }^s{∥\widehat{\mathbf{w}}\left(n\right)∥}_1\right)}^T\left(\frac{\left(\widehat{\mathbf{w}}\left(n\right)-\widehat{\mathbf{w}}(n-T_s)\right)}{T_s}\right)\hfill \\ \simeq {\left({\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1\right)}^T\left(\frac{\left(\widehat{\mathbf{w}}\left(n\right)-\widehat{\mathbf{w}}(n-T_s)\right)}{T_s}\right)=0.\hfill \end{array}$$

(8)

By referring to (7) and assuming a normalized sample time $T_s$ to 1 [24], (8) is as follows.

$${\left({\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1\right)}^T\left(\mathbf{k}\left(n\right)\xi \left(n\right)+(1-\lambda )\widehat{\gamma }(n-1)\mathbf{P}\left(n\right){\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1\right)=0.$$

(9)

Hence, $\widehat{\gamma }(n-1)$ in (9) becomes the following:

$$\widehat{\gamma }(n-1)=-\frac{{\left({\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1\right)}^T\mathbf{k}\left(n\right)\xi \left(n\right)}{\left(1-\lambda \right){\left({\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1\right)}^T\mathbf{P}\left(n\right){\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1}.$$

(10)

Eksioglu derived the upper bound of the regularization parameter $\widehat{\gamma }\left(n\right)$ in [13] as follows:

$$\overline{\gamma }\left(n\right)=2\frac{{\left({\nabla }^s{∥\widehat{\mathbf{w}}\left(n\right)∥}_1\right)}^T\mathbf{P}\left(n\right)\widetilde{ϵ}\left(n\right)}{{∥\mathbf{P}\left(n\right){\nabla }^s{∥\widehat{\mathbf{w}}\left(n\right)∥}_1∥}_2^2},$$

(11)

where $\widetilde{ϵ}\left(n\right)=\tilde{\mathbf{w}}\left(n\right)-{\mathbf{w}}_o$ and $\tilde{\mathbf{w}}\left(n\right)$ is the solution to the conventional non-regularized normal equation as $\tilde{\mathbf{w}}\left(n\right)=\mathbf{P}\left(n\right)\mathbf{r}\left(n\right)$ [25]. In (11), the true system impulse response, ${\mathbf{w}}_o$, must be known in advance for $\widetilde{ϵ}\left(n\right)=\tilde{\mathbf{w}}\left(n\right)-{\mathbf{w}}_o$.

If we assume that $x\left(n\right)$ is a white noise process, its autocorrelation matrix, $\mathsf{\Phi }\left(n\right)$, asymptotically becomes ${\sigma }_x^2\mathbf{I}$ and then $\mathbf{P}\left(n\right)=\mathsf{\Phi }{\left(n\right)}^{-1}\cong {\sigma }_x^{-2}(1-\lambda )\mathbf{I}$ [25]. Therefore, (11) asymptotically becomes the following:

$$\overline{\gamma }\left(n\right)\cong \frac{2{\sigma }_x^2}{\left(1-\lambda \right)}\frac{{\left({\nabla }^s{∥\widehat{\mathbf{w}}\left(n\right)∥}_1\right)}^T\widetilde{ϵ}\left(n\right)}{{∥{\nabla }^s{∥\widehat{\mathbf{w}}\left(n\right)∥}_1∥}_2^2}.$$

(12)

Eksioglu also showed that, if the regularization parameter $\widehat{\gamma }\left(n\right)$ is $0\le \widehat{\gamma }\left(n\right)\le \overline{\gamma }\left(n\right)$, the following inequality is satisfied.

$${∥\widehat{\mathbf{w}}\left(n\right)-{\mathbf{w}}_o∥}_2^2\le {∥\tilde{\mathbf{w}}\left(n\right)-{\mathbf{w}}_o∥}_2^2.$$

(13)

This means that the estimation error with regularization is lower than the estimation error without regularization.

In order to review the asymptotic meaning of the proposed regularization in (10), we also assume $x\left(n\right)$ is a white noise process whose autocorrelation matrix, $\mathsf{\Phi }\left(n\right)$, asymptotically becomes ${\sigma }_x^2\mathbf{I}$ and then $\mathbf{P}\left(n\right)=\mathsf{\Phi }{\left(n\right)}^{-1}\cong {\sigma }_x^{-2}(1-\lambda )\mathbf{I}$ [25]. We use the simplified gain vector k(n) in [25].

$$\mathbf{k}\left(n\right)=\mathbf{P}\left(n\right)\mathbf{x}\left(n\right),$$

(14)

and approximate the error, $\xi \left(n\right)$, in high SNR,

$$\begin{array}{c}\xi \left(n\right)=y\left(n\right)-{\mathbf{x}}^T\left(n\right)\widehat{\mathbf{w}}(n-1)\cong {\mathbf{x}}^T\left(n\right){\mathbf{w}}_o-{\mathbf{x}}^T\left(n\right)\widehat{\mathbf{w}}(n-1)\hfill \\ =-{\mathbf{x}}^T\left(n\right)\left(\widehat{\mathbf{w}}(n-1)-{\mathbf{w}}_o\right)\cong -{\mathbf{x}}^T\left(n\right)\left(\tilde{\mathbf{w}}(n-1)-{\mathbf{w}}_o\right).\hfill \end{array}$$

(15)

When we substitute (10) with (14) and (15), we can rewrite the proposed regularization in (10) as follows:

$$\widehat{\gamma }(n-1)=\frac{{\left({\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1\right)}^T\mathbf{P}\left(n\right)\mathbf{x}\left(n\right){\mathbf{x}}^T\left(n\right)\left(\widehat{\mathbf{w}}(n-1)-{\mathbf{w}}_o\right)}{\left(1-\lambda \right){\left({\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1\right)}^T\mathbf{P}\left(n\right){\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1}.$$

(16)

In (16), $\mathbf{P}\left(n\right)\mathbf{x}\left(n\right){\mathbf{x}}^T\left(n\right)$ can be asymptotically approximated as the expected value of (17) in [26].

$$E\left(\mathbf{P}\left(n\right)\mathbf{x}\left(n\right){\mathbf{x}}^T\left(n\right)\right)\cong \frac{\left(1-\lambda \right)}{\left(1-{\lambda }^{k+1}\right)}\mathbf{I}\cong \left(1-\lambda \right)\mathbf{I}.$$

(17)

Substituting (16) with (17) and $\mathbf{P}\left(n\right)\cong {\sigma }_x^{-2}(1-\lambda )\mathbf{I}$, the variable regularization factor in (16) becomes the following:

$$\widehat{\gamma }(n-1)\cong \frac{{\sigma }_x^2}{1-\lambda }\frac{{\left({∥\widehat{\mathbf{w}}(n-1)∥}_1\right)}^T\widetilde{ϵ}(n-1)}{{∥{\nabla }^s{∥\widehat{\mathbf{w}}(n-1)∥}_1∥}_2^2}=\frac{1}{2}\overline{\gamma }(n-1),$$

(18)

$$\widehat{\gamma }(n-1)=\frac{1}{2}\overline{\gamma }(n-1).$$

(19)

In (19), the proposed regularization parameter is approximately converse to the scaled optimal regularization parameter in [13]. Therefore, the newly derived regularization parameter $\widehat{\gamma }\left(n\right)$ in (19) satisfies $0\le \widehat{\gamma }\left(n\right)\le \overline{\gamma }\left(n\right)$. Therefore, the proposed regularization parameter can be used as a regularization parameter on behalf of the optimal regularization parameter. In addition, the proposed regularization parameter needs no true system parameter ${\mathbf{w}}_o$.

In terms of computational complexity, the difference between SRLS in [13] and the proposed algorithm lies only in the computational complexity of the regularization factor. Therefore, we can compare the computational complexity between the proposed regularization factor in (10) and the regularization factor in (11) from [13]. Actually, the regularization factor used in [13] is (20), which is an approximation of (11).

$$\begin{array}{c}\widehat{\gamma }(n-1)=2\left(\frac{\frac{tr\left(\mathbf{P}(n-1)\right)}{N}\left(f\left(\widehat{\mathbf{w}}(n-1)\right)-{\mathbf{w}}_o\right)}{{∥\mathbf{P}(n-1){\nabla }^{s}f\left(\widehat{\mathbf{w}}(n-1)\right)∥}_2^2}\right.\left.+\frac{{\nabla }^{s}f\left(\widehat{\mathbf{w}}(n-1)\right)\mathbf{P}(n-1)\left(\widehat{\mathbf{w}}\left(n\right)-{\widehat{\mathbf{w}}}_{RLS}(n-1)\right)}{{∥\mathbf{P}(n-1){\nabla }^{s}f\left(\widehat{\mathbf{w}}(n-1)\right)∥}_2^2}\right)\hfill \\ =2\left(\frac{\frac{tr\left(\mathbf{P}(n-1)\right)}{N}\left({∥\widehat{\mathbf{w}}(n-1)∥}_1-{∥{\mathbf{w}}_o∥}_1\right)}{{∥\mathbf{P}(n-1)\mathrm{sgn}\left(\widehat{\mathbf{w}}(n-1)\right)∥}_2^2}\right.\left.+\frac{\mathrm{sgn}\left(\widehat{\mathbf{w}}(n-1)\right)\mathbf{P}(n-1)\left(\widehat{\mathbf{w}}\left(n\right)-{\widehat{\mathbf{w}}}_{RLS}(n-1)\right)}{{∥\mathbf{P}(n-1)\mathrm{sgn}\left(\widehat{\mathbf{w}}(n-1)\right)∥}_2^2}\right),\hfill \end{array}$$

(20)

where ${\mathbf{w}}_o$ is the true system response, $f\left(\widehat{\mathbf{w}}\left(n\right)\right)={∥\widehat{\mathbf{w}}\left(n\right)∥}_1$, ${\nabla }^s{∥\widehat{\mathbf{w}}∥}_1=\mathrm{sgn}\left(\widehat{\mathbf{w}}\right)$ and ${\widehat{\mathbf{w}}}_{RLS}\left(n\right)$ is a solution from the conventional RLS algorithm [13]. And applying ${\nabla }^s{∥\widehat{\mathbf{w}}∥}_1=sgn\left(\widehat{\mathbf{w}}\right)$ to the proposed regularization factor (10), it becomes as follows:

$$\begin{array}{c}\widehat{\gamma }(n-1)=-\frac{{\nabla }^{s}f{\left(\widehat{\mathbf{w}}(n-1)\right)}^T\mathbf{k}\left(n\right)\zeta \left(n\right)}{(1-\lambda ){\nabla }^{s}f{\left(\widehat{\mathbf{w}}(n-1)\right)}^T\mathbf{P}(n-1){\nabla }^{s}f\left(\widehat{\mathbf{w}}(n-1)\right)}\hfill \\ =-\frac{\mathrm{sgn}{\left(\widehat{\mathbf{w}}(n-1)\right)}^T\mathbf{k}\left(n\right)\zeta \left(n\right)}{(1-\lambda )\mathrm{sgn}{\left(\widehat{\mathbf{w}}(n-1)\right)}^T\mathbf{P}(n-1)\mathrm{sgn}\left(\widehat{\mathbf{w}}(n-1)\right)}.\hfill \end{array}$$

(21)

The regularization factor is calculated in line 6 in Algorithm 2. When calculating the complexity of the regularization factor, it should be taken into account that the regularization factor uses the elements calculated before line 6. Taking this into account and calculating the complexity, (20) requires $2N+3$ multiplications, whereas (21) requires $N+1$ multiplications.

4. Simulation Results

For the simulation in this section, we set the same experimental conditions as in [13] (code is available at https://web.itu.edu.tr/eksioglue/pubs.htm accessed on 16 March 2021). We assume two system parameters, ${\mathbf{w}}_o$, with N = 64 taps and 256 taps, respectively. Out of the N coefficients, the only S coefficients are not zero. We generate the values of coefficients from an N(0, 1/S) distribution and randomly place the non-zero coefficients.

In the simulation, we show the system estimation results, using the proposed regularization factor that does not require the true system response information. We compare the $l_1$-RLS using the true system response [13], the $l_1$-RLS in [19], $l_1$-IWF in [21] and the $l_1$-RLS, using the proposed regularization factor selection method. In addition, we compare the estimation result from the conventional RLS. We simulate these algorithms in the sparse impulse response for S = 2, 4, 8, and 16 for the performance evaluation. Figure 1 shows the Mean Square Deviation (MSD) comparison results in the case of an order $N=64$ at SNR = 20 dB, where $\mathrm{MSD}\left(\widehat{\mathbf{w}}\left(n\right)\right)=E\left\{{∥{\mathbf{w}}_o-\widehat{\mathbf{w}}\left(n\right)∥}^2\right\}$. Figure 2 also shows the MSD comparison results in the case of an order $N=256$ at SNR = 20 dB. In addition, the MSD comparison results are summarized in

Table 1. MSD comparison in various SNR.

Channel Length	No. of Non-Zero Coefficients	SNR	20 dB	10 dB	0 dB
L = 64	S = 4	proposed method	−30.0	−20.2	−10.2
		$l_1$-RLS with true impulse response	−30.8	−20.6	−11.0
		$l_1$-RLS from [19]	−28.7	−18.7	−8.7
		$l_1$-IWF from [21]	−29.5	−19.7	−9.0
		conventional RLS	−27.7	−17.5	−7.8
	S = 16	proposed method	−28.4	−18.6	−9.1
		$l_1$-RLS with true impulse response	−28.5	−18.5	−9.3
		$l_1$-RLS from [19]	−28.5	−18.4	−8.6
		$l_1$-IWF from [21]	−29.0	−18.8	−9.2
		conventional RLS	−27.6	−17.7	−7.8
L = 256	S = 4	proposed method	−24.4	−14.3	−4.2
		$l_1$-RLS with true impulse response	−25.3	−15.5	−5.5
		$l_1$-RLS from [19]	−22.7	−12.7	−2.6
		$l_1$-IWF from [21]	−25.2	−15.2	−3.6
		conventional RLS	−21.2	−11.2	−1.4
	S = 16	proposed method	−24.2	−14.3	−4.4
		$l_1$-RLS with true impulse response	−24.6	−14.6	−5.0
		$l_1$-RLS from [19]	−22.7	−12.5	−2.6
		$l_1$-IWF from [21]	−24.8	−14.8	−3.4
		conventional RLS	−21.1	−11.1	−1.4

at SNR = 20 dB, 10 dB and 0 dB, respectively.

5. Discussion

The MSD comparison results in Figure 1 show that the estimation performance of $l_1$-RLS, using the regularization factor of the proposed method, is almost the same as the $l_1$-RLS, using the regularization factor with the true system impulse response. $l_1$-IWF also performs almost the same as the $l_1$-RLS with the true system information as the results in [21]. The results from the proposed algorithm are even better than those of the $l_1$-RLS in [19]. Predictably, the conventional RLS has the worst MSD in all cases. Figure 1 confirms that, despite having no prior knowledge of the true system response, the proposed regularization factor is comparable to the conventional regularization factor with the true system impulse response.

Figure 2 shows the MSD results for a longer-order system: $N=256$. We can see that Figure 2 has very similar results to those in Figure 1. Therefore, it can be confirmed that the proposed regularization factor calculation method works well, regardless of the system dimension.

Table 1 shows the performance of sparse RLS, using the proposed regularization factor calculation method, compared with other algorithms in various SNR situations. The performance of $l_1$-IWF is also similar to that of the proposed algorithm, but the performance of the proposed algorithm is better at low SNR. The results of Table 1 shows that, although SRLS using the proposed regularization factor does not actually utilize the impulse response of the true system, the MSD difference between SRLS with the regularization factor using information on the actual impulse response of the target system, and SRLS with the proposed regularization factor, is less than 2 dB. Therefore, it can be said that the performance of the two SRLS is very similar. In addition, as mentioned in Section 3, it is also remarkable that the computational complexity of the regularization factor can be reduced by about half of the conventional regularization factor.

6. Conclusions

In this paper, we proposed a new calculation method for a regularization factor in $l_1$-RLS, requiring no prior knowledge of the true system response. We also showed that the proposed regularization factor converges to the scaled optimal regularization factor. Therefore, we have shown that the proposed regularization factor can be used on behalf of the optimal regularization factor. The simulation results confirmed that the proposed regularization factor behaves almost the same as the conventional regularization factor with the true system response.