Fast and Robust Recursive Algorithms for Separable Nonnegative Matrix Factorization

Introduction

A hyperspectral image consists of a set of images taken at different wavelengths. It is acquired by measuring the spectral signature of each pixel present in the scene, that is, by measuring the reflectance (the fraction of the incident electromagnetic power that is reflected by a surface at a given wavelength) of each pixel at different wavelengths. One of the most important tasks in hyperspectral imaging is called unmixing. It requires the identification of the constitutive materials present in the image and estimation of their abundances in each pixel. The most widely used model is the linear mixing model: the spectral signature of each pixel results from the additive linear combination of the spectral signatures of the constitutive materials, called endmembers, where the weights of the linear combination correspond to the abundances of the different endmembers in that pixel.

where $h_{kj}\geq 0$ is the abundance of the $k$th endmember in the $j$th pixel, with $\sum_{k=1}^{r}h_{kj}=1$ $\forall j$. Defining the $m$-by-$r$ matrix $W=[w_{1}\,w_{2}\,\dots w_{k}]\geq 0$ and the $r$-by-$n$ matrix $H$ with $H_{kj}=h_{kj}$ $\forall j,k$, the equation above can be equivalently written as $M=WH$ where $M$, $W$ and $H$ are nonnegative matrices. Given the nonnegative matrix $M$, hyperspectral unmixing amounts to recovery of the endmember matrix $W$ and the abundance matrix $H$. This inverse problem corresponds to the nonnegative matrix factorization problem (NMF), which is a difficult and highly ill-posed problem .

However, if we assume that, for each constitutive material, there exists at least one pixel containing only that material (a ‘pure’ pixel), then the unmixing problem can be solved in polynomial time: it simply reduces to identifying the vertices of the convex hull of a set of points. This assumption, referred to as the pure-pixel assumption , is essentially equivalent to the separability assumption : a nonnegative matrix $M$ is called separable if it can be written as $M=WH$ where each column of $W$ is equal, up to a scaling factor, to a column of $M$. In other words, there exists a cone spanned by a small subset of the columns of $M$ containing all columns (see Section 2.1 for more details). It is worth noting that this assumption also makes sense for other applications. For example, in text mining, each entry $M_{ij}$ of matrix $M$ indicates the ‘importance’ of word $i$ in document $j$ (e.g., the number of appearances of word $i$ in text $j$). The factors $(W,H)$ can then be interpreted as follows: the columns of $W$ represent the topics (i.e., bags of words) while the columns of $H$ link the documents to these topics. Therefore,

Separability of $M$ (that is, each column of $W$ appears as a column of $M$) requires that, for each topic, there exists at least one document discussing only that topic (a ‘pure’ document).

Separability of $M^{T}$ (that is, each row of $H$ appears as a row of $M$) requires that, for each topic, there exists at least one word used only by that topic (a ‘pure’ word).

These assumptions often make sense in practice and are actually part of several existing document generative models, see and the references therein.

We focus in this paper on hyperspectral unmixing algorithms under the linear mixing model and the pure-pixel assumption, or, equivalently, to nonnegative matrix factorization algorithms under the separability assumption. Many algorithms handling this situation have been developed by the remote sensing community, see for a comprehensive overview of recent hyperspectral unmixing algorithms. Essentially, these algorithms amount to identifying the vertices of the convex hull of the (normalized) columns of $M$, or, equivalently, the extreme rays of the convex cone generated by the columns of $M$. However, as far as we know, none of these algorithms have been proved to work when the input data matrix $M$ is only approximately separable (that is, the original separable matrix is perturbed with some noise), and many algorithms are therefore not robust to noise. However, there exists a few recent notable exceptions:

Arora et al. [3, Section 5] proposed a method which requires the resolution of $n$ linear programs in $\mathcal{O}(n)$ variables ($n$ is the number of columns of the input matrix), and is therefore not suited to dealing with large-scale real-world problems. In particular, in hyperspectral imaging, $n$ corresponds to the number of pixels in the image and is of the order of $10^{6}$. Moreover, it needs several parameters to be estimated a priori (the noise level, and a function of the columns of $W$; see Section 2.4).

Esser et al. proposed a convex model with $n^{2}$ variables (see also where a similar approach is presented), which is computationally expensive. In order to deal with a large-scale real-world hyperspectral unmixing problem, the authors had to use a preprocessing, namely $k$-means, to select a subset of the columns in order to reduce the dimension $n$ of the input matrix. Their technique also requires a parameter to be chosen in advance (either the noise level, or a penalty parameter balancing the importance between the approximation error and the number of endmembers to be extracted), only applies to a restricted noise model, and cannot deal with repeated columns of $W$ in the data set (i.e., repeated endmembers).

Bittorf et al. proposed a method based on the resolution of a single convex optimization problem in $n^{2}$ variables (cf. Section 5.2). In order to deal with large-scale problems ($m~{}\sim~{}10^{6}$, $n~{}\sim~{}10^{5}$), a fast incremental gradient descent algorithm using a parallel architecture is implemented. However, the algorithm requires several parameters to be tuned, and the factorization rank has to be chosen a priori. Moreover, it would be impractical for huge-scale problems (for example for web-related applications where $n\sim 10^{9}$), and the speed of convergence could be an issue.

2 Contribution and Outline of the Paper

In this paper, we propose a new family of recursive algorithms for nonnegative matrix factorization under the separability assumption. They have the following features:

They are very fast, running in approximately $6mnr$ floating point operations, while the memory requirement is low, as only one $m$-by-$n$ matrix has to be stored.

They are extremely simple to implement and would be easily parallelized.

They do not require any parameter to be chosen a priori, nor to be tuned.

The solution does not need to be recomputed from scratch when the factorization rank is modified, as the algorithms are recursive.

A simple post-processing strategy allows us to identify outliers (Section 3).

Even if the input data matrix $M$ is not approximately separable, they identify $r$ columns of $M$ whose convex hull has large volume (Section 4.1).

To the best of our knowledge, no other algorithms share all these desirable properties. The weak point of our approach is that the bound on the noise to guarantee recovery is weaker than in ; see Section 2.4. Also, we will need to assume that the matrix $W$ is full rank, which is not a necessary condition for the approaches above . However, in practice, this condition is satisfied in most cases. At least, it is always assumed to hold in hyperspectral imaging and text mining applications, otherwise the abundance matrix $H$ is typically not uniquely determined; see Section 2.1. Moreover, in Section 5.2, our approach will be shown to perform in average better than the one proposed in on several synthetic data sets.

The paper is organized as follows. In Section 2, we introduce our approach and derive an a priori bound on the noise to guarantee the recovery of the pure pixels. In Section 3, we propose a simple way to handle outliers. In Section 4, we show that this family of algorithms generalizes several hyperspectral unmixing algorithms, including the successive projection algorithm (SPA) , the automatic target generation process (ATGP) , the successive volume maximization algorithm (SVMAX) , and the $p$-norm based pure pixel algorithm (TRI-P) . Therefore, our analysis gives the first theoretical justification of the better performances of this family of algorithms compared to algorithms based on locating pure pixels using linear functions (such as the widely used PPI and VCA algorithms) which are not robust to noise. This was, until now, only experimentally observed. Finally, we illustrate these theoretical results on several synthetic data sets in Section 5.

Robust Recursive NMF Algorithm under Separability

In this section, we analyze a family of simple recursive algorithms for NMF under the separability assumption; see Algorithm 1.

Given an input data matrix $M$ and a function $f$, it works as follows: at each step, the column of $M$ maximizing the function $f$ is selected, and $M$ is updated by projecting each column onto the orthogonal complement of the selected column.

Instead of fixing a priori the number $r$ of columns of the input matrix to be extracted, it is also possible to stop the algorithm whenever the norm of the residual (or of the last extracted column) is smaller than some specified threshold.

In Section 2.1, we discuss the assumptions on the input separable matrix $M=WH$ and the function $f$ that we will need in Section 2.2 to prove that Algorithm 1 is guaranteed to recover columns of $M$ corresponding to columns of the matrix $W$. Then, we analyze Algorithm 1 in case some noise is added to the input separable matrix $M$, and show that, under these assumptions, it is robust under any small perturbations; see Section 2.3. Finally, we compare our results with the ones from in Section 2.4.

In the remainder of the paper, we will assume that the original data matrix $M=WH$ is separable, that is, each column of $W$ appears as a column of $M$. Recall that this condition is implied by the pure-pixel assumption in hyperspectral imaging; see Introduction. We will also assume that the matrix $W$ is full rank. This is often implicitly assumed in practice otherwise the problem is in general ill-posed, because the matrix $H$ is then typically not uniquely determined; see, e.g., .

The assumption on matrix $H$ is made without loss of generality by

Permuting the columns of $M$ so that the first $r$ columns of $M$ correspond to the columns of $W$ (in the same order).

Normalizing $M$ so that the entries of each of its columns sum to one (except for its zero columns). In fact, we have that

By construction, the entries of each column of $MD_{M}^{-1}$ and $WD_{W}^{-1}$ sum to one (except for the zero columns of $M$), while the entries of each column of $(D_{W}HD_{M}^{-1})$ have to sum to one (except for ones corresponding to the zero columns of $M$ which are equal to zero) since $M=WH$.

In the hyperspectral imaging literature, the entries of each column of matrix $H$ are typically assumed to sum to one, hence Assumption 1 is slightly more general. This has several advantages:

It allows the image to contain ‘background’ pixels with zero spectral signatures, which are present for example in hyperspectral images of objects in outer space (such as satellites).

It allows us to take into account different intensities of light among the pixels in the image, e.g., if there are some shadow parts in the scene or if the angle between the camera and the scene varies. Hence, although some pixels contain the same material(s) with the same abundance(s), their spectral signature could differ by a scaling factor.

In the noisy case, it allows us to take into account endmembers with very small spectral signature as noise, although it is not clear whether relaxing the sum-to-one constraint is the best approach .

Our assumptions actually do not require the matrix $M$ to be nonnegative, as $W$ can be any full-rank matrix. In fact, after the first step of Algorithm 1, the residual matrix will typically contain negative entries.

We will also need to assume that the function $f$ in Algorithm 1 satisfies the following conditions.

Notice that, for any strongly convex function $g$ whose gradient is Lipschitz continuous and whose global minimizer is $\bar{x}$, one can construct the function $f(x)=g(\bar{x}+x)-g(\bar{x})$ satisfying Assumption 2. In fact, $f(0)=0$ while $f(x)\geq 0$ for any $x$ since $g(\bar{x}+x)\geq g(\bar{x})$ for any $x$. Recall that (see, e.g., ) a function is strongly convex with parameter $\mu$ if and only if it is convex and for any $x,y\in\text{dom}(f)$

for any $\delta\in$. Moreover, its gradient is Lipschitz continuous with constant $L$ if and only if for any $x,y\in\text{dom}(f)$

Convex analysis also tells us that if $f$ satisfies Assumption 2 then, for any $x,y$,

since $f(0)=0$ and $\nabla f(0)=0$ (because zero is the global minimizer of $f$).

The most obvious choice for $f$ satisfying Assumption 2 is $f(x)=||x||^{2}_{2}$; we return to this matter in Section 4.1.

2 Noiseless Case

We now prove that, under Assumption 1 and 2, Algorithm 1 recovers a set of indices corresponding to the columns of $W$.

where $e_{j}$ is the $j$th column of the identity matrix.

By assumption on $f$, we have $f(w)>0$ for any $w\neq 0$; see Equation (1). Hence, if $Yh=0$, we have the result since $f(Yh)=0<f(w_{i})$ for all $i$. Otherwise $Yh=\sum_{i=1}^{k}w_{i}h_{i}$ where $h_{i}\neq 0$ for at least one $1\leq i\leq k$ so that

The first inequality is strict since $h\neq e_{j}\forall j$ and $h_{i}\neq 0$ for at least one $1\leq i\leq k$, and the second follows from the fact that $\sum_{i=1}^{k}h_{i}\leq\sum_{i=1}^{r}h_{i}\leq 1$. ∎

Let the matrix $M=WH$ satisfy Assumption 1 and the function $f$ satisfy Assumption 2. Then Algorithm 1 recovers a set of indices $J$ such that $M(:,J)=W$ up to permutation.

First step. Lemma 1 applies since $f$ satisfies Assumption 2 while $W$ is full rank. Therefore, the first step of Algorithm 1 extracts one of the columns of $W$. Assume without loss of generality the last column $w_{r}$ of $W$ is extracted, then the first residual has the form

i.e., the matrix $R^{(1)}$ is obtained by projecting the columns of $M$ onto the orthogonal complement of $w_{r}$. We observe that $W^{(1)}$ satisfies the conditions of Lemma 1 as well because $W^{\prime}$ is full rank since $W$ is. This implies, by Lemma 1, that the second step of Algorithm 1 extracts one of the columns of $W^{\prime}$.

Induction step. Assume that after $k$ steps the residual has the form $R^{(k)}=[W^{*}\;\mathbf{0}_{m\times k}]H$ with $W^{*}$ full rank. Then, by Lemma 1, the next extracted index will correspond to one of the columns of $W^{*}$ (say, without loss of generality, the last one) and the next residual will have the form $R=[W^{\dagger},\mathbf{0}_{m\times(k+1)}]H$ where $W^{\dagger}$ full rank since $W^{*}$ is, and $H$ is unchanged. By induction, after $r$ steps, we have that the indices corresponding to the different columns of $W$ have been extracted and that the residual is equal to zero ($R=\mathbf{0}_{m\times r}H$). ∎

3 Adding Noise

In this section, we analyze how perturbing the input data matrix affects the performances of Algorithm 1. We are going to assume that the input perturbed matrix $M^{\prime}$ can be written as $M^{\prime}=M+N$ where $M$ is the noiseless original separable matrix satisfying Assumption 1, and $N$ is the noise with $||n_{i}||_{2}\leq\epsilon$ for all $i$ for some sufficiently small $\epsilon\geq 0$.

Given a matrix $W$, we introduce the following notations: $\gamma(W)=\min_{i\neq j}||w_{i}-w_{j}||_{2}$, $\nu(W)=\min_{i}||w_{i}||_{2}$, $\omega(W)=\min\left\{\nu(W),\frac{1}{\sqrt{2}}\gamma(W)\right\}$, and $K(W)=\max_{i}||w_{i}||_{2}$.

then, for any $0\leq\delta\leq\frac{1}{2}$,

satisfies $f^{*}\leq\max_{i}f(w_{i})-\frac{1}{2}\,\mu\,(1-\delta)\,\delta\,\omega(W)^{2}$.

$x^{*}=(1-\delta)e_{j}$ for $1\leq j\leq k$,

$x^{*}=\delta e_{i}+(1-\delta)e_{j}$ for $1\leq i,j\leq k$, or

$x^{*}=\delta e_{i}+(1-\delta)e_{j}$ for $k+1\leq i\leq r$ and $1\leq j\leq k$.

Before analyzing the different cases, let us provide a lower bound for $f^{*}$. Using Equation (1), we have

Since $(1-\delta)w_{i}$ is a feasible solution and $0\leq\delta\leq\frac{1}{2}$, this implies $f^{*}\geq\frac{\mu}{8}K(W)^{2}$. Recall that since $f$ is strongly convex with parameter $\mu$, we have

Clearly, $x^{*}\neq 0$ since $f(0)=0$ and $f(y)>0$ for all $y\neq 0$, cf. Equation (1).

$Yx^{*}=q_{i}$ for some $i$. Using Equation (1), we have

since $\nu(W)>2\sqrt{\frac{L}{\mu}}K(Q)$, a contradiction.

By strong convexity, we also have $f(w_{i})\geq\frac{\mu}{2}||w_{i}||_{2}^{2}\geq\frac{1}{2}\mu(1-\delta)||w_{i}||_{2}^{2}$. Plugging it in (3) gives

$Yx^{*}=\delta w_{i}+(1-\delta)w_{j}$ for some $i\neq j$ :

$Yx^{*}=\delta q_{i}+(1-\delta)w_{j}$ for some $i$, $j$. First, we have

In fact, $||q_{i}-w_{j}||_{2}\geq||w_{j}||_{2}-||q_{i}||_{2}\geq\nu(W)-K(Q)$. Then, using

$\nu(W)-K(Q)>\left(1-\frac{1}{2}\sqrt{\frac{\mu}{L}}\right)\nu(W)\geq\frac{1}{2}\nu(W)$ and $f(w_{j})\geq\frac{\mu}{2}(1-\delta)\nu(W)^{2}$, we obtain

For the upper bound (4), we use the fact that the gradient of $f$ is Lipschitz continuous with constant $L$

for any $||x||_{2}\leq K$, $||n||_{2}\leq\epsilon\leq K$. The second inequality follows from the fact that $\nabla f(0)=0$ and by Lipschitz continuity of the gradient: $||\nabla f(x)-0||_{2}\leq L||x-0||_{2}\leq LK$ for any $||x||_{2}\leq K$.

For the lower bound (5), we use strong convexity

for any $||x||_{2}\leq K$, $||n||_{2}\leq\epsilon\leq K$. The third inequality follows from the fact that

We can now prove the theorem which will be used in the induction step to prove that Algorithm 1 works under small perturbations of the input separable matrix.

$f$ satisfy Assumption 2, with strong convexity parameter $\mu$, and its gradient have Lipschitz constant $L$.

$\epsilon$ be sufficiently small so that $\epsilon\leq\frac{\mu{\omega}^{2}}{20KL}$.

Then the index $i$ corresponding to a column $m^{\prime}_{i}$ of $M^{\prime}$ that maximizes the function $f$ satisfies

and $\delta=\frac{10\epsilon KL}{\mu\omega^{2}}$, which implies

First note that $\epsilon\leq\frac{\mu\omega^{2}}{20KL}$ implies $\delta=\frac{10\epsilon KL}{\mu\omega^{2}}\leq\frac{1}{2}$. Let us then prove Equation (6) by contradiction. Assume the extracted index, say $i$, for which $m^{\prime}_{i}=m_{i}+n_{i}=Yh_{i}+n_{i}$ satisfies $h_{i}(l)<1-\delta$ for $1\leq l\leq k$. We have

where $w^{\prime}_{j}$ is the perturbed column of $M$ corresponding to $w_{j}$ (that is, the $j$th column of $M^{\prime}$). The first inequality follows from Lemma 3. In fact, we have $\epsilon\leq K$ since $\mu\leq L$ and $\omega\leq K$, $||m_{i}||_{2}=||Wh_{i}||_{2}\leq\max_{i}||w_{i}||_{2}=K$ (by convexity of $||.||_{2}$), and $||n_{i}||_{2}\leq\epsilon$ $\forall i$ so that $f(m^{\prime}_{i})\leq f(m_{i})+\frac{3}{2}\epsilon KL$. The second inequality is strict since the maximum is attained at a vertex with $x(l)=1-\delta$ for some $1\leq l\leq k$ at optimality (see proof of Lemma 2). The third inequality follows from Lemma 2 while the fourth follows from the fact that $||w_{j}||_{2}\leq K$ so that $f(w_{j})-\epsilon KL\leq f(w^{\prime}_{j})$ for all $j$ by Lemma 3.

We notice that, since $\delta\leq\frac{1}{2}$,

Combining this inequality with Equation (8), we obtain $f(m^{\prime}_{i})<\max_{j}f(w^{\prime}_{j})$, a contradiction since $m^{\prime}_{i}$ should maximize $f$ among the columns of $M^{\prime}$ and the $w^{\prime}_{j}$’s are among the columns of $M^{\prime}$.

To prove Equation (7), we use Equation (6) and observe that

so that $\sum_{k\neq p}\alpha_{k}\leq\delta^{\prime}\leq\delta$. Therefore,

It is interesting to relate the ratio $\frac{K(W)}{\omega(W)}$ to the condition number of matrix $W$, given by the ratio of its largest and smallest singular values $\kappa(W)=\frac{\sigma_{1}(W)}{\sigma_{r}(W)}$.

In particular, this inequality implies that if $\kappa(W)=1$ then $\frac{K(W)}{\omega(W)}=1$.

3.2 Error Bound for Algorithm 1

We have shown that, if the input matrix $M^{\prime}$ has the form

where $Q$ and $N$ are sufficiently small and the sum of the entries of each column of $H^{\prime}\geq 0$ is smaller than one, then Algorithm 1 extracts one column of $M^{\prime}$ which is close to a column of $W$; cf. Theorem 2. We now show that, at each step of Algorithm 1, the residual matrix satisfies these assumptions so that we can prove the result by induction.

We first give some useful lemmas; see and the references therein.

We can now prove the main theorem of the paper which shows that, given a noisy separable matrix $M^{\prime}=M+N=WH+N$ where $M$ satisfies Assumption 1, Algorithm 1 is able to identify approximately the columns of $W$.

and $J$ be the index set of cardinality $r$ extracted by Algorithm 1. Then there exists a permutation $P$ of $\{1,2,\dots,r\}$ such that

Let us prove the result by induction. First, let us define the residual matrix $R^{(k)}$ obtained after $k$ steps of Algorithm 1 as follows:

with $P^{(k)}=\left(I-\frac{uu^{T}}{||u||_{2}^{2}}\right)$ is the orthogonal projection performed at step 5 of Algorithm 1 where $u$ is the extracted column of $R^{(k)}$, that is, $u=r^{(k)}_{i}$ for some $1\leq i\leq n$.

for some $1\leq p\leq r-k$. Let us assume without loss of generality that $p=r-k$. The next residual $R^{(k+1)}$ has the form

$||P^{(k)}w^{(k)}_{p}||_{2}\leq\bar{\epsilon}=\epsilon\left(1+80\frac{K(W)^{2}}{\sigma_{r}^{2}(W)}\frac{L}{\mu}\right)$ since

where the first inequality follows from $P^{(k)}w^{(k)}_{p}$ being the projection of $w^{(k)}_{p}$ onto the orthogonal complement of $r^{(k)}_{i}$. Moreover,

In fact, $K(W^{(k)})\leq K(W)$ because of the orthogonal projections, while $\omega(W^{(k)})\geq\frac{1}{2}\sigma_{r}(W)$ follows from

Lemma 7 applies since $K(Q^{(k)})\leq\bar{\epsilon}$. The last inequality follows from $\bar{\epsilon}\leq\frac{\sigma_{r}(W)}{2\sqrt{r-1}}$ since

For $R^{(k+1)}$ to satisfy the same conditions as $R^{(k)}$, it remains to show that $\nu(W^{(k+1)})>2\sqrt{\frac{L}{\mu}}\bar{\epsilon}$ and $\epsilon\leq\frac{\omega(W^{(k+1)})^{2}\mu}{20K(W^{(k+1)})L}$. Let us show that these hold for all $k=0,1,\dots,r-1$ :

Since $\nu(W^{(k)})\geq\omega(W^{(k)})\geq\frac{1}{2}\sigma_{r}(W)$ (see above), $\nu(W^{(k)})>2\sqrt{\frac{L}{\mu}}\bar{\epsilon}$ is implied by $\frac{1}{2}\sigma_{r}(W)>2\sqrt{\frac{L}{\mu}}\bar{\epsilon}$, that is,

By assumption on the matrix $M^{\prime}$, these conditions are satisfied at the first step of the algorithm (we actually have that $Q^{(0)}$ is an empty matrix), so that, by induction, all the residual matrices satisfy these conditions. Finally, Theorem 2 implies that the index $i$ extracted by Algorithm 1 at the $(k+1)$th step satisfies

$\delta^{(k)}=\frac{10\epsilon K(W^{(k)})L}{\omega(W^{(k)})^{2}\mu}$, and $p=r-k$ without loss of generality. Since the matrix $H$ is unchanged between each step, this implies $m_{i}=Wh_{i}$ where $h_{i}(p)\geq 1-\delta^{(k)}$, hence

The second inequality is obtained using $m_{i}=Wh_{i}=h_{i}(p)w_{p}+\sum_{k\neq p}h_{i}(k)w_{k}$ and $\sum_{k\neq p}h_{i}(k)\leq 1-h_{i}(p)$ so that

4 Bounds for Separable NMF and Comparison with the Algorithms of Arora et al. [3] and Bittorf et al. [6]

Arora et al. identify a matrix $U$ such that

The algorithm of Bittorf et al. identifies a matrix $U$ satisfying

By Theorem 3, Algorithm 1 therefore requires

This shows that the above bounds are tighter, as they only require the noise to be bounded above by a constant proportional to $\alpha^{2}$ to guarantee an NMF with error proportional to $\frac{\epsilon_{1}}{\alpha}$. In particular, if $W$ is not full rank, Algorithm 1 will fail to extract more than $\operatorname{rank}(W)$ columns of $W$, while the value of $\alpha$ can be much larger than zero implying that the algorithms from will still be robust to a relatively large noise.

To conclude, the techniques in based on linear programming lead to better error bounds. However, there are computationally much more expensive (at least quadratic in $n$, while Algorithm 1 is linear in $n$, cf. Section 1.1), and have the drawback that some parameters have to be estimated in advance: the noise level $\epsilon_{1}$, and

the parameter $\alpha$ for Arora et al. (which is rather difficult to estimate as $W$ is unknown),

the factorization rank $r$ for Bittorf et al. In their incremental gradient descent algorithm, the parameter $\epsilon$ does not need to be estimated. However, other parameters need to be tuned, namely, primal and dual step sizes., hence the solution has to be recomputed from scratch when the value of $r$ is changed (which often happens in practice as the number of columns to be extracted is typically estimated with a trial and error approach).

Moreover, the algorithms from heavily rely on the separability assumption while Algorithm 1 still makes sense even if the separability assumption is not satisfied; see Section 4.1. Table 1 summarizes these results. Note that we keep the analysis simple and only indicate the growth in terms of $n$. The reason is threefold: (1) in many applications (such as hyperspectral unmixing), $n$ is much larger than $m$ and $r$, (2) a more detailed comparison of the running times would be possible (that is, in terms of $m$, $n$, and $r$) but is not straightforward as it depends on the algorithm used to solve the linear programs (and possibly on the parameters $\alpha$ and $\epsilon_{1}$), and (3) both algorithms are at least quadratic in $n$ (for example, the computational cost of each iteration of the first-order method proposed in is proportional to $mn^{2}$, so that the complexity is linear in $m$).

Outlier Detection

It is important to point out that Algorithm 1 is very sensitive to outliers, as are most algorithms aiming to detect the vertices of the convex hull of a set of points, e.g., the algorithms from discussed in the previous section. Therefore, one should ideally discard the outliers beforehand, or design variants of these methods robust to outliers. In this section, we briefly describe a simple way for dealing with (a few) outliers. This idea is inspired from the approach described in .

where $h^{\prime}_{i}\in\Delta^{r}$ for all $i$, hence $f^{\prime}_{i}\in\Delta^{r}$ for all $i$. Assuming $[W,T]$ has rank $r+t$ (hence $t\leq m-r$), the matrix $M$ above also satisfies Assumption 1. In the noiseless case, Algorithm 1 will then extract a set of indices corresponding to columns of $W$ and $T$ (Theorem 1). Therefore, assuming that the matrix $H^{\prime}$ has at least one non-zero element in each row, one way to identifying the outliers would be to

Extract $r+t$ columns from the matrix $M$ using Algorithm 1,

Compute the corresponding optimal abundance matrix $F$, and

Select the $r$ columns corresponding to rows of $F$ with the largest sum,

see Algorithm 2. (Note that Algorithm 2 requires to solve a convex quadratic program, hence it is computationally much more expensive than Algorithm 1.)

It is easy to check that Algorithm 2 will recover the $r$ columns of $W$ because the optimal solution $G$ computed at the second step is unique and equal to $F$ (since $[W,T]$ is full rank), hence $||G(j,$:$)||_{1}>1$ for the indices corresponding to the columns of $W$ while $||G(j,$:$)||_{1}=1$ for the outliers; see Equation (14).

In the noisy case, a stronger condition is necessary: the sum of the entries of each row of $H^{\prime}$ must be larger than some bound depending on the noise level. In terms of hyperspectral imaging, it means that for an endmember to be distinguishable from an outlier, its abundance in the image should be sufficiently large, which is perfectly reasonable.

and $J$ be the index set of cardinality $r$ extracted by Algorithm 2. If

then there exists a permutation $P$ of $\{1,2,\dots,r\}$ such that

By Theorem 3, the columns extracted at the first step of Algorithm 2 correspond to the columns of $W$ and $T$ up to error $\bar{\epsilon}$. Let then $W+N_{W}$ and $T+N_{T}$ be the columns extracted by Algorithm 1 with $K(N_{W}),K(N_{T})\leq\bar{\epsilon}$.

At the second step of Algorithm 2, the matrix $G$ is equal to

up to the permutation of its rows. It remains to show that

so that the last step of Algorithm 2 will identify correctly the columns of $W$ among the ones extracted at the first step. We are going to show that

More precisely, we are going to prove the following lower (resp. upper) bounds for the entries of the first $r$ (resp. last $t$) rows of $G$:

For $1\leq i\leq r$, $G_{ij}\geq\max\left(0,F_{ij}-2\frac{\bar{\epsilon}+\epsilon}{\sigma_{s}(B)}\right)$ for all $j$.

For $r+1\leq i\leq r+t$, $G_{ij}\leq\min\left(1,F_{ij}+2\frac{\bar{\epsilon}+\epsilon}{\sigma_{s}(B)}\right)$ for all $j$.

Therefore, assuming (a) and (b) hold, we obtain

which proves the result. It remains to prove (a) and (b). First, we have that

since $G(:,j)$ leads to the best approximation of $m^{\prime}_{j}$ over $\Delta$ (see step 2 of Algorithm 2) and $f_{j}\in\Delta$.

Then, let us prove the upper bound for the block of matrix $G$ at position $(2,1)$, that is, let us prove that

(Note that $0\leq G\leq 1$ by construction, hence some of the bounds are trivial, e.g., for the block (1,2).) The derivations necessary to obtain the bounds for the other (non-trivial) blocks are exactly the same and are then omitted here. Let $r+1\leq i\leq t$ and $1\leq j\leq r$ and denote $G_{ij}=\delta$, and let also $I=\{1,2,\dots,r+t\}\backslash\{i\}$. We have

The first inequality follows from $K(N)\leq\epsilon$, $K([N_{W},N_{T}])\leq\bar{\epsilon}$ and $G_{:j}\in\Delta^{r+t}$, while the second inequality follows from the fact that $w_{j}$ is a column of $B(:,I)$. The last inequality follows from the fact that the projection of any column of $B$ onto the subspace spanned by the other columns is at least $\sigma_{s}(B)$. Finally, using Equations (17) and (18), we have $\bar{\epsilon}+\epsilon\geq\delta\sigma_{s}(B)-\bar{\epsilon}-\epsilon$, hence $G_{ij}=\delta\leq 2\frac{\bar{\epsilon}+\epsilon}{\sigma_{s}(B)}$. ∎

Choices for f𝑓f and Related Methods

In this section, we discuss several choices for the function $f$ in Algorithm 1, and relate them to existing methods.

According to our derivations (see Theorem 3), using functions $f$ whose strong convexity parameter $\mu$ is equal to the Lipschitz constant $L$ of its gradient is the best possible choice (since it minimizes the error bounds). The only function satisfying Assumption 2 along with this condition is, up to a scaling factor, $f(x)=||x||_{2}^{2}=\sum_{i=1}^{m}x_{i}^{2}$. In fact, Assumption 2 implies $\frac{\mu}{2}||x||_{2}^{2}\leq f(x)\leq\frac{L}{2}||x||_{2}^{2}$; see Equation (1). However, depending on the problem at hand, other choices could be more judicious (see Sections 4.2 and 4.3). It is worth noting that Algorithm 1 with $f(x)=||x||_{2}^{2}$ has been introduced and analyzed by several other authors:

Successive Projection Algorithm. Araújo et al. proposed the successive projection algorithm (SPA), which is equivalent to Algorithm 1 with $f(x)=||x||_{2}^{2}$. They used it for variable selection in spectroscopic multicomponent analysis, and showed it works better than other standard techniques. In particular, they mention ‘SPA seems to be more robust than genetic algorithms’ but were not able to provide a rigorous justification for that fact (which our analysis does). Ren and Chang rediscovered the same algorithm, which was referred to as the automatic target generation process (ATGP). It was empirically observed in to perform better than other hyperspectral unmixing techniques (namely, PPI and VCA ). However, no rigorous explanation of their observations was provided. In Section 5, we describe these techniques and explain why they are not robust to noise, which theoretically justifies the better performances of Algorithm 1. Chan et al. analyzed the same algorithm (with the difference that the data is preprocessed using a linear dimensionality reduction technique). The algorithm is referred to as the successive volume maximization algorithm (SVMAX). They also successfully use Algorithm 1 as an initialization for a more sophisticated approach which does not take into account the pure-pixel assumption.

Greedy Heuristic for Volume Maximization. Çivril and Magdon-Ismail showed that Algorithm 1 with $f(x)=||x||_{2}^{2}$ is a very good greedy heuristic for the following problem: given a matrix $M$ and an integer $r$, find a subset of $r$ columns of $M$ whose convex hull has maximum volume. More precisely, unless $\mathcal{P}=\mathcal{NP}$, they proved that the approximation ratio guaranteed by the greedy heuristic is within a logarithmic factor of the best possible achievable ratio by any polynomial-time algorithm. However, the special case of separable matrices was not considered. This is another advantage of Algorithm 1: even if the input data matrix $M$ is not approximately separable, it identifies $r$ columns of $M$ whose convex hull has large volume. For the robust algorithms from discussed in Section 2.4, it is not clear whether they will be able to produce a meaningful output in that case; see also Section 5.2 for some numerical experiments.

This choice limits the impact of large entries in $x$, hence would potentially be more robust to outliers. In particular, as $\alpha$ goes to zero, $f(x)$ converges to $||x||_{1}$ while, when $\alpha$ goes to infinity, it converges to $\frac{||x||_{2}^{2}}{\alpha}$ (in any bounded set).

is strongly convex with parameter $\mu=\frac{2\alpha^{2}}{(\alpha+K)^{3}}$ and its gradient is Lipschitz continuous with constant $L=\frac{2}{\alpha}$.

hence $\mu=\frac{2\alpha^{2}}{(\alpha+K)^{3}}$, and $L=\frac{2}{\alpha}$. ∎

For example, one can choose $\alpha=K$ for which we have $\frac{L}{\mu}=2$, which is slightly larger than one but is less sensitive to large, potentially outlying, entries of $M$. Let us illustrate this on a simple example:

One can check that, for any $\epsilon\leq 0.69$, Algorithm 1 with $f(x)=||x||_{2}^{2}$ recovers the first two columns of $M$, that is, the columns of $W$. However, using $f(x)=\sum_{i}\frac{x_{i}^{2}}{1+|x_{i}|}$, Algorithm 1 recovers the columns of $W$ for any $\epsilon\leq 1.15$. Choosing appropriate function $f(x)$ depending on the input data matrix and the noise model is a topic for further research.

The condition that the gradient of $f$ must be Lipschitz continuous in Assumption 2 can be relaxed to the condition that the gradient of $f$ is continuously differentiable. In fact, in all our derivations, we have always assumed that $f$ was applied on a bounded set (more precisely, the ball $\{x\ |\ ||x||_{2}\leq K(W)\}$). Since $g\in\mathcal{C}^{1}$ implies that $g$ is locally Lipschitz continuous, the condition $f\in\mathcal{C}^{2}$ is sufficient for our analysis to hold. Similarly, the strong convexity condition can be relaxed to local strong convexity.

For $1<p\leq 2$, $f(x)$ is strongly convex with parameter $2(p-1)$ with respect to the norm $||.||_{p}$ [22, Section 4.1.1], while its gradient is locally Lipschitz continuous (see Remark 3). For $2\leq p<+\infty$, the gradient of $f(x)$ is Lipschitz continuous with respect to the norm $||.||_{p}$ with constant $2(p-1)$ (by duality), while it is locally strongly convex. Therefore, $f$ satisfies Assumption 2 for any $1<p<+\infty$ in any bounded set, hence our analysis applies. Note that, for $p=1$ and $p=+\infty$, the algorithm is not guaranteed to work, even in the noiseless case (when points are on the boundary of the convex hull of the columns of $W$): consider for example the following separable matrices

Numerical Experiments

In the first part of this section, we compare Algorithm 1 with several fast hyperspectral unmixing algorithms under the linear mixing model and the pure-pixel assumption. We first briefly describe them (computational cost and main properties) and then perform a series of experiments on synthetic data sets in order to highlight their properties. For comparisons of Algorithm 1 with other algorithms on other synthetic and real-world hyperspectral data sets, we refer the reader to since Algorithm 1 is a generalization of the algorithms proposed in ; see Section 4.

In the second part of the section, we compare Algorithm 1 with the Algorithm of Bittorf et al. .

Algorithm 1 with $f(x)=||x||_{2}^{2}$. We will only test this variant because, according to our analysis, it is the most robust. (Comparing different variants of Algorithm 1 is a topic for further research.) The computational cost is rather low: steps 3 and 5 are the only steps requiring computation, and have to be performed $r$ times. We have

Step 3. Compute the squared norm of the columns of $R$, which requires $n$ times $2m$ operations (squaring and summing the elements of each column), and extract the maximum, which requires $n$ comparisons, for a total of approximately $2mn$ operations.

Step 5. It can be compute in the following way

where computing $x^{T}=u_{j}^{T}R$ requires $2mn$ operations, $y=\frac{u_{j}}{||u_{j}||_{2}^{2}}$ $m$ operations, and $R-yx^{T}$ $2mn$ operations, for a total of approximately $4mn$ operations.

The total computational cost of Algorithm 1 is then about $6mnr$ operations, plus some negligible terms.

If the matrix $M$ is sparse, $R$ will eventually become dense which is often impractical. Therefore, $R$ should be kept in memory as the original matrix $M$ minus the rank-one updates.

It is not robust to noise. In fact, linear functions can be maximized at any vertex of the convex hull of a set of points. Therefore, in the noisy case, as soon as a column of the perturbed matrix $M^{\prime}$ is not contained in the convex hull of the columns of $W$, it can be identified as a vertex. This can occur for arbitrarily small perturbation, as will be confirmed by the experiments below.

If not enough linear functions are generated, the algorithm might not be able to identify all the vertices (even in the noiseless case). This is particularly critical in case of ill-conditioning because the probability that some vertices maximize a randomly generated linear function can be arbitrarily low.

If the input noisy data matrix contains many columns close to a given column of matrix $W$, the score of these columns will be typically small (they essentially share the score of the original column of matrix $W$), while an isolated column which does not correspond to a column of $W$ could potentially have a higher score than these columns, hence be extracted. This can for example be rather critical for hyperspectral images where there typically are many pixels close to pure pixels (i.e., columns of $M$ corresponding to the same column of $W$). Moreover, for the same reasons, PPI might extract columns of $M$ corresponding to the same column of $W$.

Vertex Component Analysis (VCA) . The first step of VCA is to preprocess the data using principal component analysis which requires $O(nm^{2}+m^{3})$ . Then, the core of the algorithm requires $O(rm^{2})$ operations (see below), for a total of $O(nm^{2}+m^{3})$ operationsThe code is available at http://www.lx.it.pt/~bioucas/code.htm.. Notice that the preprocessing is particularly well suited for data sets where $m\ll n$, such as hyperspectral images, where $m$ is the number of hyperspectral images with $m\sim 100$, while $n$ is the number of pixels per image with $n\sim 10^{6}$. The core of VCA is very similar to Algorithm 1: at each step, it projects the data onto the orthogonal complement of the extracted column. However, instead of using a strictly convex function to identify a vertex of the convex hull of $M$ (as in Algorithm 1), it uses a randomly generated linear function (that is, it selects the column maximizing the function $f(x)=c^{T}x$ where $c$ is randomly generated, as PPI does). Therefore, for the same reasons as for PPI, the algorithm is not robust. However, it solves the second and third pitfalls of PPI (see point (b) and (c) above), that is, it will always be able to identify enough vertices, and a cluster of points around a vertex are more likely to be extracted than an isolated point. Note that, in the VCA implementation, only one linear function is generated at each step which makes it rather sensitive to this choice. Finally, the two main differences between VCA and Algorithm 1 are that: (1) VCA uses a pre-processing (although we could implement a version of Algorithm 1 with the same pre-processing), and (2) VCA uses randomly generated linear functions to pick vertices of the convex hull of the columns of $M$ (which makes it non-robust to noise, and also non-deterministic).

Simplex Volume Maximization (SiVM) . SiVM recursively extracts columns of the matrix $M$ while trying to maximize the volume of the convex hull of the corresponding columns. Because evaluating the volumes induced by adding a column not yet extracted to the previously selected ones is computationally expensive, the function is approximated via a heuristic, for a total computational cost of $O(mnr)$ operations (as for Algorithm 1). The heuristic assumes that the columns of $W$ are located at the same distance, that is, $||w_{i}-w_{j}||_{2}=||w_{k}-w_{l}||_{2}$ for all $i\neq j$ and $k\neq l$. Therefore, there is not guarantee that the algorithm will work, even in the noiseless case; this will be particularly critical for ill-conditioned problems, which will be confirmed by the experiments below.

We now generate several synthetic data sets allowing to highlight the properties of the different algorithms, and in particular of Algorithm 1. We are going to consider noisy separable matrices generated as follows.

The matrix $W$ will be generated in two different ways :

The matrices $H$ and $N$ will be generated in two different ways as well :

where $\bar{w}$ is the average of the columns of $W$ (geometrically, this is the vertex centroid of the convex hull of the columns of $W$). This means that we move the columns of $M$ toward the outside of the convex hull of the columns of $W$. Hence, for any $\delta>0$, the columns of $M^{\prime}$ are not contained in the convex hull of the columns of $W$ (although the rank of $M^{\prime}$ remains equal to 20).

Finally, we construct the noisy separable matrix $M^{\prime}=WH+N$ in four different ways, see Table 2, where $W$, $H$ and $N$ are generated as described above for a total of four experiments.

For each experiment, we generate 100 matrices for 100 different values of $\delta$ and compute the percentage of columns of $W$ that the algorithms were able to identify (hence the higher the curve, the better); see Figure 1. (Note that we then have 10000 matrices generated for each experiment.)

Algorithm 1 is the most robust algorithm as it is able to identify all the columns of $W$ for the largest values of the perturbation $\delta$ for all experiments; see Table 3 and Figure 1.

In Exp. 1, PPI and SiVM perform relatively well, the reason being that the matrix $W$ is well-conditioned (see Table 2) while VCA is not robust to any noise. In fact, as explained in Section 5.1, VCA only uses one randomly generated linear function to identify a column of $W$ at each step, hence can potentially extract any column of $M$ since they all are vertices of the convex hull of the columns of $M$ (the last columns of $M$ are the middle points of the columns of $W$ and are perturbed toward the outside of the convex hull of the columns of $W$).

In Exp. 2, the matrix $W$ is well-conditioned so that SiVM still performs well. PPI is now unable to identify all columns of $M$, because of the repetition in the data set (each column of $W$ is present twice as a column of $M$)For $\delta=0$, we observe that our implementation of the PPI algorithm actually recovers all columns of $W$. The reason is that the first columns of $M$ are exactly equal to each other and that the MATLAB max(.) function only outputs the smallest index corresponding to a maximum value. This is why PPI works in the noiseless case even when there are duplicates.. VCA now performs much better because the columns of $M$ are strictly contained in the interior of the convex hull of the columns of $W$.

In Exp. 3 and 4, SiVM performs very poorly because of the ill-conditioning of matrix $W$.

In Exp. 3, as opposed to Exp. 1, PPI is no longer robust because of ill-conditioning, although more than 97% of the columns of $W$ are perfectly extracted for all $\delta\leq 10^{-3}$. VCA is not robust but performs better than PPI, and extracts more than 97% of the columns of $W$ for all $\delta\leq 10^{-2}$ (note that Algorithm 1 does for $\delta\leq 0.05$).

In Exp. 4, PPI is not able to identify all the columns of $W$ because of the repetition, while, as opposed to Exp. 2, VCA is not robust to any noise because of ill-conditioning.

Algorithm 1 is the fastest algorithm although PPI and SiVM have roughly the same computational time. VCA is slower as it uses PCA as a preprocessing; see Table 4.

These experiments also show that the error bound derived in Theorem 3 is rather loose, which can be partly explained by the fact that our analysis considers the worst-case scenario (while our experiments use either a structured noise or Gaussian noise). Recall that the value of $\epsilon$ in Theorem 3 is the smallest value such that $||n_{i}||_{2}\leq\epsilon$ for all $i$; see Equation (9). Table 2 gives the average value of the maximum norm of the columns of $N$ for each experiment. Based on these values, the first row of Table 5 shows the average upper bound for $\delta$ to guarantee recovery; see Theorem 3.

2 Comparison with the Algorithm of Bittorf et al. [6]

In this section, we compare the Algorithm of Bittorf et al. (BRRT) (Algorithm 3 in ; see also Algorithm 2 in )We do not perform a comparison with the algorithm of Arora et al. as it is not very practical (the value of $\alpha$ has to be estimated, see Section 2.4) and has already been shown to perform similarly as BRRT in . with Algorithm 1. BRRT has to solve a linear program with $\mathcal{O}(n^{2})$ variables which we solve using CVX . Therefore we are only able to solve small-scale problems (in fact, CVX uses an interior-point method): we perform exactly the same experiments as in the previous section but for $m=10$ and $r=5$ for all experiments, so that

$n=5+\binom{5}{2}=15$ for the first and third experiments (the last ten columns of $M$ are the middle points of the five columns of $W$).

$n=5+5+10=20$ for the second and fourth experiments (we repeat twice each endmember, and add 10 points in the convex hull of the columns of $W$).

The average running time for BRRT on these data sets using CVX is about two seconds while, for Algorithm 1, it is less than $10^{-3}$ seconds. (Bittorf et al. propose a more efficient solver than CVX for their LP instances. As mentioned in Section 2.4, even with their more efficient solver, Algorithm 1 is much faster for large $n$.) Figure 2 shows the percentage of correctly extracted columns with respect to $\delta$, while Table 6 shows the robustness of both algorithms.

Quite surprisingly, Algorithm 1 performs in average better than BRRT. Although BRRT is more robust in two of the four experiments (that is, it extracts correctly all columns of $W$ for a larger value of $\delta$), the percentage of columns it is able to correctly extract decreases much faster as the noise level increases. For example, Table 7 shows the maximum value of $\delta$ for which 99% percent of the columns of $W$ are correctly extracted. In that case, Algorithm 1 always performs better.

A possible explanation for this behavior is that, when the noise is too large, the condition for recovery are not satisfied as the input matrix is far from being separable. However, using Algorithm 1 still makes sense as it extracts columns whose convex hull has large volume while it is not clear what BRRT does in that situation (as it heavily relies on the separability assumption). Therefore, although BRRT guarantees perfect recovery for higher noise levels, it appears that, in practice, when the noise level is high, Algorithm 1 is preferable.

Conclusion and Further Work

In this paper, we have introduced and analyzed a new family of fast and robust recursive algorithms for separable NMF problems which are equivalent to hyperspectral unmixing problems under the linear mixing model and the pure-pixel assumption. This family generalizes several existing hyperspectral unmixing algorithms, and our analysis provides a theoretical framework to explain the better performances of these approaches. In particular, our analysis explains why algorithms like PPI and VCA are less robust against noise compared to Algorithm 1.

Many questions remain open, and would be interesting directions for further research:

Is it possible to provide better error bounds for Algorithm 1 than the ones of Theorem 3? In other words, is our analysis tight? Also, can we improve the bounds if we assume specific generative and/or noise models?

How can we choose appropriate functions $f(x)$ for Algorithm 1 depending on the input data matrix?

Can we design other robust and fast algorithms for the separable NMF problem leading to better error bounds?

Acknowledgments

The authors would like to thank the reviewers for their feedback which helped improve the paper significantly.

References