extfieldsofstudy:
Engineering
Computer Science

text:
Multi-modal and frequency-weighted tensor nuclear norm for hyperspectral image denoising

Low-rankness is important in the hyperspectral image (HSI) denoising tasks. The tensor nuclear norm (TNN), defined based on the tensor singular value decomposition, is a state-of-the-art method to describe the low-rankness of HSI. However, TNN ignores some physical meanings of HSI in tackling denoising tasks, leading to suboptimal denoising performance. In this paper, we propose the multi-modal and frequency-weighted tensor nuclear norm (MFWTNN) and the non-convex MFWTNN for HSI denoising tasks. Firstly, we investigate the physical meaning of frequency components and reconsider their weights to improve the low-rank representation ability of TNN. Secondly, we consider the correlation among two spatial dimensions and the spectral dimension of HSI and combine the above improvements to TNN to propose MFWTNN. Thirdly, we use non-convex functions to approximate the rank function of the frequency tensor and propose the NonMFWTNN to relax the MFWTNN better. Besides, we adaptively choose bigger weights for slices mainly containing noise information and smaller weights for slices containing profile information. Finally, we develop the efficient alternating direction method of multiplier (ADMM) based algorithm to solve the proposed models, and the effectiveness of our models are substantiated in simulated and real HSI datasets.

I. INTRODUCTION
H YPERSPECTRAL image (HSI) has been widely used in many fields due to its wealthy spectral information of real scenes, e.g., target detection [1], unmixing [2], classification [3] and so on. Unfortunately, due to imaging conditions and the influence of weather, observed HSIs are often accompanied by the Gaussian noise, salt and pepper noise, and stripe noise. This not only reduces the quality of HSIs, but also hinders subsequent applications. Therefore, to improve the accuracy of subsequent applications, an important research topic is to recover the clear HSI from the observed HSI in the hyperspectral processing tasks [4,5].
Since the spectral bands are the imaging results of the same scene under different wavelengths, they have global correlation or low-rankness. Based on this prior, HSI can be unfolded into matrices along the spectral mode and then processed by matrix restoration methods. The low-rank matrix recoverybased method restores the clean HSI by minimizing the ranks of unfolding matrices, which is the most popular way [4,[6][7][8]. Zhang et al. [7] propose to disassemble the noisy HSI into the clean HSI and noise information and use the matrix nuclear norm as the relaxation of the rank function. To improve Sheng Liu, Xiaozhen Xie, and Wenfeng Kong are with College of Science, Northwest A&F University, Yangling 712100, China (e-mail: liu sheng@nwafu.edu.cn; e-mail: xiexzh@nwafu.edu.cn (Corresponding author); e-mail: kksama@nwafu.edu.cn. the low-rank performance of the model, some non-convex methods are introduced into the HSI denoising task. Xie et al. [9] use the log sum function, Chen et al. [10] use the γ function, Lou et al. [11] use the L 1 − L 2 function and Gao et al. [12] use the SCAD and MCP function as the nonconvex approximation of the rank function. Low-rank matrix decomposition is also an effective way to solve the low-rank minimization problem [13,14].
However, the unfolding operator destroys the highdimensional structure of HSI. Treating HSI as a third-order tensor can protect this structure, and tensor rank minimization is the most effective way to overcome the above shortcoming. Due to the non-unique definitions of the tensor rank, different tensor decompositions and the corresponding tensor ranks are proposed, e.g., the CANDECOMP/PARAFAC (CP) decomposition [15,16], Tucker decomposition [17,18] and tensor singular value decomposition (t-SVD) [19,20]. CP rank is defined as the minimum of rank-1 tensors by CP decomposition. Whereas, it is difficult to estimate when HSI is unknown. By Tucker decomposition, a tensor can be decomposed into core tensor and modal matrices. Then, Tucker rank is approximated by the sum of modal matrix nuclear norms (SNN) [18]. However, SNN has shortcomings in preserving low-rank structure, and its results are materially sub-optimal [21,22].
Based on t-SVD, the tensor tubal rank has attracted extensive attention recently [23][24][25][26]. Its convex relaxation is the tensor nuclear norm (TNN), which is defined as the sum of the matrix nuclear norms of each frontal slice in the Fourier transformed tensor [19,21]. The minimization of TNN can be quickly solved by convex optimization algorithms, and it is effective in characterizing the low-rankness of tensor. Therefore, it is widely used in HSI denoising task [24,27]. In the definition process of TNN, the information in the image domain of HSI is transformed into frequency information in the Fourier domain. In the Fourier transformed tensor, the lowfrequency slices mainly carry the profile information, while the high-frequency slices mainly carry detail and noise information. Besides, in each frontal slice of the Fourier transformed tensor, bigger singular values mainly contain information of clean data and smaller singular values mainly contain information of noise. Meanwhile, HSIs have correlations among different modes, but TNN lacks the flexibility to explore different correlations along with different modes [25]. Thus, it is advantageous to improve the flexibility of TNN when considering the prior information about the physical meaning of frequency, correlation among different modes, and singular value.
In this paper, to take full advantage of the above prior information and improve the capability and flexibility of model, we propose the multi-modal and frequency-weighted tensor nuclear norm (MFWTNN) and non-convex MFWTNN. In the proposed MFWTNN, we explore the physical meaning of frequency to explore the low-rankness of HSI in the Fourier domain and give an adaptive calculation method of frequency weights. In addition, we consider the correlation among the spatial modes and the spectral mode. Furthermore, based on MFWTNN, we take into account the physical meaning of singular values inside frequency slices and propose the nonconvex approximation of their nuclear norms. Finally, we apply them to the HSI denoising task. The main contributions of this paper are summarized as follows: • Motivated by Tucker decomposition, we consider the correlation among spectral and spatial modes of HSI and use the mode-p permutation of tensor instead of the mode-p unfolding matrix. And then we use their weighted sum to explore the low-rankness of HSI. Furthermore, we combine the advantages of multi modes with the improvement of TNN in the frequency domain and propose MFWTNN in HSI denoising. The MFWTNN is the fusion of TNN and SNN, which inherits the advantages of t-SVD and Tucker decomposition. • In each frontal slice of the Fourier transformed tensor, we choose the log sum function as a non-convex approximation to the rank function, which can provide a better approximation to rank function compared with nuclear norms in MFWTNN. Different from the onedimensional weights of the classical non-convex models, NonMFWTNN considers the weight of frequency direction and singular value direction, which is the twodimensional weight. • According to information types in different frequency slices in the Fourier transformed tensor, we adaptively choose bigger weights for slices mainly containing noise information and smaller weights for slices containing profile information, which can depress noise more and simultaneously preserve the profile information of clean HSI better. • We develop the efficient alternating direction method of multiplier based algorithm to solve the proposed models, and obtain the best restorative performance both on the simulated and real HSI dataset in comparison to all competing HSI denoised methods.

A. Notations
In this section, we follow the work in [28,29] and give some basic notations used about matrix and tensor. We use bold upper-case letter X and calligraphic letter X to denote matrix and tensor, respectively. For a 3rd-order tensor X ∈ R n1×n2×n3 , its (i, j, k)-th component is denoted as X (i, j, k). For X , Y ∈ R n1×n2×n3 , their inner product is defined as n3 k=1 x ijk y ijk . Then the Frobenius norm of a tensor X is computed as n3 k=1 |x ijk |. The k-th frontal slice of X is represented as X (k) = X (:, :, k). For one-dimensional vector v, its fast Fourier transform is denoted asv = fft(v) ∈ C N , and inverse transformation is represented as v = ifft(v). For three-dimensional tensor X , its fast Fourier transform is represented asX = fft(X , [], 3) and its inverse operation is X = ifft(X , [], 3). The mode-p permutation of X is defined as X p = permute(X , p), p = 1, 2, 3, where the mth mode-3 slice of X p is the m-th mode-p slice of X , i.e., X (i, j, k) = X 1 (j, k, i) = X 2 (k, i, j) = X 3 (i, j, k). Also, its inverse operation is X = ipermute(X p , p).

B. Problem Formulation
A clean HSI can be treated as a third-order tensor X ∈ R n1×n2×n3 and is usually assumed to be low-rank. Corrupted by mixed noise, its observed version can be modeled as where Y, S, N ∈ R n1×n2×n3 ; S denotes the sparse noise; N denotes the Gaussian white noise. HSI denoising aims to recover the clean HSI X from the observed HSI Y in (1). Under the framework of regularization theory, it can briefly be formulated as arg min where · 1 describes the sparse noise; · F describes the Gaussian noise; Rank(·) represents the rank of unknown ideal HSI; λ and τ are non-negative parameters.

C. Tensor Nuclear Norm
Since the tensor rank function is non-convex and discrete, its rank minimization problem is NP-hard [30]. To overcome computational difficulties, a usual approach is to approximate the non-convex rank function to the convex function and then substitute it into the minimization problem. Fazel [31] proves that the matrix nuclear norm is the convex envelope rank function on the unit sphere of the spectral norm. Thereafter, the nuclear norm is used as a convex replacement of the rank function. The rank minimization problem is transformed into nuclear norm minimization, and it is widely used in low-rank matrix recover based HSI denoising tasks [7].
Furthermore, HSI can be regarded as a third-order tensor. Then, based on t-SVD [19], Lu et al. [21] propose the tensor nuclear norm (TNN) as the relaxation function of the tensor rank function. TNN is defined as whereX (k) is the k th frontal slice ofX . Due to the property of fft, we have Based on (4), frequency component [26] is defined as if i = 2, · · · , n3+1 2 (5) Fig. 1 shows the relationship between the frontal slice and the frequency component. Due to TNN can be effective in expressing the low-rank characteristics, it has received extensive attention in the HSI denoising task. However, it focuses on the low-rankness in the spectral dimension. For HSI, its low-rank structure exists with multiple dimensions [29,32]. To connect the correlation along different dimensions, Liu et al. [18] use modal matrices of the tensor in different directions to explore low-rank structures, and then use SNN as the convex relaxation of the rank function. SNN is defined as where X p is the modal matrix that unfold X along the direction of the mode-p, α p is weighting parameter of modes and 3 p=1 α p = 1. However, this way will inevitably destroy the internal structure of HSI. Zheng et al. [25] propose the weighted sum of TNN (WSTNN). They consider the lowrankness of the three modes, which are two spatial modes and one spectral mode. WSTNN is defined as where X p is the mode-p permutation of X .
III. MULTI-MODAL AND FREQUENCY-WEIGHTED TENSOR NUCLEAR NORM AND ITS NON-CONVEX APPROXIMATION.

A. Multi-modal and Frequency-weighted Tensor Nuclear Norm
In the definition of TNN (3), fft plays a core important role in the process of transforming t-products into matrix multiplication in the Fourier domain. It makes TNN to have fast calculation capabilities. However, the role of fft is not just for fast computation in the framework of t-SVD. When Fourier matrix F n is treated as a transform matrix,X is the feature tensor of X . To improve the low-rank representation ability of TNN, Wang et al. [26] propose the frequency-filtered tensor nuclear norm (FTNN), which is defined as is the sum of two matrix nuclear norms, I denotes the number of frequency bands, and α k is a pre-defined parameter assigned to the k-th frequency band, whose value depends on the prior knowledge. They treat the FTNN as a frequency filtering method. α k characterizes the filtering coefficient for the kth frequency component. However, these parameters should be data-dependent, and the pre-weight parameters cannot well reflect the real physical meaning between filter coefficient and frequency component. For the HSI, it is difficult to give n3+1 2 pre-weight parameters in advance. Therefore, we reconsider these weight parameters based on the physical properties of HSI. For the clean Washington DC Mall dataset, sparse noise is added to form clean data, and then the zerofrequency component and non-zero frequency component of the noise data are separated. Fig. 2 shows the grayscale image of the 60th band. We can find an interesting phenomenon that just using the information of the zero frequency component can get satisfactory results without any processing and the noise is more concentrated in the non-zero frequency component. As the transform operator, fft plays the role of information separation, which makes the image information and noise information as separate as possible. The distribution of the nuclear norm shown in Fig. 2(e) is consistent with the observed phenomenon. When the HSI is polluted by noise, the variation in the high-frequency component is bigger than that in the low-frequency component. There is a natural idea that we can reduce the shrinkage of the low-frequency component and increase the penalty of the high-frequency component. Therefore, base on FTNN, we redefine the frequency-weighted tensor nuclear norm (FWTNN): Our weights are different from the weights of FTNN, which are not pre-weight but data dependent. They are discussed in Section III-C.
In the definition of FWTNN, by calculating the nuclear norm along the spectral direction, fft can capture low rank of the HSI spectrum. However, the low-rankness of HSI exists not only in the spectral dimension, but also in the spatial dimensions [18,32]. Therefore, fft is insufficient in exploring intra-mode and inter-mode correlations of HSI. Besides, in the HSI denoising task, the HSI band may be polluted by stripe noise in the same direction. In a low-rank model that only considers the spectral dimension, these stripe noises will be regarded as part of the low-rank image and are hard to remove. In the frequency domain, for the same HSI, the information expressed by the permuted versions X p in different directions is even more different. Therefore, motivated by SNN [18], it is necessary to consider the low-rankness of different dimensions for HSI denoising task. Different from SNN, to better retain the high-dimensional structures of HSI, its modep permutation X p to replace its mode-p unfolding matrix X p .
Based on the above analysis, we propose the multi-modal and frequency-weighted tensor nuclear norm (MFWTNN) as follows: is the k-th frontal slice of X p and its assigned weight is w p k ; α p > 0 and 3 p=1 α p = 1.

B. Non-convex Multi-modal and Frequency-weighted Tensor Nuclear Norm
In the same frequency band, its singular values are treated equally. However, the major information of HSI, such as smooth zones and profile, is contained in the larger singular value; the noise information of HSI is contained in the smaller singular value [9][10][11]33]. To combine both the frequency and the singular value prior information, we use the log norm to perform non-convex relaxation of the nuclear norm in the frequency band to more accurately describe the tensor rank function [9]. Thus, we propose the non-convex multi-modal and frequency-weighted tensor nuclear norm (NonMFWTNN) to more accurately describe the low-rankness of HSI. It is defined as follows: where log Fig. 3 shows the mechanism of Non-MFWTNN. As can be seen from Fig. 3, NonMFWTNN takes into account the relationship between the singular values inside and outside of frequency band. The weights of classic matrix / tensor non-convex models [9,10,29,34] are one-dimensional, which only consider the relationship between the singular values inside the matrix. The way is equivalent to pulling the core singular value matrix S p in Fig. 3(c) into a column vector and using the weights as the coefficient of this column vector. The weight of MFWTNN is to make a Cartesian product of the weight vector W F in the frequency direction and the weight vector W S in the singular value direction, which are twodimensional weights. If W F = 1, it degenerates to classical one-dimensional weights.

C. Frequency-weighted parameter
In the definition of MFWTNN and NonMFWTNN, the values of the frequency-weighted parameter w p k are important. The ideal w p k can take smaller values that reduce the penalty for the low-frequency slice nuclear norm to protect the main information of HSI in the low-frequency band, and take larger values that increase the penalty for the high-frequency slice nuclear norm to fully remove the noise in the high-frequency band. It can be seen from Fig. 4(a) that the nuclear norm is opposite to the distribution characteristics of w p k in the frequency domain. It takes larger values in the low-frequency band and takes smaller values in the high-frequency band. And the value of the nuclear norm decreases as the frequency increases. However, the calculation of the nuclear norm brings high cost, which requires singular value decomposition. In a finite dimensional space, any two norms are equivalent. This shows that when transforming from the nuclear norm to the Frobenius norm, they are equal under the normalization condition. Thus, the reciprocal of the Frobenius norm of each frequency slice matrix can be used as w p k .
where δ = 10 −6 . The purpose of using log function is to shrink the scale of the value without changing the trend of the Frobenius norm to prevent the extreme influence of the maximum value. The above method of calculating weights is feasible when the clean HSI is known. However, for HSI denoising tasks, we often only have noisy HSI, and clear HSI is unknown. It can be seen from Fig. 4(b) that the weights calculated by the denoised HSI is closer to the result obtained by the clean HSI than the noisy HSI. Further, as shown in Fig. 3(e) , as the number of iterations n continues to increase, the weights calculated by the denoised HSI will get closer to the result of clean HSI. Based on this, we propose an iterative calculation method to update w p k . Assuming that the nth iteration getsX where h p k,n = 1/(log X (k,n) p 2 F + δ). C 1 is the scaling factor after frequency normalization; C 2 is a smaller constant. Affected by noise, there is less useful information from the denoised HSI in the initial iteration. If it is not corrected, the denoised result will develop in a bad direction. Since TNN is a robust model, inspired by this, we introduce a constant C 2 to improve our weights. When C 1 =0, it degenerates into TNN.
Another reason we introduced C 2 is that noise exists in all frequency bands. The previous conclusion is that there is less noise information in the low-frequency part instead of none. C 2 can achieve a slight shrinkage of all frequency bands.

IV. HSI DENOISING VIA MFWTNN AND NONMFWTNN MINIMIZATION
A. Proposed model MFWTNN uses frequency components, modal information, and NonMFWTN uses the physical meaning of singular value distribution on the basis of MFWTNN. They can provide a better approximation to the tensor rank. Then we use MFWTNN and NonMFWTNN to replace the regularization term Rank in (2) and propose the HSI denoising model as follows: arg min arg min Introducing auxiliary variables, model (14) and (15) are equivalent to arg min arg min By augmented Lagrangian multiplier method, the Lagrangian function of model (16) can be written as

B. Algorithmic optimization
ADMM is an effective framework to solve these type of minimization problems [35]. When other variables are fixed at the n-th iteration, each variable in the Lagrangian function can be updated by solving its corresponding subproblem respectively at the (n + 1)-th iteration.
For Z p , p = 1, 2, 3, their corresponding subproblems can be written as The closed-form solution of (19) obtained from theorem 1 of [26] are as follows: Similarly, when (17), the sub-problem of Z p is rewritten as The closed-form solution of (21) obtained from theorem 1 are as follows: For X , its corresponding subproblem can be reformulated as X can be updated as follows: For S, its corresponding subproblem can be reformulated as It can be solved by the soft-thresholding operator [36] as: N can be updated as follows : For multipliers Γ p and Λ, they can be updated as follows:

C. Time Complexity Analysis
The HSI denoising models based on MFWTNN and NonMFWTNN can be solved by Algorithm 1. Further, we discuss the time complexity of this algorithm. Computing Z p , p = 1, 2, 3, in both MFWTNN-based solver and NonMFWTNN based solver, has a complexity of O(n 1 n 2 n 3 log(n 1 n 2 n 3 ) + 3 i=1 max(n i , n i+1 ) min 2 (n i , n i+1 )n i+2 ), where n 4 = n 1 , n 5 = n 2 . Calculating w k has a complexity of O( 3 i=1 min (n i , n i+1 ) n i+2 ). The cost of computing X , S, and N are all O(n 1 n 2 n 3 ). Then, by calculating the complexity of the above variables, the total complexity of the proposed algorithm can be obtained as O(n 1 n 2 n 3 log (n 1 n 2 n 3 ) + 3 i=1 max (n i , n i+1 ) min 2 (n i , n i+1 )n i+2 + 3 i=1 min (n i , n i+1 )n i+2 + 3n 1 n 2 n 3 ).

V. EXPERIMENT RESULTS AND DISCUSSION
In this section, we conduct experiments on simulated and real-world HSIs to substantiate the effectiveness of the proposed MFWTNN and NonMFWTNN model for HSI denoising. For evaluating the effectiveness of the models, our models are compared with four state-of-the-art HSI denoised methods i.e.BM4D [37], LRMR [7], LRTDTV [32], 3DLogTNN [29]. Due to the BM4D methods that are only suitable for removing Gaussian noise, we first use the RPCA [21] model to remove the sparse noise, and then use it to obtain the denoised results. Before denoising all models, all pixels in the HSI band are linearly mapped to [0,1]. In all experiments, the hyper-parameters in these comparative methods are manually adjusted according to the suggestions in relevant papers. In Section V-C, the parameters are discussed in detail in our proposed models. We use visual comparison and quantitative comparison to comprehensively evaluate the performance of different denoising methods. In quantitative comparisons, we use five quantitative picture quality indices (PQIs), i.e., PSNR [38], SSIM [39], FSIM [40], ERGAS [41] and SAM [42], to assess the denoised results. The better denoised results correspond to higher values in MPSNR, MSSIM, MFSIM, and lower values in EGRAS, MSAM.

A. Simulated HSI Data Experiments
In this section, we select two HSI datasets for the simulated experiments to show the robustness of our methods. One is filmed by the reflection optical system imaging spectrometer (ROSIS-03) from Pavia City Center dataset 1 . The size is 1096×1096, with a total of 102 bands. The reason we chose a sub-block with a size of 200×200×80 as a simulation dataset is because some bands are seriously polluted by noise and have lost the meaning of reference. The other one is collected by a HSI acquisition sensor (HYDICE) from Washington DC Mall dataset 2 . The size is 1208 × 307, with a total of 191 bands. we choose a sub-block with a size of 256 × 256 × 191 as a simulation dataset.
The observed HSIs are degraded by the hybrid noises, i.e. Gaussian noise, sparse noise, and stripe noise. Thus, to better simulate the degradation mechanism of HSI, we add various intensity of hybrid noises to the clean HSI. In the simulation experiments, the hybrid noises with eight different intensity levels are added to simulation dataset band by band. Let G and P denote the variance of Gaussian white noise and percentage of impulse noise, respectively. In noise cases 1-5, the same intensity noise is added to all the bands. In noise case 1, G=0.1 and P =0.2; In noise case 2, G=0.1 and P =0.3; In noise case 3, G=0.1 and P =0.4; In noise case 4, G=0.15 and P =0.2; In noise case 5, G=0.2 and P =0.2; In noise cases 6-8, the noise intensities are different for different bands. In noise case 6, G=0.1 and P is randomly selected from 0.2 to 0.4. In noise case 7, G is randomly selected from 0.1 to 0.3 and P =0.2; In noise case 8, G is randomly selected from 0.1 to 0.3 and P is randomly selected from 0.1 to 0.3. In addition, stripe noises are added to the 54-64 bands of the Pavia City Center and the 70-100 bands of the Washington DC Mall. In this part, we show the evaluation results of all denoising models including visual and quantitative quality on Pavia City Center and Washington DC Mall datasets. Table I reports the denoised results and CPU running time by all the compared methods under the condition of the eight noise cases on simulated dataset. When compared to the other methods, our proposed methods obtain the optimal PQIs among all the denoising models in most noise cases, indicating the advantage of the proposed methods in HSI denoising. Although LRTDTV also consider the low-rankness between different modes, it is affected by the matrix unfold operation, which destroy the connection between the modes. Therefore, LRTDTV models obtain the lower PQIs. BM4D and LRTDTV use non-local similarity information and spatial smoothing information respectively. These make the denoised result excessively smooth and lack texture details. Although 3DlogTNN considers the correlation between different modes, it is still not accurate enough to represent the low-rankness of HSI. MFWTNN use the physical meaning of singular values and frequency components to improve TNN respectively. But it is not the best descriptions of the low-rankness of HSI. NonMFWTNN inherits the advantages of 3DLogTNN and MFWTNN, and at the same time improves the ability to express low-rankness of HSI in terms of frequency components and singular values. Although the time complexity of our models is higher, our models have obtained better denoised results.
In Fig. 5, we show the denoised results of the Pavia City Center and Washington DC Mall dataset in terms of PSNR and SSIM under all noisy case experiments. As this figure shows, our methods obtain the optimal PSNR and SSIM values in most bands. In Fig. 6, we show the 50th band of the denoised results of all comparison methods in case 2 on Pavia City Center dataset. Also in Fig. 7, We show the 57th band of the denoised results in case 7 on Washington DC Mall dataset. As these grayscale images show, LRMR cannot completely remove some high-intensity noises and retain more noise. For BM4D and LRTDTV, although they can remove more noise, they also lose more details. This makes the denoised result too smooth. Compared with them, our proposed models can remove more noise while retaining more details. The spectral signatures are important indicators for HSI. To further compare the denoised image quality, we show the original and denoised spectral signatures curves. Therefore, we show the spectral curve of the pixel (170, 99) under case 1 and pixel (234,171) under case 6 in Fig. 8. It can be seen that spectral signatures curves are relative to ERGAS values in table I, and our models can obtain the optimal value.

B. Real HSI Data Experiments
In this subsection, we chose the real HSI dataset in this experiment. It is AVIRIS Indian Piness dataset 3 . It is filmed by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS). The reason we chose a sub-block with a size of 145×145×220 3 https://engineering.purdue.edu/ as a real dataset is because some bands are seriously polluted by Gaussian white noise and impulse noise and have lost the meaning of reference. We evaluate our denoising model from visual evaluation and vertical average profile. In Fig.  9, we show the 105th band of the observed and denoised results in Indian dataset. It can be seen that it is seriously polluted by Gaussian and sparse noise, and the original ground features information cannot be distinguished. Although the competitive denoising methods restore the feature information of HSI, the denoised HSIs still lose detailed information or residual noise, resulting in missing or blurred texture details of HSI. Our models not only restore more feature information of HSI but also remove more local noise than the competitive methods. Therefore, these denoised HSIs have clearer texture information. In Fig. 10, we show the vertical mean profiles of the 50th band of Indian dataset. Due to the influence of hybrid noises, the vertical mean profiles of the observed HSI shows rapid fluctuations. This means that the smaller fluctuations in the vertical mean profiles, the higher quality of the denoised HSI. Obviously, the mean profile curves obtained by our model are the most stable and its fluctuation is the smallest. This result is relative to the facts reflected in Fig. 9.
C. Discussion 1) Parameter Analysis: In the NonMFWTNN model, five parameters need to be discussed, i.e., the modal weight parameter α, the regularization parameters λ, τ and frequency weight parameter C 1 , and C 2 . All experiments are done under Case 2.
The modal weight α is to control the ratio between different modes. The low-rank characteristics of the spatial dimensions are similar. To reduce the parameters, we set the weights of the two spatial dimensions to equal weights. Therefore, it is set to α = [1, 1, α 3 ]/(1 + 1 + α 3 ), which makes 3 i=1 α i = 1. We study the impact of α 3 which is selected from {0.1 − 1.0} on the robustness of the model, and its optimal value is obtained at 0.2. It can be observed from Fig. 11(a) that when α 3 takes a smaller value, the PSNR value of the denoised result is better. This is because the low-rank characteristic of the spectral mode is stronger than the spatial mode. The regularization parameter τ controls the weight of the Gaussian noise term. Its value is related to the intensity of Gaussian noise. Therefore, it is set to τ = τ n /σ, where σ is Gaussian noise intensity. We study the impact of τ n which is selected from {10 −5 − 1.6 × 10 −4 } on the robustness of the model, and its optimal value is obtained at 10 −4 . It can be observed from Fig. 11(c) that when τ ≤ 0.0001, the MPSNR of the denoised result is low and fluctuates greatly. The main reason is that when τ is small, the ability to remove Gaussian noise is weak, and there is a lot of remaining Gaussian noise. The regularization parameter λ controls the weight of the sparse noise term. Defining λ = λ s ( α1 √ max(n2,n3)n1 ).
As observed from Fig. 11(b), our model can obtain remarkable results when λ s which is selected from {0.006 − 0.016}. its optimal value is obtained at 0.011. Similar to τ , when λ is small, the noise removal effect is not good. The frequency weight parameter C 1 reflects the scaling ratio of the weight. Fig. 11(d) show the PSNR values when C 1 which is selected from {0.1 − 1.0}. It can be seen from these figures that when C 1 is small, the frequency weight becomes smaller, and the noise is more difficult to remove; when C 1 is large, the frequency weight becomes larger, and the detailed information of tensor is removed. its optimal value is obtained at 0.6. The frequency weight parameter C 2 reflects the initial weight of the iteration. Fig. 11(e) show the PSNR values when C 2 which is selected from {0.1 − 1.0}. its optimal value is obtained at 0.6. It shrinks equally of all frequency components. It inherits the advantages of TNN and ensures that noise can be removed correctly.
2) Convergence Analysis: Fig.12 shows the MPSNR value, the MSSIM value, and the Error value of the proposed model according to the iterations. As shown in Fig.12, the number of iterations increases, the value gradually stabilizes by the proposed method, which justifies the numerical convergence of the proposed method.

VI. CONCLUSION
In this paper, we proposed the multi-modal and frequencyweighted tensor nuclear norm (MFWTNN) and Non-convex MFWTNN to approximate the rank function and applied them to the HSI denoising task. The MFWTNN not only considered the correlation between two spatial modes and one spectral mode but also considered the frequency meaning of HSI in the Fourier domain. The MFWTNN is the fusion of TNN and SNN, which inherits the advantages of t-SVD and Tucker decomposition. And we gave the frequency weight adaptive calculation way. Since the calculation of weights depends on iteratively restored and observed HSI, our models are more robust to different degrees of noise. It can better characterize the low-rankness of HSI. Besides, based on the MFWTNN model, we also considered the physical meaning of the internal singular values of the frequency slice and proposed its nonconvex approximation. They powerfully improved capability and flexibility for describing low-rankness in HSIs. The ex-  periments conducted with both simulated and real HSI datasets showed that our proposed models based HSI denoising model were competitive methods to remove the hybrid noises.


ASAT: Adaptively Scaled Adversarial Training in Time Series

Adversarial training is a method for enhancing neural networks to improve the robustness against adversarial examples. Besides the security concerns of potential adversarial examples, adversarial training can also improve the generalization ability of neural networks, train robust neural networks, and provide interpretability for neural networks. In this work, we introduce adversarial training in time series analysis to enhance the neural networks for better generalization ability by taking the finance field as an example. Rethinking existing research on adversarial training, we propose the adaptively scaled adversarial training (ASAT) in time series analysis, by rescaling data at different time slots with adaptive scales. Experimental results show that the proposed ASAT can improve both the generalization ability and the adversarial robustness of neural networks compared to the baselines. Compared to the traditional adversarial training algorithm, ASAT can achieve better generalization ability and similar adversarial robustness.

Adversarial examples are generated by adding small malicious perturbations to the input data that can mislead models, which reveal the vulnerability of neural networks with respect to the input data. Adversarial training [2,5,6,7,8] is a method for enhancing neural networks that can improve both the adversarial robustness and the generalization ability at the same time. Besides the security concerns that neural network faces potential risks of adversarial example attack, researches on adversarial training focus on improving the accuracy of the neural networks [6,7,8], training robust neural networks [9,10,11] Our contributions are summarized as follows: • Rethinking existing adversarial training approaches, we propose to rescale dimensions of the perturbation according to their importance and adopt timedependent or adaptive scales with multi-step risk averaging. We propose an adaptively scaled adversarial training (ASAT) algorithm.
• We adopt ASAT on various baselines, including our newly proposed Transformer baselines. Results show that ASAT can improve both the generalization ability and the adversarial robustness of neural networks compared to the baselines. Compared to the traditional adversarial training algorithm, ASAT can achieve better generalization ability and similar adversarial robustness.
• With the proposed dimension-wise adversarial sensitivity indicator, we probe the sensitivities of different input dimensions and explain the decision bases of black-box neural networks, which provides more interpretability for neural network decisions.

Adversarial Training
In this section, we first introduce the concept of adversarial examples and adversarial training. Then, we describe the proposed dimension-wise adversarial sensitivity indicator inspired by the adversarial sensitivity.

Adversarial Examples
Let D denote the dataset, x = (x 1 , x 2 , · · · , x k ) T ∈ R k and y stand for a data input and its ground truth, f denote the neural network with the parameter vector θ and L denote the loss function. The adversarial example [1,2] can be defined as the small perturbation δ = (δ 1 , δ 2 , · · · , δ k ) T ∈ R k on the data input x that can mislead the model and cause the maximum increases in the loss function: where S denotes the constraint set, the perturbation is specified by S. For example, S = {δ : δ 2 ≤ } specifies that the perturbation should be in a k-dimensional hyperball with a radius of . Typical constraints are L p -norm bounded (1 ≤ p ≤ +∞) and S = {δ : δ p ≤ }, where the L p -norm is defined as: Since neural networks are sensitive to malicious perturbations in inputs, i.e., adversarial attacks, the adversarial risk is defined as the loss after adversarial attacks and the empirical risk is defined as the initial loss: The loss increases caused by adversarial attacks can be defined as the adversarial sensitivity: The adversarial risk or sensitivity can be defined on a dataset or a single data instance. The adversarial risk or sensitivity on the test dataset can evaluate the robustness of neural networks with respect to adversarial attacks, a lower adversarial sensitivity or risk means a better adversarial robustness.

Generating Adversarial Attacks
To generate adversarial attacks, the fast gradient methods [2], including the fast gradient sign method (FGSM) and the fast gradient method (FGM), propose to optimize Eq.(1) via optimizing the first-order Taylor Expansion term of the loss function instead when S = {δ : δ p ≤ }, namely maximizing the inner product of δ and the gradient g = ∇ x L(f (x, θ), y): where denotes element-wise product and sgn denotes element-wise sign function. FGSM adopts p = +∞ and FGM adopts p = 2.
Another representative line of existing studies for generating adversarial examples is PGD-based attacks [28,29]. Suppose S = {δ : δ p ≤ }. The projected gradient descent (PGD) algorithm optimizes Eq.(1) via a multi-step variant of fast gradient methods [2]: where δ 0 is usually set to the zero vector and u k is the k-th update in PGD. It is usually solved in Eq. (7) and {u : u p ≤ τ } and the step size is set to τ . Π S (v) is the projection function that projects v into the set S. In two common cases, L 2 -norm [8] and L +∞ -norm [28], the projection functions are: where the clip(h) function clips each dimension

Traditional Adversarial Training
Since neural networks are sensitive to adversarial examples, adversarial training [2,5] algorithms are designed to improve both the adversarial robustness and the generalization ability of neural networks. The target of traditional adversarial training is seeking to find the parameter vector θ with lower adversarial risks: [2] estimate the risk for adversarial attacks via the fast gradient methods. Other existing studies also estimate the risk via PGD-based attacks [5,6,7,8,29]. Fast gradient methods, including FGSM and FGM, can be seen as PGDbased attacks with the step number K = 1. Besides, some PGD-based adversarial training algorithms [7,8] minimize the adversarial sensitivities of multiple perturbations generated in different steps.

Dimension-wise Adversarial Sensitivity
Inspired by the adversarial sensitivity, we propose to probe the sensitivities of different input dimensions. Besides the adversarial sensitivity of the whole model, we can also define the dimension-wise adversarial sensitivity indicator, namely the sensitivity of the i-th input dimen- adv ( ) denotes the sensitivities of the i-th dimension, namely the expectation of the maximum loss change when the i-th dimension of the input changes no more than , which is defined as: The sensitivities of input dimensions reveal how the perturbations on different dimensions can cause the loss to increase. Neural networks pay more attention to dimensions with higher sensitivities and these dimensions tend to be relatively important dimensions for models to make decisions. Therefore, dimension-wise sensitivities can help reveal the decision basis of the model f and the importance of different dimensions in the prediction process. In traditional factor analysis methods, the factor loading [30] measures the importance of a factor, which can be seen as , (i = 1, 2, · · · , k), (13) namely R where g = ∇ x L(f (x, θ), y) and its i-th dimension is g i .
Namely R

Methodology
In this section, we first rethink the constraint set in adversarial training. Then, we propose to rescale different dimensions of the perturbation in the constraint set.
Last, we propose the adaptively scaled adversarial training (ASAT) algorithm.

Rescaling the Constraint in Adversarial Training
In this section, we rethink the constraint in adversarial training in time series and propose to rescale the constraint. 1 Note that L is only required to be convex and M -smooth in the neighborhood of θ instead of the entire R k .

Rethinking the Constraint in Adversarial Training
Adversarial training algorithms are widely adopted in the CV [1,2,5,6] and NLP [10,8,11] fields. The widelyadopted L p -norm constraint, S = {δ : δ p ≤ }, is invariant about different dimensions of the perturbation, which is reasonable in both the CV and NLP fields.
In the CV field, an important hypothesis of convolution neural networks is shift-invariant [25,26] and different pixels in an image can be treated of equal significance. In the NLP field, adversarial attacks can be conducted on the word embeddings and different dimensions of word embeddings are of similar significance [27,8].
However, in time series, different dimensions of inputs are not symmetrical. Therefore, directly adopting the L pnorm constraint in adversarial training in time series is not reasonable.

Rescaling the Constraint in Adversarial Training
Scales of all dimensions of the L p -norm bounded constraint are the same. However, in time series, dimensions with different timestamps may be of different significance.
The time series is usually defined as X = (v (1) , v (2) , · · · , v (T ) ), where v (t) ∈ R V is the input vector at the t-th timestamp and a smaller t means an earlier timestamp. For example, is the input vector at the first timestamp 9:00, Jan. 21, and v (2) is the input vector at the second timestamp 10:00, Jan. 21. In our work, we flatten a time series j is from the input vector with the t-th timestamp. Suppose x i is from the input vector with the t i -th timestamp, namely t i = i V , where · is the ceil function.
Suppose α = (α 1 , α 2 , · · · , α k ) T , we multiply the scale of dimension i of the perturbation with α i to rescale the radius of dimension i from to α i . As illustrated in Fig. 1, the L 2 -norm bounded constraint specifies a k-dimensional hyperball in R k (a circle in R 2 , as shown in Fig. 1a) and the L +∞ -norm bounded constraint specifies a k-dimensional hypercube in R k (a square in R 2 , as shown in Fig. 1b).
After rescaling, the hyperball is transformed into a hyperellipsoid (an ellipse in R 2 , as shown in Fig. 1c) and the hypercube is transformed into a hypercuboid (a rectangle in R 2 , as shown in Fig. 1d). The rescaling process assigns different importance to different dimensions of perturbations in adversarial training. Generally, the L p -norm constraint δ p ≤ is transformed to the constraint:

Fast Gradient Methods and PGD Algorithms after
Rescaling After rescaling, the constraint set is S = {δ : α −1 δ p ≤ }. For the fast gradient method target, namely maximizing the inner product of δ and the gradient g k−1 , the solution Eq.(7) becomes: For the projection function Π S (v) that projects v into the set S, we adopt the variants of Eq.(9) and Eq.(10): The proof and theoretical details are in Appendix.

Time-dependently Scaled and Adaptively Scaled Adversarial Training
In this section, we propose two lines of adversarial training with time-dependent or adaptive scales.

Multi-step Risk Averaging Adversarial Training
To

Time-dependently Scaled Adversarial Training
Adaptive scales α i are determined by the t i , we require α i ∈ [0, 1]. Since in time series analysis, data of closer timestamps with larger t i may be more important, We may assume α i ≤ α j when t i < t j , and max 1≤i≤k α i = 1. The scale α i decays when t i decreases.

Algorithm 1 Multi-step Risk Averaging Adversarial Training Algorithm with Dimension-wise Scales
Require: Network f ; parameters θ ∈ R k ; loss L and dataset D; constraint set S; step number K and step size τ ; optimizer O; iterations; batch size |B|.

7:
Calculate time-dependent or adaptive scales α i .
where the hyperparameter γ controls the decaying speed, and T denotes the number of different timestamps. Traditional adversarial training can be seen as adversarial training with constant scales.

Adaptively Scaled Adversarial Training
Time-dependent scales α i are determined by the i-th dimension in the data instance x i and its gradient g i = Similarly to time-dependent scales, we may assume that 0 ≤ α i ≤ 1 and max The ideal attack scales should relate to the scales of data instance, which will be analyzed later in detail in Theorem 2 and Sec. 3.3.2 that the ideal attack scales β should be proportional to the shifts between the training and test data. We may assume that the ideal attack scales or the data shifts are proportional to the scale of data instance |x i |, namely: As analyzed in Proposition 2, when conducting adaptive scales, the dimension-wise adversarial sensitivities are combined with Eq.(20), we may assume α i ∝ |xi| |gi| , namely In our implementation, we adopt: where we add a small positive real number e = 10 −6 on |x i | and |g i | for numerical stability.
As shown in Algorithm 1, we calculate adaptive scales

Theoretical Analysis
In this section, we first analyze the adversarial sensitivity under our proposed adaptively scaled adversarial attacks. Then, we reveal that the adversarial sensitivity under our proposed adaptively scaled adversarial attacks determines the generalization bound when the training and test distributions are not IID.

Adversarial Sensitivity after Rescaling
We show the adversarial sensitivity or risk under our proposed adaptively scaled adversarial attacks after rescaling in Theorem 1.

Theorem 1. Consider a neural network that is convex
and M -smooth in the neighborhood of θ, suppose g denotes the gradient, g = ∇ x L(f (x, θ), y), then its adversarial sensitivity or risk under our proposed adaptively scaled adversarial attacks where S = {δ : α −1 δ p ≤ } is:

Generalization Bound and Adversarial Sensitivity
Let z = (x, y), we may assume the training distribution P := p(z) and the test distribution Q := q(z) are two close distributions but there exists a distributional shift between , and R adv (S) is the adversarial sensitivity under adaptively scaled adversarial attacks on P. Define the l-th order distributional shift distance as: For example, when the training and test distributions are two multivariate Gaussian distributions, namely N (µ p , Σ p ) and N (µ q , Σ q ) respectively, then . When training and test distributions are close, namely µ p − µ q 2 and Σ q Σ −1 p − I 2 are bounded, Based on previous PAC-Bayes bounds [31], we provide a theoretical analysis in Theorem 2.
where S = {δ : β −1 δ p ≤ }, and the remainder term When |D| → +∞, we have Rem → 0, then the generalization bound is mainly determined by 2 . Theorem 1 and Theorem 2 reveal that the generalization bound L Q relates to the adversarial sensitivity R adv (S) under adaptively scaled adversarial attacks.
Traditional adversarial training algorithms with constant scales optimize the adversarial sensitivities under L p constrained adversarial attacks, while the proposed ASAT optimizes the adversarial sensitivities under adaptively scaled adversarial attacks, which are more precise estimations of the generalization bound. It implies that adaptive scales may help alleviate overfitting and improve the generalization ability compared to traditional adversarial training.
Since β is unknown, we adopt adaptive scales determined by data scales and gradients in Sec. 3.2.3.

Experiments
In this section, we first introduce the task and dataset details, then introduce the baseline models and experimental settings. Last, we report experimental results.  Table 1. and y is the ground truth, then these metrics are defined as: here x last denotes the volume of the last time slot, and ACC is the accuracy of whether the volume increases or decreases compared to the last time slot.

Baselines
We implement multiple moving average baselines and three neural network baselines in our experiments.

Moving average baselines
We adopt moving average baselines as our baselines, which are commonly used in technical analysis in finance.
In statistics, suppose a series of volumes v 1 , v 2 , · · · , v T , three commonly adopted estimations are: • Naive forecasting. The naive forecasting algorithm uses v T as the prediction. In our experiments, we try to adopt the volumes of yesterday or last time slot as the prediction.

Linear
The linear model is formulated as where x ∈ R 160 is the flattened vector of the concatenated input of both 12-slot and 20-day history.

LSTM
The long-short term memory (LSTM) [32] networks can capture features and long-term dependency in the entire sequences of data. Following [33], we implement two one-layer LSTM models with the global attention mechanism [34] to generate the representation vectors of previous time slots and 20-day history. Before feeding input into LSTM, we first adopt two linear layers to project the inputs into two higher dimensional spaces, respectively.
Then, we concatenate the representation vectors and feed them into a linear layer to regress the log volume. The input size, hidden size, and output size of LSTM models are 200.

Transformer
We also implement a six-layer Transformer [ of the last layer into a linear layer to regress the volume.
The input size, hidden size, and output size are 200, where the hidden states are split into 8 heads.

Settings and Choices of Hyperparameters
We train every model for 5 epochs and report the test performance on the checkpoint with the lowest valid loss. We

Experimental Results
After finding the best configurations of hyperparameters, the experimental results are shown in Table 2

Analysis
In this section, we first conduct an ablation study.
Then, we compare the adversarial robustness of neural networks under adversarial training. We also adopt ASAT to enhance existing methods. Last, we probe the sensitivities of different input dimensions and analyze the decision bases of neural networks, both on the dataset level and single instance level.

Ablation Study
We conduct an ablation study to investigate the influence of different mechanisms in our proposed ASAT.
In ASAT, as analyzed in Sec. 3.2.3, scales α i are determined by both |x i | and |g i | in Eq. (22). We try both α i = |x i |/max |x i | and α i = |g i | −1 /max |g i | −1 to verify the effectiveness of the proposed adaptive scales. The proposed ASAT algorithm minimizes the multi-step risk averaging target in Eq. (19). We try to adopt the traditional

Adversarial Robustness with ASAT
We compare the adversarial robustness of baseline models and models with adversarial training under multiple adversarial attacks on the five-minute dataset.
We first adopt multiple adversarial sensitivity indicators, including adversarial sensitivities under L 2 or L +∞ attacks and mean dimension-wise adversarial sensitivity.
Results are shown in We also visualize multiple adversarial attacks with several larger attacking scales. Experimental results are shown in Fig. 2. Similarly, we can conclude that ASAT can improve both the generalization ability and the adversar-    can improve the robustness of neural networks and achieve similar robustness to the 20-day average.

ASAT for Enhancing Existing Methods
We also implement two typical existing high-performance time-series prediction methods on the hourly dataset: Deep Auto-Regression (DeepAR) [36] and Multiple time scale LSTM (M-LSTM) [37]. The base models of both DeepAR and M-LSTM methods are LSTM models.
In DeepAR [36], we choose the Gaussian distribution as the prior distribution and adopt two linear layers to predict µ and σ of the predicted distribution respectively.    Namely, the prediction process and the training loss are: (y|µ, σ) = − log p(y|µ, σ) where h is the representation vector of the encoder, matrices U 1 , U 2 , and vectors b 1 , b 2 are learnable parameters.
µ, σ is the output of the neural networks, and y is the ground truth. In the inference phase, the prediction of the volume is µ. In M-LSTM, hidden states are split into four parts, and the time steps to propagate hidden states of four parts are 1, 2, 4, and 8, respectively.
The experimental results are in Table 5. Results show that DeepAR can enhance LSTM models in the volume prediction task, while M-LSTM cannot. ASAT can enhance these methods, which validates its effectiveness. Besides, the Transformer model enhanced with ASAT has a similar performance to the high-performance DeepAR [36] method enhanced with ASAT. The Transformer model is a high-performance time-series prediction model.

Probing Dimension-wise Sensitivities
We probe the dimension-wise sensitivities of different input dimensions in linear models, trained with three random seeds, the visualization results are shown in Fig. 3 and As shown in Fig. 3a and Fig. 3b    However, as shown in Fig. 3c, Fig. 4c and Fig. 4d, models with ASAT tend to pay close attention to some recent time slots and are sensitive to them, which accords with human intuition and is reasonable since recent time slots are more important in time series regularization. It also indicates that ASAT can alleviate overfitting since models with ASAT often focus on similar dimensions during multiple training.

Analysing Decision Basis on Individual Instance
In Eq.(12), the dimension-wise sensitivity can be defined on a dataset D or a single data instance (x, y). When  We choose the five-minute dataset, = 1, and randomly choose a data instance in the dataset to plot the dimension-wise sensitivities on the dataset and a single data instance. As shown in Fig. 5, Fig. 5a shows the dimension-wise sensitivities of a Transformer model trained with ASAT on the dataset. We can see that the model mainly concerns the volumes of nearest time slots and volumes in 20-day history. We also probe the decision basis on the data instance, whose volumes are shown in Fig. 5b.
The volumes of the 9-th nearest time slots, the 14-th nearest day in history, and the 4-th nearest day in history are relatively low compared to other time slots or days in his-tory. Both the baseline Transformer (shown in Fig. 5c) and the ASAT Transformer (shown in Fig. 5d) are sensitive to these time slots or days. They pay close attention to these abnormal volumes and the corresponding prices.
Therefore, their decisions on this data instance can capture these abnormal fluctuations. Moreover, the model with ASAT tends to be more sensitive to abnormal fluctuations, which indicates that ASAT can help the decision bases of models be more reasonable and more explainable.  [1,2,3,4,5,38,27,39], existing studies also concern improving the adversarial robustness and generalization ability of neural networks via adversarial training [2,5,6,7,8]. Adversarial training is widely adopted in both the computer vision field [2,5,6,7,40,41] and the natural language processing field [11,10,8,9,42].

Generative Adversarial Networks
Generative Adversarial Networks (GANs) [24] are also applications of adversarial training, which train both the generator model and the discriminator model with an adversarial learning target. GANs have been widely used in multiple machine learning fields for generating highquality synthetic data, such as computer vision [43,44], natural language processing [45,46], and time series synthesis [14,15]. In the finance field, several GAN models [14,15,16,17,18,19,20,21,22,23] are proposed to generate or synthesize better time-series data that can preserve temporal dynamics in time-series.

Conclusion
In this work, we propose adaptively scaled adversarial Proof. Since when |δ i | ≤ , δ =i = 0, according to Taylor series with Lagrange form of the remainder, there exists η near x, such that: Appendix A.2. Proofs of Proposition 2 Proposition 2. Consider a neural network that its loss function is convex and M -smooth in the neighborhood of θ 2 . Then, where g = ∇ x L(f (x, θ), y) and its i-th dimension is g i . Namely R (i) adv ( ) ∝ E|g i | holds approximately when is small.
Proof. Since when |δ i | ≤ , δ =i = 0, according to Taylor series with Lagrange form of the remainder, there exists η near x, such that: therefore, Appendix A.3. Proofs of Theorem 1 Theorem 1. Consider a neural network that is convex and M -smooth in the neighborhood of θ, suppose g denotes the gradient, g = ∇ x L(f (x, θ), y), then its adversarial sensitivity or risk under our proposed adaptively scaled adversarial attacks where S = {δ : α −1 δ p ≤ } is: Proof. Since L is convex and M -smooth near θ, According to Holder Inequality,   Proof. Suppose 1 p + 1 q = 1, according to the Holder Inequality, δ T g = (α −1 δ) T (α g) ≤ α −1 δ p α g q = α g q .

Appendix B. Hypothesis Tests
We also conduct student-t tests on the hourly dataset to verify that the proposed ASAT outperforms baselines statistically significantly (p < 0.05). The null hypothesis is that the outperform of ASAT is worse than the baseline.
The degree of freedom is df= 5 + 5 − 2 = 8 and the critical t-value is 1.86. Experimental results of hypothesis tests are in Table B.6. We can see that t-values of MLP and Transformer models are larger than 1.86. For the LSTM model, the improvements in ASAT are not significant.


CU-UD: text-mining drug and chemical-protein interactions with ensembles of BERT-based models

Identifying the relations between chemicals and proteins is an important text mining task. BioCreative VII track 1 DrugProt task aims to promote the development and evaluation of systems that can automatically detect relations between chemical compounds/drugs and genes/proteins in PubMed abstracts. In this paper, we describe our submission, which is an ensemble system, including multiple BERT-based language models. We combine the outputs of individual models using majority voting and multilayer perceptron. Our system obtained 0.7708 in precision and 0.7770 in recall, for an F1 score of 0.7739, demonstrating the effectiveness of using ensembles of BERT-based language models for automatically detecting relations between chemicals and proteins. Our code is available at https://github.com/bionlplab/drugprot_bcvii.

I. INTRODUCTION
Relation extraction is an important task in the biomedical domain that aims to identify the semantic relationship between the biomedical entities mentioned in the text. Information about these interactions can be helpful for other biomedical research such as drug discovery, database curation and question answering systems. In recent years, there have been lots of efforts to automatically extract the relations from the biomedical articles (e.g., (1,2)). BioCreative events have focused on different relation extraction tasks over the years. One of the tasks in BioCreative VI was ChemProt (3), which aims to promote the development and evaluation of systems that can automatically detect relations between chemical compounds and genes/proteins in PubMed abstracts. This task has been continued as a DrugProt task in BioCreative VII (4). In this paper, we discuss our participation in this task, describing ensemble approaches and incorporating multiple BERT-based language representation models.

A. Dataset
Organizers of the DrugProt track created a manually annotated corpus of 5,000 abstracts from PubMed (4). The dataset was split into a training set (3,500 abstracts), a development set (750), and a test set (750) ( Table I). Following the settings of previous BioCreative tracks, the test set consists of a large collection of records containing a subset of a total of 750 Gold Standard records that was used for evaluation purposes (3). * co-corresponding

B. Preprocessing
Since all relations in the DrugProt corpus are enclosed in one sentence, our system considers interactions only occur within a sentence. We split the raw text into individual sentences by Stanza sentence splitter (5). These sentences are passed onto the BERT framework which tokenizes them using Wordpiece tokenizer (6).
Given a sentence with multiple entity mentions marked, we treated the relation extraction problem as a multi-class classification problem. Curators of the DrugProt corpus annotated thirteen key relations of biomedical importance. A sentence instance with a marked protein and a chemical mention is considered one of these thirteen DrugProt relation classes if there is a relationship between a protein and chemical occurrences. Otherwise, this instance is negative. There may be multiple interacting entities and non-interacting entities marked in the same sentence. We consider only a single chemical-protein pair in that sentence per sentence instance. Therefore, we formulate the task as predicting one of the fourteen relation types for a sentence instance.
We consider two entity tagging schemes (Table II). In the first scheme, we anonymized target named entities in a sentence using pre-defined tags to prevent overfitting to entities In the second scheme, we directly inserted the start and end entity markers into the original sentence. For example, we used the format <PROTEIN-B> human type 12 RDH <PROTEIN-E> to enclose the protein mention "human type 12 RDH" in the sentence.

C. Model Development
We hypothesize that an ensemble system that combines results of different models could lead to better predictive performance than using a single model (1). The superior performance of ensemble models is achieved when there is enough diversity among the individual models (8). In this study, we propose two different ensemble approaches which incorporate multiple BERT-based language representation models. Here, we chose BERT as the backbone because recent studies demonstrate it achieved state-of-the-art performance on a number of natural language understanding tasks, including relation extraction task (7,9,10).
We explore various architectures and pre-trained models to achieve the goal of diversifying BERT-based models. We will discuss these techniques and the ensemble algorithms in the following subsections. Due to the space limitations, we explain only the most relevant parts of the BERT model to the DrugProt task in the following. We refer the reader to the original paper for the details.

1) Single models:
We built three different architectures. The first one uses the last layer of the [CLS] token for the relation extraction task because it aggregates the information of the whole sentence. Fig. 1a demonstrates the fine-tuning process of BERT on the relation extraction task. A sentence instance is tokenized and embedded as in the original BERT paper. Inspired by the work of Su et al. (11), we also explore two other methods of summarizing the information in the last layer: Long Short Term Memory (LSTM) on the last layer (Fig. 1b) and attention mechanism on all units of the last layer (Fig. 1c).
To obtain diversity among the models used in our ensembles, we chose three BERT-based models in our approach, BioBERT (12), PubMedBERT (13), and BioM-ELECTRA  (14,15), since they differ in some pre-training aspects as well as masking, vocabulary, and the learning objectives techniques. The pre-training of BioBERT used PubMed and PMC after the initial pre-training of BERT, whereas PubMedBERT pretraining using PubMed alone started from scratch. These two models are BERT-Base with 12 bidirectional self-attention heads. To study the effectiveness of large models for this task, we also used BioM-ELECTRAL. BioM-ELECTRAL is pretrained on PubMed abstracts with no previous pre-training and is based on the ELECTRA architecture (14).
2) Class-weighted loss function: Since we have created a negative instance for each non-interacted chemical and protein in the same sentence, we have substantially more negative examples in the training dataset. To mitigate this unbalance,  we assigned class weights inversely proportional to their respective frequencies.
3) Majority voting: Majority voting is the simplest yet effective weighting method for output fusion. We select the relation type that gets the most votes out of five individual models. If the models cannot agree on the majority vote, we assign this instance as negative.

4) Meta-learning method:
We also employ another ensemble method using multilayer perceptron (MLP) as a metalearner to learn from individual classifiers. Specifically, we use MLP to combine the outputs of each individual model in an ensemble (Fig. 2). After each individual model is fine-tuned, we extract the last output layer of the [CLS] token from each model. We then use the development set that we reserve for the ensemble learning to obtain the activation vector for the [CLS] token. Afterward, we concatenate them into one vector and pass it to another MLP. We insert a softmax layer on top of the MLP to obtain final predictions.
Individual model predictions were combined using a stacking approach involving a MLP with one hidden layer having 512 units as described in Section II-C4. Run 3 involves the same Group B and uses both stacking and the majority  While participant teams were provided with the training data set and the development data set, we combined them together to develop the models. For Runs 1 and 2, we first used 70% of the total data for training, 20% for development (tuning the hyper-parameters for the individual models) and 10% to train the ensemble models. For the final models, we trained the individual models on the original training data and ensemble models on the development data with the chosen hyper-parameters.
For Run 3, we used the same approach of 70-20-10 split, but repeated this process five times (Fig 3). Thus, we obtained a group of five models corresponding to each 70-20-10 partition. Next, we trained a separate MLP for each group of models using the 10% as described above. As a result, we have five predictions for each instance corresponding to each MLP. Finally, we combined them using majority voting to obtain the final predictions.
For Run 4, we fine-tuned five BioM-ELECTRA L models using 80% of total data in each fold and then combined their predictions using majority voting.
Run 5 involved only one model, and hence 100% of the data was used for its fine-tuning.
Our models are implemented in Tensorflow 1.14 and all the experiments are carried out on NVIDIA 2080Ti GPUs. We employ AdamW optimizer with an initial learning rate of 2e-5. The maximum epoch is set as 10. Table III and IV shows the results we obtained for 5 Runs on the development set in iteration 1 (Fig 3) and the official test set, respectively. Runs 3 and 4 have the best F1 performance, which is substantially better than the average F1 across participant teams. Interestingly, F1 score of the single large model is almost the same as those of the ensembles of  Table V shows the results of individual relation types using our best system (Run 3). Our system obtained the best performance of 0.8902 in F1 for the ANTAGONIST relation, and 0.8570 for the INHIBITOR relation. As shown in Table I, INHIBITOR is the most frequent relation in the dataset. Thus its classification is probably well trained. On the other hand, ANTAGONIST is one of the least frequent relations. Further error analysis is needed to understand why our model can effectively extract it with such few cases. Compared to results of individual relation types by other Runs, we observed that Runs 3 and 4 have a uniform increase in F1 performance for almost all individual relation types. On the other hand, Run 1 has the best F1 performance on the AGONIST type with a 2.5% increase over Run 3.

IV. CONCLUSION
In this manuscript, we describe our submission in the BioCreative VII DrugProt task. The results demonstrate that our ensemble system can effectively detect the chemicalprotein relations from biomedical literature.