A New Model for Tensor Completion: Smooth Convolutional Tensor Factorization

Tensor completion is the problem of filling-in missing parts of multidimensional data using the values of the reference elements. Recently, Multiway Delay-embedding Transform (MDT), which considers a low-dimensional space in a delay-embedded space with high expressive capability, has attracted attention as a tensor completion method. Although MDT has a high complementary performance, its computational cost is considerably high. Therefore, we propose a new model called smooth convolutional tensor factorization (SCTF) for tensor completion based on a delay-embedded space. The proposed method is small in computational complexity because of its concise model of rank-1 decomposition in the delay-embedded space, and because it does not directly perform optimization in the delay-embedded space. In addition, a smooth constraint term is assigned to the factor tensors as a prior data structure in the optimization to improve the completion accuracy further. In our experiments, we completed clipped and random missing image data, and confirmed that the proposed method achieved high completion accuracy without high computational cost.


I. INTRODUCTION
Several real-world signals such as image and video data are multidimensional, and tensors are mathematical models that represent such signals [2], [8], [14], [20], [24], [27]. Such data, modeled as tensors, are often corrupted by measurement errors and missing observations. Tensor completion is the task of filling-in the missing values of the tensor data using the values of the reference elements [24], [37], [48], [52]. Because tensor completion is an ill-posed problem, we consider the prior structure in the target tensor to narrow down the solution set. The completion value should be appropriate as per the properties of the analyzed data, and the prior structure includes smoothness [15], [48], [53], nonnegativity [17], [39], [43], sparsity [23], low-rank [14], [24], etc.
2) Completion of the Hankelized tensor using Tucker decomposition. 3) Inverse MDT of the completed tensor. This method considers a delay-embedded space with a high expressive capability and exhibits higher completion accuracy than existing methods [5], [24], [47], [48], [53]. However, MDT-Tucker has the disadvantages of considerable time requirement and space computational complexity. For example, for an N th-order tensor of average size T , if the delay window size is τ , the space complexity is O(τ N T N ) and the time complexity is O(τ N +1 T N ); thus, the complexity increases exponentially with order.
In this study, we propose a novel smooth convolutional tensor factorization (SCTF) model, which decomposes a tensor into two smooth factor tensors by convolution instead of a product. Figure 1 shows a schematic of the algorithm. This model implicitly implements tensor decomposition in the delay-embedded space, whereas optimization is performed in the original space. The model is based on the relationship between the inverse MDT of the rank-1 model and the cyclic convolution of the factor tensors. In addition, because it is a rank 1 model, the SCTF is simpler than the MDT-Tucker model, which considers the Tucker decomposition model. These properties are expected to reduce computational complexity. In addition, a smoothness constraint was imposed on the factor tensors to further narrow the solution set. Our contributions can be summarized as follows: • We have mathematically proven that tensor decomposition based on inverse MDT has sufficient representation ability in rank-1 decomposition.
• Based on the relationship of inverse MDT of rank 1 decomposition and cyclic convolution of factor tensors and the introduction of smooth prior structure into factor tensors, we proposed a new tensor completion model named smooth convolutional tensor factorization (SCTF).
• We derived a solution method of the proposed SCTF with the Majorization-Minimization (MM) algorithm [18], [30], which is expected to provide a stable optimization in which the cost function decreases monotonically. Moreover, we exploit the equivalence of cyclic convolution in the time domain and Hadamard product in the frequency domain to reduce computation time. The remainder of this paper is organized as follows: related works in Section II, a review of MDT in Section III, the proposed method in Section IV, experiments using the proposed method in Section V, and conclusions in Section VI.

A. MATHEMATICAL NOTATION
We follow the basic mathematical notation in [46]; that is, we use ⊛ for Hadamard product in this paper. In addition, we consider N matrices U n ∈ R I n ×R n (n = 1, . . . , N ) and an N -th order tensor X ∈ R R 1 ×···×R N . The all-mode product is denoted as (1)

II. RELATED WORKS
t-SVD [50] is a convolutional tensor decomposition method as well as the proposed method. It can achieve accurate tensor recovery based on group theory. t-SVD considers a new SVD for tensors by using some convolution, and the rank in t-SVD (tubal rank) is defined as the number of non-zero singular values. Since it is difficult to minimize the tubal rank directly, its convex relaxation is usually employed. The convex relaxation of tubal rank is given by the sum of singular values based on t-SVD, and it is called as the tensor VOLUME 11, 2023 67527 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. nuclear norm (TNN). Low-rank approximation in the t-SVD is substituted for a problem of minimizing TNN, which has the advantage of incorporating a global structure. TNN [50] is a typical model for tensor completion problems based on t-SVD, and PSTNN [19] is a further developed model. PSTNN suppresses the excessively low rank of the estimated tensor by considering partial sums of only small singular values in the tensor nuclear norm. RTF [9] and UTF [10] have also been proposed as models that avoid the high computational cost of these t-SVD models. RTF considers a factorization model of a low-rank tensor of small size and a dictionary (orthogonal) tensor. Since t-SVD is applied only to low-rank tensors of small size, the computational cost is lower than that of other t-SVD models. On the other hand, UTF uses the fact that the TNN is transformed to the minimum sum of the Frobenius norms of the two low-rank tensors so that the algorithm does not directly compute t-SVD in its calculation. UTF achieves very fast inference despite t-SVD model. However, these methods differ from the proposed method because it uses only the third-order tensor, and the convolution operation is performed only in the channel direction. Also, unlike the proposed model, these models do not have a smoothness term. CNNM [22], [21] is a mathematical model of nuclear norm minimization of convolutional tensors applied to image completion and time series prediction. This research shows the equivalence of nuclear norm minimization of the convolution tensor and sparse approximation in Fourier space. However, the relationship between inverse MDT and cyclic convolution, and smoothness constraints is not discussed.

III. REVIEW OF MDT
This section summarizes the Multiway Delay-embedding Transform (MDT). Note that there are two types of MDT: noncyclic MDT [46] and cyclic MDT [45]. In this study, we consider a cyclic MDT. First, we discuss the Delay-embedding Transform (DT) for one-dimensional data (vectors), which is then extended to multidimensional data. Finally, we describe the drawbacks of MDT and a method, Fast-MDT-Tucker [45], for avoiding those.

A. DELAY-EMBEDDING TRANSFORM (DT) 1) OVERVIEW OF DT
A delay-embedding Transform (DT) is the transformation of data into a high-dimensional space representing a time delay. In physics, DT has been studied by reconstructing dynamic attractors from time-series data in a delay-embedded space [32]. Mathematically, the DT converts a vector into a Hankel matrix (Hankelization) [3]. When embedding an observed signal from the original space into a highdimensional space, it is assumed that the signal is represented by a low-rank and smooth manifold in the delay-embedded space [26], [29], [41]. Figure 2 shows the results of DT of the signal generated by the Lorenz system, indicating that the transformed signal is smooth and low-dimensional in the delay-embedded space. Based on this assumption, a lowrank approximation of the Hankel matrix is used in the data analysis [3], [11].

2) MATHEMATICAL OPERATION OF DT
Following [45], the DT for an observation vector x = (x 1 , . . . , x T ) T ∈ R T with a delay window size τ is defined as where the DT operation is denoted as H τ . Because X is a Hankel matrix, DT is also called Hankelization. Each row of X is identical to the local window of vector x. Notably, this study assumes that the signal is cyclic. Figure 3a shows a concrete example of a DT operation. The DT can be considered a linear operation. By using the duplication matrix S ∈ R T τ ×T , we obtain DT can be given by where fold (V ,v) : R Vv → R V ×v is a folding operator, that is reshaping from a vector to a matrix.
The pseudo-inverse of the DT can be expressed as the average of the anti-diagonal entries. Considering a matrix X ∈ R T ×τ , the inverse DT operation H † τ can be given by where S † := (S T S) −1 S T is a Moore-Penrose pseudoinverse of S, and vec(·) represents an operation of vectorization. We note that (S T S) −1 = 1 τ I. The tth element of the inverse DT is given by Figure 3 illustrates the DT and inverse DT matrix computations.

B. MULTIWAY DELAY-EMBEDDING TRANSFORM (MDT) 1) OVERVIEW OF MDT
A Multiway Delay-embedding transform (MDT) is an embedding transform for tensor data of two or more orders. Tensor decomposition based on the MDT has various applications, particularly for tensor completion. A representative model is MDT-Tucker, which considers the Tucker decomposition model of the Hankel tensor and has been applied to image completion [46] and time-series data [36]. Another model, the HT-RPCA, was proposed in [42]. Unlike general RPCA, HT-RPCA solves the rank minimization of the tensor Hankelized by MDT, instead of the rank minimization of the matrix. This method enables anomaly detection by considering the time series. Furthermore, the TT and TR decomposition models of the Hankel tensor have been proposed and applied to image completion and time-series data, [33] [34].

2) MATHEMATICAL OPERATION OF MDT
The DT can be naturally extended to an N -th order tensor X ∈ R T 1 ×···×T N of size T = (T 1 , . . . , T N ) ∈ R N . Let us consider N duplication matrices S n ∈ {0, 1} T n τ n ×T n (n = 1, . . . , N ) with a window size τ = (τ 1 , . . . , τ N ) ∈ R N (see Equations (3)). The MDT is defined using an all-mode product and folding as follows: where fold (V ,v) : is the folding operator from an N -th order tensor to a 2N -th order tensor. Conversely, the inverse MDT is defined as where an unfolding operator from the 2N -th order tensor of an N -th order tensor.

C. FAST-MDT-TUCKER
The methods introduced in Section III-B1 have the disadvantage of considerable time requirement and space computational complexity because of the Hankelization in each mode of the tensor. Fast-MDT-Tucker [45] was proposed to improve the high computational complexity of MDT-Tucker. This method focuses on the redundant structure of the Hankel matrix and improves the time complexity to O(NT N log NT ) and the space complexity to O(T N ) using two techniques: 1) Omission of duplicate computations.
2) Equivalence of cyclic convolution in the time domain and Hadamard product in the frequency domain. In 2), the Fast-MDT-Tucker exploits the relationship between the inverse MDT and cyclic convolution. Fast-MDT-Tucker provides a fast and accurate completion; however, only low-rank priors are available.
The proposed method is also an algorithm based on the relationship between the inverse MDT and cyclic convolution and similarly avoids the issues of MDT. Note that the low-rank model of the proposed method is not a Tucker decomposition but a rank-1 decomposition. In addition, the proposed method imposes a smoothness constraint on the factor tensors.

IV. PROPOSED METHOD
The proposed method solves the optimization problem by assuming that the observation tensor can be represented by a cyclic convolution of two smooth factors of the same size (See Figure 1). We describe the key theory behind the proposed method in Section IV-A, the smoothness constraints in Section IV-B, and the formulation and algorithm in Section IV-C.

A. KEY THEORY OF PROPOSED METHOD 1) RELATIONSHIP BETWEEN THE INVERSE DT AND CYCLIC CONVOLUTION
Any rank-R matrix X ∈ R T ×τ has a singular value decomposition that can be expressed as Because the inverse DT is a linear operation, H † τ (X) can be separated into rank-1 bases: From Equations (6) and (10), the tth element of the inverse DT for a single basis u r v T r is given by Now, let us consider the matrix P = (I τ O) T ∈ R T ×τ and set v ∈ R τ to be the same vector as the dimension of u ∈ R T , i.e., zero padding operation is given bỹ Note that the sizes of u andṽ are equal and the elements ofṽ are zero after the size of the delay window τ . From Equations (11) and (12), the inverse DT of the rank-1 basis uv T is given by where * denotes a cyclic convolution operation. From Equation (13), the inverse DT of the rank-1 basis can be formulated in terms of a cyclic convolution. Eventually, from Equations (10) and (13), the inverse DT of X is 2

) SUFFICIENT REPRESENTATION ABILITY EVEN WITH A RANK-1 MATRIX
We now discuss the rank-1 representation of X. From Equation (14), the rank X denotes the number of convolutional bases. The degrees of freedom of each convolutional basis determine the representational ability of the model. In this study, we consider X to be rank-1 and show that it has sufficient representation ability for vector reconstruction.
Rank-1 matrix model X = uv T ∈ R T ×τ can generate any x ∈ R T . Let us put where 0 τ −1 is a (τ − 1)-dimensional vector of zeros and we have This suggests that the inverse DT of a matrix, even rank-1, is over-parameterized and does not work as a model. Therefore, u and v must impose constraints to narrow the solution.

3) EXTENTION TO MDT
The properties of DT discussed in Sections IV-A1 and IV-A2 can be applied to the MDT. First, we show the relationship between the inverse MDT and N -dimensional cyclic convolution. Let us consider factor tensors A ∈ R T 1 ×···×T n , B ∈ R τ 1 ×···×τ n , and we define a := vec(A) ∈ R n T n , b := vec(B) ∈ R n τ n . We assume X ∈ R T 1 ×τ 1 ×···×T N ×τ N is given by where bunfold (V ,v) : R V 1 ×v 1 ×···×V N ×v N → R n V n × n v n is the unfolding operator from an 2N -th order tensor to the block matrix. We also define bfold (V ,v) : R n V n × n v n → R V 1 ×v 1 ×···×V N ×v N as the inverse transform of bunfold (V ,v) . Using the zero padding matrix P n = (I τ n O) T ∈ {0, 1} T n ×τ n (n = 1, . . . , N ), we define a tensorB = B × 67530 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. {P} ∈ R T 1 ×···×T n of the same size as A. The inverse MDT of X = bfold (τ ,T ) (ab T ) is derived by Thus, the inverse MDT is represented by an N -dimensional cyclic convolution. Furthermore, we show that the tensor which is folded from rank-1 matrix has sufficient representation ability, as discussed in Subsection IV-A2. Rank-1 tensor model of X := bunfold (T ,τ ) (X ) = ab T ∈ R n T n × n τ n can generate any X ∈ R T 1 ×···×T n . Let us where 0 n τ n −1 is an ( n τ n −1)-dimensional vector of zeros, and then we have Because this operation is equivalent to the cyclic convolution of a tensor X with only one element, the tensor X is derived identically (See Figure 4). As discussed in Section IV-A2, the inverse MDT of an unconstrained tensor, even rank-1, is over-parameterized and does not work as a model. In this study, additional constraints were imposed on A and B (see Section IV-B).

B. SMOOTHNESS CONSTRAINTS
Because the convolution of factor tensors can represent any tensor (even rank 1 models), it is necessary to impose constraints to narrow down the candidate solutions. In this study, smoothness is used as a constraint. The reasons for introducing smoothness as a constraint are as follows.
• As shown in Figure 2, the embedded data is represented by a smooth manifold on the delay embedded space.
• The data mainly targeted in our study are images, and there are many reports that smooth constraints are effective in image completion [48], [15], [53], [47].
Note that we do not introduce smoothness for the reconstructed tensor but the factor tensors. Unlike the model which smoothens the reconstructed tensor, the proposed model enables completion without excessive smoothing. We also set the scale adjustment terms for both A and B in the optimization equation to avoid smoothing by increasing only one factor of the tensors.

C. OPTIMIZATION FORMULAS AND ALGORITHM 1) OPTIMIZATION FORMULAS
In this paper, we propose a new tensor completion model. We assume that the observed tensor Y ∈ R T 1 ×···×T N is incomplete and that some entries have no values. The projection tensor O ∈ {0, 1} T 1 ×···×T N passes the observed entries and makes the missing entries equal to zero. The entries are given by The problem involves obtaining the complete tensor A * B.
In this study, we impose a smoothness constraint on A and B. The optimization problem is then given by where is a differential filter, and A ∈ R T 1 ×···×T N , B ∈ R τ 1 ×···×τ N are factor tensors and fold V : is a folding operator from a vector to the N -th order tensor. Equation (22) evaluates the reconstruction loss in the first term. The second and third terms are smooth penalties for A and B, and the fourth and fifth terms adjust the scales of A and B. The equality constraint is for zero padding, based on Equation (12). Note that when τ = 1, Equation (22) is equivalent to QV regularization. The relaxation of optimization problem (22) for an unconstrained optimization problem including a penalty term yields the following equation: where ] ∈ R T n . i τ n serves as a penalty for B and simulates zero padding {P}. Note that we also redefine the size of B as T 1 × · · · × T N .

2) ALGORITHM FOR SOLVING OPTIMIZATION
In this study, we solved the optimization problem (23) using the majorization-minimization (MM) algorithm [18], [30]. The MM algorithm is an iterative method involving two steps. 1) Constructs a auxiliary function h(A, B|A (k) , B (k) ) for L(A, B) at A (k) , B (k) . Note,

2) Update as in
3) Update as in A conceptual diagram of the algorithm is shown in Figure 5. The MM algorithm was used because of convergence due to its monotonic convergence and ease of analytical computation.
The auxiliary function h is defined as follows: UpdateB by (32) 10: A ← IFFT(Â) 11: B ← IFFT(B) 12: Calculate L (A, B). 13: if convergence of L then break 14: end if 15: end for where and Furthermore, we exploit the property that the cyclic convolution of the time domain is a Hadamard product in the frequency domain to reduce the time complexity of the optimization problem. Because the Frobenius norm is invariant to the Fourier transform, the auxiliary function h is redefined as whereL n ,Â,B,Ẑ,Ŵ are the Fourier transform of L n , A, B, Z, W.
The final auxiliary function isĥ, which is minimized using Equations (24) and (25). To minimizeĥ, we deriveÂ and B such that ∂ĥ 67532 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.   Note that (30) is substituted by alternating the optimization with (24) and (25) because they are not satisfied simultaneously. Thus, at every optimization step, although the cost function L decreases monotonically, the auxiliary function h is not always optimal. After solving Equation (30), we obtain whereẐ * is the complex conjugate ofẐ. In summary, Equations (27) and (28)

E. RELATION TO OUR EXISTING PAPER
This paper is an extended version based on our paper [38]. Therefore, the concept of the proposed data completion model is almost identical. However, existing methods discuss only vector data (one-dimensional data) and do not target tensor data (multidimensional data). The proposed method differs significantly from existing methods in extending the model design and algorithms to multiple dimensions. Furthermore, in addition to the 1D signal completion of the existing method, this paper also experiments with image completion for RGB and MRI images (see VOLUME 11, 2023   Sections V-B1 and V-B2) and hyper-parameter sensitivities experiments (see Section V-C2).

V. EXPERIMENT A. COMPLETION OF CLIPPED DATA
This section presents the evaluation of the proposed method (SCTF) for reconstructing clipped data (declipping) as a type of completion. Clipping is an operation that uses a certain clipping level c > 0 to replace entries above c and below −c with c and −c. The clipping operation on the entry of a tensor is given by X (t 1 , t 2 , · · · , t N ) = min(c, max(−c, X 0 (t 1 , · · · , t N ))). (33) The value range of the clipped data was [−c, c].
The indices of the clipped entries were recorded and treated as missing values for reconstruction. Thus, the set of observed entries can be expressed as

1) CLIPPED SIGNAL COMPLETION (1ST ORDER tensor)
In this experiment, we evaluated SCTF using signal completion. The signals to be completed were the sine and wavelet functions with a clip level of 0.2, and both had a maximum amplitude of 1. Examples of clipping are presented in Figures 6 and 7. The original signals are shown in Figures 6a and 7a, and the signals after clipping are shown in Figures 6b and 7b. We compared SCTF with QV regularization and spline interpolation. In SCTF, we set τ 1 = 32, γ = 1.0 × 10 8 , and λ A,1 = λ B,1 = η A = η B = 0.1. Figures 6c and 7c show the signals reconstructed using QV regularization, spline interpolation, and SCTF. The signal completed by SCTF had a maximum amplitude of approximately 1, which was the best among the three methods. The spline method reconstructs the waveform signals; however, its amplitude is smaller than that of the original signal. However, the QV method failed to restore the signal. Figure 8 shows the signal-to-noise ratio (SNR) values of the declipping experiment for various clipping levels. The clipping levels were 0.8, 0.6, 0.4, and 0.2, and the parameters (λ, τ ) were adjusted for all levels. SCTF achieved significantly higher values of SNR than QV regularization and spline interpolation.

2) COMPLETION OF CLIPPED DATA (2ND ORDER tensor)
In this experiment, we evaluated SCTF by completing a clipped image of 2D-sin. The maximum value of the image (amplitude of 2D-sin) was 255, and the clip threshold c = 230. The original image is shown in Figure 9a, and the image after clipping is shown in Figure 9b. We compared SCTF with the QV model and Fast-MDT-Tucker [45]. The QV model is the model without convolution in Equation (22), and is     In addition, Table 1 shows the numerical evaluation of the completion accuracy. The completion image by SCTF has a maximum amplitude close to 255 and shows an improvement in the oscillations that occurred in Fast-MDT-Tucker, indicating the effect of the smoothing term. In fact, SCTF had the best value in PSNR. However, SCTF had a flat shape that differed from that of sine at a maximum amplitude of 255. This affects the numerical evaluation of the completion accuracy, and SCTF is worse in SSIM than in Fast-MDT-Tucker. On the other hand, the QV model fails to recover 2D-sin because it flatly completes the missing parts.

B. COMPLETION OF RANDOM MISSING DATA
This section presents the results of the completion of random missing data.

C. ANALYSIS OF ALGORITHM 1) CONVERGENCE OF ALGORITHM
Monotonic convergence is expected because the proposed algorithm uses the MM algorithm. Figure 13 shows that the objective function converges monotonically and that the algorithm works correctly.
The experimental setup was the same as that in Section V-A2; we recovered the clipped 2D-sin image. Three hyper-parameters were varied in the range τ ∈ {1, 9, 25, 81, 151, 315} and λ ∈ {0, 200, 400, 600, 800, 1000} and η ∈ {0, 200, 400, 600 , 800, 1000}. A declipping experiment was performed for all the combinations to calculate the recovery accuracy. Figure 14 shows the experiment's five-time average of the PSNR. Increasing τ improves the accuracy, whereas making it too large worsens the accuracy. For example, when τ = 1, the algorithm matches the QV regularization and recovers smoothly; however, when the delay window is smaller than the clip range, such as when τ = 9 or τ = 25, it cannot recover at all. Therefore, it is important to set τ appropriately.

VI. CONCLUSION
In this study, we proposed a new model and algorithm for tensor completion using a convolution of smooth-factor tensors. Because the proposed method corresponds to a rank-1 decomposition in the delay-embedded space, it can achieve high completion accuracy in a short computation time. In the optimization formulation, we extended the existing mathematical model based on the inverse MDT by adding a penalty term for the factor tensor corresponding to the delay-embedding width. In addition, we set smoothing constraints for the factor to narrow down the candidate solutions. As for the algorithm for solving the optimization, we employed the MM algorithm with the expectation of monotonic convergence. Our experiments mainly completed clipped and random missing image data and confirmed that the proposed method achieves high completion accuracy with low computational cost. In the experiment, we also confirmed the effect of the completion accuracy on the variation in the delay-embedding width and monotone convergence of the algorithm.