Underdetermined Blind Source Separation via Weighted Simplex Shrinkage Regularization and Quantum Deep Image Prior

1 Underdetermined Blind Source Separation via W eighted Simple x Shrinkage Re gularization and Quantum Deep Image Prior Chia-Hsiang Lin, Senior Member , IEEE , and Si-Sheng Y oung, Student Member , IEEE Abstract —As most optical satellites remotely acquire multi- spectral images (MSIs) with limited spatial resolution, multi- spectral unmixing (MU) becomes a critical signal processing technology f or analyzing the pure material spectra f or high- precision classification and identification. Unlike the widely in vestigated hyperspectral unmixing (HU) problem, MU is much more challenging as it corresponds to the underdetermined blind source separation (BSS) problem, where the number of sources is lar ger than the number of a vailable multispectral bands. In this article, we transform MU into its overdetermined counterpart (i.e., HU) by in venting a radically new quantum deep image prior (QDIP), which r elies on the virtual band- splitting task conducted on the observed MSI for generating the virtual hyperspectral image (HSI). Then, we perform HU on the virtual HSI to obtain the virtual hyperspectral sources. Though HU is overdetermined, it still suffers fr om the ill-posed issue, for which we employ the con vex geometry structure of the HSI pixels to customize a weighted simplex shrinkage (WSS) regularizer to mitigate the ill-posedness. Finally , the virtual hyperspectral sources are spectrally downsampled to obtain the desired multispectral sources. The proposed geometry/quantum- empower ed MU (GQ- µ ) algorithm can also effectiv ely obtain the spatial abundance distrib ution map for each source, where the geometric WSS regularization is adaptively and automatically controlled based on the sparsity pattern of the abundance tensor . Simulation and real-w orld data experiments demonstrate the practicality of our unsupervised GQ- µ algorithm f or the challenging MU task. Ablation study demonstrates the strength of QDIP , not achieved by classical DIP , and validates the mechanics-inspired WSS geometry regularizer . The associated code will be av ailable at https://github .com/IHCLab/GQ- mu. The Supplementary Material can be found at the provided link. Index T erms— Multispectral unmixing, hyperspectral unmix- ing, multispectral image, underdetermined blind source sep- aration, hybrid quantum-classical computing, quantum deep learning, deep image prior , con vex geometry . I . I N T R O D U C T I O N Simplex-structured matrix/tensor factorization has emer ged as a critical branch in the blind source separation (BSS) area This study w as supported by the Emer ging Y oung Scholar Program (namely , the 2030 Cross-Generation Y oung Scholars Program) of National Science and T echnology Council (NSTC), T aiwan, under Grant NSTC 114-2628-E-006- 002. W e thank the National Center for Theoretical Sciences (NCTS) and the National Center for High-performance Computing (NCHC) for providing the computing resources. (Chia-Hsiang Lin and Si-Sheng Y oung contribute equally to this work.) (Corr esponding author: Chia-Hsiang Lin.) C.-H. Lin is with the Department of Electrical Engineering, Na- tional Cheng Kung Uni versity , T ainan, T aiwan (R.O.C.) (e-mail: chiah- siang.stev en.lin@gmail.com). S.-S. Y oung is with the Institute of Computer and Communication Engineer - ing, Department of Electrical Engineering, National Cheng Kung Univ ersity , T ainan, T aiwan (R.O.C.) (e-mail: q38121509@gs.ncku.edu.tw). due to its wide applicability to solve numerous real-world problems, including the NP-hard spectrum unmixing [1]. In recent years, popular optical remote sensing satellites (e.g., Sentinel-2) capture multispectral images (MSIs) over different landscapes, enabling v arious crucial applications [2]. Ho we ver , most optical satellites remotely acquire MSIs with limited spatial resolution; therefore, multispectral unmixing (MU) becomes a critical signal processing technique to resolve the so-called mixed pixel phenomenon [3], [4]. Although the MU technique can facilitate the analysis of pure material spectra for high-precision material classification and identification, it has been rarely in vestigated compared to the widely studied hyperspectral unmixing (HU) problem [5], [6]. Both MU and HU are BSS problems, aiming at separating the mixed pixel into pure endmember spectra (sources), while MU is much more challenging as it corresponds to the underdetermined BSS. Underdetermined BSS addresses the scenario wherein the number of sources is lar ger than the number of available multispectral bands. In this article, since MU is transformed into its overde- termined counterpart (i.e., HU), we recall some HU tech- niques in the following. Numerous NP-hard HU criteria (e.g., con ve x geometry-based criteria [7] and non-negativ e matrix factorization (NMF)-based criteria [8]) ha ve been proposed in existing studies. Also, deep neural networks (DNN) ha ve attracted significant attention recently , so a variety of DNN- based algorithms have demonstrated their ef fecti veness in the HU task [9], [10]. Existing literature rarely discusses the underdetermined MU problem; the only relev ant work is [11], which tries to solve MU by extending sev eral benchmark HU algorithms, including vertex component analysis (VCA) [12], N-FINDR [13], and NMF [14]. Howe ver , the experiments in [11] consider only the o verdetermined cases, such as a scenario with P := 8 MSI bands and N := 4 materials/sources. In other words, there is no e xisting literature that fundamentally solves the underdetermined MU challenge, as far as we ha ve observed. W e remark that unlike the algebra-based NMF , both VCA and N-FINDR are geometry-based HU methods, and rely on the so-called pure-pixel assumption (PP A) [12], meaning that each source has at least one pix el exclusi vely composed of that source. Considering that PP A is often violated in remote sensing, hyperplane-based Craig-simplex-identification (HyperCSI) [5] algorithm has been proposed to conduct fast HU without relying on PP A, and the Craig criterion behind has theoretically guaranteed source identifiability [15]. T o mitigate the NP-hardness, the polynomial-time conic HU criterion has © 2026 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. This is the accepted version of the manuscript. The final, published version is av ailable at IEEE Xplore via https://doi.org/10.1109/TIP .2026.3673957. 2 also been proposed based on the L ¨ owner-John ellipsoid in recent machine learning literature [16]. Beyond image processing, we briefly revie w some under- determined BSS literature, which are mostly dev eloped for speech signals (thus unsuitable for image-type signals). First, a statistical approach based on W igner-V ille distribution (WVD) is proposed for time-frequenc y signals underdetermined BSS in [17]. The modified mixing matrix estimation is proposed to consider the ne gativ e v alue of the auto WVD of the mixture signal in [17], thereby ef fecti vely resolving the underdeter- mined system. Another representativ e approach is proposed in [18], where singular spectrum analysis is le veraged for those statistically independent sources. Howe ver , the statistically independent assumption does not seem to be true for optical satellite image signals [19]. T o understand the difficulty of the underdetermined MU, we mathematically formulate our tar get problem under the matrix/tensor factorization framew ork. Giv en a P -band MSI Z m ∈ R L 1 × L 2 × P with L = L 1 L 2 pixels, MU seeks to recov er the pure endmember matrix B ∈ R P × N and the three- way abundance tensor S ∈ R L 1 × L 2 × N from the observed MSI (i.e., Z m = S × 3 B ), where N and “ × 3 ” denote the number of sources and the tensor mode- 3 multiplication [defined in the paragraph follo wing (1)], respecti vely . In this work, we focus on the underdetermined scenario, i.e., P < N . This scenario is more practically relev ant, because a typical multispectral satellite image has a much smaller number of bands P when compared with hyperspectral satellite ones. As mentioned, we will tackle this underdetermined MU problem in its counterpart domain (i.e., HU), where we obtain the hyperspectral sources that are subsequently downsampled back to the multispectral sources. Specifically , inspired by the prism embedded within the optical satellite sensor , we de velop an ef fectiv e yet efficient approach to split the MSI into a virtual HSI Z h ∈ R L 1 × L 2 × M , where M > N denotes the number of hyperspectral bands. The MSI is first spectrally split into the hyperspectral domain by the proposed band-splitting (BSP) function, resulting in an auxiliary HSI-like variable with prism attributes. Considering the nonlinearity of the real-world prism, we employ an addi- tional deep regularization, leading to an ef ficient computing of the virtual HSI. Remarkably , this fully unsupervised approach is quite different from the existing spectral super-resolution (SSR) methods, because the mapping of MSI into HSI in typical SSR relies on a large amount of training data [20]. Be- sides, such mapping is adaptiv e for dif ferent real hyperspectral sensors (e.g., cooled CCD camera used in CA VE datasets [21]) while the approach should be general for an y MSIs captured by optical sensors. Conv ersely , we just need a virtual HSI to gain the desired overdetermined system, and hence it does not need to be associated with a real HSI sensor; such a high flexibility enables a general framework to tackle the underdetermined MU problem. W ith the obtained virtual HSI, we perform HU on it to obtain the virtual hyperspectral sources. Although HU is ov erdetermined, it still suffers from ill-posedness, for which we employ the con vex geometry structure of the HSI data pixels to design a weighted simplex shrinkage (WSS) regular - izer to mitigate this issue. Specifically , [15] demonstrates that the true spectral signatures are the v ertices of the minimum- volume data-enclosing simplex under a rather mild sufficient condition (i.e., when the data purity γ > 1 / √ N − 1 ), thereby leading to various computationally efficient quadratic regular- izers for estimating the signatures; for example, the “Center” regularizer in [22] tries to minimize the distance between each signature and the simple x center . The proposed WSS regularizer enforces the signatures to the center but using different weights, which can be adaptively and automatically controlled based on the sparsity pattern of the abundance ten- sor . W e further introduce the advanced quantum deep network (QUEEN) technique [23] rather than classical models to ad- dress the challenging abundance estimation task. Specifically, the depth of classical networks has been theoretically proven to induce a rank-diminishing issue [24]. Since multispectral endmembers are linearly dependent, which may refer to the same abundance in the low-dimensional space (to be shown in Section III-B). The classical deep layer selecting task-rele vant features in a lo w-dimensional space [25], [26] hinders the abundance estimation in underdetermined cases. In contrast, quantum machine learning (e.g., QUEENs) adopts unitary- computing for quantum feature extraction, hence preserving the rank of input data unaltered (thanks to the full rank of unitary quantum operators) [27], [28], and is able to transform the data into higher-dimensional space [29]. The unitary feature learning significantly contrib utes to solving challenging in verse problems [27], [28]. These advantages motiv ate us to incorporate QUEENs to tackle the abundance estimation. More detailed explanation for adopting QUEEN is summarized in Section II-D. T o ensure a fully unsupervised setting, we will implement QUEEN using the concept of deep image prior (DIP) [30] that considers the deep network itself as a regular - izer . Consequently, the aforementioned virtual hyperspectral sources are spectrally downsampled to obtain the desired multispectral sources. W ith a customized MU e xperimental protocol, the effecti veness of the proposed geometry/quantum- empowered MU (GQ- µ ) algorithm will be demonstrated. The remaining parts of this paper are organized as below . The details of the proposed GQ- µ algorithm, including the cri- terion design, implementation tricks, as well as the architecture of the quantum DIP (QDIP), will be collectiv ely presented in Section II. Next, the designed experimental protocol, experi- mental settings, and ev aluations are presented and discussed in Section III. Finally , we conclude this paper in Section IV. Some notations are defined as follows. Z ++ denotes the set of positi ve integers. Gi ven N ∈ Z ++ , we define I N ≜ { 1 , . . . , N } . R N , R M × N , and R M × N × L are the N - dimensional Euclidean space, ( M × N ) -dimensional real- valued matrix space, and ( M × N × L ) -dimensional real- valued 3-way tensor space, respectiv ely . conv S is the con ve x hull of S . 1 (resp., 0 ) denotes the all-one vector/matrix (resp., all-zero matrix) with the dimensionality specified in its subscript. I m is the m × m identity matrix. [ T ] l 1 ,l 2 ,...,l N represents the ( l i , l 2 , . . . , l N ) th entry of the N -way tensor T . M ≥ 0 and T ≥ 0 denote non-negativ e matrix and tensor, respectiv ely . M T and M − 1 denote the matrix transpose and in verse, respectively . vec ( M ) denotes the vector obtained 3 by sequentially stacking the columns M . | · | , ∥ · ∥ 2 , and ∥ · ∥ F are the absolute value, ℓ 2 -norm, and Frobenius norm, respectiv ely . ⊗ is the Kronecker product. Diag ( w 1 , · · · , w N ) denotes diagonal matrix with the n th diagonal entry being w n , ∀ n ∈ I N . I I . T H E P R O P O S E D M U L T I S P E C T R A L U N M I X I N G F R A M E W O R K In this section, we describe how we tackle the challenging underdetermined MU problem. Our approach is outlined as follows. Initially , we propose an unsupervised approach with BSP function and deep regularization for obtaining the virtual HSI (cf. Section II-C). W ith this virtual HSI, we transform the MU problem to its overdetermined counterpart (i.e., HU). Thus, the data-fitting model in HU can be explicitly written as a non-ne gativ e tensor factorization problem for hyperspectral source recovery . Furthermore, as the HU task still suffers from ill-posedness, we design a WSS regularizer to mitigate the ill-posedness. Also, we incorporate the quantum abundance obtained from unsupervised QDIP into the regularization term to address the abundance estimation task. Details on designing the data-fitting and re gularization terms are summarized in Section II-A. Then, the implementations of the ov erall MU algorithm (i.e., GQ- µ ) are presented in Section II-B. Finally , the quantum architecture of QDIP is presented in Section II-D. A. Criterion Design Giv en a three-way MSI tensor Z m , it can be factorized into a non-negati ve endmember matrix B and a non-negati ve three- way ab undance tensor S . Specifically , according to the well- known linear mixing model (LMM) [3], the ( ℓ 1 , ℓ 2 ) th pixel of the MSI [ Z m ] ℓ 1 ,ℓ 2 , : can be expressed as a specific linear combination of the multispectral endmembers, i.e., [ Z m ] ℓ 1 ,ℓ 2 , : = N X n =1 [ S ] ℓ 1 ,ℓ 2 ,n b n , ∀ ( ℓ 1 , ℓ 2 ) ∈ I L 1 × I L 2 , where b n ∈ R P denotes the n th multispectral endmember (i.e., the n th column of B ). W e assume that the model-order N is known a prior for MU or HU [31]. Thus, a natural MU criterion is to minimize the following data-fitting objective, DF ( B , S ) = 1 2 ∥ Z m − S × 3 B ∥ 2 F , (1) where × 3 denotes the tensor mode- 3 multiplication. More precisely, for a tensor T and a matrix M , the ten- sor mode- n multiplication (denoted by T × n M ) is to left-multiply M to each vector along the dimen- sion n within T (i.e., [ T × n M ] i 1 ,...,i n − 1 ,j,i n +1 ,...,i N ≜ P m M j,m [ T ] i 1 ,...,i n − 1 ,m,i n +1 ,...,i N ). It can be observed that (1) refers to an underdetermined system in typical MU ( P < N ), so directly solving (1) is inef fecti ve. A natural idea for solving the underdetermined issue is to gain more observations/equations, probably virtually . The idea is inspired by the recent advances of the spectral super-resolution (SSR) technologies [32], [33] that map the MSI into its hyperspectral counterpart, and we will show that the idea is feasible. According to a typical SSR model, the relation between a MSI and its M -band hyperspectral counterpart Z h ∈ R L 1 × L 2 × M can be modeled as Z m = Z h × 3 D with a spectral response matrix D ∈ R P × M [32]. Moreov er, the hyperspec- tral version of LMM can be formulated as Z h = S × 3 A , where A ≜ [ a 1 , . . . , a N ] ∈ R M × N with a n being the n th hyperspectral endmember . Thus, we have Z m = S × 3 A × 3 D , and hence (1) can be equiv alently reformulated as 1 2 ∥ Z m − S × 3 B ∥ 2 F = 1 2 ∥ Z m − S × 3 A × 3 D ∥ 2 F , (2) which enables us to estimate the hyperspectral signatures A (overdetermined) rather than the multispectral ones B (underdetermined). Accordingly, we incorporate the HSI Z h into this underdetermined system, thereby yielding a new data- fitting model for the MU task, i.e. DF ( A , S ) ≜ 1 2 ∥ Z m − S × 3 A × 3 D ∥ 2 F + 1 2 ∥ Z h − S × 3 A ∥ 2 F . (3) W e note that (3) is similar to a typical unmixing-based fusion problem [34], yet it is fundamentally different. Specifically , the MSI, HSI, spectral response, and spatial downsampling matrix are typically required for data fusion, while our MU framew ork requires only MSI as the input. The HSI Z h here is virtual, which is obtained by prism-inspired BSP function and deep regularization (cf. Section II-C), rather than real hyperspectral sensors. This design allo ws us to characterize their relations using a virtual D without considering resolution variability between Z m and Z h . Thus, Z m and Z h in the virtual framew ork are spatially-aligned and necessarily share the same abundances. Besides, the virtuality of Z h directly provides a deriv able D (directly determined by BSP) without additional estimation, while the spectral response between two real sensors is typically unav ailable (requiring estimation). Con versely , if we already hav e a real hyperspectral sensor with a real D , why would we still need to address the MU problem? In this case, we just have an ov erdetermined system, with no need to solv e the underdetermined MU problem. Thus, adopting a virtual D is actually a natural approach for the MU problem. This practical limitation (unav ailability of real hyperspectral sensors) also makes the typical data-driven SSR methods in v alid for our tar get task. Hence, we obtain Z h by the unsupervised spectral augmentation method, resulting in a highly fle xible, general, and unsupervised SSR stage, as detailed in Section II-C. Through (3), we are able to explore the spatial characteristics from Z m and spectral attributes from Z h , thereby solving the challenging underdetermined MU. Nev ertheless, the data-fitting model alone still suffers from the ill-posed issue, for which we design several regularization schemes to mitigate the ill-posedness of (3). First, based on the intrinsic properties of remote sensing images, a specific material is usually distributed only around certain regions of the image. Thus, a valid abundance tensor is often sparse [35], resulting in the ℓ 1 -norm regularization, i.e., λ 1 ∥ S ∥ 1 ≜ λ 1 N X n =1 L 2 X ℓ 2 =1 L 1 X ℓ 1 =1 | [ S ] ℓ 1 ,ℓ 2 ,n | , (4) where λ 1 ≥ 0 controls the regularization strength. Howe ver , relying solely on the sparsity feature is not suf ficient to 4 Fig. 1. Graphical illustration of the proposed weighted simplex shrinkage (WSS) regularization. Unlike existing volume-minimization types of regular - ization that blindly promote the volume minimization without considering the sparsity pattern of the abundances, the proposed WSS regularization further takes the sparsity level into consideration. Intuiti vely , when an endmember a k n corresponds to a highly sparse abundance map at the k th algorithmic iteration, this endmember is geometrically very far away from the data center c . So, the endmember requires a larger force (i.e., a larger weight w 3 ) to draw it to be closer to c at the ( k + 1) th iteration. By contrast, for the relativ ely dense abundance map (low sparsity), the corresponding endmember just needs to be slightly pushed to the center c , thus requiring a smaller weight (e.g., w 1 ). A larger weight w n will be associated with more orange arrows (denoting a higher force). regularize the challenging problem. Therefore, we further explore quantum features for the regularization of S . Once the quantum abundance feature S qu ∈ R L 1 × L 2 × N is extracted via QDIP (to be done in Section II-D), its information is incorporated into the final ab undance solution S by minimizing this regularizer , λ 2 2 ∥ S − S qu ∥ 2 F , (5) where λ 2 ≥ 0 controls the regularization strength. As dis- cussed in Section II-D, the quantum abundance feature S qu will be computed by QUEEN (rather than traditional deep networks, such as conv olution neural network (CNN)), be- cause QUEEN can extract unique quantum unitary-computing features using very light-weight network architecture (thanks to the quantum full expressibility [36, Theorem 2]). Notably , for the first time, the unsupervised version of QUEEN is employed to solve the matrix/tensor factorization problem, and experimental results evidence the effecti veness of our customized quantum regularization (cf. Section III-C). Furthermore, we propose a geometry-based WSS regular - izer for endmember matrix estimation. It is known that the hyperspectral endmember matrix A = [ a 1 , · · · , a N ] can be perfectly recovered from the vertices of the minimum-volume data-enclosing simplex conv { a 1 , · · · , a N } , under a very mild sufficient condition (i.e., when the data purity γ > 1 / √ N − 1 ) [15]. Therefore, numerous regularization schemes hav e been proposed to minimize the volume of the corresponding end- member simple x con v { a 1 , · · · , a N } . Mathematically , the sim- plex volume can be formulated as 1 ( N − 1)!  det( ¯ A T ¯ A )  1 2 , where ¯ A ≜ [ a 1 − a N , · · · , a N − 1 − a N ] [34]. Howe ver , the inherent non-con ve xity of the determinant-based simplex volume regularization may complicate the gradient descent procedure in NMF frameworks [37]. As an alternative, ex- isting literature adopts con vex surrogates, sum-of-squared distance (SSD) between the endmembers, i.e., ψ ( A ) = 1 2 P N − 1 i =1 P N j = i +1 ∥ a i − a j ∥ 2 2 [34], which enables an efficient implementation. In addition, quadratic minimum simplex vol- ume regularizations [22], [38] hold the general form, 1 2 ∥ AZ − O ∥ 2 F , (6) where Z ∈ R N × N and O ∈ R M × N are some spe- cific matrices corresponding to different surrogates, including “Center”, “Boundary”, and “TV”. For e xample, the “Center” regularization refers to ( Z , O ) ≜ ( I N , c 1 T N ) with c ≜ 1 L P ( ℓ 1 ,ℓ 2 ) [ Z h ] ℓ 1 ,ℓ 2 , : , which minimizes the simplex volume by shrinking the vertices (i.e., the columns of A ) toward the data center c (cf. Figure 1). Moreover, when ( Z , O ) ≜ ( I N − 1 N 1 N × N , 0 M × N ) , the surrogate (6) shrinks the simplex volume by minimizing the total variations between vertices, i.e., “TV” re gularization. Ho wev er , these re gularizations pro- mote volume minimization without considering the shrinkage “force” and “direction”, potentially resulting in less ef fecti ve estimations. Con versely , we design a ne w mechanics-inspired WSS geometry regularizer to address these limitations. First, we conceptually explain the ef fect of the “force”. Intuitiv ely , when an endmember a k n corresponds to a highly sparse ab undance map at the k th algorithmic iteration, the endmember is geometrically far away from the data center c . Thus, this endmember should be forced more to the center at the ( k + 1) th algorithmic iteration, as illustrated in Figure 1. In contrast, for a relativ ely dense abundance map, the corresponding endmember only needs to be slightly pushed to the center c (cf. Figure 1). Moreover, it is essential to consider the shrinkage directions. When the number of algorithmic iterations is large enough, a k n is expected to recov er the n th true simplex vertex. Howe ver , in a typical NMF problem, the endmember and abundance matrices are not jointly conv ex, which implies that the optimal endmember and abundance may not be simultaneously achiev ed. For instance, when the end- member matrix is well recovered at some iteration k (but not the ab undance matrix), the simple x formed by the endmembers tends to be excessi vely shrunk by volume regularizers when further optimizing the abundances in the upcoming algorithmic iterations. Hence, the mechanics-inspired WSS regularization is proposed to simultaneously take the force and direction into account. W e will experimentally demonstrate that the strategy of simultaneously considering force and direction (formulated belo w) does lead to a stronger regularization for the challenging MU task in Section III-C. In practice, we prefer to implement the aforementioned idea using a conv ex function like (6). Thus, we define the weighted simplex shrinkage (WSS) regularizer as WSS ( A ; w ) ≜ λ 3 2 ∥ ( A − c 1 T N ) p Diag ( w ) ∥ 2 F + λ 4 2 ∥ A − A k ∥ 2 F , (7) where λ 3 ≥ 0 and λ 4 ≥ 0 control the corresponding regularization strengths, A k denotes the updated endmember matrix at the k th iteration, and c := 1 N P N i =1 a 0 i (rather than the mean of all hyperspectral pixels, thereby a voiding potential outlier ef fects). Moreov er, the weighting vector is defined as 5 w ≜ softmax ( L ′ 1 , . . . , L ′ N ) , where L ′ n denotes the normalized sparsity lev el for the n th endmember (defined ne xt), and softmax ( · ) is the softmax operator [4]. The sparsity level for the n th endmember can be naturally defined as the reciprocal L n := ∥ [ S 0 ] : , : ,n ∥ − 1 1 = ( P L 2 ℓ 2 =1 P L 1 ℓ 1 =1 | [ S 0 ] ℓ 1 ,ℓ 2 ,n | ) − 1 , where S 0 is the initial abundance estimation (cf. Section II-B). Subsequently, the normalized sparsity level is given by L ′ n := L n / max { L 1 , . . . , L N } . T o hav e an in-depth understanding of our WSS regularization design, one may consider the gradient of our WSS regularization in the case of λ 3 := 1 and λ 4 := 1 , which can be explicitly deriv ed as, ∇ A WSS ( A ; w ) = ( A − c 1 T N ) Diag ( w ) + ( A − A k ) . (8) For the “force”, the first term in the RHS of (8) adaptiv ely shrinks the endmembers to the center based on the sparse pattern w of their corresponding ab undances. A higher sparsity lev el leads to a larger w n (i.e., the n th entry of w ), and hence the first term of (7) will adaptiv ely force a n to the center c , as desired (cf. Figure 1). Additionally , the second term in the RHS of (8) is responsible for re gularizing the “direction” to prevent the ov er-shrinkage or biased issues. Therefore, the second term of (7) encourages a more stable update of the endmember matrix A (not too far away from the previous update A k ) during the optimization procedure. In summary, the proposed regularization strategy is explic- itly summarized as REG ( A , S ) ≜ λ 1 ∥ S ∥ 1 + λ 2 2 ∥ S − S qu ∥ 2 F + WSS ( A ; w ) , resulting in the ov erall criterion, i.e., min A ≥ 0 , S ≥ 0 DF ( A , S ) + REG ( A , S ) , (9) in which the two non-negativ e constraints are to capture the natural properties of the endmembers and abundances. Once the criterion (9) is solved, its optimal solutions ( A ⋆ , S ⋆ ) directly yields the MU solutions, which are the multispectral endmember matrix B ⋆ = D A ⋆ (multispectral sources) and the ab undance tensor S ⋆ . The e xhaustiv e implementations of (9) are summarized in Section II-B. Besides, the details for computing the virtual HSI Z h and for designing the QDIP are presented in Section II-C and Section II-D, respectively . B. Algorithm Implementation In this section, we present the GQ- µ algorithm (i.e., Al- gorithm 1) to implement the criterion (9). In the MU task, we pay more attention to the ab undance tensor S because the abundance explicitly characterizes the distribution of the pure materials. Thus, we first update S k +1 , and then A k +1 . T o this purpose, we initialize A 0 and S 0 by performing HU on the virtual HSI Z h using the highly-mixed/ill-conditioned spectrum unmixing (HISUN) algorithm [16], where Z h will be computed in Section II-C. Subsequently, we use A 0 , Z h , and the proposed QDIP to obtain the quantum abundance tensor (QA T) S qu , and this will be detailed in Section II-D. Once we have S qu , we can then update the S k +1 by solving the first subproblem induced from (9), i.e., the subproblem for Algorithm 1 GQ- µ Algorithm for Solving (9) 1: Gi ven the MSI Z m . 2: Construct the virtual HSI Z h and spectral response matrix D from Z m (cf. Section II-C). Set k := 0 . 3: Initialize A k and S k by conducting HU on Z h using HISUN [16]. 4: Extract QA T S qu from Z h and A k via QDIP (cf. Section II-D). 5: r epeat 6: Update S k +1 by Algotirhm 2. 7: Update A k +1 by (17) using D , A k , S k +1 , Z m , and Z h . 8: k := k + 1 . 9: until the predefined stopping criterion is met. 10: Output multispectral endmember matrix B ⋆ := D A k and abundance tensor S ⋆ := S k . Algorithm 2 ADMM Algorithm for Solving (12) 1: Gi ven Z m , Z h , S qu , D , S k , and A k . Set q := 0 . 2: Initialize Y q := S k and V q := 0 L 1 × L 2 × N 3: r epeat 4: Update S q +1 ∈ arg min S L ( Y q , S , V q ) . 5: Update Y q +1 ∈ arg min Y L ( Y , S q +1 , V q ) . 6: Update V q +1 := V q − ( Y q +1 − S q +1 ) . 7: q := q + 1 . 8: until the predefined stopping criterion is met. 9: Output S k +1 := S q . optimizing S (with fixed A ), which can be deri ved based on (3), (4), (5) and (9) as S k +1 ∈ arg min S ≥ 0 DF ( A k , S ) + λ 1 ∥ S ∥ 1 + λ 2 2 ∥ S − S qu ∥ 2 F , (10) which will be optimized using Algorithm 2 (to be introduced next). After we have S k +1 , we subsequently update A k +1 by solving the second subproblem induced from (9), i.e., the subproblem for optimizing A (with fixed S ), which can be deriv ed from (3), (7) and (9) as A k +1 ∈ arg min A ≥ 0 DF ( A , S k +1 ) + WSS ( A ; w ) . (11) Once the abov e altering optimization con ver ges to ( A , S ) , we hav e the multispectral endmember matrix B ⋆ := D A ⋆ and the abundance tensor S ⋆ as the final MU solution, as summarized in Algorithm 1. Next, we solve (10), which is now conv ex, using the power - ful con ve x optimization solv er known as alternating direction method of multipliers (ADMM) [39]. T o handle the non- differentiable ℓ 1 -norm term in (10), we employ an auxiliary variable Y to reformulate (10) as the standard ADMM form, min S DF ( A k , S ) + λ 1 ∥ Y ∥ 1 + λ 2 2 ∥ S − S qu ∥ 2 F + I + ( S ) , s.t. Y = S , (12) where I + ( · ) is the indicator function of the non-negati ve orthant defined as I + ( S ) ≜ ( 0 , if S ≥ 0 , ∞ , otherwise. 6 According to the ADMM theory , we first write down the augmented Lagrangian of (12), i.e., L ( Y , S , V ) = DF ( A k , S ) + λ 1 ∥ Y ∥ 1 + λ 2 2 ∥ S − S qu ∥ 2 F + I + ( S ) + µ 2 ∥ Y − S − V ∥ 2 F , (13) where µ > 0 and V ∈ R L 1 × L 2 × N are the penalty parameter and the dual variable, respectively . Owing to the dual decom- posability of ADMM, it allows us to solve the non-differential term of Y separately , and the standard ADMM optimization procedure [39] is summarized in Algorithm 2. T o complete Algorithm 2, we need to solv e the two resulting optimization problems presented in Line 4 (for S q +1 ) and Line 5 (for Y q +1 ) of Algorithm 2. For updating S q +1 , the non-negati ve constraint (associated with the indicator function) leads to a slow optimization procedure, if the global optimum is pursued. T o efficiently yet approximately compute S q +1 , we first minimize L ( Y q , S , V q ) (with fixed Y q and V q ) without considering the indicator term of I + , followed by projecting the minimizer back to the non-negati ve orthant. Regarding the initialization of the auxiliary variable Y 0 , we judiciously adopt the previous-stage abundance tensor (i.e., S k ) as a warm start strategy [40], resulting in the e xplicitly updated e xpression, S q +1 ≡ S q +1 (3) = Π(( R + ( λ 2 + µ ) I N ) − 1 ( P + µ T q )) , (14) where R ≜ ( A k ) T D T D A k + ( A k ) T A k , P ≜ ( A k ) T D T Z m (3) + ( A k ) T Z h (3) + λ 2 S qu (3) , T q ≜ Y q (3) − V q (3) , and Π( · ) is the projector onto the non-negati ve orthant. Here, the subscript “ ( n ) ” denotes the tensor mode- n matricization. Specifically , for the N -way tensor G ∈ R M 1 × M 2 ×···× M N , its mode- n matricization is denoted as G ( n ) , and it maps [ G ] ℓ 1 , ··· ,ℓ N to [ G ( n ) ] ℓ n ,j ( j ≜ 1 + P N t =1 ,t  = n ( ℓ t − 1) J t and J t ≜ Q t − 1 m =1 ,m  = n M m ) . T o update Y q +1 in Line 5 of Algorithm 2 (with fixed S q +1 and V q ), we employ the renowned “soft-threshold” (shrinkage) algorithm [41] to solve the ℓ 1 -norm regularization is defined as Y q +1 = arg min Y L ( Y , S q +1 , V q ) = arg min Y λ 1 ∥ Y ∥ 1 + µ 2 ∥ Y − S q +1 − V q ∥ 2 F , = η λ 1 /µ ( S q +1 + V q ) , (15) where η c ( · ) denotes the shrinkage operator . In detail, giv en a three-way tensor T ∈ R L 1 × L 2 × N and a constant c ≥ 0 , the ( ℓ 1 , ℓ 2 , n ) th element of the output three-way tensor η c ( T ) is defined as [ η c ( T )] ℓ 1 ,ℓ 2 ,n ≜      [ T ] ℓ 1 ,ℓ 2 ,n − c, if [ T ] ℓ 1 ,ℓ 2 ,n ≥ c, 0 , if | [ T ] ℓ 1 ,ℓ 2 ,n | ≤ c, [ T ] ℓ 1 ,ℓ 2 ,n + c, if [ T ] ℓ 1 ,ℓ 2 ,n ≤ − c. The deri v ation of Algorithm 2 is thus completed, and conse- quently Line 6 of Algorithm 1 is done. The remaining task is to deriv e Line 7 of Algorithm 1, which is to update the endmember matrix A k +1 by solving (11). T o explicitly derive it, we equiv alently reformulate (11) as min a ≥ 0 1 2 ∥ z m − ( S T ⊗ D ) a ∥ 2 2 + 1 2 ∥ z h − ( S T ⊗ I M ) a ∥ 2 2 + λ 3 2 ∥ ( √ W ⊗ I M )( a − b c ) ∥ 2 2 + λ 4 2 ∥ a − a k ∥ 2 2 , (16) where z m ≜ vec ( Z m (3) ) , S ≜ S ( k +1) (3) , a ≜ vec ( A ) , z h ≜ vec ( Z h (3) ) , W ≜ Diag ( w ) , b c ≜ vec ( c 1 T N ) , and a k ≜ vec ( A k ) . For the non-negati ve constraint of the signa- ture matrix, we adopt the same trick to first ignore it (followed by projecting the solution back to the non-negati ve orthant), thereby leading to an efficient update of (16), i.e., A k +1 ≡ a k +1 = Π( J − 1 1 ( J 2 z m + J 3 z h + λ 3 J 4 b c + λ 4 a k )) , (17) where J 1 ≜ ( S S T ⊗ D T D ) + ( S S T ⊗ I M ) + λ 3 J 4 + λ 4 I N M , J 2 ≜ ( S ⊗ D T ) , J 3 ≜ ( S ⊗ I M ) , and J 4 ≜ ( W ⊗ I M ) . Ac- cording to the abov e implementations, we already completed the proposed GQ- µ algorithm (cf. Algorithm 1). C. Spectr al Augmentation In this section, we describe the approach for obtaining Z h . As emphasized before (cf. Section I), this approach should be general for any MSIs with different spectral responses, which implies that existing data-driv en SSR methods [33] are unsuit- able for our target application. Accordingly , optimization cri- teria with sophisticated regularizations are ine vitable to ensure that Z h preserves the natural properties (e.g., non-negati vity , and spectral-spatial continuity) to be a valid HSI. T o address this challenge, we develop a simple yet effecti ve approach to obtain the virtual HSI. Specifically , the natural properties of real-world HSIs originate from the prism embedded within optical systems, meaning that the HSI can be obtained by further performing the light splitting on its counterpart MSI. Con versely , MSI can be considered as a weighted summation along the spectral dimension of its hyperspectral counterpart. This relation is the fundamental concept of the SSR technique. For example, in an M -time SSR, the green spectral band of an MSI can be approximated by a weighted summation of M hyperspectral bands within the green wav elength range, while these non-negativ e weights are determined by the spectral response function (SRF). Accordingly, we judiciously dev elop the prism-inspired BSP function to obtain a BSP variable e Z h ∈ R L 1 × L 2 × M , which entails the natural properties of HSIs. In other words, the physical meaning of BSP is an unsuper- vised (and virtual) SSR. Therefore, the design of our BSP ensures that the weighted sum of the two split bands is equal to the original bands, while preserving their spectral continuity . More details are presented in the following paragraph. W ith this BSP variable, we transform the challenging SSR problem to a simpler one, i.e., the least-squares problem ∥ Z ′ h − e Z h ∥ 2 F , followed by incorporating a deep regularization (to be defined next) to compensate some non-linear effects during the BSP . Then, we deri ve the virtual HSI from the optimal solution of this deeply regularized problem [cf. (20)]. 7 For this purpose, we first describe the implementation of the BSP function as follows. Given a P -band MSI Z m , the BSP function f BSP ( · ) splits each spectral band within the MSI into τ ∈ Z ++ bands to obtain the M -band e Z h (thus, τ P = M ), i.e., e Z h = f BSP ( Z m ) . T o maintain the computational efficiency while resolving the underdetermined issue, τ can be determined as τ ⋆ := arg min e τ ∈ Z ++ e τ s.t. e τ P ≥ N . (18) In practice, representati ve multispectral satellite images in- clude the 12-band European Space Agency (ESA)’ s Sentinel-2 and the 4-band Systeme Probatoire d’Observ ation de la T erre (SPO T) (i.e., 4 ≤ P ≤ 12 ) [42], and usually the number of materials in a multispectral remote sensing image is about in the range of 4 ≤ N ≤ 7 [11]. Thus, according to (18), e τ := 2 basically guarantees the ov erdetermined property of e τ P ≥ N . Accordingly, we focus on the scenario of τ ⋆ := 2 , for which the two split HSI bands (from a single MSI band [ Z m ] : , : ,q ) can always be represented as 0 . 5[ Z m ] : , : ,q + [ e H ] : , : ,q and 0 . 5[ Z m ] : , : ,q − [ e H ] : , : ,q for some e H ∈ R L 1 × L 2 × P (to be b uilt below) according to the nature of light splitting. Specifically, the r th HSI band can be explicitly written as [ e Z h ] : , : ,r =    0 . 5  [ Z m ] : , : , r +1 2 − [ H ] : , : , r +1 2  , for odd r , 0 . 5  [ Z m ] : , : , r 2 + [ H ] : , : , r 2  , for ev en r, (19) where we define [ H ] : , : ,q := 1 4 ([ Z m ] : , : ,q +1 − [ Z m ] : , : ,q ) , ∀ q ∈ I P − 1 , for spectral continuity . For the boundary case [ H ] : , : ,P , we directly compute [ H ] : , : ,P := 1 4 ([ Z m ] : , : ,P − [ Z m ] : , : ,P − 1 ). Therefore, the spectral response matrix can be determined as D ≜ I P ⊗ 1 T τ to bridge the spectral relation of Z m = Z h × 3 D . As a result, if there is an y negati ve entry in the BSP variable e Z h , it is set to zero to maintain the non-neg ativity of a valid HSI. Note that e Z h is just a rough guess of the virtual HSI. In order to make it closer to the real situation, we refine its structure using the proximal mapping, i.e., Z h = arg min Z ′ h ∥ Z ′ h − e Z h ∥ 2 F + ϕ ( Z ′ h ) , (20) where ϕ promotes some desired structural properties of an image, such as the sparsity or the self-similarity [43]. T o implement (20) efficiently , note that it is nothing but a de- noising problem. Thus, as the plug-and-play (PnP) strategy [44] suggests that ϕ can be implicit, we solv e (20) to obtain the desired virtual HSI Z h , by denoising the image e Z h using the fast and effecti ve denoising CNN (DnCNN) function [45]. W e notice that several DnCNN variants can be selected for our implementation. For the sake of user-friendliness, we employ the default pretrained DnCNN network provided in MA TLAB’ s “Deep Learning T oolbox” to conduct a single- step denoising for resolving (20). Moreover , since this default pretrained DnCNN is a scene-adapti ve denoising network, the manual tuning for the trade-off parameter (i.e., noise lev el) is unnecessary in our implementation. Thus, there is no trade-off parameter in (20). W ith the above strategies, we can ef ficiently and effecti vely obtain Z h . Fig. 2. Graphical illustration of the proposed quantum deep image prior (QDIP), composed of a core quantum module, reshape operator, inv erse quantum collapse (QC) module, and softmax operator. First, we use the core quantum module to obtain the quantum features. Specifically , the module is developed by the quantum full-expressibility (FE) layer “ R X − X X − R Z − X X − R X ”, with trainable real-valued parameters (i.e., α i , β i , δ i , and ρ i ). W ith the provable FE, it can perform any valid quantum operators (cf. Theorem 1). Also, we use the T offoli (i.e., CCNOT) entangled layers for additional entanglement, and we employ the Pauli-Z measurement (cf. Sup- plementary T able I) to obtain the quantum statistics. Once the core quantum module holds the desired quantum states, the inverse-QC module is skillfully lev eraged to prevent the valuable quantum information from being vanished during the reading-out stage. Eventually , the softmax operator guarantees that the output of QDIP satisfies the ab undance sum-to-one property . More detailed configurations are summarized in Supplementary T able II. D. Quantum Deep Ima ge Prior This section presents the architectural design and training strategy of the proposed QDIP to extract the QA T S qu . Before proceeding, we provide the intuition for adopting the quantum network instead of its classical counterpart (e.g., CNN) to facilitate a better understanding. The network here (whether quantum or classical) aims to obtain a approximate solution and, hence, provide an informativ e prior for the proposed GQ- µ [cf. Eq. (5)]. In the underdetermined MU, the multispectral endmembers are linearly dependent (i.e., P < N ), which may refer to the same abundance in the low-dimensional space. For example, the pair of 1st and 6th abundances, the pair of 2nd and 3rd abundances, and the pair of 4th and 5th abundances obtained by the classical DIP (i.e., non-quantum DIP) appear nearly identical, as illustrated in Figure 4. Therefore, an encoder that can search for the desired abundance in the high- dimensional space potentially yields effecti ve estimations. For intuition, considering a typical conv olutional layer with kernel size of 1 × 1 , the relation between a specific C -channel input noise N ∈ R C × L and H -channel output feature F ∈ R H × L can be generally expressed as F = σ ( C N + b 1 T L ) , where C H × C , b H × 1 , and σ ( · ) are the con volutional weights (in matrix-multiplication form), bias vector , and the activ ation function, respectiv ely . T ypically , the default setting of b is an all-zeros v ector , and the ReLU acti vation function is adopted in DIP . When a con volutional network is composed of K layers, it hardly guarantees that all conv olutional kernel matrices remain full-rank during the training procedure, and may result in rank deficiency in classical con volutional networks. While the abov e analysis just sho ws an approximation for kno wing the intuition, similar conclusions of the restrictions in classical networks ha ve also been reported in the literature [24], [26], 8 T ABLE I Q UA N T I TA T I V E E V A L UATI O N S O F T H E B S S A L G O R I T H M S OV E R V A R I O U S S C E NA R I O S I N T E R M S O F ϕ E N ( ↓ ) , ϕ A B ( ↓ ) , A N D R M S E ( ↓ ) . W E U S E B O L D FAC E D A N D U N D E R L I N E D N U M B E R S T O H I G H L I G H T T H E B E S T A N D T H E S E C O N D Q UA N T I TA T I V E P E R F O R M A N C E , R E S P E C T I V E LY . B E S I D E S , T H E C O M P U TA T I O NA L T I M E WAS M E A S U R E D I N S E C O N D S ( S E C . ) . A L L Q UA N T I TA T I V E M E T R I C S A R E R E P O RT E D A S T H E M E A N O F T E N I N D E P E N D E N T R U N S ; T H E C O R R E S P O N D I N G S TA N DA R D D E V I A T I O N S A R E P R E S E N T E D I N S U P P L E M E N TA RY T AB L E I I I . Methods Metrics Boulder Ottawa Haw aii A verage VCA ϕ en ( ↓ ) 24.442 17.450 15.663 19.185 ϕ ab ( ↓ ) 71.227 75.748 54.124 67.033 RMSE ( ↓ ) 35.581 114.802 19.146 56.510 T ime (sec.) 0.009 0.016 0.009 0.011 NMF ϕ en ( ↓ ) 31.424 22.056 27.048 26.843 ϕ ab ( ↓ ) 48.784 44.366 47.464 46.871 RMSE ( ↓ ) 28.306 27.223 25.616 27.048 T ime (sec.) 14.463 14.674 14.463 14.533 MVC-NMF ϕ en ( ↓ ) 29.007 33.669 37.489 33.388 ϕ ab ( ↓ ) 53.958 51.863 62.907 56.243 RMSE ( ↓ ) 26.458 22.790 28.952 26.067 T ime (sec.) 462.671 472.424 465.071 466.722 R-CoNMF ϕ en ( ↓ ) 71.392 63.839 63.460 66.230 ϕ ab ( ↓ ) 50.608 48.348 53.670 50.875 RMSE ( ↓ ) 24.847 21.852 26.597 24.432 T ime (sec.) 44.579 39.475 40.865 41.640 UnDIP ϕ en ( ↓ ) 26.364 3.870 17.701 15.978 ϕ ab ( ↓ ) 55.474 42.573 51.956 50.001 RMSE ( ↓ ) 31.487 13.389 16.021 20.299 T ime (sec.) 57.733 58.001 57.497 57.744 MiSiCNet ϕ en ( ↓ ) 28.693 11.192 19.481 19.789 ϕ ab ( ↓ ) 60.524 49.966 65.380 58.623 RMSE ( ↓ ) 28.944 21.848 31.973 27.588 T ime (sec.) 151.434 140.646 152.659 148.246 D AAN ϕ en ( ↓ ) 19.215 21.446 16.730 19.130 ϕ ab ( ↓ ) 46.528 44.965 47.768 46.420 RMSE ( ↓ ) 25.530 21.591 22.339 23.153 T ime (sec.) 68.521 69.000 68.604 68.708 SSAFNet ϕ en ( ↓ ) 22.979 15.110 19.397 19.162 ϕ ab ( ↓ ) 59.917 41.666 68.099 56.561 RMSE ( ↓ ) 32.336 18.962 34.441 28.580 T ime (sec.) 39.043 80.182 79.030 66.085 GQ- µ ϕ en ( ↓ ) 11.885 5.491 8.191 8.522 ϕ ab ( ↓ ) 26.659 25.215 34.482 28.785 RMSE ( ↓ ) 10.069 8.033 6.477 8.193 T ime (sec.) 51.766 50.194 50.225 50.728 [46]. This explains why classical DIP hardly transfers the problem into an ov erdetermined one. Con versely , quantum machine learning (QML) provides a natural advantage in mapping classical data into a higher - dimensional space (e.g., Hilbert space) [29]. QML also excels at capturing the comple x feature interactions to facilitate more abstract feature representations via quantum entangle- ment (which is also the key behind quantum small-data machine learning) [47]. By integrating these advantages with the pro vable full expressibility (FE, as shown in Theorem 1) of our quantum network, the proposed quantum FE module enables searching for the target state throughout the higher- dimensional superposition quantum space [47] under a single- data setting. The architecture of the QDIP is illustrated in Figure 2, where the detailed configuration is summarized in Supplemen- tary T able II due to space limitations. First, the QDIP consists of the core quantum module, which is dev eloped with prov able FE, followed by the entanglement and measurement layers. The mathematically guaranteed FE ensures that the quantum FE layer can realize an y valid quantum unitary operators (cf. Theorem 1); note that all the valid quantum operators must be unitary . Accordingly , we believe that the core quantum module can represent any desired quantum image state S qu ; this is why the proposed QDIP needs no specific input, just like the con ventional DIP [30]. Once the core quantum module holds the desired quantum state, we use the in verse quantum collapse (in verse-QC) module [36] composed of transposed con volution (TCon v) blocks (cf. Supplementary T able II) to preserve the valuable quantum information from being collapsed during the quantum read-out procedure (cf. Figure 2). Finally , the softmax operator guarantees the output QA T to preserve the abundance sum-to-one nature. For completeness, the ov erall configurations of the QDIP are detailed in Supplementary T able II. Subsequently , we introduce the design of the core quantum module. W e summarize the employed quantum gates along with their corresponding unitary operators in Supplementary T able I. As illustrated in Figure 2, the core quantum module is constructed based on the FE architecture “ R X − X X − R Z − X X − R X ”, i.e., the Rotation-Ising quantum FE layer . T o exploit the quantum entanglement attribute, we also employ the T offoli (i.e., CCNOT) entangled layers for additional entanglement, followed by the Pauli-Z quantum measurement. Notably, this light-weight design of the core quantum module is well aligned with the near-term quantum computers, because ev en some of the most advanced quantum computers (e.g., the 127-qubit IBM Eagle) have very limited qubit resources [48]. The core quantum module, ev en containing just a hundred- scale qubits, can e xpress any valid quantum unitary function: Theorem 1 The trainable quantum neurons deployed in the quantum FE layer embedded within the cor e quantum module (cf. Figur e 2) can expr ess any valid quantum unitary oper- ator U , with some r eal-valued trainable network parameter s { α i , β i , δ i , ρ i } . □ The proof of Theorem 1 is provided in Supplementary Note 1 due to space limitations. Instead, we briefly discuss the key concept of quantum expressibility . F or a specific parameterized quantum circuit (PQC), such as QUEENs, its expressibility is characterized by the ability to e xplore the target quantum states in the quantum space [49, Section 3.1]. In the single-qubit PQC case, the e xpressibility is the capability for exploring the Bloch sphere. A typical quantitati ve assessment of the expressibility is the comparison between the distribution of sampling states and the uniform distribution of states [49]. A higher expressibility may not necessarily yield a better performance, while restricted expressibility indeed results in limited performance. Ho wev er , in vestigating the expressibility of all kinds of PQCs is impractical during architecture develop- ment. Therefore, we de velop the QUEEN with the theoretically guaranteed full expressibility (FE). This design ensures that our QUEEN can realize any v alid quantum operators for the challenging abundance estimation. Next, we discuss the training of QDIP . W e explicitly express the output of QDIP as e S qu := Q θ ( · ) , where Q θ ( · ) denotes the ov erall function of QDIP with θ being the network parameters. 9 Fig. 3. False-color compositions (bands 25, 12, and 6 as RGB) of NASA ’ s HSIs captured by the A VIRIS sensor for different locations and landscapes. As discussed earlier, we can obtain the desired QA T S qu using the proposed QDIP , i.e., θ ⋆ ∈ arg min θ L QDIP ( e S qu ); S qu := Q θ ⋆ ( · ) , (21) where L QDIP and θ ⋆ denote the task-specific loss function and the optimal network parameters of QDIP , respectively . Inspired by [50] and LMM, we design the loss function as L QDIP ( e S qu ) = ∥ Z h − e S qu × 3 A 0 ∥ 2 F + QREG ( e S qu ) , (22) where the regularization term QREG ( · ) is implicitly induced by QDIP network itself (cf. Figure 2) [30], while Z h and A 0 are specified in Section II-C and Section II-B, respectiv ely . The quantum neurons are unitary operators (cf. Supplementary T a- ble I) instead of affine mappings (e.g., conv olution), so QDIP provides unique quantum features for abundance estimation [cf. (5)]. W e employ the adaptive moment estimation (Adam) [51] to optimize (21)-(22), with detailed training settings summarized in Section III-A. For the first time, the quantum v ersion of the seminal idea of DIP [30] is proposed above. For a more comprehensive introduction to quantum deep learning, including the Dirac notation system, barren plateaus (BP) ef fect (a.k.a. quantum gradient vanishing) and quantum measurement/collapse, we refer readers to [36, Section II.A] and [36, Section II.B]. I I I . E X P E R I M E N TA L E V A L U A T I O N In this section, we conduct extensi ve experiments to ev aluate the performance of the proposed GQ- µ . For this purpose, we design an experimental protocol (cf. Section III-A) for the underdetermined MU task due to the lack of reference datasets. Moreov er, the experimental settings (including the proposed GQ- µ and baselines) are also described in Section III-A. Next, the qualitativ e and quantitativ e ev aluations are presented in Section III-B. Section III-C reports ablation studies to substantiate the ef fectiv eness of the proposed QDIP [cf. (5)] and WSS regularization [cf. (7)]. Section III-D presents the real-data analysis to demonstrate the practical applicability of our GQ- µ . Besides, parameter and con ver gence analysis are provided in Supplementary Note 2 due to space limitations. A. Experimental Pr otocol Design and Settings Initially , we observe that the pre vious MU studies do not consider the underdetermined issue; thus, experimental datasets or protocols for MU are una vailable. F or system- atic ev aluation, we design the experimental protocol for this challenging MU task, which is illustrated in Supplementary Figure 1 due to space limitations. Specifically , inspired by W ald’ s protocol for hyperspectral restoration or image fusion [34], we adopt a similar approach to obtain the reference endmember matrix B and abundance tensor S . In this protocol (cf. Supplementary Figure 1), we first perform an effecti ve BSS algorithm (e.g., HISUN [16]) on the reference HSI, thereby obtaining the reference hyperspectral signature A ref and the reference abundance tensor S . Subsequently, we uniformly downsample the reference hyperspectral signature A ref along the spectral dimension (to be described next) to obtain the reference endmember matrix B [52]. Accordingly, we synthesize the reference MSI, i.e., Z m = S × 3 B + e N de , where N de ∈ R L 1 × L 2 × P is the normally distributed LMM deviation tensor; and e := 1E-4 is the scaling factor determined by the energy of the LMM de viation of reference HSI. Once we hav e Z m , the estimated signature b B and abundance b S can be obtained by Algorithm 1 or other benchmarks. Eventually , the qualitativ e and quantitative ev aluations between B and b B (resp. S and b S ) can be readily conducted (cf. Supplementary Figure 1). Moreover, we also experimentally replace N de with a residual-based noise tensor . T o approximate realistic sensor noise characteristics, we first perform spectrally down- sampling on the reference HSI (using the same approach as downsampling A ref ) to approximate the noisy MSI. Then, we apply the pretrained DnCNN denoiser on the noisy MSI and treat the residual as a real-world noise tensor . W ith the quantitativ e metrics (cf. Section III-B) and default parameter settings (to be detailed later), our GQ- µ achieves an average quantitativ e performance of ϕ en = 8 . 249 , ϕ ab = 28 . 768 , and RMSE = 8 . 166 with normally distributed LMM deviation tensor . In contrast, it achieves ϕ en = 7 . 266 , ϕ ab = 28 . 080 , and RMSE = 7 . 385 with residual-based noise tensor . T o implement the designed protocol, we adopt the HSIs captured by N ASA ’ s Airborne V isible/Infrared Imaging Spec- trometer (A VIRIS) sensor [53]. These HSIs are originally composed of 224 contiguous bands, while we remove those water -vapor absorption bands (i.e., bands 1-10, 104-116, 152- 170, and 215-224), followed by using the remaining 172 clean bands as the reference HSIs for the e xperimental protocol (cf. Supplementary Figure 1). For the spectral downsampling stage of A ref , we select the spectral ranges covering 450–520, 520–600, 630–690, and 760–900 nm for synthesizing the MSIs. These w avelength ranges correspond to the Landsat TM satellite bands 1–4 and are consistent with the red- green-blue (RGB) and near-infrared (NIR) bands provided by typical multispectral imagery systems. The NIR band exhibits distinct spectral characteristics compared to RGB [54, T able 3], thereby enabling the material identification of multispectral data to some extent. Furthermore, according to the spectral continuity , the reflectances in these hyperspectral bands are averaged (for each band range) to approximate the multispectral reflectance (of each band) instead of di- rectly choosing/quantizing a specific band. These approxi- mated multispectral reflectances enable us to consider the spectral reflectance of the objects, thereby obtaining the re- alistic reference endmember matrix B ( P = 4 ). Remarkably , ev en if we hav e very fe w MSI bands, the literature often 10 Fig. 4. The qualitative comparisons over the Boulder dataset. The corresponding multispectral endmembers (i.e., M-endmember) of the n th abundance are shown in the n th subfigure within the rightmost column, where the horizontal and vertical axes represent the spectral bands and reflectance values, respectively . used the “multispectral shape” to distinguish objects in, for example, precision agriculture. That said, for MSI, we still hav e the concept of reflectance endmember signatures [11]. In this experiment, we empirically set N := 6 , and the comparisons of the multispectral endmembers within a single R OI are presented in Supplementary Figure 2. Despite the limited number of bands, the multispectral endmembers ex- hibit distinct reflectance with non-overlapping shapes, indicat- ing their material identifiability . Howe ver , when multispectral observations are insufficient for a more meticulous spectral analysis, conv erting them to their corresponding hyperspectral counterparts (cf. Supplementary Figure 3) via customized al- gorithms [55] can provide impro ved practicality . Subsequently , the descriptions of the used HSIs for data simulation are summarized below . For extensi ve ev aluations, we select three HSIs captured ov er different locations (cf. Figure 3); and all datasets have L = 256 × 256 pixels. Specifically , the first dataset (i.e., Boulder) was acquired ov er Boulder , Colorado, with a spatial resolution (SR) of 10.7 m. The second dataset (i.e., Ottawa) was captured over Ottawa, Canada, with a SR of 17.1 m; and the third dataset (i.e., Hawaii) was captured ov er Hawaii County , Hawaii, with a SR of 2.2 m. Next, we describe the benchmarks and their parametric settings. Notably, this is a pioneering study for underdeter- mined MU; thus, there is no e xisting algorithm for this chal- lenging task. Accordingly, we extend several representativ e HU algorithms to compare with the proposed GQ- µ , and the details are described below. First, we adapt the representative approaches, i.e., VCA [12] and NMF [56], to compare with our approach [11]. Specifically , VCA originally performs the dimensional reduction (DR) for HU, while we directly conduct the algorithm in the original multispectral domain because of the underdetermined issue (i.e., P < N ). As VCA only estimates the signature b B , the abundance of VCA is obtained by left-multiplying the pseudo-inv erse of the estimated sig- nature b B † on the MSI (i.e., b S = Z m × 3 b B † ), as suggested in the original VCA paper . In NMF , we adopt HISUN [16] as a warm start strategy to initialize B and S ; and we set the number of iterations to 1E+3 to improve con vergence. After NMF con ver ges, we refine the abundance by solving the non-negati ve least-squares optimization problem of the returned signature and the MSI, which yields an impro ved abundance than the returned one. Moreover , we additionally incorporate two well-known simple x-volume regularized NMF approaches, i.e., minimum volume constrained NMF (MVC- NMF) [37] and robust collaborativ e NMF (R-CoNMF) [38], with spectral augmentation to conduct the underdetermined BSS. In detail, these approaches are also designed for HU and, hence, inevitably require the DR basis matrix for the gradient descent procedure during the alternative optimization (A O). T o address the underdetermined issue (i.e., P < N ), we first apply the BSP function defined in (19), followed by adding a small amount of normally distributed perturbation to resolve the rank deficiency . This yields the full-rank M - band (with M > N ) input for these methods (i.e., MVC- NMF and R-CoNMF). After estimating M -band signatures b A (also the corresponding ab undance b S ) using these methods, we consequently obtain the estimated multispectral endmembers by performing the spectral do wnsampling, i.e., b B = D b A . For the parameter settings, we set the trade-of f parameters τ := 1E-3, tolerance factor of 1E-6, and the maximum iteration number of 5E+2 for MVC-NMF . In R-CoNMF , we 11 Fig. 5. The qualitative comparisons over the Ottawa dataset. The corresponding multispectral endmembers (i.e., M-endmember) of the n th abundance are shown in the n th subfigure within the rightmost column, where the horizontal and vertical axes represent the spectral bands and reflectance values, respectively . set the α := 1E-6 × p ( N / 1000) , β := 5E-1 × ( N / 1000) , λ t := ( N / 1000) , µ t := 1 , θ := N × 1E-3, and the number of A O to 2E+2, which are all the default settings from the official implementation. Beyond the con ventional approaches, we also employ numerous representative and state-of-the-art deep learning approaches, including HU using DIP (UnDIP) [50], minimum simplex con volutional network for deep HU (MiSiCNet) [10], deep autoencoder-based augmented network (D AAN) [9], and spatial-spectral adapti ve fusion network (SSAFNet) [57]; the corresponding experimental settings are described belo w . In UnDIP , we defaultly set the number of iterations to 3E+3; and we set the learning rate (LR) to 1E- 5 for training stability . In MiSiCNet, we set the number of iterations, LR, and the trade-off parameter to 8E+3, 1E-3, and 1E+2, respectively , across all datasets. In DAAN, the LRs of the model encoder and decoder are set to 1E-4 and 1E-2, respectiv ely , under 8E+2 training iterations. In SSAFNet, we set a small LR of 1E-4 to prevent the gradient vanishing, while the trade-of f parameters of the loss function are defaultly set to λ KL := 0 . 001 , λ SAD := 5 , and λ vol := 10 . For the proposed GQ- µ , we first set the LR and the number of iterations of the QDIP (cf. Section II-D) to 5E-2 and 2E+2, respectively , to e xtract S qu ; and the parametric settings of Algorithm 1 are summarized below . First, one can notice that the number of pixels is much larger than the number of sources (i.e., L ≫ N ), and hence we define λ † 3 ≜ λ 3 1E4 and λ † 4 ≜ λ 4 1E4 . Accordingly , we set λ 1 := 1E-3, λ 2 := 1E-2, and λ † 3 := 1E-1. As for λ † 4 , we first initialize it to 1E-2, followed by increasing λ † 4 by 1.2 times in each iteration for a stable optimization. Besides, we set the ADMM penalty parameter µ := 1 for Algorithm 2. Finally , we set the number of iterations to 5 for GQ- µ ; and set it to 20 for Algorithm 2. All benchmarks are conducted on the computational equip- ment equipped with an NVIDIA R TX 3090 GPU, AMD Ryzen 5950X CPU operating at 3.40-GHz, 24-GB RAM, as well as with Core-i7-11700 CPU and 128 GB RAM. The soft- ware en vironments for this experiment are MA TLAB R2023a, Python 3.10.9, and T orch 2.0.0. For the proposed approach, we train the QDIP on the same equipment as benchmarks, and subsequently perform the Algorithm 1 on a laptop with a Core-i7-11700 CPU 2.50-GHz and speed 24-GB RAM. B. Qualitative and Quantitative Evaluations This section presents comprehensiv e qualitative and quan- titativ e ev aluations. First, we graphically illustrate the ref- erence abundances S (resp., estimated abundances b S ) and the reference sources B (resp., estimated sources b B ) for Boulder , Ottawa, and Hawaii datasets in Figure 4, Figure 5, and Supplementary Figure 4 (due to space limitation), respec- tiv ely . In these figures, the n th row presents the qualitative comparisons for the n th abundance and the associated P -band sources, where the color legends are sho wn at the bottom of each abundance column. With these in mind, the qualitativ e ev aluations are presented as follows. As mentioned earlier, we pay more attention to the abun- dances than the endmember signatures in typical MU. Conse- quently, we first conduct qualitativ e analysis for the abundance ov er the Boulder data, with results illustrated in Figure 4. As shown in Figure 4, the proposed GQ- µ yields the most precise abundance (cf. the last abundance column of Figure 4) among all benchmarks. Conv ersely , the benchmarks exhibit unsatisfactory performance for the abundance in the Boul- 12 der data, characterized by indistinguishable abundances for distinct sources, sev ere noise-like patterns, and unexpected structural patterns. For example, this indistinguishability can be clearly observed in the 3rd and 5th abundances of VCA, the 2nd and 3rd abundances of UnDIP , the 1st and 5th abundances of D AAN, and the 1st and 6th ab undances of SSAFNet. W e further note that the qualitativ e performance of the simplex volume-based NMF approaches is severely affected by noise- like patterns; see, e.g., the 1st and 4th ab undances of MVC- NMF , and the 1st, 2nd, and 4th abundances of R-CoNMF . This may be attributed to the fact that, ev en though the spectral augmentation addresses the rank deficiency issue, the theo- retical formulations of MVC-NMF and R-CoNMF primarily consider the hyperspectral structure and fail to exploit the spatial coherence of its multispectral counterpart. In addition, we note that the ab undances of UnDIP have unexpected patterns (e.g., the 4th and 5th abundance of UnDIP), indicating the instability and unreliability of the “black-box” BSS model. Beyond abundance ev aluations, we further conduct qualitativ e comparisons for the estimated multispectral signatures. As shown in the last column of Figure 4, the proposed GQ- µ has numerous signatures well-matched with reference signatures; e.g., see the 1st, 2nd, 3rd, and 6th signatures. In comparison, the signatures estimated from the benchmarks achie ve sev ere shape and intensity discrepancies; e.g., see the 3rd and 6th signatures of R-CoNMF and the 1st signature of D AAN. From the above e valuations, both representati ve NMF and state- of-the-art deep learning HU approaches struggle with MU in the underdetermined scenarios. In contrast, the proposed GQ- µ can effecti vely tackle the underdetermined MU in its hyperspectral counterpart, i.e., HU (cf. Section II-A). Further- more, the proposed WSS and quantum regularizers effecti vely address the ill-posedness of HU. Consequently , the proposed GQ- µ outperforms the baselines on the Boulder dataset. Next, we further validate the effecti veness of our approach on the Ottawa dataset, where the qualitative results are pre- sented in Figure 5. The results indicate that both traditional and deep learning approaches remain limited in terms of performance (e.g., abundance indistinguishability) due to the underdetermined setting. For example, the 1st, 2nd, and 4th abundance of VCA, and the 2nd, 3rd, and 6th abundance of NMF , the 1st and 6th abundance of R-CoNMF , and the 2nd and 3rd abundance of SSAFNet. Although the baseline methods can partially yield relatively accurate estimations for abundance [see Figure 5(d)], the overall performance remains limited. Notably, the undesirable patterns are ag ain observed in the 3rd ab undance of UnDIP . Furthermore, the sum-to-one constraint within UnDIP may further propagate this distortion effect from the 3rd abundance to the 6th abundance. In contrast, the proposed GQ- µ yields substantially promising estimated abundances as shown in the last abundance column of Figure 5, which closely match the reference ab undances. Additionally , the proposed GQ- µ provides numerous precise estimations for signatures on the Ottawa datasets; e.g., see the 1st, 3rd, 4th, and 6th signatures. W e acknowledge that some of the signatures estimated by GQ- µ may occasionally exhibit obvious imperfections; ho we ver , severe shape and intensity deviations rarely occur, such as the 3rd signature of VCA and R-CoNMF . In summary , our GQ- µ yields the best performance compared to the baselines on the Otta wa dataset. Finally, we perform the qualitativ e comparisons for the Hawaii dataset, where the estimated abundance and signatures are illustrated in Supplementary Figure 4 due to space lim- itations. As expected, the proposed GQ- µ yields promising qualitativ e performance in abundance estimations. Nearly all estimated abundances from the proposed GQ- µ are closely matched to the reference abundance. Conv ersely , the tradi- tional benchmarks (i.e., VCA and NMF) still struggle with the underdetermined issue and fail in the abundance estimation. For example, the 1st, 3rd, and 5th abundance of VCA, and the 2nd, 3rd, 4th, and 5th ab undance of NMF . Moreo ver, the noise-like patterns are observed in the estimated abundances from the simplex-volume regularized NMF approaches; e.g., see the 1st, 2nd, and 5th abundance of MVC-NMF and R- CoNMF . Beyond these approaches, the deep learning methods suffer from the underdetermined issue and result in similar abundances for different sources. For instance, the 3rd and 5th abundance of UnDIP , the 2nd, 3rd, 5th abundance of D AAN, and the 1st, 2nd, 3rd, and 5th abundance of SSAFNet. Ac- cording to the abo ve ev aluations, the results strongly support that our GQ- µ achieves the best qualitative performance on all the datasets. Furthermore, we conduct quantitativ e e valuations to fairly ev aluate the results of MU. For this purpose, we first employ a representati ve metric for endmembers ev aluation, i.e., root- mean-square (RMS) spectral angle mapper for endmember ϕ en , whose definition is given in [5, Eq. (32)]. This metric quantitativ ely measures the RMS angle between the reference signatures and estimated signatures. W e further adopt two representativ e metrics for abundance ev aluation, i.e., RMS error (RMSE) and RMS spectral angle mapper for abundance ϕ ab , whose definitions are provided in [50, Eq. (14)] and [5, Eq. (33)], respecti vely . The RMSE and ϕ ab quantitativ ely measure the Euclidean and cosine distance, respecti vely , be- tween the reference and estimated abundance. T o mitigate the randomness of each BSS approach from biasing the quantitativ e ev aluation, the quantitati ve metrics summarized in T able I are reported as the av erage over ten independent runs. Furthermore, the standard deviations of the quantitative metrics are provided in the Supplementary T able III, and the quantitativ e ev aluations are as summarized follows. According to T able I, the proposed GQ- µ achiev es the best quantitativ e performance in both abundance and signatures on the Boulder and Hawaii datasets, aligning with the qualitativ e ev aluations. Furthermore, we note that the quantitative per- formance of the runner -up methods (i.e., UnDIP and DAAN) is substantially inferior to our approach (cf. T able I). This observation indicates that, even with the state-of-the-art deep learning techniques, MU can not be ef fecti vely addressed with- out considering its “underdetermined” nature. In contrast, we effecti vely tackle the underdetermined MU in its hyperspectral counterpart, i.e., HU. Subsequently , the ill-posedness of HU is effecti vely resolved using novel quantum and WSS regular- izers. In the proposed GQ- µ , the WSS regularizer considers both “direction” and ”force” to recover the minimum v olume data-enclosing simplex, instead of simply minimizing the data 13 volume, such as MVC-NMF and R-CoNMF . Additionally , the quantum regularizer effecti vely extracts the informativ e structure to estimate the true abundance rather than using merely the naiv e sparse prior (e.g., ℓ -1 norm regularization). These collectively yield the superior performance (cf. T able I) on the challenging MU applications. On the Ottaw a dataset, the quantitativ e results show that UnDIP is slightly better than our approach in terms of ϕ en , whereas its abundance estimation is unsatisfactory in terms of ϕ ab and RMSE, consistent with the qualitative ev aluations. By comparison, the proposed GQ- µ achieves the best performance in terms of ϕ ab and RMSE on the Ottawa data (cf. T able I). Hence, the proposed GQ- µ still yields the best a verage quantitati ve performance compared to the representative baselines. As a pioneering study of underdetermined MU, we have demonstrated the feasibility of this task based on the abov e comprehensiv e experiments. In the next section, we present an ablation study to e valuate the effecti veness of the proposed quantum and WSS regularizers. C. Ablation Study In this section, we conduct the ablation studies to ex- amine the effecti veness of the quantum and WSS regular - izers (REGs). First, we replace the proposed WSS REG [see (7)] with various types of minimum volume (MV) REGs. Specifically, giv en the hyperspectral endmember ma- trix A = [ a 1 , . . . , a N ] , the simplex volume representation can be expressed as vol ( con ( A )) = 1 ( N − 1)!  det( ¯ A T ¯ A )  1 / 2 with ¯ A ≜ [ a 1 − a N , . . . , a N − 1 − a N ] [34]. Considering the computational complexity of the determinant and potential gradient descent procedure [37], we instead adopt its con vex surrogate, sum-of-squared distances (SSD) REG [34, Eq. (6)], as the first baseline MV REG. W e additionally incorporate three additional quadratic baseline MV REGs referring to “Boundary”, “Center”, and “TV” into this comparison. The explicit formulations of these quadratic MV REGs can be found in [22, T able I]. Finally , by setting all the weights to one (i.e., w := 1 N ), the proposed WSS REG is simplified to its non-weighted counterpart to ev aluate the contribution of the weighting strategy . This baseline MV REG is denoted as “nWSS”. In summary , we comprehensiv ely compare the proposed WSS REG with fi ve different MV REGs, including SSD, Boundary , Center , TV , and nWSS. In addition, we replace the S qu in (5) with the estimated abundance b S from classical DIP , and maintain the same trade-off parameter λ 2 := 1E-2 to v alidate the effecti veness of our QDIP . Here, the widely used UnDIP is selected as the classical DIP . UnDIP first performs the geometry-based endmember extraction, fol- lowed by abundance estimation using a classical con volutional network. Accordingly , the entire implicit model prior of the classical DIP exclusiv ely influences the resulting abundance maps, which aligns with the proposed QDIP and enables a fair comparison in terms of abundance regularization. More importantly , the quantitati ve performance of UnDIP is the second best in ϕ eb and RMSE (cf. T able I). This representative classical DIP will facilitate a persuasi ve ablation study . Subsequently , we present details on the experimental set- tings of this ablation study below. In the follo wing experiment, Fig. 6. True-color composition (left) of the real Sentinel-2 MSI acquired at Northern Alaska, United States; and the 6 estimated multispectral endmembers and their corresponding abundances (right) using the proposed GQ- µ . we select one type of MV REG and DIP to regularize the proposed GQ- µ framew ork, resulting in 12 combinations (6 types of MV REGs and 2 types of DIP). Due to space limita- tions, we summarize the quantitativ e results in Supplementary T able IV , where the selected MV REG/DIP will be marked by “ " ”. Re garding the parameter settings of the trade-of f parameters for MV REGs, since the SSD, Boundary , Center , and TV REGs inv olve only a single term in their formulation, we employ the same λ † 3 := 1E-1 as the proposed GQ- µ for these MV REGs. For the nWSS and WSS REGs, we retain the same settings, λ † 3 := 1E-1 and λ † 4 := 1E-2, as presented in Section III-A. All experiments in this ablation study are based on the av eraged quantitativ e results (i.e., ϕ en , ϕ ab , and RMSE) across the three datasets (see Figure 3). Since our GQ- µ demonstrates stable results across runs (cf. Supplementary T able III), the quantitative results are reported from a single inference. W e now turn to the analysis. First, we employ the classical DIP and compare the effecti veness of the MV REGs, as presented in the first six rows of Supplementary T able IV . In these cases, we observe that the Boundary REG provides the most precise estimations for abundance and signatures. By comparison, the proposed WSS REG appears to be effecti ve only in abundance estimation with the classical DIP; and the rest of the MV REGs seem to achie ve competiti ve performance in terms of ϕ en , ϕ ab , and RMSE. Based on this observation, we experimentally assume that the proposed WSS REG can effecti vely help the NMF procedure with a strong abundance REG (e.g., QDIP) since the shrinkage weights are related to abundance. Accordingly , we replace the classical DIP with the proposed QDIP and again ev aluate the performance of each MV REGs. As demonstrated in the last 6 ro ws of Supplementary T able IV , the results show that the quantitativ e metrics for ab undance (i.e., ϕ ab and RMSE) are significantly improv ed across all MV REGs. These improv ements greatly substantiate the effecti veness and necessity of the proposed QDIP . Moreov er , the proposed GQ- µ (with WSS REG and QIDP) further achieves remarkable improv ements in both abundance and signature estimations. The quantitativ e per- formance substantially outperforms the other 11 cases. These experiments hav e accordingly demonstrated the effecti veness and necessity of our WSS REG and QDIP . D. Real-W orld MU for Sentinel-2 Data This section presents the real-data MU to demonstrate the practical applicability of our GQ- µ . In remote sensing 14 Fig. 7. True-color composition (left) of the real Sentinel-2 MSI acquired at Mackenzie River Delta, Northwest T erritories, Canada; and the 6 estimated multispectral endmembers and their corresponding abundances (right) using the proposed GQ- µ . applications, multispectral sensors typically acquire blue (B), green (G), red (R), and near-infrared (NIR) bands (e.g., SPO T - 7 [42] and RapidEye [58]), facilitating agricultural monitoring, land cover detection, and vegetation index estimation (e.g., ND VI). In this analysis, we utilize MSIs acquired by Sentinel- 2 imagery for cross-dataset validation. Specifically , Sentinel-2 provides 12 bands at multiple resolutions, where the B, G, R, and NIR bands share the highest spatial resolution [43]. W e thus extract these bands as a 4-band MSI with 128 × 128 pix els and perform MU with the proposed GQ- µ . T o preserve spatial consistency , we remove the upsample from the TConv Module 4 of QDIP (cf. Supplementary T able II), while retaining the other settings (cf. Section III-A) for the proposed GQ- µ . Owing to the lack of ground-truth abundances, the estimated abundance is rigorously ev aluated with the high-resolution Google Earth imagery and relev ant geomorphological refer- ences. The follo wing analysis demonstrates that all estimated abundances exhibit consistent spatial coherence with geomor- phological patterns. Furthermore, we validate the estimated endmembers using full Sentinel-2 data and overdetermined endmember extraction (OEE). First, the resolutions of 12 bands are unified by Sentinel-2 super -resolution via scene- adapted self-similarity (SSSS) [43], followed by conducting the OEE using HyperCSI [5]. Subsequently, the B, G, R, and NIR spectra are treated as reference endmembers, which correspond to “Reference (HU)” in Figure 6 and Figure 7. The first real dataset (cf. Figure 6) was acquired over North- ern Alaska (NA), United States, in July 2019. The curvilinear feature in the upper -left region is the Dalton Highway , as identified from Google Earth imagery . The curve of the lower - right area corresponds to the inland branch of Sagav anirktok Riv er . Besides, the distinct white blocks near the riverbank are potentially interpreted as the sedimentary packages [59]. Most of the land surface in this region is composed primar- ily of permanent permafrost, mineral, and vegetation of the tundra. The second real dataset (cf. Figure 7) was acquired ov er Mack enzie River Delta (MRD), Northwest T erritories, Canada, in July 2019. The MRD is notable for the largest concentration of pingos (i.e., small ice-cored hills formed in a piece wise manner) in the w orld [60]. In this R OI, the dark irregular patches correspond to open-water inland lakes, while the yellow patch is identified as an ice-covered water body according to the high-resolution Google Earth imagery . The surrounding uplands are ice-rich, and the unambiguous white patterns likely represent the ice-pure tabular bodies formed at the beginning of ground freezing. The MU results for the NA dataset are illustrated in Figure 6. From left to right, the 2nd, 3rd, and 5th abundances characterize the spatial distributions of the background land cov erage, sedimentary packages, and the main component of the ri ver , respecti vely . The 4th and 6th abundances appear to correspond to mixed substances since this region is composed of multiple substances, including permanent permafrost, min- erals, and ve getation of the tundra. W e further observe the yellow curve in the center of the ri ver (cf. the R OI in Figure 6); the corresponding proportion is more significant in the 1st ab undance instead of the 5th abundance. Although the 1st and 5th abundance seem similar , their endmembers are substantially distinct, indicating that they likely correspond to different materials. Moreover, in the 5th abundance, the locations of sedimentary packages and the yellow curve e xhibit lower proportions, consistent with the 1st and 5th abundance. Since they share close locations (the reflectance may be affected by the water body), it is reasonable that these locations retain a moderate proportion in the 5th abundance. Addition- ally , the multispectral endmembers estimated by our GQ- µ are well aligned with the reference endmembers, implying the precise MU results of GQ- µ . The MU results for the MRD dataset are presented in Figure 7. The 1st and 2nd abundances characterize the spatial distribution of ice-pure tabular bodies and the background grounding, respectively . The 3rd, 4th, and 5th ab undances correspond to the open-water lakes, edges of water bodies, and the ice-covered lakes, respectively . In the abundance of open water , the locations at ice-cov ered lakes retain low proportions since the ice cov erages are floating on water . The estimated endmembers from the proposed GQ- µ are again consistent with those obtained by OEE. The abo ve real MU analysis further supports the effecti veness of the proposed GQ- µ for the challenging MU. The MRD dataset is further used to inv estigate the multi- spectral identifiability . Specifically , the 4-band raw multispec- tral data is classified by K-means clustering with the clustering centers fixed as the MU-based endmember signatures (i.e., the orange curves in Figure 7). As sho wn in Supplementary Figure 6, the clustering results do exhibit a strong spatial coherence with the MRD dataset, alluding that the 4-band multispectral data still hold material identifiability to some extent. Ne v- ertheless, in Supplementary Figure 6, the ice-covered water body (cf. the red box region) and the open-water inland lakes (cf. the blue box regions) are classified into the same cluster , mainly o wing to the limited multispectral representation. In contrast, GQ- µ successfully separates the ice region and open- water lakes into two clearly distinguishable abundance maps (cf. Figure 7 and Supplementary Figure 6). In other words, while the multispectral data alone may provide preliminary material identifiability , our GQ- µ algorithm further improv es the discriminability of homogeneous materials. I V . C O N C L U S I O N S A N D F U T U R E W O R K S As current optical satellites mostly acquire MSIs with lim- ited spatial resolution, multispectral unmixing (MU) has be- come a critical signal processing technology . Hence, this work dev elops an unsupervised geometry/quantum-empowered MU 15 (GQ- µ ) algorithm for solving this underdetermined blind source separation problem (i.e., MU). W e first compute the counterpart HSI from a giv en MSI to transform the problem into a virtual o verdetermined hyperspectral domain via the unsupervised band-splitting (BSP) function. In the BSP [cf. (19)], the spectral deviation v ariable for each splitting band plays a critical role. When it is set to excessiv ely large, the resulting splitting bands may extend beyond their wavelength ranges. Conv ersely , an ov erly small deviation hardly provides sufficient spectral variation and hence f ails to yield informati ve structures. Therefore, the deviation needs to be band-adaptiv e. In our work, we find that setting the de viation as one quarter of the difference between two adjacent bands yields promising results. On the other hand, we inv ent the first quantum version of the seminal deep image prior (DIP) technique to extract specific quantum features to regularize the abundance solution, and tailor a new simplex-geometry regularizer to estimate the hyperspectral endmembers, which are then spectrally down- sampled back to the MSI space to obtain the multispectral endmembers. Ablation study also validates the mechanics- inspired WSS geometry re gularizer . The proposed quantum DIP has mathematically guaranteed FE (ensuring a light- weight QUEEN), and ablation study demonstrates the strength of quantum DIP (not achiev ed by classical DIP). All the above mechanisms can be done under a fully blind setting. In the future, we may further advance the unsupervised spectral augmentation, and hence facilitate a more precise endmember/abundance estimation in MU. In the proposed GQ- µ , the BSP function is designed to provide an initialization of virtual HSI. Then, the proximal mapping, implemented by DnCNN denoising (originally for grayscale images), refines its spatial structure to compute the virtual HSI. A natural improv ement is to consider the spectral-spatial structure si- multaneously , such as incorporating spectral regularization or a spatial-spectral denoiser . Furthermore, inspired by the basis-learning strategy [43], designing a customized update strategy for the virtual HSI during A O can potentially yield a better endmember/ab undance estimation. Regarding the com- plexity of QUEEN, a feasible direction is to explore a new unsupervised or weakly-supervised strategy for pretraining. Once trained, the quantum inference can be formulated as a multiplication of structural matrices (e.g., sparse or unitary) [36], offering a practical strategy for hardware acceleration and efficient quantum gate deployment. R E F E R E N C E S [1] C.-H. Lin, R. W u, W .-K. Ma, C.-Y . Chi, and Y . W ang, “Maximum volume inscribed ellipsoid: A new simplex-structured matrix factoriza- tion framework via facet enumeration and conve x optimization, ” SIAM Journal on Imaging Sciences , vol. 11, no. 2, pp. 1651–1679, 2018. [2] P .-W . T ang and C.-H. Lin, “Transformer-dri ven inv erse problem trans- form for fast blind hyperspectral image dehazing, ” IEEE T rans. Geo- science and Remote Sensing , vol. 62, pp. 1–14, 2024. [3] D. Hong, N. Y okoya, J. Chanussot, and X. X. Zhu, “ An augmented linear mixing model to address spectral v ariability for h yperspectral unmixing, ” IEEE Tr ans. Image Pr ocessing , vol. 28, no. 4, pp. 1923–1938, 2019. [4] S.-S. Y oung, C.-H. Lin, and Z.-C. Leng, “Unsupervised abundance matrix reconstruction transformer-guided fractional attention mechanism for hyperspectral anomaly detection, ” IEEE T rans. Neural Networks and Learning Systems , pp. 1–15, 2024. [5] C.-H. Lin, C.-Y . Chi, Y .-H. W ang, and T .-H. Chan, “ A fast hyperplane- based minimum-volume enclosing simplex algorithm for blind hyper - spectral unmixing, ” IEEE T rans. Signal Processing , vol. 64, no. 8, pp. 1946–1961, 2016. [6] R. A. Borsoi, T . Imbiriba, and P . Closas, “Dynamical hyperspectral unmixing with variational recurrent neural networks, ” IEEE T rans. on Image Pr ocessing , vol. 32, pp. 2279–2294, 2023. [7] S. Henrot, C. Soussen, M. Dossot, and D. Brie, “Does deblurring improve geometrical hyperspectral unmixing?” IEEE Tr ans. Image Pro- cessing , vol. 23, no. 3, pp. 1169–1180, 2014. [8] S. A. V av asis, “On the complexity of nonnegativ e matrix factorization, ” SIAM Journal on Optimization , vol. 20, no. 3, pp. 1364–1377, 2010. [9] Y . Su, Z. Zhu, L. Gao, A. Plaza, P . Li, X. Sun, and X. Xu, “DAAN: A deep autoencoder-based augmented network for blind multilinear hyperspectral unmixing, ” IEEE T rans. Geoscience and Remote Sensing , vol. 62, pp. 1–15, 2024. [10] B. Rasti, B. Koirala, P . Scheunders, and J. Chanussot, “MiSiCNet: Min- imum simplex conv olutional network for deep hyperspectral unmixing, ” IEEE Tr ans. Geoscience and Remote Sensing , vol. 60, pp. 1–15, 2022. [11] J. Cai, H. Chatoux, C. Boust, and et al., “Extending the unmixing methods to multispectral images, ” arXiv Pr eprint arXiv:2111.11893 , 2021. [12] J. Nascimento and J. Dias, “V ertex component analysis: A fast algorithm to unmix hyperspectral data, ” IEEE T rans. Geoscience and Remote Sensing , vol. 43, no. 4, pp. 898–910, 2005. [13] W . Xiong, C.-I. Chang, C.-C. Wu, K. Kalpakis, and H. M. Chen, “Fast algorithms to implement N-FINDR for hyperspectral endmember extrac- tion, ” IEEE Journal of Selected T opics in Applied Earth Observations and Remote Sensing , vol. 4, no. 3, pp. 545–564, 2011. [14] R. Rajabi and H. Ghassemian, “Spectral unmixing of hyperspectral imagery using multilayer NMF, ” IEEE Geoscience and Remote Sensing Letters , vol. 12, no. 1, pp. 38–42, 2014. [15] C.-H. Lin, W .-K. Ma, W .-C. Li, C.-Y . Chi, and A. Ambikapathi, “Identifiability of the simplex volume minimization criterion for blind hyperspectral unmixing: The no-pure-pixel case, ” IEEE T rans. Geo- science and Remote Sensing , vol. 53, no. 10, pp. 5530–5546, 2015. [16] C.-H. Lin and J. M. Bioucas-Dias, “Nonnegati ve blind source separation for ill-conditioned mixtures via John ellipsoid, ” IEEE Tr ans. Neural Networks and Learning Systems , vol. 32, no. 5, pp. 2209–2223, 2021. [17] S. Xie, L. Y ang, J.-M. Y ang, G. Zhou, and Y . Xiang, “Time-frequency approach to underdetermined blind source separation, ” IEEE T rans. Neural Networks and Learn. Sys. , vol. 23, no. 2, pp. 306–316, 2012. [18] H.-G. Ma, Q.-B. Jiang, Z.-Q. Liu, G. Liu, and Z.-Y . Ma, “ A novel blind source separation method for single-channel signal, ” Signal Pr ocessing , vol. 90, no. 12, pp. 3232–3241, 2010. [19] J. M. Nascimento and J. M. Bioucas-Dias, “Does independent compo- nent analysis play a role in unmixing hyperspectral data?” IEEE T rans. Geoscience and Remote Sensing , vol. 43, no. 1, pp. 175–187, 2005. [20] Y . Cai, J. Lin, Z. Lin, H. W ang, Y . Zhang, H. Pfister , R. T imofte, and L. V an Gool, “MST++: Multi-stage spectral-wise transformer for efficient spectral reconstruction, ” in Proc. IEEE Conf. on Computer V ision and P attern Recognition , New Orleans, Louisiana, 2022, pp. 745– 755. [21] F . Y asuma, T . Mitsunaga, D. Iso, and S. Nayar, “Generalized assorted pixel camera: Post-capture control of resolution, dynamic range and spectrum, ” T ech. Rep., Nov 2008. [22] L. Zhuang, C.-H. Lin, M. A. T . Figueiredo, and J. M. Bioucas-Dias, “Regularization parameter selection in minimum volume hyperspectral unmixing, ” IEEE Tr ans. Geoscience and Remote Sensing , vol. 57, no. 12, pp. 9858–9877, 2019. [23] C.-H. Lin and S.-S. Y oung, “HyperKING: Quantum-classical generative adversarial networks for hyperspectral image restoration, ” IEEE T rans. on Geoscience and Remote Sensing , vol. 63, pp. 1–19, 2025. [24] R. Feng, K. Zheng, Y . Huang, D. Zhao, M. Jordan, and Z.-J. Zha, “Rank diminishing in deep neural networks, ” Advances in Neural Information Pr ocessing Systems , vol. 35, pp. 33 054–33 065, 2022. [25] L. Lu, P . Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, “Learning nonlinear operators via deeponet based on the univ ersal approximation theorem of operators, ” Natur e Machine Intelligence , vol. 3, no. 3, pp. 218–229, 2021. [26] A. Ansuini, A. Laio, J. H. Macke, and D. Zoccolan, “Intrinsic dimension of data representations in deep neural networks, ” Advances in Neural Information Processing Systems , vol. 32, 2019. [27] C.-H. Lin, P .-W . T ang, and A. R. Huete, “Quantum feature-empowered deep classification for fast mangrove mapping, ” IEEE T rans. on Geo- science and Remote Sensing , vol. 63, pp. 1–13, 2025. 16 [28] C.-H. Lin, T .-H. Lin, and J. Chanussot, “Quantum information- empowered graph neural network for hyperspectral change detection, ” IEEE Tr ans. Geoscience and Remote Sensing , vol. 62, pp. 1–15, 2024. [29] M. C. Caro, H.-Y . Huang, M. Cerezo, K. Sharma, A. Sornborger, L. Cincio, and P . J. Coles, “Generalization in quantum machine learning from few training data, ” Natur e Comm. , vol. 13, no. 1, p. 4919, 2022. [30] D. Ulyanov , A. V edaldi, and V . Lempitsky , “Deep image prior, ” in Pr oc. IEEE Conference on Computer V ision and P attern Recognition , Salt Lake City , UT , USA, Jun. 18-23, 2018, pp. 9446–9454. [31] C.-H. Lin, C.-Y . Chi, L. Chen, D. J. Miller , and Y . W ang, “Detection of sources in non-negati ve blind source separation by minimum description length criterion, ” IEEE T rans. Neural Networks and Learning Systems , vol. 29, no. 9, pp. 4022–4037, 2018. [32] R. Hang, Q. Liu, and Z. Li, “Spectral super-resolution network guided by intrinsic properties of hyperspectral imagery , ” IEEE T rans. Ima ge Pr ocessing , vol. 30, pp. 7256–7265, 2021. [33] C.-H. Lin, S.-H. Huang, T .-H. Lin, and P . C. W u, “Metasurface- empowered snapshot hyperspectral imaging with conv ex/deep (CODE) small-data learning theory , ” Nature Comm. , vol. 14, no. 1, p. 6979, 2023. [34] C.-H. Lin, F . Ma, C.-Y . Chi, and C.-H. Hsieh, “ A conv ex optimization- based coupled nonnegativ e matrix factorization algorithm for hyperspec- tral and multispectral data fusion, ” IEEE T rans. Geoscience and Remote Sensing , vol. 56, no. 3, pp. 1652–1667, 2017. [35] J. B. Greer , “Sparse demixing of hyperspectral images, ” IEEE T rans. Image Pr ocessing , vol. 21, no. 1, pp. 219–228, 2012. [36] C.-H. Lin and Y .-Y . Chen, “HyperQUEEN: Hyperspectral quantum deep network for image restoration, ” IEEE T rans. Geoscience and Remote Sensing , vol. 61, pp. 1–20, 2023. [37] L. Miao and H. Qi, “Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization, ” IEEE T rans. on Geoscience and Remote Sensing , v ol. 45, no. 3, pp. 765–777, 2007. [38] J. Li, J. M. Bioucas-Dias, A. Plaza, and L. Liu, “Robust collaborative nonnegati ve matrix factorization for hyperspectral unmixing, ” IEEE T rans. on Geoscience and Remote Sensing , vol. 54, no. 10, pp. 6076– 6090, 2016. [39] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein et al. , “Distributed optimization and statistical learning via the alternating direction method of multipliers, ” F oundations and T r ends® in Machine learning , vol. 3, no. 1, pp. 1–122, 2011. [40] C.-Y . Chi, W .-C. Li, and C.-H. Lin, Con vex Optimization for Signal Pr ocessing and Communications: F r om Fundamentals to Applications . Boca Raton, FL, USA: CRC Press, 2017. [41] D. Donoho, “De-noising by soft-thresholding, ” IEEE Tr ans. Information Theory , vol. 41, no. 3, pp. 613–627, 1995. [42] ESA, “SPO T -7 Instruments Details, ” https://earth.esa.int/eog ateway/ missions/spot- 7#instruments- section, Accessed: Sep. 4, 2025. [43] C.-H. Lin and J. M. Bioucas-Dias, “ An explicit and scene-adapted definition of conv ex self-similarity prior with application to unsuper- vised Sentinel-2 super-resolution, ” IEEE Tr ans. Geoscience and Remote Sensing , vol. 58, no. 5, pp. 3352–3365, 2020. [44] C.-H. Lin and S.-S. Y oung, “Signal subspace identification for incom- plete hyperspectral image with applications to various in verse problems, ” IEEE Tr ans. Geoscience and Remote Sensing , vol. 62, pp. 1–16, 2024. [45] K. Zhang, W . Zuo, Y . Chen, D. Meng, and L. Zhang, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising, ” IEEE Tr ans. Image Pr ocessing , vol. 26, no. 7, pp. 3142–3155, 2017. [46] S. Recanatesi, M. Farrell, M. Advani, T . Moore, G. Lajoie, and E. Shea- Brown, “Dimensionality compression and expansion in deep neural networks, ” arXiv Pr eprint arXiv:1906.00443 , 2019. [47] Z. W ang, F . W ang, L. Li, Z. W ang, T . van der Laan, R. C. Leon, J.- K. Huang, and M. Usman, “Quantum kernel learning for small dataset modeling in semiconductor fabrication: Application to ohmic contact, ” Advanced Science , p. e06213, 2024. [48] J. Gambetta, “IBM’ s roadmap for scaling quantum technology , ” IBM Resear ch Blog (September 2020) , 2020. [49] S. Sim, P . D. Johnson, and A. Aspuru-Guzik, “Expressibility and entan- gling capability of parameterized quantum circuits for hybrid quantum- classical algorithms, ” Advanced Quantum T echnologies , vol. 2, no. 12, p. 1900070, 2019. [50] B. Rasti, B. Koirala, P . Scheunders, and P . Ghamisi, “UnDIP: Hyper- spectral unmixing using deep image prior , ” IEEE Tr ans. Geoscience and Remote Sensing , vol. 60, pp. 1–15, 2021. [51] D. P . Kingma and J. Ba, “Adam: A method for stochastic optimization, ” in Pr oc. International Conference for Learning Representations , San Diego, CA, USA, May . 7-9, 2015. [52] C.-H. Lin and J.-T . Lin, “PRIME: Blind multispectral unmixing using virtual quantum prism and conv ex geometry , ” TGRS , 2025. [53] N ASA Jet Propulsion Laboratory, “Airborne V isible/Infrared Imaging Spectrometer Data Portal, ” https://aviris.jpl.nasa.gov/dataportal/. [54] W .-K. Baek, M.-J. Lee, and H.-S. Jung, “Land cover classification from RGB and NIR satellite images using modified U-net model, ” IEEE Access , vol. 12, pp. 69 445–69 455, 2024. [55] C.-H. Lin, J.-T . Chen, Z.-C. Leng, and J.-T . Lin, “COS2A: Conv ersion from Sentinel-2 to A VIRIS hyperspectral data using interpretable algo- rithm with spectral–spatial duality , ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 63, pp. 1–16, 2026. [56] R. Hennequin, B. David, and R. Badeau, “Beta-divergence as a subclass of Bregman diver gence, ” IEEE Signal Pr ocessing Letters , vol. 18, no. 2, pp. 83–86, 2011. [57] W . Gao, J. Y ang, Y . Zhang, Y . Akoudad, and J. Chen, “SSAF-Net: A spatial-spectral adaptive fusion network for hyperspectral unmixing with endmember variability , ” IEEE T rans. on Geoscience and Remote Sensing , vol. 63, pp. 1–15, 2025. [58] ESA Planet Laboratory, “RapidEye Mission Parameters, ” https://earth. esa.int/eogatew ay/missions/rapideye, Accessed: Sep. 4, 2025. [59] R. Gillis, P . Decker , M. W artes, A. Loveland, and T . Hubbard, “Geologic map of the south-central sagav anirktok quadrangle, north slope, alaska, ” Report of Investigations , p. 4, 2014. [60] C. R. Burn, “The mackenzie delta: an archetypal permafrost landscape, ” in Geomorphological Landscapes of the W orld . Springer , 2009, pp. 1–12. Chia-Hsiang Lin (S’10-M’18-SM’24) received the B.S. degree in electrical engineering and the Ph.D. degree in communications engineering from Na- tional Tsing Hua University (NTHU), T aiwan, in 2010 and 2016, respectiv ely . From 2015 to 2016, he was a V isiting Student of V irginia T ech, Arlington, V A, USA. He is currently a Professor with the Department of Electrical Engineering, National Cheng Kung Uni- versity (NCKU), T aiwan. Before joining NCKU, he held research positions with The Chinese University of Hong K ong, HK (2014 and 2017), NTHU (2016-2017), and the Univ ersity of Lisbon (ULisboa), Lisbon, Portugal (2017-2018). He was an Assistant Professor with the Center for Space and Remote Sensing Research, National Central University , T aiwan, in 2018, and a V isiting Professor with ULisboa, in 2019. His research interests include network science, quantum computing, con vex geometry and optimization, blind signal processing, and imaging science. Dr . Lin receiv ed the Emerging Y oung Scholar A ward (The 2030 Cross- Generation Program) from National Science and T echnology Council (NSTC), from 2023 to 2027, the Future T echnology A ward from NSTC, in 2022, the Outstanding Y outh Electrical Engineer A ward from The Chinese Institute of Electrical Engineering (CIEE), in 2022, the Best Y oung Professional Member A ward from IEEE T ainan Section, in 2021, the Prize Paper A ward from IEEE Geoscience and Remote Sensing Society (GRS-S), in 2020, the T op Perfor- mance A ward from Social Media Prediction Challenge at ACM Multimedia, in 2020, and The 3rd Place from AIM Real W orld Super-Resolution Challenge at IEEE International Conference on Computer V ision (ICCV), in 2019. He receiv ed the Ministry of Science and T echnology (MOST) Y oung Scholar Fellowship, together with the EINSTEIN Grant A ward, from 2018 to 2023. In 2016, he was a recipient of the Outstanding Doctoral Dissertation A ward from the Chinese Image Processing and Pattern Recognition Society and the Best Doctoral Dissertation A ward from the IEEE GRS-S. Si-Sheng Y oung (S’23) is currently a Ph.D. student with the Intelligent Hyperspectral Computing Lab- oratory (IHCL), National Cheng Kung University (NCKU), T ainan, T aiwan. In 2023 to 2024, he received the Merit A ward from The Grand Challenge “Computing for the Fu- ture”, Miin W u School of Computing, NCKU, “Pan W en Y uan Scholarship” from Industrial T echnology Research Institute, Hsinchu, T aiwan, and “Scholar- ship Pilot Program to Culti vate Outstanding Doctoral Students” from National Science and T echnology Council (NSTC), T aipei, T aiwan. His research interests include conv ex op- timization, deep learning, anomaly detection, and imaging in verse problems.