In this subsection, we first introduce some fundamental tensor concepts and terminology. An extensive review of tensors are discussed in Kolda and Bader (2009).

For the definition of a tensor, there are several definitions at different abstraction levels (Boumal, 2023). Here, an Nth-order tensor is defined as a multidimensional array with N dimensions (N > 2), which is an extension of a matrix to higher order. Here, N represents the order of the tensor and the number of dimensions of the tensor, also known as modes. Analogous to the rows and columns of a matrix, we can create subarrays by fixing some of the given tensor's indices. A fiber of a tensor is a vector created by fixing all but one index of a particular mode (Filipović and Jukić, 2015). A third-order tensor X has row, column and tube fibers, denoted by x(:, j, k), x(i, :, k), x(i, j, :), respectively. A slice of a tensor is defined as a two-dimensional matrix, obtained by fixing all but the indices of two particular modes, e.g., a third-order tensor with horizontal, lateral, and frontal slices, denoted by x(i, :, :), x(:, j, :), and x(:, :, k), respectively.

We can also transform a tensor into a matrix via reordering the elements of an Nth-order tensor, known as matricization, unfolding or flattening. The mode-n unfolding of an Nth-order tensor X∈ℝI1×I2×⋯×IN is obtained by keeping the nth mode fixed and concatenating the slices of the rest of the modes into one long matrix, denoted by X_(n), which is an I_n × ∏_{j ≠ n}I_j matrix for n = 1, 2, ⋯ , N, where the symbol ∏ is continued product. The tensor element (i₁, i₂, …, i_N) is mapped into the element (i_n, j) of the matrix X_(n), where

The mode-n product of a tensor X∈ℝI1×I2×⋯×IN with a matrix U∈ℝJ×In is denoted as X×nU∈ℝI1×⋯×In-1×J×In+1×⋯×IN. If there exists a series of distinct modes, the order of multiplication in the mode-n product is independent, that is,

According to the mode-n unfolding of the tensor, we have the following formula

The inner product of two tensors X∈ℝI1×I2×⋯×IN and Y∈ℝI1×I2×⋯×IN with the same size is denoted as a sum of the products of the corresponding elements, that is,

By an extension of the matrix property, the Frobenius norm of the tensor X∈ℝI1×I2×⋯×IN is defined as the square root of the sum of the squares of all elements:

The distance between the tensor X and Y is defined as ||X−Y||_F.

The Kronecker product of two matrices A ∈ ℝ^I×J and B ∈ ℝ^K×L is denoted by A ⊗ B, and yields a matrix of size (IK) × (JL), defined as:

Tensor decomposition refers to expressing a tensor as a sequence of elementary operations on other simple arrays. There are many tensor decompositions. We introduce two major tensor decompositions, including CANDECOMP/PARAFAC (CP) decomposition (Cattell, 1944) and Tucker decomposition (Tucker, 1963).

For a given tensor A∈ℝI1×I2×⋯×IN, the CP decomposition is defined as the factorization of A into the minimum number of linear combinations of rank-1 tensors:

where “°” denotes the outer product of vectors, and ar(k)∈ℝIk(k=1,2,  ,R) is the factor vector.

The Tucker decomposition is based on an extension of the CP decomposition, a form of higher-order principal component analysis. The Tucker decomposition of the tensor A is expressed as:

where C∈ℝR1×R2×⋯×RN is the core tensor and denotes the degree of interaction between the different components, and Un∈ℝIn×Rn is the nth-mode factor matrix, denoting the principal components on each mode, usually orthogonal. In general, R_n is much smaller than I_n for n = 1, ⋯ , N, so the core tensor C can be considered as a compressed version of the original tensor A, which reduces the computational complexity to some extent.

The matricized form of the Tucker decomposition is as follows:

The Tucker-rank of the Nth-order tensor A∈ℝI1×I2×⋯×IN, also known as the multilinear rank (Kasai and Mishra, 2016), is denoted as an N-dimensional vector composed of the ranks of the various mode unfolding matrices of the tensor:

where rank(A_(n)) is the rank of the mode-n unfolding matrices, representing the dimension of the vector space spanned by the mode-n fiber. In contrast to the CP rank, the Tucker-rank can be easily obtained by computation, while the CP rank can only be determined empirically or by multiple experimental comparisons to search for the optimal value (Mu et al., 2014). Diagram of the Tucker decomposition of a third-order tensor is shown in Figure 2.

The LRTC aims to infer and fill missing entries of the tensor from partially observed values via utilizing the low-rank structure of the tensor, which is an extension of the low-rank matrix completion.

Let M ∈ ℝ^p×q be the observed incomplete matrix and Ω be the index set, i.e., the observed positions of the matrix elements. Recovering the matrix M based on the observed partial elements is referred to as low-rank matrix completion, that is,

where the complete matrix X ∈ ℝ^p×q is obtained via the above optimization problem with constraints.

As a broadening of the low-rank matrix completion (Equation 2), the LRTC (Ji et al., 2009; Liu et al., 2013; Yuan et al., 2019a) is formulated as the following optimization problem:

where X,T∈ℝI1×I2×⋯×IN, T is the incomplete tensor containing the missing elements, and X is the completed tensor. The optimization problem is built with a rank minimization constraint, which ensures that the output tensor has the lowest rank possible for tensor completion. In the LRTC, there is also a category of optimization based on tensor decomposition, formulated as minimizing the discrepancy between the real observations and the predicted values to complete the missing entries under the given tensor rank. That is,

where R is a bounded constraint on the rank of the X and the tensor X fulfills a specific tensor decomposition model. The traditional models are the completion models based on the CP and Tucker decompositions, e.g., the optimization model for tensor completion using CP decomposition is stated as follows:

The LRTC algorithm minimizing tensor rank is more widely used because the low-rank decomposition-based completion algorithm is impractical in some situations in reality, especially it is difficult to select an appropriate tensor rank for use in the algorithm when there are few observed elements. However, since the tensor rank is not unique, this makes the LRTC algorithm very challenging (Liu et al., 2013). Besides the CP rank and the Tucker rank in Equation (1), there are various forms of tensor rank, such as TT rank, tubal rank, TR rank, etc. Selecting one of these forms will result in different LRTC models.

In the LRTC, since the computation of the tensor CP rank is an NP-hard problem (Hillar and Lim, 2013) and the Tucker rank has the excellent property of being easy to compute, optimal completion models based on the Tucker rank of the tensor are widely used (Romera-Paredes and Pontil, 2013; Liu et al., 2019a; Yu et al., 2023). That is,

where the parameters α_i are pre-specified weights and satisfy α_i ≥ 0, ∑α_i = 1. The Tucker rank of the tensor is represented in the form of a weighted sum. Gandy et al. (2011) suggests the unweighted model of α_i = 1(i = 1, 2, ⋯ , N).

Since the rank of the matrix is non-convex, it is difficult to guarantee to finding a globally optimal solution in the optimization (3). Since the nuclear norm of the matrix is shown to be the tightest convex envelop of the matrix rank (Candes and Recht, 2012), the completion model based on the nuclear norm minimization is proposed as follows:

where the nuclear norm ‖X(i)||*=∑jσj(X(i)), and σ_j(X_(i)) denotes the j-th largest singular value of the matricization X_(i). The ∑i=1Nαi‖X(i)||* is defined as the nuclear norm of the tensor X, denoted as ||X||_*. The tensor nuclear norm is also known as the tensor schatten-q norm (Ji et al., 2009; Signoretto et al., 2010). Based on the nuclear norm, the concept of low-rank matrix completion is generalized to high-dimensional data, where the low-rank property of the tensor is often used as a necessary assumption to limit the degrees of freedom of missing entries (Kressner et al., 2014).

The tensor nuclear norm minimization completion model (TNNM) (4) can be expressed as the following:

where λ > 0 is a regularization parameter.

Since the matrix nuclear norm in the model (5) is interdependent and cannot be optimized independently, auxiliary variables for splitting need to be introduced for simplifying solution algorithm. Liu et al. (2013) presented three solution algorithms, which are simple LRTC, fast LRTC, and high-accuracy LRTC. Based on tensor decomposition, Liu (2013); Liu et al. (2014a) proposed completion models for the factor matrix nuclear norm minimization and the core tensor nuclear norm minimization. Nevertheless, the majority of TNNM algorithms involve SVDs of several large matrices in each iteration, resulting in high computing costs.

In the subsection, we describe the proposed CTNM-QR algorithm and the solution procedure in detail. Specifically, we adopt the ideas of the thin QR decomposition to update the factor matrix on the framework of the core tensor nuclear norm minimization (CTNM) and propose an efficient low-rank tensor completion algorithm.

In the CTNM model based on the Tucker decomposition (Liu, 2013; Liu et al., 2014b), the following assumption is stated:

Assumption 1. Let X∈ℝI1×I2×…×IN with rank_tc = (R₁, R₂, …, R_N) and satisfy Tucker decomposition X = C ×₁U₁ ×₂U₂ × … ×_NU_N, where the core tensor C∈ℝR1×R2×…×RN, and the orthogonal factor matrix Ui∈ℝIi×Ri. Then

where ||X||_* and ||C||_*, denotes the nuclear norm of the original tensor and the core tensor, respectively.

According to the Assumption 1, the completion model for the CTNM is:

where I_Ii is an I_i × I_i Identity matrix. Owing to the interdependence of the unfolding matrices, we introduce auxiliary tensors Vi∈ℝR1×R2×⋯×RN(i=1,2,  ,N) to facilitate variable separation. Then, the model (6) is transformed into the following equivalent form:

For the optimization problem Equation (7), we propose an efficient CTNM-QR algorithm based on alternating direction multiplication method (ADMM) to address it. The ADMM is a commonly employed optimization method that decomposes a problem into a sequence of sub-problems to calculate the optimal solution, and is adept at proficiently addressing optimization problems featuring multiple non-smooth terms in the objective function (Glowinski, 2014; Boţ and Nguyen, 2020; Han, 2022). The partial augmented Lagrangian function for the model (7) is as follows:

where Y_i(i = 1, ⋯ , N) is the Lagrange multiplier tensor and μ > 0 is the penalty parameter.

The ADMM converts the original tensor completion problem (7) into solving the five subproblems of (9a–9e) by minimizing the augmented Lagrangian function (8), which iteratively updates the corresponding parameters. Specifically, we update Vik+1,Ck+1,Uik+1,Xk+1 and Yik+1 sequentially in the order of (9a–9e) with the fixed other parameters. The steps of the algorithm solution are as follows:

That is, for Vik+1 in the k + 1th iteration, (9a) can be written as the following optimization subproblem:

The iteration for Vik+1 for i = 1, 2, ⋯ , N can be effectively addressed through the application of the singular value thresholding (SVT) operator (Cai et al., 2010). It is easy to obtain a closed solution to the problem (10) as:

where fold_i denotes the transformation of the matrix into a tensor along the ith mode, if the SVD of a matrix Z is Z = Udiag(σ)V^T, SVTτ(Z)=Udiag(max{σ-τ,0})VT, max(·, ·) is the element-based maximization operator.

For C^k+1 in the k + 1th iteration, the optimization problem (9b) with respect to C can be rewritten as:

With other variables fixed, it is easy to see that Equation (12) is a smooth and differentiable optimization problem that can be derived by first-order optimality conditions for C (Liu, 2013; Liu et al., 2014b).

For X^k+1 in the k + 1th iteration, solving X by the following subproblem:

By deriving the Karush-Kuhn-Tucker (KKT) condition for Equation (13), it is easy to obtain the update form of X as XΩk+1=TΩ and XΩck+1=(Ck+1×1U1k+1×⋯×NUNk+1)Ωc, where Ω is the index set of the observed elements, Ω^c is the complement of Ω.

For the factor matrix Uik+1∈ℝIi×Ri(i=1,2,…,N) in the k + 1th iteration, the optimization (9c) is equivalent to the following optimization problem:

This is a peculiar Tucker decomposition model in which U_i, i = 1, 2, …, N is an orthogonal matrix. HOSVD (De Lathauwer et al., 2000a; Chao et al., 2021) and HOOI (De Lathauwer et al., 2000b) are frequently used in the computation of Equation (14). Nevertheless, these approaches need alternating iterations of solving Equation (14) and involve SVDs of multiple tensor matricizations in each loop iteration, which causes the algorithms to converge quite slowly. To enhance the approaches, we adopt thin QR decomposition instead of HOOI for solving U_i, which essentially avoid the SVDs during factorization and improve the computing efficiency.

Specifically, we first isolate the ith factor matrix U_i from the Tucker decomposition since the order of multiplications in the tensor mode-i product is not relevant. Then the equivalent form of the optimization model (14) for the solution of the ith factor matrix U_i is as follows:

where the given tensor Zik=Ck+1×1U1k+1×2⋯×i-1Ui-1k+1×i+1Ui+1k⋯×NUNk. The problem (15) is nonconvex for all variables, yet it is convex for each block tensor. When the variables C,X,{Uj}j=1,j≠iN are fixed, the optimization (15) with respect to the factor matrix U_n is a smooth and differentiable optimization problem and can be solved via the first-order optimality condition. Hence the following equation holds:

where (Zik×iUi-Xk)(i) and Zi(i)k are the mode-i unfolding matrices of the tensors Zik×iUi-Xk and Zik, respectively. Then the least-squares solution without taking into account the orthogonality constraints on the factor matrices can be obtained by

Motivated by the work of Shi et al. (2015), we adopt the thin QR decomposition for the orthogonality constraints of the factor matrices, that is,

where QR(·) denotes the thin QR decomposition of the matrix. At this step, the proposed CTNM-QR algorithm merely needs to carry out QR decomposition of a matrix with the computational complexity O(∑i=1NIi×Ri2). While the computational complexity of the HOSVD algorithm in terms of U_i iterations is O(∑i=1NIi2×Πj≠iIj), the HOOI algorithm is O(∑i=1NIi2×Πj≠iRj), where R_i ≪ I_i. This indicates that the significant computational complexity associated with the SVD can be effectively mitigated by employing QR decomposition to resolve the factor matrix.

In the following, we discuss several implementations issues. To effectively solve the objective function (8), the penalty parameter μ is adaptively varied via setting the initialization of μ to μ₀, and iteratively increasing μ^k by μ^k+1 = ρμ^k, where ρ ∈ (1.0, 1.1] and μ₀ is usually a smaller constant. For the iteration among the problems (9a–9e), we principally use two convergence criteria. One is the relative change satisfaction between two neighboring recovery tensors X^k+1 and X^k:

where k is the number of iterations and tol is the tolerance error. The alternative halting condition is k ≥ K, where K is the maximum number of iterations. The iteration ends if any of the two conditions is met. We summarize the implementation of the CTNM-QR as Algorithm 1.

Let {Ck,U1k,…,UNk,V1k,…,VNk,Xk} be a sequence generated by the Algorithm 1. Although the CTNM-QR is a nonconvex ADMM algorithm, the subsequence are convergent and converge to the KKT point of the problem (6) proved by Liu et al. (2014b). In the subsection, we only illustrate the convergence of the improved part of the algorithm for the subproblem (9c).

Denote the objective function in Equation (15) as f(Un)=12‖Zn×nUn-X‖F2, where (Un)TUn=IIn(n=1,2,…,N). The function f(U_n) is convex on U_n without considering orthogonality constraints, so minimizing f(U_n) yields U~nk=X(n)k-1(Zn(n)k-1)T(Zn(n)k-1(Zn(n)k-1)T)-1. Then the following equation holds:

Denote g(Ck,Un)=12‖Ck×1U1k⋯×n-1Un-1k×n+1Un+1k-1⋯×NUNk-1×nU~n-Xk‖F2. By carrying out a thin QR decomposition of U~nk, we have U~nk=UnkQk, where Unk is an orthogonal matrix. Then we have

where C~k=Ck×nQk. Therefore, when C,X,{Ui}i=1,i≠nN is fixed, the objective function f(U_n) tends to be decreasing during the iteration process, and the sequence {U1k,…,UNk} is the convergence sequence.

In the subsection, we analyze the complexity of the proposed CTNM-QR algorithm. The running time of the CTNM-QR algorithm is primarily consumed in solving the auxiliary variables with the SVT operator and some multiplicative computations. Suppose the size of the input tensor is I₁ × I₂ × ⋯ × I_N and Tucker rank is R₁ = ⋯ = R_N = R, the time complexity of solving Equation (9a) is O(NR^N+1), the time complexity of solving X in (9d) is O(R∏i=1NIi), and the time complexity of computing U_n in Equation (16) is O(NR∏i=1NIi+R2∑iIi). Consequently, the total time complexity of the Algorithm 1 is O(T(N+1)R∏i=1NIi), where T is the total number of iterations. When the N and I_i, i = 1, 2, ⋯ , N are fixed, the time complexity of the CTNM-QR algorithm mainly dominates on R. Fortunately, the R generally takes small values due to the low-rank assumption of the tensor.

For comparison, we list the major computational complexity of several classical tensor-completion algorithms in Table 1. Although the per-iteration complexity of TT-SGD (Yuan et al., 2019b) is incredibly low, the total number of iterations required is higher and more time-consuming than the proposed CTNM-QR algorithm in practice. In addition, the Polak-Ribiere nonlinear conjugate gradient algorithm of CP-WOPT (Acar et al., 2011) also suffers from the problem of being time-consuming. Compared with the nuclear norm minimization algorithms such as SiLRTC, HaLRTC, CTNM, etc., the proposed CTNM-QR algorithm saves some computational complexity. The CTNM-QR is consistent with the generalized higher-order orthogonal iteration (gHOI) algorithm (Liu et al., 2014b) in terms of the algorithm's primary computational cost.

In the simulation studies, we generate the low rank synthetic tensor T∈ℝI1×I2×⋯×IN by the Tucker decomposition model, i.e., T = C ×₁U₁ ×₂U₂ × … ×_NU_N, where C∈ℝR1×R2×…×RN, Ui∈ℝIi×Ri(i=1,2,…,N) and their entries independently follow the standard normal distribution. This implies that the Tucker rank of the tensor T is (R₁, R₂, …, R_N).

In the concrete implementation, we set the tolerance to tol = 10⁻⁵ for all algorithms, the maximum number of iterations to K = 6 × 10⁵ for the TT-SGD, and K = 500 for the rest of the algorithms. For the FaLRTC, gHOI, CTNM, and CTNM-QR, we set the weights α_i = 1/N, i = 1, 2, ⋯ , N. The smoothing parameter of the FaLRTC is set to α_i/I_i, i = 1, …, N, and the rest of the parameters are maintained at their default values. For the CP-WOPT, the tensor rank R is set to 15. For the TT-SGD, the TT rank is set to R₁ = ⋯ = R_N−1 = 15. For the gHOI, CTNM, and CTNM-QR, we set μ0=10-3, μmax=1010, ρ = 1.05.

We consider two cases including the fourth-order tensor with the size 40 × 40 × 40 × 40 and the fifth-order tensor with the size 20 × 20 × 20 × 20 × 20. The Tucker-rank of the tensor is set to R₁ = ⋯ = R_N = 5, where the tensor orders N are 4 and 5, respectively. We consider four distinct MRs 30, 50, 70, 80% for the test tensor and replicate 100 simulations. Table 2 presents the average results (RSE and time cost) for 100 independent experiments. The outcomes reveal that the CTNM-QR method yields more accurate solutions with the less time costs, and outperforms other algorithms in the term of the RSE regardless of whether it's a higher order tensor or a lower order tensor. Although the gHOT method has a somewhat shorter running time than the CTNM-QR at the low MRs, the time costs of both methods are quite close and the RSE of the CTNM-QR is less than one of the gHOI.

We also investigate the convergence of the CTNM, gHOI, and CTNM-QR methods on synthetic tensor data with the Tucker-rank (R₁, R₂, R₃) = (10, 10, 10) and the size 100 × 100 × 100, as shown in Figure 3A, where the abscissa indicates the number of iterations and the vertical indicates the log-value of the relative change of the X^k. The results display that the relative change of the convergence of the CTNM-QR drops more quickly and smoothly than that of CTNM and gHOI. Figure 3B illustrates the relative changes of the convergence of the CTNM-QR approach at 40, 60, and 80% MRs. From the Figure 3B, we find that even at the high MR, the CTNM-QR converges fast in fewer than 100 iterations.

In this subsection, we confirm the efficacy of the CTNM-QR by performing recovery on RGB real-color images, where each image is represented as a third-order tensor. The low-rank nature of color images has been demonstrated in Figure 1, which allows excellent recovery of incomplete images utilizing the low-rank structure of the tensor. The low-rank color pictures chosen for this experiment include Flower, Starfish, Sky, Lena, and Apple. The original images are portrayed in Figure 4. The images are 500 × 500 × 3 tensors and encompass portraits and various natural landscapes. The image pixels are normalized to [0, 1]. We carry out multiple experimental comparisons for every image in order to determine the optimal rank value.

In the concrete implementation, the tolerance error of all algorithms is set to tol = 10⁻⁴, the maximum number of iterations K = 5 × 10⁵ for TT-SGD and K = 500 for the rest of the algorithms. For the FaLRTC, gHOI, CTNM and CTNM-QR, we set the weights αi=[1,1,10-3], and the smoothing parameter for the former to αi/Ii for i = 1, 2, 3. For the CTNM-QR, we set μ0=10-3, μmax=1010, ρ = 1.05. The rest of the algorithm parameters are preserved at the default settings.

For the color images, we utilize random uniform sampling to make them missing, and the average reconstruction results from 100 independent experiments, including RSE and time cost (seconds), are reported in Table 3, where MR is set to 40, 60, and 80%. The best results are highlighted in bold in the table, and it is clear that the proposed CTNM-QR method achieves the best RSE values in all images and outperforms the other competing methods in terms of CPU time in a substantial number of cases. For example, when recovering Sky with an MR of 80%, the RSE of the CTNM-QR is improved by 0.036 and 0.071 over gHOI and the original algorithm CTNM, respectively, and by 0.119 over the FaLRTC. In terms of running speed, it outperforms the CTNM and the TT-SGD by around 4 and 7 times, respectively. Furthermore, in the majority of image completions, the CP-WOPT method takes the longest. It is obvious that the CTNM-QR can efficiently mine the low-rank image spatial correlation to precisely recover the missing tensor in a reduced iteration time.

Figure 5 depicts visual comparisons of the observed images with the images following the various completion algorithms, where the missing rate of incomplete images is 60 and 80%, and the missing entries are displayed in black. From Figure 5, it can be seen that the best visibility outcomes are obtained by the CTNM-QR method (Figure 5G), and the recovered image is very clear and close to the original image. In most instances, the results of the CP-WOPT contain some fairly ambiguous parts, and this is the most acute, especially for the Lena and Apple. In the completion of missing data for images with an 80% MR, the FaLRTC and CTNM exhibit suboptimal performance, particularly on the Flower and Lena datasets. Although the results of the gHOI and TT-SGD appear to be more discernible in comparison, the finer details of the images have not been fully recovered, indicating a certain degree of information loss. It implies that the proposed CTNM-QR method can recover the images accurately and efficiently.

In the subsection, we evaluate the performance of the CTNM-QR algorithm via the MRI images completion of the brain. MRI images of the brain are acquired through the use of magnetic resonance imaging, which utilizes a strong magnetic field and harmless radio waves. It has an influential role in the fields of medicine and neuroscience, being a key tool used by physicians to diagnose diseases, study brain activity, and guide treatment programs. Here, we analyze brain MRI data with the size 236 × 236 × 180. Similarly, medical MRI data is approximated as a low-rank tensor and thus can be recovered with low-rank completion algorithms. For comparison, the tolerance value for all algorithms is set to tol = 10⁻³, the tensor Tucker-rank is set to R₁ = R₂ = R₃ = 30 for the gHOI, CTNM, and CTNM-QR, the tensor CP rank is set to 40 for the CP-WOPT, and the TT rank is set to R₁ = R₂ = 20 for the TT-SGD. We choose a slice of the MRI images for completion due to their enormous computation cost.

Table 4 presents the recovery accuracy (RSE) and time cost (seconds) of every algorithm at varying MRs, where the MRs are set to 30, 50, and 70%, respectively. The findings show that the RSE of the CTNM-QR, gHOI, and CTNM are nearly equal at MR values of 30 and 50%. However, when the MR increases to 70%, the RSE of the CTNM and gHOI dramatically rose. In terms of the runtime, the proposed CTNM-QR method typically runs faster than others, especially more noticeable at the higher MR. In addition, Figure 6 displays the recovery results of the brain MRI slice data at the 70% MR, where the CTNM-QR outperforms the other algorithms in terms of both recovered vision and efficiency.