Matrix and tensor completion using tensor ring decomposition with sparse representation

Completing a data tensor with structured missing components is a challenging task where the missing components are not distributed randomly but they admit some regular patterns, e.g. missing columns and rows or missing blocks/patches. Many of the existing tensor completion algorithms are not able to handle such scenarios. In this paper, we propose a novel and efficient approach for matrix/tensor completion by applying Hankelization and distributed tensor ring decomposition. Our main idea is first Hankelizing an incomplete data tensor in order to obtain high-order tensors and then completing the data tensor by imposing sparse representation on the core tensors in tensor ring format. We apply an efficient over-complete discrete cosine transform dictionary and sparse representation techniques to learn core tensors. Alternating direction methods of multiplier and accelerated proximal gradient approaches are used to solve the underlying optimization problems. Extensive simulations performed on image, video completions and time series forecasting show the validity and applicability of the method for different kinds of structured and random missing elements.


Introduction
In various applications of multi-dimensional data processing, the underlying data tensors or matrices are corrupted by outliers and/or incomplete due to imprecise data acquisition process, or corruption by artifacts [1,2]. The problem of recovering an original data tensor from its partially observed entries is called tensor completion and has found many applications such as in computer vision, data mining, recommender systems, image inpainting and signal reconstruction. Due to its importance and many practical applications, many algorithms have been developed to solve this problem. Most of the methods are based on low-rank matrix and tensor decomposition. For a comprehensive overview of such algorithms, we refer to [3,4]. The applicability and performance of these algorithms highly depend on the probability distribution upon which the data components are missing. Many of these algorithms work quite well only when the underlying missing components are distributed uniformly and fail to recover the data tensor when a whole part of a data tensor is missing, e.g. several sequential rows or columns of a frame are missing.
Motivated by such limitations of existing tensor completion techniques, we exploit tensor ring (TR) decomposition combined with Hankelization [5,6]. The idea of Hankelization was first proposed in [7] and recently generalized to tensors in [5] but so far has not been exploited extensively for the TR representation with sparse constraints. A key pre-processing step is a prior Hankelization procedure (to be discussed in section 3) which allows us to apply the TR representation and develop an efficient tensor completion algorithm. The original data tensor is recovered through a so-called inverse Hankelization or De-Hankelization [5]. In this paper, we propose a new approach for the tensor completion problem. Our main idea is based on applying dictionary learning technique to the block Hankelized high-order data tensor in TR format to yield sparse representations of the core tensors. More precisely, in contrast to the existing Hankelized tensor completion methods [5,8], we do not complete a Hankelized data tensor directly or using Tucker decomposition, instead we use constrained TR format. We first initialize the core tensors and then they are updated sequentially using an efficient dictionary learning procedure applied to each core tensor. The size of core tensors is smaller than the original data tensor which leads to lower-scale optimization problems.
The problem is formulated as a convex optimization problem which is solved efficiently by the alternating direction methods of multiplier (ADMM) [9]. Our main contributions in this paper are as follows: • Applying Hankelization of incomplete data in order to obtain higher order tensors and represent it in the TR format. • Exploiting the dictionary learning technique for sparse representations of the TR core tensors.
• Developing an optimised algorithm, which allows to reconstruct images corrupted by structured distortions, like rows, columns or blocks.
The rest of this paper is organized as follows. In section 2, basic notations and concepts used throughout the paper are introduced. In section 3, we describe the Hankelization technique and related concepts. The proposed method and methodology is presented in section 4. Extensive results are illustrated in section 6, and finally, a conclusion is given in section 7.

Notations and preliminaries
In this section, basic notations used throughout the paper are introduced and they are adapted from [10]. A scalar is denoted by standard lowercase letters, e.g. a, vectors or 1st-order tensors are denoted by boldface lower case letters, e.g. X, matrices or second order tensors are denoted by boldface capital letters, e.g. X and higher order tensor are denoted by bold underlined capital letters, e.g. X. Also a tensor graphically represented as, see figure 1(a). An (i 1 , i 2 , . . . , i N )th element of an Nth-order tensor X ∈ R I1×I2×···×IN is denoted by x i1,i2,...,iN = X(i 1 , i 2 , . . . , i N ). The trace and Moore-Penrose pseudoinverse of matrices are denoted by 'Tr' and †, respectively. The multi-index is defined as be an Nth-order tensor, then the n-unfolding of the tensor X is of size I i and denoted by X ⟨n⟩ whose elements are defined as (see figure 1) Similarly, the mode-n unfolding matrix of size I n × i̸ =n I i is denoted by X (n) whose components are For a graphical illustration see figure 1(b). Slices are sub-tensors generated by fixing all indices except two of them and as a result, they are matrices. For example, for a tensor X ∈ R I1×I2×I3 , the slices X(:, :, i), i = 1, 2, . . . , I 3 , are called frontal slices. The notation 'Vec' denotes the vectorization operator which stacks the columns of a matrix. The inner product of two tensors X , Y ∈ R I1×I2×···×IN is defined as ⟨X, Y⟩ = i1 i2 . . . iN a i1,...,iN b i1,...,iN and the induced Frobenius norm is ∥X∥ F = ⟨X, X⟩.

Tensor train (TT) and TR decompositions
Tensor networks (TNs) decompose higher-order tensors into sparsely interconnected small-scale low-order core tensors [10]. TT [11] and TR 3 (TR) [12][13][14][15] are two powerful and relatively simple TNs representing a higher order tensor as a train and ring (or chain) of third order core tensors, for a comprehensive study on TNs and applications in physics and machine learning, we refer to [16][17][18][19][20][21][22][23][24][25][26]. The TT and the TR decompositions are also known as matrix product state (MPS) and matrix product state with periodic boundary condition (MPS-PBC) in quantum physics [16,27]. They have found many practical applications such as machine learning [21,[28][29][30], compressing deep neural networks [31][32][33], tensor completion [34][35][36][37], and hyperspectral image super-resolution [38,39]. The memory storage requirement of these two TNs scale linearly with the order of the tensor so they break the curse of dimensionality which is a main bottleneck in handling high order data tensors [40]. The TT decomposition can be considered as a special case of the TR decomposition and because of this, throughout the paper, we mainly focus on the TR decomposition. Let X ∈ R I1×I2×···×IN be an Nth-order data tensor. The TR decomposition represents the data tensor X into a sequence of latent core tensors G (n) ∈ R Rn−1×In×Rn , n = 1, 2, . . . , N. The element-wise representation of the data tensor X in the TR format can be expressed as where G (n) (i n ) ∈ R Rn−1×Rn is the i n th lateral slice matrix of the core tensor G (n) ∈ R Rn−1×In×Rn . The expanded form of (2), is and the N-tuple (R 0 , R 1 , . . . , R N−1 ) is called TR-ranks. Note that in the TR decomposition, we have R 0 = R N and it is also shown in [15] that the TR-ranks satisfy R 0 R n ⩽ rank X ⟨n⟩ for n = 1, 2, …, N. A shorthand notation for the TR decomposition is as follows [10]: (1) , G (2) , . . . , G (N) ≫ .
As mentioned earlier, the TT/MPS decomposition is a particular case of the TR/MPS-PBC decomposition, for R 0 = R N = 1. This means that the first and last cores in the TT decomposition are matrices and rest of them are 3rd-order tensors. For graphical illustration of the TT and the TR decompositions, see figures 2(a)-(c). For a comprehensive overview on fast algorithms for the TR decomposition, see [41].

Low-rank TR completion
Let X ∈ R I1×I2×···×IN be an incomplete data tensor with observation index tensor Ω ∈ R I1×I2×···×IN , defined as ..,iN is unknown, and its complement as The operator P Ω (X) projecting the data tensor X onto the observation index tensor Ω is defined as The task of low rank TR completion, can be formulated as the following optimization problem arg min where and the goal is finding a low rank TR approximation X. If the TR-ranks (R 0 , R 1 , . . . , R N−1 ) are known in advance, then the completion algorithm is called non-adaptive otherwise it is called adaptive where the TR ranks is obtained by the algorithm automatically. Our proposed algorithm is non-adaptive. For simplicity, we assume that R n = R, ∀n.

Block Hankelization
A matrix is called a Hankel matrix if all its elements along the skew-diagonals are constant. A Hankel tensor is a generalization of Hankel matrices. A data tensor X is called a Hankel tensor if all components x i1,i2,...,iN for which the quantity i 1 + i 2 + · · · + i N is fixed, are the same. The procedure of generating a Hankel matrix/tensor from a given data vector/matrix/tensor is called Hankelization.  We first review basic concepts of Hankelization taken from [5]. From a vector x ∈ R N , we can generate a Hankel matrix. For example, for a given vector x = (x 1 , x 2 , . . . , x N ) T ∈ R N , and a window size τ , the operator H τ (x) : R N → R τ ×(N−τ +1) generates a Hankel matrix as The operator H τ (X) is also called delay embedding transform of vector X with window size τ [5,6]. Equivalently the Hankelization procedure can be expressed by a duplication matrix which provides a relationship between a Hankel matrix X H and a corresponding given vector X. Let X H , be a Hankel matrix, The operator fold (N,τ ) : R τ (N−τ +1) → R τ ×(N−τ +1) which returns back a vectorized form of a Hankel matrix to the original Hankel matrix is called folding operator, that is A sample illustration of this definition for a vector of size 8 and window size τ = 5 is presented in figure 3. It is seen that the duplication matrices includes block of shifted identity matrices. A block Hankel matrix is defined similarly, where instead of the individual entries of a matrix, its blocks are repeated along the block skew-diagonals of the matrix. More precisely, from a set of matrices, we can generate a block Hankel matrix; see figure 4 for an illustration of the block Hankelization procedure. A vector from which a Hankel matrix was constructed can be obtained through so called inverse Hankelization [5] as follows: where T is a duplication matrix. The concept of delay embedding can be generalized to higher order tensors via multi-way delay embedding [5,6,42]. The intuition behind this generalization is that the classical delay embedding can be considered as a tensor-vector multiplication i.e. X× 1 T = TX where for a vector as a first order tensor, we have only one duplication matrix T. For a matrix X ∈ R I1×I2 (which is a second order tensor), two duplication matrices Then the operator fold (N,τ ) (X× 1 T 1 × 2 T 2 ) reshapes the data matrix Y to a fourth order Hankel tensor of size For an Nth-order tensor, we should consider N duplication matrices T n ∈ {0, 1} τn(In−τn+1)×In , n = 1, 2, . . . , N and the multi-way embedded space is defined as [5] transforms an Nth-order tensor to an 2Nth-order tensor. The N-tuple (τ 1 , τ 2 , . . . , τ N ) indicates the window size of different duplication matrices corresponding to different modes. Illustration for Hankelization of a matrix and transforming it to a 4th-order tensor is shown in figure 4. The reverse procedure where the task is retrieving the original data tensor from the tensor X H was constructed, can be expressed as where unfold (I,τ ) = fold −1 (I,τ ) , and it basically transforms the 2Nth-order block Hankel tensor X H of size In the image processing community, the Hankelization procedure or equivalently delay/shift embedding is a technique that duplicates patches of an image with prescribed window sizes. It is known that a prior Hankelization procedure as a pre-processing step can capture the delay/shift-invariant features (e.g. non-local similarity) of signals/images [42]. This is an effective property to learn the hidden structure of the images which allow us to apply this technique in denoising and completion [5,42]. In particular, it is experimentally shown that for incomplete data tensors with unstructured missing patterns, this technique is quite efficient. For example, for a video with removing slices instead of pixels, the superiority of Hankelization in recovering such data tensors is shown in [5,42].
Although it is possible to use more sophisticated Hankelization, such as patch-based/block Hankelization [8], we however perform a relatively simple row Hankelization procedure due to its simplicity and efficiency. In our approach, each row of the frontal slices of the data tensor is Hankelized using a given window size and a corresponding duplication matrix. Then, from these matrices, block Hankel matrices are constructed.
Having computed all the block Hankel matrices, they are reshaped to higher order tensors, as illustrated in figure 5. By column Hankelization, we can obtain similar or the same results.
The existing algorithms are based on applying a completion technique to the Hankelized data tensor. In [5] Tucker decomposition was used, while in [8] the TT is exploited. Our work is quite different from them in the sense that we learn the core tensors of the TR decomposition of the incomplete data tensor through a dictionary learning technique and sparse representations.

Proposed optimization procedure and learning algorithm
In this section, we discuss in detail our approach. Our main idea is first Hankelizing the incomplete data tensor X into X H using row Hankelization approach and then recovering the incomplete Hankelized data tensor while imposing sparse representation constraint on the core tensors. We will experimentally show that such a formulation allows us to recover the structured missing patterns of an incomplete data tensor efficiently and achieve better performance compared to many other existing algorithms for this task. Unlike the methods in [5,8], where the completion approaches employ Tucker or the TT representations is applied to X H and then the original data tensor is reconstructed through inverse Hankelization, we consider the optimization problem (4) with block Hankelized data tensor X H (described in section 3 and figures 4 and 5) and exploit a TR decomposition as , G (2) , . . . , G where the core tensors G (n) , n = 1, 2, . . . , N, are estimated using dictionary learning. To this end, we formulate an optimization problem in which the constrained core tensors G, (n) n = 1, 2, . . . , N, are updated sequentially to minimize the cost function in (4). Moreover, inspired by [43,44], we imposed implicit sparse representation constraint on the core tensors G , (n = 1, 2, . . . , N). So we formulated the following constrained optimization problem: , G , . . . , G where G (n) (2) ∈ R In×Rn−1Rn is mode-2 matricization of the core tensor G (n) ∈ R Rn−1×In×Rn , D (n) ∈ R In×Rn−1Rn is a pre-defined dictionary, C (n) ∈ R RnRn−1×RnRn−1 is a sparse matrix and η > 0 is a regularization parameter. The method works well for quite a range of η from 0.05 to 0.3. We chose η = 0.082 to achieve best performance. It is worth mentioning that the optimization problem (8) is convex with respect to each core tensor G (n) but it is not convex with respect to the whole set of core tensors. As a result, there is not any guarantee to converge to a global minimum, however we have used a suboptimal approximation. Also, in contrast to [43], where the sparse representation of unfolding matrices are employed, we impose sparsity representation constraint on the core tensors of the TR network. This can reduce the computational complexity of the algorithm because of smaller size of the TR core tensors compared to the unfolding matrices and also improve the performance. Various dictionary learning algorithms can be used for computing the dictionary D (n) such as over-complete discrete cosine transform (ODCT) dictionary 4 [45,46], online dictionary learning [47] and group sparse dictionary learning technique [44]. We adopted for our purpose the ODCT method [46,48] because it is quite versatile and works for a broad range of data.
In the algorithm proposed in [43], the sparse representation constraints are imposed on the mode-n unfolding matrices while the nuclear norm minimization is applied on the n-unfolding matrices 5 . As a result, our proposed formulation (8) can be considered as a non-trivial extension and combination of techniques considered in [43] and [35], where on the one hand, we exploit l 1 -norm instead of nuclear norm minimization [35] and on the other hand, similar to [43], the dictionary learning technique is employed. It is worth mentioning that ignoring the sparsity constraints, we come up to the algorithm proposed in [43].
According to works by [15,35,49], the relation between the core tensor G (n) and the matricized tensor can be represented by 4 The MATLAB implementation of this dictionary can be found in https://github.com/kvdeepak/ksvd-denoising. 5 Note that mode-n unfolding and n-unfolding matrices are related to the Tucker and the TT decompositions, respectively.
where G (̸ =n) is a 3rd-order tensor of size R n × N i=1,i̸ =n I i × R n−1 obtained by merging all other core tensors whose slice matrices are computed as For brevity of notation, in the rest of the paper we define G (n) := G (n) (2) and G (̸ =n) := G (̸ =n) <2> . By employing formula (9) and taking into account the fact that the core tensors are independent, we can optimize the core tensors in optimization problem (8) independently while keeping other core tensors fixed as follows: for n = 1, 2, …, N, where The augmented Lagrangian function corresponding to the constrained optimization problem (10), can be constructed as where B is a matrix representing the Lagrangian multipliers and λ is penalty parameter. While the method works well for λ between 0.1 and 1.1. We set λ = 0.5 in our simulations. All hyper-parameters were optimized to provide the best performance. Based on Lagrangian function (11), the ADMM [9] update rules for solving (8) are straightforwardly converted to simpler optimization problems as follows: It is worth to note that the ADMM algorithm is a workhorse approach for solving a variety of optimization problems such as sparse inverse covariance selection, generalized and group lasso problems, etc. For a comprehensive study of this technique and related problems, see [9]. Moreover, this strategy has been utilized for solving many tensor completion problems: please see [43,[50][51][52][53], and for a comprehensive list of references see the review papers [3,4] and the references therein.

Solving the sub optimization problem (12)
The first optimization sub-problem (12) can be solved via the standard gradient descent method. To be more precise, the gradient of its objective function of (12) in kth iteration is where J = L G (n) k , C (n) , B and the first and second terms of (15) are the gradient of the first and second terms of the objective function of optimization problem (12), respectively. The gradient of the first term of (15) has been derived in [35] and the gradient of the second term can be derived by the following straightforward computations: and taking the gradient with respect to G (n) k .

Solving the sub-optimization problem (13)
In order to solve optimization problem (13), we applied accelerated proximal gradient (APG) method [54,55]. Assuming that the E C , we have the following cost function: The first step in the proximal gradient algorithms is considering an auxiliary proximal variable Z k , and computing a quadratic approximation of E(C (n) k ) around Z k using the following equation: where the gradient ∇ C k F (Z k ) is computed straightforwardly by expanding F (Z k ) and taking gradient with respect to C k as The parameter c > 0 in (16) is a given parameter and in our simulation we set c = D (n) 2

F
. Substituting where As a result, T(C k , Z k up to a constant and it can be minimised instead of H C (n) k , Z k . In the proximal gradient algorithms, instead of minimizing E(C (n) k ), which is a non-differentiable optimization problem, the objective function T(C (n) k , Z k ) is minimized for a sequence of proximal variables Z j k . This leads to the following minimization problem which can be solved iteratively: where Finally, we obtained the following closed form solution: where the soft thresholding operator is defined as soft (g, τ ) = sign (g) max(|g| − τ , 0). Note that in (20), the soft thresholding operator is applied in an element-wise manner to all components of the matrix W j+1 k . Algorithm 1. Algorithm for tensor ring decomposition with sparse representation (TRDSR).
Input: An incomplete data tensor X ∈ R I1×I2×···×IN , the observation index tensor Ω, TR-ranks (R0, R2, . . . , R N−1 ), for simplicity, we assume that Rn = R ∀ n, a window size τ and regularization parameter λ > 0. Output: Completed data tensor X 1 Hankelize the data tensor X and Ω as X H and Ω H 2 Initialize the TR core tensors for the Hankelized tensor X H as G (1) , G (2) , . . . , G (N) with Gaussian random tensors 3 while A stopping criterion is not satisfied do 4 for n = 1, 2, …, N do 5 Compute dictionary D (n) for mode-2 matricization of the core tensor G (n)

6
Solve optimization Problem (10) for G (n) and C (n) using the update ADMM rules (12)-(14) 7 end There are several possible options for updating the proximal variable Z j+1 k , but in this paper we used the accelerated one introduced in [54] as follows: where t 1 = 1. The whole procedure is summarized in the pseudo-code of algorithm 1. The stopping criterion which we used was that the relative error be less than a predefined tolerance or the maximum number of iterations is reached.

Discussion on computational complexity
Most of state-of-the-art tensor completion algorithms exploit the nuclear norm minimization formulation in which the computation of the singular value decomposition (SVD) is required multiple times. This increases their computational complexity and also makes them relatively slow. On the contrary, we avoid SVD computation because we do not exploit the nuclear norm minimization formulation. For the sake of simplicity, let us assume I 1 = I 2 = · · · = I N = I and R 1 = R 2 = · · · = R N = R, then the computational complexity of TR-ALS [34] and TRLRF [49,51] are O(pNR 4 I N + NR 6 ) and O(NR 2 I N + NR 6 ), respectively, where p is constant between 0 and 1. The main operation in our algorithms is matrix-matrix multiplication in (15) where we need to multiply matrices G (n) k ∈ R R 2 ×In and G . This indicates the lower complexity of our algorithm compared with others

Experiments
In this section, we evaluate the proposed algorithm referred to as TRDSR using real datasets including images, videos, and time series signals. All simulations were conducted on a laptop computer with 2.60 GHz Intel(R) Core(TM) i7-5500U processor and 16GB memory. The relative squared error (RSE) and peak signal-to-noise ratio (PSNR) are used to evaluate and compare the performance of various algorithms. The RSE is defined as where X and X are the original data tensor with all observations and the reconstructed data tensor, respectively, while the PSNR is defined as PSNR = 10 log 10 255 2 /MSE ,  and 'num(.)' denotes the number of components of the tensor X. We used the Matlab toolbox poblano toolbox 1.1 [56] for solving the underlying non-linear optimization problems. Example 1. (Image recovery) In these experiments, we used the benchmark images 'Lena' , 'Barbara' , 'House' , 'Facade' , 'Peppers' and 'Giant' as shown in figure 6. All color images are of size 256 × 256 × 3. We have tried different types of missing patterns. For example, in the case of 'Lena' image, we removed randomly 20% of pixels grouped in form of white big spots with some holes in the picture (first image) or continuous scratching the image (sixth image) both of them are considered as structured missing. For the 'Barbara' image (second image), we removed more than 40% of whole columns and rows of pixels randomly. The last two 'Pepper' and 'Giant' images have 90% of missing pixels. The missing patterns of the 'House' , 'Facade' images were also structured distortions. Using the window size τ = (252, 252, 3), the incomplete images were first Hankelized into a 5th-order tensor of size 252 × 5 × 252 × 5 × 3 and represented in TR format before applying the completion algorithms. We compared the performance of our algorithm with some recent tensor completion algorithms under various conditions and for different missing rates. The algorithms which we used were TRLRF [57], TRALS [34], TR-WOPT [35], TT-WOPT [55], CP-WOPT [58], SPC-TV [59] and MDT [5]. In order to provide a fair comparison with other algorithms, we use the same TT/TR ranks for TRLRF, TRALS, TR-WOPT, TT-WOPT, CP-WOPT algorithms and tune the hyper-parameters of each algorithm to make performance as best as possible. The MDT and the SPC-TV algorithms were able to tune their ranks automatically. The reconstructed images are displayed in figure 7 and detailed PSNR and SSIM comparison are reported in table 1. From the results, it is seen that our proposed algorithm, in general, outperforms the other algorithms in most of the cases for various scenarios and missing patterns. Also, we examined our proposed algorithms in terms of robustness to TR rank selection. To this end, we considered the 'Lena' image with 90% random missing pixels and compared the RSE. The results are reported in figure 8(a) which show that our TRDSR algorithm achieves the lowest RSE compared to the other completion methods and also it is relatively insensitive to selected TR ranks. We should mention that one of the main challenging task in tensor completion problem is a good selection of rank and a larger rank does not always necessarily provides better results. This is mainly because of the over-fitting problem where the model is not selected appropriately. Because of this issue, it is of interest to have a tensor completion algorithm which is less sensitive to rank selection. Actually, in this experiment we wanted to show that different from other algorithms, ours is less sensitive to the choice of TR-ranks and it is more stable to various TR-ranks. We also compared the computational time of our method with other existing methods.
The results are reported in figure 8(b). It is seen that our algorithm is the fastest algorithm while at the same time, it provides better performance. Finally, we applied our proposed algorithm to an MRI images.
Here, we mainly exploited four different kinds of structural missing patterns (see figure 9). In the first experiment, we scratched the data image (first image in figure 9), in the second one we removed 10% of pixels and 10% of columns and rows randomly (second image in figure 9), in the third experiment, we removed 40% of pixels and putting some spots in the image (third image in figure 9), and in the fourth experiment, we removed 40% columns and rows (fourth image in figure 9). The recovered images together with indicated corresponding PSNR index are reported in figure 9. Our computer experiments indicate that the proposed algorithm also provides good results for recovering real-life medical images. Moreover, we compared the performance of our algorithm for the MRI image with some other tensor completion algorithm. These results are reported in figure 10. Here again, our algorithms provided better performance.

Example 2.
In this experiment, we consider the video recovery task. The dataset which we used was the GunShot video [34,60] of the size 6 100 × 260 × 3 × 85. We removed 80% of voxels randomly. We Hankelized video frames into an 8th-order tensor of the size 252 × 9 × 96 × 5 × 3 × 1 × 83 × 3 for reconstruction using the window sizes τ = (252, 96, 3, 83). The reconstructed results of TR-ALS, MDT, SPC and our proposed TRDSR algorithm are displayed in figure 11. From the results, the superiority of our technique is visible.
Example 3. In this simulation, we evaluate our proposed algorithm for reconstruction of time series generated by the Lorentz system which is illustrated in figure 13. The size of the Lorentz system was 1 × 4600, which we reshaped to 3rd-order tensors of size 46 × 10 × 10. We considered two kinds of missing patterns as follows: • Removing some parts of the signal in the first and middle parts along with a Gaussian noise (see figure 13).
• Removing 20% of samples and also removing 300 samples in the middle and last parts of the signal (see figure 14).
In both above cases, we used window size τ = (5, 44,5) and Hankelized the incomplete 3rd-order tensor into a 6th-order tensor of size 5 × 6 × 44 × 3 × 5 × 6. The constructed signals using our proposed algorithm 6 Video includes 80 frames of RGB frames of size 100 × 260 × 3. 7 The video can be downloaded from www.nhk.or.jp/archives/creative/material/category-list.html?i = 22.     compared to the MDT and the SPC-TV algorithms are displayed in figures 13 and 14, respectively. Our algorithm achieved the lowest RSE and this confirms the effectiveness and applicability of our method.

Conclusion
In this paper, a new and efficient tensor completion algorithm was proposed for recovering data tensors with random and/or structured missing entries. The main idea is a prior Hankelization of the incomplete data tensor after which the data tensor is completed through learning the core tensors of the TR representation with sparse constraints. The ADMM algorithm and APG approaches were adopted to solve the underlying optimization problems. Extensive simulations show that our proposed method is capable of recovering incomplete data tensors with different types of structured and random missing elements. Our algorithm exhibits higher performance in comparison to many existing tensor decomposition methods, while providing lower computational cost in comparison to other tensor completion approaches.

Data availability statement
Data sharing is not applicable to this article as no new data were created or analysed in this study. The analyzed data are available online.