# Subband Adaptive Filtering withl1-Norm Constraint for Sparse System Identification

Authors: Young-Seok Choi
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper presents a new approach of the normalized subband adaptive filter (NSAF) which directly exploits the sparsity condition of an underlying system for sparse system identification. The proposed NSAF integrates a weightedl1-norm constraint into the cost function of the NSAF algorithm. To get the optimum solution of the weighted l1-norm regularized cost function, a subgradient calculus is employed, resulting in a stochastic gradient based update recursion of the weighted l1-norm regularized NSAF. The choice of distinct weighted l1-norm regularization leads to two versions of the l1-norm regularized NSAF. Numerical results clearly indicate the superior convergence of the l1-norm regularized NSAFs over the classical NSAF especially when identifying a sparse system.


---

## Body

## 1. Introduction

Over the past few decades, the relative simplicity and good performance of the normalized least mean square (NLMS) algorithm have made it a popular tool for adaptive filtering applications. However, its convergence performance is significantly deteriorated in case of correlated input signals [1, 2]. As a popular solution, adaptive filtering in the subband has been recently developed, which is referred to as subband adaptive filter (SAF) [3–7]. Its distinct feature is based on the property that the LMS-type adaptive filters converge faster for white input signals than colored ones [1, 2]. Thus, carrying out a prewhitening on colored input signals, it results in the accelerated convergence compared to the LMS-type adaptive filters. Recently, the use of multiple-constraint optimization criteria into formulation of a cost function has resulted in the normalized SAF (NSAF) with its computational complexity close to that of the NLMS algorithm [6, 7].In the context of a system identification, the unknown system to be identified is sparse in common scenarios, such as echo paths [8] and digital TV transmission channels [9]. Namely, the unknown system consists of many near-zero coefficients and a small number of large ones. However, the adaptive filtering algorithms suffer from poor convergence performance in case of identifying the sparse system [8]. Indeed, the capability of the NSAF is faded in a sparse system identification scenario. To deal with this issue, a variety of proportionate adaptive algorithms have been presented for NSAF, which utilize proportionate step sizes to distinct filter taps [10–12]. However, these algorithms have not exploited the sparsity condition of an underlying system.Recently, motivated by compressive sensing framework [13, 14] and the least absolute shrinkage and selection operator (LASSO) [15], a number of adaptive filtering algorithms which make use of the sparsity condition of an underlying system have been developed [16–20]. The core idea behind this approach is to incorporate the sparsity condition of underlying system by imposing an a sparsity-inducing constraint term. Adding the sparsity constraint using l0 or l1-norm constraint to the cost function makes the least relevant weights of the filter shrink to zeros. However, to the best of the author’s knowledge, adaptive filtering in subband which exploits the sparsity condition has not been studied yet.With regard, this paper presents a novel approach of the sparsity-regularized NSAFs, which incorporates the sparsity condition of the system directly into the cost function via a sparsity-inducing constraint term. This is carried out by regularizing a weightedl1-norm of the filter weights estimate to the cost function. Considering the two choices of the weighted l1-norm regularization, two stochastic gradient-based l1-norm regularized NSAF algorithms are derived. First, the l1-norm NSAF (l1-NSAF) is obtained by utilizing the identity matrix as a weighting matrix. Second, the reweighted l1-norm NSAF (l1-RNSAF) which uses the current estimate of the system as a weighted l1-norm is developed. Through numerical simulations, the resultant sparsity-regularized NSAFs have proven their superiority over the classical NSAFs, especially when the sparsity of the underlying system becomes severe.The remainder of the paper is organized as follows. Section2 introduces the classical NSAF, followed by the derivation of the proposed sparsity-regularized NSAFs in Section 3. Section 4 illustrates the computer simulation results and Section 5 concludes this study.


## 2. Conventional NSAF

Consider a desired signald(n) that arises from the system identification model
(1)d(n)=u(n)w°+v(n),
where w° is a column vector for the impulse response of an unknown system that we wish to estimate, v(n) accounts for measurement noise with zero mean and variance σv2, and u(n) denotes the 1×M input vector,
(2)u(n)=[u(n)u(n-1)⋯u(n-M+1)].Figure1 shows the structure of the NSAF, where the desired signal d(n) and input signal u(n) are partitioned into N subbands by the analysis filters H0(z),H1(z),…,HN-1(z). The resulting subband signals, di(n) and yi(n) for i=0,1,…,N-1, are critically decimated to a lower sampling rate commensurate with their bandwidth. Here, the variable n to index the original sequences and k to index the decimated sequences are used for all signals. Then, the decimated filter output signal at each subband is defined as yi,D(k)=ui(k)w(k), where ui(k) is 1×M row vector, such that
(3)ui(k)=[ui(kN),ui(kN-1),…,ui(kN-M+1)],
and w(k)=[w0(k),w1(k),…,wM-1(k)]T denotes an estimate for w° with length M. Thus the decimated subband error signal is given by
(4)ei,D(k)=di,D(k)-yi,D(k)=di,D(k)-ui(k)w(k),
where di,D(k)=di(kN) is the decimated desired signal at each subband.Figure 1
Subband structure with the analysis filters and synthesis filters and the subband desired signals, subband filter outputs, and subband error signals.In [6], the authors have formulated the Lagrangian-based multiple-constraint optimization problem, which is formulated as
(5)JNSAF(k)=∥w(k+1)-w(k)∥2+∑i=0N-1‍λi[di,D(k)-ui(k)w(k+1)],
where λi for i=0,1,…,N-1 denote the Lagrange multipliers. Solving the cost function (5), the update recursion of the NSAF algorithm is given by [6, 7]. Consider
(6)w(k+1)=w(k)+μ∑i=0N-1‍uiT(k)∥ui(k)∥2ei,D(k),
where μ is the step-size parameter.


## 3. Weightedl1-Norm Regularized NSAF

### 3.1. Derivation of the Proposed Algorithm

To reflect the sparsity condition of the true system, that is,w°, a weighted l1-norm of the filter weight estimate is regularized on the cost function of the NSAF, which is given by
(7)Jl1-NSAF(k)=∥w(k+1)-w(k)∥2+∑i=0N-1‍λi[di,D(k)-ui(k)w(k+1)]+γ∥Πw(k+1)∥1,
where ∥Πw(k+1)∥ accounts for the weighted l1-norm of the filter weight vector w(k+1) and is written as
(8)∥Πw(k+1)∥1=∑m=0M-1‍πm|wm(k+1)|,
where Π is a M×M weighting matrix whose diagonal elements are πm and other elements are equal to zero, and wm(k+1) denotes the mth tap weight of w(k+1), for m=0,1,…,M-1. In addition, γ is a positive value parameter which plays a role in compromising the error related term and the weighted l1-norm regularization in the right-hand side of (7).To find the optimal weight vectorw(k+1) which minimizes the cost function (7), the derivative of (7) with respect to w(k+1) is taken and set to zero. Note that the weighted l1-norm regularization term, that is, ∥Πw(k+1)∥1, is not differentiable at any point in case wm(k+1)=0. To address this issue, a subgradient calculus [21] is carried out.Thus, taking the derivative of (7) with respect to the weight vector w(k+1) and letting the derivative be equal to zero, it leads to
(9)w(k+1)=w(k)+12∑i=0N-1‍λiuiT(k)-γ2∇ws∥Πw(k+1)∥1,
where ∇wsf(·) denotes a subgradient vector of a function f(·) with respect to w(k+1). An available subgradient vector ∇ws∥Πw(k+1)∥1 is obtained as [21]. Consider
(10)∇ws∥Πw(k+1)∥1=ΠTsgn(Πw(k+1))=Πsgn(w(k+1)),
since Π is assumed as a diagonal matrix with positive-valued elements, where sgn(·) is a componentwise sign function defined by
(11)sgn(x)={x|x|,x≠00,elsewhere.Substituting (10) into (9) and assuming sgn[w(k+1)]≈sgn[w(k)], it is given by
(12)w(k+1)=w(k)+12∑i=0N-1‍λiuiT(k)-γ2Πsgn(w(k)).Substituting (12) into the multiple constraints of the NSAF, that is, di,D(k)=ui(k)w(k+1), i=0,1,…,N-1 and rewriting as a matrix form, it leads to
(13)Λ=2[U(k)UT(k)]-1eD(k)+γ[U(k)UT(k)]-1U(k)Πsgn(w(k)),
where Λ=[λ0,λ1,…,λN-1]T is the N×1 Lagrange vector,
(14)U(k)=[u0(k)⋮uN-1(k)],eD(k)=[e0,D(k)⋮eN-1,D(k)].By neglecting the off-diagonal elements ofU(k)UT(k) [6], the components of Λ in (13) can be simplified to
(15)λi=2ei,D(k)∥ui(k)∥2+γui(k)∥ui(k)∥2Πsgn(w(k)),
for i=0,1,…,N-1.Consequently, combining (12) and (15), the update recursion of the sparsity-regularized NSAF is given by
(16)w(k+1)=w(k)+μ∑i=0N-1‍[uiT(k)∥ui(k)∥2ei,D(k)ffffffff+12γui(k)∥ui(k)∥2Πsgn(w(k))uiT(k)]-μγ2Πsgn(w(k)),
where μ is the step-size parameter.


### 3.2. Determination of the Weightedl1-Norm Regularization

Here, by choosing the weighting matrixΠ, two versions of the sparsity-regularized NSAF are developed. First, the use of the identity matrix as the weighting matrix, that is, Π=IM, results in the following update recursion:
(17)w(k+1)=w(k)+μ∑i=0N-1‍[uiT(k)∥ui(k)∥2ei,D(k)ffffffff+12γui(k)∥ui(k)∥2sgn(w(k))uiT(k)]-μγ2sgn(w(k)),
which is referred to as the l1-norm NSAF (l1-NSAF) as an unweighted case. The l1-NSAF uniformly attracts the tap coefficients of w(k) to zero. The zero attraction process leads to the improved convergence of the l1-NSAF when the majority of entries of a system are zero; that is, a system is sparse.Second, to approximate the actual sparsity condition of an underlying system, that is,l0-norm of the system, the weights of Π are chosen inversely proportional to magnitude of the actual coefficients of the system as given by
(18)πm={1|wm|,wm≠0∞,wm=0,
where wm denotes the mth coefficients of the system, that is, w°. However, since the actual coefficients of the system are unavailable, the estimate of the current filter weights is utilized instead of the actual weights, which is referred to as the reweighting scheme [22], as follows:
(19)πm(k)=1|wm(k)|+ϵform=0,1,…,M-1,
where wm(k) denotes the mth tap weight of the w(k) and ϵ is a small positive value to avoid singularity when |wm(k)|=0. Then, the weighting matrix Π consists of the values of πm(k) as the mth diagonal elements and has a time-varying feature. Finally, the update recursion is given by
(20)w(k+1)=w(k)+μ∑i=0N-1‍[uiT(k)∥ui(k)∥2ei,D(k)ffffffff+12γu¯i(k)∥ui(k)∥2sgn(w(k))uiT(k)]-μγ2sgn(w(k))|w(k)|+ϵ,
where u¯i(k)=ui(k)Π and the vector division operation in last term accounts for a componentwise division. Then, this recursion is called the reweighted l1-norm NSAF (l1-RNSAF)Table1 lists the number of multiplications and divisions of the NSAF [6], l1-NSAF, and l1-RNSAF per iteration. As shown in Table 1, the use of l1-norm constraint leads to an acceptable increase in computation.Table 1
Computational complexity.
NSAF
l
1-NSAF
l
1-RNSAF
Multiplications
3
M
+
3
N
L
6
M
+
3
N
L
7
M
+
3
N
L
Divisions
1
2
2
+
M
/
N
(M: filter length, N: number of subbands, and L: length of the analysis filters and synthesis filters).


## 3.1. Derivation of the Proposed Algorithm

To reflect the sparsity condition of the true system, that is,w°, a weighted l1-norm of the filter weight estimate is regularized on the cost function of the NSAF, which is given by
(7)Jl1-NSAF(k)=∥w(k+1)-w(k)∥2+∑i=0N-1‍λi[di,D(k)-ui(k)w(k+1)]+γ∥Πw(k+1)∥1,
where ∥Πw(k+1)∥ accounts for the weighted l1-norm of the filter weight vector w(k+1) and is written as
(8)∥Πw(k+1)∥1=∑m=0M-1‍πm|wm(k+1)|,
where Π is a M×M weighting matrix whose diagonal elements are πm and other elements are equal to zero, and wm(k+1) denotes the mth tap weight of w(k+1), for m=0,1,…,M-1. In addition, γ is a positive value parameter which plays a role in compromising the error related term and the weighted l1-norm regularization in the right-hand side of (7).To find the optimal weight vectorw(k+1) which minimizes the cost function (7), the derivative of (7) with respect to w(k+1) is taken and set to zero. Note that the weighted l1-norm regularization term, that is, ∥Πw(k+1)∥1, is not differentiable at any point in case wm(k+1)=0. To address this issue, a subgradient calculus [21] is carried out.Thus, taking the derivative of (7) with respect to the weight vector w(k+1) and letting the derivative be equal to zero, it leads to
(9)w(k+1)=w(k)+12∑i=0N-1‍λiuiT(k)-γ2∇ws∥Πw(k+1)∥1,
where ∇wsf(·) denotes a subgradient vector of a function f(·) with respect to w(k+1). An available subgradient vector ∇ws∥Πw(k+1)∥1 is obtained as [21]. Consider
(10)∇ws∥Πw(k+1)∥1=ΠTsgn(Πw(k+1))=Πsgn(w(k+1)),
since Π is assumed as a diagonal matrix with positive-valued elements, where sgn(·) is a componentwise sign function defined by
(11)sgn(x)={x|x|,x≠00,elsewhere.Substituting (10) into (9) and assuming sgn[w(k+1)]≈sgn[w(k)], it is given by
(12)w(k+1)=w(k)+12∑i=0N-1‍λiuiT(k)-γ2Πsgn(w(k)).Substituting (12) into the multiple constraints of the NSAF, that is, di,D(k)=ui(k)w(k+1), i=0,1,…,N-1 and rewriting as a matrix form, it leads to
(13)Λ=2[U(k)UT(k)]-1eD(k)+γ[U(k)UT(k)]-1U(k)Πsgn(w(k)),
where Λ=[λ0,λ1,…,λN-1]T is the N×1 Lagrange vector,
(14)U(k)=[u0(k)⋮uN-1(k)],eD(k)=[e0,D(k)⋮eN-1,D(k)].By neglecting the off-diagonal elements ofU(k)UT(k) [6], the components of Λ in (13) can be simplified to
(15)λi=2ei,D(k)∥ui(k)∥2+γui(k)∥ui(k)∥2Πsgn(w(k)),
for i=0,1,…,N-1.Consequently, combining (12) and (15), the update recursion of the sparsity-regularized NSAF is given by
(16)w(k+1)=w(k)+μ∑i=0N-1‍[uiT(k)∥ui(k)∥2ei,D(k)ffffffff+12γui(k)∥ui(k)∥2Πsgn(w(k))uiT(k)]-μγ2Πsgn(w(k)),
where μ is the step-size parameter.


## 3.2. Determination of the Weightedl1-Norm Regularization

Here, by choosing the weighting matrixΠ, two versions of the sparsity-regularized NSAF are developed. First, the use of the identity matrix as the weighting matrix, that is, Π=IM, results in the following update recursion:
(17)w(k+1)=w(k)+μ∑i=0N-1‍[uiT(k)∥ui(k)∥2ei,D(k)ffffffff+12γui(k)∥ui(k)∥2sgn(w(k))uiT(k)]-μγ2sgn(w(k)),
which is referred to as the l1-norm NSAF (l1-NSAF) as an unweighted case. The l1-NSAF uniformly attracts the tap coefficients of w(k) to zero. The zero attraction process leads to the improved convergence of the l1-NSAF when the majority of entries of a system are zero; that is, a system is sparse.Second, to approximate the actual sparsity condition of an underlying system, that is,l0-norm of the system, the weights of Π are chosen inversely proportional to magnitude of the actual coefficients of the system as given by
(18)πm={1|wm|,wm≠0∞,wm=0,
where wm denotes the mth coefficients of the system, that is, w°. However, since the actual coefficients of the system are unavailable, the estimate of the current filter weights is utilized instead of the actual weights, which is referred to as the reweighting scheme [22], as follows:
(19)πm(k)=1|wm(k)|+ϵform=0,1,…,M-1,
where wm(k) denotes the mth tap weight of the w(k) and ϵ is a small positive value to avoid singularity when |wm(k)|=0. Then, the weighting matrix Π consists of the values of πm(k) as the mth diagonal elements and has a time-varying feature. Finally, the update recursion is given by
(20)w(k+1)=w(k)+μ∑i=0N-1‍[uiT(k)∥ui(k)∥2ei,D(k)ffffffff+12γu¯i(k)∥ui(k)∥2sgn(w(k))uiT(k)]-μγ2sgn(w(k))|w(k)|+ϵ,
where u¯i(k)=ui(k)Π and the vector division operation in last term accounts for a componentwise division. Then, this recursion is called the reweighted l1-norm NSAF (l1-RNSAF)Table1 lists the number of multiplications and divisions of the NSAF [6], l1-NSAF, and l1-RNSAF per iteration. As shown in Table 1, the use of l1-norm constraint leads to an acceptable increase in computation.Table 1
Computational complexity.
NSAF
l
1-NSAF
l
1-RNSAF
Multiplications
3
M
+
3
N
L
6
M
+
3
N
L
7
M
+
3
N
L
Divisions
1
2
2
+
M
/
N
(M: filter length, N: number of subbands, and L: length of the analysis filters and synthesis filters).


## 4.  Numerical Results

The performance of the proposed sparsity-regularized NSAFs is validated by carrying out computer simulations in a system identification scenario in which the unknown channel is randomly generated. The lengths of the unknown system areM=128 and 512 in experiments where S of them are nonzero. The nonzero filter weights are positioned randomly and their values are taken from a Gaussian distribution 𝒩(0,1/S). Here, S=4 is used in the simulations except Figure 5 in which various S values are considered. The adaptive filter and the unknown system are assumed to have the same number of taps. The input signals are obtained by filtering a white, zero-mean, Gaussian random sequence through a first-order system G(z)=1/(1-0.9z-1). The signal-to-noise ratio (SNR) is calculated by
(21)SNR=10log10(E[y2(n)]E[v2(n)]),
where y(n)=u(n)w°. The measurement noise v(n) is added to y(n) such that SNR = 10, 20, and 30dB. In order to compare the convergence performance, the normalized mean square deviation (MSD),
(22)NormalizedMSD=E[∥w°-w(k)∥2∥w°∥2],
is taken and averaged over 50 independent trials. The cosine-modulated filter banks [23] with the subband number of N=4 are used in the simulations. The prototype filter of length L=32 is used. For comparison purpose, the proportionate NSAF (PNSAF) [12] is considered, which has been developed for sparse system identification. The step-size is set to μ=0.5 for SAF algorithms except the PNSAF where the step sizes μ=0.6 (Figure 2) and μ=0.65 (Figure 6) are chosen to achieve similar steady-state MSD with the l1-RNSAF for comparison purpose. For the l1-RNSAF, ϵ=0.01 is chosen. In addition, ρ=0.05 is used for the PNSAF. The γ values are obtained by repeated trials to minimize the steady-state MSD.Figure 2
Normalized MSD curves of the NLMS, NSAF, PNSAF,l1-NSAF, and l1-RNSAF (N=4).Figure2 shows the normalized MSD curves of the NLMS, NSAF, l1-NSAF, and l1-RNSAF, in cases of N=4 and SNR=30 dB. For the l1-NSAF and l1-RNSAF, γ=3×10-5 is chosen. As shown in Figure 2, the not only l1-RNSAF outperforms the conventional NLMS, NSAF, PNSAF, and l1-NSAF, but also the l1-NSAF has better performance than other conventional algorithms, in terms of the convergence rate and the steady-state misalignment.In Figure3, to verify the effect of γ on convergence performance, the normalized MSD curves of the l1-RNSAF for different γ values are illustrated, in case of N=4 and SNR=30 dB. For different γ values (γ=1×10-4,  1×10-5,  5×10-5, and 1×10-6), the l1-RNSAF is not excessively sensitive to γ. The analysis of an optimal γ value remains a future work.Figure 3
Normalized MSD curves of thel1-RNSAF for various γ values (N=4).Next, the performance of the proposedl1-norm regularized NSAFs is compared to the original NSAF under different SNR conditions. Figure 4 depicts the normalized MSD curves of the NSAF, l1-NSAF, and l1-RNSAF under SNR=10 and 20 dB, respectively. The γ value for the l1-NSAF and l1-RNSAF is set to 5×10-5. It is clear that both the l1-NSAF and l1-RNSAF are superior to the NSAF under several SNR cases. Furthermore, the l1-RNSAF performs well compared to l1-NSAF.Figure 4
Normalized MSD curves of the NSAF,l1-NSAF, and l1-RNSAF under various SNR conditions (SNR = 10, 20, and 30 dB, N=4).Figure 5
Normalized MSD curves of the NSAF andl1-RNSAF under various sparsity conditions (S=4,8,16, and 32, N=4).Figure 6
Normalized MSD curves of the NLMS, NSAF, PNSAF,l1-NSAF, and l1-RNSAF for long system of M=512 (N=4).In Figure5, the convergence properties of the NSAF and l1-RNSAF are compared under various sparsity conditions of an underlying system. With the same length of the system, that is, M=128, different sparsity conditions (S=4, 8, 16, and 32) are considered under SNR=30 dB. The value of γ is set to 3×10-5 for the l1-RNSAF. Figure 5 shows that the NSAF is independent of the sparsity condition. On the other hand, the results indicate that the more sparse the underlying system, the better the l1-RNSAF.The comparison of performance of the NSAF,l1-NSAF, and l1-RNSAF with a long system, here, the filter length M=512, is presented in Figure 6. For the l1-NSAF and l1-RNSAF, γ=5×10-5 is chosen. A similar result of Figure 2 is observed in Figure 6.Finally, the tracking capabilities of the algorithms of a sudden change in the system are tested forN=4 and SNR=30dB. Figure 7 shows the results when the unknown system is right-shifted for 20 taps. Same value of γ  in Figure 2 is used. As can be seen, the l1-NSAF and l1-RNSAF keep track of weight change without losing the convergence rate nor the steady-state misalignment compared to the conventional NLMS, NSAF, and PNSAF. To be specific, the l1-RNSAF achieves better performance than the l1-NSAF in terms of both convergence rate and steady-state misalignment.Figure 7
Normalized MSD curves of the NSAF, PNSAF,l1-NSAF, and l1-RNSAF in case of a time-varying unknown system (N=4). The system is right-shifted for 20 taps at 1500 iterations.


## 5. Conclusion

A new family of the NSAFs which takes into account the sparsity condition of an underlying system has been presented by incorporating a weightedl1-norm constraint of filter weights in the cost function. The update recursion is obtained by employing subgradient calculus on the weighted l1-norm constraint term. Subsequently, two sparsity regularized NSAFs, that is, the unweighted l1-NSAF and l1-RNSAF have been developed. The numerical results indicate that the proposed l1-NSAF and l1-RNSAF achieve highly improved convergence performance over the conventional algorithms for sparse system identification.



---

*Source: 601623-2013-12-12.xml*