# Multiplicative Updates for NMF with $\beta$-Divergences under Disjoint Equality Constraints

Valentin Leplat, Nicolas Gillis, Jérôme Idier

► To cite this version:

Valentin Leplat, Nicolas Gillis, Jérôme Idier. Multiplicative Updates for NMF with $\beta$-Divergences under Disjoint Equality Constraints. SIAM Journal on Matrix Analysis and Applications, 2021, 42 (2), pp.730-752. 10.1137/20M1377278. hal-02993314

HAL Id: hal-02993314

https://hal.science/hal-02993314v1

Submitted on 7 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multiplicative Updates for NMF with $\beta$-Divergences under Disjoint Equality Constraints

Valentin Leplat*

Nicolas Gillis*

Jérôme Idier†

April 29, 2022

## Abstract

Nonnegative matrix factorization (NMF) is the problem of approximating an input nonnegative matrix, $V$, as the product of two smaller nonnegative matrices, $W$ and $H$. In this paper, we introduce a general framework to design multiplicative updates (MU) for NMF based on $\beta$-divergences ($\beta$-NMF) with disjoint equality constraints, and with penalty terms in the objective function. By disjoint, we mean that each variable appears in at most one equality constraint. Our MU satisfy the set of constraints after each update of the variables during the optimization process, while guaranteeing that the objective function decreases monotonically. We showcase this framework on three NMF models, and show that it competes favorably the state of the art: (1) $\beta$-NMF with sum-to-one constraints on the columns of $H$, (2) minimum-volume $\beta$-NMF with sum-to-one constraints on the columns of $W$, and (3) sparse $\beta$-NMF with $\ell_2$-norm constraints on the columns of $W$.

**Keywords:** nonnegative matrix factorization (NMF), $\beta$-divergences, disjoint constraints, simplex-structured NMF, minimum-volume NMF, sparsity

## 1 Introduction

Given a non-negative matrix $V \in \mathbb{R}_+^{F \times N}$ and a factorization rank $K \ll \min(F, N)$, nonnegative matrix factorization (NMF) aims to compute two non-negative matrices, $W$ with $K$ columns and $H$ with $K$ rows, such that $V \approx WH$ [24]. Over the last two decades, NMF has shown to be a powerful tool for the analysis of high-dimensional data. The main reason is that NMF automatically extracts sparse and meaningful features from a set of nonnegative data vectors. NMF has been successfully used in many applications such as image processing, text mining, hyperspectral imaging, blind source separation, single-channel audio source separation, clustering and music analysis; see [15, 7, 5, 29, 13, 16] and the references therein.

To compute $W$ and $H$, the most standard approach is to solve the following optimization problem

$$ \min_{W \in \mathbb{R}^{F \times K}, H \in \mathbb{R}^{K \times N}} D(V|WH) \quad \text{such that} \quad H \ge 0 \text{ and } W \ge 0, \quad (1) $$

<footnote>*Department of Mathematics and Operational Research, Faculté Polytechnique, Université de Mons, Rue de Houdain 9, 7000 Mons, Belgium. Authors acknowledge the support by the Fonds de la Recherche Scientifique - FNRS and the Fonds Wetenschappelijk Onderzoek - Vlanderen (FWO) under EOS Project no O005318F-RG47, and by the European Research Council (ERC Starting Grant no 679515, ERC Consolidator Grant no 681839). E-mails: {valentin.leplat, nicolas.gillis}@umons.ac.be.</footnote>

<footnote>†Laboratoire des Sciences du Numérique de Nantes (LS2N, CNRS UMR 6004), Ecole Centrale de Nantes, 44321 Nantes, France. E-mail: jerome.idier@ls2n.fr</footnote>
where $D(V|WH) = \sum_{f,n} d(V_{fn}|[WH]_{fn})$ with $d(x|y)$ a measure of distance between two scalars, and $A \ge 0$ means that the matrix $A$ is component-wise nonnegative. In this paper, we focus on $\beta$-NMF for which the measure of fit is the discrete $\beta$-divergence denoted $d_{\beta}(x|y)$ and defined as

$$d_{\beta}(x|y) = \begin{cases} \frac{1}{\beta(\beta-1)}(x^{\beta} + (\beta-1)y^{\beta} - \beta xy^{\beta-1}) & \text{for } \beta \in \mathbb{R} \setminus \{0, 1\}, \\ x \log \frac{x}{y} - x + y & \text{for } \beta = 1, \\ \frac{x}{y} - \log \frac{x}{y} - 1 & \text{for } \beta = 0. \end{cases}$$

For $\beta = 2$, $d_2(x|y) = \frac{1}{2}(x - y)^2$, so $D(V|WH)$ is the halved standard squared Euclidean distance between $V$ and $WH$, that is, the halved squared Frobenius norm $\frac{1}{2}\|V - WH\|_F^2$. For $\beta = 1$ and $\beta = 0$, the $\beta$-divergence corresponds to the Kullback-Leibler (KL) divergence and the Itakura-Saito (IS) divergence, respectively. The error measure should be chosen accordingly with the distribution of the noise assumed on the data. The Frobenius norm assumes i.i.d. Gaussian noise, KL divergence assumes a Poisson distribution, and the IS divergence assumes multiplicative gamma noise; see for example [11, 8, 20] and the references therein. In the NMF literature, $\beta$-divergences are the most widely used objective functions.

Most NMF algorithms developed to tackle (1) are based on iterative schemes that alternatively updates the factors $W$ and $H$. At each iteration, the minimization over one factor, $W$ or $H$, is performed with various optimization methods. For $\beta$-divergences, the most popular approach is to use multiplicative updates (MU) which were introduced for NMF in the seminal papers of Lee and Seung [24, 23]. In all applications we are aware of, $\beta$ is always chosen smaller than two. The reason is that, for $\beta > 2$, $\beta$-divergences become more and more sensitive to outliers. Already for $\beta = 2$, it is well-known that the squared Frobenius norm is sensitive to outliers. However, the case $\beta = 2$ is particular because the subproblem in $W$ and $H$ are nonnegative least squares problems, that is, convex quadratic problems with Lipschitz continuous gradient. Therefore, highly efficient schemes exist when $\beta = 2$ that outperform the MU; for example exact block coordinate descent methods [6, 17, 22, 2], or fast gradient methods [18, 19]. In this paper, we focus on the case $\beta < 2$.

In many applications, on top of the nonnegative constraints on the variables, additional constraints are needed to provide a meaningful solution. An instrumental example is the constraint that the entries in each column of $H$ sum to one; this is the so-called sum-to-one constraint that is crucial in blind hyperspectral unmixing; see Section 3. Another example is a sum-to-one constraint on the columns of $W$ along with a volume regularizer on $W$. This model leads to identifiability of the factors $W$ and $H$ under mild conditions; see Section 4. Most algorithms that deal with such equality constraints do it a posteriori with a projection onto the feasible set, or with a renormalization of the columns of $W$ and the rows of $H$ (that is, replace $W(:,k)$ and $H(k,:)$ by $\alpha_k W(:,k)$ and $H(k,:)/\alpha_k$ for some $\alpha_k > 0$), so that their product $WH$ remains unchanged, and hence $D(V|WH)$ remains unchanged. Such approaches are not ideal:

* Projection requires to perform a line-search to ensure the monotonicity of the algorithm, that is, to ensure that the objective does not increase after each iteration, which may be computationally heavy.

* Renormalization of the columns of $W$ and the rows of $H$ is only useful when each constraint applies to a column of $W$ or a row of $H$. It is not applicable for example for the sum-to-one constraint on the columns of $H$ mentioned above. Moreover, in the presence of regularization terms in the objective function, it may destroy the monotonicity of the algorithm.
Another approach is to use parametrization. However, as far as we know, it does not guarantee the monotonicity of the algorithm; see Section 3 for more details.

**Outline and contribution** In this paper, we introduce a general framework to design MU for $\beta$-NMF with disjoint linear equality constraints, and with penalty terms in the objective function. By disjoint, we mean that each variable appears in at most one equality constraint. This framework, presented in Section 2, does not resort to projection, renormalization, or parametrization. Our MU satisfy the set of constraints after each update of the variables during the optimization process, while guaranteeing that the objective function decreases monotonically. This framework works as follows:

* First, as for the standard MU for $\beta$-NMF, we majorize the objective function using a separable majorizer, that is, the majorizer is the sum of functions involving a single variable.

* Second, we construct the augmented Lagrangian for the majorizer. Because the majorizer is separable, the problem can be decomposed into independent subproblems involving only variables that occur in the same equality constraint since they are disjoint. For a fixed value of the Lagrange multipliers, we prove that the solution of these subproblems are unique, under mild conditions (Proposition 1). Moreover, it can be written in closed form via MU for specific values of $\beta$ and depending on the regularizer used (this is summarized in Table 1).

* Finally, we prove that, under mild conditions, there is a unique solution for the Lagrange multipliers so that the equality constraints are satisfied (Proposition 2). This allows us to apply the Newton-Raphson method to compute the Lagrange multipliers while guaranteeing quadratic convergence (Proposition 3).

We then showcase this framework on two NMF models, and show that it competes favorably with the state of the art:

1. A $\beta$-NMF model with sum-to-one constraints on the columns of $H$, which we refer to as simplex-structured $\beta$-NMF (Section 3), and

2. A minimum-volume $\beta$-NMF model with sum-to-one constraints on the columns of $W$ (Section 4).

Finally, Section 5 shows that the framework can be extended to the case of quadratic disjoints constraints, which we showcase on sparse $\beta$-NMF with $\ell_2$-norm constraints on the columns of $W$.

## 2 General framework to design MU for $\beta$-NMF under disjoint linear equality constraints and penalization

In this paper, we introduce a general framework to tackle $\beta$-NMF with disjoint linear equality constraints, and with penalty terms in the objective function. Let us first introduce specific notations: given a matrix $A \in \mathbb{R}^{F \times N}$ and a list of indices $K \subseteq \{(f, n) \mid 1 \le f \le F, 1 \le n \le N\}$, we denote by $A(K)$ the vector of dimension $|K|$ whose entries are the entries of $A$ corresponding to the indices within $K$. Let us introduce $K_i (1 \le i \le I)$ and $B_j (1 \le j \le J)$ to be disjoint sets of indices for the entries of $W$ and $H$, respectively, that is,

* $K_i \subseteq \{(f, k) \mid 1 \le f \le F, 1 \le k \le K\}$ for $i = 1, 2, \dots, I$,

* $B_j \subseteq \{(k, n) \mid 1 \le k \le K, 1 \le n \le N\}$ for $j = 1, 2, \dots, J$,
• $K_u \cap K_v = \emptyset$ for all $1 \le u, v \le I$ and $u \ne v$,

• $B_p \cap B_q = \emptyset$ for all $1 \le p, q \le J$ and $p \ne q$.

We now define penalized $\beta$-NMF with disjoint linear equality constraints as follows

$$ \begin{array}{ll} \min_{W \in \mathbb{R}_+^{F \times K}, H \in \mathbb{R}_+^{K \times N}} & D_{\beta}(V|WH) + \lambda_1 \Phi_1(W) + \lambda_2 \Phi_2(H) \\ \text{such that} & \alpha_i^T W (K_i) = b_i \text{ for } 1 \le i \le I, \\ & \gamma_j^T H (B_j) = c_j \text{ for } 1 \le j \le J, \end{array} \qquad (2) $$

where

• the penalty functions $\Phi_1(W)$ and $\Phi_2(H)$ are lower bounded and admit a particular upper approximation; see Assumption 1 below.

• $\lambda_1$ and $\lambda_2$ are the penalty weights (nonnegative scalars).

• $\alpha_i \in \mathbb{R}_{++}^{|K_i|}$ ($1 \le i \le I$) and $\gamma_j \in \mathbb{R}_{++}^{|B_j|}$ ($1 \le j \le J$) are vectors with positive entries. Note that if $\alpha_i$ or $\gamma_j$ contains zero entries, the corresponding indices can be removed from $K_i$ and $B_j$.

• $b_i$ ($1 \le i \le I$) and $c_j$ ($1 \le j \le J$) are positive scalars.

As for most NMF algorithms, we propose to resort to a block coordinate descent (BCD) framework to solve problem (2): at each iteration we tackle two sub-problems separately; one in $W$ and the other in $H$. The subproblems in $W$ and $H$ are essentially the same, by symmetry of the model, since transposing the relation $X \approx WH$ gives $X^T \approx H^T W^T$. Hence, we may focus on solving the subproblem in $H$ only, namely

$$ \min_{H \in \mathbb{R}_+^{K \times N}} D_{\beta}(V|WH) + \lambda_2 \Phi_2(H) \quad \text{such that} \quad \gamma_j^T H (B_j) = c_j \text{ for } 1 \le j \le J. \qquad (3) $$

In order to solve (3), we will design MU based on the majorization-minimization (MM) framework [33], which is the standard in the NMF literature; see [12] and the references therein. Let us briefly recall the high-level ideas to obtain MU via MM. Let us consider the general problem

$$ \min_{h \in \mathcal{H}} f(h). $$

Given an initial iterate $\tilde{h} \in \mathcal{H}$, MM generates a new iterate $\hat{h} \in \mathcal{H}$ that is guaranteed to decrease the objective function, that is, $f(\hat{h}) \le f(\tilde{h})$. To do so, it uses the following two steps:

• Majorization: find a function that is an upper approximation of the objective and is tight at the current iterate, which is referred to as a majorizer. More precisely find a function $g(h|\tilde{h})$ such that

$$ (i) g(\tilde{h}|\hat{h}) = f(\tilde{h}) \quad \text{and} \quad (ii) g(h|\hat{h}) \ge f(h) \text{ for all } h \in \mathcal{H}. $$

• Minimization: minimize the majorizer, that is, solve $\min_{h \in \mathcal{H}} g(h|\tilde{h})$ approximately or exactly, to obtain the next iterate $\hat{h} \in \mathcal{H}$ which is such that (iii) $g(\hat{h}|\tilde{h}) \le g(\tilde{h}|\hat{h})$. This guarantees the objective function to decrease at each step of this iterative process since

$$ f(\hat{h}) \le g(\hat{h}|\tilde{h}) \le g(\tilde{h}|\hat{h}) = f(\tilde{h}). $$
The MU for NMF are obtained using MM where the majorizer g is chosen separable, that is, $g(h|\tilde{h}) = \sum_{i=1} g_i(h_i|\tilde{h}_i)$ for some well chosen univariate functions $g_i$'s; see (4) in the next section. This choice typically makes the minimization of g admits a closed-form solution which is multiplicative, that is, it has the form $\hat{h} = \tilde{h} \odot c(\tilde{h})$ where $\odot$ is the component-wise product, and $c(\tilde{h})$ is a nonnegative vector that depends on $\tilde{h}$. We will encounter several examples later in this paper.

In summary, to derive MU for (3), we will follow the MM framework. We first provide a majorizer for the objective of (3) in Section 2. This majorizer has the property to be separable in each entry of $H$. In order to handle the equality constraints, we introduce Lagrange dual variables in Section 2.2, and explain how they can be computed efficiently. This allows us to derive general MU in Section 2.3 in the case of non-penalized $\beta$-NMF under disjoint linear equality constraints. This is showcased on simplex-structured $\beta$-NMF in Section 3. In Section 4, we will illustrate on minimum-volume KL-NMF how to derive MU in the presence of penalty terms.

## 2.1 Separable majorizer for the objective function

Let us derive a majorizer for $\Psi(H) := D_\beta(V|WH) + \lambda\Phi(H)$, that is, a function $G(H|\tilde{H})$ satisfying (i) $G(H|\tilde{H}) \ge \Psi(H)$ for all $H$, and (ii) $G(\tilde{H}|\tilde{H}) = \Psi(\tilde{H})$. Note that, to simplify the presentation, we denote $\Phi_2(H) = \Phi(H)$ and $\lambda = \lambda_2$. To do so, let us analyze each term of $\Psi(H)$ independently.

**Majorizing $D_\beta(V|WH)$** The first term $D_\beta(V|WH)$ can be decoupled into $N$ independent terms, one for each column $h_n$ of $H$, that is, $D_\beta(V|WH) = \sum_{n=1}^N D_\beta(v_n|Wh_n)$, where $v_n$ denotes the nth column of matrix $V$. Let us focus on a specific column of $H$, denoted $h \in \mathbb{R}_+^K$, and the corresponding column of $V$, denoted $v \in \mathbb{R}_+^K$. We majorize $D_\beta(v|Wh) = \sum_{f=1}^F d_\beta(v_f|(Wh)_f)$ following the methodology introduced in [12], which consists in applying a convex-concave procedure [36] to $d_\beta$, as presented in Appendix A. The resulting upper bound is given by

$$d_\beta(v_f|(Wh)_f) \le \sum_{k=1}^K \frac{w_{fk}\tilde{h}_k}{\tilde{v}_f} \tilde{d}\left(v_f|\tilde{v}_f \frac{h_k}{\tilde{h}_k}\right) + \hat{d}'(v_f|\tilde{v}_f) \sum_{k=1}^K w_{fk}(h_k - \tilde{h}_k) + \hat{d}(v_f|\tilde{v}_f), \quad (4)$$

where $w_{fk}$ denotes the entry of matrix $W$ at position $(f, k)$, $\tilde{v}_f := (W\tilde{h})_f$ denotes the fth entry of $\tilde{v}$, and $\hat{d}$ and $\tilde{d}$ are the concave and convex parts of $d$, respectively.

**Majorizing $\Phi(H)$** For the second term $\Phi(H)$, we rely on the following assumption for $\Phi$.

**Assumption 1.** The function $\Phi: \mathbb{R}_+^{K \times N} \to \mathbb{R}$ is lower bounded, and for any $\tilde{H} \in \mathbb{R}_+^{K \times N}$ there exists constants $L_{kn}$ ($1 \le k \le K, 1 \le n \le N$) such that the inequality

$$\Phi(H) \le \Phi(\tilde{H}) + \langle \nabla\Phi(\tilde{H}), H - \tilde{H} \rangle + \sum_{k,n} \frac{L_{kn}}{2} (H_{kn} - \tilde{H}_{kn})^2 \quad (5)$$

is satisfied for all $H \in \mathbb{R}_+^{K \times N}$. (Note that the constants $L_{kn}$ may depend on $\tilde{H}$, this will be the case for example in Section 4).

Let us mention two important classes of functions satisfying Assumption 1.

1. Smooth concave functions that are lower bounded on the nonnegative orthant. For such functions, we can take $L_{kn} = 0$ for all $k, n$ since they are upper approximated by their first-order Taylor approximation. Note that, in this case,

$$\nabla\Phi(\tilde{H}) \ge 0, \quad (6)$$
otherwise we would have $\lim_{y \to \infty} \Phi(H + ye_i e_j^T) = -\infty$, where $e_i$ is the $i$th unit vector, and this would contradict the fact that $\Phi$ is bounded from below. This observation will be useful in the proof of Proposition 1 and is only valid for the special case $L_{kn} = 0$ for all $k, n$.

Examples of such penalty functions include the sparsity-promoting regularizers $\Phi(H) = \|H\|_p^p = \sum_{k,n} H(k,n)^p$ for $0 < p \le 1$ since $H \ge 0$.

2. Lower-bounded functions with Lipschitz continuous gradient for which (5) follows from the descent lemma [3].

Examples of such penalty functions include any smooth convex functions; for example any quadratic penalty, such as $||AH - B||_2^2$ for some matrices $A$ and $B$ in which case $L_{kn} = \sigma_1(A)^2$ for all $k, n$. We will encounter another example later in the paper, namely logdet $(HH^\top + \delta I)$ for $\delta > 0$ which allows to minimize the volume of the rows of $H$; see Section 4 for the details (Note that we will use this regularizer for $W$).

**Majorizing $\Psi(H)$** Combining (4) and (5), we can construct a majorizer for $\Psi(H)$. Since both (4) and (5) are separable in each entry of $H$, their combination is also separable into a sum of $K \times N$ component-wise majorizers, up to an additive constant:

$$G(H|\tilde{H}) = \sum_{n=1}^{N} \sum_{k=1}^{K} g(h_{kn}|\tilde{H}) + C(\tilde{H}), \quad (7)$$

where

$$g(h_{kn}|\tilde{H}) = \sum_{f=1}^{F} \frac{w_{fk}\tilde{h}_{kn}}{\tilde{v}_{fn}} \tilde{d}\left(v_{fn}|\tilde{v}_{fn}\frac{h_{kn}}{\tilde{h}_{kn}}\right) + a_{kn}h_{kn}^2 + p_{kn}h_{kn}, \quad (8)$$

$$C(\tilde{H}) = \sum_{n=1}^{N} \sum_{f=1}^{F} \left( \hat{d}(v_{fn}|\tilde{v}_{fn}) - \sum_{k=1}^{K} \hat{d}(v_{fn}|\tilde{v}_{fn})w_{fk}\tilde{h}_{kn} \right) + a_{kn}\tilde{h}_{kn}^2,$$

with $a_{kn} = \lambda \frac{L_{kn}}{2}$, and

$$p_{kn} = \sum_{f=1}^{F} w_{fk} \hat{d}'(v_{fn}|\tilde{v}_{fn}) + \lambda \left( \frac{\partial \Phi}{\partial h_{kn}}(\tilde{H}) - L_{kn} \tilde{h}_{kn} \right).$$

## 2.2 Dealing with equality constraints via Lagrange dual variables

In the previous section, we derived a majorizer for $\Psi(H)$, $G(H|\tilde{H})$, which is separable in each entry of $H$. Without the equality constraints, we could then compute closed-form solutions to univariate problems to minimize $G(H|\tilde{H})$ to obtain the standard MU for NMF as in [12].

However, in problem (3), the entries of $H$ in the subsets $B_j$ are not independent as they are linked with the equality constraints $\gamma_j^T H(B_j) = c_j$ for $j = 1, 2, \dots, J$. In fact, to minimize the majorizer under the equality constraints, we need to solve

$$\min_{H \in \mathbb{R}_+^{K \times N}} G(H|\tilde{H}) \quad \text{such that} \quad \gamma_j^T H(B_j) = c_j \text{ for } 1 \le j \le J. \quad (9)$$

The variables in different sets $B_j$ can be optimized independently, as they do not interact in the majorizer nor in the constraints. Note that, for the entries of $H$ that do not appear in any constraints,
the standard MU [12] can be used. For simplicity, let us fix $j$ and denote $\mathcal{B} = \mathcal{B}_j$, $Q = |\mathcal{B}|$, $\mathbf{y} = H(\mathcal{B}) \in \mathbb{R}_+^Q$, $\gamma = \gamma_j \in \mathbb{R}_{++}^Q$, and $c = c_j > 0$. The problems we need to solve have the form

$$ \min_{\mathbf{y} \in \mathcal{Y}} G(\mathbf{y}|\tilde{H}), \qquad (10) $$

where $\mathcal{Y} = \{\mathbf{y} \in \mathbb{R}_+^Q | \gamma^T \mathbf{y} = c\}$ and

$$ G(\mathbf{y}|\tilde{H}) = \sum_{(k,n) \in \mathcal{B}} g(h_{kn}|\tilde{H}), \qquad (11) $$

where the component-wise majorizers $g(h_{kn}|\tilde{H})$ are defined by (8). Let us introduce a convenient notation: for $q = 1, 2, \dots, Q$, we denote by $(k(q), n(q))$ the $q$th pair belonging to $\mathcal{B}$. Hence the Lagrangian function of (11) can be written as

$$ G^{\mu}(\mathbf{y}|\tilde{H}) = G(\mathbf{y}|\tilde{H}) - \mu(\gamma^T \mathbf{y} - c) = \mu c + C(\tilde{H}) + \sum_{q=1}^{Q} g^{\mu}(y_q|\tilde{H}), \qquad (12) $$

where

$$ g^{\mu}(y_q|\tilde{H}) = g(y_q|\tilde{H}) - \mu\gamma_q y_q = \sum_{f=1}^{F} \frac{w_{fk(q)}\tilde{y}_q}{\tilde{v}_{fn(q)}} \tilde{d}' \left( v_{fn(q)}|\tilde{v}_{fn(q)} \frac{y_q}{\tilde{y}_q} \right) + a_q y_q^2 + (p_q - \mu\gamma_q)y_q, \qquad (13) $$

$$ p_q = \sum_{f=1}^{F} w_{fk(q)} \tilde{d}'(v_{fn(q)}|\tilde{v}_{fn(q)}) + \lambda \left( \frac{\partial \Phi}{\partial y_q}(\tilde{H}) - L_{k(q)n(q)}\tilde{y}_q \right), \qquad (14) $$

and $\mu \in \mathbb{R}$. Note that $G^{\mu}$ is separable, as is $G$, because the term $\gamma^T \mathbf{y}$ is linear.

Assume for now that the Lagrangian multiplier $\mu$ is known, and let us minimize $G^{\mu}(\mathbf{y}|\tilde{H})$ on $(0, \infty)^Q$. Such a problem is separable under the form of $Q$ subproblems, consisting in minimizing univariate functions $g^{\mu}(\cdot|\tilde{H})$ separately over $(0, \infty)$. We now show in Proposition 1 that, under mild conditions, each subproblem admits a unique solution over $(0, \infty)$.

**Proposition 1.** Let $q \in \{1, 2, \dots, Q\}$. Assume that $\beta < 2$ and $\tilde{y}_q, v_{fn(q)}, w_{fk(q)} > 0$ for all $f$. Moreover, when $\beta \le 1$, assume that $\mu < \frac{p_q}{\gamma_q}$ for all $q$ such that $a_q = 0$. Then there exists a unique minimizer $y_q^*(\mu)$ of $g^{\mu}(y_q|\tilde{H})$ in $(0, \infty)$.

*Proof.* According to Proposition 4 (see Appendix A), each $g^{\mu}$ is $C^{\infty}$ and strictly convex on $(0, \infty)$, so its infimum is uniquely attained in the closure of $(0, \infty)$. We have to prove that it is neither reached at 0 nor at $\infty$. On the one hand, from (13), we have

$$ (g^{\mu})'(y_q|\tilde{H}) = \sum_{f=1}^{F} w_{fk(q)} \tilde{d}' \left( v_{fn(q)}|\tilde{v}_{fn(q)} \frac{y_q}{\tilde{y}_q} \right) + 2a_q y_q + p_q - \gamma_q \mu \qquad (15) $$

and, for any $\beta < 2$ and any $x > 0$,

$$ \lim_{y \to 0^+} \tilde{d}'(x|y) = -\infty, $$
so $\lim_{y_q \to 0^+} (g^\mu)'(y_q | \tilde{H}) = -\infty$, which ensures that the infimum is not reached at 0. On the other hand,

$$ \lim_{y \to \infty} \tilde{d}'(x|y) = \begin{cases} 0 & \text{if } \beta \le 1, \\ \infty & \text{otherwise.} \end{cases} \qquad (16) $$

According to (15) and (16), the distinction must be made between two cases:

* If $a_q > 0$ or $\beta \in (1, 2)$: $\lim_{y_q \to \infty} (g^\mu)'(y_q | \tilde{H}) = \infty$, so the infimum is reached for a finite $y_q$.

* If $a_q = 0$ and $\beta \le 1$: $\lim_{y_q \to \infty} (g^\mu)'(y_q | \tilde{H}) = p_q - \gamma_q \mu$, so the same conclusion holds if $\mu < \frac{p_q}{\gamma_q}$.

□

We just proved that, under mild conditions, each $g^\mu$ has a unique minimizer over $(0, \infty)$. However we assumed that the value of $\mu$ is fixed. Now given $\mathbf{y}^*(\mu) = [y_1^*(\mu), \dots, y_Q^*(\mu)]^T$, let us show that the solution to $\gamma^T \mathbf{y}^*(\mu) = c$ is unique. The corresponding value of $\mu$, which we denote $\mu^*$, provides the minimizer $\mathbf{y}^*(\mu^*)$ of $G^\mu(\mathbf{y} | \tilde{H})$ that satisfies the linear constraint $\gamma^T \mathbf{y}^*(\mu^*) = c$. Moreover, $\mu^*$ naturally fulfills $\mu^* < \frac{p_q}{\gamma_q}$ for all $q$ when $\beta \le 1$ and $a_q = 0$, as required in Proposition 1.

**Proposition 2.** Assume that $\beta < 2$ and $\tilde{y}_q, v_{fn(q)}, w_{fk(q)} > 0$ for all $q, f$. Then the scalar equation $\gamma^T \mathbf{y}^*(\mu) = c$ in the variable $\mu$ admits a unique solution $\mu^*$ in $(-\infty, t)$, where

$$ t = \min_{1 \le q \le Q} t_q, \quad \text{where } t_q = \begin{cases} \frac{p_q}{\gamma_q} & \text{if } \beta \le 1 \text{ and } a_q = 0, \\ \infty & \text{otherwise,} \end{cases} \qquad (17) $$

so that $\mathbf{y}^*(\mu^*) \in (0, \infty)^Q$ is the unique solution to problem (10).

*Proof.* Under the conditions of Proposition 1, $g^\mu(y_q | \tilde{H})$ has a unique minimizer $y_q^*(\mu)$ for each $j$. By the first-order optimality condition, $y_q^*(\mu)$ is a solution of $(g^\mu)'(y_q | \tilde{H}) = 0$ or equivalently, by (15), a solution of $\gamma_q^{-1} g'(y_q | \tilde{H}) = \mu$ over $(0, \infty)$ where

$$ \gamma_q^{-1} g'(y_q | \tilde{H}) = \gamma_q^{-1} \sum_{f=1}^{F} w_{fk(q)} \tilde{d}' \left( v_{fn(q)} | \tilde{v}_{fn(q)} \frac{y_q}{\tilde{y}_q} \right) + 2 \frac{a_q}{\gamma_q} y_q + \frac{p_q}{\gamma_q} \qquad (18) $$

is strictly increasing on $(0, \infty)$ (since $g$ is strictly convex) and one-to-one from $(0, \infty)$ to an open interval $T_q = (t_q^-, t_q^+)$ where

$$ t_q^- = \lim_{y_q \to 0} g'(y_q | \tilde{H}) = -\infty, \qquad (19) $$

$$ t_q^+ = \lim_{y_q \to \infty} g'(y_q | \tilde{H}) = t_q. \qquad (20) $$

Moreover, $p_q \ge 0$ if $a_q = 0$ (then $L = 0$) and $\beta \le 1$ according to (6) and (14). As a consequence, $\gamma_q^{-1} g'(y_q^* | \tilde{y}_q) = \mu$ is equivalent to

$$ y_q^*(\mu) = (g')^{-1}(\gamma_q \mu), \qquad (21) $$

where $\mu \in T_q$ and $(g')^{-1}$ denotes the inverse function of $g'$.
Coming back to the multivariate problem (10), we must find a value $\mu^*$ of the Lagrangian multiplier such that the constraint $\gamma^T y^*(\mu) = c$ is satisfied. Given (21), $\mu^*$ is a solution of

$$ \sum_{q=1}^{Q} \gamma_q (g')^{-1} (\gamma_q \mu) = c. \qquad (22) $$

Each $g'(y_q|\tilde{H})$ being strictly increasing on $(0, \infty)$, $(g')^{-1}(\gamma_q \mu)$ is also strictly increasing (from $T_q$ to $(0, \infty)$), this is a direct consequence of $(f^{-1})' = \frac{1}{f' \circ f^{-1}}$ where $f$ is any strictly increasing function on some interval. Finally $\sum_{q=1}^Q \gamma_q (g')^{-1}(\gamma_q \mu)$ is strictly increasing from $\cap_{j=1}^J T_q = (-\infty, t)$ to $(0, \infty)$, with $t \ge 0$. Therefore, the solution $\mu^*$ is unique. $\square$

Proposition 2 shows that the optimal Lagrangian multiplier is the unique solution of (22). Finding the solution of (22) is equivalent to finding the root of a function $r(\mu)$. We propose here-under to use a Newton-Raphson method to compute $\mu^*$, and show that this method generates a sequence of iterates $\mu_n$ that converges towards $\mu^*$ at a quadratic speed.

**Proposition 3.** Assume that $\beta < 2$ and $\tilde{y}_q, v_{fn(q)}, w_{fk(q)} > 0$ for all $q, f$. Let

$$ r(\mu) = \sum_{q=1}^{Q} \gamma_q (g')^{-1} (\gamma_q \mu) - c $$

for $\mu \in (-\infty, t)$ where $t$ is defined in (17), and denote $\mu^*$ the unique solution of $r(\mu) = 0$. From any initial point $\mu_0 \in (\mu^*, t)$, Newton-Raphson's iterates

$$ \mu_{n+1} = \mu_n - \frac{r(\mu_n)}{r'(\mu_n)} $$

decrease towards $\mu^*$ at a quadratic speed.

*Proof.* We already know that $r$ is strictly increasing from $(-\infty, t)$ to $(0, \infty)$. Let us show that $r$ is also strictly convex. According to the third item of Proposition 4 in Appendix A, $\tilde{d}''(x|y)$ is completely monotonic, so it is strictly decreasing in $y$. Equivalently, $\tilde{d}'(x|y)$ is strictly concave in $y$, and each $g'$ is also strictly concave according to (18). Since the inverse of a strictly increasing, strictly concave function $f$ is strictly increasing and strictly convex, which is a direct consequence of $(f^{-1})'' = -\frac{f'' \circ f^{-1}}{(f' \circ f^{-1})^3}$, then each $(g')^{-1}$ is strictly convex, and finally, $r$ is strictly convex.

For any $\mu_0 \in (\mu^*, t)$, we have $r(\mu_0) > 0$, so $\mu_1 = \mu_0 - \frac{r(\mu_0)}{r'(\mu_0)} < \mu_0$. We have also $\mu_1 > \mu^*$ as a consequence of the strict convexity of $r$. By immediate recurrence, we obtain that $\mu_n$ is a decreasing series that converges towards $\mu^*$. According to [30], it converges at a quadratic speed since $|r'|$ and $|r''|$ are bounded away from 0 in $[\mu^*, \mu_0]$. $\square$

**Discussion** At this point, we have derived an optimization framework to tackle problem (10). The optimal Lagrangian multiplier value is determined before each majorization-minimization update using a Newton-Raphson algorithm. However, such a formal solution is implementable if and only if each $y_q^*(\mu)$ can be actually computed as the minimizer of $g^\mu(y_q|\tilde{H})$ in $(0, \infty)$. In some cases, computing $y_q^*(\mu)$ is equivalent to extracting the roots of a polynomial of a degree smaller or equal to four, which is possible in closed form. In other cases, we have to solve a polynomial equation of degree larger
<table><caption>Table 1: Cases where (21) can be computed in closed form. They are indicated by the degree of the corresponding polynomial equation, otherwise the symbol $\ddagger$ is used. The constants $L_{kn}$ are the one needed in Assumption 1 for the penalization functions $\Phi_1(H)$ and $\Phi_2(W)$; see (5).</caption><thead><tr><th></th><th>β ∈ (-∞, 1)\{0}</th><th>β = 0</th><th>β = 1</th><th colspan="4">β ∈ (1, 2)</th></tr><tr><th></th><th></th><th></th><th></th><th colspan="2">5/4 | 4/3 | 3/2</th><th>other</th></tr></thead><tbody><tr><th>No penalization, or<br>L<sub>kn</sub> = 0 for all k, n</th><td>1</td><td>1</td><td>1</td><td>3</td><td>4</td><td>2</td><td>‡</td></tr><tr><th>L<sub>kn</sub> > 0 for some k, n</th><td>‡</td><td>3</td><td>2</td><td>‡</td><td>‡</td><td>3</td><td>‡</td></tr></tbody></table>

than four, or even an equation that is not polynomial. Table 1 indicates the cases where a closed-form solution is available, and hence when our framework can be efficiently implemented. We observe that, without penalization or with penalization satisfying $L_{kn} = 0$ for all $k, n$ (e.g., smooth concave functions), the polynomial equation is of degree one, and hence always admit a closed form for $\beta \le 1$ and $\beta \in \{\frac{5}{4}, \frac{4}{3}, \frac{3}{2}\}$. This particular case is discussed in the next section, which we will exemplify in Section 3 with $\beta$-NMF with sum-to-one constraints on the columns of $H$. In Section 4, we will present an important example with $L_{kn} > 0$ for all $k, n$ and $\beta = 1$, namely minimum-volume KL-NMF.

## 2.3 MU for $\beta$-NMF with disjoint linear equality constraints without penalization

In this section, we derive an algorithm based on the general framework presented in the previous section to tackle the $\beta$-NMF problem under disjoint linear equality constraints without penalization, that is, problem (2) with $\lambda_1 = \lambda_2 = 0$. We consider this simplified case here as it allows to provide explicit MU for any value of $\beta < 2$; see the row ‘No penalization’ of Table 1. These updates satisfy the constraints after each update of $W$ or $H$, and monotonically decrease the objective function $D_{\beta}(V|WH)$.

Let us then consider the subproblem of (2) over $H$ when $W$ is fixed and with $\lambda_2 = 0$, that is,

$$ \min_{H \in \mathbb{R}_{+}^{K \times N}} D_{\beta}(V|WH) \quad \text{such that} \quad \gamma_j^T H(\mathcal{B}_j) = c_j \text{ for } 1 \le j \le J. \qquad (23) $$

Let us follow the framework presented above. First, an auxiliary function, which we denote $G(H|\tilde{H})$, is constructed at the current iterate $\tilde{H}$ so that it majorizes the objective for all $H$ and is defined as follows:

$$ G(H|\tilde{H}) = \sum_{f,n} \left[ \sum_k \frac{w_{fk} \tilde{h}_{kn}}{\tilde{v}_{fn}} \tilde{d} \left( v_{fn} \left| \tilde{v}_{fn} \frac{h_{kn}}{\tilde{h}_{kn}} \right. \right) \right] + \left[ \hat{d}(v_{fn}|\tilde{v}_{fn}) \sum_k w_{fk} (h_{kn} - \tilde{h}_{kn}) + \hat{d}(v_{fn}|\tilde{v}_{fn}) \right], \qquad (24) $$

where $\tilde{d}(\cdot|\cdot)$ and $\hat{d}(\cdot|\cdot)$ are given in Appendix A. Second, we need to minimize $G(H|\tilde{H})$ while imposing the set of linear constraints $\gamma_j^T H(\mathcal{B}_j) = c_j$. The Lagrangian function of $G$ is given by

$$ G^{\mu}(H|\tilde{H}) = G(H|\tilde{H}) - \sum_{j}^{J} [\mu_{j} (\gamma_{j}^{T} H(\mathcal{B}_{j}) - c_{j})], \qquad (25) $$

where $\mu_j$ are the Lagrange multipliers associated to each linear constraint $\gamma_j^T H(\mathcal{B}_j) = c_j$. We observe that $G^\mu$ in (25) is a separable majorizer in the variables $H$ of the Lagrangian function $D_\beta(V|WH)$ –
$\sum_j^J [\mu_j (\gamma_j^T H(\mathcal{B}_j) - c_j)]$. Due to the disjoitness of each subset of variables $\mathcal{B}_j$ (25), we only consider the optimization over one specific subset $\mathcal{B}_j$. The minimizer (21) of $G^\mu(H(\mathcal{B}_j) | \tilde{H}(\mathcal{B}_j))$ has the following component-wise expression:

$$H^*(\mathcal{B}_j) = \tilde{H}(\mathcal{B}_j) \odot \left( \frac{[C(\mathcal{B}_j)]}{[D(\mathcal{B}_j) - \mu_j \gamma_j]} \right)^{\eta(\beta)}, \quad (26)$$

where $C = W^T ((WH)^{(\beta-2)} \odot V)$, $D = W^T (WH)^{(\beta-1)}$, $\eta(\beta) = \frac{1}{2-\beta}$ for $\beta \le 1$, and $\eta(\beta) = \frac{1}{\beta-1}$ for $\beta \ge 2$ [12, Table 2], $A \odot B$ (resp. $\frac{[A]}{[B]}$) is the Hadamard product (resp. division) between A and B, $A^\alpha$ is the element-wise $\alpha$ exponent of A. The case $\beta \in (1, 2)$ is more difficult: we need to find a root of a function of the form $\mu + bx^{\beta-1} - cx^{\beta-2} = 0$. For example, for $\beta = \frac{3}{2}$, we have $\mu + bx^{1/2} - cx^{-1/2} = 0$. Using $y = \sqrt{x}$, and after simplifications, we obtain $\mu y + by^2 - c = 0$ leading to the positive root $x = \left( \frac{\sqrt{\mu^2+4bc-\mu}}{2b} \right)^2$.

According to Proposition 2, (26) is a well-defined update from $(0, \infty)^Q$ to itself, provided that $\mu_j$ is tuned to $\mu_j^*$. This brings us a structural guarantee that $D(\mathcal{B}_j) - \mu_j^*\gamma_j$ cannot cancel.

Finally, we need to evaluate $\mu_j^*$, which is uniquely determined on some interval $(-\infty, t)$ according to Proposition 2. This amounts to solve $\gamma_j^T H^*(\mathcal{B}_j) = c_j$. When $\beta \notin (1, 2)$, this is equivalent to find the root of the function

$$r_j(\mu_j) = \sum_{q=1}^{Q} \gamma_{j,q} \left[ \tilde{H}(\mathcal{B}_j) \odot \left( \frac{[C(\mathcal{B}_j)]}{[D(\mathcal{B}_j) - \mu_j \gamma_j]} \right)^{\eta(\beta)} \right]_q - c_j, \quad (27)$$

where $[A]_q$ denotes the q-th entry of expression A. Indeed, $r_j(\mu_j)$ is a finite sum of elementary rational functions of $\mu_j$ and each of them is an increasing, convex function in $\mu_j$ over $(-\infty, t_q)$ with $t_q = \frac{D_q(\mathcal{B}_j)}{\gamma_{j,q}}$ for each $q$. It is even completely monotone for all $\mu$ in $(-\infty, t_q)$ because $\eta(\beta) > 0$ [28]. As a consequence $r_j(\mu_j)$ is also a completely monotone, convex increasing function of $\mu_j$ in $(-\infty, t)$, where $t = \min(t_q)$. Finally, we can easily show that the function $r_j(\mu_j)$ changes of sign on the interval $(-\infty, t)$ by computing two limits at the closure of the interval. As $\mu^* \in (-\infty, t)$, the update (26) is nonnegative. To evaluate $\mu^*$, we use a Newton-Raphson method, with any initial point $\mu_0 \in (\mu^*, t)$, with a quadratic rate of convergence as demonstrated in Proposition 3. Algorithm 1 summarizes our method to tackle (2) for all the $\beta$-divergences, $\beta \notin (1, 2)$ which we refer to as disjoint-constrained $\beta$-NMF algorithm. The update for matrix $W$ can be derived in the same way, by symmetry of the problem. For $\beta \in (1, 2)$, a case-by-case analysis could be carried out for the values of $\beta$ for which the minimizer of (25) takes a closed-form expression.

**Remark 1.** As noted above, the denominators of (26) and (27) will be different from zero. This follows notably from our assumption that $(W, H) > 0$; see Propositions 1, 2 and 3. This is a standard assumption in the NMF literature: the entries of $(W, H)$ are initialized with positive entries which ensures all iterates to remain positive. This is important because the MU cannot change an entry equal to zero [26]; this is the so-called zero-locking phenomenon. This implies $C$ and $D$ in (26) and (27) are positive matrices (as long as $V$ has at least one nonzero entry per row and column). In practice, one should however be careful because some entries of $W$ and $H$ can numerically be set to zero (because of finite precision). Hence, in our implementation, we use the machine precision as a lower bound for the entries of $W$ and $H$, as recommended in [17].
**Computational cost** The computational cost of Algorithm 1 is asymptotically equivalent to the standard MU for $\beta$-NMF, that is, it requires $\mathcal{O}(FNK)$ operations per iteration. Indeed, the complexity is mainly driven by matrix products required to compute $C$ and $D$; see (26). To compute the roots of (27) corresponding to $H$ using Newton-Raphson, each iteration requires to compute $r_j(\mu_j)/r'_j(\mu_j)$ for all $j$ which requires $\mathcal{O}(KN)$ operations (when every entry of $H$ appears in a constraint). Finding the roots therefore requires $\mathcal{O}(KN)$ operations times the number of Newton-Raphson iterations. By symmetry, it requires $\mathcal{O}(KF)$ operations to compute the roots corresponding to $W$. Because of the quadratic convergence, the number of iterations required for the convergence of the Newton-Raphson method is typically small, namely between 10 to 100 in our experiments using the stopping criterion $|r(\mu_j)| \le 10^{-6}$ for all $j$. Therefore, in practice, the overall complexity of Algorithm 1 is dominated by the matrix products that require $\mathcal{O}(FNK)$ operations. The same conclusions apply to the algorithms presented in Sections 3, 4 and 5, and this will be confirmed by our numerical experiments.

<caption>**Algorithm 1** $\beta$-NMF with disjoint linear constraints</caption>

**Input:** A matrix $V \in \mathbb{R}^{F \times N}$, an initialization $H \in \mathbb{R}_+^{K \times N}$ and $W \in \mathbb{R}^{F \times K}$, a factorization rank $K$, a maximum number of iterations, maxiter, a value for $\beta \notin (1, 2)$, and the linear constraints defined by $K_i$, $\alpha_j$ and $b_i$ for $i = 1, 2, \dots, I$, and $B_j$, $\gamma_j$ and $c_j$ for $j = 1, 2, \dots, J$.

**Output:** A rank-$K$ NMF $(W, H)$ of $V$ satisfying constraints in (2).

1: **for** $i = 1$ : maxiter **do**

2: &nbsp;&nbsp;&nbsp;&nbsp;% Update of matrix $H$

3: &nbsp;&nbsp;&nbsp;&nbsp;$C \leftarrow W^T ((WH)^{(\beta-2)} \odot V)$

4: &nbsp;&nbsp;&nbsp;&nbsp;$D \leftarrow W^T (WH)^{(\beta-1)}$

5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1$ : $J$ **do**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mu_j \leftarrow \text{root}(r_j(\mu_j))$ &nbsp;&nbsp;&nbsp;&nbsp;% see Equation (27)

7: &nbsp;&nbsp;&nbsp;&nbsp;$H(B_j) \leftarrow H(B_j) \odot \left( \frac{[C(B_j)]}{[D(B_j) - \mu_j \gamma_j]} \right)^{(\eta(\beta))}$

8: &nbsp;&nbsp;&nbsp;&nbsp;**end for**

9: &nbsp;&nbsp;&nbsp;&nbsp;$B_c = \{(k, n) \mid 1 \le k \le K, 1 \le n \le N\} \setminus (\cup_j^J B_j)$ &nbsp;&nbsp;&nbsp;&nbsp;% $B_c$ is the complement of $\cup_j^J B_j$

10: &nbsp;&nbsp;&nbsp;&nbsp;$H(B_c) \leftarrow H(B_c) \odot \left( \frac{[C(B_c)]}{[D(B_c)]} \right)^{(\eta(\beta))}$

11: &nbsp;&nbsp;&nbsp;&nbsp;% Update of matrix $W$

12: &nbsp;&nbsp;&nbsp;&nbsp;$W$ is updated in the same way as $H$, by symmetry of the problem.

13: **end for**

## 3 Showcase 1: Simplex-structured $\beta$-NMF

In this section, we showcase a particularly important example of $\beta$-NMF with linear disjoint constraints and no penalization, namely, the simplex-structured matrix factorization (SSMF) problem. It is defined as follows: given a data matrix $V \in \mathbb{R}^{F \times N}$ and a factorization rank $K$, SSMF refers to the problem of computing $W$ and $H$ such that $V \approx WH$ and the columns of $H$ lie on the unit simplex, that is, the entries of each column of $H$ are nonnegative and sum to one. SSMF is a powerful tool in many applications such as hyperspectral unmixing in geoscience and remote sensing [4, 27, 1], document analysis [5], and self-modeling curve resolution [29]. We refer the reader to the recent survey [13] for more applications and details about SSMF.
To understand the underlying significance of SSMF, it is necessary to give more insights on a research topic for which important SSMF techniques were initially developed which is the blind Hyperspectral Unmixing (HU), a main research topic in remote sensing. The task of blind HU is to decompose a remotely sensed hyperspectral image into endmember spectral signatures and the corresponding abundance maps with limited prior information, usually the only known information being the number of endmembers. In this context, the columns of W correspond to the endmembers spectral signatures and the columns of H contain the proportion of the endmembers in each column of V, so the column-stochastic assumption for H naturally holds. The nonnegativity of W follows from the nonnegativity of the spectral signatures. We refer to the corresponding problem as simplex-structured nonnegative matrix factorization with the β-divergence (β-SSNMF), and is formulated as follows:

$$ \min_{W \in \mathbb{R}_+^{F \times K}, H \in \mathbb{R}_+^{K \times N}} D_{\beta}(V|WH) \quad \text{such that} \quad e^T h_j = 1 \text{ for } 1 \le j \le N, \qquad (28) $$

where e is the vector of all ones of appropriate dimension. This is particular case of (2) where

* the subsets $B_j$ correspond to the columns of H, and there is no subset $K_i$ (no constraint on W),

* $\gamma_j^T = e$ and $c_j = 1$ for $j = 1, 2, \dots, N$.

Hence Algorithm 1 can be directly applied to (28).

**Numerical experiments** Let us perform numerical experiments to evaluate the effectiveness of Algorithm 1 on the simplex-structure β-NMF problem against existing methods. To the best of our knowledge, the so-called group robust NMF (GR-NMF) algorithm<fnref n="1" /> from [14] is the most recent algorithm that is able to tackle problem (28) for the full range of β-divergences. The approach is not based on Lagrangian multipliers but introduces a change of variables for matrix H. This approach, initially used for NMF in [10], does not provide an auxiliary function for the subproblem in H and resort to a heuristic commonly used in NMF, see for example [34, 11]. Therefore there is no guarantee that the objective function is decreasing at each update of the abundance matrix, unlike Algorithm 1. We apply Algorithm 1 and GR-NMF on three widely used real hyperspectral data sets<fnref n="2" /> [37]:

* Samson: 156 spectral bands with 95×95 pixels, containing mostly 3 materials ($K = 3$), namely “Soil”, “Tree” and “Water”.

* Jasper Ridge: 198 spectral bands with 100×100 pixels, containing mostly 4 materials ($K = 4$), namely “Road”, “Soil”, “Water” and “Tree”.

* Cuprite: 188 spectral bands with 250×190 pixels, containing mostly 12 types of minerals ($K = 12$).

β-SSNMF has shown itself as a powerful one to tackle blind HU, hence this comparative study between Algorithm 1 and GR-NMF [14] focuses on the convergence aspects including the evolution of the objective function and the runtime. The algorithms are compared<fnref n="3" /> for $\beta \in \{0, \frac{1}{2}, 1, \frac{3}{2}, 2\}$. To

<footnote n="1">¹https://www.irit.fr/~Cedric.Fevotte/extras/tip2015/code.zip</footnote>

<footnote n="2">²http://lesun.weebly.com/hyperspectral-data-set.html</footnote>

<footnote n="3">³For $\beta = 3/2$, we had an error in our derivations, and use (26) with $\eta(\beta) = 1$ for $\beta \in (1, 2)$; see the discussion after (26). However, the corresponding MU always decreases the objective function values (which we were monitoring), although we do not have a theoretical justification for this. A possible approach to obtain such as result would be to come up with a majorizer of the majorizer that has a closed-form minimizer given by (26) with $\eta(\beta) = 1$ for $\beta \in (1, 2)$.</footnote>
report the results, we use the relative objective function, denoted $\bar{F}(W, H)$ and defined as<fnref n="4" />

$$ \bar{F}(W, H) = \frac{D_{\beta}(V|WH)}{D_{\beta}(V|vee e^T)}, $$

where $v = \frac{e^T V e}{F_N}$ is the average of the entries of $V$. The relative error $\bar{F}$ should be between 0 and 1: it is equal to 0 for an exact decomposition with $V = WH$, and is equal to 1 for a trivial rank-one approximation where all entries are equal to the average of the entries of $V$. This allows to meaningfully interpret the results, especially since we consider in this comparative study multiple values for $\beta$. In fact, the degree of homogeneity of the $\beta$-divergence is a function of $\beta$. For example, if all the entries of the input matrix are multiplied by 10 and keeping the same NMF solution properly scaled, the squared Frobenius error ($\beta = 2$) is multiplied by 100 while the IS-divergence ($\beta = 0$) is not affected.

As for all tests performed in this paper, the algorithms are tested on a desktop computer with Intel Core i7-8700@3.2GHz CPU and 32GB memory. The codes are written in MATLAB R2018a, and available from https://sites.google.com/site/nicolasgillis/. For all simulations, the algorithms are run for 20 random initializations of $W$ and $H$ (each entry sampled from the uniform distribution in $[0, 1]$). Table 2 reports the average and standard deviation of the runtime (in seconds) as the final value for the relative objective function over these 20 runs for a maximum of 300 iterations.

<table><caption>Table 2: Runtime performance in seconds and final value of relative objective function $\bar{F}_{\text{end}}(W, H)$ for Algorithm 1 and the GR-NMF reported for $\beta \in \{0, \frac{1}{2}, 1, \frac{3}{2}, 2\}$. The table reports the average and standard deviation over 20 random initializations with a maximum of 300 iterations for three hyperspectral data sets. A bold entry indicates the best value for each experiment.</caption><thead><tr><th rowspan="2">Algorithms</th><th colspan="2">Samson</th><th colspan="2">Jasper Ridge</th><th colspan="2">Cuprite</th></tr><tr><th>runtime (s.)</th><th>F&#x0303;<sub>end</sub>(W, H)</th><th>runtime (s.)</th><th>F&#x0303;<sub>end</sub>(W, H)</th><th>runtime (s.)</th><th>F&#x0303;<sub>end</sub>(W, H)</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center;">β = 2</td></tr><tr><td>Algorithm 1<br>GR-NMF</td><td><b>16.62</b>±<b>0.15</b><br>18.23±0.29</td><td>(<b>1.89</b>±<b>0.04</b>)10<sup>-3</sup><br>(1.91±0.05)10<sup>-3</sup></td><td><b>22.86</b>±<b>0.08</b><br>25.32±0.16</td><td>(<b>4.68</b> ± <b>0.39</b>)10<sup>-3</sup><br>(5.87 ± 1.22)10<sup>-3</sup></td><td>121.04 ± 0.62<br><b>114.27</b> ± <b>0.20</b></td><td>(<b>0.98</b> ± <b>0.06</b>)10<sup>-3</sup><br>(1.29 ± 0.07)10<sup>-3</sup></td></tr><tr><td colspan="7" style="text-align:center;">β = 3/2</td></tr><tr><td>Algorithm 1<br>GR-NMF</td><td><b>63.69</b>±<b>0.40</b><br>80.09±0.60</td><td>(<b>2.52</b> ± <b>0.78</b>)10<sup>-3</sup><br>(2.60 ± 0.63)10<sup>-3</sup></td><td><b>89.23</b> ± <b>0.30</b><br>112.72 ± 0.67</td><td>(<b>4.92</b> ± <b>0.29</b>)10<sup>-3</sup><br>(6.32 ± 1.37)10<sup>-3</sup></td><td><b>421.49</b> ± <b>2.79</b><br>508.57 ± 3.50</td><td>(<b>1.54</b> ± <b>0.07</b>)10<sup>-3</sup><br>(2.01 ± 0.09)10<sup>-3</sup></td></tr><tr><td colspan="7" style="text-align:center;">β = 1</td></tr><tr><td>Algorithm 1<br>GR-NMF</td><td><b>18.33</b> ± <b>0.08</b><br>44.78 ± 0.18</td><td>(<b>3.54</b> ± <b>0.27</b>)10<sup>-3</sup><br>(3.77 ± 0.38)10<sup>-3</sup></td><td><b>24.82</b> ± <b>0.35</b><br>62.83 ± 0.76</td><td>(<b>6.07</b> ± <b>0.21</b>)10<sup>-3</sup><br>(7.26 ± 1.50)10<sup>-3</sup></td><td><b>182.98</b> ± <b>14.14</b><br>370.25 ± 21.33</td><td>(<b>2.07</b> ± <b>0.09</b>)10<sup>-3</sup><br>(2.67 ± 0.10)10<sup>-3</sup></td></tr><tr><td colspan="7" style="text-align:center;">β = 1/2</td></tr><tr><td>Algorithm 1<br>GR-NMF</td><td><b>89.80</b> ± <b>0.65</b><br>102.21 ± 0.72</td><td>(7.21 ± 0.75)10<sup>-3</sup><br>(<b>6.93</b> ± <b>0.88</b>)10<sup>-3</sup></td><td><b>126.43</b> ± <b>0.61</b><br>141.75 ± 0.69</td><td>(<b>1.08</b> ± <b>0.10</b>)10<sup>-2</sup><br>(1.12 ± 0.13)10<sup>-2</sup></td><td>682.80 ± 3.32<br><b>642.49</b> ± <b>1.22</b></td><td>(<b>3.13</b> ± <b>0.15</b>)10<sup>-3</sup><br>(3.14 ± 0.14)10<sup>-3</sup></td></tr><tr><td colspan="7" style="text-align:center;">β = 0</td></tr><tr><td>Algorithm 1<br>GR-NMF</td><td><b>52.89</b>±<b>0.54</b><br>55.61±0.47</td><td>(4.60 ± 0.66)10<sup>-2</sup><br>(<b>4.22</b> ± <b>0.79</b>)10<sup>-2</sup></td><td><b>69.59</b>±<b>0.44</b><br>77.87±0.63</td><td>(<b>3.76</b> ± <b>0.11</b>)10<sup>-2</sup><br>(3.76 ± 0.44)10<sup>-2</sup></td><td>479.84 ± 16.02<br><b>354.65</b> ± <b>6.01</b></td><td>(4.39 ± 0.31)10<sup>-3</sup><br>(<b>3.35</b> ± <b>0.10</b>)10<sup>-3</sup></td></tr></tbody></table>

We observe that Algorithm 1 outperforms the GR-NMF in terms of runtime and final values for the relative objective function for all test cases except when $\beta = 0$ for the Samson and Cuprite data sets. In particular, for $\beta = 1$, Algorithm 1 is up to 2.5 times faster than the GR-NMF. For the Cuprite data set with $\beta = 1/2$, Algorithm 1 and GR-NMF perform similarly. We also observe that the standard

<footnote n="4">⁴For the Frobenius norm, that is, $\beta = 2$, the relative error is typically defined as $\frac{D_{\beta}(V|WH)}{D_{\beta}(V|0)}$ meaning that the trivial solution used is the all-zero matrix. However, for other $\beta$-divergences, the value of $D_{\beta}(V|0)$ might not be defined; in particular, for $\beta \le 1$ and $v_{fn} > 0$ for some $f, n$.</footnote>
deviations obtained with Algorithm 1 are in general significantly smaller for all $\beta$, except for $\beta = 0$ for the Samson and Cuprite data sets.

In the supplementary material S1, we provide figures that show the evolution of the relative objective function values with respect to iterations, and that confirm the observations above.

# 4 Showcase 2: minimum-volume KL-NMF

In this section, we showcase another important example of $\beta$-NMF with linear disjoint constraints, namely, the minimum volume NMF with the $\beta$-divergences (min-vol $\beta$-NMF) model. This model is based on the minimization of $\beta$-divergences including a penalty term promoting solutions with minimum volume spanned by the columns of the matrix $W$. It is defined as follows:

$$ \min_{W \in \mathbb{R}_+^{F \times K}, H \in \mathbb{R}_+^{K \times N}} D_{\beta}(V|WH) + \lambda \mathrm{vol}(W) \quad \text{such that} \quad W^T e = e. \qquad (29) $$

where $\lambda$ is a penalty parameter, and $\mathrm{vol}(W)$ is a function measuring the volume spanned by the columns of $W$. In [25], the authors use $\mathrm{vol}(W) = \logdet(W^T W + \delta I)$, where $\delta$ is a small positive constant that prevents $\logdet(W^T W)$ to go to $-\infty$ when $W$ tends to a rank-deficient matrix (that is, when $\mathrm{rank}(W) < K$). This model is particularly powerful as it leads to identifiability which is crucial in many applications such as in hyperspectral imaging or audio source separation [13]. Indeed, under some mild assumptions and in the exact case, authors prove in [25] that (29) is able to identify the groundtruth factors $(W^\#, H^\#)$ that generated the input data $V$, in the absence of noise. In [25], (29) is used for blind audio source separation. In a nutshell, blind audio source separation consists in isolating and extracting unknown sources based on an observation of their mix recorded with a single microphone$^5$. We have to mention that model (29) is also well suited for hyperspectral imaging as discussed in [16].

In the next subsections, we show that we can tackle the min-vol $\beta$-NMF optimization problem defined in (29) with the general framework presented in Section 2 in the case $\beta = 1$.

## 4.1 Problem formulation and algorithm

As the minimum-volume penalty of model (29) concerns matrix $W$ only, the main challenge concerns the update of $W$. Indeed, the update of $H$ is simply the one from [23]. Let us therefore consider the subproblem in $W$ for $H$ fixed:

$$ \min_{W \in \mathbb{R}_+^{F \times K}} D_{\beta}(V|WH) + \lambda \logdet(W^T W + \delta I) \quad \text{such that} \quad e^T w_i = 1 \text{ for } 1 \le i \le K. \qquad (30) $$

Compared to the general model (2), we have that

* the subsets $K_i$ correspond to the columns of $W$, and there is no subset $B_j$,

* $\alpha_i^T = e$ and $b_i = 1$ for $1 \le i \le K$.

To upper bound $\logdet(W^T W + \delta I)$ as required by (5) in Assumption 1, we majorize it using a convex quadratic separable auxiliary function provided in [25, Eq. (3.6)] and which is derived as follows.

<footnote>$^5$We invite the interested reader to watch the video https://www.youtube.com/watch?v=1BrpxvpghKQ to see the application of min-vol KL-NMF on the decomposition of a famous song from the city of Mons.</footnote>
First, the concave function logdet($Q$) for $Q > 0$ can be upper bounded using the first-order Taylor approximation: for any $\tilde{Q} > 0$,

$$ \logdet(Q) \le \logdet(\tilde{Q}) + \langle \tilde{Q}^{-1}, Q - \tilde{Q} \rangle = \langle \tilde{Q}^{-1}, Q \rangle + \text{cst}, $$

where $\text{cst}$ is some constant independent of $Q$. For any $W, \tilde{W}$, and denoting $\tilde{Q} = \tilde{W}^T \tilde{W} + \delta I > 0$, we obtain

$$ \logdet(W^T W + \delta I) \le \langle \tilde{Q}^{-1}, W^T W \rangle + \text{cst} = \text{trace}(W \tilde{Q}^{-1} W^T) + \text{cst}, $$

which is a convex quadratic and Lipschitz-smooth function in $W$. In fact, letting $\tilde{Q}^{-1} = DD^T$ be a decomposition (such as Cholesky) of $\tilde{Q}^{-1} > 0$, we have $\text{trace}(W \tilde{Q}^{-1} W) = \|WD\|_F^2$, from which (5) can be derived easily; see [25] for the details. With this and following our framework from Section 2, we obtain the Lagrangian function

$$ G^{\mu}(W|\tilde{W}) = \sum_{f} G(w_{f}|\tilde{w}_{f}) + \lambda \left( \sum_{f} \bar{l}(w_{f}|\tilde{w}_{f}) + c \right) + \mu^{T} \sum_{f} \left( w_{f} - \frac{1}{F} e \right), \quad (31) $$

where $w_f \in$ denotes the $f$-th row of $W$, $G$ is given by (24), $\bar{l}$ by [25, Eq. (3.6)] and derived as explained above, and $c$ is a constant. Let $\mu$ is the vector Lagrange multipliers of dimension $K$ associated to each linear constraint $e^T w_i = 1$. Exactly as before (hence we omit the details here), $G^\mu$ is separable and, given $\mu$, one can compute the closed-form solution:

$$ W^*(\mu) = \tilde{W} \odot \frac{\left[ \left[ \left[ C + e\mu^T \right]^2 + S \right]^{\frac{1}{2}} - (C + e\mu^T) \right]}{[D]}, \quad (32) $$

where $C = e_{F,N}H^T - 4\lambda(\tilde{W}Y^-)$, $D = 4\lambda\tilde{W}(Y^+ + Y^-)$, and $S = 8\lambda\tilde{W}(Y^+ + Y^-) \odot \left( \frac{[V]}{[\tilde{W}H]} H^T \right)$ with $Y = Y^+ - Y^- = (\tilde{W}^T \tilde{W} + \delta I)^{-1}$, $Y^+ = \max(Y, 0) \ge 0$ and $Y^- = \max(-Y, 0) \ge 0$, and $e_{F,N}$ is the $F$-by-$N$ matrix of all ones. As proved in Proposition 2, the constraint $W^*(\mu)^T e = e$ is satisfied for a unique $\mu$ in $(-\infty, t)$ where $t = \infty$ in this case. We can therefore use a Newton-Raphson method to find the $\mu_i$ with quadratic rate of convergence, see Proposition 3. Algorithm 2 summarizes our method to tackle (29).

## 4.2 Numerical experiments

In this section we compare baseline KL-NMF (that is, the standard MU), the min-vol KL-NMF from [25, Algorithm 1] that solves (29) using MU combined with line search (min-vol KL-NMF LS), and Algorithm 2 applied to the spectrogram of two monophonic piano sequences considered in [25]. The first audio sample is the first measure of “Mary had a little lamb”, a popular English song. The second audio sample corresponds to the first 30 seconds of “Prelude and Fugue No.1 in C major” from de Jean-Sebastien Bach played by Glenn Gould<fnref n="6" />. We use the following three setups:

* Setup §1: sample “Mary had a little lamb” with $K = 3$, 200 iterations.

* Setup §2: sample “Mary had a little lamb” with $K = 7$, 200 iterations.

<footnote n="6">⁶https://www.youtube.com/watch?v=ZlbK5r5mBH4</footnote>
**Algorithm 2** Min-vol KL-NMF

**Input:** A matrix $V \in \mathbb{R}^{F \times N}$, an initialization $H \in \mathbb{R}_+^{K \times N}$, an initialization $W \in \mathbb{R}^{F \times K}$, a factorization rank $K$, and a maximum number of iterations, maxiter, the parameters $\delta > 0$ and $\lambda > 0$.

**Output:** A min-vol rank-$K$ NMF $(W, H)$ of $V$ satisfying constraints in (29).

1: **for** $it = 1$ : maxiter **do**

2: % Update of matrix $H$

3: $H \leftarrow H \odot \frac{[W^T(\frac{[V]}{[WH]})]}{[W^T e_{F,N}]}$

4: % Update of matrix $W$

5: $Y \leftarrow (W^T W + \delta I)^{-1}$

6: $Y^+ \leftarrow \max(Y, 0)$

7: $Y^- \leftarrow \max(-Y, 0)$

8: $C \leftarrow e_{F,N} H^T - 4\lambda(WY^-)$

9: $S \leftarrow 8\lambda W (Y^+ + Y^-) \odot (\frac{[V]}{[WH]} H^T)$

10: $D \leftarrow 4\lambda W (Y^+ + Y^-)$

11: $\mu \leftarrow \text{root} (W^*(\mu)^T e = e)$ over $\mathbb{R}^K$ % see (32) for the expression of $W^*(\mu)$

12: $W \leftarrow W \odot \frac{[[C+e\mu^T]^2+S]^{\frac{1}{2}}-(C+e\mu^T)}{[D]}$

13: **end for**

• Setup §3: “Prelude and Fugue No.1 in C major” with $K = 16$, 300 iterations.

For each setup, the algorithms are run for the same 20 random initializations of $W$ and $H$. Table 3 reports the average and standard deviation of the runtime (in seconds) over these 20 runs. Table 4 reports the average and standard deviation of the final values for $\beta$-divergences (data fitting term) and the objective function of (29) over these 20 runs for min-vol KL-NMF LS and Algorithm 2. For this last comparison, the value for the penalty weight $\lambda$ has been chosen so that KL-NMF leads to reasonable solutions for $W$ and $H$. More precisely, the values for $\lambda$ are chosen so that the initial value of $\frac{\lambda|\logdet(W^{(0)T}W^{(0)}+\delta I)|}{D_\beta(V|WH)}$ is equal to 0.1, 0.1 and 0.022 for setup §1, setup §2 and setup §, respectively.

<table><caption>Table 3: Runtime performance in seconds of baseline KL-NMF, min-vol KL-NMF LS and Algorithm 2. The table reports the average and standard deviation over 20 random initializations.</caption><thead><tr><th rowspan="2">Algorithms</th><th colspan="3">runtime in seconds</th></tr><tr><th>setup §1</th><th>setup §2</th><th>setup §3</th></tr></thead><tbody><tr><td>baseline KL-NMF</td><td>0.53±0.03</td><td>0.45±0.02</td><td>4.32±0.30</td></tr><tr><td>min-vol KL-NMF LS [25]</td><td>3.79±0.13</td><td>2.39±0.30</td><td>10.19±1.28</td></tr><tr><td>Algorithm 2</td><td>0.58±0.03</td><td>0.66±0.03</td><td>4.80±0.38</td></tr></tbody></table>

We observe that the runtime of Algorithm 2 is close to the baseline KL-NMF algorithm which confirms the negligible cost of the Newton-Raphson steps to compute $\mu^*$ as discussed in Section 2.3. On the other hand, since no line search is needed, we have a drastic acceleration from 2x to 7x
<table><caption>Table 4: Final values for $D_β$ and the penalized objective $Ψ$ from (29) obtained with min-vol KL-NMF LS and Algorithm 2. The table reports the average and standard deviation over 20 random initializations for three experimental setups. A bold entry indicates the best value for each experiment.</caption><thead><tr><th></th><th></th><th>min-vol KL-NMF LS [25]</th><th>Algorithm 2</th></tr></thead><tbody><tr><td rowspan="2">setup #1</td><td>$D_{β,end}$</td><td>(3.52 ± 0.03)10<sup>3</sup></td><td>(<b>2.31</b> ± <b>0.01</b>)10<sup>3</sup></td></tr><tr><td>$Ψ_{end}$</td><td>(4.17 ± 0.03)10<sup>3</sup></td><td>(<b>3.08</b> ± <b>0.01</b>)10<sup>3</sup></td></tr><tr><td rowspan="2">setup #2</td><td>$D_{β,end}$</td><td>(3.54 ± 0.03)10<sup>3</sup></td><td>(<b>1.77</b> ± <b>0.02</b>)10<sup>3</sup></td></tr><tr><td>$Ψ_{end}$</td><td>(4.42 ± 0.04)10<sup>3</sup></td><td>(<b>2.87</b> ± <b>0.02</b>)10<sup>3</sup></td></tr><tr><td rowspan="2">setup #3</td><td>$D_{β,end}$</td><td>(7.77 ± 0.23)10<sup>3</sup></td><td>(<b>4.67</b> ± <b>0.08</b>)10<sup>3</sup></td></tr><tr><td>$Ψ_{end}$</td><td>(9.14 ± 0.20)10<sup>3</sup></td><td>(<b>6.50</b> ± <b>0.06</b>)10<sup>3</sup></td></tr></tbody></table>

compared to the backtracking line-search procedure integrated in min-vol KL-NMF LS [25]. Moreover, we observe in Table 4 that Algorithm 2 outperforms min-vol KL-NMF LS in terms of final values for the data fitting term and objective function values, with lower standard deviations.

# 5 Extension to quadratic disjoints constraints

Our general framework presented in Section 2 applies to $β$-NMF under disjoint linear equality constraints with penalty terms satisfying Assumption 1; see problem (2). We have showcased our approach on $β$-SSNMF in Section 3 and on min-vol KL-NMF under sum-to-one constraints on the columns of $W$ in Section 4. In this section, we show that the same framework can be extended to other simple constraints, namely disjoint quadratic constraints.

We consider sparse $β$-NMF for $β = 1$ where the rows of $H$ are penalized with the $ℓ_1$ norm and each column of $W$ have a fixed $ℓ_2$ norm. We show that MU satisfying the set of constraints can be derived which we apply on blind HU.

## 5.1 Problem formulation and algorithm

In this section we consider the following model involving quadratic disjoints constraints, that we refer to as hyperspheric-structured sparse $β$-NMF:

$$ \min_{W \in \mathbb{R}_+^{F \times K}, H \in \mathbb{R}_+^{K \times N}} D_{\beta}(V|WH) + \sum_{k=1}^{K} \lambda_k \|H(k, :)_1\|_1 \quad \text{such that} \quad e^T w_j^{(2)} = \rho \text{ for } 1 \le i \le K, \quad (33) $$

where $λ_k$ is a penalty weight to control the sparsity of the $k$-th row of $H$, and the quadratic constraints require the columns of $W$ to lie on the surface of a hyper-sphere centered at the origin with radius $\sqrt{\rho} > 0$. Without this normalization, the $ℓ_1$-norm regularization would make $H$ tends to zero and $W$ grows to infinity.

As done before, we update $W$ and $H$ alternatively. We tackle the subproblem in $H$ with $W$ fixed based on the MU developed in [11] and guaranteed to decrease the objective function:

$$ H = \tilde{H} \odot \frac{[W^T (V \odot [W\tilde{H}]^{(\beta-2)})]}{[W^T [W\tilde{H}]^{(\beta-1)} + \lambda e^T]}, \quad (34) $$
where $\lambda \in \mathbb{R}_+^K$ is the vector of penalty weights. It remains to compute an update for $W$. To do so, we use the convex separable auxiliary function $G$ from [12] constructed at the current iterate $\widetilde{W}$, from which we obtain, as before, the Lagrangian function

$$G^{\mu}(W|\widetilde{W}) = \sum_{f} G(w_{f}|\widetilde{w}_{f}) + \sum_{k} \lambda_{k} \|H(k, : )\|_{1} + \boldsymbol{\mu}^{T} \sum_{f} \left( w_{f}^{(2)} - \frac{1}{F} \rho e \right), \quad (35)$$

where $\boldsymbol{\mu} \in \mathbb{R}^K$ is the vector of Lagrange multipliers associated to the constraint $e^T W^{(2)} = \rho e^T$. Exactly as before (hence we omit the details here), given $\boldsymbol{\mu}$, one can obtain a closed-form solution:

$$W^{\star}(\boldsymbol{\mu}) = \frac{\left[ \left[ \left[ C \right]^{2} + 8 (e\boldsymbol{\mu}^{T}) \odot S \right]^{\frac{1}{2}} - C \right]}{[4e\boldsymbol{\mu}^{T}]}, \quad (36)$$

where $C = e_{F,N}H^T$ and $S = \widetilde{W} \odot \left( \frac{[V]}{[\widetilde{WH}]} H^T \right)$. Let us now write the expression of the quadratic constraint $\sum_f (W^*(\boldsymbol{\mu})_{f,i})^2 - \rho = 0$ for one specific column of $W$, say the $i$-th:

$$r_i(\mu_i) := \sum_f (W_{f,i}^*(\mu_i))^2 - \rho = \sum_f \left( \frac{\sqrt{(C_{f,i})^2 + 8\mu_i S_{f,i} - C_{f,i}}}{4\mu_i} \right)^2 - \rho = 0. \quad (37)$$

Computing the Lagrangian multiplier $\mu_i$ to satisfy the constraint requires computing the roots of the functions $r_i(\mu_i)$. We can show that each $W_{f,i}^*(\mu_i)$ (36) is a monotone decreasing, nonnegative convex function over $(0, +\infty)$. Therefore $\sum_f (W_{f,i}^*(\mu_i))^2$ is also monotone decreasing and convex in $\mu_i$ over $(0, +\infty)$. Indeed, let $g: \mathbb{R}_+ \to \mathbb{R}_+$ be a monotone decreasing, nonnegative convex function. If $g$ is twice-differentiable, then $(g^2)'' = 2(g')^2 + 2gg'' \ge 0$ since $g, g'' \ge 0$ and $(g^2)' = 2g'g \le 0$ since $g \ge 0$, $g' \le 0$ by hypothesis. Now we can conclude that $r_i(\mu_i)$ is a monotone decreasing convex function over $(0, +\infty)$. Moreover, using Hospital's rule, we have:

$$\lim_{\mu_i \to 0^+} \sum_f (W_{f,i}^*(\mu_i))^2 - \rho = +\infty \quad \text{and} \quad \lim_{\mu_i \to +\infty} \sum_f (W_{f,i}^*(\mu_i))^2 - \rho = -\rho < 0,$$

since $\rho > 0$. Therefore, the root of $r_i(\mu_i)$ is unique over $(0, +\infty)$. We use a Newton-Raphson method to solve the problem. Algorithm 3 summarizes our method.

## 5.2 Numerical experiments

In this section, we perform numerical experiments to evaluate the effectiveness of Algorithm 3 on the HU problem. To the best of our knowledge, sparse $\beta$-NMF$^7$ from [32] is the most recent algorithm that is able to tackle problem (33) for the KL-divergence by integrating the $\ell_2$-normalization for each update of matrix $W$. This approach is similar to that of [14] for $\beta$-SSNMF, that is, it uses parametrization, and resort to a heuristic with no guarantee on the decrease of the objective function. We refer to this algorithm as $\beta$-SNMF.

We apply Algorithm 3 and $\beta$-SNMF [32] to the three real hyperspectral datasets detailed in Section 3. This comparative study focuses on the convergence aspects including the evolution of the

<footnote>$^7$http://www.jonathanleroux.org/software/sparseNMF.zip</footnote>
**Algorithm 3** Hyperspheric-structured sparse KL-NMF

**Input:** A matrix $V \in \mathbb{R}^{F \times N}$, an initialization $H \in \mathbb{R}_+^{K \times N}$, an initialization $W \in \mathbb{R}^{F \times K}$, a factorization rank $K$, a maximum number of iterations, maxiter, a weight vector $\lambda > 0$.

**Output:** A sparse rank-$K$ NMF ($W, H$) of $V$ satisfying constraints in (33).

1: **for** $it = 1$ : maxiter **do**

2: % Update of matrix $H$

3: $H \leftarrow H \odot \frac{[W^T (V \odot [WH])^{(\beta-2)}]}{[W^T [WH]]^{(\beta-1)} + \lambda e^T}$

4: % Update of matrix $W$

5: $C \leftarrow e_{F,N} H^T$

6: $S \leftarrow W \odot (\frac{[V]}{[WH]} H^T)$

7: **for** $j = 1$ : K **do**

8: $\mu_i \leftarrow \text{root}(r_i(\mu_i))$ over $(0, +\infty)$ % See Equation (37)

9: **end for**

10: $W \leftarrow \frac{[[C]^2 + 8(e\mu^T) \odot S]^\frac{1}{2} - C]}{[4\mu e^T]}$

11: **end for**

objective function and the runtime; we refer the interested reader to the Supplementary Material S2 for qualitative result on the ability of sparse $\beta$-NMF to decompose such images. For all simulations, the algorithms are ran for 20 random initializations of $W$ and $H$, the entries of the penalty weight $\lambda$ has been set to 0.1, 0.05 and 0.05 for Samson, Jasper Ridge and Cuprite data sets, respectively. In order to fairly compare both algorithms, $\rho$ has been set to 1 as $\beta$-SNMF considers a $\ell_2$-normalization for the columns of $W$, and the entries of the weight vector $\lambda$ in Algorithm 3 have the same values as $\beta$-SNMF requires to use the same values for all rows of $H$. Table 5 reports the average and standard deviation of the runtime (in seconds) as the final value for the objective function over these 20 runs for a maximum of 300 iterations. Figure 1 displays the objective function values.

<table><caption>Table 5: Runtime performance in seconds and final value of objective function $\Phi_{\text{end}}(W, H)$ for Algorithm 3 and $\beta$-SNMF. The table reports the average and standard deviation over 20 random initializations with a maximum of 300 iterations for three hyperspectral data sets. A bold entry indicates the best value for each experiment.</caption><thead><tr><th rowspan="2">Algorithms</th><th colspan="2">Samson data set</th><th colspan="2">Jasper Ridge data set</th><th colspan="2">Cuprite data set</th></tr><tr><th>runtime (sec)</th><th>Φ<sub>end</sub>(W, H)</th><th>runtime (sec)</th><th>Φ<sub>end</sub>(W, H)</th><th>runtime (sec)</th><th>Φ<sub>end</sub>(W, H)</th></tr></thead><tbody><tr><td>Algorithm 3</td><td>11.07±0.19</td><td>(2.68±0.00)10<sup>3</sup></td><td>15.67±0.17</td><td>(4.65 ± 0.00)10<sup>3</sup></td><td>70.16 ± 0.85</td><td>(2.12 ± 0.00)10<sup>3</sup></td></tr><tr><td>β-SNMF [32]</td><td><b>7.63±0.13</b></td><td>(2.68±0.00)10<sup>3</sup></td><td><b>10.98±0.18</b></td><td>(4.71 ± 0.00)10<sup>3</sup></td><td><b>51.86 ± 0.74</b></td><td>(2.18 ± 0.00)10<sup>3</sup></td></tr></tbody></table>

According to Table 5 (top row), we observe that Algorithm 3 outperforms the heuristic from [32] in terms of final value for the objective functions while $\beta$-SNMF shows lower runtimes. Additionally, based on Figure 1, we observe that Algorithm 3 converges on average faster than $\beta$-SNMF for all the data sets, in terms of iterations. However, $\beta$-SNMF has a lower computational cost per iteration. Thus, we complete the comparison between both algorithms by imposing the same computational time: we run Algorithm 3 for 300 iterations, record the computational time and run $\beta$-SNMF for the
<figure><img src="image_2.png" /><figcaption>Figure 1: Averaged objective functions over 20 random initializations obtained for Algorithm 3 with 300 iterations (red line with circle markers), and the heuristic $\beta$-SNMF from [32] (black dashed line).</figcaption></figure>

same amount of time. Table 6 reports the average and standard deviation of the final value for the

<table><caption>Table 6: Final value of objective function values $\Phi_{end}(W, H)$ for Algorithm 3 and the heuristic from [32]. The table reports the average and standard deviation over 20 random initializations for an equal computational time that corresponds to 300 iterations of Algorithm 3. A bold entry indicates the best value for each experiment.</caption><thead><tr><th>Algorithms</th><th>Samson data set<br>$\Phi_{end}(W, H)$</th><th>Jasper Ridge data set<br>$\Phi_{end}(W, H)$</th><th>Cuprite data set<br>$\Phi_{end}(W, H)$</th></tr></thead><tbody><tr><td>Algorithm 3<br>$\beta$-SNMF [32]</td><td>(2.68±0.00)10<sup>3</sup><br>(2.68±0.00)10<sup>3</sup></td><td>(4.65 ± 0.00)10<sup>3</sup><br>(4.66 ± 0.00)10<sup>3</sup></td><td>(2.12 ± 0.00)10<sup>3</sup><br>(2.15 ± 0.00)10<sup>3</sup></td></tr></tbody></table>

objective function over 20 runs in this setting. Figure 1 (bottom row) displays the objective function w.r.t. time for the three data sets. On this comparison, Algorithm 3 and the heuristic from [32] perform similarly although Algorithm 3 has slightly better final objective function values. However, keep in mind that only Algorithm 3 is theoretically guaranteed to decrease the objective function.

# 6 Conclusion

In this paper we have presented a general framework to solve penalized $\beta$-NMF problems that integrates a set of disjoint constraints on the variables; see the general formulation (2). Using this framework, we showed that we can derive algorithms that compete favorably with the state of the art for a wide variety of $\beta$-NMF problems, such as the simplex-structured NMF and the minimum-volume
β-NMF with sum-to-one constraints on the columns of W. We have also shown how to extend the framework to non-linear disjoints constraints, with application to a sparse β-NMF model for β = 1 where each column of W lie on a hyper-sphere.

Further works will focus on the possible extension of the methods to non-disjoints constraints. The non-disjoint constraints will lead to roots finding problems of polynomial equations in the Lagrangian multipliers for which we hope to find conditions that ensure the uniqueness of the solution.

Another interesting direction of research would be to apply our framework to other NMF models. For example, in probabilistic latent semantic analysis/indexing (PLSA/PLSI), the model is the following: given a nonnegative matrix V such that $e^T V e = 1$ (this can be assumed w.l.o.g. by dividing the input matrix by $e^T V e$), solve

$$ \max_{W \ge 0, H \ge 0, s \ge 0} \sum_{f,n} v_{fn} \log(W \operatorname{Diag}(s)H)_{fn} \text{ such that } W^T e = e, He = e, s^T e = 1. $$

This model is equivalent to KL-NMF [9], with the additional constraint that $e^T W H e = e^T X e$, and hence our framework is applicable to PLSA/PLSI. Such constraints have also applications in soft clustering contexts; see [35].

## Acknowledgment

We would like to thank the Associate Editor and the reviewers for taking the time to carefully read the paper and for the useful feedback that helped us improve the paper. We also thank Arthur Marmin for identifying an error in our derivations when $\beta \in (1, 2)$ (indicated in red color in this version of the manuscript).

## A Convexity, concavity and complete monotonicity for a convex-concave decomposition of the discrete β-divergence

The discrete β-divergence can always be expressed as the sum of convex, concave, and constant terms. In Table 7 we introduce a convex-concave decomposition of the β-divergence which slightly differ from the one given in [12, Table 1] (by the fact that ours contains no constant term $\tilde{d}$) as given in Table 7.

<table><caption>Table 7: Proposed concave-convex decomposition of the discrete β-divergence.</caption><thead><tr><th>Decomposition<br>d<sub>β</sub> = d̃ + d̂</th><th>β ∈ (-∞, 1) \ {0}</th><th>β = 0</th><th>β = 1</th><th>β ∈ (1, 2)</th><th>β ∈ [2, +∞)</th></tr></thead><tbody><tr><td>d̃(x|y)</td><td>1/(1-β)x y<sup>β-1</sup></td><td>x/y</td><td>-x log y</td><td>(1/β)y<sup>β</sup> - (1/(β-1))xy<sup>β-1</sup></td><td>(1/β)y<sup>β</sup></td></tr><tr><td>d̂(x|y)</td><td>(1/β)y<sup>β</sup> - (1/(β(1-β)))x<sup>β</sup></td><td>log y/x - 1</td><td>y + x log x - x</td><td>(1/(β(β-1)))x<sup>β</sup></td><td>- (1/(β-1))x y<sup>β-1</sup> + (1/(β(β-1)))x<sup>β</sup></td></tr></tbody></table>

In Table 7, $y \in (0, \infty)$, $\beta$ is real valued and $x \in (0, \infty)$. Further, $\beta$ and $x$ are considered as parameters, $d_\beta$, $\hat{d}$ and $\tilde{d}$ being handled as univariate functions of $y$.

Let us now recall the definition of a complete monotonic function $f$:

**Definition 1.** A function $f$ is said to be completely monotonic (c.m.) on an interval $I$ if $f$ has derivatives of all orders on $I$ and $(-1)^n f^{(n)}(x) \ge 0$ for $x \in I$ and $n \ge 0$.
We can now introduce the properties of concavity, convexity and monotonicity for our convex-concave formulation of the discrete $\beta$-divergence:

**Proposition 4.** Given $\tilde{d}(\cdot|\cdot)$ and $\hat{d}(\cdot|\cdot)$ as defined above, we have that

1. $\tilde{d}(x|y)$ is $C^\infty$ and strictly convex on $(0, \infty)$ for $x > 0$ and $\beta \in \mathbb{R}$;

2. $\hat{d}(x|y)$ is concave for $x > 0$ and $\beta \in \mathbb{R}$;

3. for all $\beta < 2$, $\tilde{d}''(x|y)$ and $\hat{d}''(x|y)$ are c.m.

*Proof.* The proof is straightforward, given that $\tilde{d}(x|y)$ and $\hat{d}(x|y)$ linearly combine $C^\infty$ functions on $(0, \infty)$, and that in the same interval,

* $\log y$ is strictly concave;

* $y^\nu$ is strictly convex for all $\nu \in (-\infty, 0) \cup (1, \infty)$, and strictly concave for all $\nu \in (0, 1)$;

* $y^\nu$ is c.m. for all $\nu < 0$.

□

According to the first two items of Proposition 4, $\tilde{d}$ and $\hat{d}$ indeed yield a convex-concave decomposition of the $\beta$-divergence, which is a variant of [12, Table 1]. Let us remark that the successive minimization of an upper approximation of this convex-concave decomposition following the methodology presented in [12] yields to the usual multiplicative update scheme.

## References

[1] M. ABDOLALI AND N. GILLIS, *Simplex-structured matrix factorization: Sparsity-based identifiability and provably correct algorithms*, arXiv preprint arXiv:2007.11446, (2020).

[2] A. ANG AND N. GILLIS, *Accelerating nonnegative matrix factorization algorithms using extrapolation*, Neural computation, 31 (2019), pp. 417–439.

[3] D. P. BERTSEKAS, *Nonlinear Programming*, Athena Scientific, Belmont, MA, second ed., 1999.

[4] J. M. BIOCAS-DIAS, A. PLAZA, N. DOBIGEON, M. PARENTE, Q. DU, P. GADER, AND J. CHANUSSOT, *Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches*, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 5 (2012), pp. 354–379.

[5] E. C. CHI AND T. G. KOLDA, *On tensors, sparsity, and nonnegative factorizations*, SIAM Journal on Matrix Analysis and Applications, 33 (2012), pp. 1272–1299.

[6] A. CICHOCKI, R. ZDUNEK, AND S.-I. AMARI, *Hierarchical ALS algorithms for nonnegative matrix and 3D tensor factorization*, in Lecture Notes in Computer Science, Vol. 4666, Springer, 2007, pp. 169–176.

[7] A. CICHOCKI, R. ZDUNEK, A. H. PHAN, AND S.-I. AMARI, *Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-Way Data Analysis and Blind Source Separation*, John Wiley & Sons, 2009.

[8] O. DIKMEN, Z. YANG, AND E. OJA, *Learning the information divergence*, IEEE Transactions on Pattern Analysis and Machine Intelligence, 37 (2015), pp. 1442–1454.

[9] C. DING, T. LI, AND W. PENG, *On the equivalence between non-negative matrix factorization and probabilistic latent semantic indexing*, Computational Statistics & Data Analysis, 52 (2008), pp. 3913–3927.
[10] J. EGGERT AND E. KORNER, *Sparse coding and NMF*, in IEEE International Joint Conference on Neural Networks, vol. 4, 2004, pp. 2529–2533 vol.4.

[11] C. FÉVOTTE, N. BERTIN, AND J.-L. DURRIEU, *Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis*, Neural computation, 21 (2009), pp. 793–830.

[12] C. FÉVOTTE AND J. IDIER, *Algorithms for nonnegative matrix factorization with the β-divergence*, Neural computation, 23 (2011), pp. 2421–2456.

[13] X. FU, K. HUANG, N. D. SIDIROPOULOS, AND W.-K. MA, *Nonnegative matrix factorization for signal and data analytics: Identifiability, algorithms, and applications.*, IEEE Signal Process. Mag., 36 (2019), pp. 59–80.

[14] C. FÉVOTTE AND N. DOBIGEON, *Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization*, IEEE Transactions on Image Processing, 24 (2015), pp. 4810–4819.

[15] N. GILLIS, *The why and how of nonnegative matrix factorization*, in Regularization, Optimization, Kernels, and Support Vector Machines, J. Suykens, M. Signoretto, and A. Argyriou, eds., Machine Learning and Pattern Recognition, Chapman & Hall/CRC, Boca Raton, Florida, 2014, ch. 12, pp. 257–291.

[16] N. GILLIS, *Nonnegative Matrix Factorization*, SIAM, Philadelphia, 2020.

[17] N. GILLIS AND F. GLINEUR, *Accelerated multiplicative updates and hierarchical ALS algorithms for non-negative matrix factorization*, Neural computation, 24 (2012), pp. 1085–1105.

[18] N. GUAN, D. TAO, Z. LUO, AND B. YUAN, *NeNMF: An optimal gradient method for nonnegative matrix factorization*, IEEE Transactions on Signal Processing, 60 (2012), pp. 2882–2898.

[19] L. T. K. HIEN, N. GILLIS, AND P. PATRINOS, *Inertial block proximal methods for non-convex non-smooth optimization*, in International Conference on Machine Learning, 2020, pp. 5671–5681.

[20] D. HONG, T. G. KOLDA, AND J. A. DUERSCH, *Generalized canonical polyadic tensor decomposition*, SIAM Review, 62 (2020), pp. 133–163.

[21] P. HOYER, *Non-negative matrix factorization with sparseness constraints*, J. Mach. Learn. Res., 5 (2004), p. 1457–1469.

[22] J. KIM, Y. HE, AND H. PARK, *Algorithms for nonnegative matrix and tensor factorizations: A unified view based on block coordinate descent framework*, Journal of Global Optimization, 58 (2014), pp. 285–319.

[23] D. LEE AND H. SEUNG, *Algorithms for non-negative matrix factorization*, in Proceedings of the 13th International Conference on Neural Information Processing Systems, NIPS, MIT Press Cambridge, 2000, pp. 535–541.

[24] D. D. LEE AND H. S. SEUNG, *Learning the parts of objects by non-negative matrix factorization*, Nature, 401 (1999), p. 788.

[25] V. LEPLAT, N. GILLIS, AND A. M. S. ANG, *Blind audio source separation with minimum-volume beta-divergence NMF*, IEEE Transactions on Signal Processing, 68 (2020), pp. 3400–3410.

[26] C.-J. LIN, *Projected gradient methods for nonnegative matrix factorization*, Neural computation, 19 (2007), pp. 2756–2779.

[27] W.-K. MA, J. M. BIOUCAS-DIAS, T. CHAN, N. GILLIS, P. GADER, A. J. PLAZA, A. AMBIKAPATHI, AND C. CHI, *A signal processing perspective on hyperspectral unmixing: Insights from remote sensing*, IEEE Signal Processing Magazine, 31 (2014), pp. 67–81.

[28] K. MILLER AND G. SAMKO, *Completely monotonic functions*, Integral Transforms and Special Functions, 12 (2001), p. 389–402.

[29] K. NEYMEYR AND M. SAWALL, *On the set of solutions of the nonnegative matrix factorization problem*, SIAM Journal on Matrix Analysis and Applications, 39 (2018), pp. 1049–1069.
[30] J. ORTEGA AND W. RHEINBOLDT, *Iterative Solution of Nonlinear Equations in Several Variables*, Academic Press, New York, NY, 1970.

[31] Y. QIAN, S. JIA, J. ZHOU, AND A. ROBLES-KELLY, *Hyperspectral unmixing via $l_{1/2}$ sparsity-constrained nonnegative matrix factorization*, IEEE Transactions on Geoscience and Remote Sensing, 49 (2011), pp. 4282-4297.

[32] J. L. ROUX, F. J. WENINGER, AND J. R. HERSHEY, *Sparse NMF – half-baked or well done?*, tech. rep., Mitsubishi Electric Research Laboratories (MERL), 2015.

[33] Y. SUN, P. BABU, AND D. PALOMAR, *Majorization-minimization algorithms in signal processing, communications, and machine learning*, IEEE Transactions on Signal Processing, 65 (2017), pp. 794-816.

[34] T. VIRTANEN, *Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria*, IEEE Transactions on Audio, Speech, and Language Processing, 15 (2007), pp. 1066-1074.

[35] Z. YANG, J. CORANDER, AND E. OJA, *Low-rank doubly stochastic matrix decomposition for cluster analysis*, The Journal of Machine Learning Research, 17 (2016), pp. 6454-6478.

[36] A. L. YUILLE AND A. RANGARAJAN, *The concave-convex procedure*, Neural Computation, 15 (2003), pp. 915-936.

[37] F. ZHU, *Hyperspectral unmixing: ground truth labeling, datasets, benchmark performances and survey*, arXiv preprint arXiv:1708.05125, (2017).
# Supplementary Material

This supplementary materials provide additional numerical experiments. In S1, we show the evolution of the error as a function of the iterations for Algorithm 1 and GR-NMF [14] for the tests performed in Section 3. In S2, we provide qualitative results obtained with Algorithm 3 on hyperspectral images.

## S1. Evolution of the objective function for β-SSNMF

Figure 2 displays the evolution of the relative objective function values, that is, $\frac{D_{\beta}(V|WH)}{D_{\beta}(V|vee e^T)}$, of β-SSNMF for Algorithm 1 and GR-NMF [14] on the experiments described in Section 3. As mentioned in the paper, Algorithm 1 performs better than GR-NMF [14], except for $\beta = 0$.

## S2. Qualitative results obtained with Algorithm 3

In the following we report qualitative results obtained with Algorithm 3 applied to three HS real data sets, that are Samson, Jasper and Urban data sets. The first two data sets are detailed in Section 3. The Urban data set contains 162 spectral bands with 307×307 pixels with mostly six endmembers. Note that Cuprite data set is replaced by the Urban data set since endmembers for Cuprite correspond to chemical components which are more difficult to interpret visually while endmembers for Urban data sets are more easily interpretable.

As mentioned earlier, $\lambda_k$ enables to control sparsity of the k-th row of $H$. Given a row $H(k, :)$ $\in \mathbb{R}_+^N$ of $H$, a meaningful way to measure its sparsity is to consider the following measure [21]:

$$ \text{sp}(H(k, :)) = \frac{\sqrt{N} - \frac{\|H(k, :)\|_1}{\|H(k, :)\|_2}}{\sqrt{N} - 1} \in [0, 1]. \quad (38) $$

During the numerical experiments, we observed that Algorithm 3 gives better results when the initial values for $\lambda$ are low and progressively increased. During a specified interval of iterations $[it_{min}, it_{max}]$, the sparsity of the current iterate is measured by using equation (38), and the entries of $\lambda$ are dynamically updated (increased with a rate $\alpha > 1$) to achieve a desired sparsity level sp. The dynamic update of the weight vector to reach the desired levels of sparsity has been activated in the iterations intervals $[1, 150]$, $[1, 150]$ and $[1, 75]$ for Samson, Jasper and Urban, respectively. We report here the abundance maps of each end-member for two levels of average target sparsity that are 0.25 and 0.5. For all the simulations, the weight vector $\lambda$ has been initialized to $0.05e$, and the algorithm was run for 300 iterations.

We fix the number of endmembers to 3, 4 and 6 respectively for Samson, Jasper Ridge and Urban data sets, these values are commonly considered in the HS community [37]. Figures 3 to 5 picture the abundance estimation for the three data sets for the two levels of sparsity.

In order to validate the results obtained for the abundances of the endmembers, we display in Figures 6, 7 and 8 the ground truth results obtained in [37]. Note that the grayscale used in [37] is the complementary of the one used in Figures 3 to 5.

We observe that the abundance estimation gets significantly more accurate when the level of average sparsity is higher. For the Samson and Jasper Ridge data sets, the abundances for the endmembers are nicely estimated while five endmembers over six are well estimated for the Urban data set. The “Roof” is divided into “Roof1” and “Roof2/shadow” [31, 37]. In our simulations, it seems that the sixth endmember corresponds to some shadows with a small residual of “Grass”, while the “Roof” is not split into two groups.
<figure><img src="image_3.png" /><figcaption>Figure 2: Averaged relative objective function values over 20 random initializations obtained for Algorithm 1 (red line with circle markers) and the GR-NMF (black dashed line) applied to the three data sets detailed in the text for 300 iteration. The comparison is performed for different values of $\beta$, from top to bottom: $\beta = 2$, $\beta = 3/2$, and $\beta = 1$. Logarithmic scale for y axis.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Samson data set (a) and results ((b) and (c)) for the Abundance maps estimated using Algorithm 3 for the three endmembers: #1 Tree, #2 Soil and #3 Water. Two average sparsity levels considered: 0.25 (b) and 0.5 (c).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Jasper Ridge data set (a) and results ((b) and (c)) for the Abundance maps estimated using Algorithm 3 for the four endmembers: #1 Road, #2 Tree, #3 Water and #4 Soil. Two average sparsity levels are considered: 0.25 (b) and 0.5 (c).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Urban data set (a) and results ((b) and (c)) for the Abundance maps estimated using Algorithm 3 for the six endmembers: #1 Soil, #2 Tree, #3 Grass, #4 Roof, #5 Road/Asphalt and #6 Roof2/shadows. Two average sparsity levels are considered: 0.25 (b) and 0.5 (c).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Baseline abundances for the endmembers obtained for Samson data extracted from [37]: #1 Soil, #2 Tree and #3 Water.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Baseline abundances for the endmembers obtained for Jasper Ridge data extracted from [37]: #1 Road, #2 Soil, #3 Water and #4 Tree.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Baseline abundances for the endmembers obtained for Urban data extracted from [37]: #1 Asphalt, #2 Grass, #3 Tree, #4 Roof1, #5 Roof2/Shadow and #6 Soil.</figcaption></figure>