Single-channel source separation problems arise when a number of sources emit signals that are mixed and recorded by a single sensor and one is interested in estimating the original sources of the signals based on the recorded mixture. This problem is ill-posed. Several different model channel source separation methods have been used. One method is the autoregressive model (AR). This model captures temporal correlations at the source, as was shown in [2,3]. In these papers it was proved that, for a single channel mixture of stationary (AR) sources, the (AR) coefficients can be uniquely identified and the sources separated. For non-stationary (AR) sources an adaptive sliding -window was introduced to update the process. A complete study of (AR) models can be found in [4,5].

Let $M(m,n)$ denote the collection of $m\times n$ matrices with nonnegative entries. For a given matrix $V\in M(m,n)$, the Nonnegative Matrix Factorization procedure ( $NMF$ ) is used to find matrices $W\in M(m,r)$ and $H\in M(r,n)$ such that $V=WH$. Matrix factorization has many applications. For example it is used in source separation and dimensionality reduction or clustering. In using (NMF) it is necessary to compute

An excellent survey on (NMF) and matrix factorization could be found in [15].
Other factorization methods used are Vector Quantization (VQ), Principal Component Analysis (PCA) and Independent Component Analysis (ICA). These can be written in the form $V\approx WH$. The differences between these methods and (NMF) are due to the different constraints placed on factoring matrices. In the (VQ) method the columns of $H$ are constrained to be unary vectors (i.e., all components are zero except for one element equal to 1 ). In the (PCA) procedure the columns of $W$ and the rows of $H$ must be orthogonal. In the (ICA) procedure the rows of $H$ are maximally statistically independent.

A major problem with the (PCA) procedure is that it allows the basis vectors to have both positive and negative components and the data are represented by linear combinations of these vectors. In some applications, the presence of negative

components contradicts the physical reality. For example, the pixels in a gray scale image must have nonnegative entries, so any image with negative intensities would have no reasonable interpretation. Also STFT magnitude is given by nonnegative quantities. The ( $NMF$ ) procedure was developed in an attempt to address this problem. (See $[1,6,7,11,18,19]$ ).

Two iterative (NMF) methods are in use - Alternative Least Squares (ALS) and that of Lee-Seung (LS), the reader may consult [8-13].

Consider two functions $F:\mathbb{R}^{M\times N}\rightarrow\mathbb{R}$ and $G:\mathbb{R}^{M\times N}\times\mathbb{R}^{M\times N}\rightarrow\mathbb{R}$.
Definition 1. [9] One says that $G\left(H,H^{\prime}\right)$ is an auxiliary function for $F(H)$ if the conditions

are satisfied.
In [9] the following quantity is considered as the cost function

where $1\leqslant i\leqslant m,1\leqslant a\leqslant r,V_{i}$ is the "ith", ( $1,n$ ) dimension row from $V,H_{a}$ is the " $a$ th", $(1,n)$ row from $H$ and $H=\left(H_{a}\right)_{1\leqslant a\leqslant r}$.
Lemma 2. [9] If $G$ is an auxiliary function, then $F$ is nonincreasing under the update

The Taylor expansion for $F$ in (3), leads to

Leaving both $H$ and $H_{n}$ as variables one obtains as an auxiliary function for $F$,

The quantity $\left(W_{n}^{T}W_{n}\right)$ is a positive semidefinite matrix. Moreover, as noted in [9], the difference matrix between,

and

remains positive semidefinite. (In Matlab, "./" denotes the pointwise division between matrices and "diag" the diagonal of a Matrix, seen as vector). Note that we need it twice, in order to keep the dimensionality right. Hence, a similar quantity as the above Taylor expansion, (4),

satisfies the conditions of an auxiliary function of the above $F$, see [9]. In order to make an upgrade each step for $H_{n+1}$, the arg min is involved. It is know that under appropriate conditions over $\Phi,\Psi$ the quadratic form $F$, (with $\Phi$ a positive defined matrix),

attains its minimum at

Thus, by setting (6) into (7),

we obtain $\bar{x}=H_{n+1}-H_{n}$ and therefore the new $H_{n+1}$ is

This leads to the Lee-Seung iterative (multiplicative) method.

More specific, for (8), the matrix $H_{n+1}$ is given by

Basically, the rule is to reduce the "distance" (or the cost function) by choosing at each step, appropriate $\eta_{n}$, respectively $\gamma_{n}$,

The heart of each iterative method consists in the choice of such $\eta_{n}$ and $\gamma_{n}$. Specifically, at step $n$, in [9], each matrix is given pointwisely by

where $(\cdot)_{ab}$, gives the location, (row and column), within the matrix. In (9) set(10). This leads to Hadamard or pointwise multiplication, denoted by ".". Eventually, the following iteration method is obtained

It was reported in [11,13], that the convergence result from [9], actually, does not provide enough conditions for the convergence of (11). A new more stable iteration based on (11)was presented and its convergence was study.

Remark 3. "Lin’s modification" for (11), from [11], consists in adding a "small" positive quantity (i.e. $\delta=10^{-9}$ ), such that the Lee-Seung iteration becomes:

For this new method, $\eta_{n}$ and $\gamma_{n}$, were set to be:

In (9) set

and alternatively,

The matrix V is generated by real audio data.

to obtain,

The problem with such iteration is that one or both of $\left(W_{n}^{T}W_{n}\right)^{-1}$ and $\left(H_{n+1}H_{n+1}^{T}\right)^{-1}$ can be negative. In order to solve this problem one can consider the projection onto the nonnegative orthant, denoted by $P_{+}[\cdot]$. The above Alternating Least Squares iteration becomes

The convergence for such iteration is described in [7]. We remark that a quantity such as $\left(A^{T}A\right)^{-1}A$ is similar to the projection operator obtained for the least squares method, (see for example [14]). This method was reported (see also our experiments) to be very fast but unstable, see also [17]. In [17], as well as here, the aim is to obtain a new iteration which has the speed convergence of (ALS) and the stability of (LS).

Typical convergence and divergence behaviors are indicated in Fig. 1. V is the spectrogram of an audio signal A3_whistle.wav ¹, the frame length of the FFT was set to 512 . Hence, the dimension of $\mathbf{V}$ was $m=512,n=174$. Both $\mathbf{W}$ and H were randomly initialized. The rank $r$ was set to 7 and $50\%$ overlap between the windows was used for generating the spectrogram. All the three algorithms were applied to decompose the music notes from the audio signal.

In (9) insert $\eta=\operatorname{in}v\left(W^{-1}W\right)$, to obtain the iteration

A "general" iteration method for such " $H_{n+1}$ " as in (9) would have the following structure

or,

where $K\left(X_{n}\right)$ is a surrogate for $A\left(W_{n}\right)$ which may be different at each step. We introduce $K$ to be a nonnegative matrix "close enough" to $A\left(W_{n}\right)$ so as to avoid the errors introduced by using a projection onto the positive orthant. Note that $Y_{n}$ is different from $W_{n}$ at each step. But, if $Y_{n}$ is close enough to $W_{n}$, then $X_{n}$ behaves analogously to $H_{n}$. We introduced within (16 and (17) two "general" iterations. The study will cover, in this matter, all possible choices for those $K,N,A$ and $B$ within ( 16 and (17). Appropriate settings for $K$ and $N$ from (16) will lead to (LS) iteration. Eventually we will compare it with (ALS), if $A$ and $B$ are well chosen within (17).

Remark 4. Set $A\left(W_{n}\right)^{-1}=\left(\operatorname{diag}\left(\operatorname{diag}\left(W_{n}^{T}W_{n}H_{n}^{T}\cdot/H_{n}^{T}\right)\right)\right)^{-1}$ and $B\left(H_{n+1}\right)^{-1}=\left(\left(\operatorname{diag}\left(\operatorname{diag}\left(W_{n}H_{n+1}H_{n+1}^{T}\cdot/W_{n}\right)\right)\right)\right)^{-1}$ to obtain the Lee-Seung iteration (11). As we can see at each step $A\left(W_{n}\right)$ is changing, even within (15); set $A\left(W_{n}\right)^{-1}=\left(W_{n}^{T}W_{n}\right)^{-1}$ and $B\left(H_{n+1}\right)^{-1}=\left(H_{n+1}H_{n+1}^{T}\right)^{-1}$, to obtain the (ALS) iteration.

Our main purposes are the following: first to be able to mix (LS) and (ALS) iterations in order to obtain a better one. Second, we show that "structurally" the two algorithms are not very different and we provide the mathematical background for such hybrid to converge.

Recall the following Lemma:
Lemma 5. [16] Let $\left\{a_{n}\right\}$ be a nonnegative sequence that satisfies

where $w\in(0,1)$ and $\mathbf{M}>0$ are fixed numbers and $\left\{\sigma_{n}\right\}$; is a nonnegative sequence which converges to zero. Then $\lim_{n\rightarrow\infty}a_{n}=0$.
The result remains true provided that the coefficient of $a_{n}$ stays within an interval in ( 0,1 ).
Proposition 6. Let $\left\{a_{n}\right\}$ be a nonnegative sequence that satisfies

where $\mathbf{M}>0$ is fixed number, $\left\{\lambda_{n}\right\}$ and $\left\{\sigma_{n}\right\}$; are nonnegative sequences such that $\left\{\lambda_{n}\right\}\subset(0,\Lambda)$, for some $\Lambda<1$ and $\lim_{n\rightarrow\infty}\sigma_{n}=0$. Then $\lim_{n\rightarrow\infty}a_{n}=0$.

Proof. Note that, for each $n\in\mathbb{N},\lambda_{n}\leqslant\Lambda$. By defining $(1-w)=\max_{n}\left(1-\lambda_{n}\right)$, the result follows from Lemma 5 .

Remark 7. If $\lambda_{n}>1$ for each $n$, the Proposition 6 fails. As an example, choose $a_{n}=n,\lambda_{n}=2$ and $\sigma_{n}=0$ for each $n$. Then it follows that $n+1=a_{n+1}\leqslant\lambda_{n}a_{n}+\sigma_{n}M=2n$. Proposition 6 remains true if $\lambda_{n}>1$, for only a finite subset of $\mathbb{N}$.

For sake of simplicity, through out this paper, we shall consider the sup - norm for all matrices involved. Within Matlab one use "max (max (…))" command.

Theorem 8. If iteration (17) converges i.e. $\lim_{n\rightarrow\infty}H_{n}=H^{*},\lim_{n\rightarrow\infty}W_{n}=W^{*}$ and $\lim_{n\rightarrow\infty}Y_{n}=W^{*}$, and there exists $\lambda\in(0,1)$ and $\mathbf{M}>0$ such that for each step $n$, we have the following relations satisfied

then iteration (16) is also convergent; i.e. $\lim_{n\rightarrow\infty}X_{n}=H^{*}$.

Proof. Define $M_{n}=Y_{n}-W_{n}$. Note that,

where $W_{n}$ is obtained from (17). Using (17) and (19) one obtains

Thus,

We shall consider here the sup - sup-norm such that $\left\|I_{r,r}\right\|=1$. Denote by

Note that one has $\lambda_{n}\in(0,\Lambda)$ and $\sigma_{n}\rightarrow 0$, therefore from Proposition 6, we obtain $\lim_{n\rightarrow\infty}a_{n}=\lim_{n\rightarrow\infty}\left\|H_{n}-X_{n}\right\|=0$, that is $\lim_{n\rightarrow\infty}\left\|X_{n}\right\|=H^{*}$. Using the inequality,

it follows that $\lim_{n\rightarrow\infty}X_{n}=H^{*}$.

Remark 9. Using Matlab it is easy to verify each step, of the condition of (18) by using (max (max (eye (r,r)-inv (diag (diag $\left.\left.\left.\left.\left(\left(A^{**}A^{*}B\right)./B\right)\right)\right)^{*}\left(A^{**}A\right)\right)\right)$, where size $(A)=(m,r)$ and size $(B)=(r,n)$. The second condition of (18), simply demands an upper bound for those matrices involved in the process.

In (17), set $A\left(W_{n}\right)=\left(W_{n}^{T}W_{n}\right)^{-1}$, (respectively, $B\left(H_{n+1}\right)=\left(H_{n+1}H_{n+1}^{T}\right)^{-1}$ ) to obtain the (ALS) method. The above Theorem leads to the following result.

Corollary 10. If the iteration (15) converges i.e. $\left(\lim_{n\rightarrow\infty}H_{n}=H^{*},\lim_{n\rightarrow\infty}W_{n}=W^{*}\right.$ and $\left.\lim_{n\rightarrow\infty}Y_{n}=W^{*}\right)$, and there exists a $\lambda\in(0,1)$ and an $\mathbf{M}>0$ such that for each step $n$, the following relations are satisfied

then the iteration (16) is also convergent; i.e. $\lim_{n\rightarrow\infty}X_{n}=H^{*}$.