∎

11institutetext: ✉ Nicolas Nadisic
<EMAIL_ADDRESS>
1 Department of Mathematics and Operational Research, University of Mons,
Mons, Belgium

# A Homotopy-based Algorithm for Sparse Multiple Right-hand Sides Nonnegative
Least Squares ††thanks: NN and NG acknowledge the support by the European
Research Council (ERC starting grant No 679515), and by the Fonds de la
Recherche Scientifique - FNRS and the Fonds Wetenschappelijk Onderzoek -
Vlanderen (FWO) under EOS project O005318F-RG47.

Nicolas Nadisic1 Arnaud Vandaele1 Nicolas Gillis1

(Received: date / Accepted: date)

###### Abstract

Nonnegative least squares (NNLS) problems arise in models that rely on
additive linear combinations. In particular, they are at the core of
nonnegative matrix factorization (NMF) algorithms. The nonnegativity
constraint is known to naturally favor sparsity, that is, solutions with few
non-zero entries. However, it is often useful to further enhance this
sparsity, as it improves the interpretability of the results and helps
reducing noise. While the $\ell_{0}$-“norm”, equal to the number of non-zeros
entries in a vector, is a natural sparsity measure, its combinatorial nature
makes it difficult to use in practical optimization schemes. Most existing
approaches thus rely either on its convex surrogate, the $\ell_{1}$-norm, or
on heuristics such as greedy algorithms. In the case of multiple right-hand
sides NNLS (MNNLS), which are used within NMF algorithms, sparsity is often
enforced column- or row-wise, and the fact that the solution is a matrix is
not exploited. In this paper, we first introduce a novel formulation for
sparse MNNLS, with a matrix-wise $\ell_{0}$ sparsity constraint. Then, we
present a two-step algorithm to tackle this problem. The first step uses a
homotopy algorithm to produce the whole regularization path for all the
$\ell_{1}$-penalized NNLS problems arising in MNNLS, that is, to produce a set
of solutions representing different tradeoffs between reconstruction error and
sparsity. The second step selects solutions among these paths in order to
build a sparsity-constrained matrix that minimizes the reconstruction error.
We illustrate the advantages of our proposed algorithm for the unmixing of
facial and hyperspectral images.

###### Keywords:

nonnegative least squares sparsity nonnegative matrix factorization homotopy
algorithm

††journal: Preprint

##  1 Introduction

Nonnegative least squares (NNLS) problems arise in many applications where
data points can be represented as additive linear combinations of meaningful
components (Lee and Seung, 1997). For instance,

  * •

In facial images, the faces are the nonnegative linear combination of facial
features such as eyes, noses and lips (Lee and Seung, 1999).

  * •

In hyperspectral images, the spectral signature of a pixel is the nonnegative
linear combination of the spectral signature of the materials it contains
(Bioucas-Dias et al, 2012).

NNLS problems are also at the core of many alternating approaches to solve
nonnegative matrix factorization (NMF) (Gillis, 2014b). The standard NNLS
problem using least squares can be formulated as follows

$\min\limits_{x\in\mathbb{R}^{r}}\frac{1}{2}\|Ax-b\|_{2}^{2}\quad\text{ such
that }\quad x\geq 0,$ (1)

where $A\in\mathbb{R}^{m\times r}$, and $b\in\mathbb{R}^{m}$.

###  1.1 Sparsity and NNLS

The nonnegativity constraint is known to naturally produce sparse solutions,
that is, solutions with few non-zero entries (Foucart and Koslicki, 2014).
Sparsity often improves the interpretability of the results by modelling data
points as combinations of only a few components. For example, in hyperspectral
unmixing, that is, the task of identifying materials in a hyperspectral image,
sparsity means that a pixel contains only a few materials. Unfortunately, the
sparsity of the solution to a NNLS problem is not guaranteed in general,
whereas controlling it can be helpful in many applications. For this reason,
numerous techniques have been developed to favor sparsity.

A natural sparsity measure is the $\ell_{0}$-“norm”, defined as the number of
non-zero entries in a given vector, $\|x\|_{0}=|\\{i:x_{i}\neq 0\\}|$. Given a
positive integer $k$, a vector $x$ is $k$-sparse if $\|x\|_{0}\leq k$. A
sparsity-constrained variant of Problem (1), referred to as _$k$ -sparse
NNLS_, is the following

$\min\limits_{x\in\mathbb{R}^{r}}\frac{1}{2}\|Ax-b\|_{2}^{2}\quad\text{ such
that }\quad x\geq 0\;\text{ and }\;\|x\|_{0}\leq k.$ (2)

Albeit intuitive and easily interpretable, the constraint of $k$-sparsity is
hard to enforce in practice, because the $\ell_{0}$-“norm” is non-convex and
non-smooth.

For this reason, practitioners often prefer other sparsity-inducing
approaches, such as $\ell_{1}$-regularization. The $\ell_{1}$-norm, defined as
$\|x\|_{1}=\sum_{i=1}^{r}|x_{i}|$, is a convex surrogate of the
$\ell_{0}$-“norm”, it is therefore easier to optimize in practice while being
a good sparsity measure. The $\ell_{1}$-regularization consists in penalizing
the solution in the objective function of (1), leading to the following
problem, referred to as $\ell_{1}$-NNLS,

$\min\limits_{x\in\mathbb{R}^{r}}\frac{1}{2}\|Ax-b\|_{2}^{2}+\lambda\|x\|_{1}\quad\text{
such that }\quad x\geq 0.$ (3)

Note that, without the nonnegativity constraint, this is the well-known LASSO
model (Tibshirani, 1996). Problem (3) can be thought of as the weighted-sum
form of a multiobjective problem, where the objectives are minimizing the
reconstruction error $\frac{1}{2}\|Ax-b\|_{2}^{2}$ on one hand, and minimizing
the $\ell_{1}$-norm of the solution $\|x\|_{1}$ on the other hand. The former
measures the quality of the solution, while the latter favors its sparsity.
The parameter $\lambda$ controls the trade-off between the two objectives.

Despite its popularity, this technique suffers from several drawbacks. In
particular, there is no explicit relation between the parameter $\lambda$ and
the sparsity of the solution, hence choosing an appropriate value for
$\lambda$ can be tricky, and often involves a tedious trial-and-error process.
Also, the $\ell_{1}$-penalty introduces a bias. Although there exist
theoretical guarantees for support recovery, such as the _Exact Recovery
Condition_ (Tropp, 2006), they are restrictive and often not realistic in
practice.

To overcome these issues, _homotopy_ methods have been introduced that
generate the full _regularization path_ of a given $\ell_{1}$-NNLS problem,
that is, the set of the solutions for all possible values of $\lambda$. They
allow the user to choose the relevant solution within this path, each solution
representing a different trade-off between sparsity and reconstruction error.
The first homotopy method has been introduced by Osborne et al (2000), for
sparse least squares with no nonnegativity constraint. Variants have been
developed by Efron et al (2004) and Donoho and Tsaig (2008). Kim et al (2013)
introduced a variant to deal specifically with the nonnegativity constraint.
We detail the homotopy method for NNLS in Section 3.

###  1.2 Sparsity and multiple right-hand sides NNLS

In many applications, we have to deal with multiple right-hand sides NNLS
problems (MNNLS), that is, problems of the form

$\min\limits_{H}\frac{1}{2}\|M-WH\|_{2}^{2}\quad\text{ such that }\quad H\geq
0,$ (4)

where $M\in\mathbb{R}^{m\times n}$, $W\in\mathbb{R}^{m\times r}$, and
$H\in\mathbb{R}^{r\times n}$. We note $M(:,j)$ the $j$-th column of the matrix
M. Problem (4) can be decomposed into $n$ NNLS subproblems of the form (1),
where $H(:,j)$, $W$, and $M(:,j)$ correspond to $x$, $A$, and $b$,
respectively. For example, in the unmixing of a hyperspectral image, every
column $M(:,j)$ represents a pixel, and the corresponding column $H(:,j)$
represents its composition, in terms of the abundances of the $r$ materials
whose spectral signatures are the columns of $W$. In this work, for the sake
of conciseness, we focus only on the (sparse) optimization of $H$, but all
concepts and algorithms can be applied symmetrically on $W$.

This is related to the nonnegative matrix factorization (NMF) problem, of the
form

$\min\limits_{W,H}\frac{1}{2}\|M-WH\|_{F}^{2}\quad\text{ such that }\quad
W\geq 0\text{ and }H\geq 0,$ (5)

in which we aim at finding the factors $W$ and $H$, given $M$ and a
factorization rank $r$. The usual optimization scheme for NMF consists in
alternatively optimizing one factor while fixing the other, which is
equivalent to solving MNNLS subproblems.

To encourage sparsity in MNNLS, one can apply a sparse NNLS model column-wise,
leading to

$\min\limits_{H}\frac{1}{2}\|M-WH\|_{2}^{2}\quad\text{ such that }\quad H\geq
0\text{ and }\|H(:,j)\|_{0}\leq k\text{ for all }j.$ (6)

However, in some applications, setting the sparsity parameter $k$ is tricky,
as the relevant value can vary for different columns. For example, in
hyperspectral unmixing, pixels will be composed of different numbers of
materials. Therefore, one can consider a more global approach, such as

$\min\limits_{H}\frac{1}{2}\|M-WH\|_{2}^{2}\quad\text{ such that }\quad H\geq
0\text{ and }\|H\|_{0}\leq q,$ (7)

where $\|H\|_{0}=\sum_{j}\|H(:,j)\|_{0}$ and $q$ is a matrix-wise sparsity
parameter, hence enforcing an _average_ sparsity $q/n$ for the columns of $H$.

###  1.3 Contribution and outline of the paper

In this work, we present an approximate method to solve Problem (7). Our
approach consists in using a homotopy method to compute a regularization path
for every column of $H$, and then exploit these paths with a dedicated
optimization algorithm to build a solution to Problem (7).

The paper is organized as follows. In Section 2, we present existing
approaches for sparse NNLS and sparse MNNLS . In Section 3, we detail the
homotopy principle and algorithm. In Section 4, we present the main
contribution of this work, that is, an algorithm to solve Problem (7). We
illustrate the effectiveness of our proposed algorithm with experiments on
real-life facial and hyperspectral image datasets in Section 5, and conclude
in Section 6.

##  2 Related work

We discussed in Section 1.1 the methods based on the $\ell_{1}$-norm, and
their limitations. In this section, we first discuss the existing works that
solve the $\ell_{0}$-constrained problem, as defined in Problem (2). We then
review the most common approaches for sparse MNNLS.

###  2.1 $k$-sparse NNLS

Greedy algorithms are one of the most popular approaches for solving Problem
(2). They start with an empty support ($x_{i}=0$ for all $i$), and select
components one by one to enrich the support, until the target sparsity $k$ is
reached. The selection of a component is done greedily, by choosing the
component maximizing the decrease of the residual error. Orthogonal versions
of these methods, such as _Orthogonal Least Squares_ (OLS) (Chen et al, 1989)
and _Orthogonal Matching Pursuit_ (OMP) (Pati et al, 1993) make sure a
component can be selected only once. Nonnegative variants have been proposed
recently; see for example Nguyen et al (2019) and the references therein. In
general, these algorithms do not give the globally optimal solution. Again,
theoretical recovery guarantees exist, but they are restrictive (Soussen et
al, 2013).

Other works were proposed for solving globally the $\ell_{0}$-constrained
problem, although not specifically with nonnegativity constraints. Notably,
Bienstock (1996) proposed a reformulation of Problem (2) to a mixed-integer
program (MIP), and a branch-and-cut algorithm to solve it, using continuous
relaxation at every iteration. Several contributions extended this approach,
for example under different constraints, or using new methods to solve the
continuous relaxation subproblem; see for example Bertsimas and Shioda (2009);
Bourguignon et al (2015) and the references therein. Mhenni et al (2019)
further extended it, using a homotopy algorithm to compute the relaxations.

###  2.2 Sparse MNNLS

Most approaches that tackle sparse MNNLS were actually introduced in the
context of sparse NMF. Since its very introduction by Lee and Seung (1999),
NMF is appreciated for the sparsity of the produced factors. A variety of
works have been proposed to further enhance this sparsity, making _sparse NMF_
one the most popular variants of NMF. Many authors worked on the
$\ell_{1}$-penalized formulation, notably Hoyer (2002); Eggert and Korner
(2004); Kim and Park (2007); Cichocki et al (2008); Gillis (2012).

To avoid the issues linked to the $\ell_{1}$-penalty, Hoyer (2004) introduced
a more explicit sparsity measure, based on the ratio between the
$\ell_{1}$-norm and the $\ell_{2}$-norm, and considered an NMF variant with a
column-wise constraint on this measure. Other works considered the
$\ell_{0}$-“norm” formulation, that can be decomposed into a series of
$k$-sparse NNLS subproblems. We can cite Aharon et al (2005); Morup et al
(2008); Peharz and Pernkopf (2012). Cohen and Gillis (2019) proposed a method
that solves _exactly_ the $k$-sparse NNLS subproblems, using a bruteforce
approach. Nadisic et al (2020) extended this work by replacing the bruteforce
subroutine by a dedicated branch-and-bound algorithm.

##  3 The homotopy algorithm

In this section, we review the homotopy method, which is central in our
proposed algorithm for solving Problem (7).

Given an $\ell_{1}$-NNLS problem, the homotopy algorithm computes sequentially
all optimal solutions for the different values of $\lambda$. In a nutshell, it
uses the KKT conditions (necessary conditions for optimality) to first find
the value $\lambda_{\max}$ for which the optimal solution of $\ell_{1}$-NNLS
is $x=0$ for any $\lambda\geq\lambda_{\max}$, and then to compute the next
smaller values of $\lambda$ for which the support (that is, the set of non-
zero entries) of the optimal solution changes (one zero entry becomes non-
zero, or the other way around). This is similar in spirit to active-set
methods, akin to the simplex algorithm for linear programming. These values of
$\lambda$ are called _breakpoints_ , between which the support of the optimal
solution does not change. We denote the first breakpoint $\lambda_{\max}$, and
the following ones $\lambda_{2}$, $\lambda_{3}$, $\ldots$; see Figs. 2 and 2
for an illustration.

$\lambda$$x_{i}=0$$x_{i}$$\lambda=0$$\lambda_{4}$$\lambda_{3}$$\lambda_{2}$$\lambda_{max}$$x_{1}$$x_{2}$$x_{3}$
Figure 1: Example of a solution path of a homotopy algorithm, depending on
$\lambda$, for a problem of 3 variables. Vertical dotted lines correspond to
breakpoints.

$\|Ax-b\|_{2}^{2}$0$\|x\|_{1}$$\lambda_{max}$$\|b\|_{2}^{2}$$\lambda_{2}$$\lambda_{3}$$\lambda_{4}$$\lambda=0$
Figure 2: Trade-off between the $\ell_{1}$-norm and the error
$\|Ax-b\|_{2}^{2}$ corresponding to the solution path in Fig. 2.

###  3.1 Optimality conditions

The homotopy algorithm uses the first-order optimality conditions, that is,
the KKT conditions, to determine the breakpoints and the supports of the
corresponding solutions. Because $x$ is nonnegative, we have
$\|x\|_{1}=\sum_{i}x_{i}=e^{T}x$, where $e$ is the vector of all ones whose
dimension will be clear from context. Therefore, the $\ell_{1}$-NNLS problem
can be written as follows

$\min\limits_{x\geq 0}f(x),\text{ where
}f(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}+\lambda e^{T}x.$ (8)

The KKT conditions are $x\geq 0$,

$\displaystyle\nabla
f(x)=\underbrace{A^{\top}A}_{P}x-\underbrace{A^{\top}b}_{\ell}+\lambda e$
$\displaystyle\geq 0,$ (9) $\displaystyle x_{i}(A^{\top}Ax-A^{\top}b+\lambda
e)_{i}$ $\displaystyle=0\text{ for all }i.$ (10)

Eq. 10 is the complementary condition; for every entry of $x$, either the
entry itself or the corresponding gradient entry is equal to zero. To simplify
the notation, let us define $P=A^{\top}A$ and $\ell=A^{\top}b$.

As $f$ is a convex function and the feasible set contains a Slater point
(e.g., $x=e$), the KKT conditions are necessary and sufficient. Therefore, any
solution $x$ satisfying them is optimal. Suppose we know the optimal support,
that is, the set $\mathcal{K}$ such that $x(\mathcal{K})>0$ and
$x(\bar{\mathcal{K}})=0$, where
$\bar{\mathcal{K}}=\\{1,2,\ldots,n\\}\setminus\mathcal{K}$. Then

$\displaystyle x(\mathcal{K})>0$ $\displaystyle\Rightarrow
P(\mathcal{K},\mathcal{K})x(\mathcal{K})-\ell(\mathcal{K})+\lambda e=0$
$\displaystyle\Rightarrow
x(\mathcal{K})=P(\mathcal{K},\mathcal{K})^{-1}(\ell(\mathcal{K})-\lambda
e)\geq 0,$ ($\mathcal{C}1$)

and

$Px-\ell+\lambda e\geq 0\Rightarrow
P(\bar{\mathcal{K}},\mathcal{K})x(\mathcal{K})-\ell(\bar{\mathcal{K}})+\lambda
e\geq 0.$ ($\mathcal{C}2$)

Replacing $x(\mathcal{K})$ in ($\mathcal{C}2$) by ($\mathcal{C}1$), we have

$P(\bar{\mathcal{K}},\mathcal{K})[P(\mathcal{K},\mathcal{K})^{-1}(\ell(\mathcal{K})-\lambda
e)]-\ell(\bar{\mathcal{K}})+\lambda e\geq 0.$

Let us simplify the notation. Let

$a_{\mathcal{K}}=P(\mathcal{K},\mathcal{K})^{-1}\ell(\mathcal{K}),\text{ and
}b_{\mathcal{K}}=P(\mathcal{K},\mathcal{K})^{-1}e.$ (11)

We can rewrite ($\mathcal{C}1$) as

$a_{\mathcal{K}}-\lambda b_{\mathcal{K}}\geq 0.$ ($\mathcal{C}1b$)

Note that the dimension of $a_{\mathcal{K}}$ and $b_{\mathcal{K}}$ is the
cardinality of $\mathcal{K}$. Let

$c_{\mathcal{K}}=P(\bar{\mathcal{K}},\mathcal{K})a_{\mathcal{K}}-\ell(\bar{\mathcal{K}}),\text{
and }d_{\mathcal{K}}=P(\bar{\mathcal{K}},\mathcal{K})b_{\mathcal{K}}-e.$

We can rewrite ($\mathcal{C}2$) as

$c_{\mathcal{K}}-\lambda d_{\mathcal{K}}\geq 0.$ ($\mathcal{C}2b$)

Note that the dimension of $c_{\mathcal{K}}$ and $d_{\mathcal{K}}$ is the
cardinality of $\bar{\mathcal{K}}$. Moreover, Equations ($\mathcal{C}1b$) and
($\mathcal{C}2b$) are linear in $\lambda$: Given $\mathcal{K}$, we can easily
compute $a_{\mathcal{K}}$, $b_{\mathcal{K}}$, $c_{\mathcal{K}}$,
$d_{\mathcal{K}}$.

###  3.2 Algorithm description

The goal of the homotopy algorithm is to compute breakpoints and their
corresponding supports. It starts with an empty support, corresponding to the
zero vector for any $\lambda\geq\lambda_{max}$, and iteratively adds or
removes entries to the support while decreasing the value of $\lambda$.

The first step to build the regularization path is to find the first
breakpoint $\lambda_{\max}$, that is, the minimum value of $\lambda$ such that
the solution is the zero vector. If the optimal solution is $x=0$, that is,
$\mathcal{K}=\varnothing$ and $\bar{\mathcal{K}}=\\{1,2,\ldots,n\\}$, then
from Eq. 9 we have $\lambda\geq\max\limits_{i}\ell_{i}$. Therefore, the first
breakpoint is

$\lambda_{\max}=\max\limits_{i}\ell_{i}=\max\limits_{i}(A^{\top}b)_{i}=\max\limits_{i}A(:,i)^{T}b.$

The index $i_{1}=\arg\max\limits_{i}\ell_{i}$ is the first to enter the
support111If two of more columns maximize the value of $\ell_{i}$, we have to
pick one to start the homotopy. If we always pick the one with smallest index,
then the algorithm behaves normally, except that at the next iteration we will
have $\lambda_{j+1}=\lambda_{j}$. This is similar to Bland’s rule for the
simplex algorithm. , thus for $\lambda_{2}\leq\lambda<\lambda_{\max}$ we have
$\mathcal{K}=\\{i_{1}\\}$.

From a given support $\mathcal{K}_{j}$, the next breakpoint
$\lambda_{j+1}\leq\lambda_{j}$ is the largest value of $\lambda$ that violates
one of the conditions ($\mathcal{C}1$) or ($\mathcal{C}2$). If
($\mathcal{C}1$) is violated then a variable will leave the support, that is,
a positive entry will become zero. Denoting $k^{*}$ the index of this entry,
we have $\mathcal{K}_{j+1}=\mathcal{K}_{j}\setminus\\{k^{*}\\}$. If
($\mathcal{C}2$) is violated then a variable will enter the support, that is,
a zero entry will become positive,
$\mathcal{K}_{j+1}=\mathcal{K}_{j}\cup\\{k^{*}\\}$.

Let us consider ($\mathcal{C}1$). We have $a_{\mathcal{K}}(k)-\lambda
b_{\mathcal{K}}(k)\geq 0$ for all $k$ so

$\lambda_{j+1}\geq\max_{\\{k|b_{\mathcal{K}}(k)<0\\}}\frac{a_{\mathcal{K}}(k)}{b_{\mathcal{K}}(k)}.$

Similarly, for ($\mathcal{C}2$) we have $c_{\mathcal{K}}(k)-\lambda
d_{\mathcal{K}}(k)\geq 0$ for all $k$ so

$\lambda_{j+1}\geq\max_{\\{k|d_{\mathcal{K}}(k)<0\\}}\frac{c_{\mathcal{K}}(k)}{d_{\mathcal{K}}(k)}.$

Therefore,

$\lambda_{j+1}=\max\Bigg{(}\underbrace{\max_{\\{k|b_{\mathcal{K}}(k)<0\\}}\frac{a_{\mathcal{K}}(k)}{b_{\mathcal{K}}(k)}}_{\text{Case
}1},\underbrace{\max_{\\{k|d_{\mathcal{K}}(k)<0\\}}\frac{c_{\mathcal{K}}(k)}{d_{\mathcal{K}}(k)}}_{\text{Case
}2}\Bigg{)}.$

The algorithm is detailed formally in Algorithm 1. Note that, inside the
algorithm loop, once a support $\mathcal{K}$ is identified, getting the
corresponding optimal solution is straightforward. From Eq. 11, if
$a_{\mathcal{K}}$ is nonnegative, then the unbiased optimal solution of the
NNLS problem is the vector $x^{*}$ such that $x^{*}(\bar{\mathcal{K}})=0$ and
$x^{*}(\mathcal{K})=a_{\mathcal{K}}$. If $a_{\mathcal{K}}$ has negative
entries, then we can compute the unbiased solution with a standard NNLS
solver, with $x^{*}(\bar{\mathcal{K}})=0$ and
$x^{*}(\mathcal{K})=\arg\min\limits_{x\geq
0}\|P(\mathcal{K},\mathcal{K})x-l(\mathcal{K})\|_{2}^{2}$.

1

Input: $A\in\mathbb{R}^{m\times r}$, $b\in\mathbb{R}^{m}$

Output: Breakpoints $\lambda_{j}$, corresponding supports $\mathcal{K}_{j}$
and solutions $x^{*}_{j}$, for all $j$

2

3$i_{1}\leftarrow\arg\max\limits_{i}\ell_{i}$

4 $K\leftarrow\\{i_{1}\\}$ ; $\lambda_{1}=\lambda_{\max}=\ell_{i_{1}}$,
$j\leftarrow 1$

5 while  _ $\lambda_{j}>0$ _ do

6 $j\leftarrow j+1$

7 Compute $a_{K}$, $b_{K}$, $c_{K}$, $d_{K}$

8
$\lambda_{C1}\leftarrow\max_{\\{k|b_{\mathcal{K}}(k)<0\\}}\frac{a_{\mathcal{K}}(k)}{b_{\mathcal{K}}(k)}$,
and $k^{*}_{C1}\leftarrow$ corresponding argmax

9
$\lambda_{C2}\leftarrow\max_{\\{k|d_{\mathcal{K}}(k)<0\\}}\frac{c_{\mathcal{K}}(k)}{d_{K}(k)}$,
and $k^{*}_{C2}\leftarrow$ corresponding argmax

10 if  _$\lambda_{C1}\geq\lambda_{C2}$_ then

11 $\lambda_{j}\leftarrow\lambda_{C1}$

12 $\mathcal{K}\leftarrow\mathcal{K}\setminus\\{k^{*}_{C1}\\}$

13 else if  _$\lambda_{C1} <\lambda_{C2}$_ then

14 $\lambda_{j}\leftarrow\lambda_{C2}$

15 $\mathcal{K}\leftarrow\mathcal{K}\cup\\{k^{*}_{C2}\\}$

16

17 $x^{*}_{j}(\bar{\mathcal{K}})\leftarrow 0$

18 if  _$a_{\mathcal{K}} >0$_ then

19 $x^{*}_{j}(\mathcal{K})\leftarrow a_{K}$

20 else

21 $x^{*}_{\mathcal{K}}=\arg\min\limits_{x\geq
0}\|P(\mathcal{K},\mathcal{K})x-\ell(\mathcal{K})\|_{2}^{2}$

22

23

Algorithm 1 Homotopy algorithm for sparse NNLS

###  3.3 Computational cost

As explained by Kim et al (2013), the time complexity of one iteration of the
homotopy algorithm is the same as one iteration of the standard active-set
algorithm (Lawson and Hanson, 1995), that is, $\mathcal{O}(r^{3})$ operations.
It is dominated by the computation of $a_{\mathcal{K}}$, that entails solving
a linear system in at most $r$ variables. The number of iterations equals the
number of breakpoints, which is in practice similar to those of the active-
set. In the worst case, active-set methods might require an exponential number
of iterations, up to $\mathcal{O}(2^{r})$, as the simplex algorithm for linear
programming. However, in practice, we have observed that it typically requires
much less iterations, of the order of $\mathcal{O}(r)$. In particular, when
$P^{-1}$ is diagonally dominant, we have $b_{\mathcal{K}}>0$, so
($\mathcal{C}1$) is never violated. As a result, when $\lambda$ decreases, no
positive entry of $x$ becomes zero. We only add entries to the support, so the
homotopy algorithm will be done in at most $r$ iterations. In practice, even
when this condition is not met, we have observed that adding entries to the
support happens far more often than removing entries.

As the homotopy algorithm solves a series of $\ell_{1}$-penalized NNLS
problems, there exist conditions under which it is guaranteed to recover the
correct supports, that is, the supports of the solutions of the corresponding
$\ell_{0}$-constrained NNLS problems; see Itoh et al (2017) and the references
therein.

##  4 Proposed algorithm for sparse MNNLS

In this section, we present the main contribution of our work, that is, how we
leverage the NNLS homotopy algorithm to solve Problem (7). First of all, given
the data matrix $M$ and the dictionary $W$, we run the homotopy algorithm for
every column of $M$, that is, with input $A=W$ and $b=M(:j)$ for all $j$. This
gives us the set of supports of the optimal solutions of the
$\ell_{1}$-regularized NNLS for every column of $M$, for all values of
$\lambda$. From these sets, we build a cost matrix $\mathcal{C}$ where each
column represents a column $j$ of $H$, each row represents a $k$-sparsity
between $0$ and $r$, and each entry is the error $\|M(:,j)-WH(:,j)\|_{2}^{2}$
of the $k$-sparse solution of column $j$. When this is done, we need to select
one solution per column to build the solution matrix $H$, that is, we need to
choose the sparsity level for each column of $H$. This selection is a
combinatorial problem, similar to an assignment problem. It can be tackled in
many different ways. For instance, one can trivially select the $k$-sparse
solution for every column, thus reducing to Problem (6). In the following, we
propose a global selection strategy that solves the selection subproblem
exactly.

We call our algorithm222The name stands for “shamans is a Homotopy-based
Algorithm for Matrix Approximation under Nonnegativity and Sparsity” shamans,
and detail it in Algorithm 2. Note that we do not use the values of the
breakpoints $\lambda$; we only use the regularization path, that is, the
solutions $x^{*}$ and their supports $\mathcal{K}$. Let us explain how shamans
proceeds.

1

Input: Data matrix $M\in\mathbb{R}^{m\times n}$, dictionary
$W\in\mathbb{R}^{m\times r}$, sparsity target $q\in\mathbb{N}$

Output: $H^{*}\in\mathbb{R}^{r\times n}$, $H^{*}=\arg\min\limits_{H\geq
0}\|M-WH\|_{F}^{2}\text{ s.t. }\|H\|_{0}\leq q+r-1$

2

3Init matrix $\mathcal{C}\in\mathbb{R}^{r\times n}$, $\forall i,\forall
j,\mathcal{C}(i,j)\leftarrow\infty$

4 Init solutions $Sol_{i,j}$

5foreach  _ $j\in\\{1,\ldots,n\\}$ _ do

6 $\\{K_{t},x^{*}_{t}\\}_{\forall t}\leftarrow$ homotopy($W$,$M(:,j)$)

7 foreach  _ $t$ _ do

8 $k\leftarrow|K_{t}|$

9 $err\leftarrow\|M(:,j)-Wx^{*}_{t}\|_{2}^{2}$

10 if  _$err <\mathcal{C}(k,j)$ _ then

11 foreach  _$i\in\\{k,\ldots,r\\}$ _ do

12 $\mathcal{C}(i,j)\leftarrow err$

13 $Sol_{i,j}\leftarrow x^{*}_{t}$

14

15

16

17

18Init cursors, $\forall j,Cur_{j}\leftarrow 0$

19 Init matrix $\Delta\mathcal{C}\in\mathbb{R}^{r\times n}$

20 Init matrix $\mathcal{G}\in\mathbb{R}^{r\times n}$

21

22foreach  _$i\in\\{1,\ldots,r\\},j\in\\{1,\ldots,n\\}$ _ do

23 $\Delta\mathcal{C}(i,j)\leftarrow\mathcal{C}(i-1,j)-\mathcal{C}(i,j)$

24

25foreach  _$i\in\\{1,\ldots,r\\},j\in\\{1,\ldots,n\\}$ _ do

26 $\mathcal{G}(i,j)\leftarrow(\sum_{l=1}^{i}\Delta\mathcal{C}(l,j))/i$

27

28$it\leftarrow 0$

29 while  _$it <q$_ do

30 $i^{*},j^{*}\leftarrow\arg\max\mathcal{G}$

31 $it\leftarrow i^{*}-Cur_{j^{*}}$

32 $Cur_{j^{*}}=i^{*}$

33 foreach  _$i\in\\{i^{*}+1,\ldots,r\\}$ _ do

34
$\mathcal{G}(i,j^{*})\leftarrow\frac{\sum_{l=i^{*}+1}^{i}\Delta\mathcal{C}(l,j^{*})}{(i-i^{*})}$

35

36

37foreach  _$j\in\\{1,\ldots,n\\}$ _ do

38 $H^{*}(:,j)=Sol_{Cur_{j},j}$

39

Algorithm 2 shamans

First, between Algorithms 2 and 2, the homotopy algorithm is used to compute
$k$-sparse solutions for the columns of $H$, as explained above. On Algorithm
2, we begin to loop over the columns of $M$. On Algorithm 2, we call the
homotopy algorithm and store the resulting solutions and their supports.
Between Algorithm 2 and Algorithm 2, we build the cost matrix $\mathcal{C}$
such that every column $\mathcal{C}(:,j)$ corresponds to the subproblem
$\min\limits_{x\geq 0}\|M(:,j)-Wx\|_{2}^{2}$, that is, to the computation of
the column $H(:,j)$, and every row corresponds to a value of $k$, between 0
and $r$. An entry of $\mathcal{C}$ thus represents the reconstruction error of
the $k$-sparse solution for the column $H(:,j)$ obtained with the homotopy.
The corresponding $Sol_{k,j}$ stores this solution. The loop on Algorithm 2
allows, when updating a $k$-sparse solution, to also update the
$k^{\prime}$-sparse solutions, with $k^{\prime}$ greater than $k$. This is
necessary so that, if a $k$-sparse solution has a lower reconstruction error
than a $k^{\prime}$-sparse solution with $k^{\prime}>k$ (although this case is
unlikely to happen), then the entry corresponding to $k^{\prime}$ will contain
the $k$-sparse solution. It is also useful when, for a given column and for
some $k$, a $k$-sparse solution does not exist because the homotopy algorithm
does not necessarily generate $k$-sparse solutions for all $k=0,1,2,\dots,r$.
Note that, with $k<k^{\prime}$, a $k$-sparse vector $x_{k}$ is also
$k^{\prime}$-sparse, because $\|x_{k}\|_{0}\leq k<k^{\prime}$.

The selection step, between Algorithms 2 and 2, begins with the initialization
of the cursors. There is one cursor $Cur_{j}$ for every column $H(:,j)$, and
it indicates the number of non-zero entries selected in the current solution
(that is, the $k$-sparsity of the currently selected solution) for this
column. We begin with all cursors at zero, meaning that the $0$-sparse
solution (the zero vector) is selected for each column. We then increment one
cursor at every iteration by selecting new non-zero entries until the sparsity
target is reached. We first compute the matrix $\Delta\mathcal{C}$ (Algorithms
2 and 2). An entry $\Delta\mathcal{C}(k,j)$ represents, for the column
$H(:,j)$, the error decrease when going from the $(k-1)$-sparse solution to
the $k$-sparse one. An example of such a matrix is given in Section 4.1. From
$\Delta\mathcal{C}$, we then build a matrix $\mathcal{G}$ (Algorithms 2 and
2), that is essentially $\Delta\mathcal{C}$ with entries divided by the
corresponding $k$. Every entry of $\mathcal{G}$ then represents the error
decrease of a solution, divided by the corresponding cost in sparsity budget.
Finally, at every iteration, shamans selects the maximum entry in
$\mathcal{G}$ (Algorithm 2), guaranteeing to have the largest possible
decrease of the error $\|M-WH\|_{F}^{2}$ for a fixed budget of additional non-
zero entries in $H$ (at each step, there will be between $1$ and $r$ new non-
zero entries in $H$). Note that this operation is very cheap in practice, as
we can implement it by maintaining a sorted list of the maximum entry of each
column, and thus avoid to recompute the maximum over the whole matrix. We use
its coordinates to update the cursor (Algorithm 2). After every iteration, the
column of $\mathcal{G}$ of the last selected solution must be updated
(Algorithms 2 and 2). If the last selected solution is $k$-sparse, then all
the entries corresponding to $k^{\prime}$-sparse solutions with $k^{\prime}<k$
are put to zero in the corresponding column of $\mathcal{G}$, so they will not
be selected during the next iterations. The entries corresponding to
$k^{\prime}$-sparse solutions with $k^{\prime}>k$ are recomputed from
$\Delta\mathcal{C}$ to take into account the new value of the cursor. In other
words, the corresponding column of $\mathcal{G}$ maintains the reduction in
the error $\|M(:,j)-WH(:,j)\|_{2}^{2}$ if we further add non-zero entries in
$H(:,j)$. Finally, on Algorithms 2 and 2, we build the matrix $H^{*}$ with the
final cursors. In Section 4.1, we illustrate shamans with a toy example.

Finally, although motivated by nonnegative problems, our approach is
generalized easily to problems with no nonnegativity constraint. In shamans,
the selection step is not specific to nonnegativity, and the homotopy step for
sparse NNLS can be easily replaced with a standard homotopy algorithm for
sparse LS.

###  4.1 Running example for shamans

Let us consider the following example, with the $5$-by-$6$ data matrix $M$ and
the $5$-by-$4$ dictionary $W$:

$M=\begin{pmatrix}0.89&1.21&0.73&0.8&0.06&0.02\\\
0.65&0.97&1.17&0.23&0.36&0.27\\\ 1.06&1.63&1.27&0.76&0.49&0.15\\\
0.98&1.41&1.32&0.59&0.51&0.2\\\
1.01&1.66&1.57&0.57&0.56&0.29\end{pmatrix},\quad
W=\begin{pmatrix}0.8&0.07&0.1&0.81\\\ 0.07&0.51&0.78&0.4\\\
0.77&0.92&0.4&0.76\\\ 0.47&0.9&0.51&0.7\\\ 0.58&0.9&0.87&0.59\end{pmatrix}.$

We compute $H$ with shamans with a constraint of average column-wise sparsity
$k=3$, that is, a matrix-wise constraint $q=3\times 6=18$. Table 1 shows the
output of the homotopy algorithm, that is, the regularization path containing
the solutions for all possible values of $\lambda$ for each column.

Table 1: Regularization paths corresponding to $\min\limits_{x\geq
0}\|M(:,j)-Wx\|_{2}^{2}+\lambda\|x\|_{1}$, computed with the homotopy
algorithm. Columns $2$ to $5$ are not shown.

Solutions $x^{*}_{t}$ | $\begin{pmatrix}0.0\\\ 0.0\\\ 0.0\\\ 0.0\end{pmatrix}$ | $\begin{pmatrix}0.0\\\ 1.15\\\ 0.0\\\ 0.0\end{pmatrix}$ | $\begin{pmatrix}0.0\\\ 0.34\\\ 0.0\\\ 1.05\end{pmatrix}$ | $\begin{pmatrix}0.0\\\ 0.19\\\ 0.21\\\ 1.04\end{pmatrix}$ | $\begin{pmatrix}0.2\\\ 0.16\\\ 0.28\\\ 0.85\end{pmatrix}$
---|---|---|---|---|---
$\lambda$ | 3.16 | 2.75 | 0.25 | 0.06 | 0.0
Error ${\|M(:,1)-Wx^{*}_{t}\|_{2}^{2}}$ | 4.31 | 0.67 | 0.02 | 0.0 | 0.0
$k$ | 0 | 1 | 2 | 3 | 4
(a) Column $1$

[…]

Solutions $x^{*}_{t}$ | $\begin{pmatrix}0.0\\\ 0.0\\\ 0.0\\\ 0.0\end{pmatrix}$ | $\begin{pmatrix}0.0\\\ 0.26\\\ 0.0\\\ 0.0\end{pmatrix}$ | $\begin{pmatrix}0.0\\\ 0.03\\\ 0.31\\\ 0.0\end{pmatrix}$
---|---|---|---
$\lambda$ | 0.71 | 0.37 | 0.0
Error $\|M(:,6)-Wx^{*}_{t}\|_{2}^{2}$ | 0.21 | 0.03 | 0.0
$k$ | 0 | 1 | 2
(b) Column $6$

From these paths, we can trivially get the solution to the unpenalized problem
$\min\limits_{H\geq 0}\|M-WH\|_{F}^{2}$ by selecting the last solution for
each column. We can also trivially select column-wise the $k$-sparse solution.
To select the best solution with respect to a matrix-wise sparsity constraint,
shamans needs extra steps. Tables 2 and 3 show respectively the corresponding
matrices $\mathcal{C}$, and $\Delta\mathcal{C}$. Table 4 shows the
corresponding matrix $\mathcal{G}$ at different steps of the selection
process.

Table 2: Matrix $\mathcal{C}$, where every entry $\mathcal{C}(k,j)$ represents the error of the $k$-sparse solution of the column $j$, $\mathcal{C}(k,j)=\|M(:,j)-WSol_{k,j}\|_{2}^{2}$. | $j$ | 1 | 2 | 3 | 4 | 5 | $n=6$
---|---|---|---|---|---|---|---
$k$ |  |  |  |  |  |  |
0 |  | 4.31 | 9.81 | 7.70 | 1.95 | 0.94 | 0.21
1 |  | 0.67 | 1.24 | 0.59 | 0.03 | 0.01 | 0.03
2 |  | 0.02 | 0.09 | 0.24 | 0.0 | 0.0 | 0.0
3 |  | 0.0 | 0.05 | 0.0 | 0.0 | 0.0 | 0.0
$r=4$ |  | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0
Table 3: Matrix $\Delta\mathcal{C}$, where every entry $\Delta\mathcal{C}(k,j)$ represents the error decrease by going in column $j$ from the $(k-1)$-sparse solution to the $k$-sparse one, $\Delta\mathcal{C}(k,j)=\mathcal{C}(k-1,j)-\mathcal{C}(k,j)$. | $j$ | 1 | 2 | 3 | 4 | 5 | $n=6$
---|---|---|---|---|---|---|---
$k$ |  |  |  |  |  |  |
1 |  | 3.64 | 8.57 | 7.11 | 1.92 | 0.93 | 0.19
2 |  | 0.65 | 1.14 | 0.35 | 0.03 | 0.01 | 0.03
3 |  | 0.01 | 0.05 | 0.24 | 0.0 | 0.0 | 0.0
$r=4$ |  | 0.0 | 0.05 | 0.0 | 0.0 | 0.0 | 0.0
Table 4: Matrix $\mathcal{G}$ during the first three iterations of the
selection step of shamans. Every entry represents the error decrease of a
solution divided by its cost in sparsity budget (the number by which we need
to increment the cursor of its column in order to select this solution). Every
entry is initialized so that
$\mathcal{G}(k,j)=(\sum_{l=1}^{k}\Delta\mathcal{C}(l,j))/i$. When an entry is
selected, we update the other entries of this column to account for the new
error decrease and cursor value. Below, the boxed entry is the maximum of the
matrix, so it is the entry to be selected next.

| $j$ | 1 | 2 | 3 | 4 | 5 | $n=6$
---|---|---|---|---|---|---|---
$k$ |  |  |  |  |  |  |
1 |  | 3.64 | $None$ | 7.11 | 1.92 | 0.93 | 0.19
2 |  | 2.15 | 4.86 | 3.73 | 0.97 | 0.47 | 0.11
3 |  | 1.44 | 3.25 | 2.57 | 0.65 | 0.31 | 0.07
$r=4$ |  | 1.08 | 2.45 | 1.93 | 0.49 | 0.24 | 0.05
(a) Initialization

| $j$ | 1 | 2 | 3 | 4 | 5 | $n=6$
---|---|---|---|---|---|---|---
$k$ |  |  |  |  |  |  |
1 |  | 3.64 | 0.0 | $None$ | 1.92 | 0.93 | 0.19
2 |  | 2.15 | 1.14 | 3.73 | 0.97 | 0.47 | 0.11
3 |  | 1.44 | 0.6 | 2.57 | 0.65 | 0.31 | 0.07
$r=4$ |  | 1.08 | 0.41 | 1.93 | 0.49 | 0.24 | 0.05
(b) After one iteration

| $j$ | 1 | 2 | 3 | 4 | 5 | $n=6$
---|---|---|---|---|---|---|---
$k$ |  |  |  |  |  |  |
1 |  | $None$ | 0.0 | 0.0 | 1.92 | 0.93 | 0.19
2 |  | 2.15 | 1.14 | 0.35 | 0.97 | 0.47 | 0.11
3 |  | 1.44 | 0.6 | 0.29 | 0.65 | 0.31 | 0.07
$r=4$ |  | 1.08 | 0.41 | 0.20 | 0.49 | 0.24 | 0.05
(c) After two iterations

With a $3$-sparse column-wise selection, the result is $H_{ks}$ shown in (12),
for a relative reconstruction error $\frac{\|M-WH\|_{F}}{\|M\|_{F}}=4.50\%$.
With shamans with $q=18$, that is, with average 3-sparse columns is $H_{Sh}$
shown in (12), and has a relative reconstruction error of $0.73\%$:

$H_{ks}=\begin{pmatrix}0.0&0.0&0.0&0.49&0.0&0.0\\\
0.19&0.42&0.38&0.0&0.45&0.03\\\ 0.21&0.38&0.88&0.0&0.16&0.31\\\
1.04&1.38&0.75&0.49&0.02&0.0\end{pmatrix},H_{Sh}=\begin{pmatrix}0.2&0.79&0.06&0.49&0.0&0.0\\\
0.16&0.32&0.37&0.0&0.46&0.03\\\ 0.28&0.66&0.9&0.0&0.16&0.31\\\
0.85&0.61&0.69&0.5&0.0&0.0\end{pmatrix}.$ (12)

###  4.2 Optimality and computational cost of shamans

In the following two paragraphs, we discuss the optimality and the
computational cost of shamans.

#### Optimality

Interestingly, although shamans is a greedy algorithm as it performs an
optimal selection at each iteration, it is globally optimal. The reason is
that the columns of $H$ do not interact with each other in the objective
function (see Section 1.2), so the greedy solution cannot possibly be
improved.

###### Proposition 1

The selection step of shamans is guaranteed to select the best solution $H$
with average $(q+p)/n$-sparse columns among the set of solutions returned by
the homotopy algorithm, where $p\in[0,r-1]$, meaning that $\|H\|_{0}\leq q+p$.

Indeed, shamans selects the best solution in terms of reducing the error
$\|M-WH\|_{F}^{2}$ for a fixed sparsity budget, since the largest entry of
$\mathcal{G}$ is selected, while the columns of $H$ do not interact with each
other in the objective. Note that shamans adds non-zero entries in $H$ as long
as $\|H\|_{0}<q$, and it adds at each iteration at most $r$ non-zero entries
in a column (a column could change from the zero vector to the unconstrained
NNLS solution with potentially $r$ non-zero entries). In this particular case,
the solution produced by shamans could reach a sparsity of $\|H\|_{0}=q+r-1$.
However, this is not an important drawback for three reasons:

  * •

Only the last updated column can lead to the $q$-sparsity constraint to be
violated. Moreover, since $n\gg r$, we have $(q+r)/n\approx q/n$ and this does
not play an important role in practice.

  * •

If the constraint $\|H\|_{0}\leq q$ is critical, it suffices to keep the
penultimate solution generated by shamans.

  * •

In most cases, we have observed that each iteration of shamans only increases
the number of non-zero entries in $H$ by one; see for example Table 4.
However, there are some pathological cases where shamans may increase that
number by more than one.

Note that the homotopy step is approximate, so shamans does not solve Problem
(7) optimally in general. What shamans does optimally is to select a
$q$-sparse solution $H$ from a set of solutions provided (in our case, by the
homotopy algorithm).

#### Computational cost

The selection step consists in building $\Delta\mathcal{C}$ in
$\mathcal{O}(rn)$, building $\mathcal{G}$ in $\mathcal{O}(rn)$, then iterating
$q$ times to select a solution. As we maintain updated a sorted list to avoid
recomputing the maximum at each iteration, this is done in
$\mathcal{O}(q\log(n))$. Assuming $q$ is in the order of $rn$, the cost of the
selection step is hence dominated by the cost of the homotopy step, that
requires at least $O(r^{3}n)$ operations to solve $n$ NNLS problems via
homotopy in $r$ variables; see Section 3.3.

##  5 Experiments

In this section, we study the performance of shamans on the unmixing of 7
datasets: 3 faces datasets and 4 hyperspectral images.

###  5.1 Data

For the faces datasets, each column of $M$ corresponds to a pixel and each row
to an image (that is, $M(i,j)$ is the intensity of pixel $j$ in image $i$). It
is well-known that NMF will extract facial features as the rows of matrix $H$
(Lee and Seung, 1999). As no groundtruth is available, we first compute $W$
with SNPA (Gillis, 2014a), an algorithm for separable NMF, setting the
factorization rank $r$ as in the literature. We then compute $H$ with our
sparsity-enhancing method. Imposing sparsity on $H$ means that we require that
only a few pixels are contained in each facial feature (Hoyer, 2004). We
consider the 3 widely used face datasets CBCL333Downloaded from http://poggio-
lab.mit.edu/codedatasets, Frey444Downloaded from
https://cs.nyu.edu/~roweis/data.html, and Kuls555Downloaded from
http://www.robots.ox.ac.uk/.

Similarly as for the faces, a hyperspectral image is an image-by-pixel matrix
where each image corresponds to a different wavelength. The columns of $W$
represent the spectral signature of the pure materials (also called
endmembers) present in the image (Bioucas-Dias et al, 2012), and we use the
ground truth $W$ of Zhu (2017). We compute $H$, whose columns represent the
abundance of materials in each pixel. It makes sense to impose $H$ to be
sparse as most pixels contain only a few endmembers (Ma et al, 2013). We
consider the 4 widely used datasets666Downloaded from
http://lesun.weebly.com/hyperspectral-data-set.html, Jasper, Samson, Urban,
and Cuprite. The characteristics of the considered datasets are summarized in
Table 5.

Table 5: Summary of the datasets studied in this work. $M\in\mathbb{R}^{m\times n}$ and $W\in\mathbb{R}^{m\times r}$. Dataset | Type | m | n | r
---|---|---|---|---
CBCL | Faces | 2429 | $19\times 19=361$ | 49
Frey | Faces | 1965 | $20\times 28=560$ | 36
Kuls | Faces | 20 | $64\times 64=4096$ | 5
Jasper | Hyperspectral | 198 | $100\times 100=1000$ | 4
Samson | Hyperspectral | 156 | $95\times 95=9025$ | 3
Urban | Hyperspectral | 162 | $307\times 307=94249$ | 6
Cuprite | Hyperspectral | 188 | $250\times 191=47750$ | 12

###  5.2 Methods

We compare the performances of shamans with three other approaches:

  * •

Homotopy with no sparsity constraint, noted _H no-s_. It consists in computing
the homotopy path for each column and keep the best solution.

  * •

Homotopy with a column-wise $k$-sparse constraint, noted _H $k$-s_. It
consists in computing the homotopy path for each column and keep the
corresponding $k$-sparse solution.

  * •

NNOLS (Nguyen et al, 2019). This greedy algorithm is one of the best
performing heuristic for $k$-sparse NNLS.

To the best of our knowledge, no other algorithm tackles the problem of MNNLS
with a matrix-wise $q$-sparsity constraint.

We choose the parameter $k$ by trial-and-error. Unless stated otherwise, we
define the sparsity parameter of shamans as $q=k\times n$, which is equivalent
to an average column-wise $k$-sparsity constraint.

The homotopy algorithm and shamans are implemented in Julia, and ran with
Julia 1.5. For NNOLS, we use the Matlab implementation provided by the
authors, and run it in Matlab R2019b. All experiments are made on a computer
with a processor Intel Core i5-2520M @2.50GHz. The code and experiments
scripts are provided in an online
repository777https://gitlab.com/nnadisic/giant.jl.

For every dataset, we run the four methods, and measure the average sparsity
of the given solutions (defined as number of entries superior to $10^{-3}$
divided by the number of columns), the relative reconstruction error
$\frac{\|M-WH\|_{F}}{\|M\|_{F}}$, and the computing time, for which we measure
the median over 10 runs. Note that, for a given dataset, a given algorithm
always gives the same output.

###  5.3 Results and interpretation

The results of our experiments are shown in Table 6. Some abundance maps for
Kuls, Jasper, and Urban are shown respectively in Figs. 3, 4 and 5. Images
from the other experiments are available as a supplementary material.

Table 6: Results of the experiments, for the unmixing of faces and hyperspectral datasets. Time is in seconds, relative error is in percent, and sparsity is the average number of entries above $10^{-3}$. Jasper is processed once with all algorithms for $k=q/n=2$, and once with shamans for $q/n=1.8$ (which is not possible with other algorithms). |  | H no-s | H $k$-s | shamans | NNOLS
---|---|---|---|---|---
CBCL | Time | 0.78 | 0.77 | 0.84 | 0.26
$r=49$ | Rel error | 12.04 | 16.19 | 13.22 | 12.60
$k=3$ | Sparsity | 6.53 | 2.69 | 3.0 | 2.71
Frey | Time | 1.31 | 1.32 | 1.46 | 0.86
$r=36$ | Rel error | 19.35 | 23.22 | 20.75 | 20.30
$k=6$ | Sparsity | 12.17 | 5.51 | 6.0 | 5.53
Kuls | Time | 0.25 | 0.21 | 0.30 | 1.87
$r=5$ | Rel error | 19.05 | 20.13 | 19.12 | 19.20
$k=3$ | Sparsity | 3.44 | 2.86 | 3.0 | 2.84
Jasper | Time | 0.63 | 0.63 | 0.74 | 2.25
$r=4$ | Rel error | 5.71 | 6.99 | 5.72 | 6.09
$k=2$ | Sparsity | 2.27 | 1.79 | 2.0 | 1.81
Jasper | Time |  |  | 0.65 |
$r=4$ | Rel error |  |  | 5.95 |
$q/n=1.8$ | Sparsity |  |  | 1.8 |
Samson | Time | 0.33 | 0.33 | 0.45 | 1.96
$r=3$ | Rel error | 3.3 | 3.34 | 3.30 | 3.34
$k=2$ | Sparsity | 2.14 | 1.84 | 2.0 | 1.83
Urban | Time | 6.67 | 6.55 | 10.73 | 22.04
$r=6$ | Rel error | 7.67 | 8.62 | 7.83 | 8.14
$k=2$ | Sparsity | 2.61 | 1.90 | 2.0 | 1.90
Cuprite | Time | 17.83 | 20.27 | 22.31 | 36.30
$r=12$ | Rel error | 1.74 | 2.37 | 2.01 | 2.04
$k=4$ | Sparsity | 6.60 | 3.93 | 4.0 | 3.93

(a) Homotopy with no sparsity constraint

(b) Homotopy with $k$-sparse selection

(c) shamans

(d) NNOLS

Figure 3: Abundance maps of the unmixing of the Kuls facial images (that is,
reshaped rows of $W$) with different algorithms. The extracted components
correspond to different light directions.

(a) H no-s

(b) H $k$-s

(c) shamans

(d) NNOLS

Figure 4: Abundance maps of the sixth endmember (that is, the sixth row of $H$
reshaped) from the unmixing of the hyperspectral image Urban by different
algorithms.

(a) H no-s

(b) shamans with $k=q/n=2$

(c) shamans with $k=q/n=1.8$

(d) NNOLS

Figure 5: Abundance maps of the second endmember (that is, the second row of
$H$ reshaped) from the unmixing of the hyperspectral image Jasper by different
algorithms.

The results show that, without sparsity constraint, the resulting factors are
already quite sparse. The column-wise sparse methods (H $k$-s and NNOLS)
produce factors with an average sparsity below the sparsity threshold $k$,
meaning that many columns actually have a sparsity level lower than the
threshold. H $k$-s gives relatively bad results in terms of error, comparing
with the unconstrained solution, while shamans is able to enforce sparsity
while limiting the loss in error. H $k$-s is as fast as the unconstrained
method, as expected. shamans is only slightly slower, meaning that the
selection step takes less time than the homotopy step, as expected in Section
4.2. These methods are faster than NNOLS for all datasets expect CBCL and
Frey.

NNOLS outperforms our approach in terms of error for CBCL and Frey. This means
that the column-wise sparsity is a good assumption for these datasets, and
also that our homotopy step does not give good enough solutions. However,
shamans outperforms NNOLS for all the other datasets. This verifies our
assumption that columns have different sparsity levels, a situation in which
our approach is especially interesting. Another advantage of the matrix-wise
constraint over the column-wise one is the possibility of tuning the parameter
more precisely. This is illustrated with the result for Jasper with $q/n=1.8$
(that is, we require an average $k$-sparsity with $k=1.8$), for which we reach
an average sparsity similar to that of NNOLS with $k=2$, and still get a lower
error.

Visually, in Fig. 3, the unconstrained method already gives good results, but
some features are not well separated. With the $k$-sparse selection, it fails
to retrieve the features and sometimes produces noisier images. The results of
NNOLS are better, but lack spatial coherence (some regions are pixelated; see
for example the chin on the last image in Figure 3). shamans produces better
separated features, and a more consistent spatial coherence.

In Fig. 4, we focus on the sixth endmember for the Urban dataset. It mainly
shows pixels from rooftops. With no sparsity constraint, the image is very
noisy and includes pixels composed of other materials. Visually, NNOLS gives
almost no improvement, and worsen the spatial coherence (some regions are
pixelated). H $k$-s removes too many pixels and most of the information for
that endmember is lost. shamans produces a better solution, with less noise,
and a more distinct separation.

In Fig. 5, we focus on the second endmember for the Jasper dataset,
corresponding to water. Without sparsity constraint, they are mixed with
pixels of the road on the right. NNOLS with $k=2$ and shamans with $k=q/n=2$
fail to separate these pixels. shamans with $k=q/n=1.8$ performs an almost
perfect separation, and eliminates the non-water pixels while sharpening the
edge between water and non-water pixels. In these images, it is clear that
some pixels are composed of only one endmember (for instance, water) while
others are mixed (trees and grass mixed). This is the kind of setting where
shamans can bring significant improvement.

##  6 Conclusion

In this paper, we have focused on the sparse multiple right-hand sides NNLS
problem with a matrix-wise $\ell_{0}$-constraint. We introduced shamans
(Algorithm 2), a two-step algorithm to tackle this problem, that first
computes for each column the regularization path of the $\ell_{1}$-penalized
subproblem via a homotopy algorithm, and then applies a selection strategy to
build the solution matrix $H$ from these paths. Our method solves exactly the
selection subproblem. We illustrated the advantages of shamans for the
unmixing of facial and hyperspectral images, for which it generally
outperformed the state-of-the-art heuristic NNOLS in both reconstruction error
and computing time. shamans is particularly useful for datasets for which the
columns of $H$ have different sparsity levels. Also, it is almost as fast as a
standard active-set algorithms for NNLS with no sparsity constraint, and thus
can be applied to large datasets.

###### Acknowledgements.

We thank Maxime De Wolf for his technical contribution to the implementation
of our methods. We thank T.T. Nguyen and his co-authors for making the codes
of nonnegative greedy algorithms available online under a free software
license.

## References

  * Aharon et al (2005) Aharon M, Elad M, Bruckstein AM (2005) K-SVD and its non-negative variant for dictionary design. In: Wavelets XI, International Society for Optics and Photonics
  * Bertsimas and Shioda (2009) Bertsimas D, Shioda R (2009) Algorithm for cardinality-constrained quadratic optimization. Computational Optimization and Applications 43(1):1–22
  * Bienstock (1996) Bienstock D (1996) Computational study of a family of mixed-integer quadratic programming problems. Mathematical programming 74(2):121–140
  * Bioucas-Dias et al (2012) Bioucas-Dias JM, Plaza A, Dobigeon N, Parente M, Du Q, Gader P, Chanussot J (2012) Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 5(2):354–379
  * Bourguignon et al (2015) Bourguignon S, Ninin J, Carfantan H, Mongeau M (2015) Exact sparse approximation problems via mixed-integer programming: Formulations and computational performance. IEEE Transactions on Signal Processing 64(6):1405–1419
  * Chen et al (1989) Chen S, Billings SA, Luo W (1989) Orthogonal least squares methods and their application to non-linear system identification. International Journal of control 50(5):1873–1896
  * Cichocki et al (2008) Cichocki A, Phan AH, Caiafa C (2008) Flexible HALS algorithms for sparse non-negative matrix/tensor factorization. In: IEEE Workshop on Machine Learning for Signal Processing, pp 73–78
  * Cohen and Gillis (2019) Cohen JE, Gillis N (2019) Nonnegative Low-rank Sparse Component Analysis. In: 2019 IEEE International Conference on Acoustics, Speech and Signal Processing, pp 8226–8230
  * Donoho and Tsaig (2008) Donoho DL, Tsaig Y (2008) Fast Solution of $\ell_{1}$-Norm Minimization Problems When the Solution May Be Sparse. IEEE Trans Inform Theory 54(11):4789–4812
  * Efron et al (2004) Efron B, Hastie T, Johnstone I, Tibshirani R, et al (2004) Least angle regression. The Annals of statistics 32(2):407–499
  * Eggert and Korner (2004) Eggert J, Korner E (2004) Sparse coding and NMF. In: IEEE International Joint Conference on Neural Networks, vol 4, pp 2529–2533
  * Foucart and Koslicki (2014) Foucart S, Koslicki D (2014) Sparse recovery by means of nonnegative least squares. IEEE Signal Processing Letters 21(4):498–502
  * Gillis (2012) Gillis N (2012) Sparse and unique nonnegative matrix factorization through data preprocessing. Journal of Machine Learning Research 13:3349–3386
  * Gillis (2014a) Gillis N (2014a) Successive Nonnegative Projection Algorithm for Robust Nonnegative Blind Source Separation. SIAM Journal on Imaging Sciences 7(2):1420–1450
  * Gillis (2014b) Gillis N (2014b) The why and how of nonnegative matrix factorization. Regularization, Optimization, Kernels, and Support Vector Machines 12(257):257–291
  * Hoyer (2002) Hoyer PO (2002) Non-negative sparse coding. In: IEEE Workshop On Neural Networks for Signal Processing, pp 557–565
  * Hoyer (2004) Hoyer PO (2004) Non-negative matrix factorization with sparseness constraints. Journal of machine learning research 5:1457–1469
  * Itoh et al (2017) Itoh Y, Duarte MF, Parente M (2017) Perfect Recovery Conditions For Non-Negative Sparse Modeling. IEEE Transactions on Signal Processing 65(1):69–80
  * Kim and Park (2007) Kim H, Park H (2007) Sparse non-negative matrix factorizations via alternating non-negativity-constrained least squares for microarray data analysis. Bioinformatics 23(12):1495–1502
  * Kim et al (2013) Kim J, Ramakrishnan N, Marwah M, Shah A, Park H (2013) Regularization Paths for Sparse Nonnegative Least Squares Problems with Applications to Life Cycle Assessment Tree Discovery. In: IEEE 13th International Conference on Data Mining, pp 360–369
  * Lawson and Hanson (1995) Lawson CL, Hanson RJ (1995) Solving least squares problems. Society for Industrial and Applied Mathematics
  * Lee and Seung (1997) Lee DD, Seung HS (1997) Unsupervised learning by convex and conic coding. In: Advances in neural information processing systems, pp 515–521
  * Lee and Seung (1999) Lee DD, Seung HS (1999) Learning the parts of objects by non-negative matrix factorization. Nature 401(6755):788–791
  * Ma et al (2013) Ma WK, Bioucas-Dias JM, Chan TH, Gillis N, Gader P, Plaza AJ, Ambikapathi A, Chi CY (2013) A signal processing perspective on hyperspectral unmixing: Insights from remote sensing. IEEE Signal Processing Magazine 31(1):67–81
  * Mhenni et al (2019) Mhenni R, Bourguignon S, Ninin J (2019) Global optimization for sparse solution of least squares problems. https://hal.archives-ouvertes.fr/hal-02066368/
  * Morup et al (2008) Morup M, Madsen KH, Hansen LK (2008) Approximate L0 constrained non-negative matrix and tensor factorization. In: 2008 IEEE International Symposium on Circuits and Systems, IEEE, pp 1328–1331
  * Nadisic et al (2020) Nadisic N, Vandaele A, Gillis N, Cohen JE (2020) Exact Sparse Nonnegative Least Squares. In: ICASSP 2020, pp 5395 – 5399
  * Nguyen et al (2019) Nguyen TT, Idier J, Soussen C, Djermoune EH (2019) Non-Negative Orthogonal Greedy Algorithms. IEEE Transactions on Signal Processing pp 1–16
  * Osborne et al (2000) Osborne MR, Presnell B, Turlach BA (2000) A new approach to variable selection in least squares problems. IMA journal of numerical analysis 20(3):389–403
  * Pati et al (1993) Pati YC, Rezaiifar R, Krishnaprasad PS (1993) Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition. In: Proceedings of 27th Asilomar Conference on Signals, Systems and Computers, pp 40–44
  * Peharz and Pernkopf (2012) Peharz R, Pernkopf F (2012) Sparse nonnegative matrix factorization with l0-constraints. Neurocomputing 80:38–46
  * Soussen et al (2013) Soussen C, Gribonval R, Idier J, Herzet C (2013) Joint $k$-Step Analysis of Orthogonal Matching Pursuit and Orthogonal Least Squares. IEEE Transactions on Information Theory 59(5):3158–3174
  * Tibshirani (1996) Tibshirani R (1996) Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological) 58(1):267–288
  * Tropp (2006) Tropp JA (2006) Just relax: convex programming methods for identifying sparse signals in noise. IEEE Transactions on Information Theory 52(3):1030–1051
  * Zhu (2017) Zhu F (2017) Hyperspectral unmixing: ground truth labeling, datasets, benchmark performances and survey. arXiv preprint arXiv:170805125

## Supplementary material for the paper “A Homotopy-based Algorithm for Sparse
Multiple Right-hand Sides Nonnegative Least Squares”

In this document, we provide the abundance maps resulting from our
experiments, that could not be included in the paper:

  * •

CBCL in Figure 6,

  * •

Frey in Figure 7,

  * •

Jasper in Figure 8,

  * •

Samson in Figure 9,

  * •

Urban in Figure 10, and

  * •

Cuprite in Figure 11.

As explained in the paper, the abundance maps generated by shamans tends to be
less noisy and less pixelated compare to the other sparse MNNLS algorithms.

(a) Homotopy with no sparsity constraint

(b) Homotopy with $k$-sparse selection

(c) shamans

(d) NNOLS

Figure 6: Abundance maps of the unmixing of the CBCL facial images (that is,
reshaped rows of $W$) with different algorithms.

(a) Homotopy with no sparsity constraint

(b) Homotopy with $k$-sparse selection

(c) shamans

(d) NNOLS

Figure 7: Abundance maps of the unmixing of the Frey facial images (that is,
reshaped rows of $W$) with different algorithms.

(a) Homotopy with no sparsity constraint

(b) Homotopy with $k$-sparse selection

(c) shamans with $q=20000$

(d) shamans with $q=18000$

(e) NNOLS

Figure 8: Abundance maps (that is,reshaped rows of $H$) from the unmixing of
the hyperspectral image Jasper by different algorithms.

(a) Homotopy with no sparsity constraint

(b) Homotopy with $k$-sparse selection

(c) shamans

(d) NNOLS

Figure 9: Abundance maps (that is, reshaped rows of $H$) from the unmixing of
the hyperspectral image Samson by different algorithms.

(a) Homotopy with no sparsity constraint

(b) Homotopy with $k$-sparse selection

(c) shamans

(d) NNOLS

Figure 10: Abundance maps (that is, reshaped rows of $H$) from the unmixing of
the hyperspectral image Urban by different algorithms.

(a) Homotopy with no sparsity constraint

(b) Homotopy with $k$-sparse selection

(c) shamans

(d) NNOLS

Figure 11: Abundance maps (that is, reshaped rows of $H$) from the unmixing of
the hyperspectral image Cuprite by different algorithms.