Low-Rank Tensor Decomposition With Smooth and Sparse Regularization for Hyperspectral and Multispectral Data Fusion

The fusion of hyperspectral and multispectral images is an effective way to obtain hyperspectral super-resolution images with high spatial resolution. A hyperspectral image is a datacube containing two spatial dimensions and a spectral dimension. The fusion methods based on non-negative matrix factorization need to reshape the three-dimensional data in matrix form, which will result in the loss of data structure information. Owing to the non-uniqueness of tensor rank and noise inference, there is a lot of redundant information in the spatial and spectral subspaces of tensor decomposition. To address the above problems, this article incorporates smooth and sparse regularization into low-rank tensor decomposition to reformulate a fusion method, in which the logarithmic sum function is adopted to eliminate the effect of redundant information and shadows in both spatial and spectral domains. Moreover, the total-variation-based regularizer is employed to vertically smooth the spectral factor matrix to suppress the noise. Then, the alternating direction multiplier method, as well as the conjugate gradient approach, is utilized to design a set of efficient algorithms by complexity reduction. The experimental results demonstrate that the proposed method can yield better performance than the state-of-the-art benchmark algorithms in most cases, which also verifies the effectiveness of incorporated regularizers in low signal-to-noise ratio environments for hyperspectral super-resolution images.


I. INTRODUCTION
Ahyper-spectral image (HSI) is a three-dimensional datacube with abundant spatial and spectral information, which is widely used to accurately identify substance and classification in remote sensing [1]- [3]. Generally speaking, HSIs have high spectral resolution but low spatial resolution owing to their mutual restriction. Therefore, how to obtain a hyperspectral super-resolution image is of vital importance in practical applications. In recent years, as an economical and effective way, data fusion of low spatial resolution HSI (LR-HSI) and high spatial resolution multispectral image (HR-MSI) has become a research hotspot in remote sensing The associate editor coordinating the review of this manuscript and approving it for publication was Wenming Cao .
[4]- [6]. Many related research results have been obtained on data fusion of HSIs and MSIs, which can be roughly divided into four categories, including panchromatic sharpening (pansharpening), deep neural network, subspace-based fusion, and tensor decomposition.
Pansharpening is extended to fuse the LR-HSI and HR-MSI by the spectral extension and projection [6], although it was originally developed for the fusion of a panchromatic image and an MSI, whose representative methods include component substitution (CS) [7], [8] and multi-resolution analysis (MRA) [9]. For instance, the CSbased pansharpening method [8] was proposed to sharpen the MSI by adding the spatial details obtained from the difference of a panchromatic image and a synthetic intensity component. The MRA-based method was designed to generate more spatial details from each multispectral band by multiplying the difference between a panchromatic image and its low-pass version through a filter, in which the spatial information can be gained in the wavelet or Laplace pyramid domain. However, spectral distortion is usually caused by ignoring the inherent correlation between spatial and spectral signatures, although some methods (i.e. GLP and SFIM) work very efficiently [6]. With the development of deep learning, the hyperspectral and multispectral data or features are trained by building a deep neural network model, in which the network hyperparameters are determined for the reconstructed super-resolution image by relevant constraints [10]- [13]. For example, considering spatial smoothing and spectral correlation, a hyperspectral image super-resolution based on convolutional neural network (CNN) [10] was proposed to alleviate the distortion. Yang et al. [11] proposed a fusion method based on deep CNN, which used the spectral correlation to fuse the two branches of the LR-HSI and HR-MSI separately. However, these methods need a training library for data or features, and involve a large number of hyperparameters adjustments.
The subspace-based fusion methods include Bayesian inference and spectral unmixing [5]. As mentioned above, data fusion aims at integrating the observed LR-HSI and HR-MSI in the subspaces to obtain the hyperspectral super-resolution image in the original space. Conversely, the observed LR-HSI and HR-MSI can be generated by spatially and spectrally degrading the target image, respectively, which is called the observation model. Thus, the subspacebased fusion method is usually an ill-posed inverse problem, which needs prior knowledge or regularization to limit the range of solution. A Bayesian-based method was proposed by adding the priors (eg. Guassian prior) to handle the fusion problem in the wavelet domain or principal component subspace [14]. Wei et al. [15] decomposed the images on a set of dictionaries by redefining a proper posterior distribution, and proposed a Bayesian-based sparse representation (BSR) algorithm. Incorporating the total variation regularizer, Simoes et al. [16] developed the hyperspectral super-resolution (HySure), which was solved by the split augmented Lagrangian shrinkage algorithm. However, since these methods generally involve large computation, complexity reduction is required to simplify the Bayesian-based fusion approaches, such as a Sylvester equation-based fast fusion [17].
Another subspace method is based on spectral unmixing, and its representative work is called couple non-negative matrix factorization (CNMF) [18]. From the viewpoint of the observational model, CNMF, like Bayesian inference, was also an inverse problem. As is well known, incorporating the regularization of endmember and/or abundance signature is intended to eliminate the ill-posedness. Sparsitypromoting regularizers were used to yield high-quality fused images [19]- [22]. Furthermore, joint spatial and spectral regularization can improve the fused performance [23], [24]. Wu et al. [25] formulated hyperspectral super-resolution as a global-local rank-regularized least-squares problem, assuming that the spectral-spatial matrices associated with the whole image and the local areas of the image had low-rank structures. Veganzones et al. [26] partitioned the image into patches, and proposed a method to solve the data fusion problem independently for each patch based on the fact that the real data HSIs are locally of low rank. Yang et al. [24] employed total variation and signature-based regularization, in which the horizontal and vertical difference matrices were constructed separately for spatial smoothing. However, besides the ill-posedness, the spectral-unmixing-based fusion methods need to convert the datacube of HSI and MSI in matrix form, which leads to the loss of three-dimensional structure information.
From the perspective of tensor decomposition, the data of LR-HSI and HR-MSI can be regarded as three-dimensional tensors. To take full advantage of data structure information, it is a better way to utilize the tensor decomposition to fuse the HSI and MSI [27]- [31]. Tensor models include canonical polyadic (CP) decomposition, Tucker decomposition, tensor ring, and others. Owing to the low-rank property, tensor decomposition has been widely used in remote sensing, such as hyperspectral image restoration [32], image smooth and denoising [33], [34], feature extraction [35], image compression and reconstruction [36], [37]. For instance, Xue et al. imposed the nonlocal low-rank regularization on tensor decomposition for hyperspectral image denoising [34]. Certainly, there are also some related research works on tensor-based image fusion. For example, Kanatsoulis et al. [27] proposed a coupled tensor factorization framework (called STEREO) to overcome the aforementioned problems, which needed little information of the degradation operators. Xu et al. [28] brought forward a nonlocal tensor decomposition model for HSI-MSI fusion by exploring the relationship between the LR-HSI and the HR-MSI via a coupled CP decomposition. Li et al. [29] redefined the fusion problem as an estimation of the core tensor and dictionaries of the three modes, and proposed a coupled sparse tensor factorization (CSTF) method. Dian et al. [38] proposed a novel HSI super-resolution method based on non-local sparse tensor factorization (called NLSTF) via Tucker decomposition to further exploit the non-local spatial self-similarities of the HSI. In order to preserve spatial-spectral structures in the LR-HSIs and the HR-MSIs effectively, Zhang et al. [30] put forward an image fusion method based on spatial-spectral-graph-regularized low-rank tensor decomposition to preserve the spatial correlation and the spectral structure in the fused images by deriving two graphs in the spatial and spectral domains, respectively. Dian et al. [39] developed a subspace-based low-tensor method with multi-rank regularization, which fully exploited the spectral correlations and non-local similarities in the HR-HSI. He et al. [40], [41] developed a hyperspectral superresolution based on a coupled tensor ring model to decompose a higher-order tensor into a series of three-dimensional tensors. Xu et al. [42] proposed a super-resolution method using a VOLUME 8, 2020 high-order coupled tensor ring for the fusion of LR-HSI and HR-MSI, which maintained spectral information and a core tensor in the tensor ring to reconstruct the high-resolution image. However, coupled tensor decomposition methods often face a problem of the redundant information caused by the non-uniqueness of tensor rank, and their performance decreases rapidly in the low signal-to-noise ration (SNR) environment.
To address the effect of redundant information and noise, this article proposes a novel fusion method based on low-rank Tucker tensor decomposition to fuse the LR-HSI and HR-MSI for the hyperspectral super-resolution image. Our main contributions in this article include the aspects dicussed below. First of all, on the basis of the tensor observation model, we propose a low-rank tensor fusion method with total variation and sparse regularization, in which the logarithmic sum function is incorporated to strengthen the low-rank property for removing the redundancy of factor matrices. Then, because of the gradual change in endmember curves, total variation is utilized to suppress the noise by smoothing the spectral factor matrix vertically. Next, containing the weights of three-dimensional factor matrices, the core tensor is generally sparse, since Tucker decomposition can be regarded as a higher-order form of singular value decomposition (SVD). Finally, a set of efficient algorithms is carefully designed to reduce the computational complexity using the conjugate gradient (CG) and alternating direction multiplier method (ADMM) [43]. Three real datasets captured by different sensors are adopted to conduct the exprimental tests, showing that the proposed method can significantly enhance the spatial resolution of reconstructed images.
The remainder of this paper is organized as follows: Section II develops a tensor-based observation model. Section III reformulates a regularized fusion problem. Section IV designs the solving algorithms and simplifies the solutions to reduce the computational complexity. Section V conducts the experimental tests and evaluates the performance of the proposed algorithm. Section VI draws a conclusion. In addition, Some key notations are used in the subsequent sections. R, R n and R m×n denote the set of real number, n-vector and m × n matrices, respectively. R + , R n + and R m×n + represent the set of non-negative real number, nvector and m × n matrices, respectively. · 1 stands for l 1norm. · F denotes Frobenius norm. I stands for the identity matrix with proper dimensions. vec(C) stands for the vector of a tensor C. ⊗ represents the Kronecker product.
[·] + stands for the orthogonal projection onto the non-negative orthant of Euclidean space. represents the componentwise inequality operation. × n denotes the mode-n product. · stands for the simplified form of mode-n product. · * denotes the nuclear norm.

II. TENSOR OBSERVATION MODEL
Under the framework of convex optimization, this article proposes a low-rank tensor fusion method with smooth and sparse regularization. Firstly, the tensor observation model is established via Tucker decomposition. Then, the smooth and sparse regularization is imposed on the factor matrices and core tensor, respectively.
In this article, the target HR-HSI (called the reference image) is expressed as Z ∈ R N w ×N h ×N a , in which N w , N h and N a denote the sizes of width, height and spectral dimensions, respectively; Y ∈ R n w ×n h ×N a represents the observed LR-HSI, in which n w , n h and N a represent the sizes of horizontal, vertical and spectral dimensions, respectively; X ∈ R N w ×N h ×n a represents the observed HR-MSI with the sizes of N w , N h and n a . Image fusion aims at fusing the observed LR-HSI data Y and HR-MSI data X to reconstruct the target HSI with high spatial and spectral resolution. Conversely, the data Y and X can be regarded as the spatial and spectral degradation versions of the reference image Z, respectively. So the tensor-based observation model is presented as where P 1 ∈ R n w ×N w , P 2 ∈ R n h ×N h and P 3 ∈ R n a ×N a denote the degraded matrices of three dimensions with the downsampling factor of spatial dimensions r = N w n w = N h n h . According to Tucker model, the reference image Z can be decomposed as where W ∈ R N w ×r w , H ∈ R N h ×r h and A ∈ R N a ×r a denote the factor matrices of width, height and spectral dimensions, respectively; C ∈ R r w ×r h ×r a is the core tensor with the sizes of r w , r h , r a ; E z is the residual. If equation (3) is substituted into equations (1) and (2), then the observed tensor model can be given as

III. REFORMULATION OF FUSION PROBLEM
According to the least square criterion for equations (4) and (5), a fusion optimization problem based on coupling tensor decomposition is defined as From the perspective of tensor decomposition, the spatial factor matrices in problem (6) do not satisfy the non-negativity and sum-to-one constraints; the core tensor contains the coefficients of each factor matrix, which is usually sparse. Furthermore, because the non-uniqueness of tensor rank generally result in the redundant information, this article proposes a fusion method based on low-rank tensor decomposition to eliminate the redundancy and suppress the shadow and/or noise, in which smooth and sparse regularizers are imposed on the spectral factor matrix and the core tensor. The novel fusion problem can be reformulated as where φ w (W), φ h (H) and φ a (A) denote the low-rank regularizers for the three-dimensional matrices; φ d (A) is intended to regularize the vertical smoothness of the spectral factor matrix (endmember matrix); φ c (C) is the sparse regularizer of the core tensor; λ w , λ h , λ a , λ d and λ c are the weights.

A. REGULARIZATION EXPRESSIONS
In this subsection, the expressions of the regualrizers are defined and then reformulated in matrix form on the purpose of efficient computation.

1) SPARSE REGULARIZATION
The sparse regularization of core tensor C denotes the number of non-zero elements, which is written as C 0 . Since the solution of the l 0 -norm is an NP-hard problem, it can be determined approximately by the l 1 -norm. Hence, as the surrogate of the l 0 -norm, the l 1 -norm of the core tensor is defined as where |c ijk | is the element of C.

2) SPECTRAL SMOOTH REGULARIZATION
The endmember signature curve describes the spectral reflectance of the endmember, which is usually a smooth (or piecewise smooth) curve owing to the gradual change [44]. Therefore, the total variation is utilized to smooth the spectral signature. We have where γ denotes the vertical-neighbourhood set, and a j (m) is the mth element of the jth endmember a j . In order to facilitate matrix operations, equation (9) should be reformulated in matrix form as where D d is the first-order vertical difference matrix.

3) LOW-RANK REGULARIZATION
The low-rank regularizer is imposed separately on three factor matrices W, H and A. Taking the factor matrix W as an example, the definition and solution method of logarithmic low rank are given below. As is well known, the rank of W is equal to the number of its positive singular value, which is non-convex and hard to solve directly. In general, the nuclear norm is used to calculate the approximate rank of a matrix. On that basis, this article adopts the logarithmic sum function as the surrogate of low-rank property to suppress the small singular values caused by noise or interference. We take the factor matrix W as an example to illustrate the solving procedure. First of all, the nuclear norm of W is defined as where σ i (W) denotes the ith singular value of W. From equation (11), low-rank regularization is to constrict the matrix rank by minimizing its singular value. Furthermore, the bigger singular values in the factor matrix represents the primary components, while the smaller elements are easily affected by noise or redundant information. Thus, the logarithmic sum function [39], [45] can suppress the noise or secondary components carried in the small singular values, given in equation (12) as where ε is a small positive number to avoid the logarithm of zero. Specifically speaking, each low-rank regularizer of φ w (W) and φ h (H) can suppress the effect of spatial redundant information, and the term φ h (A) is employed to remove the spectral shadow.

4) SOLVING PROCEDURE OF LOW-RANK REGULARIZATION
Given a matrix W 0 ∈ R m×n (m < n) and a weight α, the optimization problem of nuclear norm regularization is defined as where W 0 = Udiag(σ 1 , σ 2 , . . . , σ n )V T . Thus, the closedform solution of W is where When the low-rank regularizer is the logarithmic sum function (12), the optimization problem can be defined as where 0 < α, 0 < ε < min( √ α, (α/σ 1 )). The closed-form solution of equation (15) is where

IV. PROPOSED ALGORITHM
This article proposes a novel fusion method based on low-rank tensor decomposition with total variation and sparse regularization (called LRTVS). Incorporating the regularizers (8), (9) and (12), the non-convex multi-variable fusion problem (7) can be solved via alternating optimization (AO), which is decoupled into single-variable convex subproblems. We have In sum, the complete procedure of LRTVS for solving the problem is summarized in Algorithm 1, in which the factor matrices W and H are initialized via the dictionary-update KSVD method (DU-KSVD) [46], and the initial value of the factor matrix A is obtained via the successive projection algorithm (SPA) [47].

A. ESTIMATION OF FACTOR MATRIX W
Since there are two data fitting terms based on the Frobenius norm and a low-rank regularizer in the objective function of equation (17), it is difficult to solve directly. Hence, ADMM is employed to split the variable of W by adding equality constraints. Equation (17) is reformulated as In general, the tensor needs to be converted into a matrix to be computed. According to mode product, the subproblem of (21) can be unfolded in matrix form along mode 1 (horizontal dimension): (1) , V w = (C× 2 H× 3 (P 3 A)) (1) . The augmented Lagrangian function of equation (22) can be presented as where β w1 is a dual variable, and η is the weight of the augmented Lagrangian function. The primal and dual variables W, W 1 and β w1 can be updated in turn until convergence and are written as where j represents the number of iterations. We take the Lagrangian function (23) into equation (24a) and (24b) separately, and we can get the scaled expressions of two original variables as is a convex quadratic programming problem, but it involves large-scale matrix calculations. Thus, according to the least square criterion, a Sylvester equation of W can be obtained as which can be solved via CG for the solution of W.
2) UPDATE W 1 Equation (25b) is a low-rank decomposition problem. Using equation (16), we can get the closed-form solution of W 1 as 129846 VOLUME 8, 2020 where The dual variable β w1 can be updated by equation (24c).
To sum up, the solving procedure of factor matrix W can be completed by alternately iterating the primal and dual variables, as presented in Algorithm 2.
In the same way as W, problem (28) is unfolded in matrix form along mode 2 (vertical dimension) as where (2) . The ADMM and CG are utilized to obtain the solution of each variable. The iteration process is shown as The Sylvester equation of H is expressed as (30), which can be computed via CG.
2) UPDATE H 1 According to equation (16), the closed-form solution of H 1 is presented as where q h i = Q α,ε (σ h i ), and U h diag(σ h 1 , σ h 2 , . . . , σ h n ) V T h is the SVD form of the matrix H j+1 − β h1 j /η.
In sum, the solving algorithm of H is similar to W, so it can also be solved by Algorithm 2.

C. ESTIMATION OF FACTOR MATRIX A
Besides two data-fitting terms in (33), the low-rank and total variation regularizers are also introduced to reformulate the subproblem as Unfolding equation (33) in matrix form along mode 3 (spectral dimension), the variables can be split by adding three equality constraints. Thus, equation (33) can be reshaped as where U a = (C × 1 (P 1 W) × 2 (P 2 H)) (3) , V a = (C × 1 W × 2 H) (3) . In the same way, equation (34) is solved via the ADMM and CG methods to obtain the solution of all variables, which are updated in turn as follows:

1) UPDATE A
Because of the high complexity of matrix operations, the CG method can be used to solve the Sylvester equation (35).

2) UPDATE THE OTHER VARIABLES
where q a i = Q α,ε (σ a i ), and U a diag(σ a 1 , σ a 2 , . . . , σ a n ) V T a is the SVD of A j+1 − β a1 j /η. The complete solver of factor matrix A is displayed in Algorithm 3.

D. ESTIMATION OF CORE TENSOR C
The sparse regularizer is imposed on the core tensor C, and the tensor form of equation (20) is given as If c = vec(C), then c 1 is equivalent to C 1 . After splitting the variables by the equality constraints, problem (37) The closed-form solutions of the primal and dual variables can be updated iteratively by the following equations as where β 1 and β 2 are dual variables, B 1 = A⊗(P 2 H)⊗(P 1 W), It is found that the bottlenecks of equation (39a) and (39b) are the computational complexities of (B T 1 B 1 + 2η I NL ) −1 and (B T 2 B 2 + η I NL ) −1 , respectively. Tensor operators and the matrix structure of B 1 and B 2 are utilized to simplify the solutions. Hence, the optimal computation of (B T 1 B 1 + 2η I NL ) −1 in equation (39a) is rewritten as where U Bi and i (i = 1, 2, 3) are unitary matrices and diagonal matrices by the SVD decomposition of W T P T 1 P 1 W, H T P T 2 P 2 H and A T A, respectively. In the same way, the simplified expression of (B T 2 B 2 + η I NL ) −1 in (39b) is reformulated as whereŪ Bi and¯ i (i = 1, 2, 3) are unitary matrices and diagonal matrices by the SVD decomposition of W T W, H T H and A T P T 3 P 3 A, respectively. The solving procedure of core tensor C is shown in Algorithm 4.

V. EXPERIMENTS AND PERFORMANCE EVALUATION
In this article, Wald's protocol [48] is adopted to design the experiments for the performance evaluation of the LRTVS algorithm, which is shown in Figure 1. This simulation flowchart not only completes the registration of two observed images, but also provides a comparative reference for performance evaluation. Three real datasets, termed as the reference images Z, are degraded spatially and spectrally to generate the observed HSIs Y and MSIs X that act as the input data for fusion algorithms. Then, Algorithm 1 is taken to integrate the observed data Y and X , in which Algorithms 2, 3 and 4 are called to implement the factor matrices W, H, A and the core tensor C until convergence for reconstructing the HR-HSI Z. Thus, the similarity between Z and Z is measured by such quality evaluation metrics [29], [44] as degree of distortion (DD), reconstruction signal-to-noise ratio (RSNR), spectral angle mapper (SAM), root mean squared error (RMSE), erreur relative globale adimensionnelle de synthè se (ERGAS), and structural similarity index for measuring image quality (SSIM). The larger the value of RSNR and SSIM, the better the performance. However, the smaller the values of SAM, DD, RMSE and ERGAS, the less distortion. In addition, six state-of-the-art methods were taken as the benchmark, including CNMF [18], BSR [15], HySure [16], CSTF [29], NLSTF [38], STEREO [27]. Because the pansharpening methods have higher efficiency and much lower performance than that of the aforementioned coupling-based algorithms, some algorithms (such as GLP and SFIM) are not adopted as the benchmark methods. And the deep-learningbased approaches can also be regarded as the baseline method owing to the different running environments. All the methods were tested on a laptop with Intel CPU Core-3210 with 2.5 GHz speed and 16 GB RAM.

A. DATASETS
In this article, three datasets captured by different sensors are adopted for the performance tests. The first dataset was taken by the reflective optics system imaging spectrometer (ROSIS) sensor over Pavia University in northern Italy [49], with 115 spectral bands and spatial resolution of 1.3 meters. After removing the effect of water vapor absorption, 103 spectral bands were reserved from 430 nm to 860 nm. At the same time, the corresponding MSI in the same region was acquired by the multispectral sensor (IKONOS) [50], with the four bands covering 445-516, 516-595, 632-698 and 757-853 nm. Thus, the spectral downsampling matrix is P 3 ∈ R 4×103 . The second dataset was collected by the AVIRIS (Airborne Visible/Infrared Imaging Spectrometer) over the Moffett, California area [51], in which 183 bands were retained. The corresponding MSI was captured by Landsat TM with six bands (Landsat 1-5 and 7) ranging from 400 nm to 2500 nm [16], [23], [52]. The spectral downsampling matrix is P 3 ∈ R 6×183 . The third dataset was obtained by HYDICE (Hyperspectral digital imagery collection experiment) over Washington DC Mall [53]. After preprocessing, 191 spectral bands (from 400 to 2500 nm) were retained, corresponding to the multispectral sensor Landsat TM band 1-5 and 7. The spectral downsampling matrix is P 3 ∈ R 6×191 to generate the observed multispectral data X .

B. PARAMETER SELECTION AND ANALYSIS
In this part, the experimental tests are conducted to determine the parameters, including the weights of regularizers, the dimensionality of the core tensor, the iteration numbers of ADMM and the SNR.

1) REGULARIZATION PARAMETERS
The low-rank regularization in this article is imposed on the horizontal, vertical and spectral factor matrices with logarithm sum function, in which the weights are λ w , λ h and λ a . Because of symmetry, the weights of two spatial factor matrices satisfy λ w = λ h . Moreover, considering the gradual properties of the spectral signature [44], the total variation is adopted to vertically smooth the spectral factor matrix with the coefficient λ d . The weight of the core tensor is λ c . Generally speaking, the regularization weight is set to balance the relationship between data-fitting terms and the regularizer in the objective function [55], which has a strong correlation with the noise level of the observed image. For example, the corresponding coefficient can be set to 0.001 when the SNR is 30 dB. To determine the optimal weights, we test the RMSE with respect to each weight for three datasets, as shown in sub-figures 2 (a)∼(d). It can be seen from these figures that the RMSE performance fluctuates with the change of the weight for three datasets, while other parameters are set to be a constant value (for instance, 0 or 0.001). Therefore, we select the parameter values with the minimum  error RMSE of three datasets as λ w = λ h = 0.1, λ a = 0.5, λ d = 0.1, λ c = 0.0005.

2) DIMENSIONALITY OF THE CORE TENSOR
A hyperspectral image can be factorized into three factor matrices (i.e. W, H and A) and a core tensor C via Tucker tensor decomposition, in which the core tensor contains the coefficients of factor matrices. Owing to the effect of noise and shadow in the images, the rank of factor matrices is inaccurate. When the dimensionality is smaller than the groundtruth, it will result in serious distortion. Otherwise, it will lead to redundant information in the factor matrices. Generally, a larger dimensionality of the factor matrix is set, and then low-rank regularization is integrated to suppress the effect of noise or shadow. Therefore, the parameters r w , r h , r a of factor matrices are also the dimensionalities of the core tensor, whose RMSE performance is shown in Figure 3. We can observe that the RMSE is larger when the values of each dimensionality are smaller. However, with the increase in the weight, the RMSE performance is gradually stable. For instance, the RMSE performance is stable when the values of spatial dimensionalities r w and r h are greater than 120, and the dimensionality of spectral factor matrix r a is no less than 10 for three datasets. Since the dimensionality of the spectral factor matrix represents the number of endmembers, it is robust in the mode-3 expression [18], [24]. In other words, as long as it is greater than the groundtruth, the fusion method shows a relatively stable performance, which can be interpreted as the shadow effect. Therefore, the dimensionalities of the core tensor in this article are set to r w = r h = 140, and r a = 15.

3) STOPPING RULES OF CONVERGENCE
As shown in Algorithm 1, the proposed LRTVS method decomposes the problem (7) into alternating iterations of three factor matrices and a core tensor. The stopping rule satisfies the relative difference threshold between the successive updates of the objective function F 0 (W, H, A, C) is less than 0.001 [44], [56], as shown in equation (43). The experimental results show that the increase in iterations has no significant effect on the convergence in the ADMM-based algorithms, so they need not run exhaustively. Therefore, the number of ADMM iterations is set to 20.

4) SNR SETTINGS
The proposed method promotes the vertical smoothing based on total variation for mode-3 factor matrix A. Thus, in order           difference images can only reflect the performance of partial bands, and the overall performance difference between the fusion algorithms can be distinguished by the tables and curves below.
The overall performance of all methods under three noise conditions for three datasets is shown in tables 1, 2 and 3. We can see that the performance of all algorithms declines with an increase in noise level (i.e. the reduction of SNRs). For the Pavia University and Moffett datasets, the perfor- mance of the proposed LRTVS is obviously better than the benchmark algorithms. Compared with CSTF, the performance RSNR of the LRTVS algorithm is improved by 2∼5dB, the RMSE decreases by more than 20%, and the performance of other metrics is also improved significantly. For the Washington DC dataset, except for the metric ERGAS, the performance of the proposed LRTVS is better than that of the benchmark methods. Especially, the proposed LRTVS performs significantly better than CSTF, although both of them are based on Tucker decomposition.
To test the performance of different bands, the curves of RSNR, RMSE and ERGAS are drawn in figures 7, 8 and 9. We can see from figures 7 and 8 that the performance of the proposed LRTVS is much better than that of the benchmark algorithms for Pavia University and Moffett datasets. Especially, compared with the CSTF, the LRTVS's performance has been greatly improved, and its curve is relatively smooth. In figure 9 of the Washington DC dataset, the proposed LRTVS has obvious advantages over other benchmark algorithms in the band range from 500 nm to 1000 nm, while the BSR and CNMF algorithms perform well in the other bands. In conclusion, the proposed LRTVS method in this article performs better than the baseline algorithms with different bands in most cases.
The target image is unfolded along mode 3 (spectral dimension) to obtain the factor matrix A, which represents the endmember signature, as does the NMF-based method. Using the gradual change in endmember signature curves, the total variation regularization is adopted to vertically smooth the spectral factor matrix, which can suppress the noise effect in the reconstructed images. Because both CSTF and LRTVS are based on Tucker decomposition, they are selected to verify the validity of the smooth regularization. For the Washington DC dataset, the spectral curves in the pixel with spatial coordinates of (50, 50) are shown in Figure 10. When the reflectance of the endmember is large among eight major endmembers, the spectral curves of the LRTVS and CSTF algorithms are almost identical, as displayed in sub-figures (a) and (b). Conversely, when the reflectance value is small, the spectral curve of the CSTF algorithm changes dramatically by the noise, while the curve of the LRTVS is relatively smooth, as shown in sub-figures (e) and (f). Thus, the curves illustrate that the total variation regularizer is effective to remove the noise.

D. COMPUTATIONAL COMPLEXITY ANALYSIS
The proposed LRTVS method involves the calculation of W, H, A, and C. Firstly, the Sylvester equations of factor matrices W, H and A are solved by using the conjugate gradient method, whose computational complexities are O(n 2 w N 2 w ), O(n 2 h N 2 h ) and O(n 2 a N 2 a ), respectively. Using the vector-matrix operator, they can be reduced to O(n 2 w N w +n w N 2 w ), O(n 2 h N h + n h N 2 h ) and O(n 2 a N a + n a N 2 a ). For the core tensor C, the complexities of equations (39a) and (39b) are both equal to O(n 3 w n 3 h n 3 a ), which can be simplified as O(n 2 w n h n a +n w n 2 h n a + n w n h n 2 a ) by Kronecker operators. In addition, the runtime of the algorithm is an important measure of complexity. The runtime of the proposed LRTVS and benchmark algorithms in this article is presented in table 4. We can observe that the runtime of NLSTF and CSTF algorithms is less, and the proposed LRTVS is still time-consuming, which needs to be further improved.

VI. CONCLUSION
Hyperspectral and multispectral data are three-dimensional data cubes, which need to be reshaped into matrices in the NMF-based fusion methods. To address the structure loss of hyperspectral and multispectral data cubes and suppress the effect of noise and redundant information, this article proposes a low-rank decomposition-based fusion method with total variation and sparse regularization to achieve the fusion of HSIs and MSIs for hyperspectral super-resolution images. Furthermore, the multi-variable fusion problem is decoupled into multiple single-variable subproblems via the AO and mode-n product, which is solved by the ADMM and conjugate gradient methods. In addition, tensor operators are used to simplify the closed-form solutions so as to reduce the computational complexity. The experimental results illustrate that the proposed LRTVS method can improve the performance of RSNR by 2∼5 dB and drop the RMSE by more than 22%. In conclusion, the proposed algorithm can not only improve the spatial resolution of the reconstructed image, but also greatly enhance the anti-noise aperformance, which fully verifies the validity of this regularization method.