\section{Experiments}
\label{sec:xp}


The code and data are available online\footnote{\url{https://gitlab.com/nnadisic/ssnmf}}.

All experiments have been performed on a personal computer with an i5
processor, with a clock frequency of 2.30GHz.
All algorithms are single-threaded.
All are implemented in Matlab, except the sparse NNLS solver \arbo{}, which is implemented in C++ with a Matlab MEX interface. 

As far as we know, no algorithm other than \brassens{} can tackle SSNMF with more than one interior point (see Section~\ref{sec:relwork}) hence comparisons with existing works are unfortunately limited.
For example, separable NMF algorithms can only identify the exterior vertices; see Section~\ref{sec:intro}.
However, we will compare \brassens{} to SNPA on a real multispectral image in Section~\ref{sec:xp-real}, to show the advantages of the SSNMF model over separable NMF.
In Section~\ref{sec:xp-synth}, we illustrate the correctness and efficiency of \brassens{} on synthetic data sets. 


\subsection{Synthetic Data Sets}\label{sec:xp-synth}

In this section, we illustrate the behaviour of \brassens{} in different experimental setups. 
The generation of a synthetic data set is done as follows: for a given number of dimensions $m$, number of vertices $r$, number of data points $n$, and data sparsity $k$, we generate matrices $W \in \mathbb{R}_+^{m \times r}$ such that the last $r-m$ columns of $W$ are linear combinations of the first $m$ columns, and $H \in \mathbb{R}_+^{r \times n}$ such that $H = [ I_r , H' ]$ and $\|H(:,j)\|_0 \leq k$ for all $j$.
We use the uniform distribution in the interval [0,1] to generate random numbers (columns of $W$ and columns of $H$), and then normalize the columns of $W$ and $H$ to have unit $\ell_1$ norm.
We then compute $M = W H$.
This way, the matrix $M$ is $k$-sparse $r$-separable, with $r$ vertices, of which $r-m$ are interior vertices (in fact, the first $m$ columns of $W$ are linearly independent with probability one as they are generated randomly). 
We then run \brassens{} on $M$, with the parameter $k$, and no noise-tolerance.
For a given setup, we perform 30 rounds of generation and solving, and we measure the median of the running time and the median of the number of candidates extracted by kSSNPA.
This number of candidates corresponds to $| \mathcal{J}' |$ in \cref{alg:ssnmf}, that is, the number of interior points selected by kSSNPA as potential interior vertices.
Note that a larger number of candidates only results in an increased computation time, and does not change the output of the algorithm which is guaranteed to extract all vertices (\cref{th:recovery}).

\cref{fig:xpnvar} shows the behaviour of \brassens{} when $n$ varies, with fixed $m=3$, $k=2$, and $r=5$.
To the best of our knowledge, this case is not handled by any other algorithm in the literature.
Both the number of candidates and the run time grow slower than linear.

The irregularities in the plot are due to the high variance between runs.
Indeed, if vertices are generated in a way that some segments between vertices are very close to each other, \brassens{} typically selects more candidates before identifying all columns of $W$.  
\begin{figure}[ht!] 
    \centering
    \input{img/xp-nvar}
    \caption{Results for \brassens{} on synthetic data sets for different values of $n$, with fixed $m=3$, $k=2$, and $r=5$ with 3 exterior and 2 interior vertices. The values showed are the medians over 30 experiments.}\label{fig:xpnvar}
\end{figure}

In \cref{tab:synth2} we compare the performance of \brassens{} for several sets of parameters. 
The number of candidates grows relatively slowly as the dimensions $(m,n)$ of the problem increase, showing the efficiency of the screening performed by kSSNPA. 
However, the run time grows rather fast when the dimensions $(m,n)$ grow.  
This is because, not only the number of NNLS subproblems to solve increase, but also their size. 
\setlength{\tabcolsep}{6pt}
\begin{table}[ht]
\centering
\caption{Results for \brassens{} on synthetic data sets (median over 30 experiments).} 
\label{tab:synth2}
\begin{tabular}{c|c|c|c|c|c}
m & n  & r  & k & Number of candidates & Run time in seconds \\ \hline
3 & 25 & 5  & 2 & 5.5                  & 0.26                \\
4 & 30 & 6  & 3 & 8.5                  & 3.30                \\
5 & 35 & 7  & 4 & 9.5                  & 38.71               \\
6 & 40 & 8  & 5 & 13                   & 395.88              \\

\end{tabular}
\end{table}
In all cases, as guaranteed by Theorem~\ref{th:recovery}, \brassens{} was able to correctly identify the columns of $W$. Again, as far as we know, no existing algorithms in the literature can perform this task. 

To summarize, our synthetic experiments show the efficiency of the screening done by kSSNPA, and the capacity of \brassens{} to handle medium-scale data sets.


\subsection{Blind Multispectral Unmixing}\label{sec:xp-real}

A multispectral image is an image composed of various wavelength ranges, called \emph{spectral bands}, where every pixel is described by its \emph{spectral signature}.
This signature is a vector representing the amount of energy measured for this pixel in every considered spectral band.
Multispectral images usually have a small number of bands (between 3 and 15).
These bands can be included or not in the spectrum of visible light.
In the NMF model, if the columns of $M$ are the $n$ pixels of the image, then its rows represent the $m$ spectral bands.

The unmixing of a multispectral image consists in identifying the different materials present in that image.
When the spectral signatures of the materials present in the image are unknown, it is referred to as \emph{blind unmixing}.
The use of NMF for blind unmixing of multispectral images relies on the linear mixing model, that is, the assumption that the spectral signature of a pixel is the linear combination of the spectral signatures of the material present in this pixel.
This corresponds exactly to the NMF model, which is therefore able to identify both the materials present in the image ($W$) and the proportions/abundances of materials present in every pixel ($H$); see~\cite{bioucas2012hyperspectral,ma2013signal} for more details. 

Let us apply \brassens{} to the unmixing of the well-known Urban satellite image~\cite{fyzhu2014hyperspectraldata}, composed of 
$309 \times 309$ pixels. The original cleaned image has 162 bands, but we only keep 3 bands, namely the bands 2, 80, and 133 -- these were obtained by selecting different bands with SPA applied on $M^T$ -- to obtain a data set of size $3 \times \num{94249}$. The question is: can we still recover materials by using only 3 bands? (The reason for this choice is that this data set is well known and the ground truth is available, which is not the case of most multispectral images with only 3 bands.)  
We first normalize all columns of $M$ so that they sum to one.
Then, we run \brassens{} with a sparsity constraint $k = 2$ (this means that we assume that a pixel can be composed of at most 2 materials, which is reasonable for this relatively high resolution image) and a noise-tolerance threshold of 4\%; this means that we stop SNPA and kSSNPA when 
$\|M - M(:,\mathcal{J}) H\|_F \leq 0.04 \|M\|_F$. 
\brassens{} extracts 5 columns of the input matrix.  
For comparison, we run SNPA with $r = 5$. 
Note that this setup corresponds to underdetermined blind unmixing, because $m = 3 < r = 5$. 
It would not be possible to tackle this problem using standard NMF algorithms (that would return a trivial solution such as $M = I_3 M$). 
It can be solved with SNPA, but SNPA cannot identify interior vertices. 

SNPA extracts the 5 vertices in \num{3.8} seconds.
\brassens{} extracts 5 vertices, including one interior vertex, in \num{33} seconds. 
The resulting abundance maps are showed in \cref{fig:mapbrassens}.
They correspond to the reshaped rows of $H$, hence they show which pixel contains which extracted material (they are more easily interpretable than the spectral signatures contained in the columns of $W$).  
The materials they contain are given in \cref{tab:unmixing}, using the ground truth from~\cite{zhu2017hyperspectral}.
We see that \brassens{} produces a better solution, as the materials present in the image are better separated: the first three abundance maps of \brassens{} are sparser and correspond to well-defined materials. The last two abundances maps of SNPA and of \brassens{} are similar but extracted in a different order. 
The running time of \brassens{} is reasonable, although ten times higher than SNPA. 


\begin{figure}[ht!]
    \centering
    \includegraphics[width=\textwidth]{img/mapsnpa.png}
    \vspace{0.1cm}\\ 
    \includegraphics[width=\textwidth]{img/mapbrassens.png} 
    \caption{Abundances maps of materials (that is, reshaped rows of $H$) extracted by SNPA (top) and by \brassens{} (bottom) in the Urban image with only 3 spectral bands.}\label{fig:mapbrassens}
\end{figure}

\begin{table}[ht!] 
\centering
\caption{Interpretation of the unimixing results from Figure~\ref{fig:mapbrassens}.}\label{tab:unmixing}
\begin{tabular}{c|c|c}
Image & Materials extracted by SNPA & Materials extracted by \brassens{} \\ \hline
1     & Grass + trees + roof tops   & Grass + trees                   \\
2     & Roof tops 1                 & Roof tops 1                     \\
3     & Dirt + road + roof tops     & Road                            \\
4     & Dirt + grass                & Roof tops 1 and 2 + road        \\
5     & Roof tops 1 + dirt + road   & Dirt + grass                   
\end{tabular}
\end{table}
{"config": "arxiv", "file": "2006.07553/60xp.tex"}