Accurate and Computational Efficient Joint Multiple Kronecker Pursuit for Tensor Data Recovery

This paper addresses the problem of tensor completion from limited samplings. Generally speaking, in order to achieve good recovery result, many tensor completion methods employ alternative optimization or minimization with SVD operations, leading to a high computational complexity. In this paper, we aim to propose algorithms with high recovery accuracy and moderate computational complexity. It is shown that the data to be recovered contains structure of Kronecker Tensor decomposition under multiple patterns, and therefore the tensor completion problem becomes aKronecker rank optimization one, which can be further relaxed into tensor Frobenius-norm minimization with a constraint of a maximum number of rank-1 basis or tensors. Then the idea of orthogonalmatching pursuit is employed to avoid the burdensome SVD operations. Based on these, two methods, namely iterative rank-1 tensor pursuit and joint rank-1 tensor pursuit are proposed. Their economic variants are also included to further reduce the computational and storage complexity, making them effective for large-scale data tensor recovery. To verify the proposed algorithms, both synthesis data and real world data, including SAR data and video data completion, are used. Comparing to the single pattern case, when multiple patterns are used, more stable performance can be achieved with higher complexity by the proposed methods. Furthermore, both results from synthesis and real world data shows the advantage of the proposed methods in term of recovery accuracy and/or computational complexity over the state-of-the-art methods. To conclude, the proposed tensor completion methods are suitable for large scale data completion with high recovery accuracy and moderate computational complexity.


Introduction
The problem of data recovery from limited number of observed entries, referred to as matrix completion or tensor recovery problem, had attracted significant attention in the pass decade and had been applied in various fields such as recommendation system [1,2], Bioinformatics data [3], computer vision [4,5] and image data [6,7].
In practice, the matrix data we observed is often of low rank structure. Therefore, the matrix completion problem can be solved by finding a low rank matrix to approximate the data matrix to be recovered. However, the process of rank minimization is NP-hard [8]. In order to relax the non-convex NP-hard problem to a traceable and convex one, [9,10] proposed to minimize the matrix nuclear norm instead of matrix rank. Based on this concept, many algorithms, such as the singular values thresholding (SVT) [11] and singular value projection (SVP) [12] approaches. Unfortunately, a great number of singular value decomposition (SVD) operations or iterations are required for these methods, leading to high computational and storage complexity.
In real world systems, the observed data is sometimes of high dimensional structure. The recovery of these tensor data, which can be regarded as a generalization of that of matrix completion to the high dimensional case, is becoming more and more important. It is assumed that the observed data tensor contains a low rank structure and thus the missing entries can be recovered by the minimizing the tensor rank. In fact, several tensor decomposition methods and definitions of tensor rank are employed for the completion problem. Reference [13] proposed to minimize the number of rank-1 tensors from CANDECOMP/PARAFAC (CP) decomposition to recover the original tensor data. However, the decomposition is very costly, and it might fail to give reasonable results in some cases because of the ill-posedness property [14]. Other than the CP decomposition, the Tucker decomposition that yields a core tensor and corresponding subspace matrices of each dimension, are employed for tensor completion [15,16]. However, these methods require the knowledge of tensor rank, and might be not applicable in some applications.
Recently, some new types of tensor decomposition methods are proposed. By operating the 3-dimensional (3-D) tensors as 2-D matrices using the tensor product, the tensor singular value decomposition (t-SVD) is proposed [17]. Then the corresponding tensor rank and tensor nuclear norm [18] are defined and used for tensor data recovery. This method can transform the decomposition of one 3-D tensor to slice-wise matrix ones. However, for tensor with number of dimension higher than 3, it will fold all the dimensions higher than 3 into the third one, resulting in a lost of higher order information. Other researchers, proposed to employ the Kronecker tensor decomposition (KTD) based on the Kronecker structure of general R-D tensor to solve the tensor completion problem [19]. This approach employed the idea that folding an R-D tensor into a highorder tensor and then unfolding it into a matrix by general unfolding [20]. By reducing the tensor computation to matrix computation, a similar matrix rank can be defined, and good recovery results can be obtained. However, the computational complexity of the KTD based methods are still unsatisfying because of the SVD operation of large matrices.
Furthermore, in order to improve the efficiency of data recovery, some orthogonal rank-1 matrix and tensor pursuit approaches had been proposed [21,22]. These methods transfer the minimization of the matrix rank or tensor rank to an iterative rank-1 pursuit problem, and achieve a moderate data recovery result with low computational complexity. In this work, we will propose two novel tensor completion methods based on the minimization of tensor Kronecker rank as well as technique of matching pursuit, and apply the tensor data recovery methods to several real world data completion applications.
The remainder of the paper is organized as follows. Section 2 presents the notations and definitions. In Section 3, the tensor recovery algorithms are developed in details. Experiments are provided in Section 4 and conclusions are drawn in Section 5.

Notations and Definitions
The notations used in this paper is first defined. Scalars, vectors, matrices and tensors are denoted by italic, bold lower-case, bold upper-case and bold calligraphic symbols, respectively.
The Frobenius-norm of a tensor A is defined as .., m R ) . Furthermore, by stacking the first F-th dimensions of an R-D tensor A ∈ C M 1 ×M 2 ×···×M R one by one into one dimension and the remaining dimensions into the other [23], the tensor A can be unfolded into a matrix . . , K r and l r = 1, 2, . . . , L r . This folding procedure in fact folds the r-th and (r + R)-th dimensions of B.
We can now go on to define the Kronecker tensor product [19] of two tensors A ∈ This unfolding can be operated by first compute Y = fold {K} (X ) ∈ C K 1 ×···×K R ×L 1 ×···×L R and then calculate Y = {Y} 1→R . Then the Kronecker tensor decomposition (KTD) [19] can be defines as the Definition 2.3.
, and a g, f g 2 = b g, f g 2 = 1 the α g, f g is the amplitude of the f g -th term, and the K g is the g-th pattern. This factorization can be calculated using truncated SVD or non-negative matrix factorization [19] on X g for g = 1, 2, . . . , G. It is worth noting that the word 'pattern' here refers to using different shape to explore different structure information of the tensor data.
Following this definition, the tensor rank of X can be defined as the total number of Kronecker terms the joint KTD approach generates, which is G g=1 F g . Furthermore, when G = 1 and F G = 1, it is called Kronecker rank-1 tensor under the pattern {K G } [19].

Main Result
The task is to recover the R-D tensor X ∈ C M 1 ×M 2 ×···×M R from the partly observed tensor X , denoted as X where is the 1/0 set indicating the location of the observed entries. By writing the tensor X as a linear combination of several Kronecker rank-1 tensors under a given set of patterns {K 1 , K 2 , . . . , K G }, we get: where M g, f is the f-th Kronecker rank-1 tensor with the g-th pattern K g and θ g, f is the magnitude of M g, f for g = 1, 2, . . . , G and f = 1, 2, . . . , F. Note that here we assume that the ranks , the low rank tensor completion problem becomes: where θ 0 denotes the number of nonzero elements of θ. For noisy scenario, the optimization where is a small value. The (3) can be transformed into: where F is the maximum number of Kronecker rank. We can solve this problem by a greedy matching pursuit type method.

Iterative Multiple Kronecker Tensor Pursuit
In this Iterative multiple Kronecker Tensor Pursuit (IKTP) method, we update one Kronecker rank-1 basis tensor M n and one magnitude θ n with one pattern K h ∈ {K 1 , K 2 , . . . , K G }, where h = (n − 1) %G + 1 in an iteration. Therefore, the problem (4) can be rewritten as Supposing that after the (k − 1)-th iteration, (k − 1) Kronecker rank-1 basis tensors M 1 , M 2 , . . . , M k−1 and their weights θ k−1 = [θ 1 , θ 2 , . . . , θ k−1 ] T were obtained. In the k-th iteration, the target is to pursue a new Kronecker rank-1 basis tensor M k from the residual term where u k = v T k = 1. The u k and v k are pair of top left and right singular vectors of R k , which can be solved using the power method [23,24]. Then the tensor M k can be folded form By reshaping the tensor X and , the optimal solution θ k of (7) is calculated as Then we can go on to the (k + 1)-th iteration until all N ranks are solved, and this iterative multiple Kronecker Tensor Pursuit (IKTP) method is now summarized in Tab. 1.
Step 2: update (u k , v k ) from R k according to (6).
Step 4: Until a stopping criterion is satisfied. And give Y as the output.

Joint Multiple Kronecker Tensor Pursuit
The proposed IKTP approach updates one Kronecker rank-1 basis tensor according to one of the G patterns alternatively. When the number of ranks to be recovered is large, a lot of iteration is required, making the proposed approach highly computational inefficient. To overcome this problem, we propose a Joint multiple Kronecker Tensor Pursuit (JKTP) method that updates all the G Kronecker rank-1 basis tensors in one iteration under the given patterns.
According to (4), suppose that after the (k − 1)-th iteration, the Kronecker rank-1 basis tensors M 1, 1 , . . . , M G, 1 , M 1, 2 , . . . , M G, 2 , . . . , M 1, k−1 , . . . , M G, k−1 and their weights Under the given patterns {K 1 , K 2 , . . . , K G }, we define R k g = unfold {Kg} R k g and M g, k = unfold {Kg} M g, k for g = 1, 2, . . . , G, and proposed to solve the M g, k , g = 1, 2, . . . , G, by: where u g, k = v g, k = 1. The u g, k and v g, k are pair of top left and right singular vectors of R k g , which can be solved using the power method [25,26]. The new M g, k can be available by computing Then the optimal solution θ k of (10) is given by
Step 4: Until a stopping criterion is satisfied. And give Y as the output.

Economic Algorithms
In each iteration, the proposed IKTP algorithm will track all pursued bases and save them in the memory, leading to a high requirement of storage space. Furthermore, when the number of iteration is large, the LS solution (8) will compute the inverse of a large matrix, making it be rather computationally unattractive. To overcome these problems, an economic updating scheme which update α = [α 1 , α 2 ] instead of the magnitude term θ can be used in Algorithm 1, where the α can be figured by And the updated recovered tensor becomes Similarly, the economic updating technique can be used to update α = [α 1 , α 2 , . . . , α G+1 ] for the JKTP approach instead of (10), where α is calculated as And the recovered tensor can be obtained by

Performance Evaluation
In this section, we conduct experiments based on both visual data and SAR imaging data completion. The proposed IKTP and JKTP methods as well as they economic realization, IKTPecon and JKTP-econ, are evaluated, and the state-of-the-art algorithms, including Tucker [15], KTD [19], R1TP [21] and R1MP [22] are employed for comparison. For all algorithms, the stop criterion is met if the two conditions are reached: X k − X k−1 F < ε = 10 −6 or a maximum number of i max is satisfied. All the experiments are performed using MATLAB running on Inter(R) Core(TM) i7-8700@3.2 GHz for 100 Monte Carlo trials.

Image Recovery
First, the problem of image recovery by JKTP, JKTP-econ, IKTP, IKTP-econ, Tucker [15], KTD [19], and R1TP [22] methods from randomly sampled entries is tackled. Six images of dimensions 256 × 256 × 3 and 1024 × 1024 × 3 that separated into 2 groups are used for testing and they are shown in Fig. 1, we can infer that the completion result of the images from the first group might be better than that from the second one, as the pixels in the image form the former group is more related to each other. , and the observation rate is set to be 0.2. Note that in order to achieve a good performance, the patterns are suggested to be set as square as possible. By writing 'JKTP-i,' 'IKTP-i,' 'JKTP-econ-i' and 'IKTP-econ-i,' i = 1, 2, 3, 4, we mean that the proposed approach will use i patterns out if all 4 patterns for data recovery, and the patterns used are randomly selected in each Monte Carlo trial. The PSNR results of the two groups images data under numbers of iteration in Figs. 2 and 3, respectively. It is shown that the JKTP and JKTP-econ algorithms can converge within 15 iterations and the IKTP and IKTP-econ approaches will converge within 25 iterations when the n umber of patterns is larger than 3. In the remainder image and video experiments, the maximum numbers of iteration are set to be 20 and 30 for these methods, respectively, and we will use at least 3 patterns for proposed algorithms.  Fig. 1, and thus their recovery results are not shown here. It is shown that when more pattern used, better performance the proposed approached can achieve. Tucker and R1MP methods as well as the KTD approach with i = 1, 2, 3, 4 patterns are also included for comparison. It is shown that the JKTP and JKTP-econ methods give inferior results. For i ≥ 2, all the proposed algorithms can give better performance than the Tucker and R1MP approaches. Furthermore, the computational complexity of the proposed methods are much smaller than KTD, making them more attractive in dealing image data.   14.75, 9.51, 1.57 × 10 3 , 17.54 and 77.13 s. It is shown that the proposed methods run slightly slower than the R1MP methodology because the R1MP is matrix based rank-1 pursuit method while the proposed methods are tensor based ones. Comparing to the other methods, the computational complexity of proposed algorithms are much smaller, indicating the efficiency of the proposed approaches especially when the dimensions of the data to be recovered is large.

Video Inpainting
In this part, we test the algorithms JKTP, JKTP-econ, IKTP, IKTP-econ, Tucker [15], R1MP [19] and R1TP [21] under 'Gun Shooting' video [25,26] completion from randomly sampled entries, and the first and last frames of the video are shown in Fig. 6. Note that the KTD approach is not include for comparison because of its high computational complexity. The video is of dimensions 200 × 520 × 3 × 80.  Fig. 7 and the computational complexity results are shown in Fig. 8. Note that for the proposed methods, '−3' means that only the first 3 patterns are used, while '−6' means that all 6 patterns are employed in data recovery. Generally speaking, the proposed methods give the best PSNR performance when the number of frames used is higher than 20. It is worth noting that although the proposed methods with 3 patterns perform worse than R1TP method when the number of frames used is less than 20, when more frames are used, the recovery performance of the proposed approaches increase rapidly, and far outperform the other methods.

High Resolution SAR Imaging
In the end, we evaluate the algorithms with under-sampled high resolution SAR imaging data [27,28]. We employ the model shown in Fig. 9 to conduct our experiments as follows. First, scan the test scene to generate the data X by using the SAR system, and the scanning parameters follow those in [27]. Note that the true image I can be computed from the data X by the FFT process. Furthermore, the data X can be sampled randomly under an observation rate to obtain M . Finally, the tensor recovery methods can be applied for the recovery of the original data asX , and the recovered image can be computed asÎ. The mean square error (MSE) between X andX is used to evaluate the performance of the proposed four methods, R1MP [27], Tucker [18] and KTD [22] methods. The tested high resolution SAR imaging scene with the size 6200 × 12000 is shown in Fig. 10a.  Considering that the size 6200 × 12000 of the SAR imaging data is too large to be computed effectively by one computer, we will not complete the whole SAR data in one tune. We complete small subdata whose data dimensions are 620 × 50, 620 × 200 and 620 × 400 in one time and repeat the process for 2400, 600 and 300 times, respectively, for the whole image. The patterns for the proposed methods with different completed data size are K 50 The situations of different observation percentages ranging from 10% to 30% are tested and a maximum number of 40 iterations is used. The SAR imaging scene sampled under the observation percentage 10% is shown in Fig. 10b, and the SAR imaging scenes with the size of 620 × 200 recovered by all the algorithms under 10% observation rate are shown in Figs. 10c-10f. The completion results of the four methods including MSE performances and CUP times are listed in Tab. 4. Generally speaking, the proposed methods can give the best performance with moderate computational complexity among all the algorithms. tensor based algorithms, the results also show the advantage of the proposed methods in terms of computational complexity.
Funding Statement: The work described in this paper was supported in part by the Foundation of Shenzhen under Grant JCYJ20190808122005605, and in part by National Science Fund for Distinguished Young Scholars under grant 61925108, and in part by the National Natural Science Foundation of China (NSFC) under Grant U1713217 and U1913203.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.