# Guaranteed Best Sparse Solutions for Spectral Unmixing

Nils Foix-Colonier, Sébastien Bourguignon, Frédéric Schmidt, Mehdi Latif

► To cite this version:

Nils Foix-Colonier, Sébastien Bourguignon, Frédéric Schmidt, Mehdi Latif. Guaranteed Best Sparse Solutions for Spectral Unmixing. 15th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), 2025, Barcelona, Spain. hal-05319871

HAL Id: hal-05319871

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

Submitted on 17 Oct 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
# GUARANTEED BEST SPARSE SOLUTIONS FOR SPECTRAL UNMIXING

Nils Foix-Colonier¹, Sébastien Bourguignon¹, Frédéric Schmidt²,³, Mehdi Latif¹

¹ Nantes Université, École Centrale Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France

² Université Paris-Saclay, CNRS, GEOPS Orsay, 91405, France

³ Institut Universitaire de France, Paris, France

## ABSTRACT

Cardinality-constrained sparse spectral unmixing can be solved using Branch-and-Bound algorithms, provided that the number of reference endmembers and the cardinality constraint are reasonably small. However, focusing solely on the best solution may not always be the most relevant approach, especially in the presence of high correlation between endmembers: solutions close to the optimal one—in terms of objective function—but with different supports (activated endmembers) may offer better interpretability. We propose a Branch-and-Bound algorithm that returns a given number of supports which are guaranteed to provide the lowest least-squares errors among all possible ones. The interest of such an approach is illustrated on numerical simulations, where chances for retrieving the ground truth among the few best solutions are shown to increase compared to considering only the optimal one, especially in the noisiest settings. A Python implementation is made available.

**Index Terms—** Spectral unmixing, sparsity, exact optimization, Branch-and-Bound

## 1. INTRODUCTION

Spectral unmixing of remote sensing data can be used to estimate the composition of scenes observed by a hyperspectral sensor. In planetary science, for instance, the unmixing of satellite hyperspectral data is a lever for characterizing the mineralogical composition of the Martian surface, otherwise difficult to access [1]. Numerous unmixing algorithms exist [2], and those that take into account sparsity constraints (i.e., the solution is composed of only a few spectra) are particularly relevant as they incorporate a physical constraint: in one pixel, i.e., over a limited area on the ground, only few minerals are observed, and therefore few spectra are mixed. Although sparse spectral unmixing has already been studied with approximate or heuristic methods (for example [3, 4, 5]), sparsity is rarely taken into account exactly, e.g., by imposing a maximum number of non-zero components. Indeed, although it can significantly improve solutions [6, 7], the $l_0$-"norm" constraint makes the problem difficult, and therefore potentially time-consuming to solve.

Branch-and-Bound methods were shown to enable exact optimization of cardinality-constrained linear least squares problems. First, a reformulation as a mixed-integer problem, solved by generic solvers, made it possible to obtain exact solutions, with a computational cost drastically reduced compared to exhaustive search [8]. An example of such implementation exists in the case of spectral unmixing, where non-negativity and sum-to-one constraints are taken into account [6]. Since then, it has been shown, for generic sparse optimization problems, that a dedicated Branch-and-Bound implementation is likely to strongly reduce the computation time [9]. An recent adaptation proposed a Branch-and-Bound algorithm tailored to sparse spectral unmixing [7].

Such $l_0$-norm solution, although optimal, may not be the most relevant to the user, as it may lack interpretability. In particular, in hyperspectral imaging, due to noise and possibly high correlation between the reference spectra, solutions obtained at neighbor locations may lack coherence. In this paper, we generalize the work in [7], by proposing a Branch-and-Bound algorithm which returns a set of solutions, corresponding to the set of supports (non-zero coefficients indices) that achieve the best cardinality-constrained least-squares fit.

The remainder of this paper is organized as follows. Section 2 introduces the problem of cardinal-constrained spectral unmixing and the associated notations. Section 3 describes the proposed algorithm in detail and provides a short pseudocode. Section 4 presents simulation results, along with a link to the code and data used. Finally, Section 5 concludes and discusses directions for future research.

## 2. CARDINALITY-CONSTRAINED UNMIXING

We consider that the reflectance spectrum measured in a pixel, $y \in \mathbb{R}^N$, is derived from a linear mixture of the spectra of the components present in the surface covered by the pixel [10]. The matrix of reference spectra (the dictionary) is denoted by $A \in \mathbb{R}^{N \times P}$ and $x \in \mathbb{R}^P$ is abundance vector, whose non-zero coefficients activate and weight the corresponding columns of $A$. Finally, an additive term $\varepsilon$ includes model errors and noise, yielding the model:

$$y = Ax + \varepsilon \quad (1)$$
Finding the abundance vector $\boldsymbol{x}$ from the measurement $\boldsymbol{y}$ and the dictionary $\mathbf{A}$ is an inverse problem, for which solutions can be found by maximizing the likelihood. Assuming Gaussian white noise, it amounts to minimizing the least-squares error between $\boldsymbol{y}$ and $\mathbf{A}\boldsymbol{x}$. This approach generally gives very poor results, due to the high correlation of the columns of $\mathbf{A}$ (if the dictionary is large, many spectra resemble each other). Physical constraints are therefore usually added on the abundance vector, in order to obtain more meaningful solutions Since the coefficients of $\boldsymbol{x}$ are relative areas of the corresponding endmembers, we can enforce non-negativity and sum-to-one constraints (the latter equals the $\ell_1$-norm due to the non-negativity). The resulting problem, often called FCLS (Fully Constrained Least Squares) in the literature [11], reads:

$$ \hat{\boldsymbol{x}}^{\text{FCLS}} = \arg \min_{\boldsymbol{x} \in [0,1]^P} \frac{1}{2} \|\boldsymbol{y} - \mathbf{A}\boldsymbol{x}\|_2^2 \quad \text{s.t.} \quad \|\boldsymbol{x}\|_1 = 1 \quad (2) $$

However, in practice, the solution to the FCLS problem (2), which has zeros due to non-negativity, tends to disperse the non-zero abundances over too many spectra of the dictionary, and may also fail in detecting some [3, 4, 6]. It is then interesting to take into account a stronger sparsity assumption, which rigorously constrains the $\ell_0$ “norm” of $\boldsymbol{x}$, that is, the number of non-zero coefficients, to $K \in \mathbb{N}^*$:

$$ \boldsymbol{x}^* = \arg \min_{\boldsymbol{x} \in [0,1]^P} \frac{1}{2} \|\boldsymbol{y} - \mathbf{A}\boldsymbol{x}\|_2^2 \quad \text{s.t.} \quad \|\boldsymbol{x}\|_1 = 1, \quad \|\boldsymbol{x}\|_0 \le K \quad (3) $$

Solving (3) in a guaranteed way is far more complex than (2), as the additional discrete constraint makes the problem combinatorial. Specific algorithms are needed to explore a non-continuous space, namely the space of solution supports. A Branch-and-Bound algorithm was recently proposed to that purpose in [7]. In the following section, we show that the Branch-and-Bound framework enables one to generalize this approach to the search for the set of best supports.

### 3. PROPOSED ALGORITHM

Let $\mathbb{S}^K$ be the set of supports with a cardinality limited to $K$:

$$ \mathbb{S}^K = \{S \subset [1, P] \mid |S| \le K\} \quad (4) $$

Let $\mathbf{A}_S$ denote the matrix $\mathbf{A}$ restricted to the columns with indices in $S$. We define the abundance vector solving the problem on a given support $S$:

$$ \boldsymbol{x}_S^* = \arg \min_{\boldsymbol{x}_S \in [0,1]^{|S|}} \frac{1}{2} \|\boldsymbol{y} - \mathbf{A}_S \boldsymbol{x}_S\|_2^2 \quad \text{s.t.} \quad \|\boldsymbol{x}_S\|_1 = 1 \quad (5) $$

Each support $S$ thus provides “its best solution”, measured by the value of the objective function at $\boldsymbol{x}_S^*$, denoted by $f^*$:

$$ f^*(S) = \frac{1}{2} \|\boldsymbol{y} - \mathbf{A}_S \boldsymbol{x}_S^*\|_2^2 \quad (6) $$

This allows for the comparison and ranking of the supports of $\mathbb{S}^K$: $f^*(S_1) \le f^*(S_2) \le \dots \le f^*(S_{|\mathbb{S}^K|})$. Hence, problem (3) comes down to finding the best sparse support:

$$ S_1 = \arg \min_{S \in \mathbb{S}^K} f^*(S) \quad (7) $$

In this work, we aim to obtain the $m$ best supports with at most $K$ components:

$$ \mathbb{S}^{*m} = \left\{ \begin{array}{l} S_i, i \in [1, m] \mid f^*(S_1) \le \dots \le f^*(S_m), \\ \forall S \in \mathbb{S}^K \setminus \mathbb{S}^{*m}, f^*(S) \ge f^*(S_m) \end{array} \right\} \quad (8) $$

To achieve this, the discrete search space (set of supports) is explored. The initial problem (equation (3), root node) is divided into two complementary sub-problems (left and right children), based on a decision: the nullity or not of a coefficient of $\boldsymbol{x}$, the $j$-th for example. On the one hand, we obtain a sub-problem (e.g., left child) where the $j$-th spectrum in $\mathbf{A}$ is included in the mixture—the abundance vector therefore has support $S = \{j\}$, with a cost of 1 for the $\ell_0$ constraint. On the other hand, we obtain a sub-problem (right child) where this spectrum is excluded from the solution (empty support for now). The other variables are undecided, and we denote the set of their (remaining) indices by $R = [1, P] \setminus \{j\}$. Repeating iteratively such division process, we create a binary tree of sub-problems, representing the set of possible supports. At a node $\mathcal{N}$, the cardinality-constrained sub-problem is:

$$ \boldsymbol{x}_{\mathcal{N}}^* = \left\{ \begin{array}{ll} \arg \min_{\boldsymbol{x}_{S \cup R} \in [0,1]^{|S \cup R|}} & \frac{1}{2} \|\boldsymbol{y} - \mathbf{A}_{S \cup R} \boldsymbol{x}_{S \cup R}\|_2^2 \\ \text{s.t.} \quad \|\boldsymbol{x}_{S \cup R}\|_1 = 1, \quad \|\boldsymbol{x}_R\|_0 \le K - |S| & \end{array} \right. \quad (9) $$

At each node, bounds are computed on the optimal value of the objective function (9). An upper bound $\text{UB}$ is computed by considering a particular solution. A lower bound $\text{LB}$ is computed, that should be smaller than the lowest value taken by the objective function over all the sub-problems derived from $\mathcal{N}$. Such bounding procedure is the key point which avoids exploring the entire tree, while preserving the mathematical guarantee of the optimality of $\mathbb{S}^{*m}$. Indeed, if the lower bound of a node is greater than the largest $\text{UB}$ value in the list of the $m$ best solutions (the worse solution temporarily accepted in $\mathbb{S}^{*m}$), then it is proved that the solutions for all descendant nodes will not provide a better solution, so the branch can be pruned.

To compute $\text{UB}$, we choose to impose $\boldsymbol{x}_R = \mathbf{0}$, which amounts to solving an FCLS problem on the support $S$:

$$ \text{UB} = f^*(S) \quad (10) $$

By design, through the abundance vectors corresponding to the minimizers associated with each support $S$, we enumerate the set of sparse candidates for problem (3). We then build a list of predefined maximum size $m$, which receives the solutions found as we explore the tree, updated in order to keep only the $m$ best ones.
Usually, the lower bound is derived from a convex relaxation of problem (9). One can show that, on $[0, 1]^{\mathbb{R}}$, the convex relaxation of the constraint $||x_R||_0 \le K - |S|$ is the $\ell_1$ constraint: $||x_R||_1 \le K - |S|$. We omit its proof due to lack of space. In our case, as the $\ell_1$ norm is set to 1, this $\ell_1$-norm relaxation becomes useless, except when $|S| = K$: then, the node is a leaf, that does not need to be divided: $x_R = 0$ and problem (9) is simply solved by the upper bound in equation (10). The best possible convex relaxation is therefore:

$$ LB = f^*(S \cup R) \qquad (11) $$

The tree exploration strategy, and the way nodes are divided into sub-problems, remain to be chosen. Let $x_{S \cup R}^*$ denote the minimizer of (11). We choose the branching index $j \in R$ with the largest value in $x_{S \cup R}^*$, which favors components more likely to be part of the optimal sets [9]. Thus, the left child—noted $\mathcal{N}^{(\ell)}$—corresponds to the sub-problem:

$$ S^{(\ell)} = S \cup \{j\} \quad \text{and} \quad R^{(\ell)} = R \setminus \{j\} \qquad (12) $$

Whereas for the right node, $\mathcal{N}^{(d)}$, we have:

$$ S^{(r)} = S \quad \text{and} \quad R^{(r)} = R \setminus \{j\} \qquad (13) $$

This node creation is called **divide** in Algorithm 1 pseudo-code, which summarizes the proposed method. By construction, the right node inherits its parent's upper bound. We therefore consider only left nodes as possible solutions.

<caption>**Algorithm 1:** BnB-m, returns the m best supports</caption>

**Input:** A, y, K, m
1 &nbsp; $\mathcal{N}$ = Node(S = $\emptyset$, R = [1, P])
2 &nbsp; $L$ = $\mathcal{N}$ # stack initialized with root
3 &nbsp; $L_{\text{sol}}$ = $\emptyset$ # the list yielding $\mathbb{S}^{*m}$
4 &nbsp; $\text{UB}^{*m} = +\infty$
5 &nbsp; **while** $L \neq \emptyset$ **do**
6 &nbsp;&nbsp;&nbsp;&nbsp; $\mathcal{N}$ = pop from $L$ # implicitly load S and R
7 &nbsp;&nbsp;&nbsp;&nbsp; LB = computed with equation (11)
8 &nbsp;&nbsp;&nbsp;&nbsp; **if** $\text{LB} \le \text{UB}^{*m}$ **then**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\text{UB}$ = computed with equation (10)
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $\text{UB} \le \text{UB}^{*m}$ and $\mathcal{N}$ is left node **then**
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $|L_{\text{sol}}| < m$ **then**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; add $\mathcal{N}$ to $L_{\text{sol}}$
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **else**
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; remove last element of $L_{\text{sol}}$
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; add $\mathcal{N}$ to $L_{\text{sol}}$
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; sort $L_{\text{sol}}$ on the UB of the nodes
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\text{UB}^{*m} = \text{UB}$ of the $m^{\text{th}}$ element of $L_{\text{sol}}$
18 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $|S| < K$ and $|R| > 0$ **then**
19 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\{\mathcal{N}^{(\ell)}, \mathcal{N}^{(r)}\} = \text{divide}(\mathcal{N})$
20 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; push $\mathcal{N}^{(r)}$ on L
21 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; push $\mathcal{N}^{(\ell)}$ on L
22 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; # sub-optimal branch cut otherwise
23 &nbsp; **return** $L_{\text{sol}}$

Created nodes are added to a stack, storing the left node last, in order to process it first. This corresponds to a depth-first exploration strategy which, in combination with the choice of branching variable, will prioritize the exploration of solutions that are assumed to be better.

Most of the computation time is devoted to determining the values for UB and LB at each explored node (equations (10) and (11), respectively). These problems, formulated as quadratic programs under positivity and sum-to-one constraints, are passed to the Clarabel [12] interior-point solver, which proved to be the most efficient one among the various open-source software packages evaluated.

## 4. SIMULATIONS AND RESULTS

The proposed algorithm is compared with other methods on synthetic data. A Python implementation, as well as the data used in this section, are available here: https://gitlab.univ-nantes.fr/ls2n-sims/bnb-whispers-2025. Linear mixtures $A\hat{x}$ were created using a dictionary A with $N = 113$ spectral bands and a maximum of 410 endmembers, taken from the United States Geological Survey [13] and adapted to simulate the OMEGA/Mars Express spectrometer observing Mars [1]. The experimental protocol consists in generating 100 unmixing problems per parameter triplet (K, P, SNR), for $K \in \{3, 4, 5\}$ mixed spectra, randomly selected from a dictionary of $P \in \{50, 100, 150, 200\}$ spectra composed of randomly selected columns of A, with a signal-to-noise ratio SNR $\in \{30, 35, 40, 50, 55, 60\}$ (dB), where SNR = $10 \log(\|A\hat{x}\|_2^2/(N\sigma^2))$. Abundances are drawn uniformly in $\{[0.05, 1]^K \|\|\hat{x}\|_1 = 1\}$, and Gaussian white noise is added, corresponding to the desired SNR. The SNR scale may seem unusual, but it is a consequence of the non-negativity of all quantities. Thus, 30 dB corresponds to fairly strong noise.

Each unmixing problem is solved using four sparse methods: “KFCLS” consists in taking the K largest abundances from the FCLS solution<fnref n="1" /> obtained by solving (2). “Backward” is a deflation heuristic [4], iteratively removing the smallest components in the FCLS solution until the desired level of sparsity is reached. “BnB” corresponds to the monosolution version of the proposed method: $m = 1$. “BnB-10” is the proposed method, which returns a set of $m = 10$ solutions. For a given support, amplitudes $\hat{x}$ minimizing the least-squares fit (6) are considered. Accuracy of a given solution $\hat{x}$ is measured using two metrics: SF (Support Found) equals 1 when the supports of $\hat{x}$ and $\hat{x}$ are equal, and 0 otherwise, and the SRE (Signal to Reconstruction Error ratio) is defined by:

$$ \text{SRE (dB)} = 10 \log_{10} \left( \frac{\|\hat{x}\|_2^2}{\|\hat{x} - \hat{x}'\|_2^2} \right) \qquad (14) $$

<footnote n="1">¹With sum-to-one $\hat{x}$, $\ell_1$-norm-penalized solutions (e.g. obtained by SuNSAL [3]) would achieve their best performance when the $\ell_1$ norm is set to 1, which is the standard FCLS solution (more details at the paper GitLab).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Percentage of success (SF) in finding the initial support, for various parameters (K, P, SNR).</figcaption></figure>

<table><caption>**Table 1.** Comparative performance of four methods, including the proposed algorithm, averaged over *K* and *P*.</caption><thead><tr><th>SNR</th><th>Metric</th><th>KFCLS</th><th>Greer</th><th>BnB</th><th>BnB-10<sup>(*)</sup></th></tr></thead><tbody><tr><td rowspan="2">30 dB</td><td>SF (%)</td><td>10.0</td><td>14.6</td><td>28.2</td><td>55.0</td></tr><tr><td>SRE (dB)</td><td>7.2</td><td>8.7</td><td>13.1</td><td>-</td></tr><tr><td rowspan="2">35 dB</td><td>SF (%)</td><td>19.2</td><td>27.4</td><td>52.4</td><td>77.8</td></tr><tr><td>SRE (dB)</td><td>11.8</td><td>14.2</td><td>22.0</td><td>-</td></tr><tr><td rowspan="2">40 dB</td><td>SF (%)</td><td>33.2</td><td>44.6</td><td>74.8</td><td>91.4</td></tr><tr><td>SRE (dB)</td><td>18.3</td><td>21.9</td><td>31.6</td><td>-</td></tr><tr><td rowspan="2">45 dB</td><td>SF (%)</td><td>48.8</td><td>61.8</td><td>88.0</td><td>97.5</td></tr><tr><td>SRE (dB)</td><td>26.0</td><td>30.7</td><td>39.5</td><td>-</td></tr><tr><td rowspan="2">50 dB</td><td>SF (%)</td><td>63.4</td><td>79.7</td><td>94.5</td><td>99.1</td></tr><tr><td>SRE (dB)</td><td>34.7</td><td>40.8</td><td>46.0</td><td>-</td></tr></tbody></table>

Table 1 compares performance metrics averaged over the 100 instances and over *K* and *P*. As expected, the global minimizer BnB always performs better than its competitors KFCLS and Backward. More importantly, it enhances the benefit of returning the *m* best solutions, which enables us to find a set containing the true support in cases where the global minimizer does not find it, and this is even more impactful as the SNR decreases. For example, for SNR = 30 dB, the true composition is found in only 28% of cases by BnB, the probability to find it among the 10 best supports, obtained with BnB-10, increases to 55%.

More detailed results are given in Figure 1 as a function of *K*, *P* and SNR. The interest brought by returning 10 solutions

<figure><img src="image_3.png" /><figcaption>Fig. 2. Percentage of success in having the initial support in $L_{sol}$, with $m \in [1, 10]$, on examples with $K = 5$ and $P = 50$.</figcaption></figure>

and not only the global optimizer is all the more significant as the problem complexity increases (SNR, but also *K* and *P*). This was rather expected, since the probability that the global optimizer retrieves the true support subsequently decreases. Note that for the most difficult problems ($K = 5, P = 200, SNR = 30$ dB), the true support is rarely found among the 10 best ones. In such cases, *m* should be substantially increased.

Figure 2 shows the probability that the true support is found by BnB-*m* as a function of *m* and SNR, for *K* = 5 and *P* = 50. For almost noise-free data (SNR = 50 dB), BnB already finds the true support in 89% of cases, while BnB-*m* finds it in 99% of the cases as soon as *m* = 3. As the SNR decreases, the gap increases, showing a regime where increasing *m* by a few units considerably increases the probability to retrieve the true support among the returned solutions.

<footnote>(*) BnB-10 method results are not straightforwardly comparable: the probability displayed is that the initial support is in the set of 10 solutions.</footnote>
<table><caption>**Table 2.** Resolution times (s) for BnB-10, averaged over 100 instances (number of instances requiring more than 1 hour).</caption><thead><tr><th>K</th><th colspan="3">3</th><th colspan="3">4</th><th colspan="3">5</th></tr><tr><th>SNR</th><th>30 dB</th><th>40 dB</th><th>50 dB</th><th>30 dB</th><th>40 dB</th><th>50 dB</th><th>30 dB</th><th>40 dB</th><th>50 dB</th></tr></thead><tbody><tr><th rowspan="4">P</th><td>50</td><td>5.1 (0)</td><td>4.5 (0)</td><td>4.5 (0)</td><td>7.2 (0)</td><td>6.3 (0)</td><td>5.6 (0)</td><td>9.0 (0)</td><td>9.7 (0)</td><td>8.6 (0)</td></tr><tr><td>100</td><td>39.7 (0)</td><td>40.4 (0)</td><td>38.9 (0)</td><td>111 (0)</td><td>69.9 (0)</td><td>35.8 (0)</td><td>303 (3)</td><td>141 (0)</td><td>76.7 (0)</td></tr><tr><td>150</td><td>149 (0)</td><td>60.2 (0)</td><td>78.7 (0)</td><td>470 (4)</td><td>339 (1)</td><td>333 (2)</td><td>935 (14)</td><td>567 (9)</td><td>409 (5)</td></tr><tr><td>200</td><td>294 (1)</td><td>237 (0)</td><td>176 (0)</td><td>748 (8)</td><td>417 (2)</td><td>336 (3)</td><td>1618 (26)</td><td>1136 (16)</td><td>808 (12)</td></tr></tbody></table>

Finally, Table 2 gives the computation times of BnB-10 for various triplets (*K*, *P*, SNR), measured on a single core with an Intel® Xeon® Gold 6126 CPU (computing times could be easily halved with more recent CPUs). More comparisons between BnB and BnB-10 computation times are available at the GitLab link. Note that the performance of the proposed method comes at a higher computational cost compared to BnB, and of course of more usual estimators, and this cost highly depends on the problem complexity. For example, for (*K* = 3, *P* = 50, SNR = 30 dB), the average resolution time is 1 second for BnB and 5.1 seconds for BnB-10, and for (*K* = 4, *P* = 100, SNR = 30 dB), these times increase to 59 and 111 seconds, respectively. Finally, for the hardest problems (large dictionaries and low SNR), a vast number of candidate solutions have close objective function values and are harder to discriminate. Consequently, computing times are highly dispersed, with instances exceeding the time limit of 1 hour. For the most complex cases (*K* = 5, *P* = 200, SNR = 30 dB) with BnB-10, this occurs in 26% of instances, and the average computing time is about 27 minutes.

## 5. CONCLUSION AND FURTHER WORKS

There is a clear interest in finding solution sets by exploiting the possibilities offered by the exploration of the support space with Branch-and-Bound algorithms. The same methodology could be developed for other sparsity-enhancing problems, e.g., without non-negativity and sum-to-one constraints. This paper did not discuss how such a solution set can be exploited in practice, which is left for future work. Expert interpretation could discriminate between the obtained solutions based on “external” knowledge that was not incorporated in the problem formulation. Most of all, a direct perspective of this work is the analysis of hyperspectral data cubes, where taking into account spatial regularity in estimated abundances should enable to search for consistent sparse mixtures while preserving global optimality properties.

## 6. REFERENCES

[1] Frédéric Schmidt, Gaspard Salomon, Sébastien Bourguignon, Joanna Gurgurewicz, and Daniel Mège, “Mineralogy analysis using linear unmixing under group constraint,” in 2023 13th

Workshop on Hyperspectral Imaging and Signal Processing: Evolution in Remote Sensing (WHISPERS), 2023, pp. 1–5.

[2] José M. Bioucas-Dias, Antonio Plaza, Nicolas Dobigeon, Mario Parente, Qian Du, Paul Gader, and Jocelyn Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches,” *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, vol. 5, no. 2, pp. 354–379, 2012.

[3] Marian-Daniel Iordache, José M. Bioucas-Dias, and Antonio Plaza, “Sparse Unmixing of Hyperspectral Data,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 49, no. 6, pp. 2014–2039, 2011.

[4] John B. Greer, “Sparse demixing of hyperspectral images,” *IEEE Transactions on Image Processing*, vol. 21, no. 1, pp. 219–228, 2012.

[5] Lucas Drumetz, Travis R. Meyer, Jocelyn Chanussot, Andrea L. Bertozzi, and Christian Jutten, “Hyperspectral image unmixing with endmember bundles and group sparsity inducing mixed norms,” *IEEE Transactions on Image Processing*, vol. 28, no. 7, pp. 3435–3450, 2019.

[6] Ramzi Ben Mhenni, Sébastien Bourguignon, Jordan Ninin, and Frédéric Schmidt, “Spectral Unmixing with Sparsity and Structuring Constraints,” in 2018 9th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), 2018.

[7] Mehdi Latif, Gwen Samain, Nils Foix-Colonier, and Sébastien Bourguignon, “Branch-and-bound algorithm for exact $\ell_0$-norm sparse spectral unmixing,” in EUSIPCO 2025, Palermo, Italy, 2025, EURASIP.

[8] Dimitris Bertsimas, Angela King, and Rahul Mazumder, “Best subset selection via a modern optimization lens,” *The Annals of Statistics*, vol. 44, no. 2, pp. 813–852, 2016.

[9] Ramzi Ben Mhenni, Sébastien Bourguignon, and Jordan Ninin, “Global optimization for sparse solution of least squares problems,” *Optimization Methods and Software*, vol. 37, no. 5, pp. 1740–1769, 2022.

[10] N. Keshava and J.F. Mustard, “Spectral unmixing,” *IEEE Signal Processing Magazine*, vol. 19, no. 1, pp. 44–57, 2002.

[11] Daniel C. Heinz and Chein-I-Chang, “Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 39, no. 3, pp. 529–545, 2001.

[12] Paul J. Goulart and Yuwen Chen, “Clarabel: An interior-point solver for conic programs with quadratic objectives,” 2024.

[13] Roger N. Clark et al., “USGS digital spectral library splib05a,” *US Geological Survey, Digital Data Series*, vol. 231, 2003.