Successive Projection Algorithm Robust to Outliers

TL;DR

This paper introduces Robust SPA, an improved version of the successive projection algorithm, which is more resilient to outliers and considers data fitting, enhancing nonnegative matrix factorization applications like hyperspectral unmixing.

Contribution

The paper proposes Robust SPA (RSPA), a novel variant of SPA that handles outliers effectively and incorporates data fitting into index selection, with proven robustness in low-noise conditions.

Findings

RSPA outperforms SPA in the presence of outliers.

RSPA maintains robustness in low-noise environments.

Effective on synthetic and hyperspectral data.

Abstract

The successive projection algorithm (SPA) is a fast algorithm to tackle separable nonnegative matrix factorization (NMF). Given a nonnegative data matrix $X$, SPA identifies an index set $\mathcal{K}$ such that there exists a nonnegative matrix $H$ with $X \approx X(:,\mathcal{K})H$. SPA has been successfully used as a pure-pixel search algorithm in hyperspectral unmixing and for anchor word selection in document classification. Moreover, SPA is provably robust in low-noise settings. The main drawbacks of SPA are that it is not robust to outliers and does not take the data fitting term into account when selecting the indices in $\mathcal{K}$. In this paper, we propose a new SPA variant, dubbed Robust SPA (RSPA), that is robust to outliers while still being provably robust in low-noise settings, and that takes into account the reconstruction error for selecting the indices in $\mathcal{K}$. We illustrate the effectiveness of RSPA on synthetic data sets and hyperspectral images.

Tables

Table 1. Table 1 : Relative error ( 2 ) and, in brackets, computational time in seconds of SPA, RSPA(10,1,4) and RSPA(20,1,4) applied on three hyperspectral images.

	Urban	San Diego	Cuprite
SPA	9.58 (1.4)	12.62   (2.9)	1.83 (1.9)
RSPA(10,1,4)	7.65 (31)	6.63   (64)	1.78  (52)
RSPA(20,1,4)	6.66 (59)	6.03 (124)	1.83  (83)

Equations

$$f_i(x)=||P_i{P}_{i-1}\dots P_{1}x|{|}_2^2\text{ for }1\le i\le d,$$

$$f_{i+1}(R(:,k^{\prime }(i)))=\beta f_{i+1}(R(:,k(i)))\text{ with }\beta >1.$$

$${\alpha }^{\ast }=1-\sqrt{1-\frac{\beta ||x|{|}_2^2-||y|{|}_2^2}{\beta (u^{T}x{)}^2-(u^{T}y{)}^2}}.$$

$$||(I-\alpha u{u}^T)y|{|}_2^2=\beta ||(I-\alpha u{u}^T)x|{|}_2^2.$$

$$\underset{H\ge 0}{\min }\frac{||X-X(:,\mathcal{K})H|{|}_F}{||X|{|}_F}$$

Full text

Keywords: nonnegative matrix factorization, hyperspectral unmixing, pure-pixel search, successive projection algorithm, outliers.

1 Introduction

Given a nonnegtive data matrix $X\in\mathbb{R}^{m\times n}_{+}$ and a factorization rank $r$, NMF looks for nonnegative matrices $W\in\mathbb{R}^{m\times r}_{+}$ and $H\in\mathbb{R}^{r\times n}_{+}$ such that $WH\approx X$. NMF can be used for example in image analysis, document classification and hyperspectral unmixing; see, e.g., [4, 5] and the references therein. However, NMF is NP-hard in general [12]. Separable NMF is an NMF variant where it is assumed that $W=X(:,\mathcal{K})$ for some index set $\mathcal{K}$ of size $r$. This means that the basis matrix $W$ is contained within the data set. This assumption makes sense in several applications including document classification [2] and hyperspectral unmixing [10]. In hyperspectral unmixing, the separability assumption is known as the pure-pixel assumption and requires that for each material present in the image (also called endmember), at least one pixel contains only that material. The pure-pixel assumption has been used for a long time in the literature [10], but it is only rather recently that separable NMF algorithms with provable guarantees have been proposed [3]. Among these algorithms, the successive projection algorithm (SPA) is one of the most popular ones: it is very fast, simple to implement and robust to noise [9]. It was first introduced in [1] and has been rediscovered many times; see the discussion in [7]. However, SPA has two main drawbacks: (1) SPA is very sensitive to outliers, and (2) SPA does not take directly the data fitting term into account to select the indices in $\mathcal{K}$. SPA can be made robust to outliers either by properly preprocessing the data set and removing the outliers, or using a proper post-processing of the index set $\mathcal{K}$ [9]. However, these approaches do not alleviate the second drawback of SPA. Moreover, it would be useful to have an SPA variant robust to outliers in case these pre- and/or post-processings fail to identify all outliers.

In this paper, we propose a new variant of SPA that is robust to outliers and takes directly the data fitting term into account to select the indices in $\mathcal{K}$. Moreover, this variant retains the good properties of SPA: it is fast, simple to implement and robust to noise. The paper is organized as follows. In Section 2, we recall how SPA works and its properties. In Section 3, we present our new SPA variant, dubbed robust SPA (RSPA), that is robust to outliers and takes the data fitting term into account in the selection step. In Section 4, we illustrate the effectiveness of this new approach on synthetic data sets and hyperspectral images.

2 The successive projection algorithm

Algorithm 1 gives the pseudocode of SPA. SPA sequentially identifies indices in $\mathcal{K}$ using two steps: at iteration $k$, given the current residual matrix $R$,

• •

Selection step: add the index $k$ to $\mathcal{K}$ that maximizes $f(R(:,k))$ where $f$ is a given function (see Assumption 2).

• •

Projection step: project the residual $R$ onto the orthogonal complement of $R(:,k)$ (step 5 of Algorithm 1).

Let us define the class of matrices for which SPA will provably identify a subset $\mathcal{K}$ such that there exists a nonnegative matrix $H$ with $X=X(:,\mathcal{K})H$.

Assumption 1**.**

The separable matrix $X\in\mathbb{R}^{m\times n}$ can be written as $X=WH=W[I_{r},H^{\prime}]\Pi$, where $W\in\mathbb{R}^{m\times r}$ has rank $r$, $H\in\mathbb{R}^{r\times n}_{+}$, $I_{r}$ is the identity matrix of size $r$, $\Pi$ is a permutation matrix, and the sum of the entries of each column of $H^{\prime}$ is at most one.

Let us make a few remarks:

• •

SPA is similar to vertex component analysis (VCA) [11]: the main difference is in the selection step where VCA uses a linear function, which is not robust to noise.

• •

If the sum-to-one-constraint on the columns of $H^{\prime}$ is not satisfied by the input matrix $X$, it can be obtained by normalizing each column of $X$ to have uni $\ell_{1}$ norm [9].

• •

In the absence of noise, separable NMF is equivalent to identifying the vertices of a set of points; see [10, 5].

In the presence of bounded noise, SPA will be able to identify $\mathcal{K}$ such that $X(:,\mathcal{K})\approx W$ (up to permutation); see [9] where error bounds are provided. For this result to hold, the function $f$ must satisfy the following assumption.

Assumption 2**.**

The function $f:\mathbb{R}^{m}\to\mathbb{R}_{+}$ is strongly convex with parameter $\mu_{f}>0$, its gradient is Lipschitz continuous with constant $L_{f}$, and its global minimizer is the all-zero vector with $f(0)=0$.

The standard variant of SPA uses $f(x)=||x||_{2}^{2}$, and is the most robust to noise according to the analysis in [9] since the error bound depends on the ratio $L_{f}/\mu_{f}$. This ratio is the conditioning of the function $f$ and denoted $\kappa_{f}=L_{f}/\mu_{f}\geq 1$. For $f(x)=||x||_{2}^{2}$, we have $\kappa_{f}=1$. However, SPA remains robust in low-noise settings as long as Assumption 2 is satisfied. For example, one can choose any quadratic function $f=\frac{1}{2}x^{T}Qx$ where $Q$ is positive definite, and we have $\kappa_{f}=\frac{\lambda_{\max}(Q)}{\lambda_{\min}(Q)}$. Moreover, since the analysis of SPA is sequential, the analysis still holds if one chooses different functions $f$ to select the index to put in $\mathcal{K}$ at each step of SPA, as long as they satisfy Assumption 2.

3 Robust SPA

The main contribution of this paper is to leverage the flexibility of SPA in choosing the function $f$ in order to make SPA robust to outliers by taking into account the residual error during the selection step. Algorithm 2, which we refer to as robust SPA (RSPA), is our proposed robust variant of SPA. It only differs from SPA in the selection step.

Let us explain the selection step of RSPA described in Algorithm 3.

As opposed to SPA that simply picks the column of the current residual $R$ that maximizes $f$ (step 3 of SPA), RSPA uses $d$ well-chosen quadratic functions $f_{i}(x)$ $(1\leq i\leq d)$. Each function $f_{i}$ will correspond to a candidate column of $R$ with index $k(i)=\operatorname{argmax}_{j}f_{i}(R(:,j))$. Among these candidates, RSPA will select the one with the smallest residual after projection onto its orthogonal complement. To measure the norm of the residual, we use the $\ell_{p}$ norm of the vector containing the $\ell_{2}$ norms of the columns of the residual, but many other measures could be used. We have observed that using $p=1$ works well, as it is less sensitive to large entries; see Section 4. As long as the functions $\{f_{i}\}_{i=1}^{d}$ satisfy Assumption 2, RSPA is guaranteed to be robust in low-noise settings. Moreover, this selection step will be more robust to outliers because outliers lead in general to a smaller decrease in the residual since they are less/not correlated with the data points.

It remains to explain how the functions $\{f_{i}\}_{i=1}^{d}$ are generated. Note that one could be tempted to generate them randomly but this will most likely still lead to the identification of outliers like in SPA; in particular if an outlier has a very large norm. For example, if one uses quadratic functions $x^{T}P^{T}Px=||Px||_{2}^{2}$ where the entries of $P$ are randomly generated, an outlier with a large norm will most likely also have a large value for $||Px||_{2}^{2}$. Hence we generate $\{f_{i}\}_{i=1}^{d}$ so that the associated data points maximizing them are well spread in the data set. To do so, we define

[TABLE]

where $P_{1}=I_{m}$ and $P_{i+1}=I_{m}-\alpha_{i}u_{i}u_{i}^{T}$ for some well chosen $\alpha_{i}\in(0,1)$ and $u_{i}$ with unit $\ell_{2}$ norm. Since $\alpha_{i}\in(0,1)$ and $||u_{i}||_{2}=1$, the matrices $\{P_{i}\}_{i=1}^{d}$ are positive definite (all eigenvalues are equal to one except for one which is equal to $1-\alpha_{i}$) hence $\{f_{i}\}_{i=1}^{d}$ are strongly convex functions. Note that $f_{1}(x)=||x||_{2}^{2}$ hence, for $d=1$, RSPA is equivalent to SPA with $f(x)=||x||_{2}^{2}$. The matrices $\{P_{i}\}_{i=1}^{d}$ are chosen such that there is a diversification of the data points maximizing the functions $\{f_{i}\}_{i=1}^{d}$. Our strategy is described in Algorithm 3 and works sequentially as follows: Given an input matrix $R$ (the residual after $k$ steps of RSPA), for $i=1,2,\dots,d$:

step 3. Identify the index $k(i)$ that correspond to the data point $R(:,k(i))$ that maximizes $f_{i}(x)$.

steps 4-6. Let $u_{i}=R(:,k(i))/||R(:,k(i))||_{2}$, and let $R_{i}=(I-u_{i}u_{i}^{T})R$ be the projection of $R$ onto the orthogonal complement of $u_{i}$. The $\ell_{p}$ norm of the $\ell_{2}$ norms of the columns of $R_{i}$, denoted $e(i)$, will allow us to select the index among $\{k(i)\}_{i=1}^{d}$ such that this norm is minimized (step 11).

steps 7-9. Identify the index $k^{\prime}(i)$ corresponding to the column of $R_{i}$ with largest norm. We choose $\alpha_{i}$ such that

[TABLE]

The value of $\alpha_{i}$ that achieves this is given in Lemma 1 (see below). This guarantees that, at the next step, $f_{i+1}$ will not identify $k(i)$ again since there is at least one data point with value $\beta$ times larger for $f_{i+1}$. Note that we simplify the computation of $f_{i+1}(R(:,k))$ by introducing the matrix $Y$ initialized as $Y=R$ and updated at each step as $Y\leftarrow P_{i+1}Y$ (step 9) so that $f_{i+1}(R(:,k))=||Y(:,k)||_{2}^{2}$ for all $k$.

Lemma 1**.**

Let $x$ and $y$ be two non-zero vectors not multiple of one another, $u=x/||x||_{2}$ and $||x||_{2}>||y||_{2}$. For $\beta>1$,

[TABLE]

is the unique solution for $\alpha\in(0,1)$ of

[TABLE]

Proof.

Using $||(I-\alpha uu^{T})z||_{2}^{2}=||z||_{2}^{2}-2\alpha(u^{T}z)^{2}+\alpha^{2}(u^{T}z)^{2}$, the solution of the above problem is a root of $\Delta\alpha^{2}-2\Delta\alpha+\beta||x||_{2}^{2}-||y||_{2}^{2}$ where $\Delta=\beta(u^{T}x)^{2}-(u^{T}y)^{2}$. Note that $u^{T}x=||x||_{2}>u^{T}y$ since $u^{T}y\leq||y||_{2}<||x||_{2}$ hence $\Delta>\beta||x||_{2}^{2}-||y||_{2}^{2}$. We obtain $\alpha^{*}=\frac{\Delta-\sqrt{\Delta^{2}-\Delta(\beta||x||_{2}^{2}-||y||_{2}^{2})}}{\Delta}\in(0,1)$. ∎

Computational cost

It can be checked that SPA runs in $O(mnr)$ operations [9], while RSPA runs in $O(mnrd)$. The main computational cost lies in the selection and projection steps, each in $O(mn)$ operations.

4 Numerical experiments

We write RSPA($d$,$p$,$\beta$) to refer to RSPA with parameters ($d$,$p$,$\beta$). The code is available from

https://sites.google.com/site/nicolasgillis/code.

4.1 Synthetic data sets

Let $r=10$, $n=1000$ and the value of $m$ is varied from 10 to 50. Each entry of $W\in\mathbb{R}^{m\times r}$ is generated randomly using the uniform distribution $\mathcal{U}(0,1)$ in the interval $[0,1]$. Each entry of $H\in\mathbb{R}^{r\times n-r}$ is generated in the same way, but then each column of $H$ is normalized so that Assumption 1 holds: $H(:,j)\leftarrow H(:,j)/||H(:,j)||_{1}$. Finally, we take $X=W[I_{r},H]$ to which we add 10 outliers whose entries are generated randomly using the normal distribution $\mathcal{N}(0,1)$ of mean zero and variance 1. For each value of $m$, we generate 100 such matrices. Figure 1 reports the percentage of correctly identified columns of $W$ by SPA and by RSPA with various combinations of the parameters, with $d\in\{10,20,40,80\}$, $p\in\{1,2\}$ and $\beta\in\{2,4,8\}$.

We observe the following:

• •

When $m$ is small, no algorithm is able to recover all columns of $W$. The reason is that the outliers and the columns of $W$ are less separated (they are linearly dependant for example when $m=r=10$).

• •

SPA performs very poorly. The reason is that the norm of the outliers is larger than that of the columns of $W$; in fact, $\mathbb{E}(x^{2})=1/3$ for $x\sim\mathcal{U}(0,1)$ while $\mathbb{E}(x^{2})=1$ for $x\sim\mathcal{N}(0,1)$. Note that any other greedy algorithm such as VCA [11] or the successive nonnegative projection algorithm (SNPA) [6] would fail as well.

• •

RSPA does not perform well when $p=2$ as it is more sensitive to large entries in the residual hence to outliers. We have also tried $p=0.5$ and it performed similarly as $p=1$.

• •

RSPA performs best for $\beta=4$. The parameter $\beta$ should not be chosen too large as it makes $P_{i}$ ill conditioned since $\alpha_{i}$ will be close to 1, nor too small as it does not provide a good diversification.

• •

RSPA does not perform well when $d$ is too small ($d\leq 20$): in that case, RSPA is not able to avoid outliers.

To summarize, RSPA performs best when $d\geq 40$, $p=1$ and $\beta=4$: more than 99% of the columns of $W$ are correctly identified for $m\geq 25$. Ideally, $d$ should be of the order of the number of outliers. In fact, if $d$ is smaller than the number of outliers, the diversification procedure could only identify outliers hence fail. This is particularly crucial when the outliers have a larger norm than the inliers (see also the experiment on the San Diego hyperspectral image). However, choosing $d$ too large makes RSPA slower as it runs in $O(mnrd)$ operations.

4.2 Hyperspectral images

Let us compare SPA and RSPA on three widely used hyperspectral images (HSIs) that are described for example in [8]:

• •

Urban: 162 spectral bands, $307\times 307$ pixels and 6 endmembers. It contains a few outliers that correspond to materials present in small proportions.

• •

San Diego: 158 spectral bands, $400\times 400$ pixels and 8 endmembers. It contains quite a few outliers (see Figure 2).

• •

Cuprite: 188 spectral bands, $250\times 191$ pixels and 15 endmembers. It does not contain outliers.

Table 1 reports the relative approximation error

[TABLE]

for the index sets $\mathcal{K}$ extracted by SPA, RSPA(10,1,4) and RSPA(20,1,4).

We observe the following:

• •

For Urban, RSPA variants allow to slightly reduce the relative error compared to SPA. The gain is appreciable (from 9.58% to 7.65% for RSPA(10,1,4) and to 6.66% for RSPA(20,1,4)) but not significant because the outliers are endmembers corresponding to materials present in small proportions and sharing similarities with the main endmembers.

• •

For San Diego, RSPA allows a significant reduction of the relative error; from 12.62% to 6.63% for RSPA(10,1,4) and to 6.03% for RSPA(20,1,4). Figure 2 displays the 8 endmembers extracted by each algorithm. We observe on Figure 2 that SPA identifies 5 outliers, RSPA(10,1,4) only 2, and RSPA(20,1,4) none. This confirms our observations made on synthetic data sets: the parameter $d$ should be chosen properly so as to allow RSPA to avoid extracting outliers. Note however that the relative errors of RSPA(10,1,4) and RSPA(20,1,4) are relatively close because the outliers have a very large norm.

• •

For Cuprite, SPA and RSPA provide comparable results because of the absence of outliers.

• •

In terms of computational time, RSPA is between $2d$ to $3d$ times slower than SPA. This is expected since RSPA requires $O(d)$ times more operations than SPA.

5 Conclusion

We have proposed a new variant of SPA, namely Robust SPA (RSPA), which is robust to outliers by taking into account the residual error to identify important columns in the data set, while remaining robust in low-noise settings. We have illustrated the effectiveness of RSPA on synthetic data sets and hyperspectral images. A similar enhancement could be brought to other greedy separable NMF algorithms, such as VCA and SNPA. Further work includes a thorough analysis of the behavior of RSPA under different choices of the parameters in various conditions, as well as a rigorous robustness analysis of RSPA with explicit error bounds depending on the noise level and the number of outliers.