Student Member, IEEE Jie Chen
email: chenjieg@sina.com

Jose Carlos

Senior Member, IEEE M Bermudez

Member, IEEE Paul Honeine
email: paul.honeine@utt.fr

C Richard
email: cedric.richard@unice.fr

J.-C M Bermudez
email: j.bermudez@ieee.org

Variants of non-negative least-mean-square algorithm and convergence analysis

Due to the inherent physical characteristics of systems under investigation, non-negativity is one of the most interesting constraints that can usually be imposed on the parameters to estimate. The Non-Negative Least-Mean-Square algorithm (NNLMS) was proposed to adaptively find solutions of a typical Wiener filtering problem but with the side constraint that the resulting weights need to be non-negative. It has been shown to have good convergence properties. Nevertheless, certain practical applications may benefit from the use of modified versions of this algorithm. In this paper, we derive three variants of NNLMS. Each variant aims at improving the NNLMS performance regarding one of the following aspects: sensitivity of input power, unbalance of convergence rates for different weights and computational cost. We study the stochastic behavior of the adaptive weights for these three new algorithms for nonstationary environments. This study leads to analytical models to predict the first and second order moment behaviors of the weights for Gaussian inputs. Simulation results are presented to illustrate the performance of the new algorithms and the accuracy of the derived models.

I. INTRODUCTION

Optimization of a cost function given a set of constraints is a common objective in signal processing estimation problems. The constraints are usually imposed by system specifications which provide a priori information on the feasible set of solutions. The solution of estimation problems under constraints poses special problems for online applications. Common realtime signal processing restrictions on computational complexity and memory requirements tend to rule out several good solutions to the constrained optimization problem.

Non-negativity is one of the most commonly stated constraints. It is often imposed on the parameters to estimate in order to avoid physically absurd and uninterpretable results. Non-negativity constraints have been used for image deblurring [START_REF] Benvenuto | Nonnegative leastsquares image deblurring: improved gradient projection approaches[END_REF], deconvolution of system impulse response estimation [START_REF] Lin | Bayesian regularization and nonnegative deconvolution for room impulse response estimation[END_REF] and audio processing [START_REF] Cont | Realtime multiple pitch and multipleinstrument recognition for music signals unsing sparse non-negative constraints[END_REF]. Another similar problem is the non-negative matrix factorization (NMF), which is now a popular dimension reduction technique used in many applications [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF]. This problem is closely related to blind deconvolution, and has found direct application in neuroscience [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF] and in hyperspectral imaging [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF]. Separation of non-negative mixture of non-negative sources has also been considered in [START_REF] Plumbley | Algorithms for nonnegative independent component analysis[END_REF], [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a bayesian approach and MCMC sampling[END_REF].

Over the last fifteen years, a variety of methods have been developed to tackle non-negative least-square (NNLS) problems. Active set techniques for NNLS use the fact that if the set of variables which activate constraints is known, then the solution of the constrained least-square problem can be obtained by solving an unconstrained one that includes only inactive variables. The active set algorithm of Lawson and Hanson [START_REF] Lawson | Solving Least Squares Problems[END_REF] is a batch resolution technique for NNLS problems. It has become a standard among the most frequently used methods. In [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF], Bro and De Jong introduced a modification of the latter, called Fast NNLS, which takes advantage of the special characteristics of iterative algorithms involving repeated use of non-negativity constraints. Projected gradient algorithms [START_REF] Rosen | The gradient projection method for nonlinear programming. part 1: Linear constraints[END_REF], [START_REF] Calamai | Projected gradient methods for linearly constrained problems[END_REF], [START_REF] Barzilai | Two-point step size gradient methods[END_REF], [START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF] form another class, which is based on successive projections onto the feasible region. In [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF], Lin used this kind of algorithm for NMF problems. Low memory requirements and simplicity make algorithms in this class attractive for large scale problems. Nevertheless, they are characterized by slow convergence rate if not combined with appropriate step size selection. The class of multiplicative algorithms is very popular for dealing with NMF problems [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], [START_REF] Lin | On the convergence of multiplicative update algorithms for nonnegative matrix factorization[END_REF]. Particularly efficient updates were derived in this way for a large number of problems involving non-negativity constraints [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. However, these algorithms require batch processing, which is not suitable for online system identification problems. In [START_REF] Chen | Non-negative least-mean square algorithm[END_REF], the problem of online system identification under non-negativity constraints on the parameters to estimate was investigated. An LMS-type adaptive algorithm, called Non-Negative Least-Mean-Square (NNLMS) was proposed to solve the Wiener problem under the constraint that the resulting weights need to be non-negative. It was based on the stochastic gradient descent approach combined with a fixed point iteration which converges to a solution satisfying the Karush-Kuhn-Tucker conditions. The stochastic behavior of this algorithm was also analyzed in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF].

In this paper, we extend the work of [START_REF] Chen | Non-negative least-mean square algorithm[END_REF] and derive useful variants of the NNLMS algorithm. Each of these variants is derived to improve the NNLMS properties in some sense. A normalized algorithm is proposed to reduce the NNLMS performance sensitivity to the input power value. An exponential algorithm is proposed to improve the balance of weight convergence rates. Compared to NNLMS, the new algorithm leads to faster convergence of the weights in the active set (weights for which the inequality constraint is satisfied with the equal sign). Finally, a sign-based algorithm is proposed to reduce implementation cost in critical real-time applications.

The rest of this paper is organized as follows. Section II reviews the system identification problem under non-negative constraints and the NNLMS algorithm. Section III motivates and introduces the NNLMS variants. In Section IV and Section V, the transient behavior of each of these variants is analyzed. Analytical models are derived for the mean weight and for the mean-square error behavior. The accuracy of these models is illustrated through simulations. A final example compares the performance of the proposed algorithms with those of NLMS and Projected NLMS in solving an unconstrained non-negative parameter estimation problem.

II. REVIEW OF NON-NEGATIVE LEAST-MEAN-SQUARE ALGORITHM

Consider the estimation problem depicted in Fig. 1. The unknown system is characterized by real-valued observations In certain applications, inherent physical characteristics of systems under investigation impose a non-negativity constraint on the estimate α of the system parameters. Therefore, the problem of identifying the optimum non-negative model can be formalized as follows

y(n) = α * x(n) + z(n), (1)
z(n) y(n) x(n) α(n) Algo. + + + - e(n) α *
α o = arg min α J(α) subject to α o i ≥ 0, ∀i, (2)
where J(α) is a continuously differentiable and strictly convex cost function in R N , and α o is the solution to the constrained optimization problem.

A. A fixed-point iteration scheme

To solve the problem (2), let us consider its Lagrangian function Q(α, λ) given by [START_REF] Boyd | Convex Optimization[END_REF]

Q(α, λ) = J(α) -λ α,
where λ is the vector of non-negative Lagrange multipliers. The Karush-Kuhn-Tucker conditions must necessarily be satisfied at the optimum defined by α o , λ o , namely,

∇ α Q(α o , λ o ) = 0 (3) α o i [λ o ] i = 0, i = 1, . . . , N (4)
where ∇ α stands for the gradient operator with respect to α.

Using ∇ α Q(α, λ) = ∇ α J(α)λ, these equations can be combined into the following expression

α o i [-∇ α J(α o )] i = 0, i = 1, . . . , N (5)
where the extra minus sign is just used to make a gradient descent of J(α) apparent. To solve equation ( 5) iteratively, two important points have to be noticed. First, D(-∇ α J(α)) is also a descent direction of J(α) if D is a symmetric positive definite matrix. Second, equations of the form ϕ(u) = 0 can be solved with a fixed-point iteration algorithm by considering the problem u = u + ϕ(u) under some conditions on function ϕ. Implementing this strategy with equation ( 5) leads to the component-wise gradient descent algorithm

α i (n + 1) = α i (n) + η i (n)f i (α(n)) α i (n)[-∇ α J(α(n))] i (6
) where η i (n) a positive step size required to get a contraction scheme and to control the convergence rate, and 6) is an arbitrary positive function of α(n). Some criteria J(α) are defined only for parameter vectors α with positive entries, e.g., the Itakura-Saito distance and the Kullback-Leibler divergence. If necessary, this condition can be managed by an appropriate choice of the step size parameter. Let us assume that α i (n) ≥ 0. Non-negativity of α i (n + 1) in ( 6) is guaranteed if

f i (α(n)) > 0 is the i-th entry of a diagonal matrix D(n). Function f i (α(n)) in (
1 + η i (n)f i (α(n))[-∇ α J(α(n)] i ≥ 0. (7)
If [∇ α J(α(n))] i ≤ 0, condition ( 7) is clearly satisfied and non-negativity does not impose any restriction on the step size.

Conversely, if [∇ α J(α(n))] i > 0, non-negativity of α i (n + 1) holds if 0 ≤ η i (n) ≤ 1 f i (α(n)) [∇ α J(α(n))] i . ( 8
)
Using a single step size η(n) in (0, η max (n)] for all entries of α so that

η max (n) = min i 1 f i (α(n)) [∇ α J(α(n))] i , i = 1, . . . , N (9)
the update equation ( 6) can be written in vector form as

α(n + 1) = α(n) + η(n) d(n) ( 10
)
where the ith entry of the weight adjustment direction d(n) defined as follows

[d(n)] i = f i (α(n)) α i (n)[-∇ α J(α(n))] i (11) is a descent direction of J(α) because f i [α(n)] α i (n) ≥ 0.
It should be noted that condition [START_REF] Plumbley | Algorithms for nonnegative independent component analysis[END_REF] on the step size η(n) guarantees the non-negativity of α(n) for all n, but does not ensure algorithm stability.

B. The Non-Negative Least-Mean-Square algorithm

Let us now consider the mean-square error criterion J mse (α) = E{[y(n)α x(n)] 2 } to be minimized with respect to α so that

α o = arg min α E{[y(n) -α x(n)] 2 } subject to α o i ≥ 0, ∀i (12)
where the non-negativity constraint applies only to the optimum solution because J mse (α) is defined for α with entries α i that can be positive or negative. The gradient of J mse (α) with respect to α is ∇ α J(α) = 2 (R x αr xy ), where R x is the autocorrelation matrix of x(n) and r xy is the correlation vector between x(n) and y(n). Following a stochastic gradient approach, the second-order moments R x and r xy are replaced in [START_REF] Lawson | Solving Least Squares Problems[END_REF] by the instantaneous estimates x(n) x (n) and y(n) x(n), respectively. Then, choosing f i (α(n)) = 1/2 for all i in (6), a fixed positive step-size η, defining

D x (n) = diag{x(n)} and D α (n) = diag{α(n)}, and noting that D α (n) x(n) = D x (n) α(n) leads to the stochastic update given by α(n + 1) = α(n) + η e(n) D x (n) α(n) (13)
where

e(n) = y(n) -α (n) x(n).
A detailed study of this algorithm, named Non-Negative LMS (NNLMS) can be found in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF].

III. VARIANTS OF THE NON-NEGATIVE LEAST-MEAN-SQUARE ALGORITHM

A. Normalized NNLMS

A direct extension of the original algorithm is the Normalized NNLMS. Conditioned on α(n), the product e(n)D x (n) in (13) has dimension of signal power. Thus, η is inversely proportional to signal power. Hence, setting a constant value for η leads to different weight updates for different signal power levels. This is the same sensitivity to signal power verified in the LMS algorithm. A popular way to address this limitation is to normalize the weight update by the input vector squared 2 -norm which yields the Normalized NNLMS update equation

α N (n + 1) = α N (n) + η x (n) x(n) e(n) D x (n) α N (n) (14)
Like in Normalized LMS (NLMS) algorithm, adding a small positive regularization parameter to the denominator x (n) x(n) may be necessary to avoid numerical difficulties when x (n) x(n) becomes very small. The resulting -Normalized NNLMS will then be

α N (n + 1) = α N (n) + η x (n) x(n) + e(n) D x (n) α N (n) (15
) where we maintained the notation α N (n) because ( 14) is a particular case of (15) for = 0. From now on, we refer to [START_REF] Barzilai | Two-point step size gradient methods[END_REF] simply as the Normalized NNLMS algorithm.

B. Exponential NNLMS

Each component α i (n) in the update term of (13) can be viewed as a distinct variable step size adjustment along the ith axis 1 . Hence, each component of α(n) will have a different convergence rate in general. Specifically in the case of weights in the active set (those that tend to zero in steady-state), the convergence rate will progressively reduce in time, becoming very small near steady-state. To alleviate this convergence rate unbalance, we introduce the Exponential NNLMS algorithm.

To achieve a faster convergence for the adaptive coefficients as they get close to zero we propose the use of

f i (α E (n)) = α γ-1
Ei (n) in [START_REF] Lawson | Solving Least Squares Problems[END_REF], with parameter γ chosen in order to attract small values of α Ei (n) towards zero. This leads to the ith weight update equation

α Ei (n + 1) = α Ei (n) + η e(n) x i (n) α γ Ei (n). (16)
For 0 < γ < 1, the ith weight update in [START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF] becomes larger than that in [START_REF] Rosen | The gradient projection method for nonlinear programming. part 1: Linear constraints[END_REF] when |α Ei (n)| < 1, thus accelerating convergence towards a null steady-state coefficient value.

The condition for α Ei (n + 1) ≥ 0 given α Ei (n) ≥ 0 can be easily determined from [START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF] as

η ≤ 1 e(n) α γ-1 Ei (n) , ∀ i, n. (17)
This condition, however, is not useful for design purposes, since it requires a priori knowledge of the algorithm behavior. We then propose a modified version of the update equation ( 16) that allows for instantaneous negative values of α Ei (n). The problem with real and negative instantaneous values of α Ei (n) is that it may lead to a complex value for α γ Ei for 0 < γ < 1. To obtain always real values for α γ Ei we propose to use γ = p/q with p and q odd integers and 0 < p < q. The oddness of p and q guarantees that sgn(α γ Ei (n)) = sgn(α Ei (n)). Then, the real solution for α γ Ei can be obtained by calculating sgn{α Ei } |α Ei | γ . This leads to the following weight update equation for the Exponential NNLMS algorithm in vector form:

α E (n + 1) = α E (n) + η e(n) D x (n) α (γ) E (n) (18a)
with the ith component of α

(γ) E (n) defined as [α (γ) E (n)] i = sign{α Ei (n)} |α Ei (n)| γ . ( 18b
)
As in the gamma correction used in image processing, an exponent in the range 0 < γ < 1 reduces the dynamic range of each α Ei (n). Large values of α Ei (n) will be compressed towards 1 and small values of α Ei (n) will be increased to prevent from stalling convergence. When γ = 1, the update equation degenerates into the NNLMS algorithm [START_REF] Rosen | The gradient projection method for nonlinear programming. part 1: Linear constraints[END_REF]. Using γ > 1 is generally not recommended, as it tends to spread the vector component values.

C. Sign-Sign NNLMS

Like Sign-Sign LMS, which has been included in the CCITT standard for adaptive differential pulse code modulation [START_REF]Itu-t recommendation g[END_REF], the motivation for introducing a Sign-Sign NNLMS algorithm is its computational simplicity and its robustness against disturbances [START_REF] Koike | Analysis of the sign-sign algorithm based on gaussian distributed tap weights[END_REF]. Replacing the input regressor vector and the estimation error in the update term with their signs reduces computation time and dynamic range requirements by replacing multiplications with shifts in real-time implementations. The Sign-Sign NNLMS algorithm is given by

α S (n + 1) = α S (n) + η sgn{e(n)} sgn{D x (n)} α S (n) (19)
After the two signs are evaluated, the ith component update is given by

α Si (n + 1) = α Si (n) ± η α Si (n) ( 20
)
where the sign before η is determined by sgn{e(n)x i (n)}. The step-size η is usually selected as a power of 2 -1 , say η = 2 -m for some integer m > 0. In this case, Equation ( 20) can be efficiently implemented using shift-add operations. Moreover, the non-negativity constraint will be always satisfied if α S is initialized with a positive vector and 0 < η < 1. Table I compares the computation complexities of NNLMS and its three variants described above. The rightmost column describes the anticipated property of each algorithm, to be verified in the following.

IV. MEAN WEIGHT BEHAVIOR

Convergence in the mean sense of the NNLMS algorithm [START_REF] Rosen | The gradient projection method for nonlinear programming. part 1: Linear constraints[END_REF] has been studied in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF] for a stationary environment. We now study the stochastic behavior of the NNLMS variants introduced in Section III for fixed step sizes and for a time variant unconstrained solution given by

α * (n) = α * o (n) + ξ(n) (21)
where α * o (n) is a deterministic time-variant mean and ξ(n) is zero-mean with covariance matrix Ξ = σ 2 ξ I and independent of any other signal. This simple model provides some information on how the performances of the proposed algorithms are affected by a time variant optimal solution which consists of a deterministic trajectory and a random perturbation. The analysis using more elaborate non-stationarity models such as the random walk model or the autoregressive model [START_REF] Sayed | Adaptive Filters[END_REF] leads to mathematically intractable situations. This is due to the extra multiplication of the weight update term by a function of α(n) in ( 13), ( 15), ( 18) and [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF], as compared to the LMS algorithm. For the random walk model, the recursive equation for the covariance matrix of the adaptive weight vector becomes a function of the optimal weight covariance matrix, which becomes unbounded as time progresses [START_REF] Sayed | Adaptive Filters[END_REF]. For the autoregressive model, a nonlinear term given by the product of the weight error vector and the optimal weight update makes it impossible to determine a recursive adaptive weight vector covariance matrix equation in the state-space form [START_REF] Solo | Adaptive Signal Processing Algorithms: Stability and Performance[END_REF]. The model ( 21) leads to a tractable analysis and still permits inferences about the behavior of the algorithms in randomly time variant environments by varying the power

σ 2 ξ of ξ(n).
Inferences on the ability of the algorithm to track mean weight variations are also possible but require a different model run for each type of mean time variation of α * o (n) to be investigated.

To conserve space and to simplify notation without ambiguity, from now on we use the generic notations α(n) and α i (n) whenever the given expression is valid for all the algorithms under study. Notations α N , α E and α S will be used only for expressions which are specific to the corresponding algorithm. The same notational observation applies to any vector or matrix when referring to a specific algorithm.

For the analyses that follow, we shall define the weight error vector with respect to the unconstrained solution α * (n) as

ṽ(n) = α(n) -α * (n). ( 22
)
and the weight error vector with respect to the mean uncon-

strained solution α * o (n) as v(n) = α(n) -α * o (n). (23)
The two vectors are related by ṽ

(n) = v(n) -ξ(n).

A. Statistical assumptions

The following analysis is performed for x(n) and z(n) zero-mean stationary Gaussian and for z(n) white and statistically independent of any other signal. We assume in the subsequent mean weight behavior analysis that the input and weight vectors are statistically independent, according to the Independence Assumption [START_REF] Sayed | Adaptive Filters[END_REF]. This assumption is typical in the study of adaptive algorithms. It is sometimes used for simplification and frequently required for mathematical tractability. The simulation results will show that the resulting analytical models have low sensitivity to this assumption, as they accurately predict the behavior of the three algorithms.

B. Normalized NNLMS algorithm

Using [START_REF]Itu-t recommendation g[END_REF] with the appropriate subscript in [START_REF] Barzilai | Two-point step size gradient methods[END_REF] and

e(n) = y(n) -α N (n) x(n) = z(n) -(v N (n) -ξ(n)) x(n) yields v N (n + 1) = v N (n) + η x (n) x(n) + z(n) D x (n) v N (n) + η x (n) x(n) + z(n) D x (n) α * o (n) - η x (n) x(n) + D x (n) α * o (n) x (n) v N (n) - η x (n) x(n) + D x (n) v N (n) v N (n) x(n) + η x (n) x(n) + D x (n) v N (n) ξ (n) x(n) + η x (n) x(n) + D x (n) α * o (n) ξ (n) x(n)-∆ N (n). ( 24
)
where

∆ N (n) = α * o (n + 1) -α * o (n)
is a deterministic vector proportional to the derivative of the mean unconstrained optimal solution.

Taking the expected value of ( 24) and noting that the expectations of the second, third, sixth and seventh terms on

E{v N (n + 1)} = E{v N (n)} -ηE 1 x (n) x(n) + D x (n) α * o (n) x (n) v N (n) -ηE 1 x (n)x(n) + D x (n)v N (n)v N (n)x(n) -∆ N (n). (25)
Using the independence assumption, the second expectation in the r.h.s. of ( 25) can be written as

E 1 x (n) x(n) + D x (n) α * o (n) x (n) v N (n) =D α * o (n) E x(n) x (n) x (n) x(n) + E{v N (n)} (26)
Evaluation of the first expected value in the r.h.s. of (26) requires approximations. Each numerator element is given by x(ni)x(nj). The random part of the denominator is given by

N -1 k=0 x 2 (n -k).
A common approximation that works well for reasonably large N is to neglect the correlation between these two variates, as the latter tends to vary much slower than the former [START_REF] Samson | Fixed point error analysis of the normalized ladder algorithms[END_REF], [START_REF] Almeida | A statistical analysis of the affine projection algorithm for unity step size and autoregressive inputs[END_REF]. Moreover, given its slow variation we approximate x (n) x(n) by its mean value N σ 2

x , which is reasonable for large values of N . Using these approximations yields

E 1 x (n) x(n) + D x (n) α * o (n) x (n) v N (n) ≈ 1 N σ 2 x + D α * o (n) R x E{v N (n)}. (27)
Using again

x (n) x(n) ≈ N σ 2
x and removing it from the expected value, the ith component of the second expectation in the r.h.s. of ( 25) is

D x (n) v N (n) v N (n) x(n) i = N j=1 x(n -i + 1) v Ni (n) v Nj (n) x(n -j + 1). ( 28
)
Taking the expectation, using the independence assumption and defining

K N (n) = E{v N (n) v N (n)} we obtain E{D x (n) v N (n) v N (n) x(n)} i = N j=1 r x (j -i) [K N (n)] ij =[R x K N (n)] ii (29)
which yields

E{D x (n) v N (n) v N (n) x(n)} = diag{R x K N (n)}
where diag{•} denotes the vector of diagonal entries in the matrix. Hence, [START_REF] Solo | Adaptive Signal Processing Algorithms: Stability and Performance[END_REF] becomes

E{v N (n + 1)} = I - η N σ 2 x + D α * o (n) R x E{v N (n)} - η N σ 2 x + diag{R x K N (n)}-∆ N (n). ( 30
)
This recursion for E{v N (n)} requires a model for K N (n). A recursive model will be derived for

K N (n) in Section V, see (44)
. That model can be used along with [START_REF] Price | A useful theorem for nonlinear devices having gaussian inputs[END_REF] to predict the mean weight behavior of the Normalized NNLMS algorithm. Nevertheless, we have found that a sufficiently accurate and more intuitive mean behavior model can be obtained by neglecting the weight error fluctuations and using the following separation approximation

K N (n) ≈ E{v N (n)} E{v N (n)}. (31)
This approximation has been successfully used in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF] to study the mean behavior of the NNLMS algorithm. A discussion about the validity of the approximation can be found in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF].

Extensive simulation results have shown that this approximation achieves adequate accuracy in modeling the mean behavior of the adaptive weights. We thus obtain the following model

E{v N (n + 1)} = I - η N σ 2 x + D α * o (n) R x E{v N (n)} - η N σ 2 x + diag R x E{v N (n)} E{v N (n)} -∆ N (n). (32)

C. Exponential NNLMS algorithm

Using [START_REF]Itu-t recommendation g[END_REF] with the appropriate subscript in [START_REF] Lin | On the convergence of multiplicative update algorithms for nonnegative matrix factorization[END_REF],

e(n) = y(n) -α E (n) x(n) = z(n) -(v E (n) -ξ(n)) x(n), and considering that sign{α Ei } |α Ei | γ is equal to the real solution of α γ
Ei , the Exponential NNLMS weight error update equation can be written as

v E (n + 1) =v E (n) + η e(n) D x (n) v E (n) + α * o (n) (γ) =v E (n) + η z(n) D x (n) v E (n) + α * o (n) (γ) -η D x (n) v E (n) + α * o (n) (γ) v E (n) x(n) +η ξ (n)x(n) v E (n) + α * o (n) (γ) -∆ E (n). (33)
where

(v E (n) + α * o (n) (γ) is a real vector.
The nonlinear term

(v E (n) + α * o (n)(n)) (γ)
on the r.h.s. complicates the evaluation of the expected value of (33) because the statistics of the weight error vector are unknown. We have again found out that using a zero-th order approximation of v E (n) is sufficient to provide a reasonably good model for the mean weight error behavior. Thus, we make

(v Ei (n) + α * oi (n)) γ ≈ (E{v Ei (n)} + α * oi (n)) γ . ( 34
)
Using (34) in [START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF], taking the expected value and considering the statistical properties of z(n) and ξ(n) yields

E{v E (n + 1)} ≈E{v E (n)} -η E{D x (n) E{v E (n)}+α * o (n) (γ) • x (n)v E (n)}-∆ E (n) = I N -η D r (n)R x E{v E (n)}-∆ E (n) (35)
where

I N is the N × N identity matrix and D r (n) is an N × N diagonal matrix defined as D r (n) = diag{r(n)} with r(n) being the N × 1 vector whose ith component is r i (n) = (E{v Ei (n)} + α * oi (n)) γ .
It is simple to verify that this model collapses to the NNLMS model derived in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF] for p = q = 1.

D. Sign-Sign NNLMS algorithm

The statistical analysis of the Sign-Sign NNLMS algorithm behavior is complicated by the fact that the weight update term is discontinuous in both the input vector x(n) and the error e(n) [START_REF] Eweda | Transient and tracking performance bounds of the sign-sign algorithm[END_REF]. To make that analysis tractable, we consider the case of input signal x(n) zero-mean and Gaussian [START_REF] Koike | Analysis of the sign-sign algorithm based on gaussian distributed tap weights[END_REF], [START_REF] Eweda | Transient and tracking performance bounds of the sign-sign algorithm[END_REF].

Using [START_REF]Itu-t recommendation g[END_REF] with the appropriate subscript in [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF] and

e(n) = y(n)-α S (n) x(n) = z(n)-(v S (n)-ξ(n))
x(n), the Sign-Sign NNLMS weight error update equation can be written as

v S (n + 1) =v S (n)+η sgn{z(n)-v S (n)x(n)+ξ (n)x(n)} • sgn{D x (n)} (v S (n) + α * o (n))-∆ S (n). ( 36
)
Note that, unlike the former two variants, the non-stationarity effect appears in the weight error update equation (36) as a nonlinear function of ξ(n). The ith component of (36) is given by

v Si (n + 1) = v Si (n) + η sgn z(n) -v S (n)x(n) + ξ (n)x(n) sgn x i (n) (v Si (n) + α * oi (n))-∆ Si (n). (37)
To determine the expected value of (37), we first note that it has been demonstrated in [START_REF] Eweda | Transient and tracking performance bounds of the sign-sign algorithm[END_REF], [START_REF] Dasgupta | Sign-sign LMS convergence with independent stochastic inputs[END_REF] using Price's theorem [START_REF] Price | A useful theorem for nonlinear devices having gaussian inputs[END_REF] that

E sgn{θ 1 } sgn{θ 2 } = 2 π sin -1 E{θ 1 θ 2 } σ θ1 σ θ2 (38)
for θ 1 and θ 2 two zero-mean jointly Gaussian variables with variances σ 2 θ1 and σ 2 θ2 , respectively. Then, noting that z(n)v S (n)x(n)+ξ (n)x(n) and x i (n) are zero-mean Gaussian when conditioned on v S (n) and ξ(n), we assume them jointly Gaussian 2 and use the result in (38) to obtain 2 As x(n) and z(n) are independent and both Gaussian,

[x(n), z(n)] is jointly Gaussian. When conditioned on v S (n) and ξ(n), [x (n), z(n) - v S (n)x(n) + ξ (n)x(n)] is jointly Gaussian as a linear transformation of [x (n), z(n)]. E sgn{z(n) -v S (n) x(n)+ξ (n)x(n)} • sgn{x i (n)}|v S (n), ξ(n) ≈ 2 π sin -1 - R i v S (n) -R i ξ(n) σ x σ e|vS(n),ξ(n) (39)
where R i the i-th column of R and σ 2 e|vS(n),ξ(n) is the variance of e(n) when conditioned on v S (n) and ξ(n) . Now, since sin -1 (•) is a nonlinear function and the distribution of its argument is unknown, we proceed as we did for the Exponential NNLMS algorithm and replace the nonlinear function by its zero-th order approximation

E sgn{z(n) -v S (n) x(n) + ξ (n)x(n)} sgn{x i (n)} ≈ 2 π sin -1 - R i E{v S (n)} σ x σ e|E{vS(n)},Ξ (40)
with

σ e|E{vS(n)},Ξ = σ 2 z + tr R x E{v S (n)}E{v S (n)} + trace{R x Ξ} (41)
Taking the expected value of (37), using the results (39) and (40) and expressing the result in vector form yields the mean weight error vector behavior model

E{v S (n + 1)} = I N + η D p (n) E{v S (n)} + η D p (n) α * o (n)-∆ S (n) (42)
where D p (n) is the N × N diagonal matrix D p (n) = diag{p(n)} with p(n) being the N × 1 vector whose ith entry is given by (40).

V. SECOND MOMENT ANALYSIS

We now study the behavior of the second-order moments of the adaptive weights for the three algorithms proposed in Section III. The analysis is performed under the same statistical hypotheses used in the previous section. The following additional assumptions are used in the subsequent analysis: A1 : The input vector x(n) is Gaussian. A2 : The weight error vector v(n) is statistically independent of x(n) x (n). The reasoning for this approximation has been discussed in detail in [START_REF] Minkoff | On the unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF]. A3 : The statistical dependence of v(n) v (n) and v(n) is neglected, following the same reasoning valid for A2, see [START_REF] Minkoff | On the unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF]. A4 : v(n) and (x (n) v(n)) 2 are statistically independent given A2. This is because (x (n) v(n)) 2 is a linear combination of the entries in v(n) v (n). Hence, this approximation follows the same reasoning discussed in [START_REF] Minkoff | On the unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF] to justify A2. These assumptions are typical in the study of adaptive algorithms. They are sometimes used for simplification and sometimes required for mathematical feasibility. The simulation results will show that these assumptions lead to analytical models which are accurate enough in predicting the behavior of the algorithms for design purposes.

The excess means square estimation error (EMSE) is given by ζ(n) = E ṽ (n)x(n)x (n)ṽ(n) . Using the relation between ṽ(n) and v(n), the properties of ξ(n), and noting from ( 24), ( 33) and (36) that v(n) and ξ(n) are independent, we can write ζ(n) as

ζ(n) = E (v(n) -ξ(n)) x(n)x (n)(v(n) -ξ(n)) = trace{R x K(n)} + trace{R x Ξ} (43) with K(n) = E v(n)v (n) .
The term trace{R x Ξ} is the contribution of the random non-stationarity of the system to the EMSE. In the following, we derive recursive models for K(n) for each of the algorithms.

A. Normalized NNLMS algorithm

Post-multiplying ( 24) by its transpose, taking the expectation, using the approximation x (n) x(n) ≈ N σ 2

x in (24), defining η = η/(N σ 2

x + ), using assumptions A1-A4 and proceeding as in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF] leads to

K N (n + 1) = K N (n) -η P 1N (n)K N (n) + K N (n)P 1N (n)+P 5N (n)+P 5N (n) + η2 P 6N (n) + P 7N (n) + P 7N (n) + P 8N (n) + η2 σ 2 z P 2N (n) + P 3N (n) + P 3N (n) + P 4N (n) + η2 (P 9N (n) + P 10N (n) + P 11N (n) + P 11N (n)) + K ∆N (n) (44)
with

P 1N (n) = E{D x (n) α * o (n) x(n)} = D α * o (n) R x . (45) P 2N (n) = E{D x (n) α * o (n) α * o (n) D x (n)} = D α * o (n) R x D α * o (n). (46)
P 3N (n) = E{D x (n) v N (n) α * o (n) D x (n)} ≈ E{D vN (n)} R x D α * o (n). (47)
P 4N (n) = E{D x (n) v N (n) v N (n) D x (n)} ≈ R x • K N (n) (48
) where • denotes the so-called Hadamard entry-wise product,

P 5N (n) = E{v N (n) x (n) v N (n) v N (n) D x (n)} ≈ K N (n) R x E{D vN (n)}. ( 49
) P 6N (n) =E{D x (n) α * o (n) x (n) v N (n) v N (n) x(n) α * o (n) D x (n)} =D α * o (n) Q N (n) D α * o (n). ( 50
)
where the matrix

Q N (n) is defined by Q N (n) = D α * o (n) (2 R x K N (n) R x + trace{R x K N (n)} R x ), P 7N (n) =E{D x (n) α * o (n)x (n)v N (n) x (n) v N (n) v N (n) D x (n)} ≈D α * o (n) Q N (n) E{D vN (n)}. ( 51
) P 8N (n) =E{D x (n) v N (n) v N (n) x(n) x (n) v N (n) v N (n) D x (n)} =Q N (n) • K N (n). (52)
P 9N (n) =E{ξ (n)x(n)ξ (n)x(n)D x (n) • v N (n)v N (n)D x (n)} (53)
P 10N (n) =E{ξ (n)x(n)ξ (n)x(n)D x (n) • α * o (n)α * o (n)D x (n)} (54)
and

P 11N (n) =E{ξ (n)x(n)ξ (n)x(n)D x (n) • v N (n)α * o (n)D x (n)} (55)
In obtaining (44), it was considered that the products of the last two term of ( 24) by the other terms lead to zero mean values due to the properties of ξ(n).

Expected values P 1N (n) through P 8N (n) correspond to the terms of the weight error vector recursive equation derived for the NNLMS algorithm in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF] with η substituted for η. Thus, we use the results from [START_REF] Chen | Non-negative least-mean square algorithm[END_REF] and indicate their values directly in (45) through (52). We now derive expressions for P 9N (n) through P 11N (n). These terms convey the effect of the random part of the environment non-stationarity.

Computing (i, j)th entry of P 9N (n) yields

[P 9N ] ij (n) = E k l ξ k (n)ξ l (n)x k (n)x l (n)x i (n)v Ni (n)v Nj (n)x j (n) = k l E {ξ k (n)ξ l (n)} E v Ni (n)v Nj (n) • E {x k (n)x l (n)x i (n)x j (n)} (56) As E {ξ k (n)ξ l (n)} = 0 only for k = l, [P 9N ] ij (n) = σ 2 ξ k [K N (n)] ij E x 2 k (n)x i (n)x j (n) . Using the Gaussian moment factorizing theorem yields k E x 2 k (n)x i (n)x j (n) = [R x ] ij [R x ] kk +[R x ] ik [R x ] jk = R x trace{R x } + 2 R x R x ij
. This enables us to write the result in matrix form

P 9N (n) = σ 2 ξ K N (n) • (R x trace{R x } + 2 R x R x ) (57)
Similarly, we have

P 10N (n) = σ 2 ξ α * o (n) α * o (n) • (R x trace{R x } + 2 R x R x ) (58)
P 11N (n) =σ 2 ξ E {v N (n)} α * o (n) • (R x trace{R x } + 2 R x R x ) (59)
The last term K ∆N (n) conveys the effect of deterministic variation of the mean of system weights. Observing the terms multiplied with ∆ N (n), we have

K ∆N (n) =∆ N (n)∆ N (n) -∆ N (n)(E {v N (n + 1)} + ∆ N (n)) -(E {v N (n + 1)} + ∆ N (n))∆ N (n) = -∆ N (n)∆ N (n) -∆ N (n)E {v N (n + 1)} -E {v N (n + 1)} ∆ N (n) (60)

B. Exponential NNLMS algorithm

The second order moment analysis of the Exponential NNLMS algorithm requires an improvement on approximation (34) for the nonlinearity [START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF]. We use instead the following first order approximation for the realvalued solution of (v Ei (n

(v E (n) + α * o (n)) (γ) in
) + α * oi (n)) γ : v Ei (n) + α * oi (n)) γ ≈(E {v Ei (n)} + α * oi (n)) γ + γ g(E {v Ei (n)}) • (E {v Ei (n)} + α * oi (n)) γ-1 (v Ei (n) -E {v Ei (n)}) (61)
where

g(E {v Ei (n)}) =1 -u E {v Ei (n)} + α * oi (n) + ε -u E {v Ei (n)} + α * oi (n) -ε (62)
with u(•) being the unit step function and ε a small constant. The reason to include the gate function

g about E {v Ei (n)} + α * oi (n) in the regular Taylor series is that the derivative of (v Ei (n) + α * oi (n)) γ tends to infinity if v Ei (n) approaches -α * oi (n). It is simple to verify that lim vE i (n)→-α * o i (n) g(E {v Ei (n)}) (E {v Ei (n)} + α * oi (n)) γ-1 = 0.
The zero-th order approximation is sufficient about the point where the function is equal to zero. With this new approximation, the term on v Ei (n) in (61) will include moments of the weight error vector which are necessary to proper modeling its fluctuations.

To use vector notation, we define two deterministic vectors r(n) and s(n) whose ith entries are respectively

r i =(E {v Ei (n)} + α * oi (n)) γ -γg(E {v Ei (n)})(E {v Ei (n)} + α * oi (n)) γ-1 E {v Ei (n)} s i =γ g(E {v Ei (n)}) (E {v Ei (n)} + α * oi (n)) γ-1
respectively. We define also the corresponding diagonal matrices D r (n) = diag{r(n)} and D s (n) = diag{s(n)}. With these new definitions, the linear approximation can be written in vector form as

(v E (n) + α * o (n) (γ) ≈ r(n) + D s (n) v E (n) (63)
Post-multiplying (33) by its transpose, using (63), taking the expected value, using assumptions A1-A4 and defining matrix D vE (n) = diag{v E (n)} yields, after simple algebraic manipulations as done in [START_REF] Chen | Non-negative least-mean square algorithm[END_REF]

K E (n + 1) =K E (n) -η P 1E (n) K E (n) + K E (n) P 1E (n) -η P 5E (n)K E (n) + K E (n)P 5E (n) + η 2 P 6E (n) + P 7E (n) + P 7E (n) + P 8E (n) + η 2 σ 2 z P 2E (n) + P 3E (n) + P 3E (n) + P 4E (n) + η 2 P 9E (n) + K ∆E (n) (64)
with the eight moments P 1E (n) -P 8E (n) given by

P 1E (n) = E D x (n) r(n) x (n) = D r (n) R x .
(65)

P 2E (n) = E{D x (n)r(n)r (n)D x (n)} ≈ D r (n)R x D r (n). (66) P 3E (n) = E{D x (n) r(n) v E (n) D s (n) D x (n)} ≈ D r (n) R x E{D vE (n)} D s (n) (67) P 4E (n) = E{D x (n) D s (n) v E (n) v E (n) D s (n) D x (n)} ≈ D s (n) R x • K E (n) D s (n). ( 68
) P 5E (n) = v E (n) x(n) D x (n) D s (n) = diag R x E{v E (n)} D s (n) (69) P 6E (n) =E{v E (n) x(n) D x (n) r(n) r (n) D x (n) x (n) v E (n)} ≈D r (n) Q E (n) D r (n). ( 70
)
where the matrix Q E (n) in above equations is defined by

Q E (n) = 2 R x K E (n) R x + trace{R x K E (n)} R x P 7E (n) = E{v E (n)x(n)D x (n)r(n)v E (n)D s (n)D x (n)x (n)v E (n)} ≈ D r (n) Q E (n) E{D vE (n)} D s (n). ( 71
)
and

P 8E (n) =E{v E (n) x(n) D x (n) D s (n)v E (n) v E (n) • D s (n) D x (n)x (n) v E (n)} ≈D s (n) Q E (n) • K E (n) D s (n). (72)
The expectation P 9E (n) conveys the non-stationarity effects and is given by

P 9E (n) =E (ξ (n)x(n)) 2 v E (n) + α * o (n) (γ) • v E (n) + α * o (n) (γ) (73)
Using the first order approximation (63) and simple manipulations yields

P 9E (n) = trace{ΞR x } r(n)r (n) + D s (n)K E (n)D s (n) + r(n)E v E (n) D s (n) + D s (n)E {v E (n)} r (n) (74)
The last term K ∆E (n) is obtained in the same form of (60)

K ∆E (n) = -∆ E (n)∆ E (n) -∆ E (n)E {v E (n + 1)} -E {v E (n + 1)} ∆ E (n) (75)

C. Sign-Sign NNLMS algorithm

Using the weight error vector definition v

S (n) = α S (n) -α * o (n) in (19) and sgn{e(n)}sgn{D x (n)} = sgn{D x (n)e(n)} yields v S (n + 1) =v S (n) + η sgn{D x (n) e(n)} v S (n) + η sgn{D x (n) e(n)} α * o (n)-∆ S (n) (76)
Post-multiplying (76) by its transpose, taking the expected value and rearranging the terms leads to

K S (n + 1) = K S (n) + η P 1S (n) + P 1S (n) + η P 2S (n) + P 2S (n) + η 2 P 3S (n) + P 4S (n) + P 4S (n) + P 5S (n) +K ∆S (n) ( 77
)
where

P 1S (n) = E{v S (n) α * o sgn{D x (n) e(n)}} (78) P 2S (n) = E{v S (n) v S (n) sgn{D x (n) e(n)}} (79) P 3S (n) = E{sgn{D x (n)e(n)}α * o (n)α * o • sgn{D x (n)e(n)}} (80) P 4S (n) = E{sgn{D x (n)e(n)}α * o (n)v S (n) • sgn{D x (n)e(n)}} (81) P 5S (n) = E{sgn{D x (n)e(n)}v(n)v S (n) • sgn{D x (n)e(n)}} (82)
These expected values are calculated in the following for x(n) Gaussian.

Expected value P 1S (n):

Using the properties of statistical expectation P 1S (n) can be written as

P 1S (n) =E v v S (n) α * o (n) E sgn{D x (n) e(n)|v S (n), ξ(n)} . (83)
The conditional expectation in (83) is given by (39), which must be approximated. Approximation (40) for the ith element of (39) is too simple to predict the weight error fluctuations. A more suitable approximation is given by a first order Taylor series expansion:

2 π sin -1 - R i v S (n) σ x σ e|vS(n),ξ(n) ≈ 2 π sin -1 - R i E{v S (n)} σ x σ e|E{vS(n)},Ξ - 2 π R i (v S (n) -E{v S (n)}) σ x σ e|E{vS(n)},Ξ 1 - R i E{v S (n)}) σ x σ e|E{vS(n)},Ξ 2 =q i (n) + s i (n)v S (n). (84)
where the scalar q i (n) and the vector s i (n) are deterministic variables defined respectively as

q i (n) = 2 π sin -1 - R i E{v S (n)} σ x σ e|E{vS(n)},Ξ + 2 π R i E{v S (n)} σ x σ e|E{vS(n)},Ξ) 1 - R i E{v S (n)}) σ x σ e|E{vS(n)},Ξ 2 (85)
s i (n) = - 2 π R i σ x σ e|E{vS(n)},Ξ 1 - R i E{v S (n)}) σ x σ e|E{vS(n)},Ξ 2
(86) with σ e|E{vS(n)},Ξ defined in (41).

Using (84) in (83) and defining the N × 1 vector q(n) with ith element q i (n) and the N ×N matrix S(n) with ith column s i (n), (83) becomes, after simple manipulations,

P 1S (n) ≈E{v S (n)} α * o (n) diag{q(n)} + K S (n) S(n) D α * o (n) (87)
Expected value P 2S (n): Similar to P 1S (n), we first express P 2S (n) in the form

P 2S (n) = E v {v S (n) v S (n) E{sgn{D x (n) e(n)|v S (n)}}} (88
) Then, using (39) and (84) we obtain

P 2S (n) ≈K S (n) diag{q(n)} + E{v S (n) v S (n) diag{S (n)v S (n)}}. (89)
The (i, j)th element of the expectation in (89) is given by

E{[v S (n)v S (n)diag{S (n)v(n)}] ij } = N k=1 S kj E{v S k (n)v Si (n)v Sj (n)}. (90)
Evaluation of the third order moment in (90) requires further approximation, as the distribution of v S (n) is unknown. We assume that the distribution of v(n) can be approximated by a Gaussian distribution about its mean value. Then, using the properties of Gaussian variables [START_REF] Miller | Multidimensional Gaussian Distributions[END_REF] and defining the centered variable vp (n) = v p (n) -E{v p (n)} for p = i, j, k, we have

E{v i (n)v j (n)v k (n)} ≈E{v i (n)v j (n)} E{v k (n)} + E{v i (n)v k (n)} E{v j (n)} +E{v j (n)v k (n)} E{v i (n)} + E{v i (n)}E{v j (n)}E{v k (n)} =K ij (n) E{v k (n)} + K ik (n) E{v j (n)} + K jk (n) E{v i (n)} -2 E{v i (n)}E{v j (n)}E{v k (n)} (91)
which completes the derivation of (89).

Expected value P 3S :

The (i, j)-th entry of matrix P 3S (n) is given by

[P 3S (n)] ij =E{sgn{x i (n) e(n)} [α * o (n) α * o (n)] ij sgn{x j (n) e(n)} } =E{[α * o (n) α * o (n)] ij sgn{x i (n) x j (n)}} (92) Using (38), E{sgn{x i (n) x j (n)}} = (2/π) sin -1 ([R x ] ij /σ 2 x
) and

[P 3S (n)] ij = [α * o (n) α * o (n)] ij 2 π sin -1 [R x ] ij σ 2 x (93)
Finally, expressing the result in the matrix form yields

P 3S (n) = 2 π (α * o (n) α * o (n)) • T (94)
where the (i, j)th elements of the N × N matrix T is given by Expected values P 4S (n) and P 5S (n): Using the same reasoning and approximations used to evaluate P 1S (n) to P 3S (n) yields

[T ] ij = sin -1 ([R x ] ij /σ 2 x ).
P 4S (n) = 2 π (α * o (n) E{v S (n)}) • T (95)
P 5S (n) = 2 π (K S (n)) • T (96)
The last term K ∆S (n) writes

K ∆S (n) = -∆ S (n)∆ S (n) -∆ S (n)E {v S (n + 1)} -E {v S (n + 1)} ∆ S (n) (97)

VI. STEADY-STATE BEHAVIOR

The recursive transient models derived in Section IV and in Section V are nonlinear in the weight error correlation matrix. Hence, it is not possible to derive closed form expressions for their steady-state values. Nevertheless, these models can also provide information about the steady-state behavior of the algorithms in the limit as n → ∞. Fig. 2 shows examples of plots of the steady-state EMSE for the three new algorithms as functions of both the step size and the nonstationarity parameter σ 2 ξ . These plots were obtained from the theoretical models for an unknown response of order 10 with samples drawn from a uniform distribution in [0,1].

VII. SIMULATION RESULTS AND DISCUSSION

We now present simulation examples to illustrate the properties of the three algorithms and the accuracy of the derived models. The parameters for these examples were chosen to illustrate several properties of the three algorithms while conserving space. Similar results have been obtained using a variety of parameter sets. For all examples, N = 31. The unknown stationary system is defined as α * (stat.) oi = 0.9 -0.05 i, i = 0, . . . , 18 -0.01 (i -18) i = 19, . . . , 31

For the non-stationary case, we consider an unknown response defined by

α * (nonstat.) oi (n) =α * (stat.) oi + |α * (stat.) oi | 10 sin 2 π T n + 2 π i -1 N + ξ i (n) (99)
where the period T of the deterministic sinusoidal component was set to 2500. ξ(n) is a zero-mean Gaussian random vector with correlation matrix σ 2 ξ I with σ 2 ξ = 5 × 10 -4 . The input signal is a first-order AR process given by x(n) = 0.5 x(n -1) + w(n), with w(n) i.i.d. zero-mean Gaussian with variance σ 2 w , adjusted to obtain the desired input power σ 2

x = 1. The noise z(n) is zero-mean i.i.d. Gaussian with variance σ 2 z = 10 -2 . The adaptive weights in α i (0) were all initialized at 10/N for all realizations. The step size was always set to η = 0.005 for all but the normalized variant. For the latter we used η = 0.005 N σ 2

x , which leads to an equivalent step size η = 0.005. Monte Carlo simulations were obtained by averaging 100 runs.

1) Example 1: Figs. 3 and4 show the results for the Normalized NNLMS. The parameter was set to 0. Blue curves show simulation results and red curves show the theoretical predictions from [START_REF] Miller | Multidimensional Gaussian Distributions[END_REF]. Fig. 3 is for σ 2

x = 1 and Fig. 4 is for σ 2 x = 0.5. It can be verified that the model [START_REF] Miller | Multidimensional Gaussian Distributions[END_REF] accurately predicts the algorithm behavior, and that normalization has made the algorithm performance basically independent of the input power.

2) Example 2: Fig. 5 illustrates the results for the Exponential NNLMS algorithm. The parameter (p, q) = (5, 7) was used. Compared with Fig. 3, these figures clearly show that the coefficients that tend to zero in steady-state had their convergence rate significantly improved by the Exponential NNLMS algorithm. Also, the accuracy of the theoretical model (35) can be verified.

3) Example 3: Fig. 6 illustrate the result of the Sign-Sign NNLMS under stationary and nonstationary environment. These figures illustrate the accuracy of the model (42). It is also clear that the Sign-Sign NNLMS coefficients converge much slower than those for the NNLMS algorithm as expected.

A. Second moment behavior

We now illustrate the EMSE behavior of the NNLMS variants. Again, blue curves were obtained from Monte Carlo simulation and red curves show the theoretical predictions.

1) Example 4: Fig. 7 illustrates the EMSE behavior of the Normalized NNLMS algorithm. The accuracy of model (44) can be easily verified. Two more curves are added to each plot to illustrate the effect of the random non-stationarity parameter σ 2 ξ . Random perturbations with different variances were also added to α * (stat.) oi (Fig. 7(a)) and to the nonstationary case (Fig. 7(b)). The light blue (dash-dot) lines show the theoretical EMSE for σ 2 ξ = 10 -3 , while the green (dash) lines show the theoretical EMSE for σ 2 ξ = 5 × 10 -3 . These curves illustrate the expected extra EMSE due to tracking of the random optimal solution variations. Simulation curves coincide with the theoretical ones, but are not shown to preserve the visibility of the other curves.

2) Example 5: Fig. 8 illustrates the EMSE behavior of the Exponential NNLMS algorithm. The blue and red curves show again the simulation results and the accurate theoretical predictions using (64) for (p, q) = (5, 7). The light blue (dash-dot) and the green (dash) curves show the theoretical predictions of the EMSE behavior for (p, q) = (3, 5) and (p, q) = (1, 1) (original NNLMS), respectively. The simulation results agree with these curves but are not shown for clarity. These results confirm that the Exponential NNLMS algorithm accelerates the convergence of the adaptive weights when compared to NNLMS.

3) Example 6: Fig. 9 illustrates EMSE behavior of the Sign-Sign NNLMS algorithm. Once more, the red curves and blue curves illustrate the accuracy of the model (77). The green (dashed) curves show the performance of the original NNLMS in the same conditions. These curves illustrate the slower convergence rate of Sign-Sign NNLMS when compared to NNLMS, the price paid for a reduced computational complexity.

B. A comparative example

This example compares the performance of the NNLMS algorithm and its variants with that of unconstrained algorithms in solving the unconstrained solution problem of identifying an unknown weight vector α * with non-negative coefficients. This is an interesting application, as in this case the unconstrained algorithm will converge in the mean to the optimal solution. Though the problem description may guarantee that the optimal weights are positive, often in practice one do not have accurate information about the number of coefficients in the optimal solution. A common approach is to set the adaptive filter with a sufficient number of coefficients, usually larger than the actual unknown number. This examples illustrates the performance of the different algorithms in this case.

Consider a non-negative unknown optimal solution with i = 1, . . . , 10 and adaptive filters with N = 30 coefficients. Five algorithms were tested: NLMS [START_REF] Sayed | Adaptive Filters[END_REF], Projected Gradient NLMS [START_REF] Theodoridis | Adaptive learning in a world of projections[END_REF], Normalized NNLMS, Exponential NNLMS and Sign-Sign NNLMS. In Projected Gradient NLMS, the coefficients which activate the non-negativity constraints are projected into the feasible region, i.e. set to 0, at each iteration. The input signal was given by x(n) = 0.5 x(n -1) + w(n) with σ 2 w = 3/4 so that σ 2 x = 1. The initial weights α(0) were drawn from the uniform distribution U(0, 1). The additive noise z(n) was i.i.d. Gaussian with σ 2 z = 0.1. The step sizes were chosen for each algorithm by experimentation so that all would reach approximately the same steady-state EMSE with the value of 2 × 10 -3 . The step sizes were η = 0.035 for both NLMS and Projected Gradient NLMS, η = 0.8750 for Normalized NNLMS, η = 0.022 for Exponential NNLMS and η = 0.007 for Sign-Sign NNLMS. Figure 10 shows the EMSE evolution for the five algorithms (Monte Carlo simulation averaged over 100 realizations). Figure 11 shows the estimated weights for a single realization of the input signal at n = 15000. Although the unconstrained NLMS algorithm is able to converge to the optimal solution in the mean sense, it does not provide a good estimation of the zero-valued coefficients in a single realization. NNLMS-type algorithms, including the Sign-Sign algorithm (which has not even converged to the steady-state at n = 15000) do a better job in determining the support of the actual response.

α * i = exp(-0.6 i) (100)

VIII. CONCLUSION

Many real-life systems require non-negative coefficients when their physical behavior is parameterized. In such cases, a non-negativity constraint should be imposed on the parameters to estimate in order to avoid physically absurd and uninterpretable results. The Non-Negative Least-Mean-Square (NNLMS) algorithm has been recently proposed to solve such a constrained Wiener problem online. In this paper, we proposed three variants of NNLMS, each addressing a different issue that may affect NNLMS under given circumstances. The performances of the Normalized NNLMS, Exponential NNLMS and Sign-Sign NNLMS algorithms were studied for nonstationary environments. The optimal unconstrained solution was modeled by a time-variant mean plus a random fluctuation. The derived analytical models were shown to accurately predict both the mean and the mean-square behavior of the algorithms. Their performances were compared and their advantages in potential applications discussed. Future research efforts will further explore these NNLMS variants properties and apply them in practical situations where efficient adaptive solutions to non-negative filtering problem are required.

Fig. 1 .

 1 Fig. 1. Adaptive system under study

Fig. 2 .

 2 Fig. 2. Steady-State EMSE behavior for the Normalized NNLMS, Exponential NNLMS and Sign-Sign NNLMS algorithms.

1 Fig. 3 .

 13 Fig.3. Evolution of the coefficients α i (n) for the normalized NNLMS algorithm in stationary and nonstationary environments for σ 2x = 1. Theoretical curves are from[START_REF] Miller | Multidimensional Gaussian Distributions[END_REF].

5 Fig. 4 .

 54 Fig.[START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. Evolution of the coefficients α i (n) for the normalized NNLMS algorithm in stationary and nonstationary environments for σ 2 x = 0.5. Theoretical curves are from[START_REF] Miller | Multidimensional Gaussian Distributions[END_REF].

Fig. 5 .

 5 Fig. 5. Evolution of the coefficients α i (n) for the Exponential NNLMS algorithm in stationary and nonstationary environments. Theoretical curves are from (35).

Fig. 6 .

 6 Fig.6. Evolution of the coefficients α i (n) for the Sign-Sign NNLMS algorithm in stationary and nonstationary environments. Theoretical curves are from (42).

Fig. 7 .Fig. 8 .

 78 Fig. 7. Evolution of the EMSE for the Normalized NNLMS algorithm in stationary and nonstationary environment. Light blue dash-dot line and green dashed line show the theoretical results for σ 2 ξ = 10 -3 and 5 × 10 -3 , respectively.

Fig. 9 .

 9 Fig. 9. Evolution of the EMSE for the Sign-Sign NNLMS algorithm in stationary and nonstationary environment. Theoretical evolution of original NNLMS is represented by the green dashed line.

Fig. 10 .

 10 Fig. 10. EMSE (db) for the four algorithms compared.

  Sign-Sign NNLMS algorithm.

Fig. 11 .

 11 Fig. 11. Weights estimated by NLMS, Projected NLMS, Normalized NNLMS, and Exponential NNLMS at n = 15000 for a single realization. Real weights are marked by ×. The NNLMS variants determine clearly the support of the response.

Note that Eq. (13) can be written as α(n + 1) = α(n) + ηDα(n)x(n)e(n). Hence, Dα(n) multiplies the ith gradient component e(n)x(n -i + 1) of the classical LMS algorithm by α i (n).

This work has been partly supported by CNPq grant No. 305377/2009-4.