## 1 Introduction  
In modern industry, the increasing complexity of mechanical equipment assigns great significance to the normal operation of every single part. An unexpected failure could lead to the instability or even breakdown of the mechanical system [1; 2]. In engineering practice, incipient fault prognosis is expected to provide early warning before the occurrence of catastrophic failure. Evidently, vibration signals are effective tools for condition monitoring and fault diagnosis, mainly because of their revelations of dynamic information and their convenience for measurements. In [PERSON]'s study, various measurement systems for evaluating the vibrations of rolling bearings were thoroughly compared [3]. Generally, due to the complexity of the transmission, the vibration signal with incipient fault often incorporates a large amount of the interference apart from the fault characteristic components. The interference may include the harmonic components that are related to operational status and stochastic noise, such as electromagnetic interference. Thus, the feature of impulse fault is inundated and covered unobservable in the data. To date, plenty of studies have been carried out to identify the impulse fault components. In the theory of probability and statistics, the kurtosis, as a measure for the shape of data distribution, is conducive to capture the outliers (e.g., abnormal impulses) in the data [4]. [PERSON] et al. further investigated that peculiarity, enforcing the impulses in the signal, in an improved algorithm named the maximum correlated kurtosis deconvolution (MCKD), to indicate the periodic fault transients [5]. MCKD has received growing attention, and it can be categorized as a single-dimensional method. With respect to two-dimensional methods like time-frequency analysis, the time-frequency distribution (TFD) obtained by short time Fourier transform (STFT) [6], Wigner distribution, wavelet transform [7], etc., provides an effective tool for analyzing non-stationary vibration signals, as it incorporates information from both the time-domain and frequency-domain. Also, it is interesting to note that relevant studies based on matrix analysis emerged. Given the non-negativity of TFD, [PERSON] et al. combined the time-frequency analysis with non-negative matrix factorization for feature extraction in [8]. [PERSON] et al. utilized singular value decomposition (SVD) to identify relevant information from background noise, which was combined with envelope analysis to detect a fault feature of rolling bearings in [9]. [PERSON] et al. introduced the ratios of neighboring singular values to distinguish different beating fault types [10]. Specifically, SVD resorts to the Hankel matrix to obtain a two-dimensional representation for a time series, which signals the trend of the high-dimension analysis. The study of [PERSON] and [PERSON] reports a combination of time-frequency analysis and manifold learning for signature enhancement [11], in which phase space reconstruction (PSR) is used for constructing a three-dimensional TFD manifold. It provides finesse for enhancing the impulse feature, but inconvenience appears, as the TFDs must be changed into vectors for manifold learning. The techniques mentioned are admittedly efficacious for feature extraction in two-dimensional representations; however, the applications are limited when confronted with higher-order data. For data in a multi-dimensional representation, such as 3-way RGB images, it is efficient to directly use multi-way arrays rather than matrices for data processing. Tensors are the multi-way extension of matrices and provide a pertinent representation of natural data with multi-dimensional formation [12]. Tensor factorization, which recovers meaningful hidden structure or principal components from the multi-way array, was gradually popularized and widely applied in various fields, such as brain science [13; 14], computer vision [15], speaker recognition [16], and user preference and modeling [17].