# Hyperspectral Image Fusion via Logarithmic Low-rank Tensor Ring
Decomposition

Jun Zhang1,2, Lipeng Zhu3, Chao Wang4, Shutao Li1,∗

††footnotetext:  ∗Corresponding author.
1College of Electrical and Information Engineering, Hunan University, Changsha
410082, Hunan, China
2College of Science, Nanchang Institute of Technology, Nanchang 330099,
Jiangxi, China
3Jiangxi Province Key Laboratory of Water Information Cooperative Sensing and
Intelligent Processing, Nanchang Institute of Technology, Nanchang 330099,
Jiangxi, China
4Department of Statistics and Data Science, Southern University of Science and
Technology, Shenzhen 518055, Guangdong, China
Email addresses<EMAIL_ADDRESS>(J. Zhang<EMAIL_ADDRESS>(L.
Zhu<EMAIL_ADDRESS>(C. Wang<EMAIL_ADDRESS>(S. Li).

Abstract. Integrating a low-spatial-resolution hyperspectral image (LR-HSI)
with a high-spatial-resolution multispectral image (HR-MSI) is recognized as a
valid method for acquiring HR-HSI. Among the current fusion approaches, the
tensor ring (TR) decomposition-based method has received growing attention
owing to its superior performance on preserving the spatial-spectral
correlation. Furthermore, the low-rank property in some TR factors has been
exploited via the matrix nuclear norm regularization along mode-2. On the
other hand, the tensor nuclear norm (TNN)-based approaches have recently
demonstrated to be more efficient on keeping high-dimensional low-rank
structures in tensor recovery. Here, we study the low-rankness of TR factors
from the TNN perspective and consider the mode-2 logarithmic TNN (LTNN) on
each TR factor. A novel fusion model is proposed by incorporating this LTNN
regularization and the weighted total variation which is to promote the
continuity of HR-HSI in the spatial-spectral domain. Meanwhile, we have
devised a highly efficient proximal alternating minimization algorithm to
solve the proposed model. The experimental results indicate that our method
improves the visual quality and exceeds the existing state-of-the-art fusion
approaches with respect to various quantitative metrics.

Keywords: Hyperspectral image fusion, tensor ring decomposition, logarithmic
function, tensor nuclear norm, proximal alternating minimization.

##  1 Introduction

Hyperspectral imaging technology emerges from the fusion of imaging and
spectral technologies. In contrast to conventional imaging methods,
hyperspectral imaging captures numerous minuscule spectral regions across the
electromagnetic spectrum, offering a wealth of spectral information [1, 2].
Consequently, the exploitation of hyperspectral data proves particularly
advantageous for precise feature classification [3, 4], denoising [5],
unmixing [6], and more. Unfortunately, due to limitations inherent in the
imaging techniques employed, hyperspectral images often suffer from low
spatial resolution, constraining their practical utility [7, 8]. On the other
hand, multispectral imaging sensors yield higher spatial resolution images
[9]. Hence, a viable strategy for achieving hyperspectral super-resolution
involves fusing a low-resolution HSI (LR-HSI) with a high-resolution MSI (HR-
MSI) of the same scene.

In recent years, significant efforts have been devoted to exploring HSI-MSI
fusion techniques. These methods can be systematically classified into three
overarching categories: deep learning-based [11, 12, 10], matrix
decomposition-based [14, 15, 13, 16, 17], and tensor decomposition-based
methods [18, 19, 20, 21]. Matrix decomposition-based methods operate under the
assumption that each spectral vector within HR-HSI can be formulated as a
linear combination of distinct spectral features. Consequently, estimating the
desired HR-HSI involves the computation of both the spectral basis and its
corresponding coefficient matrix. Particularly, by imposing prior constraints
on the coefficient matrix, some effective fusion approaches have been proposed
[14, 15, 16, 17]. In [14], the spatial smoothness of the fused image was
improved by applying the total variation (TV) constraint to the coefficient
matrix. In pursuit of promoting local low-rank characteristics within images,
Dian et al. [15] introduced an HSI super-resolution technique utilizing
nuclear norm regularization. On this foundation, to better retain crucial
information in the image, Zhang et al. [17] innovatively integrated a weighted
vector with the nuclear norm regularization. Notably, the TV regularization
was incorporated to effectively preserve the piecewise smoothness in the
spatial domain. Although the vectorization of individual bands in these
methods does not alter the values within the original 3D data, it disrupts the
inherent structure of HSI, resulting in a deterioration of fusion performance.

Tensor representation can effectively preserve the inherent structural
information of HSI. Nowadays, tensor factorization, which includes Canonical
polyadic (CP) decomposition [18, 22], Tucker decomposition [19, 23, 24],
tensor-train (TT) decomposition [20, 25], and tensor ring (TR) decomposition
[26, 27], has become a popular method to deal with the HSI-MSI fusion problem.
As presented in [18], the CP decomposition was employed to partition the HR-
HSI into three parts, followed by iterative updates for each part.
Nonetheless, using the CP decomposition presents challenges in locating the
optimal latent tensors. To overcome this problem, based on the Tucker
decomposition, Li et al. [19] presented a coupled sparse tensor factorization
approach for HSI-MSI fusion. Specifically, they postulated that the HR-HSI
could be formulated as a core tensor multiplied by the dictionaries
corresponding to the three modes. On this basis, Xu et al. [23] introduced a
unidirectional TV-based approach. Especially, they applied the $\ell_{1}$-norm
to the core tensor so as to capture sparsity, while employing unidirectional
TV on three dictionaries to characterize piecewise smoothness. In [24], Wang
et al. proposed a novel HSI-MSI fusion approach using spectral unmixing and
nonlocal low-rank tensor decomposition. They harnessed nonlocal Tucker
decomposition to capture both spatial-spectral correlations and spatial self-
similarity within HSI, augmenting it with spectral regularization technique to
mitigate spectral distortions. Although the aforementioned methods based on
Tucker decomposition excel in preserving the spectral-spatial correlation of
HSI, they face the challenge of dimensionality catastrophe.

Fortunately, compared to CP decomposition and Tucker decomposition, tensor
network-based approaches such as TT decomposition [20, 25] and TR
decomposition [26, 27, 28, 29, 31, 32, 21, 30] have a stronger ability to mine
the intrinsic structure of data. In [20], Dian et al. utilized non-local
similarity to group similar patches within an image, leading to the creation
of multiple 4D tensors. They further applied a low TT rank constraint to these
4D tensors, and thus acquired better fusion performance. However, the TT
decomposition requires the imposition of rank-1 constraint on the border
factors, resulting in limitations in its expressiveness and flexibility. To
address the problem caused by the TT decomposition, the TR decomposition model
was initially proposed for color image and video recovery [26]. TR can be seen
as the linear combination of a set of TT representations [27]. It inherits the
main properties of TT. Furthermore, owing to the rotation invariance property
of latent factors, TR exhibits greater flexibility in its model constraints
compared to TT. In [28], a TV-regularized TR completion model was proposed for
HSI reconstruction. In this work, the TR decomposition was employed to explore
low-rank tensor information, and TV was utilized to promote spatial
smoothness. In addition, TR decomposition was successfully applied to HSI
denoising [33] and neural network compression [34]. Recently, TR decomposition
has also been developed for HSI-MSI fusion [31, 32, 21, 30]. He et al. [21]
proposed a coupled TR decomposition approach, which simultaneously learns the
TR factors of the HR-HSI among a pair of HSI and MSI. This approach is very
different from traditional methods since the aim of fusion is transformed into
estimating the TR factors. Based on this TR framework, Chen et al. [30]
further imposed weighted TV regularization to constrain the TR factors so as
to capture its spatial-spectral continuity. Although these approaches in the
TR decomposition framework achieve relatively satisfactory fusion results,
they don’t take into account the low-rank characteristics of HR-HSI. In other
words, the low-rank property of TR factors is ignored in the TR decomposition
framework.

Owing to the low-rank nature of HSIs, the nuclear norm regularization has been
widely used in HSI denoising and fusion [35, 36, 15, 17]. In [29], Yuan et al.
applied the matrix nuclear norm to the rank-mode unfoldings of the TR factors
and proposed a TR low-rank factors model for tensor completion. In [21], a
coupled TR factorization model with nuclear norm regularization was proposed
for HSI-MSI fusion. Specifically, to characterize the low-rank property of HSI
along the spectral dimension, the authors utilized matrix nuclear norm
regularization to the third core tensor along mode-2. Although these two
approaches acquired acceptable performance, the core tensors were all unfolded
to mode-n matrices, neglecting the correlation between different modes. In
contrast, the tensor nuclear norm (TNN) based on tensor-singular value
decomposition (t-SVD) can effectively maintain the inherent low-rank structure
of the tensor [38, 32, 37]. By introducing TNN into all 3D TR factors, Xu et
al. [32] presented a TR and subspace decompositions-based model for HSI-MSI
fusion to preserve the internal tensor structure naturally. However, TNN
treats each singular value equally, ignoring the physical significance of the
different singular values and generating a biased approximation of tensor low-
rank. To address this problem, the use of non-convex functions instead of the
$\ell_{1}$-norm to approximate the rank function can more accurately extract
the low-rank components from the degenerate HSI [39, 40]. In this way, it
allows larger singular values to be penalized less, effectively preserving the
main information of the image data. In approximating the nuclear norm, the
log-sum of singular values has demonstrated superior performance when compared
to the sum of singular values. Thus, in [40], Dian et al. introduced an HSI-
MSI fusion technique that employs subspace-based logarithmic TNN (LTNN)
regularization.

Inspired by the success of LTNN in subspace-based HSI-MSI fusion, in this
study, we expand its applicability to fusion based on the TR decomposition
framework. Specifically, we develop a novel model for HSI-MSI fusion by
combining LTNN with weighted TV regularization. To exploit the intrinsic low-
rank structure of HSIs, we directly apply the LTNN regularization to each of
the TR factors in the framework of TR decomposition. In addition, we introduce
a weighted TV regularization to preserve the continuity of the fused image in
the spatial-spectral domain. We devise a proximal alternating minimization
algorithm to address the proposed model. Our proposed approach exceeds the
current state-of-the-art fusion approach by increasing the PSNR value by over
0.69 and reducing the ERGAS value by 6.9$\%$.

The main contributions can be summarized as follows:

  * •

We propose a novel fusion model under TR decomposition and utilize TNN for TR
factors. To the best of our knowledge, this is the first work in HSI-MSI
fusion that investigates the low-rankness of TR factors from the perspective
of TNN.

  * •

We address the skew effect of TNN and introduce the LTNN method to capture the
low-rankness among these TR factors. Additionally, we incorporate the weighted
TV to ensure the continuity of the fused image.

The remainder of the paper is structured as follows. In Section 2, we
introduce the related notations and the fusion framework based on the TR
decomposition, and two previous fusion models. In Section 3, we propose a
novel HSI-MSI fusion model by combining logarithmic tensor nuclear norm with
the weighted TV and present its optimization algorithm. Experimental results
and convergence analysis are shown in Section 4. Finally, some conclusions and
outlook of this paper are given in Section 5.

##  2 Preliminaries

In this section, we will first introduce concepts and notations related to
tensor, and the definition of tensor ring (TR) decomposition. Subsequently,
the problem formulation concerning the fusion of HSI and MSI is outlined.
Finally, we will present existing fusion models related to TR decomposition.

###  2.1 Notations

In this paper, we primarily adopt tensor notations consistent with those found
in existing literature [40, 30]. Lowercase and uppercase letters are used for
denoting scalars, i.e., $n$, $N\in\mathbb{R}$. Bold lowercase is adopted to
indicate vectors, i.e., $\mathbf{x}\in\mathbb{R}{{}^{N}}$. Correspondingly,
bold uppercase letters denote matrices, i.e.,
${\mathbf{X}}\in\mathbb{R}{{}^{M\times N}}$. The calligraphic letters are
employed to represent tensors of order $N\ (N\geq 3)$, i.e.,
$\mathcal{X}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times\cdots\times{I_{N}}}}$,
where ${I_{i}}$ is the dimension of the $i$th mode. Especially, scalars,
vectors, and matrices can be seen as the zero-order, first-order, and second-
order tensors, respectively. $\mathcal{X}({i_{1}},{i_{2}},\cdots,{i_{N}})$ or
${x_{{i_{1}},{i_{2}},\cdots,{i_{N}}}}$ stands for the elemental value of
$\mathcal{X}$ at position $({i_{1}},{i_{2}},\cdots,{i_{N}})$.
${\mathbf{X}_{(n)}}\in\mathbb{R}{{}^{{I_{n}}\times{I_{1}}{I_{2}}\cdots{I_{n-1}}{I_{n+1}}\cdots{I_{N}}}}$
and
${\mathbf{X}}_{<n>}\in\mathbb{R}^{{I_{n}}\times{I_{n+1}}\cdots{I_{N}}{I_{1}}\cdots{I_{n-1}}}$
represent the first and second mode-$n$ expansion matrices of tensor
$\mathcal{X}$, respectively. ${\mathrm{Fold}}_{n}(\cdot)$ is the inverse
operator of mode-$n$ matricization.
${\left\|\mathcal{X}\right\|_{F}}=\sqrt{\left\langle{\mathcal{X},\mathcal{X}}\right\rangle}$
denotes the Frobenius-norm of tensor $\mathcal{X}$, while its $\ell_{1}$-norm
is defined by
${\left\|\mathcal{X}\right\|_{1}}=\sum\nolimits_{{i_{1}},{i_{2}},\cdots,{i_{N}}}{\left|{x_{{i_{1}},{i_{2}},\cdots,{i_{N}}}}\right|}$.
The mode-$n$ product of the matrix
${\mathbf{D}}\in\mathbb{R}{{}^{N\times{I_{n}}}}$ and $\mathcal{X}$ is defined
as $\mathcal{Y}=\mathcal{X}\times_{n}{\mathbf{D}}$, where
$\mathcal{Y}\in\mathbb{R}{{}^{{I_{1}}\times\cdots\times{I_{n-1}}\times
N\times{I_{n+1}}\times\cdots\times{I_{N}}}}$.

The objective of TR decomposition is to express a high-dimensional tensor
using circular multilinear products involving a sequence of third-order factor
tensors, commonly referred to as TR factors [26]. The TR decomposition of
$\mathcal{X}$ is expressed as a search for $N$ latent third-order core tensors
$\mathcal{G}=\left\\{\mathcal{G}^{(1)},\mathcal{G}^{(2)},\cdots,\mathcal{G}^{(N)}\right\\}$,
where
$\mathcal{G}^{(n)}\in{\mathbb{R}^{{R_{n}}\times{I_{n}}\times{R_{n+1}}}},\
n=1,2,\cdots,N.$ Here, $R_{1}=R_{N+1}$ and $R=[{R_{1},\cdots,R_{N}}]$ is
called the TR rank of tensor $\mathcal{X}$. This decomposition has a stronger
ability to mine the intrinsic structure of data than other decomposition
approaches. Therefore, they have found extensive utilization in tasks related
to third-order image processing. Here the TR decomposition of the tensor
$\mathcal{X}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times\cdots\times{I_{N}}}}$
can be expressed as

$\mathcal{X}({{i_{1}},{i_{2}},\cdots,{i_{N}}})=\mathrm{Tr}({\mathbf{G}^{(1)}({i_{1}}){\mathbf{G}^{(2)}}({i_{2}})\cdots{\mathbf{G}^{(N)}}({i_{N}})}),$
(2.1)

where $\mathrm{Tr}(\cdot)$ signifies the matrix trace operation, and
$\mathbf{G}^{(n)}({i_{n}})\in{\mathbb{R}^{{R_{n}}\times{R_{n+1}}}}$ means the
$i_{n}$th lateral slicing matrix of the core tensor $\mathcal{G}^{(n)}$. To
simplify the expression, we denote the TR decomposition of a $N$-order tensor
by

$\mathcal{X}=\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\dots,\mathcal{G}^{(N)}).$

The multilinear product of ${\mathcal{G}^{(n)}}$ and ${\mathcal{G}^{(n+1)}}$
is represented by
${\mathcal{G}^{(n,n+1)}}\in{\mathbb{R}^{{R_{n}}\times{I_{n}}{I_{n+1}}\times{R_{n+2}}}}$.
Specifically,

$\mathcal{G}^{(n,n+1)}(({j_{n}}-1){I_{n}}+{i_{n}})={\mathcal{G}^{(n)}}({i_{n}})\times{\mathcal{G}^{(n+1)}}({j_{n}}),$
(2.2)

for ${i_{n}}=1,\cdots,{I_{n}},\ {j_{n}}=1,\cdots,{I_{n+1}}$. In addition,
$\overleftarrow{\mathcal{X}}_{n}\in\mathbb{R}{{}^{{I_{n}}\times{I_{n+1}}\times\cdots\times{I_{N}}\times{I_{1}}\times\cdots\times{I_{n-1}}}}$
is a tensor derived from $\mathcal{X}$ by cyclically shifting its dimensions
by $n$ steps. Thus, the following equation can be obtained:

$\overleftarrow{\mathcal{X}}_{n}=\mathbf{\Phi}(\mathcal{G}^{(n)},\mathcal{G}^{(n+1)},\cdots,\mathcal{G}^{(N)},\mathcal{G}^{(1)},\cdots,\mathcal{G}^{(n-1)}).$
(2.3)

Based on this property, the second mode-n expansion matrix of the tensor
$\mathcal{X}$ can be written as

$\mathbf{X}_{<n>}=\mathbf{G}_{(2)}^{(n)}{\left({\mathbf{G}_{<2>}^{(\neq
n)}}\right)^{T}},$ (2.4)

where $\mathcal{G}^{(\neq
n)}\in{\mathbb{R}^{{R_{n+1}}\times\prod^{N}_{i=1,i\neq
n}{I_{i}}\times{R_{n}}}}$ is a subchain tensor acquired by fusing all $N-1$
tensors excluding the $n$th core tensor.

Next, the other crucial property is introduced. Namely, given that the TR
decomposition of $\mathcal{X}$ is
$\mathcal{X}=\mathbf{\Phi}\left({\mathcal{G}^{(1)}},{\mathcal{G}^{(2)}},\cdots,{\mathcal{G}^{(N)}}\right)$,
the operation of mode-$n$ multiplication on $\mathcal{X}$ is expressed as
follows:

$\mathcal{X}\times_{n}\mathbf{D}=\mathbf{\Phi}({\mathcal{G}^{(1)}},{\mathcal{G}^{(2)}},\cdots,{\mathcal{G}^{(n)}}{\times_{2}}{\mathbf{D}},\cdots,{\mathcal{G}^{(N)}}).$
(2.5)

###  2.2 Problem formulation

In actual scenarios, HSI and MSI constitute third-dimensional tensor data,
where the initial two modes represent the spatial dimensions of the image,
while the third mode corresponds to the spectral dimension. Within the context
of HSI fusion, a pair of images is provided, namely an LR-HSI denoted as
$\mathcal{Y}\in\mathbb{R}{{}^{w\times h\times B}}$ and an HR-MSI termed
$\mathcal{Z}\in\mathbb{R}{{}^{W\times H\times b}}$, both captured within the
same scene. HSI-MSI fusion focuses on restoring an HR-HSI
$\mathcal{X}\in\mathbb{R}{{}^{W\times H\times B}}$ from $\mathcal{Y}$ and
$\mathcal{Z}$. Here, $W>w$ and $H>h$ are the spatial sizes and $B>b$ stands
for the number of spectral bands. Specifically, the procedures for obtaining
LR-HSI and HR-MSI can be viewed as the downsampling spatially and spectrally,
respectively.

The separability applies to both the point spread function (PSF) of the
hyperspectral sensor and the downsampling matrices for the width and height
modes. As a result, we have

$\mathcal{Y}=\mathcal{X}{\times_{1}}{\mathbf{U}_{1}}{\times_{2}}{\mathbf{U}_{2}},\
\mathcal{Z}=\mathcal{X}{\times_{3}}{\mathbf{U}_{3}},$ (2.6)

where ${{\mathbf{U}}_{1}}\in\mathbb{R}^{w\times W}$ and
${\mathbf{U}_{2}}\in\mathbb{R}^{h\times H}$ are the separable degradation
operators along the width and height modes, respectively. Separable operators
are shown to offer more advantages in computing and optimization [41]. And
${\mathbf{U}_{3}}\in\mathbb{R}^{b\times B}$ represents the spectral
downsampling matrix.

By utilizing the simplified expression of the TR decomposition, it becomes
possible to depict the HR-HSI as

$\mathcal{X}=\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)}),$
(2.7)

where $\mathcal{G}^{(1)}\in{\mathbb{R}^{{R_{1}}\times W\times{R_{2}}}}$,
$\mathcal{G}^{(2)}\in{\mathbb{R}^{{R_{2}}\times H\times{R_{3}}}}$ and
$\mathcal{G}^{(3)}\in{\mathbb{R}^{{R_{3}}\times B\times{R_{1}}}}$ denote TR
factors associated with two spatial dimensions and one spectral dimension,
respectively. $R=[R_{1},R_{2},R_{3}]$ signifies the TR rank of the HR-HSI.
According to equations (2.6) and (2.7), it is easy to get the following
expressions for the LR-HSI and HR-MSI:

$\mathcal{Y}=\mathbf{\Phi}(\mathcal{G}^{(1)}\times_{2}{\mathbf{U}_{1}},\mathcal{G}^{(2)}\times_{2}{\mathbf{U}_{2}},\mathcal{G}^{(3)})$
(2.8)

and

$\mathcal{Z}=\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)}\times_{2}{\mathbf{U}_{3}}).$
(2.9)

In this way, the fusion problem is transformed into the estimation of three TR
factors by making use of the LR-HSI $\mathcal{Y}$ and the HR-MSI
$\mathcal{Z}$.

###  2.3 Related models

As previously mentioned, the aim of HSI-MSI fusion is to calculate the three
TR factors via (2.8) and (2.9). With the downsampling matrices
$\mathbf{U}_{1}$, ${{\mathbf{U}}_{2}}$ and $\mathbf{U}_{3}$ known, the three
TR factors can be solved by [21]

$\min\limits_{\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)}}\frac{1}{2}\left\|{\mathcal{Y}-\mathbf{\Phi}(\mathcal{G}^{(1)}\times_{2}{\mathbf{U}_{1}},{\mathcal{G}^{(2)}}\times_{2}{\mathbf{U}_{2}},\mathcal{G}^{(3)})}\right\|_{F}^{2}+\frac{\lambda}{2}\left\|{\mathcal{Z}-\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)}\times_{2}\mathbf{U}_{3})}\right\|_{F}^{2},$
(2.10)

where the constant ${\lambda>0}$ is employed to equilibrium the two terms
concerning data fidelity. Compared with the traditional matrix/tensor
decomposition-based approaches, TR decomposition generates a more accurate
representation for high-dimensional tensor data such as HR-HSI. Nevertheless,
the traditional TR decomposition model (2.10) only captures the spatial-
spectral high-correlation of HR-HSI. Furthermore, deriving the TR factors
directly from the model (2.10) poses an issue of instability. Therefore,
additional prior information should be integrated to acquire a stable solution
[30].

In order to design an effective regularization of the tensor factors, it is
necessary to explore the intrinsic relationship between tensor factors and HR-
HSI [42]. In terms of matrix theory, one can derive from (2.4) that each
column of $\mathbf{X}_{<n>}$ can be formulated as a linear combination of all
columns of the factor $\mathbf{G}_{(2)}^{(n)}$. Therefore, drawing from the
principle that continuous bases have the capacity to depict continuous data,
all columns of three factors $\\{\mathbf{G}_{(2)}^{(n)}\\}_{n=1}^{3}$ should
be continuous. For this purpose, the work presented in [30] introduced an
approach within the TR decomposition framework. The weighted total variation
(TV) with updating the weights was presented to constrain the TR factors so as
to promote their continuity. Specifically, the resulting model can be
expressed as

$\begin{split}\min\limits_{\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)}}&\frac{1}{2}\left\|{\mathcal{Y}-\mathbf{\Phi}(\mathcal{G}^{(1)}\times_{2}\mathbf{U}_{1},\mathcal{G}^{(2)}\times_{2}\mathbf{U}_{2},\mathcal{G}^{(3)})}\right\|_{F}^{2}+\frac{\lambda}{2}\left\|\mathcal{Z}-\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)}\times_{2}{\mathbf{U}_{3}})\right\|_{F}^{2}\\\
&+\alpha\sum_{n=1}^{3}\left\|\mathcal{W}^{(n)}\odot(\mathcal{G}^{(n)}\times_{2}\mathbf{D})\right\|_{1},\end{split}$
(2.11)

where $\mathcal{W}^{(n)}\in{\mathbb{R}^{{R_{n}}\times{I_{n}}\times{R_{n+1}}}}$
stands for a non-negative weighted tensor, and
$\mathbf{D}\in\mathbb{R}^{{I_{n}}\times{I_{n}}}$ is a first-order difference
square matrix. Here, ${I_{n}}$ refers to the second dimension of
$\mathcal{G}^{(n)}$.

##  3 The proposed method

In this section, we first propose a novel HSI fusion model by combining
logarithmic nuclear norm with the weighted TV. Then an efficient numerical
algorithm is designed.

###  3.1 Proposed model

To delve deeper into the latent information held within the TR factors
$\mathcal{G}^{(n)}$ ($n=1,2,3$), we introduce a low-rank constraint to each of
these factors. In [29], Yuan et al. applied the matrix nuclear norm to the
rank-mode unfoldings of all TR factors for tensor completion. In [21], the
nuclear norm regularization along model-2 of only the third TR factor was
proposed for HSI-MSI fusion. Due to the three-dimensional structure of each
$\mathcal{G}^{(n)}$, utilizing a matrix nuclear norm on its unfolded
representation could potentially disrupt its inherent high-order structural
correlation, consequently impacting the ultimate reconstruction outcomes. To
address this, we instead apply a tensor nuclear norm (TNN) constraint to every
individual TR factor within the framework of TR decomposition.

Nowadays, the TNN-based HSI-MSI fusion has attracted increasing attention
owing to its advantages in preserving high-dimensional low-rank structures
[38, 32, 37]. Nevertheless, the general TNN method merely represents a skewed
approximation of tensor low-rank and does not effectively enhance the
solution’s low-rank characteristics. Therefore, we design a logarithmic tensor
nuclear norm (LTNN) and apply it to each TR factor $\mathcal{G}^{(n)}$. The
proposed LTNN provides a more accurate low-rank approximation, since the
logarithmic function produces an approximation that is more closely aligned
with the $\ell_{0}$-norm compared to the $\ell_{1}$-norm. Particularly, we
employ the model-2 tensor singular value thresholding (t-SVT) to handle
singular values differently. To clarify, we apply the milder reduction to the
significant singular values, retaining the primary data components, while
implementing a more substantial reduction to the smaller ones to mitigate
random errors. Specifically, our proposed HSI-MSI fusion model can be
generated as

$\begin{split}\min\limits_{\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)}}&\frac{1}{2}\left\|{\mathcal{Y}-\mathbf{\Phi}(\mathcal{G}^{(1)}\times_{2}\mathbf{U}_{1},\mathcal{G}^{(2)}\times_{2}\mathbf{U}_{2},\mathcal{G}^{(3)})}\right\|_{F}^{2}+\frac{\lambda}{2}\left\|{\mathcal{Z}-\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)}\times_{2}\mathbf{U}_{3})}\right\|_{F}^{2}\\\
&+\sum_{n=1}^{3}\left(\alpha\left\|\mathcal{W}^{(n)}\odot(\mathcal{G}^{(n)}\times_{2}\mathbf{D})\right\|_{1}+\beta\left\|\mathcal{G}^{(n)}\right\|_{\mathrm{LTNN}}\right),\end{split}$
(3.1)

where $\beta$ is a regularization parameter and

$\left\|\mathcal{G}\right\|_{\mathrm{LTNN}}=\frac{1}{S}\sum\limits_{i=1}^{S}{\sum\limits_{j}{\log({\sigma_{j}}({\overline{\mathbf{G}}}^{[i]})+\varepsilon)}}.$
(3.2)

Here, $S$ represents the size of the second mode of tensor $\mathcal{G}$,
$\varepsilon>0$ is an extremely small positive value, and
${\sigma_{j}}({\overline{\mathbf{G}}}^{[i]})$ signifies the $j$th singular
value of the $i$th lateral slice of tensor
$\overline{\mathcal{G}}=\mathrm{fft}\left(\mathcal{G},[\enspace],2\right)$.
Note that our LTNN is different from the one in [40] since the corresponding
Fourier transform is along the second dimension instead of the third one.

###  3.2 Algorithm development

In this subsection, we introduce an iterative scheme for solving the proposed
model (3.1). Even though this model exhibits joint non-convexity with respect
to the three factors $\mathcal{G}^{(n)}$ (where $n=1,2,3$), it demonstrates
convexity for each individual separable variable. Thus, we utilize the
proximal alternating minimization (PAM) framework [43], ensuring the
convergence of the solution to a stationary point of the objective function.

For the sake of convenience, we consider the objective function of (3.1) as
$f(\mathcal{G}^{(1)},\mathcal{G}^{(2)},\mathcal{G}^{(3)})$. Additionally, the
proximal term is introduced to update each factor within the PAM framework. In
this way, we can solve the proposed model through alternately iterating three
subproblems as below:

$\left\\{\begin{aligned}
&\mathcal{G}^{(1),k+1}=\arg\min\limits_{\mathcal{G}^{(1)}}f(\mathcal{G}^{(1)},\mathcal{G}^{(2),k},\mathcal{G}^{(3),k})+\frac{\eta}{2}\left\|\mathcal{G}^{(1)}-\mathcal{G}^{(1),k}\right\|_{F}^{2},\\\
&\mathcal{G}^{(2),k+1}=\arg\min\limits_{\mathcal{G}^{(2)}}f(\mathcal{G}^{(1),k+1},\mathcal{G}^{(2)},\mathcal{G}^{(3),k})+\frac{\eta}{2}\left\|\mathcal{G}^{(2)}-\mathcal{G}^{(2),k}\right\|_{F}^{2},\\\
&\mathcal{G}^{(3),k+1}=\arg\min\limits_{\mathcal{G}^{(3)}}f(\mathcal{G}^{(1),k+1},\mathcal{G}^{(2),k+1},\mathcal{G}^{(3)})+\frac{\eta}{2}\left\|\mathcal{G}^{(3)}-\mathcal{G}^{(3),k}\right\|_{F}^{2},\end{aligned}\right.$
(3.3)

where $\eta>0$ is a parameter. And
$\frac{\eta}{2}\left\|{\cdot}\right\|_{F}^{2}$ signifies the proximal term. In
the following, we describe the solution of each subproblem in detail.

First, we discuss the solution of the $\mathcal{G}^{(1)}$-subproblem. This
subproblem can be expressed as

$\begin{split}\min\limits_{\mathcal{G}^{(1)}}&\frac{1}{2}\left\|{\mathcal{Y}-\mathbf{\Phi}(\mathcal{G}^{(1)}\times_{2}{\mathbf{U}_{1}},\mathcal{G}^{(2),k}\times_{2}\mathbf{U}_{2},\mathcal{G}^{(3),k})}\right\|_{F}^{2}+\frac{\lambda}{2}\left\|{\mathcal{Z}-\mathbf{\Phi}({\mathcal{G}^{(1)},\mathcal{G}^{(2),k},\mathcal{G}^{(3),k}\times_{2}{\mathbf{U}}_{3}})}\right\|_{F}^{2}\\\
&+\alpha{\left\|{\mathcal{W}^{(1)}\odot(\mathcal{G}^{(1)}\times_{2}{\mathbf{D}})}\right\|_{1}}+\beta{\left\|\mathcal{G}^{(1)}\right\|_{\mathrm{LTNN}}}+\frac{\eta}{2}\left\|{\mathcal{G}^{(1)}-\mathcal{G}^{(1),k}}\right\|_{F}^{2}.\end{split}$
(3.4)

To efficiently solve (3.4), we employ the alternating direction method of
multipliers (ADMM) [44, 45]. The ADMM is a fast approach for handling the
minimization of non-differentiable and high-order functionals. To begin, the
unconstrained problem (3.4) can be transformed into the following constrained
optimization problem by introducing two auxiliary variables $\mathcal{R}_{1}$
and $\mathcal{V}_{1}$:

$\begin{split}\min\limits_{\mathcal{G}^{(1)},\mathcal{R}_{1},\mathcal{V}_{1}}&\frac{1}{2}\left\|\mathcal{Y}-\mathbf{\Phi}(\mathcal{G}^{(1)}\times_{2}\mathbf{U}_{1},\mathcal{G}^{(2),k}\times_{2}\mathbf{U}_{2},\mathcal{G}^{(3),k})\right\|_{F}^{2}+\frac{\lambda}{2}\left\|\mathcal{Z}-\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2),k},\mathcal{G}^{(3),k}\times_{2}\mathbf{U}_{3})\right\|_{F}^{2}\\\
&+\alpha\left\|\mathcal{W}^{(1)}\odot\mathcal{R}_{1}\right\|_{1}+\beta\left\|\mathcal{V}_{1}\right\|_{\mathrm{LTNN}}+\frac{\eta}{2}\left\|\mathcal{G}^{(1)}-\mathcal{G}^{(1),k}\right\|_{F}^{2},\\\
\hbox{s.t.}\quad&\mathcal{R}_{1}=\mathcal{G}^{(1)}\times_{2}\mathbf{D},\
\mathcal{V}_{1}=\mathcal{G}^{(1)}.\end{split}$ (3.5)

Here we can obtain its augmented Lagrangian function:

$\begin{split}L_{\mu}&(\mathcal{G}^{(1)},\mathcal{R}_{1},\mathcal{V}_{1};\mathcal{M}_{1},\mathcal{N}_{1})\\\
&=\frac{1}{2}\left\|\mathcal{Y}-\mathbf{\Phi}(\mathcal{G}^{(1)}\times_{2}\mathbf{U}_{1},\mathcal{G}^{(2),k}\times_{2}\mathbf{U}_{2},\mathcal{G}^{(3),k})\right\|_{F}^{2}+\frac{\lambda}{2}\left\|\mathcal{Z}-\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2),k},\mathcal{G}^{(3),k}\times_{2}\mathbf{U}_{3})\right\|_{F}^{2}\\\
&+\frac{\eta}{2}\left\|\mathcal{G}^{(1)}-\mathcal{G}^{(1),k}\right\|_{F}^{2}+\alpha\left\|\mathcal{W}^{(1)}\odot\mathcal{R}_{1}\right\|_{1}+\frac{\mu}{2}\left\|\mathcal{R}_{1}-\mathcal{G}^{(1)}\times_{2}\mathbf{D}+\frac{\mathcal{M}_{1}}{\mu}\right\|_{F}^{2}\\\
&+\beta\left\|\mathcal{V}_{1}\right\|_{\mathrm{LTNN}}+\frac{\mu}{2}\left\|\mathcal{V}_{1}-\mathcal{G}^{(1)}+\frac{\mathcal{N}_{1}}{\mu}\right\|_{F}^{2},\end{split}$
(3.6)

where $\mathcal{M}_{1}$ , $\mathcal{N}_{1}$ are the Lagrangian multipliers and
$\mu>0$ is a constant. Then the update of $\mathcal{G}^{(1)}$ can be acquired
by iteratively calculating the following three subproblems and updating the
two Lagrangian multipliers:

a) The $\mathcal{R}_{1}$-subproblem is written as

$\min\limits_{\mathcal{R}_{1}}\alpha\left\|\mathcal{W}^{(1)}\odot\mathcal{R}_{1}\right\|_{1}+\frac{\mu}{2}\left\|\mathcal{R}_{1}-\mathcal{G}^{(1)}\times_{2}\mathbf{D}+\frac{\mathcal{M}_{1}}{\mu}\right\|_{F}^{2}.$
(3.7)

This is a weighted ${\ell_{1}}$-norm minimization, and its closed-form
solution is achieved by the soft-shrinkage operator, that is

$\mathcal{R}_{1}=\mathrm{soft}(\mathcal{J}_{1},\frac{\alpha}{\mu}\mathcal{W}^{(1)}).$
(3.8)

Here,
$\mathcal{J}_{1}=\mathcal{G}^{(1)}\times_{2}\mathbf{D}-\left(\mathcal{M}_{1}/\mu\right)$,
soft($a$, $b$)=sign($a$) $\ast$ $\max$($|a|$$-$$b$, 0). Following the
reweighting approach [46], the weight coefficient is assigned as

$\mathcal{W}^{(1)}=1/\left(\left|\mathcal{J}_{1}\right|+\varsigma\right),$
(3.9)

where $\varsigma>0$ is a small constant used to eliminate singularity.

b) The $\mathcal{G}^{(1)}$-subproblem is formulated as

$\displaystyle\min_{\mathcal{G}^{(1)}}$
$\displaystyle\frac{1}{2}\left\|\mathcal{Y}-\mathbf{\Phi}(\mathcal{G}^{(1)}\times_{2}\mathbf{U}_{1},\mathcal{G}^{(2),k}\times_{2}\mathbf{U}_{2},\mathcal{G}^{(3),k})\right\|_{F}^{2}+\frac{\lambda}{2}\left\|\mathcal{Z}-\mathbf{\Phi}(\mathcal{G}^{(1)},\mathcal{G}^{(2),k},\mathcal{G}^{(3),k}\times_{2}\mathbf{U}_{3})\right\|_{F}^{2}$
(3.10) $\displaystyle\quad+$
$\displaystyle\frac{\eta}{2}\left\|\mathcal{G}^{(1)}-\mathcal{G}^{(1),k}\right\|_{F}^{2}+\frac{\mu}{2}\left\|\mathcal{R}_{1}-\mathcal{G}^{(1)}\times_{2}\mathbf{D}+\frac{\mathcal{M}_{1}}{\mu}\right\|_{F}^{2}+\frac{\mu}{2}\left\|\mathcal{V}_{1}-\mathcal{G}^{(1)}+\frac{\mathcal{N}_{1}}{\mu}\right\|_{F}^{2}.$

For convenience, let
$\mathbf{P}_{1}=\left(\left(\mathcal{G}^{(2),k}\times_{2}\mathbf{U}_{2}\right)\mathcal{G}^{(3),k}\right)_{<2>}^{T}$
and
$\mathbf{P}_{2}=\left(\mathcal{G}^{(2),k}\left(\mathcal{G}^{(3),k}\times_{2}\mathbf{U}_{3}\right)\right)_{<2>}^{T}$.
Utilizing the matrix representation of TR decomposition (2.4), the above
problem (3.10) can be rewritten as

$\displaystyle\min_{\mathbf{G}_{(2)}^{(1)}}$
$\displaystyle\frac{1}{2}\left\|\mathbf{Y}_{<1>}-\mathbf{U}_{1}\mathbf{G}_{(2)}^{(1)}\mathbf{P}_{1}\right\|_{F}^{2}+\frac{\lambda}{2}\left\|\mathbf{Z}_{<1>}-\mathbf{G}_{(2)}^{(1)}\mathbf{P}_{2}\right\|_{F}^{2}+\frac{\eta}{2}\left\|\mathbf{G}_{(2)}^{(1)}-\mathbf{G}_{(2)}^{(1),k}\right\|_{F}^{2}$
(3.11)
$\displaystyle+\frac{\mu}{2}\left\|\mathbf{R}_{\mathbf{1(2)}}-\mathbf{D}\mathbf{G}_{(2)}^{(1)}+\frac{\mathbf{M}_{1(2)}}{\mu}\right\|_{F}^{2}+\frac{\mu}{2}\left\|\mathbf{V}_{1(2)}-\mathbf{G}_{(2)}^{(1)}+\frac{\mathbf{N}_{1(2)}}{\mu}\right\|_{F}^{2}.$

Its Euler-Lagrange equation is

$\displaystyle\mathbf{U}_{1}^{T}\mathbf{U}_{1}\mathbf{G}_{(2)}^{(1)}\mathbf{P}_{1}\mathbf{P}_{1}^{T}+\lambda\mathbf{G}_{(2)}^{(1)}\mathbf{P}_{2}\mathbf{P}_{2}^{T}+\eta\mathbf{G}_{(2)}^{(1)}+\mu\mathbf{D}^{T}\mathbf{D}\mathbf{G}_{(2)}^{(1)}+\mu\mathbf{G}_{(2)}^{(1)}=\mathbf{U}_{1}^{T}\mathbf{Y}_{<1>}\mathbf{P}_{1}^{T}$
(3.12)
$\displaystyle+\lambda\mathbf{Z}_{<1>}\mathbf{P}_{2}^{T}+\eta\mathbf{G}_{(2)}^{(1),k}+\mu\mathbf{D}^{T}\left(\mathbf{R}_{1(2)}+\frac{\mathbf{M}_{1(2)}}{\mu}\right)+\mu\left(\mathbf{V}_{1(2)}+\frac{\mathbf{N}_{1(2)}}{\mu}\right).$

As a general Sylvester equation, this linear system is solved using the
conjugate gradient (CG) method. Subsequently, the TR factor
$\mathcal{G}^{(1)}$ is acquired by folding the solution
${\mathbf{G}}_{(2)}^{(1)}$:
$\mathcal{G}^{(1)}=\operatorname{Fold}_{2}(\mathbf{G}_{(2)}^{(1)})$.

c) The $\mathcal{V}_{1}$-subproblem can be expressed as

$\min\limits_{\mathcal{V}_{1}}\beta\left\|\mathcal{V}_{1}\right\|_{\mathrm{LTNN}}+\frac{\mu}{2}\left\|\mathcal{V}_{1}-\mathcal{G}^{(1)}+\frac{\mathcal{N}_{1}}{\mu}\right\|_{F}^{2}.$
(3.13)

According to the work in [47], its solution can be given by

$\displaystyle\mathcal{V}_{1}=\mathcal{D}_{\mathrm{LTNN}}^{\beta/\mu}\left(\mathcal{G}^{(1)}-\frac{\mathcal{N}_{1}}{\mu}\right)=\mathcal{U}\mathcal{S}_{\mathrm{LTNN}}^{\beta/\mu,\varepsilon}\mathcal{V}^{T},$
(3.14)

where
$\mathcal{G}^{(1)}-\frac{\mathcal{N}_{1}}{\mu}=\mathcal{U}\mathcal{S}\mathcal{V}^{T},\mathcal{S}_{\mathrm{LTNN}}^{\beta/\mu,\varepsilon}=\mathrm{ifft}(\overline{\mathcal{S}}_{\mathrm{LTNN}}^{\beta/\mu,\varepsilon},[\enspace],2)$.
The thresholding operator
$\overline{\mathcal{S}}_{\mathrm{LTNN}}^{\beta/\mu,\varepsilon}$ is defined as

$\overline{\mathcal{S}}_{\mathrm{LTNN}}^{\beta/\mu,\varepsilon}(i_{1},i_{2},i_{3})=\left\\{\begin{array}[]{ll}0,&\mathrm{if}\
c_{2}\leq 0,\\\ \mathrm{sign}(\overline{\mathcal{S}}\
(i_{1},i_{2},i_{3}))\left(\frac{c_{1}+\sqrt{c_{2}}}{2}\right),&\mathrm{if}\
c_{2}>0.\end{array}\right.$ (3.15)

Among them, $\overline{\mathcal{S}}=\mathrm{fft}(\mathcal{S},[\enspace],2)$,
$c_{1}=|\overline{\mathcal{S}}(i_{1},i_{2},i_{3})|-\varepsilon$,
$c_{2}=c_{1}^{2}-4(\beta/\mu-\varepsilon|\overline{\mathcal{S}}(i_{1},i_{2},i_{3})|)$.
In this way, the thresholding operator shrinks less for larger singular values
and more for smaller ones [48].

d) The update of Lagrangian multipliers $\mathcal{M}_{1}$, $\mathcal{N}_{1}$
unfolds in the subsequent manner:

$\left\\{\begin{aligned}
&{\mathcal{M}_{1}}\leftarrow{\mathcal{M}_{1}}+\mu({\mathcal{R}_{1}}-{\mathcal{G}^{(1)}}\times_{2}{\mathbf{D}}),\\\
&{\mathcal{N}_{1}}\leftarrow{\mathcal{N}_{1}}+\mu({\mathcal{V}_{1}}-{\mathcal{G}^{(1)}}).\end{aligned}\right.$
(3.16)

Through iteratively updating $\mathcal{R}_{1}$, $\mathcal{G}^{(1)}$,
$\mathcal{V}_{1}$, $\mathcal{M}_{1}$, and $\mathcal{N}_{1}$, the first factor
$\mathcal{G}^{(1),k+1}$ could be acquired.

Second, we focus on solving the $\mathcal{G}^{(2)}$-subproblem. This
subproblem can be rewritten as

$\begin{split}\min\limits_{\mathcal{G}^{(2)}}&\frac{1}{2}\left\|{\mathcal{Y}-\mathbf{\Phi}({\mathcal{G}^{(1),k+1}\times_{2}{\mathbf{U}_{1}},\mathcal{G}^{(2)}\times_{2}{\mathbf{U}_{2}},\mathcal{G}^{(3),k}})}\right\|_{F}^{2}+\frac{\lambda}{2}\left\|{\mathcal{Z}-\mathbf{\Phi}({\mathcal{G}^{(1),k+1},\mathcal{G}^{(2)},\mathcal{G}^{(3),k}\times_{2}\mathbf{U}_{3}})}\right\|_{F}^{2}\\\
&+\alpha{\left\|\mathcal{W}^{(2)}\odot({\mathcal{G}^{(2)}\times_{2}\mathbf{D}})\right\|_{1}}+\beta{\left\|\mathcal{G}^{(2)}\right\|_{\mathrm{LTNN}}}+\frac{\eta}{2}\left\|\mathcal{G}^{(2)}-\mathcal{G}^{(2),k}\right\|_{F}^{2}.\end{split}$
(3.17)

Similarly, we can obtain the process of optimizing $\mathcal{G}^{(2)}$. By
employing the ADMM, the update for each subproblem is outlined as

$\mathcal{R}_{2}=\mathrm{soft}(\mathcal{J}_{2},\frac{\alpha}{\mu}\mathcal{W}^{(2)}),$
(3.18)
$\displaystyle\mathbf{U}_{2}^{T}\mathbf{U}_{2}\mathbf{G}_{(2)}^{(2)}\mathbf{Q}_{1}\mathbf{Q}_{1}^{T}+\lambda\mathbf{G}_{(2)}^{(2)}\mathbf{Q}_{2}\mathbf{Q}_{2}^{T}+\eta\mathbf{G}_{(2)}^{(2)}+\mu\mathbf{D}^{T}\mathbf{D}\mathbf{G}_{(2)}^{(2)}+\mu\mathbf{G}_{(2)}^{(2)}=\mathbf{U}_{2}^{T}\mathbf{Y}_{<2>}\mathbf{Q}_{1}^{T}$
(3.19)
$\displaystyle+\lambda\mathbf{Z}_{<2>}\mathbf{Q}_{2}^{T}+\eta\mathbf{G}_{(2)}^{(2),k}+\mu\mathbf{D}^{T}\left(\mathbf{R}_{2(2)}+\frac{\mathbf{M}_{2(2)}}{\mu}\right)+\mu\left(\mathbf{V}_{2(2)}+\frac{\mathbf{N}_{2(2)}}{\mu}\right),$
$\mathcal{V}_{2}=\mathcal{D}_{\mathrm{LTNN}}^{\beta/\mu}\left(\mathcal{G}^{(2)}-\frac{\mathcal{N}_{2}}{\mu}\right),$
(3.20) $\left\\{\begin{aligned}
&\mathcal{M}_{2}\leftarrow\mathcal{M}_{2}+\mu(\mathcal{R}_{2}-\mathcal{G}^{(2)}\times_{2}\mathbf{D}),\\\
&\mathcal{N}_{2}\leftarrow{\mathcal{N}_{2}}+\mu(\mathcal{V}_{2}-\mathcal{G}^{(2)}),\end{aligned}\right.$
(3.21)

where $\mathcal{R}_{2}$, $\mathcal{V}_{2}$ and $\mathcal{M}_{2}$,
$\mathcal{N}_{2}$ are auxiliary variables and Lagrangian multipliers,
respectively. In addition, some variables in the aforementioned equations are
defined as

$\left\\{\begin{aligned}
&\mathcal{J}_{2}=\mathcal{G}^{(2)}\times_{2}\mathbf{D}-(\mathcal{M}_{2}/\mu),\
\mathcal{W}^{(2)}=1/(\left|\mathcal{J}_{2}\right|+\varsigma),\\\
&\mathbf{Q}_{1}=(\mathcal{G}^{(3),k}(\mathcal{G}^{(1),k+1}\times_{2}\mathbf{U}_{1}))^{T}_{<2>},\
\mathbf{Q}_{2}=({(\mathcal{G}^{(3),k}\times_{2}\mathbf{U}_{3})\mathcal{G}^{(1),k+1}})^{T}_{<2>},\\\
&\mathcal{G}^{(2)}-\frac{\mathcal{N}_{2}}{\mu}=\mathcal{U}\mathcal{S}\mathcal{V}^{T}.\end{aligned}\right.$
(3.22)

Finally, we consider the $\mathcal{G}^{(3)}$-subproblem, which can be written
as

$\displaystyle\min\limits_{\mathcal{G}^{(3)}}$
$\displaystyle\frac{1}{2}\left\|{\mathcal{Y}-\mathbf{\Phi}(\mathcal{G}^{(1),k+1}\times_{2}{\mathbf{U}_{1}},\mathcal{G}^{(2),k+1}\times_{2}{\mathbf{U}_{2}},\mathcal{G}^{(3)})}\right\|_{F}^{2}+\alpha{\left\|{\mathcal{W}^{(3)}\odot({\mathcal{G}^{(3)}}\times_{2}{\mathbf{D}})}\right\|_{1}}+\beta{\left\|\mathcal{G}^{(3)}\right\|_{\mathrm{LTNN}}}$
(3.23)
$\displaystyle+\frac{\lambda}{2}\left\|{\mathcal{Z}-\mathbf{\Phi}(\mathcal{G}^{(1),k+1},\mathcal{G}^{(2),k+1},\mathcal{G}^{(3)}\times_{2}{\mathbf{U}_{3}})}\right\|_{F}^{2}+\frac{\eta}{2}\left\|{\mathcal{G}^{(3)}-\mathcal{G}^{(3),k}}\right\|_{F}^{2}.$

Similarly, the update of each subproblem is as follows:

$\mathcal{R}_{3}=\mathrm{soft}(\mathcal{J}_{3},\frac{\alpha}{\mu}\mathcal{W}^{(3)}),$
(3.24)
$\displaystyle\lambda\mathbf{U}_{3}^{T}\mathbf{U}_{3}\mathbf{G}_{(2)}^{(3)}\mathbf{T}_{2}\mathbf{T}_{2}^{T}+\mathbf{G}_{(2)}^{(3)}\mathbf{T}_{1}\mathbf{T}_{1}^{T}+\eta\mathbf{G}_{(2)}^{(3)}+\mu\mathbf{D}^{T}\mathbf{D}\mathbf{G}_{(2)}^{(3)}+\mu\mathbf{G}_{(2)}^{(3)}=\mathbf{Y}_{<3>}\mathbf{T}_{1}^{T}$
(3.25)
$\displaystyle+\lambda\mathbf{U}_{3}^{T}\mathbf{Z}_{<3>}\mathbf{T}_{2}^{T}+\eta\mathbf{G}_{(2)}^{(3),k}+\mu\mathbf{D}^{T}\left(\mathbf{R}_{3(2)}+\frac{\mathbf{M}_{3(2)}}{\mu}\right)+\mu\left(\mathbf{V}_{3(2)}+\frac{\mathbf{N}_{3(2)}}{\mu}\right),$
$\mathcal{V}_{3}=\mathcal{D}_{\mathrm{LTNN}}^{\beta/\mu}\left(\mathcal{G}^{(3)}-\frac{\mathcal{N}_{3}}{\mu}\right),$
(3.26) $\left\\{\begin{aligned}
&{\mathcal{M}_{3}}\leftarrow{\mathcal{M}_{3}}+\mu(\mathcal{R}_{3}-\mathcal{G}^{(3)}\times_{2}\mathbf{D}),\\\
&{\mathcal{N}_{3}}\leftarrow{\mathcal{N}_{3}}+\mu(\mathcal{V}_{3}-\mathcal{G}^{(3)}).\end{aligned}\right.$
(3.27)

Here, $\mathcal{R}_{3}$, $\mathcal{V}_{3}$ and $\mathcal{M}_{3}$,
$\mathcal{N}_{3}$ serve as auxiliary variables and Lagrangian multipliers,
respectively. Additionally, some variables in the aforementioned equations are
defined as

$\left\\{\begin{aligned}
&\mathcal{J}_{3}=\mathcal{G}^{(3)}\times_{2}{\mathbf{D}}-\left({\mathcal{M}_{3}}/\mu\right),\
\mathcal{W}^{(3)}=1/\left({\left|{\mathcal{J}_{3}}\right|+\varsigma}\right),\\\
&\mathbf{T}_{1}=({(\mathcal{G}^{(1),k+1}\times_{2}{\mathbf{U}_{1}})(\mathcal{G}^{(2),k+1}\times_{2}{\mathbf{U}_{2}})})_{<2>}^{T},\
\mathbf{T}_{2}=(\mathcal{G}^{(1),k+1}{\mathcal{G}^{(2),k+1}})^{T}_{<2>},\\\
&\mathcal{G}^{(3)}-\frac{\mathcal{N}_{3}}{\mu}=\mathcal{U}\mathcal{S}\mathcal{V}^{T}.\end{aligned}\right.$
(3.28)

We summarize the algorithm for solving the proposed HSI-MSI fusion model (3.1)
in Algorithm 3.1. Due to the non-convex nature of this optimization problem, a
well-chosen initialization will result in a favorable solution. In
experiments, we decompose LR-HSI $\mathcal{Y}$ and HR-MSI $\mathcal{Z}$ using
TR-SVD proposed in the literature [49] to initialize the variables
$\left\\{\mathcal{G}^{(1),0},\mathcal{G}^{(2),0},\mathcal{G}^{(3),0}\right\\}$.

Algorithm 3.1 Optimization of the proposed model (3.1)

0: LR-HSI $\mathcal{Y}$, HR-MSI $\mathcal{Z}$,
$\mathbf{U}_{1},\mathbf{U}_{2},\mathbf{U}_{3}$

1: Parameters: $\mathrm{R}=[\mathrm{R}_{1},\mathrm{R}_{2},\mathrm{R}_{3}]$,
$\lambda,\alpha,\beta,\eta,\mu$, and kMax

2: Initialization:
$\mathcal{G}^{(1),0},\mathcal{G}^{(2),0},\mathcal{G}^{(3),0}$, and $k=0$

3: while $k<$kMax do

4: Initialization: $\mathcal{G}^{(1)}=\mathcal{G}^{(1),k},\
\mathcal{V}_{1}=\mathcal{M}_{1}=\mathcal{N}_{1}=0$

5: while not coverged do

6: Update $\mathcal{R}_{1}$ via (3.8)

7: Update $\mathcal{G}^{(1)}$ via (3.12)

8: Update $\mathcal{V}_{1}$ via (3.14)

9: Update $\mathcal{M}_{1},\ \mathcal{N}_{1}$ via (3.16)

10: end while

11: Output: $\mathcal{G}^{(1),k+1}$

12: Initialization: $\mathcal{G}^{(2)}=\mathcal{G}^{(2),k},\
\mathcal{V}_{2}=\mathcal{M}_{2}=\mathcal{N}_{2}=0$

13: while not coverged do

14: Update $\mathcal{R}_{2}$ via (3.18)

15: Update $\mathcal{G}^{(2)}$ via (3.19)

16: Update $\mathcal{V}_{2}$ via (3.20)

17: Update $\mathcal{M}_{2},\ \mathcal{N}_{2}$ via (3.21)

18: end while

19: Output: $\mathcal{G}^{(2),k+1}$

20: Initialization: $\mathcal{G}^{(3)}=\mathcal{G}^{(3),k},\
\mathcal{V}_{3}=\mathcal{M}_{3}=\mathcal{N}_{3}=0$

21: while not coverged do

22: Update $\mathcal{R}_{3}$ via (3.24)

23: Update $\mathcal{G}^{(3)}$ via (3.25)

24: Update $\mathcal{V}_{3}$ via (3.26)

25: Update $\mathcal{M}_{3},\ \mathcal{N}_{3}$ via (3.27)

26: end while

27: Output: $\mathcal{G}^{(3),k+1}$

28: Assign $k$ by $k+1$

29: If
$\|\mathbf{\Phi}(\mathcal{G})^{k}-\mathbf{\Phi}(\mathcal{G})^{k-1}\|_{F}/\|\mathbf{\Phi}(\mathcal{G})^{k}\|_{F}\leq\epsilon$,
stop the iteration

30: end while

30: HR-HSI
$\mathcal{X}=\mathbf{\Phi}(\mathcal{G}^{(1),k},\mathcal{G}^{(2),k},\mathcal{G}^{(3),k})$

##  4 Experiments

In this section, we commence by providing an overview of the datasets utilized
in the conducted experiments. Following that, we introduce five prominent
quantitative metrics that have been employed to evaluate the quality of the
fusion outcomes. Furthermore, we delve into an analysis of the impact of
various parameters on the experimental results, subsequently presenting a list
of the identified optimal parameters. Finally, we show the experimental
results obtained by the proposed method (named LogLRTR) and provide a
convergence analysis. Specifically, to establish the superiority of our
approach, we compare it with three advanced methods, including coupled sparse
tensor factorization (CSTF) [19], coupled TR factorization (CTRF) [21], and
factor smoothed tensor ring decomposition (FSTRD) [30]. All experiments were
performed on a desktop computer equipped with a 2.70 GHz Intel(R) Core(TM)
i5-11400U CPU and 16 GB RAM using MATLAB (R2018a).

###  4.1 Datasets

The University of Pavia dataset [50] was acquired using the ROSIS sensor
during an aerial campaign conducted over the University of Pavia, Italy. The
HSI encompasses $115$ spectral bands, encompassing an image dimension of
$610\times 340$ pixels. After excluding bands with low signal-to-noise ratio
(SNR) and delineating the specific subregion, we select the upper-left 256×256
spatial pixels containing 93 spectral bands as our HR-HSI. Employing a
Gaussian blur (with a size of $7\times 7$ and a standard deviation of 2), we
simulated the low-resolution HSI. Subsequently, the image underwent a
downsampling process with a reduction factor of 4. The HR-MSI consisting of
four bands, is generated using the IKONOS-like reflectance spectral response
filter [51].

The CAVE dataset [52], which comprises 32 indoor HSIs captured by a
generalized combined pixel camera within real scenes, has become a prevalent
choice for HSI-MSI fusion [15, 40]. Each HSI within this dataset has
dimensions of $512\times 512\times 31$, with $512\times 512$ indicating
spatial pixels and 31 signifying spectral bands. These 31 spectral bands cover
the wavelength range of $400\sim 700$ nm in 10 nm intervals. For our testing
requirements, we have selected three distinctive HSI scenes (Superballs,
Flowers, and Peppers) as the primary real-world dataset. To derive the
corresponding LR-HSI, we apply a symmetric Gaussian blur with a standard
deviation of 2, utilizing a $7\times 7$ kernel to convolve the HR-HSI. During
the downsampling process, we reduce the spatial resolution by selecting every
8 pixels in both spatial dimensions for each HSI band. To generate the HR-MSI,
we employ spectral downsampling on the HR-HSI from the CAVE dataset, utilizing
the spectral response matrix associated with the Nikon D700 camera.

In addition, we add Gaussian noise with SNRs of 25dB and 30dB to the LR-HSI
and the HR-MSI, respectively.

###  4.2 Quantitative indicators

In order to fully assess the quality of HR-HSI obtained by fusing LR-HSI and
HR-MSI, five quality metrics are used in our study, including the average peak
signal-to-noise ratio (PSNR) [21], structural similarity index (SSIM) [53],
relative dimensionless global error in synthesis (ERGAS) [54], spectral angle
mapper (SAM) [55], and universal image quality index (UIQI) [56]. Among these
metrics, higher values of SSIM, PSNR, and UIQI indicate superior quality,
while lower values of ERGAS and SAM correspond to enhanced quality measures.
Before initiating the simulation, the pixel values of the HR-HSI are
normalized within the range [0, 1]. During the computation of quality metrics,
all image elements are subsequently rescaled to span the range [0, 255].

###  4.3 Parameters analysis

Within our devised method, there are six regularization parameters $\lambda$,
$\alpha$, $\beta$, $\mathrm{R}_{1}$, $\mathrm{R}_{2}$ and $\mathrm{R}_{3}$, as
well as two algorithmic parameters $\eta$, $\mu$. Notably, $\mathrm{R}_{1}$,
$\mathrm{R}_{2}$ and $\mathrm{R}_{3}$ make up the TR rank $\mathrm{R}$.
Initially, it is imperative to establish a sensible choice for the TR rank.
Subsequently, a stepwise adjustment is conducted for the remaining five
parameters. While fine-tuning these parameters, our approach began with
informed estimations for their initial values, relying on our prior expertise.
We iteratively fine-tuned the parameters using a trial-and-error methodology,
systematically searching for the optimal parameter configuration that yielded
the highest peak signal-to-noise ratio (PSNR) value.

Analysis of TR rank: In experiments, we find that the TR rank is usually
needed into a pattern that sets the two sides of the rank smaller and the
middle larger. To simplify the complexity of the parameters analysis, we
choose $\mathrm{R}_{3}=\mathrm{R}_{1}$ to maintain the consistency of the
spatial TR core tensor. $\mathrm{R}_{1}$ and $\mathrm{R}_{3}$ undergo
adjustments within the range of 2 to 6, while $\mathrm{R}_{2}$ is varied
within the interval of 50 to 300. Fig. 1 shows the PSNR curves versus TR rank
$[\mathrm{R}_{1},\mathrm{R}_{2},\mathrm{R}_{1}]$ for the reconstructed Pavia
University dataset. From this figure, we find that a larger value should be
chosen for both $\mathrm{R}_{1}$ and $\mathrm{R}_{2}$ to obtain satisfactory
results.

Figure 1: Sensitivity analysis of TR rank.

Specifically, Table 1 lists the optimal parameter sets chosen for the LogLRTR.

Table 1: Selected parameter sets for the proposed method. Image | $\mathrm{R}$ | $\lambda$ | $\alpha$ | $\beta$ | $\eta$ | $\mu$
---|---|---|---|---|---|---
Pavia University | [6,300,6] | 0.01 | $1\cdot 10^{-5}$ | 3 | 20 | 0.2
Flowers | [5,300,5] | 1.3 | $7\cdot 10^{-4}$ | 1 | 0.8 | 0.04
Superballs | [5,300,5] | 1 | $3\cdot 10^{-5}$ | 0.8 | 0.7 | 0.04
Peppers | [5,300,5] | 1.2 | $3\cdot 10^{-4}$ | 0.3 | 0.9 | 0.08

###  4.4 Experimental results

The numerical results of various fusion approaches for “Pavia University”,
“Flowers”, “Superballs”, and “Peppers” are displayed in Table 2, with the
optimal values highlighted in bold. As shown in Table 2, LogLRTR outperforms
other comparative methods in all five evaluation indicators. This suggests
that our approach substantially improves both the spatial and spectral
characteristics of the fused HR-HSI in the Pavia and CAVE datasets. In
particular, across the four images, LogLRTR exhibits a PSNR advantage over
FSTRD by 0.76 dB, 0.69 dB, 0.76 dB, and 0.72 dB, respectively.

Table 2: Quantitative evaluation of four different fusion methods.

Image | Pavia University | Flowers |
---|---|---|---
Method | PSNR | SSIM | ERGAS | SAM | UIQI | PSNR | SSIM | ERGAS | SAM | UIQI
Best Values | $+\infty$ | 1 | 0 | 0 | 1 | $+\infty$ | 1 | 0 | 0 | 1
CSTF | 39.462 | 0.963 | 1.605 | 3.119 | 0.983 | 42.640 | 0.963 | 1.445 | 14.350 | 0.718
CTRF | 39.797 | 0.971 | 1.555 | 2.922 | 0.986 | 41.934 | 0.952 | 3.995 | 18.832 | 0.682
FSTRD | 40.272 | 0.973 | 1.468 | 2.717 | 0.987 | 43.867 | 0.971 | 1.210 | 11.979 | 0.728
LogLRTR | 41.035 | 0.977 | 1.346 | 2.472 | 0.989 | 44.557 | 0.976 | 1.127 | 10.839 | 0.746
Image | Superballs | Peppers |
Method | PSNR | SSIM | ERGAS | SAM | UIQI | PSNR | SSIM | ERGAS | SAM | UIQI
Best Values | $+\infty$ | 1 | 0 | 0 | 1 | $+\infty$ | 1 | 0 | 0 | 1
CSTF | 43.809 | 0.970 | 1.780 | 9.424 | 0.762 | 43.754 | 0.974 | 1.466 | 8.557 | 0.799
CTRF | 43.088 | 0.959 | 1.995 | 12.000 | 0.753 | 42.305 | 0.959 | 1.782 | 11.876 | 0.760
FSTRD | 44.690 | 0.972 | 1.619 | 8.948 | 0.765 | 45.596 | 0.984 | 1.247 | 6.687 | 0.842
LogLRTR | 45.458 | 0.976 | 1.474 | 8.546 | 0.791 | 46.317 | 0.986 | 1.101 | 6.110 | 0.850

In addition to the numerical analysis, we depict the PSNR and UIQI versus band
curves for the various methods in Fig. 2. This graphical representation offers
compelling evidence that LogLRTR attains higher PSNR and UIQI values across
almost all spectral bands when compared to other approaches. Particularly, the
“Pavia University” and “Flowers” images display substantial enhancements. In
summation, these findings prove the superior performance of our proposed
method.

To further validate the ability of LogLRTR on maintaining spatial structures,
we show the fused and error images acquired by different methods for each test
image. The error image effectively captures the disparities between the fused
image and the ground truth reference image. In addition, to facilitate visual
comparison, we mark and enlarge meaningful regions from each fused the false
color image.

(a) Pavia University

(b) Flowers

(c) Superballs

(d) Peppers

Figure 2: The plots of PSNR (top) and UIQI (bottom) values versus the spectral
band for four test images.

In Fig. 3, we show the false color images of “Pavia University” generated
using bands 83, 38 and 34 obtained through different fusion methods, as well
as the corresponding error images. In the magnified sub-region of the false
color image, the CSTF blurs the image details, while the roofs of the other
methods are closer to the original image, proving once again the drawback of
Tucker decomposition that does not sufficiently exploit the intrinsic
information of the HSI. From the second row of Fig. 3, it can be more clearly
observed that LogLRTR achieves a smaller error, indicating that our proposed
method exhibits an advantage in preserving most of the spatial details.

(a) CSTF

(b) CTRF

(c) FSTRD

(d) LogLRTR

(e) Ground truth

Figure 3: The first row presents the false color images formed by bands 83, 38
and 34 for “Pavia University”, and the second row displays the corresponding
error images.

Additionally, we assess the reconstructed images of “Flowers”, “Superballs”
and “Peppers” in the CAVE dataset and the visualisation of their error maps to
verify the superiority of LogLRTR. Fig. 4 provides the false color images
formed by bands 28, 12 and 5 in the “Flowers” fusion images. Specifically, the
enlarged regions show that some important details are lost in the results of
CSTF and CTRF. The proposed method and FSTRD perform better in reconstructing
high-resolution details. In other words, more detailed edge structures are
preserved through these two methods. In addition, we clearly observe from the
corresponding error images that LogLRTR produces less error than FSTRD.

The false color images consisting of bands 29, 13 and 5 of “Superballs” and
the corresponding residual images are displayed in Fig. 5. The magnified
regions of the fused images show that LogLRTR generates smoother result and
less parallax compared to other approaches. This conclusion is also supported
by the fact that LogLRTR retains relatively more edge details in the error
image.

(a) CSTF

(b) CTRF

(c) FSTRD

(d) LogLRTR

(e) Ground truth

Figure 4: The first row gives the false color images formed by bands 28, 12
and 5 for “Flowers”, and the second row shows the corresponding error images.

Fig. 6 presents the false color images comprising of the 29th, 13th and 5th
bands of “Peppers” reconstructed by various fusion approaches. The magnified
regions show an evidence that the previous three methods generate artifacts on
the surface of the pepper, while the fused image acquired by LogLRTR is closer
to the ground truth image. The error images clearly illustrate that CTRF
produces less satisfactory fusion result, confirming its suboptimal
performance on the CAVE dataset. Conversely, LogLRTR showcases a notable
improvement, especially for the pepper in the upper left corner.

(a) CSTF

(b) CTRF

(c) FSTRD

(d) LogLRTR

(e) Ground truth

Figure 5: The first row shows the false color images formed by bands 29, 13
and 5 for “Superballs”, and the second row presents the corresponding error
images.

(a) CSTF

(b) CTRF

(c) FSTRD

(d) LogLRTR

(e) Ground truth

Figure 6: The first row presents the false color images formed by bands 29, 13
and 5 for “Peppers”, and the second row displays their corresponding error
images.

###  4.5 Ablation experiments

To highlight the efficacy of the three TR factor LTNN terms, we carry out
ablation experiments. This involves disabling in turn the spectral and spatial
TR factors in the proposed method, providing additional insights into their
individual contributions. More specifically, we ban the use of the one
spectral TR factor LTNN term and the two spatial TR factor LTNN terms, which
are respectively termed as LogLRTR$\cdot$ban$\cdot$spe and
LogLRTR$\cdot$ban$\cdot$spa. The numerical results of ablation experiments
conducted on “Flowers” and “Superballs” by using FSTRD, the proposed LogLRTR
and its variants are shown in Table 3. We easily observe that LogLRTR
surpasses the other three approaches, and the performance of the proposed
method benefits significantly from the inclusion of the three TR factor LTNN
terms.

Table 3: Comparison of quantitative indices from ablation experiments. Image | Method | PSNR | SSIM | ERGAS | SAM | UIQI
---|---|---|---|---|---|---
Flowers | FSTRD | 43.867 | 0.971 | 1.210 | 11.979 | 0.728
| LogLRTR$\cdot$ban$\cdot$spe | 44.172 | 0.973 | 1.162 | 11.586 | 0.733
| LogLRTR$\cdot$ban$\cdot$spa | 44.148 | 0.971 | 1.210 | 12.091 | 0.729
| LogLRTR | 44.557 | 0.976 | 1.127 | 10.839 | 0.746
Superballs | FSTRD | 44.690 | 0.972 | 1.619 | 8.948 | 0.765
| LogLRTR$\cdot$ban$\cdot$ spe | 45.196 | 0.975 | 1.497 | 8.404 | 0.781
| LogLRTR$\cdot$ban$\cdot$ spa | 45.269 | 0.975 | 1.504 | 8.571 | 0.770
| LogLRTR | 45.458 | 0.976 | 1.474 | 8.546 | 0.791

###  4.6 Convergence

Finally, we empirically demonstrate the convergence of the newly developed
algorithm. Specifically, all four images are tested by Algorithm 3.1. The
plots of relative error and PSNR versus the number of iterations are displayed
in Fig. 7. The curves depicted in Fig. 7(a) reveal a consistent trend of the
relative errors converging towards zero. Furthermore, the PSNR values grow
quickly as the number of iterations increases, eventually reaching a stable
state, as demonstrated in Fig. 7(b). These observations collectively
constitute robust evidence of the algorithm’s convergence.

Figure 7: The plots of relative error
($\left\|\mathcal{X}^{k}-\mathcal{X}^{k-1}\right\|_{F}/\left\|\mathcal{X}^{k}\right\|_{F}$)
(left) and PSNR (right) against the number of iterations for four images.

##  5 Conclusion

In this paper, we introduced a newly developed hyperspectral image fusion
model that integrates logarithmic tensor nuclear norm (LTNN) with weighted TV
of TR factors. To better preserve the correlation in the spatial-spectral
domain, the TR decomposition was utilized to represent the HR-HSI. Using this
TR framework as a foundation, the objective of HSI-MSI fusion is reframed as
the estimation of the three TR factors. To enhance the continuity of TR
factors for improved maintenance of the smooth structure, the weighted total
variation with weight updates was used. In addition, we employed the mode-2
LTNN to constrain each TR factor for delving deeper into the low-rank
characteristics. Numerical experiments showed that our proposed approach is
superior to existing HSI-MSI fusion methods. In future work, we intend to
develop some adaptive techniques for selecting the appropriate TR rank.

## Acknowledgments

This work was partly supported by the Major Program of the National Natural
Science Foundation of China (No. 61890962), the Natural Science Foundation of
Jiangxi Province (20232BAB201017), the Science Foundation for Post Doctorate
of China (2020M672484), and Jiangxi Postgraduate Innovation Special Fund
Project (YC2022-S987).

## References

  * [1] Li S, Kang X, Fang L, et al. Pixel-level image fusion: A survey of the state of the art[J]. Information Fusion, 2017, 33: 100-112.
  * [2] Qureshi R, Uzair M, Khurshid K, et al. Hyperspectral document image processing: Applications, challenges and future prospects[J]. Pattern Recognition, 2019, 90: 12-22.
  * [3] Zheng J, Feng Y, Bai C, et al. Hyperspectral image classification using mixed convolutions and covariance pooling[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 59(1): 522-534.
  * [4] Hong D, Wu X, Ghamisi P, et al. Invariant attribute profiles: A spatial-frequency joint feature extractor for hyperspectral image classification[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(6): 3791-3808.
  * [5] Li J Y, Xie J Y, Luo Y S, et al. H2TF for Hyperspectral Image Denoising: Where Hierarchical Nonlinear Transform Meets Hierarchical Matrix Factorization[J]. IEEE Geoscience and Remote Sensing Letters, 2023, 20: 5506705.
  * [6] Huang J, Huang T Z, Zhao X L, et al. Nonlocal tensor-based sparse hyperspectral unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 59(8): 6854-6868.
  * [7] Kanatsoulis C I, Fu X, Sidiropoulos N D, et al. Hyperspectral super-resolution: A coupled tensor factorization approach[J]. IEEE Transactions on Signal Processing, 2018, 66(24): 6503-6517.
  * [8] Dian R, Fang L, Li S. Hyperspectral image super-resolution via non-local sparse tensor factorization[C]. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017: 5344-5353.
  * [9] Wei Q, Bioucas-Dias J, Dobigeon N, et al. Multiband image fusion based on spectral unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(12): 7236-7249.
  * [10] Ran R, Deng L J, Jiang T X, et al. GuidedNet: A general CNN fusion framework via high-resolution guidance for hyperspectral image super-resolution[J]. IEEE Transactions on Cybernetics, 2023, 53(7): 4148-4161.
  * [11] Deng L J, Vivone G, Paoletti M E, et al. Machine Learning in Pansharpening: A benchmark, from shallow to deep networks[J]. IEEE Geoscience and Remote Sensing Magazine, 2022, 10(3): 279-315.
  * [12] Deng S Q, Deng L J, Wu X, et al. PSRT: Pyramid shuffle-and-reshuffle transformer for multispectral and hyperspectral image fusion[J]. IEEE Transactions on Geoscience and Remote Sensing, 2023, 61: 1-15.
  * [13] Liu J, Wu Z, Xiao L, et al. A truncated matrix decomposition for hyperspectral image super-resolution[J]. IEEE Transactions on Image Processing, 2020, 29: 8028-8042.
  * [14] Simões M, Bioucas‐Dias J, Almeida L B, et al. A convex formulation for hyperspectral image superresolution via subspace-based regularization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 53(6): 3373-3388.
  * [15] Dian R, Li S, Fang L, et al. Hyperspectral image super-resolution via local low-rank and sparse representations[C]. IEEE International Geoscience and Remote Sensing Symposium, 2018: 4003-4006.
  * [16] Zhang J, Liu Z, Ma M. Hyperspectral image fusion with a new hybrid regularization[J]. Computational and Applied Mathematics, 2022, 41(6): 241.
  * [17] Zhang J, Lu J, Wang C, et al. Hyperspectral and multispectral image fusion via superpixel-based weighted nuclear norm minimization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2023, 61: 5521612.
  * [18] Kanatsoulis C I, Fu X, Sidiropoulos N D, et al. Hyperspectral super-resolution: A coupled tensor factorization approach[J]. IEEE Transactions on Signal Processing, 2018, 66(24): 6503-6517.
  * [19] Li S, Dian R, Fang L, et al. Fusing hyperspectral and multispectral images via coupled sparse tensor factorization[J]. IEEE Transactions on Image Processing, 2018, 27(8): 4118-4130.
  * [20] Dian R, Li S, Fang L. Learning a low tensor-train rank representation for hyperspectral image super-resolution[J]. IEEE Transactions on Neural Networks and Iearning Systems, 2019, 30(9): 2672-2683.
  * [21] He W, Chen Y, Yokoya N, et al. Hyperspectral super-resolution via coupled tensor ring factorization[J]. Pattern Recognition, 2022, 122: 108280.
  * [22] Xu Y, Wu Z, Chanussot J, et al. Nonlocal coupled tensor CP decomposition for hyperspectral and multispectral image fusion[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 58(1): 348-362.
  * [23] Xu T, Huang T Z, Deng L J, et al. Hyperspectral image superresolution using unidirectional total variation with tucker decomposition[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2020, 13: 4381-4398.
  * [24] Wang K, Wang Y, Zhao X L, et al. Hyperspectral and multispectral image fusion via nonlocal low-rank tensor decomposition and spectral unmixing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(11): 7654-7671.
  * [25] Yang J H, Zhao X L, Ji T Y, et al. Low-rank tensor train for tensor robust principal component analysis[J]. Applied Mathematics and Computation, 2020, 367: 124783.
  * [26] Zhao Q, Zhou G, Xie S, et al. Tensor ring decomposition[J]. ArXiv Preprint ArXiv:1606.05535, 2016.
  * [27] Wang W, Aggarwal V, Aeron S. Efficient low rank tensor ring completion[C]. Proceedings of the IEEE International Conference on Computer Vision, 2017: 5697-5705.
  * [28] He W, Yokoya N, Yuan L, et al. Remote sensing image reconstruction using tensor ring completion and total variation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(11): 8998-9009.
  * [29] Yuan L, Li C, Mandic D, et al. Tensor ring decomposition with rank minimization on latent space: An efficient approach for tensor completion[C]. Proceedings of the AAAI Conference on Artificial Intelligence, 2019, 33(01): 9151-9158.
  * [30] Chen Y, Zeng J, He W, et al. Hyperspectral and Multispectral Image Fusion Using Factor Smoothed Tensor Ring Decomposition[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 1-17.
  * [31] Xu Y, Wu Z, Chanussot J, et al. Hyperspectral images super-resolution via learning high-order coupled tensor ring representation[J]. IEEE Transactions on Neural Networks and Iearning Systems, 2020, 31(11): 4747-4760.
  * [32] Xu H, Qin M, Chen S, et al. Hyperspectral-multispectral image fusion via tensor ring and subspace decompositions[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 8823-8837.
  * [33] Chen Y, He W, Yokoya N, et al. Nonlocal tensor-ring decomposition for hyperspectral image denoising[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 58(2): 1348-1362.
  * [34] Wang W, Sun Y, Eriksson B, et al. Wide compression: Tensor ring nets[C]. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018: 9329-9338.
  * [35] Xue J, Zhao Y, Liao W, et al. Nonlocal low-rank regularized tensor decomposition for hyperspectral image denoising[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(7): 5174-5189.
  * [36] Cao C, Yu J, Zhou C, et al. Hyperspectral image denoising via subspace-based nonlocal low-rank and sparse factorization[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2019, 12(3): 973-988.
  * [37] Lu C, Feng J, Chen Y, et al. Tensor robust principal component analysis with a new tensor nuclear norm[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2019, 42(4): 925-938.
  * [38] Hu W, Tao D, Zhang W, et al. The twist tensor nuclear norm for video completion[J]. IEEE Transactions on Neural Networks and Learning Systems, 2016, 28(12): 2961-2973.
  * [39] Ding M, Huang T Z, Zhao X L, et al. Tensor completion via nonconvex tensor ring rank minimization with guaranteed convergence[J]. Signal Processing, 2022, 194: 108425.
  * [40] Dian R, Li S. Hyperspectral image super-resolution via subspace-based low tensor multi-rank regularization[J]. IEEE Transactions on Image Processing, 2019, 28(10): 5135-5146.
  * [41] Rivenson Y, Stern A. Compressed imaging with a separable sensing operator[J]. IEEE Signal Processing Letters, 2009, 16(6): 449-452.
  * [42] Wang Y, Peng J, Zhao Q, et al. Hyperspectral image restoration via total variation regularized low-rank tensor decomposition[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2017, 11(4): 1227-1243.
  * [43] Attouch H, Bolte J, Redont P, et al. Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality[J]. Mathematics of Operations Research, 2010, 35(2): 438-457.
  * [44] Yang J H, Zhao X L, Ma T H, et al. Remote sensing images destriping using unidirectional hybrid total variation and nonconvex low-rank regularization[J]. Journal of Computational and Applied Mathematics, 2020, 363: 124-144.
  * [45] Liu Y Y, Zhao X L, Zheng Y B, et al. Hyperspectral image restoration by tensor fibered rank constrained optimization and plug-and-play regularization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2021, 60: 1-17.
  * [46] Candes E J, Wakin M B, Boyd S P. Enhancing sparsity by reweighted $\ell_{1}$ minimization[J]. Journal of Fourier Analysis and Applications, 2008, 14: 877-905.
  * [47] Xie Q, Zhao Q, Meng D, et al. Multispectral images denoising by intrinsic tensor sparsity regularization[C]. IEEE Conference on Computer Vision and Pattern Recognition, 2016: 1692-1700.
  * [48] Zheng Y B, Huang T Z, Zhao X L, et al. Mixed noise removal in hyperspectral image via low-fibered-rank regularization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 58(1): 734-749.
  * [49] Zhao Q, Sugiyama M, Yuan L, et al. Learning efficient tensor representations with ring-structured networks[C]. IEEE International Conference on Acoustics, Speech and Signal Processing, 2019: 8608-8612.
  * [50] Dell’Acqua F, Gamba P, Ferrari A, et al. Exploiting spectral and spatial information in hyperspectral urban data with high resolution[J]. IEEE Geoscience and Remote Sensing Letters, 2004, 1(4): 322-326.
  * [51] Wei Q, Dobigeon N, Tourneret J Y. Fast fusion of multi-band images based on solving a Sylvester equation[J]. IEEE Transactions on Image Processing, 2015, 24(11): 4109-4121.
  * [52] Yasuma F, Mitsunaga T, Iso D, et al. Generalized assorted pixel camera: postcapture control of resolution, dynamic range, and spectrum[J]. IEEE transactions on image processing, 2010, 19(9): 2241-2253.
  * [53] Wang Z, Bovik A C, Sheikh H R, et al. Image quality assessment: from error visibility to structural similarity[J]. IEEE Transactions on Image Processing, 2004, 13(4): 600-612.
  * [54] Wald L. Quality of high resolution synthesised images: Is there a simple criterion?[C]. Third conference” Fusion of Earth data: merging point measurements, raster maps and remotely sensed images”, 2000: 99-103.
  * [55] Yuhas R H, Goetz A F H, Boardman J W. Discrimination among semi-arid landscape endmembers using the spectral angle mapper (SAM) algorithm[C]. JPL, Summaries of the Third Annual JPL Airborne Geoscience Workshop, 1992: 147–149.
  * [56] Wang Z, Bovik A C. A universal image quality index[J]. IEEE Signal Processing Letters, 2002, 9(3): 81-84.