Student Member, IEEE, Tales Imbiriba, Member, IEEE Ricardo Augusto Borsoi

Senior Member, IEEE José Carlos Moreira Bermudez

J C M Bermudez
email: j.bermudez@ieee.org.

C Richard
email: cedric.richard@unice.fr.

A Blind Multiscale Spatial Regularization Framework for Kernel-based Spectral Unmixing

Abstract-Introducing spatial prior information in hyperspectral imaging (HSI) analysis has led to an overall improvement of the performance of many HSI methods applied for denoising, classification, and unmixing. Extending such methodologies to nonlinear settings is not always straightforward, specially for unmixing problems where the consideration of spatial relationships between neighboring pixels might comprise intricate interactions between their fractional abundances and nonlinear contributions. In this paper, we consider a multiscale regularization strategy for nonlinear spectral unmixing with kernels. The proposed methodology splits the unmixing problem into two sub-problems at two different spatial scales: a coarse scale containing low-dimensional structures, and the original fine scale. The coarse spatial domain is defined using superpixels that result from a multiscale transformation. Spectral unmixing is then formulated as the solution of quadratically constrained optimization problems, which are solved efficiently by exploring their strong duality and a reformulation of their dual cost functions in the form of root-finding problems. Furthermore, we employ a theory-based statistical framework to devise a consistent strategy to estimate all required parameters, including both the regularization parameters of the algorithm and the number of superpixels of the transformation, resulting in a truly blind (from the parameters setting perspective) unmixing method. Experimental results attest the superior performance of the proposed method when comparing with other, state-of-theart, related strategies.

Index Terms-Hyperspectral data, multiscale, spatial regularization, nonlinear unmixing, kernel methods.

I. INTRODUCTION

Modern remote sensing greatly relies on hyperspectral (HS) image analysis to retrieve information about surface materials in many applications such as agriculture, surveillance and space exploration [START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF]. Specifically, reflectance measures can provide detailed information about the spectral signature of pure materials present on the surface of a target scene and their proportion for each pixel. Applications often, but not exclusively, associated with remote sensing trade poor spatial resolution for high spectral resolution due to physical limitations of imaging sensors and to the distance between the sensor and the target scene. Therefore, the measured reflectance of a given pixel is usually a mixture of the pure spectral signatures of materials existing in the corresponding area. Spectral unmixing (SU) consists of extracting the pure component spectral signatures and their proportions (or abundances) for each pixel. The literature presents many mixing models to explain the observed reflectance as a mathematical function of the pure spectral components. The simplest form of such models is the linear mixing model (LMM) which confines the observed reflectance vectors into a convex hull whose extremities are the pure component spectral signatures, therefore, called endmembers. The LMM is effective in accurately modelling mixtures occurring in scenes where the materials of interest cover a large area with respect to the pixel size [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF]. It however disregards more complex mixing phenomena such as non-linearity [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF], [START_REF] Imbiriba | Nonparametric detection of nonlinearly mixed pixels and endmember estimation in hyperspectral images[END_REF] and spectral variability [START_REF] Somers | Endmember variability in spectral mixture analysis: A review[END_REF]- [START_REF] Borsoi | Deep generative models for library augmentation in multiple endmember spectral mixture analysis[END_REF], which often results in estimation errors being propagated throughout the unmixing process [START_REF] Thouvenin | Hyperspectral unmixing with spectral variability using a perturbed linear mixing model[END_REF].

Nonlinear interactions between materials occur in many scenes where there is complex radiation scattering among several endmembers, such as in some vegetation areas [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF]. In such situations, nonlinear mixing models must be considered [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF], [START_REF] Heylen | A review of nonlinear hyperspectral unmixing methods[END_REF]. Several nonlinear SU strategies have been proposed in the literature, which can be roughly divided between model-based and model-free methods. Most modelbased nonlinear SU algorithms assume that the mixing process that occurs in the scene is known a priori [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF], [START_REF] Heylen | A review of nonlinear hyperspectral unmixing methods[END_REF]- [START_REF] Altmann | Supervised nonlinear spectral unmixing using a postnonlinear mixing model for hyperspectral imagery[END_REF]. However, real mixing mechanisms can be very complex and prior knowledge about them is seldom available in practice. This led to the consideration of more flexible model-free nonlinear SU, which employ more flexible nonlinear mixing models that are able to represent generic functions. Prominent model-free strategies include the estimation of abundances as posterior class probabilities of a nonlinear classifier [START_REF] Li | Abundance estimation from hyperspectral image based on probabilistic outputs of multi-class support vector machines[END_REF], the use of graph-based approximate geodesic distances [START_REF] Heylen | Non-linear spectral unmixing by geodesic simplex volume maximization[END_REF], [START_REF] Heylen | A distance geometric framework for nonlinear hyperspectral unmixing[END_REF],

and kernel-based algorithms [START_REF] Wu | A kernel spatial complexity-based nonlinear unmixing method of hyperspectral imagery[END_REF]- [START_REF] Altmann | Nonlinear spectral unmixing of hyperspectral images using gaussian processes[END_REF]. Kernel-based methods provide non-parametric representations of functional spaces that are able to model arbitrary nonlinear mixtures [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF], [START_REF] Heylen | A review of nonlinear hyperspectral unmixing methods[END_REF], [START_REF] Wu | A kernel spatial complexity-based nonlinear unmixing method of hyperspectral imagery[END_REF]- [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF], [START_REF] Ammanouil | Nonlinear unmixing of hyperspectral data with vector-valued kernel functions[END_REF]. This flexibility, allied to a good experimental performance has led to the wide application of kernel methods.

Despite the good results obtained with kernel-based unmixing methods [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF], most algorithms fail to explore the high spatial regularity associated to many real world scenes. This property can be leveraged to improve the conditioning of the unmixing problem. Spatial regularization has already been shown to improve the performance of linear [START_REF] Shi | Incorporating spatial information in spectral unmixing: A review[END_REF], [START_REF] Imbiriba | A low-rank tensor regularization strategy for hyperspectral unmixing[END_REF] and sparse [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF], [START_REF] Feng | Adaptive spatial regularization sparse unmixing strategy based on joint MAP for hyperspectral remote sensing imagery[END_REF] SU, as well as spectral-variability-aware SU [START_REF] Imbiriba | Generalized linear mixing model accounting for endmember variability[END_REF], [START_REF] Drumetz | Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability[END_REF]- [START_REF] Hong | Sulora: Subspace unmixing with low-rank attribute embedding for hyperspectral data analysis[END_REF], which is closely connected to nonlinear SU [START_REF] Drumetz | Relationships between nonlinear and space-variant linear models in hyperspectral image unmixing[END_REF]. However, spatial information has seldom been enforced in nonlinear unmixing algorithms, partly due to the challenges associated with more complex observation models. For instance, a spatial clustering approach was used in [START_REF] Tang | Integrating spatial information in the normalized p-linear algorithm for nonlinear hyperspectral unmixing[END_REF] to divide the image into different groups of pixels. SU was then performed using the P-linear mixing model in a Bayesian framework with a unique set of regularization parameters for each group. In [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF] a Total Variation (TV) regularization was introduced in a regression-based kernel unmixing method [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF], and a variable splitting approach was then used to solve the resulting optimization problem.

The TV regularization has been widely used in many HS imaging tasks since it promotes smooth image reconstructions while still allowing for sharp discontinuities [START_REF] Iordache | Total variation spatial regularization for sparse hyperspectral unmixing[END_REF]. However, TV regularization is not the most effective approach to extract spatial information from hyperspectral images. Regularization strategies exploiting nonlocal redundancy in images were recently considered for SU [START_REF] Wang | Centralized collaborative sparse unmixing for hyperspectral images[END_REF], [START_REF] Yao | Nonconvex-sparsity and nonlocal-smoothness-based blind hyperspectral unmixing[END_REF], leading to a better abundance estimation performance at the expense of an increase in computational complexity. Other works exploited the manifold structure in the hyperspectral data by using a graph-based regularization to connect the abundances of pixels that are similar with respect to some metrics [START_REF] Lu | Manifold regularized sparse NMF for hyperspectral unmixing[END_REF], [START_REF] Ammanouil | A graph laplacian regularization for hyperspectral data unmixing[END_REF].

In [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF], a multiscale spatial regularization approach was proposed for sparse spectral unmixing. The multiscale approach led to improved results and smaller computational complexity when compared to TV regularization. The unmixing problem was split into two simpler problems in different image domains defined using a multiscale transformation. This transformation groups image pixels into contiguous regions using (over)-segmentation strategies such as the superpixel decomposition [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF]. The multiscale regularization strategy was later extended in [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF] to consider SU accounting for spectral variability. Despite the excellent results obtained with spatial regularization strategies, their performance usually depends on the careful selection of regularization parameters. This is specially important in multiscale strategies, which require a larger number of parameters. Determination of parameters for spatial regularization methods remains a challenging problem, and works applied to HSI are rare [START_REF] Feng | Adaptive spatial regularization sparse unmixing strategy based on joint MAP for hyperspectral remote sensing imagery[END_REF], [START_REF] Song | Regularization parameter estimation for non-negative hyperspectral image deconvolution[END_REF].

In this paper, we propose a new multiscale spatial regularization approach for kernel-based nonlinear unmixing. Building upon the ideas proposed in [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF], we employ a multiscale representation to divide the unmixing problem into two simpler problems in different scales. Though based on the same principle used in [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF], devising kernel-based mixing models in multiple scales is more challenging than in the linear case. Moreover, we address the parameter adjustment problem differently from what has been done in previous multiscale SU formulations in [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF], [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF]. In this work we reformulate the SU problem at multiple scales by statistically characterizing not only the algorithm reconstruction error in both scales, but also the inter-scale interaction between the abundances and the nonlinear mixing contributions across the coarse and fine image domains. This formulation leads to physically motivated constraints which are leveraged to devise the Blind Multiscale Unmixing Algorithm for Nonlinear Mixtures (BMUA-N), in which all the parameters are determined automatically from the observed data. Thus, the proposed strategy benefits from an improved quality without the need for ad hoc parameter adjustment such as in TV-based works.

We formulate the resulting unmixing problem as a sequence of two optimization problems with quadratic equality constraints. These non-convex problems are solved by reformulating the dual problems in the form of low-dimensional root finding problems, which can be solved in very few iterations using a multidimensional bisection algorithm. Moreover, we are able to prove that, under mild conditions, strong duality holds for these optimization problems, which guarantees the optimality of this approach. Simulations with synthetic and real datasets illustrate the effectiveness of the proposed methodology in producing piecewise smooth solutions while preserving sharp discontinuities existing in the image. This leads to more accurate unmixing results when compared to TV-based strategies, with less computational complexity and without the need for ad hoc parameter adjustment.

This manuscript is organized as follows. In Section II we discuss the main concepts related to regression-based kernel unmixing. In Section III we present the proposed kernelbased multiscale unmixing strategy. The automatic parameter setting methodology is presented in Section IV and the solution for the proposed optimization problems is discussed in Section V. In Section VI we propose a method for designing the multiscale transformation to yield spectral homogeneity. Section VII discusses the computational complexity of the proposed method. Experimental results are presented and discussed in Section VIII. They are followed by concluding remarks in Section IX.

II. KERNEL-BASED UNMIXING

A. Kernel-based mixture model

In this section we review the standard kernel-based mixture model introduced in [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF] and discuss the main theoretical aspects of kernel machines. As in [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF], we assume that each L-band observed pixel y n ∈ R L in an HSI can be modeled as a function of the endmember spectra as follows:

y n, = ψ an ( m ) + e n, , = 1, . . . , L , (1)
where y n, is the -th entry of vector y n , m ∈ R 1×P is the -th row of the endmember matrix M ∈ R L×P with P spectral signatures of pure materials in the scene. Function ψ an is an unknown nonlinear function defining the interactions between endmember spectra parameterized by their fractional abundances a n ∈ R P . e n, includes the observation noise and modeling errors. The problem that arises is to find a functional ψ an that can accurately represent the different and complex types of light-endmember interactions often occurring in real scenes. Since the type of nonlinearity is rarely known in practice, a popular solution is to search for kernel-based smooth function representations whose parameters can be learned directly from the data [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF], [START_REF] Imbiriba | Nonparametric detection of nonlinearly mixed pixels and endmember estimation in hyperspectral images[END_REF].

In [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF], [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF] the authors considered a semi-parametric kernel-based model consisting of a linear trend parameterized by the abundance vector plus an additive nonlinear fluctuation. The model, which allows the quantification of the abundance vectors during the unmixing process, is given by

ψ an ( m ) = m a n + ψ n ( m ) , (2)
with ψ n : R P → R being an arbitrary smooth function belonging to a Reproducing Kernel Hilbert Space (RKHS) denoted by H and defined over a nonempty compact set M ⊂ R P . This assumption allows for the kernel machinery (i.e., via the kernel trick) to obtain accurate solutions to the unmixing problem. The theory of positive definite kernels emerged from the study of positive definite integral operators [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF], and was further generalized in the study of positive definite matrices [START_REF] Moore | On properly positive hermitian matrices[END_REF]. It has been established that, to every positive definite function κ(•, •) : M × M → R, defined over a non-empty compact set M ⊂ R P , there corresponds one and only one class of realvalued functions on M forming a Hilbert space H endowed with a uniquely defined inner product •, • H , and admitting κ as a reproducing kernel (r.k.) [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. Space H is called a RKHS if its evaluation functional δ m is a linear and continuous (or equivalently bounded) functional for every m ∈ M, thus, admitting κ as its unique kernel. As a consequence of the Riesz representation theorem [46, p. 188], κ(•, m) is the representer of evaluation of any functional ψ ∈ H, such that the reproducing property

ψ( m) = ψ, κ(•, m) H (3)
holds, for all ψ ∈ H and all m ∈ M. Furthermore, since κ(•, m) ∈ H for all m, m ∈ M we also have

κ( m, m ) = κ(•, m), κ(•, m ) H . (4)
The RKHS H is then formed by a class of functions generated by all functions of the form ψ(•) = j α j κ(•, m j ), with norm defined by ψ 2 H = i j α i α j κ( m i , m j ). In the context of machine learning, kernel methods are often related with the concept of building a high dimensional feature space H, and a mapping

Φ : M -→ H m -→ Φ( m) , (5)
with inner product defined as κ( m, m ) = Φ( m), Φ( m ) H . If κ is a r.k. of H, then H is a RKHS and also a feature space of κ with Φ( m) = κ(•, m). In this case Φ is called the canonical feature map [47, p. 120]. This leads to the so-called "kernel trick" allowing one to compute inner products of data mapped into higher, or even infinite, dimensional feature spaces by evaluating a real function κ( m i , m j ) in the input space.

Although the literature proposed a variety of kernel functions elaborated during the past two decades of intense research activity [START_REF] Vapnik | The nature of statistical learning theory[END_REF]- [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], in this work we restrain ourselves to the polynomial kernel due to its intimate relation with multiple scattering phenomena known to exist in the interaction between light and the materials in the scene. Thus, the polynomial kernel is given by

κ( m i , m j ) = ( m i m j + c) d , ( 6
)
where d is the polynomial degree and c ≥ 0 is a real number. Due to relevant findings reported in [START_REF] Somers | Nonlinear hyperspectral mixture analysis for tree cover estimates in orchards[END_REF] concerning the order of multiple reflection models and the good results obtained in [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF], in this paper we assume d = 2 and c = 1 in all simulations.

B. LS-SVR-based unmixing

In [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF] the authors proposed to solve the unmixing problem accounting for the model in ( 1)-( 2) by considering a multikernel generalization of standard least-squares support vector regression (LS-SVR) methods [START_REF] Suykens | Least Squares Support Vector Machines[END_REF]. The resulting optimization problem is given by

(â n , ψn ) = arg min an, ψn 1 2 a n 2 + ψ n 2 H + 1 µ ξ n 2 2 (7)
subject to

a n ≥ 0 , 1 a n = 1 , ξ n, = y n, -a n m -ψ n ( m ), = 1, . . . , L ,
where ân and ψn are the estimated abundance vector and nonlinear function for the n-th pixel. Problem ( 7) is solved using standard dual formulation based on the Lagrangian [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF]. Although problem [START_REF]Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF] presents an effective way of modeling both the linear trend and the nonlinear mixing occurring in a given pixel, it fails to impose any smooth structure over the abundance estimation within neighboring pixels.

Standard regularization approaches such as the TV have considered an additional term to the cost function in problem [START_REF]Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF] that penalizes spatial discontinuities in the abundance maps [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF]. However, as discussed in the introduction, this strategy is not the most effective in exploring spatial information contained in the image. Besides, it introduces an additional parameter that must be carefully tuned in order to achieve good performance.

In the following, we present a multiscale formulation for the nonlinear mixing model in ( 1)-( 2) that enables us to better exploit the spatial regularity in HSIs, leading to improved results when compared to the TV-based strategy. Compared to linear unmixing, the application of a multiscale formulation to kernel-based nonlinear unmixing leads to specific challenges that need to be addressed. This is specially true in the present work, as quadratic equality constraints need to be reformulated to allow for the automatic determination of the parameters, and thus for a blind algorithm as detailed in Section IV.

III. A MULTISCALE NONLINEAR MIXING MODEL

Traditional regularization approaches (e.g., Tikhonov or TV) introduce spatial regularity by promoting similarity between abundances at fixed spatial neighborhoods. Recently, more flexible approaches emerged exploring irregular and datadependent image structures that can generally be described under a graph or manifold regularization framework [START_REF] Lu | Manifold regularized sparse NMF for hyperspectral unmixing[END_REF], [START_REF] Ammanouil | A graph laplacian regularization for hyperspectral data unmixing[END_REF], [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. In this case, a graph is first constructed to represent the similarity between pairs of pixels in the HSI using distance metrics that can be either defined explicitly [START_REF] Stevens | Spectral-densitybased graph construction techniques for hyperspectral image analysis[END_REF] or derived indirectly using image (over)-segmentation methods such as superpixels, ultrametric contour maps (UCM), or binary partition trees (BPT) [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]- [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral unmixing-based binary partition tree representation[END_REF]. Afterwards, abundances corresponding to pixels that are similar according to the graph can be constrained to have similar proportions in the SU problem, which preserves the geometric structure found in the HSI [START_REF] Lu | Manifold regularized sparse NMF for hyperspectral unmixing[END_REF], [START_REF] Ammanouil | A graph laplacian regularization for hyperspectral data unmixing[END_REF], [START_REF] Wang | Spatial group sparsity regularized nonnegative matrix factorization for hyperspectral unmixing[END_REF]. Despite providing a significant amount of flexibility, graph regularization can still lead to computationally costly SU problems. More importantly though, this approach does not provide a clear insight about the multiscale abundance interactions in a way that could motivate the design of the algorithms, which makes the selection of the regularization parameters difficult in practical scenarios.

Motivated by the results in [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF], [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF], we propose to introduce spatial information into SU by representing this problem separately in two spatial scales, which significantly reduces the computational complexity of SU. Moreover, this also makes the inter-scale interaction between the abundances explicit in the resulting optimization problem, which is essential for the theoretically principled parameter design strategy presented in Section IV. Specifically, we divide the SU problem into two consecutive steps. First, we represent the nonlinear mixing process in an approximation (coarse) spatial scale (C) which preserves relevant inter-pixel spatial contextual information. Pixels in the coarse spatial scale can be unmixed independently from each other. The recovered coarse abundance maps are then mapped back to the original image domain (D) and used as prior information to regularize the second unmixing process applied to the original image to promote spatial dependency between neighboring pixels.

A. Unmixing in the coarse scale

Denote the HSI and the abundance map for all pixels by Y = [y 1 , . . . , y N ] and A = [a 1 , . . . , a N ], respectively. We consider a dimensionality reduction transformation W ∈ R N ×K , K < N , constructed based on relevant contextual inter-pixel information present in the observed image Y , that maps both the HSI and the abundance map to the approximation domain. The transformed matrices are given by

Y C = Y W , A C = AW , (8)
where

Y C = [y C1 , . . . , y C K ] ∈ R L×K and A C = [a C1 , . . . , a C K ] ∈ R P ×K
are, respectively, the HSI and the abundance matrix in the coarse approximation scale.

Various methods can be used to construct the transformation W . Specifically, W must group pixels that are spatially adjacent and spectrally similar and must respect image borders by not grouping pixels corresponding to different image structures or features. Following the same approach as in [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF], [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF], we consider the superpixel decomposition of the image Y for the transformation W . Besides satisfying the criteria outlined above, multiscale decompositions based on superpixel algorithms have shown excellent performance in SU considering both sparsity [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF] and variability of the endmembers [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF]. Superpixel algorithms group image pixels into different spatially compact neighborhoods with similar spectral information [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF], decomposing the image into a set of contiguous homogeneous regions whose size and regularity are controlled by adjusting a set of parameters. Furthermore, the superpixel decomposition can be computed very efficiently by employing low-cost algorithms such as the SLIC [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. The transformation W is thus constructed such that Y W computes the superpixel decomposition of the image Y , and returns the averages of all pixels inside each superpixel region. The effect of the transformation W on the Cuprite HSI is illustrated in Figure 1. Note that, besides the superpixel decomposition, other image (over)-segmentation strategies can also be used to construct W , such as e.g., UCM or BPT [START_REF] Arbelaez | Boundary extraction in natural images using ultrametric contour maps[END_REF], [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral unmixing-based binary partition tree representation[END_REF], which can also provide hierarchical (multiscale) representations of HSIs.

Considering the nonlinear observation model ( 1), the transformed image in (8) leads to an equivalent mixing model in the coarse spatial domain, which is given by

y C i, = 1 |N i | n∈Ni m a n + ψ n ( m ) + e n, = m a Ci + ψ Ci ( m ) + 1 |N i | n∈Ni e n, , (9)
where y C i, is the -th entry of y Ci , N i is the set of indexes of the pixels contained within the i-th superpixel, | • | denotes the cardinality of a set, and a Ci and ψ Ci , given by

a Ci = 1 |N i | n∈Ni a n , ψ Ci ( m ) = 1 |N i | n∈Ni ψ n ( m ) ,
for i = 1, . . . , K, denote the fractional abundances and nonlinear contributions at the coarse spatial scale.

Following the observation model ( 9), the SU problem at the coarse spatial scale can be formulated using the LS-SVR framework presented in Section II-B, which leads to the following optimization problem:

{â Ci , ψCi } = arg min {a C i ,ψ C i ,ξ C i } 1 2 K i=1 ψ Ci 2 H + a Ci 2 2 (10) subject to a Ci ≥ 0 , 1 a Ci = 1 , i = 1, . . . , K , ξ Ci = y Ci -M a Ci -ψ Ci (M ), i = 1, . . . , K , 1 K K i=1 ξ Ci 2 2 = C 0 ,
where âCi and ψCi are the estimated abundance vector and nonlinear function for the i-th superpixel and ψ Ci (M ) = [ψ Ci ( m 1 ), . . . , ψ Ci ( m L )] . Parameter C 0 is a positive constant that constrains the reconstruction error of the algorithm, and operates in an analogous manner to a regularization. Differently from [START_REF]Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF], in [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] we choose to limit the reconstruction error using an equality constraint instead of directly adding ξ Ci 2 2 , i = 1, . . . , K as a regularizer in the cost function. The rationale is that, unlike the regularization parameter µ in [START_REF]Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF], the constant C 0 has a clear physical interpretation. This characteristic is exploited in the next section where we present a proper methodology for automatically setting C 0 .

B. Unmixing in the image domain

The abundance maps estimated at the coarse spatial scale, denoted by A C = [a C1 , . . . , a C K ], can be used to regularize the original unmixing problem. To this end, we convert the abundance map from the coarse approximation domain C back to the original image domain D as

A D = A C W * . ( 11
)
Matrix W * ∈ R K×N is a conjugate transformation to W , and takes the image from the coarse domain C back to the original (uniform) image domain. This is performed by attributing the value a Ci to all pixels in A D that lie within the i-th superpixel. Thus, A D can be viewed as a coarse version of A in the original image domain. After computing A D using [START_REF] Heylen | A review of nonlinear hyperspectral unmixing methods[END_REF], the SU problem for all pixels is given by

{â n , ψn } = arg min {an,ψn,ξ n } 1 2 N n=1 ψ n 2 H (12) subject to a n ≥ 0, 1 a n = 1 , n = 1, . . . , N , ξ n = y n -M a n -ψ n (M ) , n = 1, . . . , N , 1 N N n=1 ξ n 2 2 = C 1 , 1 N N n=1 a n -a Dn 2 2 = C 2 ,
where ân and ψn are the estimated abundance vector and nonlinear function for the n-th pixel in the original image domain,

ψ n (M ) = [ψ n ( m 1 ), . . . , ψ n ( m L )] ,
and C 1 , C 2 are positive constants that constrain the reconstruction error and the abundance variability across scales. Again, as in [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF], we use equality constraints instead of additive penalty terms in the cost function due to the easier interpretation of C 1 and C 2 when compared to regularization parameters. This improved interpretability is exploited in the next section to provide a methodology for automatically adjusting these constants. We note that the nonlinear equality constraints make the optimization problems [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] and ( 12) non-convex. An efficient algorithm will be proposed in Section V to address this issue.

IV. DETERMINING THE REGULARIZATION CONSTANTS

A significant challenge in regularized unmixing algorithms consists in determining the regularization constants. The proposed formulation requires the selection of C 0 , C 1 and C 2 in problems [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] and [START_REF] Mustard | Photometric phase functions of common geologic minerals and applications to quantitative analysis of mineral mixture reflectance spectra[END_REF]. Although most works assume that such constants can be determined empirically, several frameworks have been proposed to select regularization parameters in (ill-posed) inverse problems, such as the generalized cross validation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF], risk estimation for reconstruction error minimization [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF], [START_REF] Ammanouil | ADA-PT: An adaptive parameter tuning strategy based on the weighted stein unbiased risk estimator[END_REF] or the L-curve [START_REF] Belge | Efficient determination of multiple regularization parameters in a generalized L-curve framework[END_REF]. However, determining these parameters blindly for a given problem is usually difficult and computationally intensive. This motivates the consideration of external a priori information whenever possible, which allows the application of conceptually simpler and more reliable solutions.

Classical approaches in this latter category are the Chisquared method [START_REF] Hall | Common structure of techniques for choosing smoothing parameters in regression problems[END_REF]- [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF] and the closely related discrepancy principle [START_REF] Heinz Werner Engl | Regularization of Inverse Problems[END_REF]. The underlying idea behind these methods is to evaluate the statistical properties of the data reconstruction term in the cost function when the estimated solution is equal to the desired (i.e., true) parameters. For an observation model such as [START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF], evaluating the true parameters in the cost function of [START_REF]Deep generative endmember modeling: An application to unsupervised spectral unmixing[END_REF] would result in the residuals ξ n being equal to the observation noise e n -which, under a Gaussianity hypothesis, makes the reconstruction error Chi-squared-distributed, giving the technique its name. The regularization constant is then selected so that the estimated solution yields a residual with the same statistical properties of the observation noise, which are assumed to be known [START_REF] Hall | Common structure of techniques for choosing smoothing parameters in regression problems[END_REF]- [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF]. This can be performed using equality constrained optimization problems such as in [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] and [START_REF] Mustard | Photometric phase functions of common geologic minerals and applications to quantitative analysis of mineral mixture reflectance spectra[END_REF] [START_REF] Hunt | The application of constrained least squares estimation to image restoration by digital computer[END_REF], which immediately translates the choice of constants C 1 , C 2 and C 3 as the problem of determining the statistical properties of the corresponding equality terms.

In this section, we will extend this strategy to the present SU problem by determining the statistical averages of the reconstruction error in both scales and the inter-scale abundance variation. Although in principle this requires knowledge about the true abundance solutions we want to estimate (which historically limited the applicability of the Chi-square technique to more complex problems like ( 12)), we will be able to provide a theoretically sound and yet simple strategy for its reliable estimation.

More precisely, constant C 0 in (10) reflects the average noise power in the coarse image scale. Constants C 1 and C 2 in (12) reflect, respectively, the average noise power in the detail scale and the average energy of the differences between the fractional abundances and their estimates in the coarse domain. Under the Chi-squared framework, these constants are defined in terms of statistical means as

C 0 = E 1 K K i=1 y Ci -M a Ci -ψ Ci (M ) 2 2 = E 1 K K i=1 e Ci 2 2 , (13)
C 1 = E 1 N N n=1 y n -M a n -ψ n (M ) 2 2 = E 1 N N n=1 e n 2 2 , (14)
C 2 = E 1 N N n=1 a n -a Dn 2 2 , (15)
where E{•} is the expected value operator with respect to the true distribution of the parameters within the brackets, and e n and e Ci are the errors related to the fine and coarse domains of the n-th pixel and i-th superpixel, respectively. In this case, the quadratic constraints in problems ( 10) and ( 12) can be seen as approximations to the above expected values. Constants C 0 and C 1 depend directly on the noise level and modeling errors represented by e n . We write

e n = e 0,n + e ψ,n , (16)
where e 0,n is a Gaussian noise term with zero mean and covariance Σ e , and e ψ,n represents modeling errors. Before proceeding, we make the following assumptions:

A1) The additive noise is spatially uncorrelated, i.e., E{e 0,n e 0,m } = 0, ∀ n = m. A2) The noise and modeling errors e 0,n and e ψ,n are uncorrelated. A3) The modeling errors e ψ,n are assumed to be spatially correlated and approximately constant within each superpixel, that is:

e ψ,n ≈ e ψ,m , ∀ m, n ∈ N i , i = 1, . . . , K . ( 17
)
This hypothesis is motivated by the spatial smoothness of both the abundances and the nonlinear contributions in the mixing model [START_REF] Altmann | Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis[END_REF], [START_REF] Ammanouil | Spatial regularization for nonlinear unmixing of hyperspectral data with vectorvalued kernel functions[END_REF]. A4) The expected value of the modeling error's norm is the same for all pixels, and is represented as

E e ψ,n 2
A5) The noise covariance matrix is the same for all image pixels, i.e., E e 0,i e 0,i = E e 0,j e 0,j ,

1 ≤ i, j ≤ N = Σ e . (19)
A6) The vectors a n -a Dn , M † (ψ n (M ) -ψ Cn (M ) and

M † e n -1
|Nn| i∈Nn e i ) are mutually uncorrelated and zero-mean. We denote the average size of each superpixel by S = N/K.

In the following, we will evaluate the expectations in ( 13)-( 15) in order to provide well-founded means to select the constants C 0 , C 1 and C 2 .

A. Determining the constant C 1

Using hypothesis A1, A2, A4 and A5, constant C 1 can be computed as:

C 1 = E 1 N N n=1 e 0,n + e ψ,n 2 2
(A2) = 1 N N n=1 E e 0,n 2 2 + E e ψ,n 2 2
(A4,A5)

= tr{Σ e } + σ 2 e,ψ . (20)
B. Determining the constant C 0 Constant C 0 can be derived in a similar way by assuming A1-A4:

C 0 = E 1 K K i=1 e Ci 2 2 = E 1 K K i=1 1 |N i | n∈Ni e 0,n + e ψ,n 2 2 = E 1 K K i=1 1 |N i | 2 n∈Ni e 0,n 2 2 + n∈Ni e ψ,n 2 2 + 2 n∈Ni m∈Ni e 0,n , e ψ,m (A2) = 1 K K i=1 1 |N i | 2 n∈Ni E e 0,n 2 2 + 1 K K i=1 1 |N i | 2 E n∈Ni e ψ,n 2 2 . (21)
The modeling errors are not uncorrelated and zero mean in each superpixel. By approximating K -1 K i=1 |N i | -1 S -1 and using hypothesis A3, C 0 can be approximated as

C 0 1 S tr{Σ e } + 1 K K i=1 1 |N i | 2 n∈Ni E |N i | e ψ,n 2 2 1 S tr{Σ e } + σ 2 e,ψ . (22)
Note that this shows that the modeling errors are more significant relative to the noise at the coarse spatial scale. This is because the contribution of the noise is reduced by a factor of S whereas the modeling errors retain the same energy.

C. Determining the constant C 2

The constant C 2 is slightly more challenging to compute than the previous ones. We denote the left pseudo-inverse of M as M † ∈ R P ×L , that is, M † M = I P . Then, we can consider the following quantity:

E 1 N N n=1 M † (y n -y Dn ) 2 = E 1 N N n=1 M † M a n + e n + ψ n (M ) -M a Dn -ψ C * n (M ) - 1 |N n | i∈Nn e i 2 . (23)
The expectation on the right hand side of ( 23) can be computed by using assumption A6. This assumption is reasonable since each of the terms in A6 comprises fluctuations between the coarse and fine spatial scales, which can be expected to be of zero mean. Thus, by noting that

E M † (M a n -M a Dn ) 2 = E a n -a Dn 2 (24)
and using A6, we can write [START_REF] Ammanouil | Nonlinear unmixing of hyperspectral data with vector-valued kernel functions[END_REF] as

E 1 N N n=1 M † (y n -y Dn ) 2 = E 1 N N n=1 a n -a Dn 2 + E 1 N N n=1 M † (ψ n (M ) -ψ Cn (M ) 2 + E 1 N N n=1 M † e n - 1 |N n | i∈Nn e i ) 2 . (25)
In the following, we will expand the last term on the right hand side of [START_REF] Imbiriba | A low-rank tensor regularization strategy for hyperspectral unmixing[END_REF]. Using A1-A4, the summand for for each pixel can be written as

E M † e n - 1 |N n | i∈Nn e i 2 = E M † e 0,n + e ψ,n - 1 |N n | i∈Nn e 0,i + e ψ,i 2
(A2) = E M † e 0,n - 1 |N n | i∈Nn e 0,i 2 + E M † e ψ,n - 1 |N n | i∈Nn e ψ,i 2 . ( 26
)
Using hypothesis A3 and equation ( 17), the second term of ( 26) can be approximated as

E M † e ψ,n - 1 |N n | i∈Nn e ψ,i 2 ≈ 0 . (27)
This can be intuitively reasoned by accounting for the spatial correlation of the modeling errors, where e ψ,n in each pixel is very similar to the average of the modeling errors in its respective superpixel. This leads to

E M † e n - 1 |N n | i∈Nn e i ) 2 E M † e 0,n - 1 |N n | i∈Nn e 0,i 2 = E M † e 0,n 1 - 1 |N n | -M † 1 |N n | i∈Nn\{n} e 0,i 2
(A1) = E M † e 0,n 1 - 1 |N n | 2 + E M † 1 |N n | i∈Nn\{n} e 0,i 2
(A1) = E M † e 0,n 1 - 1 |N n | 2 + 1 |N n | 2 i∈Nn\{n} E M † e 0,i 2 (A5) = E M † e 0,n 2 1 - 1 |N n | 2 + |N n | -1 |N n | 2 E M † e 0,n 2 = E M † e 0,n 2 |N n | -1 |N n | , (28)
where "\" denotes the set difference operator. Since the noise statistics are spatially invariant, see A5,

E M † e 0,n 2 = E tr{M † e n e n (M † ) } = M † Σ 1/2 e 2 F (29)
and, by approximating the superpixel sizes by their average value (i.e. |N n | S), we can approximate ( 26) as

E 1 N N n=1 M † e n - 1 |N n | i∈Nn e i ) 2 M † Σ 1/2 e 2 F S -1 S . (30)
Finally, approximating the expectations with respect to a n and ψ n by their instantaneous values and using the estimates âDn and ψCn obtained as solutions to the optimization problem [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF], equation ( 25) can be approximated as

1 N N n=1 a n -a Dn 2 1 N N n=1 M † (y n -y Dn ) 2 - 1 N N n=1 M † (ψ n (M ) -ψ Cn (M ) 2 -M † Σ 1/2 e 2 F S -1 S . (31)

D. The updated fine scale optimization problem

Using these results, we can substitute the second quadratic equality constraint of problem ( 12) by [START_REF] Hong | An augmented linear mixing model to address spectral variability for hyperspectral unmixing[END_REF], resulting in the following problem

arg min {an,ψn,ξ n } 1 2 N n=1 ψ n 2 H (32) subject to a n ≥ 0 , 1 a n = 1 , n = 1, . . . , N , ξ n = y n -M a n -ψ n (M ) , n = 1, . . . , N , 1 N N n=1 ξ n 2 2 = C 1 , 1 N N n=1 a n -a Dn 2 2 = C Y -C E - 1 N N n=1 M † (ψ n (M ) -ψ Cn (M ) 2 ,
where constants C Y and C E are defined as

C Y = 1 N N n=1 M † (y n -y Dn ) 2 , C E = M † Σ 1/2 e 2 F S -1 S . (33)
Defining

ξ ψ,n = M † ψ n (M ) -ψ Cn (M ) , n = 1, . . . , N (34)
and multiplying the quadratic constraints by N/2, we can represent problem (32) equivalently as

arg min {an,ψn,ξ n ,ξ ψ,n } 1 2 N n=1 ψ n 2 H (35) subject to a n ≥ 0 , 1 a n = 1 , n = 1, . . . , N , ξ n = y n -M a n -ψ n (M ) , n = 1, . . . , N , 1 2 N n=1 ξ n 2 2 = N 2 C 1 , 1 2 N n=1 a n -a Dn 2 2 + ξ ψ,n 2 = N 2 (C Y -C E ) , ξ ψ,n = M † ψ n (M ) -ψ Cn (M ) , n = 1, . . . , N .
Problem [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF] can now be used instead of [START_REF] Mustard | Photometric phase functions of common geologic minerals and applications to quantitative analysis of mineral mixture reflectance spectra[END_REF] to perform unmixing in the original spatial scale.

V. SOLVING THE OPTIMIZATION PROBLEMS

The quadratic equality constraints in [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] and ( 35) make the optimization problems non-convex. Furthermore, the functional form of the variables ψ Ci and ψ n makes them hard to be optimized in their primal form. Thus, we resort to a Lagrangian relaxation and solve the dual optimization problem, which is concave and finite-dimensional [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF], [START_REF] Boyd | Convex optimization[END_REF]. Although the non-convexity of the constraints implies the possibility of a non-zero duality gap, we will show in Section V-C that strong duality holds under mild conditions for problems [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] and [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF]. Thus, this approach incurs no loss of performance.

A. The coarse scale dual problem

The Lagrangian of ( 10) is given by

J C = K i=1 1 2 ψ Ci 2 H + a Ci 2 2 + µ 0 2 ξ Ci 2 2 -C 0 + λ Ci (1 a Ci -1) -γ Ci a Ci -β Ci ξ Ci -y Ci + M a Ci + ψ Ci (M ) , (36)
where µ 0 , β Ci , λ Ci and γ Ci ≥ 0 are the Lagrange multipliers.

The optimality conditions with respect to the primal variables are given by

a * Ci = M β Ci + γ Ci -λ Ci 1 , ψ * Ci = L =1 β Ci, κ(•, m ) , ξ * Ci = 1 µ 0 β Ci . (37)
By replacing these solutions into the Lagrangian in [START_REF] Wang | Centralized collaborative sparse unmixing for hyperspectral images[END_REF], we can derive the dual optimization problem, which is given by

max µ0 max ω C K i=1 ω Ci B C (µ 0 ) ω Ci + c Ci ω Ci - µ 0 2 C 0 subject to γ Ci ≥ 0 , i = 1, . . . , K , (38)
where ω C = ω 1 , . . . , ω K , ω Ci = β Ci , γ Ci , λ Ci is the vector of dual variables, and B C and c Ci are given by

B C (µ 0 ) = - 1 2   K + 1 µ0 I + M M M -M 1 M I -1 -1 M -1 P   , (39)
c Ci = y Ci 0 -1 . ( 40
)

B. The fine scale dual problem

The Lagrangian of ( 35) is given by

J D = N n=1 1 2 ψ n 2 H + µ 1 2 ξ n 2 2 -C 1 + λ n (1 a n -1) + µ 2 2 a n -a Dn 2 2 -C Y + ξ ψ,n 2 + C E + µ 3,n M † (ψ n (M ) -ψ Cn (M ) -ξ ψ,n -β n ξ n -y n + M a n + ψ n (M ) -γ n a n , (41)
where µ 1 , µ 2 , µ 3,n , β n , γ n and λ n are the Lagrange multipliers. Differentiating the optimality conditions with respect to the primal variables and equating the result to zero we obtain

a * n = a Dn + 1 µ 2 M β n + γ n -λ n 1 , ψ * n = L =1 β n, κ(•, m ) - L =1 [µ 3,n ] M † κ(•, m ) , ξ * n = 1 µ 1 β n , ξ * ψ,n = 1 µ 2 µ 3,n , (42)
where [ • ] denotes the -th position of a vector. Substituting the solution to the primal problem in the Lagrangian, we obtain the following dual problem

max µ1,µ2 max ω N n=1 ω n B(µ 1 , µ 2 )ω n + c n ω n - N 2 µ 1 C 1 + µ 2 C Y -µ 2 C E ( 43
)
subject to γ n ≥ 0 , n = 1, . . . , N ,

where ω = ω 1 , . . . , ω N is a vector containing the dual variables, with entries given by ω n = β n , µ 3,n , γ n , λ n . The terms B(µ 1 , µ 2 ) and c n are defined in [START_REF] Moore | On properly positive hermitian matrices[END_REF].

Although [START_REF] Lu | Manifold regularized sparse NMF for hyperspectral unmixing[END_REF] and [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF] being Lagrangian dual problems implies that they are concave with respect to all variables [START_REF] Boyd | Convex optimization[END_REF], they are still nonlinear and thus computationally intensive to solve given their large dimension. However, the cost function of (43) (resp. ( 38)) becomes quadratic when µ 1 and µ 2 (resp. µ 0 ) are fixed. This will allow us to propose an efficient algorithm in Section V-D to solve these problems.

C. Strong duality of the optimization problems

Since optimization problems [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] and ( 35) are non-convex, it does not immediately follows that the Lagrangian duality gap is zero. This means that, unless shown otherwise, the optimal solutions to the dual problems in ( 43) and ( 38) can be different from those of ( 10) and [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF]. Fortunately, building upon results from non-convex optimization in [START_REF] Tuy | Generalized S-lemma and strong duality in nonconvex quadratic programming[END_REF] we can show that strong duality holds for this problem. This is formalized in the following result: Theorem 1. Suppose that the variables µ * 1 and µ * 2 (resp. µ * 0 ) that solve problem (43) (resp. [START_REF] Lu | Manifold regularized sparse NMF for hyperspectral unmixing[END_REF]) are strictly positive. Then, strong duality holds for problem (12) (resp. [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF]).

B(µ 1 , µ 2 ) = - 1 2       K + 1 µ1 I + 1 µ2 M M -K M † 1 µ2 M -1 µ2 M 1 -M † K 1 µ2 I + M † K M † 0 0 1 µ2 M 0 1 µ2 I -1 µ2 1 -1 µ2 1 M 0 -1 µ2 1 1 µ2 P       , ω n =     β n µ 3,n γ n λ n     , c n = -a Dn M + y n -ψ Cn (M ) M † -a Dn 1 a Dn -1 . (44)
Proof. The proof builds upon the results in [START_REF] Tuy | Generalized S-lemma and strong duality in nonconvex quadratic programming[END_REF]Theorem 6]. Due to space limitations, it is relegated to the supplemental material, also available in [START_REF] Borsoi | A blind multiscale spatial regularization framework for kernel-based spectral unmixing[END_REF].

This shows that the proposed Lagrangian relaxation strategy can achieve the same solution to the original problems, which was the case in all our experiments.

D. An efficient solution to the Lagrangian dual problem

In order to devise an efficient algorithm for solving problems ( 38) and ( 43), we first note that the purpose of the maximization with respect to µ 0 , µ 1 and µ 2 is to ensure that the quadratic equality constraints in the primal problems [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] and ( 35) are satisfied. Thus, we will attempt to write the dual problem in an equivalent form that will allow us to exploit this property to obtain a simpler solution.

The optimality conditions for the coarse scale dual problem [START_REF] Lu | Manifold regularized sparse NMF for hyperspectral unmixing[END_REF] with respect to µ 0 are obtained by differentiating the cost function and setting it equal to zero, which gives

g 0 (µ 0 ; ω) = 0 , (45)
where function g 0 is defined as

g 0 (µ 0 ; ω) = 1 µ 2 0 K i=1 β Ci 2 2 -K C 0 . (46)
Similarly, for the fine scale dual problem (43) the optimality conditions with respect to µ 1 and µ 2 will be given by

g 1 (µ 1 , µ 2 ; ω) = 0 , g 2 (µ 1 , µ 2 ; ω) = 0 , (47)
where functions g 1 and g 2 are defined as

g 1 (µ 1 , µ 2 ; ω) = 1 µ 2 1 N n=1 β n 2 2 -N C 1 , g 2 (µ 1 , µ 2 ; ω) = 1 µ 2 2 N n=1 M β n + γ n -λ n 1 2 2 + µ 3,n 2 2 -N C Y -C E . (48)
Let us define the following functions

ω C (µ 0 ) = arg max ω C : γ C i ≥0 K i=1 ω Ci B C (µ 0 )ω Ci + c Ci ω Ci (49)
for the coarse scale problem, and

ω(µ 1 , µ 2 ) = arg max ω : γ n ≥0 N n=1 ω n B(µ 1 , µ 2 )ω n + c n ω n (50)
for the fine scale problem. By substituting (49) in g 0 and (50) in g 1 and g 2 , the optimal µ 0 , µ 1 and µ 2 can be found by solving two systems of equations, one for the coarse scale problem, given by find µ 0 such that g 0 (µ

0 ; ω C ) = 0 , (51)
and another for the fine scale problem, given by

find µ 1 , µ 2 such that g 1 (µ 1 , µ 2 ; ω) = 0 , g 2 (µ 1 , µ 2 ; ω) = 0 , (52)
where we omitted the dependency of functions ω C and ω on µ 0 , µ 1 and µ 2 for notational simplicity. This consists in maximizing the inner optimization problems in [START_REF] Lu | Manifold regularized sparse NMF for hyperspectral unmixing[END_REF] and ( 43) (with respect to ω C or ω) such that the quadratic constraints of the primal problems are satisfied. Although many techniques can be used to solve (51) (resp. ( 52)), one must note that evaluating functions ω C (µ 0 ) (resp. ω(µ 1 , µ 2 )) is computationally expensive. Thus, to have an efficient solution we resort to a bisection strategy, which is a robust algorithm that converges to approximate solutions to these problems with relatively few function evaluations.

Although the solution of (51) using the conventional bisection algorithm is straightforward, the multidimensional case is less clear. Thus, we present it in the remaining of this section. The multidimensional bisection algorithm relies on the Poincaré-Miranda theorem, which states that if a set of multivariate functions change sign in an interval for any of its coordinates, then there is at least one root, common to all such functions, within that interval [START_REF] Galván | The multivariate bisection algorithm[END_REF], [START_REF] Bachrathy | Bisection method in higher dimensions and the efficiency number[END_REF]. This condition can be used to verify whether a given region in the function's domain contains a zero or not.

Thus, by defining a search space and dividing it in two parts along one of the coordinates, we can test to see in which half the root is contained. By performing this operation alternately along each of the function coordinates, we can get arbitrarily close to the root [START_REF] Galván | The multivariate bisection algorithm[END_REF], [START_REF] Bachrathy | Bisection method in higher dimensions and the efficiency number[END_REF]. This procedure is detailed in Algorithm 1. In all our experiments, we ran Algorithm 1 for ten iterations or until the relative variation of the parameters became smaller than a tolerance factor = 0.1. The final Blind Multiscale Unmixing Algorithm for Nonlinear spectral unmixing (BMUA-N) is presented in Algorithm 2.

VI. DETERMINING THE NUMBER OF SUPERPIXELS

A parameter of fundamental importance in the design of the proposed multiscale transform W is the number of superpixels K, or, equivalently, the average size of each superpixel S = N/K. The purpose of the multiscale transform is to Algorithm 1: Bi-dimensional bisection algorithm Input : Functions g1, g2 : R 2 → R. Output: The estimated root (ar, br). 1 Define an initial rectangle containing the root R = {(a1, b1), (a2, b1), (a1, b2), (a2, b2)}, a1 < a2, b1 < b2 ; 2 while Stopping criteria is not satisfied do 3 Compute centers: (ac, bc) = (a2 -a1)/2, (b2 -b1)/2 ; 4 Divide the search space in two and check for the root: R = {(a1, b1), (ac, b1), (a1, b2), (ac, b2)} ; group semantically/spectrally similar pixels, which are then averaged (within each superpixel) to constitute the coarse scale image, capturing spatial correlation and reducing the influence of noise. From this definition, the desired average superpixel size could be intuitively defined as the largest value of S such that the superpixels are still spectrally homogeneous.

In order to evaluate the homogeneity of the superpixels, we consider the distribution of the singular values of the sets of pixels within each superpixel, which are ordered in the form of matrices (matricized). Thus, the j-th matricized superpixel Y j can be written as

Y j = y I1 , . . . , y I |N j | , {I 1 , . . . , I |Nj | } ⊆ N j , (53)
for j = 1, . . . , K. Denote the singular values of Y j by ρ j,1 , ρ j,2 , . . . , ρ j,|Nj | , ordered from the largest to the smallest magnitude. The homogeneity of the j-th superpixel can then be assessed using the ratio between the two largest singular values of Y j , which intuitively evaluates how close Y j is to being a rank-1 matrix. This measure has already been successfully employed to detect heterogeneous superpixels in HS segmentation [START_REF] Yi | Low-dimensional enhanced superpixel representation with homogeneity testing for unmixing of hyperspectral imagery[END_REF]. We then define the average homogeneity of all superpixels Hom(K) as a function of the number K of superpixels as

Hom(K) = 1 K K j=1 |ρ j,1 | |ρ j,2 | , (54)
where we assume that ρ j,2 exists and is nonzero. Thus, K can be selected using the following simple criterion:

K = max j subject to Hom(j) ≥ (1 -ε) max v Hom(v) , K min ≤ j ≤ K max , (55)
where we restrict the number of superpixels to be within a prescribed interval [K min , K max ].

VII. COMPUTATIONAL COMPLEXITY ANALYSIS

The computational complexity of the proposed algorithm depends mainly on the two bisection procedures in Algorithm 2. Each iteration of the bisection method in the coarse domain problem involves solving K quadratic problems (QPs) in L + P + 1 variables (one for each superpixel), whereas each iteration of the bisection method in the original image scale involves solving N QPs in L + 2P + 1 variables (one for each pixel). Since the bisection method reduces the search domain by half at each iteration, it converges linearly [START_REF] Bachrathy | Bisection method in higher dimensions and the efficiency number[END_REF]. Nevertheless, this is sufficient to achieve a reasonable approximation of the optimal parameters µ * 0 , µ * 1 and µ * 2 in relatively few iterations (≤ 10 in our experiments). To see how this compares to other spatially regularized methods, consider for instance the TV-based nonlinear SU algorithm in [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF]. This algorithm employs a variable splitting procedure that leads to an iterative algorithm. At each iteration, N QPs in L+P +1 variables (one for each pixel) and P linear systems in N variables (one for each endmember) must be solved, which is comparable to our method. This illustrates how the proposed separation of the SU problem in two spatial scales maintains a computational complexity that is competitive with other algorithms, even when the estimation of the parameters is considered.

VIII. RESULTS

In this section, we evaluate the performance of the proposed method using both synthetic and real datasets. The BMUA-N is compared with the fully constrained least squares (FCLS), with the unregularized K-Hype [START_REF]Nonlinear estimation of material abundances in hyperspectral images with 1 -norm spatial regularization[END_REF], with the TV-based K-Hype (K-Hype-TV) [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF], with the CDA-NL [START_REF] Halimi | Hyperspectral unmixing in presence of endmember variability, nonlinearity, or mismodeling effects[END_REF] and with the NDU [START_REF] Ammanouil | Nonlinear unmixing of hyperspectral data with vector-valued kernel functions[END_REF] algorithms. The performances were evaluated using the Root Mean Squared Error (RMSE) between the estimated abundance maps (RMSE A ) and between the reconstructed images (RMSE Y ). The RMSE between a true, generic matrix X and its estimate X is defined as

RMSE X = 1 N X X -X 2 F , (56)
where N X denotes the number of elements in the matrix X.

For the proposed method, the noise covariance matrix Σ e was estimated using the residual method described in [START_REF] Roger | Principal components transform with simple, automatic noise adjustment[END_REF], [START_REF] Mahmood | Modified residual method for the estimation of noise in hyperspectral images[END_REF], and the superpixel sizes were selected using the strategy detailed in Section VI, with K min = N/8, K max = N/170 and ε = 0.1. The polynomial kernel described in (6) was used with d = 2 for all kernel-based non-linear SU algorithms.

A. Synthetic data sets

To compare the performance of the different algorithms quantitatively, we created two synthetic datasets with spatially correlated abundance maps, namely Data Cube 0 (DC0), with 70 × 70 pixels, and Data Cube 1 (DC1), with 50 × 50 pixels. Both datasets were constructed using three spectral signatures with 224 bands extracted from the USGS Spectral Library. The synthetic abundance maps are displayed in the first row of Figs. 2 and3. The reflectance values were generated using two nonlinear mixture models, namely, the bilinear mixing model (BLMM), defined as

y n = M a n + P -1 i=1 P j=i+1 a n,i a n,j m i • m j + e n , (57)
with • being the Hadamard product, and the post-nonlinear mixing model (PNMM), defined as where a fixed exponential value of 0.7 has been applied to the LMM result. Finally, white Gaussian noise with signal to noise ratios (SNR) of both 20 and 30dB was added to all datasets. The parameters for each algorithm were either fixed or selected based on a grid search performed for each dataset, with search ranges defined based on those ranges discussed I. The proposed BMUA-N method outperformed the competing algorithms for almost all cases, except for one where its result was very close to that of the TV-based solution. This is despite the fact that the parameters of K-Hype, K-Hype-TV and NDU were selected through a grid search procedure. The abundance maps provided by the nonlinear SU methods for both datacubes are displayed in Figures 2 and3 for illustrative purposes for the case of the BLMM with an SNR of 20dB. The FCLS results were not displayed for the sake of space since they were significantly worse than those of the other algorithms. It can be seen that the BMUA-N results better approximates the ground truth, and even though the K-Hype-TV solution is smoother for DC2 its mean results are farther from the true values. Moreover, although the NDU achieves relatively good abundance reconstructions for the BLMM, the results are noisy and not as good for the PNMM. The reconstruction errors of the nonlinear SU algorithms, also shown in Table I, were similar and significantly lower than those of the FCLS. The execution times of the BMUA-N algorithm, shown in Table II, are about 3.5 times higher than those of K-Hype-TV and half those of the NDU. Thus, the complexity of the BMUA-N is still on the same order of the other state of the art algorithms even though no significant parameter tuning is necessary.

y n = (M a n ) 0.7 + e n , (58)

B. Experiments with real data

For the simulations with real data we consider the Cuprite, the Urban and the Jasper Ridge datasets, which were captured by the AVIRIS instrument and originally had 224 bands. Water absorption and low SNR bands were removed before processing, resulting in 188 bands for the Cuprite image, 162 bands for the Urban image and 198 bands for the Jasper Ridge image. Previous works indicate that 14 endmembers are present at the Cuprite mining field [START_REF] Imbiriba | Low-rank tensor modeling for hyperspectral unmixing accounting for spectral variability[END_REF], [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF], [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF], while the Jasper Ridge HI is known to have four predominant endmemers [START_REF] Imbiriba | Low-rank tensor modeling for hyperspectral unmixing accounting for spectral variability[END_REF]. For the Urban HI, we consider a smaller subscene (shown in Figure 4) containing three endmembers to allow for an easier evaluation.

The reconstructed abundance maps of the nonlinear SU algorithms for the three datasets are presented in Figures 6 and7. For the Cuprite dataset, four endmembers were sewhose distribution could be clearly distinguished in the scene [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF]. The reconstructed abundance maps of the nonlinear SU algorithms are presented in Figure 5, where it can be seen that, except for the case of the CDA-NL, the results for all algorithms are generally compatible and agree with previous studies of this scene [START_REF] Imbiriba | Low-rank tensor modeling for hyperspectral unmixing accounting for spectral variability[END_REF], [START_REF] Borsoi | A data dependent multiscale model for hyperspectral unmixing with spectral variability[END_REF], [START_REF] Nascimento | Vertex Component Analysis: A fast algorithm to unmix hyperspectral data[END_REF]. Nevertheless, a careful analysis reveals that the BMUA-N results, displayed in the bottom row of the figure, show smoother abundance reconstructions without compromising image details and discontinuities. For the Urban HI, the last row of Figure 6 shows that the abundances estimated by the BMUA-N generally contain stronger components for the road, tree and grass endmembers at the regions where these materials appear more prominently, leading to a better separation between the different endmember classes according to Figure 4. The abundance maps estimated by all algorithms for the Jasper Ridge HI, shown in Figure 7, generally agree with the distribution of the corresponding materials observed in Figure 4. A careful analysis of the BMUA-N results shows that it has stronger vegetation and water components when compared to K-Hype and K-Hype-TV. Moreover, although the abundances estimated by CDA-NL and NDU contain slightly clearer results for the road endmember, they also show considerably more confusion between ground and water.

The reconstruction errors for all datasets are shown in Table III. Except for K-Hype, the nonlinear SU algorithms generally achieved much smaller reconstruction errors when compared to the FCLS. Although the results of BMUA-N were slightly higher than the TV-based solution, as can also be observed in the synthetic experiments, small variations in the reconstruction error do not necessarily correlate with better abundance reconstructions. The execution times of the proposed method were again around 3.5 times higher than that of the TV-based solution but only about half of those of the NDU, which indicates that it scales favorably with larger image sizes.

IX. CONCLUSIONS

In this paper, a blind multiscale unmixing strategy was proposed for nonlinear kernel-based mixing models. Based on the concept of a multiscale regularization strategy recently introduced in [START_REF] Borsoi | A fast multiscale spatial regularization for sparse hyperspectral unmixing[END_REF], we were able to efficiently capture image spatial information by splitting the nonlinear mixing process between two image scales, one containing the coarse, lowdimensional image structures and another representing the original image domain. Furthermore, we employed a theorybased statistical framework to devise a consistent strategy to automatically select the regularization parameters of the proposed algorithm and of the multiscale transformation. This resulted in a truly blind (from the parameter setting perspective) multiscale regularization framework. The unmixing problem was formulated using quadratically constrained optimization problems, for which efficient solutions were obtained by exploring their strong duality and a reformulation of their dual representations as root-finding problems. Simulation results with both synthetic and real data indicate that the proposed strategy leads to a consistent performance improvement when compared to the classical Total Variation regularization, even though no parameter adjustment is necessary.

Figure 1 .

 1 Figure 1. Cuprite image (left) and its superpixel decomposition (right).

2 = σ 2 e

 2 ,ψ , n = 1, . . . , N .

5then 7 a2 12 Evaluate

 712 Evaluate g1 and g2 at the four vertices of R ; 6 if the sign of both g1 and g2 is not constant at all vertices of R the other dimension R = {(a1, b1), (a2, b1), (a1, bc), (a2, bc)} ; g1 and g2 at the four vertices of R ;13 if the sign of both g1 and g2 is not constant at all vertices of R then 14 b2 = bc ; (ar, br) = (a2 -a1)/2, (a2 -a1)/2 ; 20 return (ar, br); Algorithm 2: BMUA-N Input : Y , M , σ 2 e,ψ , the number of superpixels K and multiscale decomposition matrix W . Output: The estimated abundance matrix A. 1 Estimate the noise covariance matrix Σe from Y ; 2 Compute the constants C0, C1, CY and CE using equations (22), (20), and (33) ; 3 Compute Y C = Y W ; 4 Find AC by solving (10) using the procedures detailed in Sections (V-A) and (V-D); 5 Compute AD using (11); 6 Find A by solving (35) using the procedures detailed in sections (V-B) and (V-D); 7 return A;

Figure 2 .

 2 Figure 2. Abundance maps estimated by all algorithms for the data cube DC1.

Figure 3 .

 3 Figure 3. Abundance maps estimated by all algorithms for the data cube DC2.

Figure 4 . 1 N N n=1 y n 2 2 .

 412 Figure 4. Subscene of the Urban (left) and Jasper Ridge (right) HIs.

Figure 5 .

 5 Figure 5. Estimated abundance maps for four endmembers of the Cuprite image.

Figure 6 .

 6 Figure 6. Estimated abundance maps for the Urban HI subscene.

  FCLS K-Hype K-Hype-TV CDA-NL NDU BMUA-

Figure 7 .

 7 Figure 7. Estimated abundance maps for the Jasper Ridge HI.

Table I QUANTITATIVE

 I RESULTS FOR DATA CUBES DC1 AND DC2.

			DC1 data cube
			BLMM	PNMM
	SNR	Method	RMSE A	RMSE Y	RMSE A	RMSE Y
		FCLS	0.2587	0.1143	0.1657	0.1038
		K-Hype	0.0575	0.0816	0.0972	0.0769
	20dB	K-Hype-TV CDA-NL	0.0371 0.0730	0.0814 0.0821	0.0800 0.1708	0.0766 0.0789
		NDU	0.0483	0.0736	0.1263	0.0802
		BMUA-N	0.0326	0.0820	0.0730	0.0774
		FCLS	0.2591	0.0836	0.1633	0.0731
		K-Hype	0.0346	0.0258	0.0765	0.0242
	30dB	K-Hype-TV CDA-NL	0.0323 0.0485	0.0258 0.0265	0.0757 0.1625	0.0243 0.0290
		NDU	0.0336	0.0256	0.1208	0.0323
		BMUA-N	0.0325	0.0257	0.0734	0.0243
			DC2 data cube
			BLMM	PNMM
	SNR	Method	RMSE A	RMSE Y	RMSE A	RMSE Y
		FCLS	0.1718	0.1055	0.1554	0.1033
		K-Hype	0.0723	0.0774	0.1165	0.0758
	20dB	K-Hype-TV CDA-NL	0.0557 0.0648	0.0771 0.0780	0.1037 0.1871	0.0758 0.0781
		NDU	0.0518	0.0710	0.1363	0.0794
		BMUA-N	0.0490	0.0771	0.1009	0.0759
		FCLS	0.1714	0.0750	0.1553	0.0731
		K-Hype	0.0503	0.0244	0.0994	0.0240
	30dB	K-Hype-TV CDA-NL	0.0501 0.0491	0.0244 0.0250	0.0992 0.1803	0.0240 0.0289
		NDU	0.0422	0.0243	0.1318	0.0321
		BMUA-N	0.0393	0.0245	0.0902	0.0241

Table II AVERAGE

 II EXECUTION TIME (IN SECONDS) OF THE ALGORITHMS.

		DC1	DC2	Cuprite	Urban	Jasper Ridge
	FCLS	0.59	0.30	11.66	0.34	1.37
	K-Hype	4.54	2.40	43.29	1.42	75.17
	K-Hype-TV	43.27	22.51	405.81	17.72	74.08
	CDA-NL	29.26	15.57	3295.99	27.90	137.63
	NDU	263.30	141.04	2260.45	74.32	483.34
	BMUA-N	144.60	79.03	1380.21	45.21	240.12

Table III RECONSTRUCTION

 III ERRORS (RMSE Y ) FOR THE REAL DATASETS.

This work has been supported by the National Council for Scientific and Technological Development (CNPq) under grants 304250/2017-1, 409044/2018-0, 141271/2017-5 and 204991/2018-8, by the Foundation for Research Support of the State of Rio Grande do Sul (FAPERGS) under grant 19/2551-0001844-4, and by the Brazilian Education Ministry (CAPES) under grant PNPD/1811213. The authors would like to thank Lucas Drumetz and his collaborators for providing part of the data used in the experimental section of the manuscript. R.A. Borsoi is with the Department of Electrical Engineering,

APPENDIX A SUPPLEMENTAL MATERIAL: PROOF OF THEOREM 1

For the sake of brevity, we provide a demonstration for Theorem 1 only for the problem [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF]. The demonstration for problem [START_REF] Ray | Nonlinear spectral mixing in desert vegetation[END_REF] follows the same ideas and is thus straightforward. The first part of the proof consists in using the same ideas as in the representer theorem (presented in Section II) to obtain a finite dimensional representation of [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF] in its primal form. Note that, since H is formed by a class of functions defined as linear combinations of κ(•, m), m ∈ M, we can write any candidate ψ n ∈ H, n = 1, . . . , N for solving [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF] equivalently as

where ψ ⊥ n is orthogonal to m ), = 1, . . . , L, i.e., κ(•, m ), ψ ⊥ n H = 0, and β n, are the linear combination coefficients. This allows us to write each position of the vector ψ n (M ) as

for j = 1, . . . , L. Thus,

Note that ψ Cn (M ) in ( 35) is a constant vector since it was estimated in the previous problem. The objective function can be similarly written as

Since the constraints do not depend on ψ ⊥ n , the second term of ( 62) is irrelevant to the problem, and will be equal to zero for any solution to the optimization problem. Thus, we can rewrite [START_REF] Chen | Nonlinear estimation of material abundances in hyperspectral images with L1-norm spatial regularization[END_REF] as a finite dimensional optimization problem

This derivation was just to show that we can write the primal version of the optimization problem in finite dimension. The finite dimensional Lagrangian is given by

Now, we can write [START_REF] Hall | Common structure of techniques for choosing smoothing parameters in regression problems[END_REF] equivalently by splitting each quadratic equality constraints in two quadratic inequality ones, one being convex and the other non-convex. This results in the following problem:

A is an affine manifold that encapsulates the equality constraints, functions f (x), g i (x) and h i,j (x), for all i and j are quadratic and defined as

Note that the constraints involving g 3 (x) and g 4 (x) are nonconvex, while the remaining ones are all convex. Note also that this problem does not satisfy the Slater condition, what imposes difficulties in the analysis and precludes the consideration of works that rely on this hypothesis (e.g., [START_REF] Ai | Strong duality for the CDT subproblem: a necessary and sufficient condition[END_REF], [START_REF] Bomze | Copositive relaxation beats Lagrangian dual bounds in quadratically and linearly constrained quadratic optimization problems[END_REF]). The Lagrangian of this problem is Since the constraints involving g 1 (x) and g 3 (x), and g 2 (x) and g 4 (x) are linearly dependent, we have

Thus, we can define new variables v 1 = (u 1 -u 3 ) ∈ R and v 2 = (u 2 -u 4 ) ∈ R and rewrite (72) as

where u = [u 5 , . . . , u N P +4 ] . Since we can always find nonnegative u j , j = 1, . . . , 4 satisfying v 1 = (u 1 -u 3 ) and v 2 = (u 2 -u 4 ), the optimization of the Lagrangian in ( 75) is equivalent to the optimization of the original one in (64) (subject to the corresponding nonnegativity constraints). Thus, this means that strong duality of problem [START_REF] Galatsanos | Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[END_REF] implies that problem (63) has strong duality too.

In order to study strong duality of ( 65), we resort to the following theorem, originally stated in [START_REF] Tuy | Generalized S-lemma and strong duality in nonconvex quadratic programming[END_REF].

and assume that the concave function l(u) : u → inf x∈A L(x, u) attains its maximum at a point u * ∈ R N P +4 + such that L(x, u * ) is strictly convex on the set A . More precisely, ∃u * ∈ arg max

where Hess(•) is the Hessian operator. Then, strong duality holds.

In order to apply this to our problem, we note that since h n,p (x) are linear functions, it is sufficient that u 1 -u 3 > 0 and u 2 -u 4 > 0 for L(x, u) to be strictly convex. Moreover, under the change of variables µ 1 = (u 1 -u 3 ) ∈ R and µ 2 = (u 2 -u 4 ) ∈ R, maximization of l(u) in Theorem (2) becomes equivalent to the dual problem (43) (due to the equivalence we showed between ( 41) and ( 75)). Since we assumed that the optimal solutions µ * 1 and µ * 2 to (43) are strictly positive, this means that u * 1 -u * 3 > 0 and u * 2 -u * 4 > 0 and thus L(x, u * ) is strictly convex. Therefore, strong duality holds.

APPENDIX B SUPPLEMENTAL MATERIAL: ALTERNATIVE IMAGE RECONSTRUCTION ERROR METRICS

In addition to the RMSE, we considered three other quality measures typically employed with hyperspectral images to evaluate reconstructions of real data sets: the spectral angle mapper (SAM), the spectral information divergence (SID) [S1], and a spectral-spatial similarity measure that considers the spatial neighborhood of each pixel, named IPD [S2]. The results for the real datasets are shown in Table IV. These results do not lead to a significantly different qualitative conclusion regarding the best solution.