Hyperspectral Band Selection based on Generalized 3DTV and Tensor CUR Decomposition ^††thanks: This research is partially supported by the NSF grant DMS-1941197.

Hyperspectral imaging (HSI) captures spectral information from ground objects across hundreds of narrow spectral bands. This method yields richer visual data than traditional RGB imagery and has found diverse applications in fields like remote sensing, biomedicine, agriculture, and art conservation. However, HSI data is characterized by its high dimensionality, leading to issues like spectral redundancy and increased computational demands. To address these challenges, it is crucial to apply dimensionality reduction and redundancy elimination strategies. The process of hyperspectral band selection is designed to mitigate these issues by choosing a subset of the original spectral bands, ensuring the retention of critical spectral information.

In general, HSI band selection can be implemented in an either supervised or unsupervised manner depending on the availability of label information. In this work, we focus on unsupervised band selection which does not require prior knowledge of labeled data. Many methods have been developed based on local or global similarity comparison. For example, the enhanced fast density-peak-based clustering (E-FDPC) algorithm [1] aims to identify density peaks to select the most discriminative bands, and the fast neighborhood grouping (FNGBS) method [2] utilizes spatial proximity and spectral similarity to group and select representative bands efficiently. Another notable approach is the similarity-based ranking strategy with structural similarity (SR-SSIM) [3], which prioritizes bands based on their contribution to the overall structural similarity, thus preserving the integrity of the spectral-spatial structures in HSI data.

Graph-based unsupervised band selection methods are gaining popularity for their ability to capture underlying graph similarities. Recently, the Marginalized Graph Self-Representation (MGSR) [4] works to generate segmentation of the HSI through superpixel segmentation, thereby capturing the spatial information across different homogeneous regions and encoding this information into a structural graph. To accelerate the processing, the Matrix CUR decomposition-based approach (MHBSCUR) [5] has been developed. This graph-based method stands out for leveraging matrix CUR decomposition, which enhances computational efficiency. However, graph-based strategies entail extra computational overhead to construct graphs. Furthermore, most band selection approaches necessitate converting data into a matrix instead of directly applying algorithms on the HSI tensor. To address these issues, we propose a tensor-based HSI band selection algorithm that employs tensor CUR decomposition alongside a generalized 3D total variation (3DTV), which enhances the smoothness of the traditional 3DTV by applying the $\ell_{p}$-norm to derivatives.

Matrix CUR decompositions are efficient low-rank approximation methods that use a small number of actual columns and rows of the original matrix. Recent literature has extended the matrix CUR decomposition to the tensor setting. Initially the CUR decomposition was extended to tensors through the single-mode unfolding of 3-mode tensors, as discussed in early works [6]. Subsequent research [7], introduced a more comprehensive tensor CUR variant that encompasses all tensor modes. More recently, the terminology has been refined to include specific names for these decompositions, with [8, 9] introducing the terms Fiber and Chidori CUR decompositions. While these initial tensor CUR decomposition methods utilize the mode-$k$ product, a t-CUR decomposition has also been developed under the t-product framework [10].

In this paper, we focus on the t-product based CUR decomposition, which will be utilized to expedite a crucial phase in our proposed algorithm. This approach allows for the preservation of multi-dimensional structural integrity while facilitating a computationally efficient representation of the original tensor. This innovative combination of CUR and 3DTV in our algorithm leverages the inherent multidimensional structure of HSI data, making our method particularly suitable for dealing with large-scale datasets while ensuring robust band selection performance.

The organization of the remainder of this paper is as follows: Fundamental concepts and definitions in tensor algebra together with G3DTV are introduced in Section 2. The proposed band selection method is presented in Section 3. Section 4 demonstrates the performance of the proposed method through several numerical experiments on benchmark remote sensing datasets. Conclusions and future work are discussed in Section 6.

In this section, we introduce fundamental definitions and notation about third-order tensors. Consider ${\mathcal{A}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, ${\mathcal{B}}\in\mathbb{R}^{n_{2}\times l\times n_{3}}$, and the $k$th frontal slice denoted by ${\mathcal{A}}_{k}={\mathcal{A}}(:,:,k)$. We will use $\mathcal{O}$ to denote the zero tensor.

Note that the proximal operator $\operatorname{prox}_{\lambda\left\lVert\cdot\right\rVert_{1}^{p}}({\mathcal{Z}})$ when $p=1,2,3,4$ can be computed using Algorithms 1 and 2 in [11].

In addition, to describe the high order smoothness of 3D images, we propose a novel 3D total variation regularization.

Consider a hyperspectral data tensor ${\mathcal{Y}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ with $n_{3}$ spectral bands, each composed of $n_{1}\times n_{2}$ spatial pixels. To select bands, we aim to decompose ${\mathcal{Y}}$ into a sum of low-rank and spatial-spectral smooth tensor ${\mathcal{B}}$ and a sparse tensor ${\mathcal{S}}$ via the following model

$$\min_{\text{rank}_{t}({\mathcal{B}})\leq r,{\mathcal{S}}}\frac{1}{2}\left\lVert{\mathcal{B}}+{\mathcal{S}}-{\mathcal{Y}}\right\rVert_{F}^{2}+\lambda_{1}\left\lVert{\mathcal{S}}\right\rVert_{1}+\lambda_{2}\left\lVert\mathcal{B}\right\rVert_{G3DTV}.$$

(2)

Here the third term is the G3DTV of $\mathcal{B}$ as defined in Definition 2.11, i.e., $\left\lVert\mathcal{B}\right\rVert_{G3DTV}=\sum_{i=1}^{3}\left\lVert\nabla_{i}{\mathcal{B}}\right\rVert^{p}_{1}$. The parameter $\lambda_{1}$ is the regularization parameter for controlling the sparsity of the outlier tensor $S$, and $\lambda_{2}>0$ controls the spatial-spectral smoothness of ${\mathcal{B}}$. The parameter $p$ in G3DTV is an integer, which is set as 2 throughout our experiments. In order to apply the ADMM framework to minimize (2), we introduce the auxiliary variables ${\mathcal{X}}_{i}$ and rewrite (2) as an equivalent form

	$\displaystyle\min_{\begin{subarray}{c}\text{rank}_{t}({\mathcal{B}})\leq r\\ {\mathcal{S}},{\mathcal{X}}\end{subarray}}$	$\displaystyle\frac{1}{2}\left\lVert{\mathcal{B}}+{\mathcal{S}}-{\mathcal{Y}}\right\rVert_{F}^{2}+\lambda_{1}\left\lVert{\mathcal{S}}\right\rVert_{1}+\lambda_{2}\sum_{i=1}^{3}\left\lVert{\mathcal{X}}_{i}\right\rVert^{p}_{1}$		(3)
		$\displaystyle\text{s.t.}\quad{\mathcal{X}}_{i}=\nabla_{i}{\mathcal{B}}\quad\text{for}\quad i=1,2,3.$

We introduce an indicator function to take care of the constraints. Let $\Pi=\{{\mathcal{X}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}\,|\,\text{rank}_{\text{t-SVD}}({\mathcal{X}})\leq r\}$. The indicator function $\chi_{\Pi}$ is defined as $\chi_{\Pi}({\mathcal{X}})=0$ if ${\mathcal{X}}\in\Pi$ and $\infty$ otherwise. Then the augmented Lagrangian reads as

	$\displaystyle\mathcal{L}=$	$\displaystyle\frac{1}{2}\left\lVert{\mathcal{B}}+{\mathcal{S}}-{\mathcal{Y}}\right\rVert_{F}^{2}+\lambda_{1}\left\lVert{\mathcal{S}}\right\rVert_{1}+\chi_{\Pi}({\mathcal{B}})+\lambda_{2}\sum_{i=1}^{3}\left\lVert{\mathcal{X}}_{i}\right\rVert^{p}_{1}$
		$\displaystyle+\frac{\beta}{2}\sum_{i=1}^{3}\left\lVert\nabla_{i}{\mathcal{B}}-{\mathcal{X}}_{i}+\widetilde{{\mathcal{X}}}_{i}\right\rVert_{F}^{2},$

where $\widetilde{{\mathcal{X}}}$ is a dual variable and $\beta>0$ is the penalty parameter. The resulting algorithm can be described as follows:

$$\left\{\begin{aligned} {\mathcal{B}}\leftarrow&\operatorname*{argmin}_{\text{rank}_{t}({\mathcal{B}})\leq r}\frac{1}{2}\left\lVert{\mathcal{B}}+{\mathcal{S}}-{\mathcal{Y}}\right\rVert_{F}^{2}\\ &+\frac{\beta}{2}\sum_{i=1}^{3}\left\lVert\nabla_{i}{\mathcal{B}}-{\mathcal{X}}_{i}+\widetilde{{\mathcal{X}}}_{i}\right\rVert_{F}^{2},\\ {\mathcal{S}}\leftarrow&\operatorname*{argmin}_{{\mathcal{S}}}\frac{1}{2}\left\lVert{\mathcal{B}}+{\mathcal{S}}-{\mathcal{Y}}\right\rVert_{F}^{2}+\lambda_{1}\left\lVert{\mathcal{S}}\right\rVert_{1},\\ {\mathcal{X}}_{i}\leftarrow&\operatorname*{argmin}_{{\mathcal{X}}_{i}}\lambda_{2}\left\lVert{\mathcal{X}}_{i}\right\rVert^{p}_{1}+\frac{\beta}{2}\left\lVert\nabla_{i}{\mathcal{B}}-{\mathcal{X}}_{i}+\widetilde{{\mathcal{X}}}_{i}\right\rVert_{F}^{2},\\ \widetilde{{\mathcal{X}}}_{i}\leftarrow&\widetilde{{\mathcal{X}}}_{i}+\nabla_{i}{\mathcal{B}}-{\mathcal{X}}_{i},\end{aligned}\right.$$

where ${\mathcal{X}}_{i}$ and $\widetilde{{\mathcal{X}}}_{i}$ are updated for $i=1,2,3$. Applying the ADMM algorithm requires solving three subproblems at each iteration. Specifically, we aim to minimize $\mathcal{L}$ with respect to ${\mathcal{B}}$, ${\mathcal{S}}$, and ${\mathcal{X}}_{i}$. A common approach for updating ${\mathcal{B}}$ is to use the skinny t-SVD, however this can be costly when the size of the tensor is large. To address this issue, we utilize the t-CUR decomposition and update ${\mathcal{B}}$ via

	$\displaystyle{\mathcal{B}}^{j+1}$	$\displaystyle=\operatorname*{argmin}_{\text{rank}_{t}({\mathcal{B}})\leq r}\frac{1}{2}\left\lVert{\mathcal{B}}+{\mathcal{S}}^{j}-{\mathcal{Y}}\right\rVert_{F}^{2}$
		$\displaystyle+\frac{\beta}{2}\sum_{i=1}^{3}\left\lVert\nabla_{i}{\mathcal{B}}-{\mathcal{X}}^{j}_{i}+\widetilde{{\mathcal{X}}}^{j}_{i}\right\rVert_{F}^{2}:=f({\mathcal{B}}).$

By applying gradient descent with step size $\tau>0$, we obtain the updating scheme for $\mathcal{B}$ as ${\mathcal{B}}^{j+1}={\mathcal{B}}^{j}-\tau\nabla f({\mathcal{B}}^{j})$ where the gradient is calculated as $\nabla f({\mathcal{B}}^{j})={\mathcal{B}}^{j}+{\mathcal{S}}^{j}-{\mathcal{Y}}+\beta\sum_{i=1}^{3}\nabla_{i}^{T}(\nabla_{i}{\mathcal{B}}^{j}-{\mathcal{X}}^{j}_{i}+\widetilde{{\mathcal{X}}}^{j}_{i})$. Then the t-CUR decomposition is updated as

$$\left\{\begin{aligned} {\mathcal{C}}\leftarrow&{\mathcal{C}}-\tau\nabla f({\mathcal{B}}^{j})(:,J,:),\\ \mathcal{R}\leftarrow&\mathcal{R}-\tau\nabla f({\mathcal{B}}^{j})(I,:,:),\\ {\mathcal{U}}\leftarrow&\tfrac{1}{2}({\mathcal{C}}+\mathcal{R}),\end{aligned}\right.$$

where $I$ and $J$ are the respective row and column index sets. Next we calculate ${\mathcal{U}}^{\dagger}$ using Definition 2.6

$\displaystyle{\mathcal{B}}^{j+1}={\mathcal{C}}*{\mathcal{U}}^{{\dagger}}*\mathcal{R}.$

(4)

By fixing other variables, we update ${\mathcal{S}}$ and ${\mathcal{X}}_{i}$ as

	$\displaystyle S^{j+1}$	$\displaystyle=\operatorname{prox}_{\frac{\lambda_{1}}{\beta}\left\lVert\cdot\right\rVert_{1}}({\mathcal{Y}}-{\mathcal{B}}^{j+1})$		(5)
	$\displaystyle{\mathcal{X}}_{i}^{j+1}$	$\displaystyle=\operatorname{prox}_{\frac{\lambda_{2}}{\beta}\left\lVert\cdot\right\rVert^{p}_{1}}(\nabla_{i}{\mathcal{B}}^{j+1}+\widetilde{{\mathcal{X}}}^{j}_{i}),\quad i=1,2,3.$		(6)

Here $\operatorname{prox}_{\frac{\lambda_{1}}{\beta}\left\lVert\cdot\right\rVert_{1}}$ is the soft thresholding operator and $\operatorname{prox}_{\frac{\lambda_{2}}{\beta}\left\lVert\cdot\right\rVert^{p}_{1}}$ is the $\ell_{1}^{p}$ proximal operator defined in (1). The convergence condition is defined as: ${\left\lVert{\mathcal{B}}^{j+1}-{\mathcal{B}}^{j}\right\rVert_{\infty}<\varepsilon}$, where $\varepsilon$ is a predefined tolerance. Finally, we apply a classifier such as k-means on ${\mathcal{B}}^{j+1}$ to find the desired $k$ clusters. The fiber indices of the bands closest to the cluster centroids are stored in the set $Q$. Thus the corresponding bands from the original tensor ${\mathcal{Y}}$ represent the desired band subset.

The main computational cost of Algorithm 1 arises from updating ${\mathcal{B}}$ using the t-product. The per iteration complexity is $O(2n_{3}\log(n_{3})+n_{1}n_{2}\max(s_{r},s_{3})n_{3})$. Thus a computational trade-off exists when applying the t-CUR decomposition.

In this section, we discuss parameter selection for each test dataset and provide general guidelines. Table 1 lists the optimal parameter range for each parameter and dataset. As the number of selected bands increases, the optimal choices for $\lambda_{1}$ and $\lambda_{2}$ are stable, whereas parameters $\beta$ and $\tau$ must be more carefully tuned. Importantly, there exist multiple combinations of parameter choices that may lead to optimal overall accuracy.

In addition to the optimal parameter ranges outlined in Table 1, it is important to consider the effects of noise on parameter selection. Hyperspectral data often contains noise, commonly modeled as Gaussian noise, can greatly impact the performance of band selection algorithms. To account for the noise, the parameters may need to be adjusted to accommodate the deviation from the original signal introduced by the noise.

Considering a scenario where Gaussian noise with various standard deviation levels is introduced into the dataset we recommend the following guidelines for selecting 15 bands from the Indian Pines dataset. The parameters $\lambda_{1}$ and $\lambda_{2}$, which control the trade-off between the sparsity term and the 3DTV regularization term, may require adjustment to counteract the introduction of noise. For the noisy case, $\beta=1$ and $\tau=10^{-4}$ consistently yield optimal results across noise levels. Typically, a higher level of noise may necessitate a larger value of $\lambda_{2}$ as it enforces stricter smoothness in ${\mathcal{B}}$. The parameter $\lambda_{1}$ may be chosen slightly higher than typical settings in noise-free environments. The optimal parameters for $\sigma=1,2,3$ are presented in 2.

The inherent high dimensionality and redundancy of HSI data necessitates effective strategies for band selection to ensure computational efficiency and data interpretability. In this work, we propose a novel tensor-based band selection approach that utilizes the G3DTV regularization to preserve high-order spatial-spectral smoothness, and the tensor CUR decomposition to improve processing efficiency while preserving low-rankness. Our method distinguishes itself by maintaining the tensor structure of HSI data, avoiding the often inefficient conversion to matrix form and directly addressing the spectral redundancy problem in the tensor domain. This approach not only preserves the multidimensional structure of the data but also enhances computational scalability, making it ideal for handling large-scale datasets in practical applications. In future work, we aim to further improve the efficiency of our algorithm by developing an accelerated version. Furthermore, we plan to explore various importance sampling schemes in the t-CUR decomposition and study the impact of gradient tensor sparsity on the band selection accuracy.