Clémence Prévost

Valentin Leplat

C L Émence

DATA FUSION AND UNMIXING WITH THE REGULARIZED NON-NEGATIVE BLOCK-TERM DECOMPOSITION: JOINT PROBLEMS, BLIND APPROACH AND AUTOMATIC MODEL ORDER SELECTION

Keywords: Nonnegative tensor factorization, block-term decomposition, β-divergence, minimum volume regularization, automatic model order selection, blind spectral unmixing, hyperspectral super-resolution

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

In [START_REF] Kanatsoulis | Hyperspectral Super-Resolution: A Coupled Tensor Factorization Approach[END_REF][START_REF] Kanatsoulis | Hyperspectral Super-Resolution: Combining Low Rank Tensor and Matrix Structure[END_REF], the HSR problem was formulated as a coupled canonical polyadic decomposition, while a coupled Tucker decomposition was used in [START_REF] Prévost | Hyperspectral Super-Resolution with Coupled Tucker Approximation: Identifiability and SVD-based algorithms[END_REF]. However, the factors of these decompositions lacked physical interpretation, thus these methods could not perform unmixing.

The block-term decomposition (BTD) had been successfully used to perform unmixing on the SRI in [START_REF] Qian | Matrix-vector nonnegative tensor factorization for blind unmixing of hyperspectral imagery[END_REF]. Motivated by its usefulness, and by the success of tensorbased HSR, it has since then grown very popular to perform HSR [START_REF] Ding | Hyperspectral super-resolution via interpretable block-term tensor modeling[END_REF][START_REF] Guo | Multispectral and hyperspectral image fusion based on regularized coupled non-negative block-term tensor decomposition[END_REF][START_REF] Liu | Coupled tensor block term decomposition with superpixel-based graph laplacian regularization for hyperspectral super-resolution[END_REF][START_REF] Zhang | Hyperspectral super-resolution: A coupled nonnegative block-term tensor decomposition approach[END_REF]. In [START_REF] Prévost | Hyperspectral super-resolution accounting for spectral variability: Coupled tensor ll1-based recovery and blind unmixing of the unknown super-resolution image[END_REF], the BTD was used for joint fusion and unmixing in the presence of spectral variability.

Successfully solving both the HSR and unmixing problems faces several challenges.

First, the observation models and assorted algorithms must enforce non-negativity constraints. These constraints, when applied on the BTD factors, ensure that they can be recast in the linear mixing model as endmembers and abundance maps. The vast majority of algorithms for tensor/matrix non-negative factorization (NTF/NMF) rely on Block Coordinate Descent (BCD) schemes, see e.g., [START_REF] Bro | Least squares algorithms under unimodality and nonnegativity constraints[END_REF][START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF][START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF][START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF].

Second, in remote sensing, observed images often contain highly similar or nearly collinear materials. For instance, consider the scenario where mineral components within a rock sample or materials like road, soil, and sand are primarily composed of silica. This property of the images may lead unmixing algorithms to recover rankdeficient mixing factors. This is particularly challenging when using the BTD, since its numerical implementation is very sensitive to the initialization, as highlighted in [START_REF] Prévost | Hyperspectral super-resolution accounting for spectral variability: Coupled tensor ll1-based recovery and blind unmixing of the unknown super-resolution image[END_REF].

The third and most important challenge consists in correctly estimating the number of materials underlying the images. In practice, this value is often unknown.

An incorrect estimation may negatively affect the performance of the fusion and unmixing algorithms, and may also not ensure the uniqueness of the solutions to these problems. In the absence of groundtruth reference, the number of materials can be estimated prior to processing owing to the extensive spectral information available in the HSI [START_REF] Bioucas-Dias | Hyperspectral subspace identification[END_REF][START_REF] Fu | Self-dictionary sparse regression for hyperspectral unmixing: Greedy pursuit and pure pixel search are related[END_REF]. Some recent works were devoted to the task of estimating the number of factors in BTD models using e.g., penalized optimization [START_REF] Rontogiannis | Online rank-revealing block-term tensor decomposition[END_REF], Bayesian inference [START_REF] Giampouras | Block-term tensor decomposition model selection and computation: The bayesian way[END_REF] or autoregressive models [START_REF] Zarzoso | Parameter estimation in block term decomposition for noninvasive atrial fibrillation analysis[END_REF].

In this paper, a family of coupled optimization problems based on the BTD is introduced. The case of joint HSR and unmixing in remote sensing exemplifies the interest for these optimization problems. Other fields of applications could be envisioned for the considered model, such as audio signal processing [START_REF] Leplat | Blind audio source separation with minimum-volume beta-divergence nmf[END_REF], biomedical imaging [START_REF] Krieg | Diffusion tensor imaging fiber tracking using navigated brain stimulation-a feasibility study[END_REF] or graph signal processing [START_REF] Kanatsoulis | Tex-graph: Coupled tensor-matrix knowledgegraph embedding for covid-19 drug repurposing[END_REF]. Using β-divergences, various noise statistics present within the data can be accounted for, see [START_REF] Gillis | Nonnegative Matrix Factorization[END_REF] for a detailed overview of the topic. A family of simple, efficient and flexible algorithms is developed, that are able to handle the three above challenges within a unified framework.

Previously, the short paper [START_REF] Prévost | Nonnegative block-term decomposition with the β-divergence: Joint data fusion and blind spectral unmixing[END_REF] provided a brief overview of the optimization problem and partial simulations. This paper introduces a new family of regularized and constrained non-negative BTD problems and conduct additional simulations, including cases with collinear materials, unknown spatial degradation, and estimation of the model. Our main contributions are the following.

' We develop a family of tensor-based algorithms utilizing multiplicative updates tailored for the β-divergence. Among others, these algorithms allow the estimation of the degradation operators responsible for mapping the HSI and MSI to the unknown SRI. This fundamental principle, initially introduced in [START_REF] Leplat | Multi-resolution beta-divergence nmf for blind spectral unmixing[END_REF] through coupled matrix models, is extended into a coupled tensor framework.

' We incorporate a minimum-volume regularization into our optimization procedures. This regularization is designed to facilitate the computation of an interpretable and meaningful nonnegative BTD, while concurrently strengthening the identifiability of the decomposition model. It also addresses the challenge of rank deficiency and automating the critical step of determining This manuscript is for review purposes only.

the number of factors, commonly referred to as "model order selection" in the matrix case. In this paper, we use this terminology when discussing the ability of a method to identify the factor count within the decomposition. It is imperative not to conflate this concept with the "order" of the input tensor, defined as its number of dimensions. While the advantages of the minimumvolume regularization were previously discussed in [START_REF] Leplat | Blind audio source separation with minimum-volume beta-divergence nmf[END_REF], this previous work exclusively applied this regularization to uncoupled matrix models. ' To circumvent the numerical sensitivity of the BTD, a robust way to initialize the algorithms based on matrix multiplicative updates is introduced.

' The proposed methods compete favorably with the state-of-the-art for joint fusion and unmixing of synthetic and semi-real data sets including three noise statistics: Gaussian noise, Poisson noise and multiplicative Gamma noise.

Finally, we showcase the good performance of our methods in challenging cases: partially unknown degradation operators, almost collinear materials, and estimation of the number of endmembers. This paper is organized as follows. Subsection 2.1 introduces the background on tensor algebra, the observation model and its assumptions. Section 3 describes the optimization problems that we consider. Section 4 contains the proposed algorithms and the suggested initialization strategy. Section 5 presents the numerical experiments on a series of synthetic and semi-real datasets.

Background and problem formulation.

2.1. Background on tensor algebra. The following notations [START_REF] Comon | Tensors: A brief introduction[END_REF][START_REF] Kolda | Tensor Decompositions and Applications[END_REF] are used: lower (a) or uppercase (A) plain font for scalars, boldface lowercase (a) for vectors, boldface uppercase (A) for matrices and calligraphic (A) for tensors. The elements of vectors, matrices and tensors are denoted as a i , A i,j and A i1,...,i N , respectively.

The transpose of a matrix A is denoted by A T . Notation I N is used for the N ˆN identity matrix and 0 LˆK for the L ˆK matrix of zeros. Notation 1 L denotes an all-ones vector of size L ˆ1. For a matrix X, the notation X ě 0 means that X is entry-wise non-negative. Symbols b and d denote the Kronecker and Khatri-Rao products, respectively. The Hadamard (element-wise) product is denoted by d. The operator vec for the standard column-major vectorization of a matrix or a tensor.

Each dimension of a tensor is called a mode, and the number of dimensions is called order. Definition 2.1 introduces the BTD with ranks pL r , L r , 1q, that will be used to build the model. The main advantage of this decomposition is to link the low-rank factors to high-resolution abundance matrices and spectral signatures used in blind spectral unmixing of the unknown SRI. Definition 2.1 (Block-term decomposition). An order-3 tensor X P R IˆJˆK admits a block-term decomposition (BTD) with ranks pL r , L r , 1q

(L r L r 1-BTD) if X " R ÿ r"1 `Ar B T r ˘b c r , (2.1)
where b denotes the outer product, A r P R IˆLr , B r P R JˆLr , and c r P R K , for r P t1, . . . , Ru. Moreover, we denote A " rA 1 , . . . , A R s P R Iˆř r Lr , B " rB 1 , . . . , B R s P R Jˆř r Lr and C " rc 1 , . . . , c R s P R KˆR . Property 2.2 recalls the unfolding formulae for the L r L r 1-BTD, that will be helpful for building our algorithms. Property 2.2 (Tensor unfoldings). Using the above notation, the unfoldings of This manuscript is for review purposes only. a tensor X admitting an L r L r 1-BTD as above can be expressed as

X p1q " A pC d p Bq T , X p2q " B pC d p Aq T , X p3q " C rpA 1 d B 1 q 1 L1 , . . . , pA R d B R q 1 L R s T ,
where d p denotes the partition-wise Khatri-Rao products defined as follows:

Cd p A " rc 1 b A 1 , . . . , c R b A R s.
I 1 ă I 2 , J 1 ă J 2 and K 2 ă K 1 :
indices I , J denote the spatial dimensions whereas K denote the spectral ones ( " 1, 2). In order to ease the notation, we assume that

I " I 2 , J " J 2 et K " K 1 .
Assumption 2.3 (Structure of the SRI). In the noiseless case, the tensor Y admits a BTD with ranks (L r , L r , 1) such that

Y " R ÿ r"1 `Ar B T r ˘b c r . (2.2)
Under nonnegativity constraints, the vectors c r in Equation (2.2) can be viewed as the endmembers associated to the R constitutive materials of Y, while the matrices A r B T r " S r P R IˆJ represent the corresponding abundance maps. where A r P R IˆL and B r P R JˆL admit rank L r for all r P t1, . . . , Ru.

This assumption is reasonable, since an upper bound on the reconstruction error of such matrices by (2.3) can be obtained [START_REF] Bousse | A tensor-based method for large-scale blind source separation using segmentation[END_REF]. Furthermore, this assumptions will serve to link the linear mixing model to the L r L r 1-BTD.

Let S " rvectS 1 u, . . . , vectS R us P R IJˆR be the matrix containing the vectorized abundance maps of each material and C " rc 1 , . . . , c R s P R LˆR the matrix whose columns are the spectral signatures. The transposed third-mode unfolding of Equation (2.2) [START_REF] Ding | Hyperspectral super-resolution via interpretable block-term tensor modeling[END_REF][START_REF] Zhang | Hyperspectral super-resolution: A coupled nonnegative block-term tensor decomposition approach[END_REF] reads

Y p3qT " SC T P R IJˆK , (2.4)
which is the linear mixing model (LMM) for the SRI Y under nonnegativity constraints. Using Assumption 2.4 the block-term structure (2.2) can thus be viewed as tensor format for the LMM, under low-rank constraints of the abundance maps.

Observation model.

As done in previous works (see [START_REF] Wang | Tensor decompositions for hyperspectral data processing in remote sensing: A comprehensive review[END_REF] and references therein), the following model providing the links between Y and its two degraded versions Y 1 and Y 2 is considered. The tensors Y and (Y 1 ,Y 2 ) are such that

$ ' ' & ' ' % Y 1 « R ř r"1 ´P1 A r pP 2 B r q T ¯b c r , Y 2 « R ř r"1 `Ar B T r ˘b P 3 c r , (2.5)
This manuscript is for review purposes only.

which is a coupled L r L r 1-BTD. The tensors Y 1 and Y 2 are degraded from Y using linear downsampling operators P 1 P R I1ˆI , P 2 P R J1ˆJ et P 3 P R K2ˆK , assumed to be full-rank. In a remote sensing framework, the matrix P 3 P R K2ˆK contains the spectral response functions for each band of the MSI sensor. The spatial degradation matrices P 1 P R I1ˆI and P 2 P R J1ˆJ perform Gaussian blurring and downsampling along each spatial dimension, i.e. we suppose that the spatial degradation operation is separable, as in the commonly used Wald's protocol [START_REF] Wald | Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images[END_REF]. The approximately equal symbols in Equation (2.5) account for the presence of noise during the degradation process.

3. Optimization problems. Differently from those works, fusion between an HSI with an MSI with high spatial resolution allows one to seek for abundance maps at a higher spatial resolution. Furthermore, in remote sensing, Y 1 and Y 2 can be acquired at different time instants. The different acquisition conditions can result in e.g., variations in atmospheric, seasonal or illumination conditions [START_REF] Imbiriba | Generalized linear mixing model accounting for endmember variability[END_REF]. Moreover, the specificities of the imaging devices might change after launch of installation in an aircraft, due to e.g., outgassing, aging of components, or misalignment. These uncertainties motivate the need for more flexible models, capable of estimating one or several degradation matrices P i (i P t1, . . . , 3u).

Thus jointly solving the data fusion and blind unmixing problems consists in finding the pL r L r 1q factors tA r B T r u R r"1 , C, and possibly the degradation matrices P i , under the assumption of (2.5), subject to the constraints tA r B T r u R r"1 P X A,B , C P X C , P i P X P,i for i P t1, . . . , 3u , (3.2) where X T Ď R M ˆN denotes the feasible (convex) set for matrix T of size M ˆN .

Non-negative tensor optimization problem.

Trying to minimize the approximation errors in (2.5) leads, for instance, to minimizing the following general cost function:

(3.3) Φ "D β ˜Y1 } R ÿ r"1 ´P1 A r pP 2 B r q T ¯b c r ¸`λD β ˜Y2 } R ÿ r"1 `Ar B T r ˘b P 3 c r Ŗ ÿ r"1 γ pABq,r Ψ pABq,r pA r B T r q `γC Ψ C pCq `3 ÿ i"1 γ Pi Ψ Pi pP i q
subject to the constraints in (3.2). The scalar λ is a positive penalty parameter controlling the weights for the data fitting terms associated to each observation in the cost function Φ, while γ p.q are positive parameters controlling the weights of the regularization functions Ψp.q potentially applied to all the factors of the L r L r 1-BTD. Such regularization functions are used to promote solutions with desired and meaningful structures, such as low-rank, sparsity and minimum-volume to cite a few.

For a tensor Y P R IˆJˆK , the beta-divergence is defined as:

(3.4) D β ˜Y} R ÿ r"1 `Ar B T r ˘b c r ¸" ÿ i,j,k
d β `pYq i,j,k } `pA r q i,: pB r q T :,j ˘b pc r q k ˘, with d β px}yq the β-divergence between the two scalars x and y. For β " 2, this amounts to the standard squared Euclidean distance since d 2 px}yq " 1 2 px ´yq 2 . For β " 1 and β " 0, the β-divergence corresponds to the Kullback-Leibler (KL) divergence and the Itakura-Saito (IS) divergence, respectively. For NMF and NTF models, the data fitting term should be chosen depending on the noise statistic assumed in the generative model of the data, see [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music analysis[END_REF][START_REF] Févotte | Nonlinear hyperspectral unmixing with robust nonnegative matrix factorization[END_REF][START_REF] King | Optimal cost function and magnitude power for nmf-based speech separation and music interpolation[END_REF][START_REF] Leplat | Multi-resolution beta-divergence nmf for blind spectral unmixing[END_REF] and references therein for more details.

Section 4 presents our Algorithms to tackle specific instances of the general family of problems given in Equation (3.3). Nonnegative and non-regularized. For this case, γ p.q " 0 for the regularization functions Ψp.q, and the feasible sets X p.q in (3.2) are the nonnegative orthants of appropriate dimensions. The optimization problem defined here will be dubbed as "β-pL r , L r , 1q-NBTD".

Nonnegative minimum-volume. The following family of optimization problems is considered:

(3.5) min tSru,tPiu,C D β ˜Y1 } R ÿ r"1 ´P1 A r pP 2 B r q T ¯b c r ¸`λD β ˜Y2 } R ÿ r"1 `Ar B T r ˘b P 3 c r γvolpCq s.t. S r " A r B T r ě 0 @r, P i ě 0 @i, c r P ∆ K @r, where ∆ K " tx P R K `| ř K k"1 pxq k " 1u
, and volpCq is a function measuring the volume spanned by the columns of C. In [START_REF] Leplat | Blind audio source separation with minimum-volume beta-divergence nmf[END_REF], the authors use volpCq " log detpC T C `δIq where δ is a small positive constant that prevents log detpC T Cq from going to ´8

when C tends to a rank deficient matrix (that is, when rankpCq ă R). 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 [START_REF] Fu | Nonnegative matrix factorization for signal and data analytics: Identifiability, algorithms, and applications[END_REF]. Finally, the feasible sets for factors tA r B T r u R r"1 and P i are the nonnegative orthant of appropriate size, while the feasible set for factor C is the set of column-stochastic matrices. The problem defined in (3.5) will be dubbed as "min-vol β-pL r , L r , 1q-NBTD".

Section 4 explains how to tackle the two variants of optimization problems detailed above. In particular, the optimization problem (3.5) is handled with the general framework presented in [START_REF] Leplat | Multiplicative updates for nmf with β-divergences under disjoint equality constraints[END_REF] with a special focus on the Kullback-Leibler (KL) divergence, that is when β " 1.

4.

Algorithms. This section addresses the β-pL r , L r , 1q-NBTD problem. Subsequently, we will proceed to present our algorithmic solution for tackling the "min-vol KL-pL r , L r , 1q-NBTD" problem outlined in (3.5).

4.1.

Nonnegative and non-regularized. Most nonnegative tensor decomposition algorithms are based on an iterative scheme that alternatively updates one factor at the time with the others kept fixed, as it will be adopted in this paper.

The goal in this section is to derive an algorithm to solve the nonnegative and nonregularized version of (3.3) based on the multiplicative updates (MU) rule. Let us consider the subproblem in A (with the others fixed) after unfolding along the first mode following Property 2.2:

(4.1) min Aě0 D β pY p1q 1 }P 1 ApC d p P 2 Bq T q `λD β pY p1q 2 }ApP 3 C d p Bq T q.
To tackle this problem, we follow the standard majorization-minimization (MM) framework [START_REF] Sun | Majorization-minimization algorithms in signal processing, communications, and machine learning[END_REF] and the results given by [25, Lemma 2]. Given the current iterate

This manuscript is for review purposes only.

Ã, let us pose H 1 " pC d p P 2 Bq T and H 2 " pP 3 C d p Bq T , we obtain the following update:

(4.2)

A " Ã d ¨"P T 1 ˆ´P 1 ÃH 1 ¯.pβ´2q d Y p1q 1 ˙HT 1 `λ ˆ´Ã H 2 ¯.pβ´2q d Y p1q 2 ˙HT 2  " P T 1 ´P1 ÃH 1 ¯.pβ´1q H T 1 `λ ´ÃH 2 ¯.pβ´1q H T 2  ‹ ‹ ' .γpβq
where A d B (resp. rAs rBs ) is the Hadamard product (resp. division) between A and B, A p.αq is the element-wise α exponent of A, γpβq " 1 2´β for β ă 1, γpβq " 1 for β P r1, 2s and γpβq " 1 β´1 for β ą 2 [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the βdivergence[END_REF]. The subproblems in B and C can be solved similarly and their closed form expressions can be found in Appendix A.

As opposed to the majority of state-of-the-art methods, the proposed algorithms are also able to estimate the degradation matrices P i for i P t1, ..., 3u. These updates can be derived based on the classical MU associated to the matrix model X " UV T [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the βdivergence[END_REF]. For P 1 , we are interested in solving D β pY p1q 1 }P 1 ApC d p P 2 Bq T q. By posing V T " ApC d p P 2 Bq T , we derive:

(4.3) P 1 Ð r P 1 d ¨"´r P 1 V T ¯.pβ´2q d Y p1q 1  V " r P 1 V T ı .pβ´1q V ‹ ‹ ' .γpβq
. Similar rationale has been followed for the updates of P 2 and P 3 , see Appendix A for more details. Loop 2: All the factors, including the matrices P i , are alternatively updated. The maximum number of iterations for Loop 2 is i2. For the HSR problem, the operators P i for i P t1, ..., 3u are usually known and therefore the parameter i2 is set to zero.

Loop 2 is considered in the case we have partial knowledge or uncertainties on one of more downsampling operators, similarly as done in [START_REF] Leplat | Multi-resolution beta-divergence nmf for blind spectral unmixing[END_REF] with a matrix model. This case will be later referred to as "semi-blind".

The Algorithm is stopped when the relative change of the cost function Φ from (3.3) is below some given threshold κ, or when the maximum number of iterations is reached.

Initialization: Many options are available to initialize factors pA, B, Cq. In this paper, an efficient way to initialize the low-rank factors of the L r L r 1-BTD is proposed.

For factor C, VCA [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF] is performed on the HSI Y 1 , thus extracting high-resolution spectral information. Then, the matrix S of vectorized abundance maps is obtained by solving the following inverse problem under nonnegative constraints:

(4.4) min Sě0 P 3 CS T « Y p3q 2 .
Initialization of A r and B r for r P t1, . . . , Ru is performed based on the classical MU S r " A r B T r [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the βdivergence[END_REF] with a maximum of j1 iterations.

In the semi-blind case, i.e., when one or several matrices P i are unknown, they are initialized similarly using (4.3) with a maximum of j2 iterations. The initialization procedure is summarized in Algorithm 4.2.

This manuscript is for review purposes only.

Algorithm 4.1 MU for β-(L r ,L r ,1)-NBTD Require: Input tensors Y 1 ě 0, Y 2 ě 0, initializations A, B, C ě 0, downsampling operators P i ě 0 for i P t1, ..., 3u; R, ranks tL r u R r"1 , maximum number of iterations i1 and i2, a threshold 0 ă κ ! 1, and a weight λ ą 0. Ensure: An approximate solution to (3.3) under constraints (3.2)

1: % Loop 1: update of matrices A, B and C only 2: i Ð 0, Φ 0 " 1, Φ 1 " 0. Compute S using result of (4.4); 3: % Non-blind case: update of A r , B r only 4: j Ð 0, Φ 0 " 1, Φ 1 " 0. 5: while j ă j1 and ˇˇΦ j ´Φj`1 Φ j ˇˇą κ do 6:

Update A r , B r using the classical MU

7:

Compute the objective function Φ j`1 8: end while 9: % Blind case: Update A r , B r , and P i for i P t1, . . . , 3u 10: j Ð 0 11: while j ă j2 and ˇˇΦ

j ´Φj`1 Φ j ˇˇą κ do 12:
Update A r , B r , and P i using matrix-based MU

13:

Compute the objective function Φ j`1 14: end while It is majorized using a convex quadratic separable auxiliary function provided in [24, Eq. (3.6)] and which is derived as follows. First, the concave function log detpQq for Q ą 0 can be upper bounded using the first-order Taylor approximation: for any

r Q ą 0, log detpQq ď log detp r Qq `x r Q ´1, Q ´r Qy " x r Q ´1, Qy `cst,
where cst is some constant independent of Q. For any C, r C, and denoting r

Q " r C T r C `δI ą 0, we obtain (4.6) log detpC T C `δIq ď A r Q ´1, C T C E `cst " tracepC r Q ´1C T q `cst,
which is a convex quadratic and Lipschitz-smooth function in C, and where r C will be practically chosen as the current iterate for C. With this and following the framework from [START_REF] Leplat | Multiplicative updates for nmf with β-divergences under disjoint equality constraints[END_REF], the Lagrangian function is built as follows (4.7)

G µ `C| r C ˘" K ÿ k G Y1 pc k |r c k q `λ K ÿ k G Y2 pc k |r c k q `γ ˜K ÿ k l pc k |r c k q `cμ T K ÿ k ˆck ´1 K e ˙,
where c k denotes the k-th row of C, G is the separable majorizer for β-divergence proposed in [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the βdivergence[END_REF] (one for each β-divergence D β in (4.5)), l is a convex and separable majorizer 1 for the convex quadratic (4.

(4.8) C ‹ pµq " r C d " " " Q `eµ T ‰ .2 `Rı . 1 2 ´`Q `eµ T ˘ rDs , where Q " E K,I1J1 S T ´4γ `r CU ´˘`λP T 3 E K2,IJ H T , D " 4γ r C pU ``U ´q, and
R " 8γ r C pU ``U ´q d ˆ"Y p3q 1 ı r r CH1s H T 1 `λP T 3 " Y p3q 2 ı rP3 r CH2s H T 2 ˙with U " U `´U ´" 1 tight at the current iterate r C, that is ř K k l pr c k |r c k q " tracep r C r Q ´1 r C T q
This manuscript is for review purposes only.

`r C T r C `δI ˘´1 , U `" maxpU, 0q ě 0 and U ´" maxp´U, 0q ě 0, and E K,I1J1 is the K-by-I 1 J 1 matrix of all ones. As proved in [26, Proposition. 2], the constraint e T C ‹ pµq " e T is satisfied for a unique µ P R R . Therefore a Newton-Raphson method can be used to find the µ with quadratic rate of convergence; see [START_REF] Leplat | Multiplicative updates for nmf with β-divergences under disjoint equality constraints[END_REF]Proposition. 3].

Algorithm 4.3 summarizes our method to tackle (3.5).

Algorithm 4.3 MU for min-vol KL-(L r ,L r ,1)-NBTD Require: Input tensors Y 1 ě 0, Y 2 ě 0, initializations A, B, C ě 0, downsampling operators P i ě 0 for i P t1, ..., 3u; R, ranks tL r u R r"1 , maximum number of iterations i1 and i2, a threshold 0 ă κ ! 1, and weights λ, γ ą 0. Ensure: An approximate solution to (3.5)

1: % Loop 1 2: i Ð 0, Φ 0 " 1, Φ 1 " 0. 3: while i ă i1 and ˇˇΦ i ´Φi`1 Φ i ˇˇą κ do 4:
% Update of matrices A, B U `Ð max pU, 0q

U ´Ð max p´U, 0q 10: Comments on computational complexity: The computational complexity of Algorithm 4.1, designed to solve the β-(L r , L r , 1)-NBTD problem, exhibits asymptotic equivalence to the standard MU for β-NMF, after unfoldings. This implies that it requires OpIJK ˆřR r L r q operations per iteration. The main driver of this complexity lies in the aforementioned matrix products. One might question the influence of the computation step of Algorithm 4.3, which aims to calculate the optimal Lagrangian multipliers, thereby ensuring the column-stochasticity constraints of matrix C are satisfied. In [START_REF] Leplat | Multiplicative updates for nmf with β-divergences under disjoint equality constraints[END_REF], the authors demonstrate the quadratic convergence of the procedure This manuscript is for review purposes only.

Q Ð E K,I1J1 H T 1 ´4γ `CU ´˘`λP T 3 E K2,IJ H T 2 11: D Ð 4γC pU ``U ´q 12: R Ð 8γC pU ``U ´q d ˆ"Y p3q 1 ı rCH1s H T 1 `λP T 3 " Y p3q 2 ı rP3CH2s H T 2 13: µ Ð root `eT C ‹ pµq " e T ˘over R
employed to compute the Lagrangian multipliers. As a result, the number of iterations required for this process remains low in practical scenarios. Consequently, the overall complexity of Algorithm 4.3 is predominantly determined by the computationally intensive matrix products mentioned earlier. For more comprehensive information, we refer the reader to [START_REF] Leplat | Multiplicative updates for nmf with β-divergences under disjoint equality constraints[END_REF].

Parallelization: Some of the most computationally intensive steps of the proposed algorithm can be easily ran onto a parallel computation platform. Indeed, the complexity of the MU given in Equation (4.2), for instance, is mainly driven by the matrix products in which matrices A and H are involved. On MATLAB for example, one can easily take of advantage of a GPU compatible with CUDA libraries by simply transforming usual arrays into GPU arrays and significantly speed up the Algorithm.

Comments on convergence guarantees: In practice, the updates of factors A, B, C and P i are obtained from the element-wise maximum between the matrix updates, that correspond to the closed form expression of the minimizer of the majorization built at the current iterate [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the βdivergence[END_REF][START_REF] Leplat | Multi-resolution beta-divergence nmf for blind spectral unmixing[END_REF], and a small positive scalar (here we choose the Matlab machine epsilon). These modified updates aim at establishing convergence guarantee to stationary points within the Block Successive Minimization Methods (BSUM) framework [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF].

Numerical experiments.

All tests are preformed using Matlab R2021a on a laptop Intel CORE i7-11800H CPU @2.30GHz 16GB RAM with GeForce RTX3060

GPU. The code is available online at https://github.com/cprevost4/bLL1 NBTD.

Test setup.

The proposed algorithms were compared to several tensor methods designed for solving the HSR problem, namely STEREO and Blind-STEREO [START_REF] Kanatsoulis | Hyperspectral Super-Resolution: A Coupled Tensor Factorization Approach[END_REF][START_REF] Kanatsoulis | Hyperspectral super-resolution via coupled tensor factorization: Identifiability and algorithms[END_REF], SCOTT and BSCOTT [START_REF] Prévost | Coupled tensor low-rank multilinear approximation for hyperspectral super-resolution[END_REF], CT-STAR and CB-STAR [START_REF] Borsoi | Coupled tensor decomposition for hyperspectral and multispectral image fusion with inter-image variability[END_REF], SCLL1 [START_REF] Ding | Hyperspectral super-resolution via interpretable block-term tensor modeling[END_REF] and CNN-BTD-Var [START_REF] Prévost | Hyperspectral super-resolution accounting for spectral variability: Coupled tensor ll1-based recovery and blind unmixing of the unknown super-resolution image[END_REF]. Among them, SCLL1 and CNN-BTD-Var were based on the LL1-BTD model 2 , therefore they are able to solve the unmixing problem 3 .

Several matrix-based approaches were also benchmarked: CNMF [START_REF] Yokoya | Coupled Nonnegative Matrix Factorization Unmixing for Hyperspectral and Multispectral Data Fusion[END_REF], FUSE [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF],

HySure [START_REF] Simoes | A convex formulation for hyperspectral image superresolution via subspace-based regularization[END_REF], SFIM [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF] and MR-β-NMF [START_REF] Leplat | Multi-resolution beta-divergence nmf for blind spectral unmixing[END_REF]. Being based on coupled nonnegative matrix factorization, CNMF and MR-β-NMF were able to perform joint fusion and unmixing. In particular, MR-β-NMF was based on multiplicative updates with the β-divergence, but did not include a volume-regularizing constraint. The ranks and regularization parameters were chosen according to the original works.

The proposed algorithms were initialized using Algorithm 4.2 with a maximum of 500 iterations and a threshold κ " 10 ´7. All methods were limited to a maximum of 1000 iterations, encompassing both loops, and utilized a fixed value of κ " 10 ´7.

The groundtruth SRI Y was compared to the estimated SRI p Y obtained by the algorithms. The main performance metric used in comparisons was the Peak Signalto-Noise ratio (PSNR) [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF]:

(5.1) PSNR " 10log 10 ˜}Y} 2 F } p Y ´Y} 2 F ¸.
In addition to PSNR, we considered different metrics [START_REF] Yokoya | Hyperspectral and multispectral data fusion: A comparative review of the recent literature[END_REF] described below:

(5.2) CC " 1 IJK ˜K ÿ k"1 ρ ´Y:,:,k , p Y :,:,k ¯¸,
where ρp¨, ¨q is the Pearson correlation coefficient between the estimated and original

2 These works all considered the simpler model with L 1 " . . . " L R . 3 The performance of SCLL1 for blind spectral unmixing was not assessed in the original work.

spectral slices;

(5.3) ERGAS " 100 d

g f f e 1 IJK K ÿ k"1 } p Y :,:,k ´Y:,:,k } 2 F µ 2 k ,
where µ 2 k is the mean value of p Y :,:,k . ERGAS represents the relative dimensionless global error between the SRI and the estimate, which is the root mean-square error averaged by the size of the SRI. We also used Spectral Angle Distance (SAD):

(5.4) SAD " 1 R R ÿ r"1 arccos ˆcT r p c r }c r } 2 }p c r } 2 ˙,
which computes the spectral angle distance between original and estimated spectra, and can be used to assess unmixing performance as well. Performance for recovery of the abundance maps was assessed using the root mean-squared error between reference S and estimate p S:

(5.5) RMSE " 1 R R ÿ r"1 g f f e 1 IJ IJ ÿ d"1 ´pS r q d ´pp S r q d q ¯2.
The computational time for each algorithm was given by the tic and toc functions of Matlab. Finally, the compression rate (C.R.) of each method was computed, i.e., the ratio between the number of elements in the unknown SRI and the number of unknown parameters in the algorithms. The higher this value, the less parameters need to be recovered by the model.

Degradation model.

We considered synthetic and semi-real datasets for which a reference SRI is available. The HSI was obtained by spatial degradation of Y using P 1 and P 2 while the MSI was obtained by spectral degradation of Y with P 3 according to model (2.5).

For spatial degradation, we followed the commonly used Wald's protocol [START_REF] Wald | Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images[END_REF].

The matrices P 1 , P 2 were computed with a separable Gaussian blurring kernel of size q " 9. Downsampling was performed along each spatial dimension with a ratio d " 4

between the SRI and HSI, as in previous works [START_REF] Kanatsoulis | Hyperspectral Super-Resolution: A Coupled Tensor Factorization Approach[END_REF]- [START_REF] Ding | Hyperspectral super-resolution via interpretable block-term tensor modeling[END_REF].

For the spectral degradation matrix P 3 , the spectral response functions of the Sentinel-2 instrument4 was used. It spans the electromagnetic spectrum from 412nm to 2022nm and produces a 10-band MSI (K 2 " 10) corresponding to the wavelengths 433-453nm (atmospheric correction), 458-522nm (soil, vegetation), 543-577nm (green peak), 650-680nm (maximum chlorophyll absorption), 698-712nm (red edge), 733-747nm (red edge), 773-793nm (leaf area index, edge of NIR), 785-900nm (leaf area index), 855-875nm (NIR plateau), 935-955nm (water vapour absorption). The spectral degradation matrix P 3 acts as a selection-weighting matrix, utilizing the shared spectral bands between the SRI and the MSI datasets.

As done in [START_REF] Leplat | Multi-resolution beta-divergence nmf for blind spectral unmixing[END_REF], three noise statistics were considered: Gaussian Noise, Poisson noise and Multiplicative Gamma noise. Gaussian and Poisson noise were added to the observations in order to yield 30dB SNR. For Gamma noise, we considered a distribution of mean 1 and variance 0.05. Since both observations are subject to the same level of noise, we fixed the weight λ " 1.

Datasets.

Experiments were run on synthetic and semi-real datasets described below.

Synthetic dataset: A dataset where the SRI admits an exact L r L r 1-BTD was created, in order to perform simulations in a controlled environment.

We considered R " 4 spectral signatures c r (r P t1, . . . , Ru) obtained from the Jasper Ridge reference data 5 . The SRI Y P R IˆJˆK (I " J " 120, K " 173) was split into 36 equal blocks in the spatial dimensions. We set L 1 " L 2 " 3 and L 3 " L 4 " 6.

Each abundance map S r (r P t1, . . . , Ru) was a block matrix with L r blocks of size I Lr ˆJ Lr . Each block in the parcel map was a patch composed of entries equal to one. At most one material is active in each block, thus the so-called "pure pixel assumption" was valid. Figure 1 depicts the abundance maps for the four materials in the synthetic data set. The reference SRI is computed as follows:

Y " R ÿ r"1 S r b c r .
Semi-real datasets: The first semi-real dataset we considered is built upon the Jasper Ridge reference SRI Y P R 100ˆ100ˆ173 . It includes four materials: road, soil, water and vegetation. We chose R " 4, and L 1 " 15, L 2 " 8, L 3 " 20, L 4 " 13.

A second semi-real dataset was considered and based on the Ivanpah Playa [START_REF] Mahalingam | Reflectance based vicarious calibration of hysis sensors and spectral stability study over pseudo-invariant sites[END_REF] reference SRI Y P R 95ˆ95ˆ156 This dataset is composed of four, very similar materials: road, solar panels, light sand and darker sand. Thus it constitutes a challenging example for which unmixing may lead to rank-deficient solutions.

Reconstruction results.

In order to assess the performance of our approaches for data fusion, the quality metrics obtained with each method over 5 trials are reported in the Tables 123 The proposed algorithms yield better metrics than most baseline methods. Moreover, they provide excellent reconstruction results for, e.g., the Ivanpah Playa dataset.

The compression rate is higher than that of matrix-based approaches (i.e., less information was needed to reconstruct the SRI). It is however lower than that of other tensor-based algorithms, but provide significantly better results than, e.g., CT-STAR and CNN-BTD-Var for the semi-real datasets.

As it will be highlighted in the following sections, the good unmixing performance of the minimim-volume regularized approach come at the cost of a tradeoff in the reconstruction performance. The proposed algorithms however remain very competitive for reconstruction of the three considered SRI. This manuscript is for review purposes only. One can observe that the baseline algorithm SCLL1 provide rather poor results, due to its incapacity to account for different noise statistics and different L r values.

Conversely, the spectral signatures and abundance maps are correctly estimated by the proposed algorithms for the three considered datasets. Algorithm 4.3 generally yields better estimation of the spectra and less artifacts on the abundance maps. Even though Ivanpah Playa is a notoriously difficult dataset for unmixing, Algorithm 4.3 provides a good estimation of the materials in the scene. 5.6. Performance in a semi-blind scenario. The performance of our algorithms was then evaluated in a semi-blind scenario, i.e., when the spatial degradation matrices P 1 and P 2 are unknown. This scenario is likely to occur, e.g., when the HSI and MSI are acquired at different times or in the presence of motion blur. The proposed algorithms were compared to two semi-blind tensor algorithms: BSTEREO and BSCOTT, that do not estimate P 1 and P 2 .

We considered the Jasper Ridge dataset and the three noise statistics. Table 4 presents the reconstruction metrics obtained for the benchmarked algorithms. The proposed algorithms generally provides good reconstruction, especially with additive Gaussian noise. One can observe that using Algorithm 4.3 with β " 1 improves the overall reconstruction performance. Although the computation time is higher than that of the other tensor algorithms, it is significantly smaller than the one of the MR-β-NMF approach.

Figure 5 shows a portion of the reference and estimated spatial degradation (that This manuscript is for review purposes only. This test considered the synthetic dataset with 4 endmembers. We deliberately chose R " 6 for our algorithms to simulate an overestimation of the number of materials. The first 4 columns of C were initialized using VCA while the remaining two were initialized randomly. We chose L 1 " L 2 " L 5 " L 6 " 3 and L 3 " L 4 " 6. new Algorithm, dubbed as "min-vol β-(L r ,L r ,1)-NBTD", was introduced to efficiently solve the regularized problems. The experiments show that the minimum-volume regularization approach improves the unmixing performance, at the cost of a small loss in the reconstruction performance. More importantly, the "min-vol β-(L r ,L r ,1)-NBTD" algorithm correctly estimates the number of materials in the datasets, additionally to perform high-quality unmixing. The minimum-volume regularization approach also previously showed its interest for automatic model order selection in audio signal processing. Future works will be devoted to the mathematical analysis of this mechanism. , with H 2 " rpA 1 d B 1 q 1 L1 , . . . , pA R d B R q 1 L R s T and H 1 " pP 2 b P 1 qH 2 .

The updates for P 2 (resp. P 3 ) are obtained by substituting P 1 by P 2 (resp.

This manuscript is for review purposes only.

Assumption 2 . 4 (

 24 Structure of S r ). Matrices S r are low-rank matrices, i.e., S r « A r B T r P R IˆJ , (2.3)

3. 3 .

 3 Considered variants of the problem. In this paper, two specific instances of the optimization problem (3.3) subject to constraints in (3.2) are considered.

Algorithm 4 .

 4 1 summarizes the proposed method to tackle the β-(L r ,L r ,1)-NBTD optimization problem (3.3). It consists in two optimization loops: Loop 1: A, B and C only are alternatively updated with downsampling matrices fixed for a maximum of i1 iterations. P i for i P t1, ..., 3u kept fixed to obtain good estimates for A, B and C.

15 : 1 }CH 1 q `λD β pY p3q 2 }P 3

 151123 return A " rA 1 , . . . , A R s, B " rB 1 , . . . , B R s and C 4.2. Nonnegative minimum-volume. The primary approach for addressing the "min-vol KL-pL r , L r , 1q-NBTD" problem in equation (3.5) is analogous to the methodology employed in Section 4.1, specifically the cyclic block majorization minimization (BMM) framework. This framework involves cyclically updating each block of variables with the others kept fixed, with each update accomplished by minimizing a majorizer constructed using the current iterate. The MU for updating the factors This manuscript is for review purposes only. tA r , B r u R r"1 and tP i u 3 i"1 for Problem (3.5) are the same as the ones presented in the previous section dedicated to tackling the β-(L r ,L r ,1)-NBTD problem, since the subproblems in theses factors are identical. Therefore, this subsection focuses on solving the sub-problem of (3.5) that concerns the minimum-volume penalty and columnstochasticity constraints, specifically the updating factor C. The recent optimization framework proposed in [26] is used to derive efficient MU for the sub-problem in C, which is defined as follows: CH 2 q `γ log detpC T C `δIq s.t. e T C " e T , C ě 0. with H 1 " pP 1 A d vec P 2 Bq T , and H 2 " pA d vec Bq T . An upper bound for the term log detpC T C`δIq is required by [26, Assumption. 1].

  6) and given by[START_REF] Leplat | Blind audio source separation with minimum-volume beta-divergence nmf[END_REF] Eq. (3.6)], and c is a constant. Let µ be the vector of Lagrange multipliers of dimension R associated to each linear constraint e T c r " 1. One can easily show that G µ is separable w.r.t. each component of c k and, given µ, one can compute the closed-form solution (the details are omitted):

Fig. 1 .

 1 Fig. 1. Abundance maps for the four materials in the synthetic dataset.

  below. The two best metrics of each columns are shown in bold. For comparison with Algorithm 4.3, we only show the results for the KL divergence, i.e., for β " 1 along with Poisson noise. The results for other noise statistics are available in supplementary materials.

  only the results for our algorithms and SCLL1 in the scenario of Poisson additive noise are shown. The other results, including the other baseline algorithms and noise statistics, are available as supplementary materials. Figures2-4show the reference and estimated spectra (first column), reference abundance maps (second column) and the estimated maps (third column).

Fig. 2 .

 2 Fig. 2. Unmixing on the synthetic dataset with Poisson noise, (a) Algorithm 4.1, (b) Algorithm 4.3, (c) SCLL1.

  MethodCC (Ò) SAD (Ó) RMSE (Ó) ERGAS (Ó) PSNR (dB) (Ò) Time (sec) (Ó) C.R.

Fig. 3 .

 3 Fig. 3. Unmixing on the Jasper Ridge dataset with Poisson noise, (a) Algorithm 4.1, (b) Algorithm 4.3, (c) SCLL1.

Fig. 4 .

 4 Fig. 4. Unmixing on the Ivanpah Playa dataset with Poisson noise, (a) Algorithm 4.1, (b) Algorithm 4.3, (c) SCLL1.

  is, P 2 b P 1 ). A single slice of the true HSI and the tensor constructed by mode product of the SRI with the estimated spatial degradation is also shown. The proposed algorithms recovers correctly the structure of the spatial degradation operators. Moreover, the HSI obtained from the estimated degradation is coherent with the reference HSI, and Algorithm 4.3 provides slightly better results.

Figure 6

 6 Figure 6 shows our unmixing results in the case β " 1. The spectra and abundance maps are correctly estimated by the proposed algorithms, even without prior knowledge on the spatial degradation matrix. Resorting to Algorithm 4.3 offers better estimation of the spectral signatures and less artifacts due to the low-rank assumption on the abundance maps.

5. 7 .

 7 Estimation of the number of endmembers. This subsection highlights the capabilities of Algorithm 4.3 to accurately retrieve the number of endmembers underlying an image.

Figure 7

 7 Figure 7 presents the unmixing results obtained with Algorithm 4.1 and Algo-

Fig. 5 .

 5 Fig. 5. Estimation of the spatial degradation and comparison with the real HSI, (a) Algorithm 4.1 (β " 1), (b) Algorithm 4.3 (β " 1).

Fig. 6 .

 6 Fig. 6. Semi-blind unmixing, (a) Algorithm 4.1 (β " 1), (b) Algorithm 4.3 (β " 1).

rithm 4 .

 4 3 with β " 1. Algorithm 4.1 only recovers the first three spectra correctly, while the other three do not match the reference spectral signatures. The estimated abundance maps seem like mixtures of several reference abundance maps, thus they are not estimated correctly. Algorithm 4.3 recovers with high accuracy the four spectral signatures and abundance maps of interest. The other two materials are almost constant and show very low magnitude. These numerical evidences support the capacity of Algorithm 4.3 to estimate the correct number of materials, hence performing automatic model order selection, on top of the reconstruction and unmixing tasks.

6 .

 6 Conclusion. This paper proposes a family of coupled tensor-based optimization problems. Using the beta-divergence, various noise statistics can be accounted for. The use of the BTD allows the proposed algorithms to jointly solve the reconstruction and unmixing tasks in remote sensing, even in challenging cases such as a semi-blind scenario or almost collinear materials. We introduced a new family of penalized optimization problems with a focus on the minimum-volume regularization. A This manuscript is for review purposes only.

Fig. 7 .

 7 Fig. 7. Unmixing results with overestimation of the rank R, (a) Algorithm 4.1 (β " 1), (b) Algorithm 4.3 (β " 1).

  Appendix A. Detailed MU updates. This appendix contains the closedform expressions for the updates of B and C, similarly to (4.2). (A.1) B " B d ¨"P

P 3 2 )

 32 and defining V T " BpC d p P 1 Aq T (resp. V T " CpA d vec Bq T ) in (4.3).

  2.2. Assumptions. Let us consider two tensors Y 1 P R I1ˆJ1ˆK1 and Y 2 P R I2ˆJ2ˆK2 . The observed tensors Y 1 , Y 2 are degraded versions of the same tensor Y P R IˆJˆK . In remote sensing, the tensors Y 1 and Y 2 respectively denote the HSI and MSI, and Y denotes the unknown SRI we intend to recover. Thus we assume

  3: while i ă i1 and ˇˇΦ

			i ´Φi`1 Φ i	ˇˇą κ do
	4:	Update A, B and C sequentially; see Equations (4.2)
	5:	Compute the objective function Φ i`1
	6: end while
	7: % Loop 2 : update of A, B, C and P i
	8: i Ð 0 9: while i ă i2 and ˇˇΦ i ´Φi`1 Φ i	ˇˇą κ do
	11:	Compute the objective function Φ i`1
	12: end while 13: return Ŷ "	ř R r"1 `Ar B T r ˘b c

10: Update A, B, C and P i for i P t1, ..., 3u sequentially; see Equations (4.2) and (4.3) r Algorithm 4.2 Initialization of Algorithm 4.1 Require: Input tensors Y 1 ě 0, Y 2 ě 0, downsampling operators P i ě 0 for i P t1, ..., 3u; R, ranks tL r u R r"1 , maximum number of iterations j1 and j2, a threshold 0 ă κ ! 1. Ensure: Initial values A, B, C 1: Initialize C using VCA on Y 1 ; 2:

Table 1

 1 5.5. Unmixing results. For the unmixing task, the reference and estimated spectral signatures and abundance maps are depicted. Due to space limitations, Reconstruction of the synthetic dataset with Poisson noise.

	Method	CC (Ò) SAD (Ó) RMSE (Ó) ERGAS (Ó) PSNR (dB) (Ò) Time (sec) (Ó) C.R. (Ò)
	STEREO	0.999	2.82	8.882e-3	1.347	35.31	1.573	120
	BSTEREO	0.995	7.98	2.134e-2	2.809	29.93	1.379	101
	SCOTT	0.999	2.77	1.122e-2	1.518	34.30	2.087	234
	BSCOTT	0.999	1.47	1.146e-2	1.601	34.11	0.204	159
	SCLL1	0.987	5.70	3.609e-2	5.434	24.46	33.57	266
	CT-STAR	0.999	0.25	5.272e-3	0.713	39.71	0.102	242
	CB-STAR	0.999	1.25	6.408e-3	0.991	38.02	20.14	138
	CNN-BTD-Var	0.997	5.30	1.260e-2	2.952	32.76	3.302	266
	CNMF	0.999	1.09	7.313e-3	1.364	36.20	1.872	43
	FUSE	0.994	2.77	1.470e-2	2.813	30.92	0.146	6
	HySure	0.998	3.04	1.375e-2	2.174	32.35	17.48	43
	SFIM	0.990	4.96	2.544e-2	3.817	28.07	0.281	6
	MR-β(=1)-NMF	0.998	2.22	7.026e-2	1.581	37.08	586.77	43
	Alg. 4.1 (β = 1)	0.999	1.19	1.065e-2	0.776	35.69	23.42	497
	Alg. 4.3 (β = 1)	0.999	2.96	8.703e-3	2.007	36.88	30.19	497
	Method	CC (Ò) SAD (Ó) RMSE (Ó) ERGAS (Ó) PSNR (dB) (Ò) Time (sec) (Ó) C.R. (Ò)
	STEREO	0.989	3.01	1.378e-2	2.617	31.10	1.597	92
	BSTEREO	0.989	3.05	1.435e-2	2.479	30.69	1.673	82
	SCOTT	0.991	3.30	1.584e-2	2.361	30.20	0.180	191
	BSCOTT	0.981	3.55	2.282e-2	3.191	28.11	0.213	127
	SCLL1	0.979	5.23	2.983e-2	3.656	25.46	11.387	199
	CT-STAR	0.886	8.63	5.328e-2	7.204	19.53	0.092	456
	CB-STAR	0.992	2.89	1.371e-2	2.269	31.09	7.592	116
	CNN-BTD-Var	0.920	8.63	4.152e-2	6.078	21.45	1.130	199
	CNMF	0.999	0.90	6.000e-3	1.493	27.96	1.302	42
	FUSE	0.984	2.88	2.066e-2	3.227	28.15	0.276	8
	HySure	0.998	2.05	1.196e-2	1.778	33.55	13.734	42
	SFIM	0.989	3.71	2.083e-2	2.832	29.54	0.320	8
	MR-β(=1)-NMF	0.964	5.63	2.111e-2	6.514	27.81	61.268	42
	Alg. 4.1 (β " 1)	0.990	3.16	1.712e-2	2.657	29.58	42.241	87
	Alg. 4.3 (β = 1)	0.986	3.95	2.144e-2	2.88	27.58	19.77	87

Table 2

 2 Reconstruction metrics on the Jasper Ridge dataset, Poisson noise.

Table 3

 3 Reconstruction metrics on the Ivanpah Playa dataset, Poisson noise.

				1.421e-2	1.119	32.07	1.640	92
	BSTEREO	0.979	2.84	1.897e-2	1.477	29.90	1.499	81
	SCOTT	0.979	1.80	1.629e-2	1.265	30.58	0.189	196
	BSCOTT	0.870	5.45	3.946e-2	3.320	23.82	0.115	142
	SCLL1	0.974	5.12	6.996e-2	3.985	22.05	2.838	86
	CT-STAR	0.948	1.30	2.257e-2	1.595	27.01	0.053	442
	CB-STAR	0.988	1.24	1.215e-2	0.969	32.31	5.955	196
	CNN-BTD-Var	0.973	1.11	1.827e-2	1.295	28.42	2.727	86
	CNMF	0.988	1.12	8.598e-3	0.766	34.96	1.931	42
	FUSE	0.991	1.05	9.621e-3	0.825	34.02	0.233	8
	HySure	0.939	3.71	2.583e-2	2.083	27.15	13.852	42
	SFIM	0.963	2.08	1.605e-2	1.289	30.02	0.258	8
	MR-β(=1)-NMF	0.983	1.82	1.415e-2	1.573	32.26	152.81	42
	Alg. 4.1 (β " 1)	0.983	1.75	1.489e-2	1.176	31.43	46.472	125
	Alg. 4.3 (β = 1)	0.984	1.74	1.502e-2	1.172	31.33	19.16	125

Table 4

 4 Semi-blind reconstruction on the Jasper Ridge dataset.

  T 2 ˆ´P 2 BH 1 " pC d p P 1 Aq T and H 2 " pP 3 C d p Aq T . For C, we have

				"	¯.pβ´2q 2 ´P2 BH 1 d Y 1 p2q P T ¯.pβ´1q ˙HT 1 H T 1	`λ ˆ´B H 2 `λ ´BH 2 ¯.pβ´1q ¯.pβ´2q H T d Y 2 	p2q 2	˙HT 2		' .γpβq ‹ ‹	,
	where H 1 (A.2)
	C " C d	¨"ˆ´C	H 1	´CH 1 ¯.pβ´2q d Y ¯.pβ´1q p3q 1 ˙HT H T 1 `λP T 3 ´P3 CH 2 1 `λP T 3 ˆ´P 3 CH 2 ¯.pβ´1q ¯.pβ´2q H T 2 d Y " 	p3q 2	˙HT 2		' ‹ ‹ .γpβq

This manuscript is for review purposes only.

available for download at https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/ document-library/-/assetpublisher/Wk0TKajiISaR/content/sentinel-2a-spectral-responses.

Available for download at http://lesun.weebly.com/hyperspectral-data-set.html.This manuscript is for review purposes only.

˚Submitted to the editors November 6, 2023. This work is an extension of [33]. This work was partly supported by the ANR project "Chaire IA Sherlock" ANR-20-CHIA-0031-01 hold by P. Chainais, as well as by the national support within the programme d'investissements d'avenir ANR-16-IDEX-0004 ULNE and Région HDF. VL acknowledges the support by Ministry of Science and Higher Education grant No. 075-10-2021-068.