Rank-Adaptive Tensor Completion Based on Tucker Decomposition

Abstract

Tensor completion is a fundamental tool to estimate unknown information from observed data, which is widely used in many areas, including image and video recovery, traffic data completion and the multi-input multi-output problems in information theory. Based on Tucker decomposition, this paper proposes a new algorithm to complete tensors with missing data. In decomposition-based tensor completion methods, underestimation or overestimation of tensor ranks can lead to inaccurate results. To tackle this problem, we design an alternative iterating method that breaks the original problem into several matrix completion subproblems and adaptively adjusts the multilinear rank of the model during optimization procedures. Through numerical experiments on synthetic data and authentic images, we show that the proposed method can effectively estimate the tensor ranks and predict the missing entries.

1. Introduction

Tensors, as a higher-order generalization of vectors and matrices, can preserve the original structure and the latent characteristics of multi-dimensional data, such as images and videos. To store and process these data, vectorization or matricization may break the adjacency relation and lead to the loss of crucial information in practical applications. Therefore, tensor analysis of data has gathered increasing attention due to its better performance in finding hidden structures and capturing potential features. Due to data missing from a transmission or insufficient collection, sometimes we have to deal with incomplete data, which increases the operation cost and the difficulty of multi-dimensional data analysis. Therefore, tensor completion, also referred to as tensor recovery, plays a key role in processing and analyzing incomplete multi-dimensional data. In tensor completion, the prediction of unknown data from a few observed entries is achieved based on the correlation between different parts of data or, simply, the hidden low-rank structure. Many real-world data lie in a latent low dimensional space since they may only relate to a small number of contributions; thus, we can use tensor completion to infer the missing data. It has already been widely applied in many scientific fields, including surrogate construction [1], image and video recovery [2,3,4], traffic data completion [5] and the MIMO (multi-input multi-output) problems in information theory [6].

The tensor completion methods can be roughly classified into model-based and non-model-based ones. Extending the idea from existing matrix completion algorithms [7,8], non-model-based methods directly formulate the tensor completion task as a rank minimization problem subject to linear constraints, and then relax it to a nuclear norm minimization problem, such as LRTC (low rank tensor completion) [9]. However, model-based methods achieve the goal differently. Based on a decomposition model with a given rank, these methods minimize the error on the known entries by adjusting the model parameters, such as CP-Wopt (CP-weighted optimization) [10], Tucker–Wopt (Tucker-weighted optimization) [5], ${\mathrm{M}}^2\mathrm{SA}$ (multilinear subspace analysis with missing values) [11] and so on.

The computation cost of non-model-based methods is closely related to the size of the tensor, leaving little room for reducing the time complexity. On the contrary, this can be more flexibly controlled in model-based methods by altering the model size. Nonetheless, there still exists a big challenge when choosing the model’s rank. Underestimation of the tensor rank makes it harder to reach our desired accuracy, whereas overestimation may lead to a poor generalization on the unobserved part of entries. Therefore, rank-adaptive methods are needed and worth studying. In [1], a greedy scheme is designed to find the fittest rank by gradually increasing the number of factors in the CP decomposition model. In [12], from an overestimated rank, a systematic approach to reducing the CP ranks to optimal ones is developed.

Focusing on tensor completion based on Tucker decomposition, we propose a novel rank-adaptive tensor completion method and verify its efficiency through experiments on synthetic and real practical data. The rest of the paper is organized as follows. Section 2 formulates the problem setting. Section 3 presents Tucker decomposition and the related works. We propose our rank-adaptive tensor completion based on the Tucker decomposition (RATC-TD) approach in Section 4 and verify it through numerical experiments in Section 5. Finally, we conclude our work in Section 6.

2. Notations and the Tensor Completion Problem

Before summarizing the related work and presenting our algorithm, this section first introduces some basic definitions that will be used in the following sections, and then states our problem settings.

2.1. Notations

We follow the notations in [13]. Let lowercase letters denote scalars, e.g., a; let lowercase bold letters denote vectors, e.g., $\mathbf{a}$; let bold capital letter denote matrices, e.g., $\mathit{U}$; let capital fraktur letter denote high-order tensors, e.g., $\mathcal{T}$. The ith element of a vector $\mathbf{a}$ is denoted as $a_i$; the element of matrix $\mathit{U}$ with index $(i,j)$ is denoted by $u_{ij}$; the element of N-way tensor $\mathcal{T}$ with index $(i_1,i_2,\cdots ,i_N)$ is denoted by $t_{i_1{i}_2\cdots i_N}$. A tensor can be viewed as a multi-dimensional array, and the number of dimensions is generally called the order of the tensor. For example, matrices are second-order tensors, while vectors are first-order tensors. A mode-n fiber is a vector extracted from a tensor by fixing all indices except the nth one. Consider a matrix as an example; mode-1 fibers refer to the columns, and mode-2 fibers refer to the rows. The mode-n unfolding of a tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ is to rearrange the mode-n fibers as the columns of the resulting matrix, denoted as ${\mathit{T}}_{\left(n\right)}\in {\mathbb{R}}^{I_n\times I_1\cdots I_{n-1}{I}_{n+1}\cdots I_N}$. Naturally, the reverse of unfolding, the process of rearranging ${\mathit{T}}_{\left(n\right)}\in {\mathbb{R}}^{I_n\times I_1\cdots I_{n-1}{I}_{n+1}\cdots I_N}$ into $\mathcal{T}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$, is called folding.

The n-mode matrix product is the multiplication of a tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ with a matrix $\mathit{U}\in {\mathbb{R}}^{J\times I_n}$, denoted as $\mathcal{X}=\mathcal{T}{\times }_n\mathit{U}$, where $\mathcal{X}\in {\mathbb{R}}^{I_1\times \cdots \times I_{n-1}\times J\times I_{n+1}\times \cdots \times I_N}$. It can be regarded as multiplying each mode-n fiber of $\mathcal{T}$ with $\mathit{U}$. The matricized form is given by

$${\mathit{X}}_{\left(n\right)}=\mathit{U}{\mathit{T}}_{\left(n\right)},$$

(1)

and the element-wise operations can be obtained by

$$x_{i_1\cdots i_{n-1}j{i}_{n+1}\cdots i_N}=\sum _{i_n=1}^{I_n}{t}_{i_1\cdots i_{n-1}{i}_n{i}_{n+1}\cdots i_N}{u}_{j{i}_n}.$$

(2)

The norm of a tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ is a higher-order analog of the matrix Frobenius norm, defined by

$${\parallel \mathcal{T}\parallel }_F=\sqrt{\sum _{i_1=1}^{I_1}\sum _{i_2=1}^{I_2}\cdots \sum _{i_N=1}^{I_N}{t}_{i_1{i}_2\cdots i_N}^2}.$$

(3)

The inner product of two same sized tensors $\mathcal{T},\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ is defined by

$$\langle \mathcal{T},\mathcal{X}\rangle =\sum _{i_1=1}^{I_1}\sum _{i_2=1}^{I_2}\cdots \sum _{i_N=1}^{I_N}{t}_{i_1{i}_2\cdots i_N}{x}_{i_1{i}_2\cdots i_N}.$$

(4)

2.2. The Tensor Completion Problem

Tensor completion is the problem of predicting missing entries through a few sampled entries of a tensor based on the assumption that it has a low-rank structure. Supposing we have a partially observed tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times \cdots \times I_N}$, the set of indices where entries are observed is denoted as $\mathsf{\Omega }$. Let ${\mathcal{P}}_{\mathsf{\Omega }}$ be an operator on ${\mathbb{R}}^{I_1\times \cdots \times I_N}$, such that ${\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{X}\right)$ is equal to $\mathcal{X}$ on $\mathsf{\Omega }$ and is equal to 0 otherwise. The low-rank tensor completion problem can be formulated as

$$\begin{array}{cc}\hfill \underset{\mathcal{X}}{\min }& \mathrm{rank}\left(\mathcal{X}\right)s.t.{\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{X}\right)={\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right).\hfill \end{array}$$

(5)

Although the formulation seems easy, how to measure the rank of a tensor is still an open question. In previous research, the CP rank [1,12] and the Tucker rank [9,14,15] were mostly considered. After the tensor-train (TT) decomposition was proposed in [16], tensor completion based on the TT rank was studied in [3,17,18]. The low-tubal rank tensor completion problem was studied in [19,20,21]. In addition, the combination of the CP and the Tucker ranks was investigated in [22].

Since matrices are second-order tensors, matrix completion can be viewed as a special case of tensor completion. Using the matrix nuclear norm to approximate the rank of a matrix is studied for matrix completion in [23]. Liu et al. [9] generalized this to the tensor scheme by defining the nuclear norm of tensors and proposing the following low-rank tensor completion problem:

$$\begin{array}{cc}\hfill \underset{\mathcal{X}}{\min }& {\parallel \mathcal{X}\parallel }_{*}s.t.{\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{X}\right)={\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right).\hfill \end{array}$$

(6)

The nuclear norm of a tensor is defined by ${\parallel \mathcal{X}\parallel }_{*}:={\sum }_{n=1}^N{\parallel {\mathit{X}}_{\left(n\right)}\parallel }_{*}$, where $\parallel {\mathit{X}}_{\left(n\right)}{\parallel }_{*}$ is the nuclear norm of the matrix ${\mathit{X}}_{\left(n\right)}$. The nuclear norm can be computed by the sum of the singular values. However, the computation involves SVD (singular value decomposition) of the unfolding matrix ${\mathit{T}}_{\left(n\right)}$ in each iteration, which is expensive when the size of $\mathcal{T}$ is large. To control the computation cost, this problem is reformulated based on tensor decomposition [5,11]. Particularly, the low-rank tensor completion problem based on Tucker decomposition can be written as

$$\begin{array}{cc}\hfill \underset{\mathcal{G},{\mathit{A}}_1\cdots ,{\mathit{A}}_N}{\min }& {\parallel \mathcal{T}-\mathcal{X}\parallel }_{\mathsf{\Omega }}^2\hfill \\ \hfill s.t.\mathcal{X}=\mathcal{G}{\times }_1{\mathit{A}}_1\cdots {\times }_N{\mathit{A}}_N,& {\mathit{A}}_n^T{\mathit{A}}_n=I^{R_n\times R_n},n=1,\cdots ,N,\hfill \end{array}$$

(7)

where $\mathcal{X}$ is the completed goal tensor, $(R_1,\cdots ,R_N)$ is a predefined multilinear rank, ${\mathit{A}}_n\in {\mathbb{R}}^{I_n\times R_n}$ and ${\parallel ·\parallel }_{\mathsf{\Omega }}$ denote the norm on the observation elements, i.e., $\parallel {\mathcal{P}}_{\mathsf{\Omega }}{(·)\parallel }_F$. Methods are actively developed to solve this problem, including ${\mathrm{M}}^2\mathrm{SA}$ [11], gHOI (generalized higher-order orthogonal iteration) [24] and Tucker–Wopt [5]. However, these algorithms typically require the multilinear ranks given a priori, both underestimating and overestimating the ranks can result in inaccurate completion results. To handle this issue, we in this work propose a rank-adaptive tensor completion approach without requiring the multilinear ranks given a priori.

3. The Tucker Decomposition and Its Computation

We review Tucker decomposition and its related algorithms in this section.

3.1. Tucker Decomposition

Tucker decomposition was first introduced by Tucker [25]. For a Nth-order tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$, its Tucker decomposition has the form

$$\begin{array}{c}\hfill \mathcal{T}=\mathcal{G}{\times }_1{\mathit{A}}_1{\times }_2{\mathit{A}}_2\cdots {\times }_N{\mathit{A}}_N,\end{array}$$

(8)

where $\mathcal{G}\in {\mathbb{R}}^{R_1\times \cdots \times R_N}$ is a core tensor and ${\mathit{A}}_n\in {\mathbb{R}}^{I_n\times R_n}$, $n=1,2,\cdots ,N$, are factor matrices. For simplicity, we assume all the factor matrices have orthogonal columns. $R_n$ is the rank of the mode-n unfolding matrix ${\mathit{T}}_{\left(n\right)}$. As described in [13], Tucker decomposition can be viewed as a higher-order generalization of PCA (principal component analysis) in some sense. In Tucker decomposition, the columns of the factor matrix represent the components in each mode, and each element of the core tensor characterizes the level of interaction between the components in different modes.

The tuple $(R_1,\cdots ,R_N)$ is called Tucker rank or multilinear rank. However, practically, we prefer to use the truncated version of Tucker decomposition, where we set $R_n<\mathrm{rank}\left({\mathit{T}}_{\left(n\right)}\right)$ for one or more n. Given the truncation $(R_1,\cdots ,R_N)$, the approximation of the truncated Tucker decomposition of a tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times \cdots \times I_N}$ can be described as

$$\begin{array}{cc}\hfill \underset{\mathcal{G},{\mathit{A}}_1\cdots ,{\mathit{A}}_N}{\min }& {\parallel \mathcal{T}-\mathcal{X}\parallel }_F^2\hfill \\ \hfill s.t.\mathcal{X}=\mathcal{G}{\times }_1{\mathit{A}}_1\cdots {\times }_N{\mathit{A}}_N,& {\mathit{A}}_n^T{\mathit{A}}_n=I^{R_n\times R_n},n=1,\cdots ,N,\hfill \end{array}$$

(9)

where $\mathcal{G}\in {\mathbb{R}}^{R_1\times \cdots \times R_N}$, ${\mathit{A}}_n\in {\mathbb{R}}^{I_n\times R_n}$ is the nth factor matrix, and $I^{R_n\times R_n}$ is an identity matrix of size $R_n\times R_n$.

3.2. The Higher-Order Orthogonal Iteration (HOOI) Algorithm

There are several approaches to computing the truncated Tucker decomposition of a tensor. One of the most popular methods is the higher-order orthogonal iteration (HOOI) algorithm [26], also referred to as the alternative least square (ALS) algorithm for Tucker decomposition. It is an alternative iterating method, where we fix all except one factor matrix each time, and then minimize the objective function in (9). Specifically, given initial guesses $\{{\mathit{A}}_n^{\left(0\right)}:n=1,\cdots ,N\}$, at the kth iteration, for each $n=1,\cdots ,N$, we fix all the factor matrices except ${\mathit{A}}_n$, and then find the optimal solution to the subproblem

$$\begin{array}{c}\hfill \underset{\mathcal{G},{\mathit{A}}_n^{\left(k\right)}}{\min }{\parallel \mathcal{T}-\mathcal{X}\parallel }_F^2\\ \hfill s.t.\mathcal{X}=\mathcal{G}{\times }_1{\mathit{A}}_1^{\left(k\right)}\cdots {\times }_n{\mathit{A}}_n^{\left(k\right)}{\times }_{n+1}{\mathit{A}}_{n+1}^{(k-1)}\cdots & {\times }_N{\mathit{A}}_N^{(k-1)},{\left({\mathit{A}}_n^{\left(k\right)}\right)}^T{\mathit{A}}_n^{\left(k\right)}=I^{R_n\times R_n},\hfill \end{array}$$

(10)

where $\mathcal{G}\in {\mathbb{R}}^{R_1\times \cdots \times R_N}$, ${\mathit{A}}_n^{\left(k\right)}\in {\mathbb{R}}^{I_n\times R_n}$. Denoting $\mathcal{B}:=\mathcal{G}{\times }_n{\mathit{A}}_n^{\left(k\right)}$, this optimization problem can be simplified as

$$\begin{array}{c}\hfill \underset{\mathcal{B}}{\min }{\parallel \mathcal{T}-\mathcal{X}\parallel }_F^2\\ \hfill s.t.\mathcal{X}=\mathcal{B}{\times }_1{\mathit{A}}_1^{\left(k\right)}\cdots {\times }_{n-1}{\mathit{A}}_{n-1}^{\left(k\right)}{\times }_{n+1}{\mathit{A}}_{n+1}^{(k-1)}\cdots & {\times }_N{\mathit{A}}_N^{(k-1)},\mathrm{rank}\left({\mathit{B}}_{\left(n\right)}\right)\le R_n,\hfill \end{array}$$

(11)

where $\mathcal{B}\in {\mathbb{R}}^{R_1\times \cdots R_{n-1}\times I_n\times R_{n+1}\cdots R_N}$ and ${\mathit{B}}_{\left(n\right)}$ is the mode-n unfolding matrix form of $\mathcal{B}$. It can be considered as a constrained least squares problem [27], and thus the solution to (10) can be easily obtained by the following steps:

• Compute ${\mathcal{B}}^{*}=\mathcal{T}{\times }_1{\left({\mathit{A}}_1^{\left(k\right)}\right)}^T\cdots {\times }_{n-1}{\left({\mathit{A}}_{n-1}^{\left(k\right)}\right)}^T{\times }_{n+1}{\left({\mathit{A}}_{n+1}^{(k-1)}\right)}^T\cdots {\times }_N{\left({\mathit{A}}_N^{(k-1)}\right)}^T$;
• Unfold ${\mathcal{B}}^{*}$ at the nth mode to obtain ${\mathit{B}}_{\left(n\right)}^{*}$, then perform a truncated rank-$R_n$ singular value decomposition ${\left({\mathit{B}}_{\left(n\right)}^{*}\right)}_{\left[R\right]}=\mathit{U}\mathsf{\Sigma }{\mathit{V}}^T={\mathit{B}}_{\left(n\right)}$;
• Compute ${\mathit{A}}_n^{\left(k\right)}=\mathit{U}$ and ${\mathbf{G}}_{\left(n\right)}=\mathsf{\Sigma }{\mathit{V}}^T$.

We iteratively solve these subproblems until a given stopping criterion is reached, i.e., the difference between the solutions at two adjacent iterations is small enough, or the value of the objective function decreases very slightly. The complete procedure of HOOI is presented in Algorithm 1.

Algorithm 1 The high-order orthogonal iteration (HOOI) algorithm [26]

Input: Tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ and truncation $(R_1,R_2,\cdots ,R_N)$. Output: Core tensor $\mathcal{G}\in {\mathbb{R}}^{R_1\times \cdots \times R_N}$, and factor matrices ${\mathit{A}}_n\in {\mathbb{R}}^{I_n\times R_n}$ for $n=1,2,\cdots ,N$. 1: Initialize ${\mathit{A}}_n^{\left(0\right)}\in {\mathbb{R}}^{I_n\times R_n}$ for $n=1,2,\cdots ,N$ using HOSVD. 2: $k\leftarrow 0$. 3: while not converge do 4:     $k\leftarrow k+1$. 5:     for $n=1,2,\cdots ,N$ do 6:         $\mathcal{B}\leftarrow \mathcal{T}{\times }_1{\left({\mathit{A}}_1^{\left(k\right)}\right)}^T\cdots {\times }_{n-1}{\left({\mathit{A}}_{n-1}^{\left(k\right)}\right)}^T{\times }_{n+1}{\left({\mathit{A}}_{n+1}^{(k-1)}\right)}^T\cdots {\times }_N{\left({\mathit{A}}_N^{(k-1)}\right)}^T$. 7:         ${\mathit{B}}_{\left(n\right)}\leftarrow$ mode-n unfolding matrix of $\mathcal{B}$. 8:         $\mathit{U},\mathsf{\Sigma },{\mathit{V}}^T\leftarrow$ truncated rank-$R_n$ SVD of ${\mathit{B}}_{\left(n\right)}$. 9:         ${\mathit{A}}_n^{\left(k\right)}\leftarrow \mathit{U}$. 10:     end for 11: end while 12: $\mathcal{G}\leftarrow \mathsf{\Sigma }{\mathit{V}}^T$ folding at mode-n.

3.3. The Rank-Adaptive HOOI Algorithm

HOOI requires the multilinear rank $(R_1,\cdots ,R_N)$ given a priori, which is hard to be determined in practice. Instead of (9), Xiao and Yang [27] consider the following form of the low multilinear rank approximation problem:

$$\begin{array}{c}\hfill \underset{\mathcal{G},{\mathit{A}}_1\cdots ,{\mathit{A}}_N}{\min }(R_1,R_2,\cdots ,R_N){s.t.\parallel \mathcal{X}-\mathcal{T}\parallel }_F^2<ϵ{\parallel \mathcal{T}\parallel }_F^2,\\ \hfill \mathcal{X}=\mathcal{G}{\times }_1{\mathit{A}}_1\cdots {\times }_N{\mathit{A}}_N,{\mathit{A}}_n^T{\mathit{A}}_n=I^{R_n\times R_n},n=1,\cdots ,N,\end{array}$$

(12)

where $ϵ$ is a given tolerance. A rank-adaptive HOOI algorithm is proposed in [27], which adjusts the truncation $R_n$ for dimension n in the HOOI iterations by

$$\begin{array}{c}\hfill R_n^{\left(k\right)}=\underset{R}{\mathrm{argmin}}\parallel {\mathit{B}}_{\left(n\right)}-{\left({\mathit{B}}_{\left(n\right)}\right)}_{\left[R\right]}{\parallel }_F^2\le {\parallel \mathcal{B}\parallel }_F^2-(1-ϵ){\parallel \mathcal{T}\parallel }_F^2,\end{array}$$

(13)

where ${\left({\mathit{B}}_{\left(n\right)}\right)}_{\left[R\right]}$ is the best rank-R approximation of ${\mathit{B}}_{\left(n\right)}$. For the full SVD of ${\mathit{B}}_{\left(n\right)}=\mathit{U}\mathsf{\Sigma }{\mathit{V}}^T$, it can be calculated by ${\left({\mathit{B}}_{\left(n\right)}\right)}_{\left[R\right]}={\mathit{U}}_{:,1:R}{\mathsf{\Sigma }}_{1:R,1:R}{\mathit{V}}_{:,1:R}^T$. In [27], it is proven that (13) is a local optimal strategy for updating $R_n$, i.e., the optimal solution of the following problem:

$$\begin{array}{c}\hfill \underset{\mathcal{B}}{\min }\mathrm{rank}\left({\mathit{B}}_{\left(n\right)}\right){s.t.\parallel \mathcal{X}-\mathcal{T}\parallel }_F^2<ϵ{\parallel \mathcal{T}\parallel }_F^2,\\ \hfill \mathcal{X}=\mathcal{B}{\times }_1{\mathit{A}}_1^{\left(k\right)}\cdots {\times }_{n-1}{\mathit{A}}_{n-1}^{\left(k\right)}{\times }_{n+1}{\mathit{A}}_{n+1}^{(k-1)}\cdots {\times }_N{\mathit{A}}_N^{(k-1)}.\end{array}$$

(14)

Details of the above procedure are summarized in Algorithm 2.

Algorithm 2 The rank-adaptive HOOI algorithm [27]

Input: Tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$, error tolerance $ϵ$, initial guess of factor matrices ${\mathit{A}}_n^{\left(0\right)}\in {\mathbb{R}}^{I_n\times R_n^{\left(0\right)}}$ for $n=1,2,\cdots ,N$ and initial truncation $(R_1^{\left(0\right)},R_2^{\left(0\right)},\cdots ,R_N^{\left(0\right)})$. Output: Truncation $(R_1^{\left(k\right)},R_2^{\left(k\right)},\cdots ,R_N^{\left(k\right)})$, core tensor $\mathcal{G}\in {\mathbb{R}}^{R_1^{\left(k\right)}\times \cdots \times R_N^{\left(k\right)}}$ and factor matrices ${\mathit{A}}_n^{\left(k\right)}\in {\mathbb{R}}^{I_n\times R_n^{\left(k\right)}}$ for $n=1,2,\cdots ,N$. 1: ${\mathcal{G}}^{\left(0\right)}\leftarrow \mathcal{T}{\times }_1{\left({\mathit{A}}_1^{\left(0\right)}\right)}^T\cdots {\times }_N{\left({\mathit{A}}_N^{\left(0\right)}\right)}^T$. 2: $k\leftarrow 0$. 3: while not converge do 4:     $k\leftarrow k+1$. 5:     for all $n\in \{1,2,\cdots ,N\}$ do 6:         $\mathcal{B}\leftarrow \mathcal{T}{\times }_1{\left({\mathit{A}}_1^{\left(k\right)}\right)}^T\cdots {\times }_{n-1}{\left({\mathit{A}}_{n-1}^{\left(k\right)}\right)}^T{\times }_{n+1}{\left({\mathit{A}}_{n+1}^{(k-1)}\right)}^T\cdots {\times }_N{\left({\mathit{A}}_N^{(k-1)}\right)}^T$. 7:         ${\mathit{B}}_{\left(n\right)}\leftarrow$ mode-n unfolding matrix of $\mathcal{B}$. 8:         $\mathit{U},\mathsf{\Sigma },{\mathit{V}}^T\leftarrow$ full SVD of ${\mathit{B}}_{\left(n\right)}$. 9:         $R_n^{\left(k\right)}\leftarrow$ minimum R such that ${\displaystyle \sum _{r>R}}{\mathsf{\Sigma }}_{r,r}^2{<\parallel \mathcal{B}\parallel }_F^2-(1-ϵ){\parallel \mathcal{T}\parallel }_F^2$. 10:         ${\mathit{A}}_n^{\left(k\right)}\leftarrow {\mathit{U}}_{:,1:R_n^{\left(k\right)}}$. 11:     end for 12:     $\mathcal{G}\leftarrow {\mathsf{\Sigma }}_{1:R_n^{\left(k\right)},1:R_n^{\left(k\right)}}{\mathit{V}}_{:,1:R_n^{\left(k\right)}}^T$ folding at mode-n. 13: end while

4. A Rank-Adaptive Tensor Recovery Scheme

This section proposes a rank-adaptive tensor completion scheme based on truncated Tucker decomposition (RATC-TD).

Analogous to the rank-adaptive HOOI algorithm, instead of (7), we consider a different form of the low-rank tensor completion problem:

$$\begin{array}{c}\hfill \underset{\mathcal{G},{\mathit{A}}_1\cdots ,{\mathit{A}}_N}{\min }(R_1,R_2,\cdots ,R_N)s.t.{\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{X}\right)={\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right),\\ \hfill \mathcal{X}=\mathcal{G}{\times }_1{\mathit{A}}_1\cdots {\times }_N{\mathit{A}}_N,{\mathit{A}}_n^T{\mathit{A}}_n=I^{R_n\times R_n},n=1,\cdots ,N.\end{array}$$

(15)

Note that if the data are noisy, the constraint must be relaxed to ${\parallel \mathcal{X}-\mathcal{T}\parallel }_{\mathsf{\Omega }}^2<ϵ{\parallel \mathcal{T}\parallel }_{\mathsf{\Omega }}^2$, where $ϵ$ is a given tolerance of the relative error on the observation part between the original tensor and the completed one.

Using an alternative optimization technique, problem (15) can be divided into N subproblems, which will be detailedly represented in Section 4.1. To address those subproblems, the singular value thresholding (SVT) algorithm [7] is introduced in Section 4.2. The entire algorithm is summarized in Section 4.3.

4.1. Alternative Optimization in Tensor Completion

We solve (15) by an alternative iteration method. For each $n=1,\cdots ,N$, initialize $R_n^{\left(0\right)}=I_n$, and the initial guess ${\mathit{A}}_n^{\left(0\right)}$ can be either randomly set or given by HOSVD (higher-order SVD) [13] of ${\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right)$. At the kth iteration, for each $n=1,\cdots ,N$, by fixing all the factor matrices except ${\mathit{A}}_n$, we solve the subproblem

$$\begin{array}{cc}\hfill \underset{\mathcal{G},{\mathit{A}}_n^{\left(k\right)}}{\min }\left(R_n^{\left(k\right)}\right)& s.t.{\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{X}\right)={\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right),\hfill \\ \hfill \mathcal{X}=\mathcal{G}{\times }_1{\mathit{A}}_1^{\left(k\right)}\cdots {\times }_n{\mathit{A}}_n^{\left(k\right)}{\times }_{n+1}{\mathit{A}}_{n+1}^{(k-1)}& \cdots {\times }_N{\mathit{A}}_N^{(k-1)},{\left({\mathit{A}}_n^{\left(k\right)}\right)}^T{\mathit{A}}_n^{\left(k\right)}=I^{R_n\times R_n}.\hfill \end{array}$$

(16)

Denoting $\mathcal{B}:=\mathcal{G}{\times }_n{\mathit{A}}_n^{\left(k\right)}$, (16) is equivalent to

$$\begin{array}{c}\hfill \underset{\mathcal{B}}{\min }\mathrm{rank}\left({\mathit{B}}_{\left(n\right)}\right)s.t.{\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{X}\right)={\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right),\\ \hfill \mathcal{X}=\mathcal{B}{\times }_1{\mathit{A}}_1^{\left(k\right)}\cdots {\times }_{n-1}{\mathit{A}}_{n-1}^{\left(k\right)}{\times }_{n+1}{\mathit{A}}_{n+1}^{(k-1)}\cdots {\times }_N{\mathit{A}}_N^{(k-1)},\end{array}$$

(17)

where ${\mathit{B}}_{\left(n\right)}$ is the mode-n unfolding matrix form of $\mathcal{B}$. However, this rank minimization problem is NP-hard. Its tightest convex relaxation is

$$\begin{array}{c}\hfill \underset{\mathcal{B}}{\min }{\parallel {\mathit{B}}_{\left(n\right)}\parallel }_{*}s.t.{\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{X}\right)={\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right),\\ \hfill \mathcal{X}=\mathcal{B}{\times }_1{\mathit{A}}_1^{\left(k\right)}\cdots {\times }_{n-1}{\mathit{A}}_{n-1}^{\left(k\right)}{\times }_{n+1}{\mathit{A}}_{n+1}^{(k-1)}\cdots {\times }_N{\mathit{A}}_N^{(k-1)},\end{array}$$

(18)

where $\parallel {\mathit{B}}_{\left(n\right)}{\parallel }_{*}$ is the nuclear norm of the matrix ${\mathit{B}}_{\left(n\right)}$, which can be computed by the sum of the singular values. The matrix form of this problem can be written as [27]

$$\begin{array}{c}\hfill \underset{\mathit{B}}{\min }{\parallel \mathit{B}\parallel }_{*}s.t.{\mathit{P}}_{\mathsf{\Omega }}\left(\mathit{B}{\mathit{M}}^T\right)={\mathit{P}}_{\mathsf{\Omega }}\left({\mathit{T}}_{\left(n\right)}\right),\\ \hfill \mathit{M}={\mathit{A}}_N\cdots \otimes {\mathit{A}}_{n+1}\otimes {\mathit{A}}_{n-1}\cdots \otimes {\mathit{A}}_1,\end{array}$$

(19)

where ${\mathit{T}}_{\left(n\right)}$ is the mode-n unfolding matrix form of $\mathcal{T}$, and ⊗ is the Kronecker product. In this work, we use the singular value thresholding (SVT) algorithm [7] to obtain a proximity solution to (19).

4.2. Solving the Subproblems Using SVT

Following the notation used in [7], we first introduce the singular value thresholding operator. Consider the skinny singular value decomposition of a matrix $\mathit{X}\in {\mathbb{R}}^{n_1\times n_2}$ of rank r,

$$\mathit{X}=\mathit{U}\mathsf{\Sigma }{\mathit{V}}^T,\mathsf{\Sigma }=\mathrm{diag}\left({\left\{{\sigma }_i\right\}}_{1\le i\le r}\right),$$

(20)

where $\mathit{U}$ and $\mathit{V}$ are, respectively, $n_1\times r$ and $n_2\times r$ matrices with orthonormal columns, and the singular values ${\sigma }_i>0$ for $i=1,\cdots ,r$. Given the shrinkage threshold $\tau >0$, the singular value thresholding operator ${\mathcal{D}}_{\tau }$ is defined by

$${\mathcal{D}}_{\tau }\left(\mathit{X}\right):=\mathit{U}{\mathcal{D}}_{\tau }\left(\mathsf{\Sigma }\right){\mathit{V}}^T,{\mathcal{D}}_{\tau }\left(\mathsf{\Sigma }\right)=\mathrm{diag}\left({\left\{\max ({\sigma }_i-\tau ,0)\right\}}_{1\le i\le r}\right).$$

(21)

According to the deduction of [7], for each $\tau >0$ and $\mathit{Y}\in {\mathbb{R}}^{n_1\times n_2}$, the singular value shrinkage operator (21) obeys

$${\mathcal{D}}_{\tau }\left(\mathit{Y}\right)=\underset{\mathit{X}}{\mathrm{argmin}}\left\{{\tau \parallel \mathit{X}\parallel }_{*}+\frac{1}{2}{\parallel \mathit{X}-\mathit{Y}\parallel }_F^2\right\}.$$

(22)

The SVT algorithm utilizes the singular value thresholding operator ${\mathcal{D}}_{\tau }$ and its property (22) to handle the subproblem (19). Consider the proximal problem for $\tau >0$,

$$\begin{array}{c}\hfill \underset{\mathit{B}}{\min }{\tau \parallel \mathit{B}\parallel }_{*}+\frac{1}{2}{\parallel \mathit{B}\parallel }_F^2\\ \hfill s.t.{\mathit{P}}_{\mathsf{\Omega }}\left(\mathit{B}{\mathit{M}}^T\right)={\mathit{P}}_{\mathsf{\Omega }}\left({\mathit{T}}_{\left(n\right)}\right),\end{array}$$

(23)

the solution of which converges to that of (19) as $\tau \to \infty$. This optimization problem can be solved by a Lagrangian multiplier method known as the Uzawa’s algorithm [28].

Given the Lagrangian of (23)

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{B},\mathit{Y})={\tau \parallel \mathit{B}\parallel }_{*}+\frac{1}{2}{\parallel \mathit{B}\parallel }_F^2+\langle \mathit{Y},{\mathit{P}}_{\mathsf{\Omega }}({\mathit{T}}_{\left(n\right)}-\mathit{B}{\mathit{M}}^T)\rangle ,\end{array}$$

(24)

where $\mathit{Y}$ has the same size as ${\mathit{T}}_{\left(n\right)}$, the optimal ${\mathit{B}}^{*}$ and ${\mathit{Y}}^{*}$ should satisfy

$$\mathcal{L}({\mathit{B}}^{*},{\mathit{Y}}^{*})=\underset{\mathit{B}}{\inf }\underset{\mathit{Y}}{\sup }\mathcal{L}(\mathit{B},\mathit{Y})=\underset{\mathit{Y}}{\sup }\underset{\mathit{B}}{\inf }\mathcal{L}(\mathit{B},\mathit{Y}).$$

(25)

Starting with ${\mathit{Y}}^{\left(0\right)}=\mathbf{0}$, Uzawa’s algorithm finds the saddle point $({\mathit{B}}^{*},{\mathit{Y}}^{*})$ through an iterative procedure given by

$$\left\{\begin{array}{c}{\mathit{B}}^{\left(k\right)}={\mathrm{argmin}}_{\mathit{B}}\mathcal{L}(\mathit{B},{\mathit{Y}}^{(k-1)})\hfill \\ {\mathit{Y}}^{\left(k\right)}={\mathit{Y}}^{(k-1)}+{\delta }_k{\mathit{P}}_{\mathsf{\Omega }}({\mathit{T}}_{\left(n\right)}-{\mathit{B}}^{\left(k\right)}{\mathit{M}}^T),\hfill \end{array}\right.$$

(26)

where ${\left\{{\delta }_k\right\}}_{k\ge 1}>0$ are scalar step sizes. The sequence $\left\{{\mathit{B}}^{\left(k\right)}\right\}$ converges to the unique solution to (23). The update of $\mathit{Y}$ is actually based on the gradient descent method if we note that

$${\partial }_{\mathit{Y}}\mathcal{L}(\mathit{B},\mathit{Y})={\mathit{P}}_{\mathsf{\Omega }}({\mathit{T}}_{\left(n\right)}-\mathit{B}{\mathit{M}}^T).$$

(27)

Now, we have to compute the minimizer of (26). Observe that the factor matrices $\{{\mathit{A}}_k:k=1,\cdots ,N\}$ have orthogonal columns. Based on the orthogonal invariance property of the Frobenius norm, we have

$${\parallel \mathcal{B}\parallel }_F^2={\parallel \mathcal{B}{\times }_1{\mathit{A}}_1^{\left(k\right)}\cdots {\times }_{n-1}{\mathit{A}}_{n-1}^{\left(k\right)}{\times }_{n+1}{\mathit{A}}_{n+1}^{(k-1)}\cdots {\times }_N{\mathit{A}}_N^{(k-1)}\parallel }_F^2,$$

(28)

which in the matrix form gives that

$${\parallel \mathit{B}\parallel }_F^2={\parallel \mathit{B}{\mathit{M}}^T\parallel }_F^2.$$

(29)

Utilizing this property,

$$\begin{array}{cc}\hfill & {\mathrm{argmin}}_{\mathit{B}}{\tau \parallel \mathit{B}\parallel }_{*}+\frac{1}{2}{\parallel \mathit{B}\parallel }_F^2+\langle \mathit{Y},{\mathit{P}}_{\mathsf{\Omega }}({\mathit{T}}_{\left(n\right)}-\mathit{B}{\mathit{M}}^T)\rangle \hfill \\ \hfill =& {\mathrm{argmin}}_{\mathit{B}}{\tau \parallel \mathit{B}\parallel }_{*}+\frac{1}{2}{\parallel \mathit{B}{\mathit{M}}^T\parallel }_F^2+\langle \mathit{Y},{\mathit{P}}_{\mathsf{\Omega }}({\mathit{T}}_{\left(n\right)}-\mathit{B}{\mathit{M}}^T)\rangle \hfill \\ \hfill =& {\mathrm{argmin}}_{\mathit{B}}{\tau \parallel \mathit{B}\parallel }_{*}+\frac{1}{2}{\parallel \mathit{B}{\mathit{M}}^T-{\mathit{P}}_{\mathsf{\Omega }}\left(\mathit{Y}\right)\parallel }_F^2\hfill \\ \hfill =& {\mathrm{argmin}}_{\mathit{B}}{\tau \parallel \mathit{B}\parallel }_{*}+\frac{1}{2}{\parallel \mathit{B}-{\mathit{P}}_{\mathsf{\Omega }}\left(\mathit{Y}\right)\mathit{M}\parallel }_F^2.\hfill \end{array}$$

(30)

According to (22), the optimal ${\mathit{B}}^{*}$ is given by ${\mathcal{D}}_{\tau }\left({\mathit{P}}_{\mathsf{\Omega }}\left(\mathit{Y}\right)\mathit{M}\right)={\mathcal{D}}_{\tau }\left(\mathit{Y}\mathit{M}\right)$ since ${\mathit{P}}_{\mathsf{\Omega }}\left(\mathit{Y}\right)=\mathit{Y}$ for all $k\ge 0$. Therefore, Uzawa’s algorithm finally takes the form

$$\left\{\begin{array}{c}{\mathit{B}}^{\left(k\right)}={\mathcal{D}}_{\tau }\left({\mathit{Y}}^{(k-1)}\mathit{M}\right)\hfill \\ {\mathit{Y}}^{\left(k\right)}={\mathit{Y}}^{(k-1)}+{\delta }_k{\mathit{P}}_{\mathsf{\Omega }}({\mathit{T}}_{\left(n\right)}-{\mathit{B}}^{\left(k\right)}{\mathit{M}}^T),\hfill \end{array}\right.$$

(31)

also referred to as the shrinkage iterations in SVT.

To obtain an approximation to the solution of (19), we choose a large enough $\tau$ and perform the iterations (31) until the stopping criteria $\parallel \mathit{T}-{\mathit{B}}^{\left(k\right)}{\mathit{M}}^T{\parallel }_{\mathsf{\Omega }}^2<ϵ{\parallel \mathit{T}\parallel }_{\mathsf{\Omega }}^2$ are reached, starting with ${\mathit{Y}}^{\left(0\right)}=\mathbf{0}\in {\mathbb{R}}^{I_1\times \cdots \times I_N}$.

The overall process of solving the subproblem (18) is shown in Algorithm 3.

Algorithm 3 The SVT algorithm for solving (18)

Input: Set of observed indices $\mathsf{\Omega }$, tensor with observed entries ${\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right)\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$, fixed factor matrices ${\mathit{A}}_m\in {\mathbb{R}}^{I_m\times R_m}$ for $m=\{1,2,\cdots ,N\}/n$, error tolerance $ϵ$, shrinkage threshold $\tau$ and scalar step sizes ${\left\{{\delta }_k\right\}}_{k\ge 1}$. Output: Optimized $\mathcal{B}$. 1: ${\mathit{P}}_{\mathsf{\Omega }}\left(\mathit{T}\right)\leftarrow$ mode-n unfolding matrix of ${\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right)$. 2: $\mathit{M}\leftarrow {\mathit{A}}_N\cdots \otimes {\mathit{A}}_{n+1}\otimes {\mathit{A}}_{n-1}\cdots \otimes {\mathit{A}}_1$. 3: Initialize ${\mathit{B}}^{\left(0\right)}\leftarrow \left({\mathit{P}}_{\mathsf{\Omega }}\left(\mathit{T}\right)\right)\mathit{M}$,   ${\mathit{Y}}^{\left(0\right)}\leftarrow \mathbf{0}$. 4: $k\leftarrow 0$. 5: while$\parallel \mathit{T}-{\mathit{B}}^{\left(k\right)}{\mathit{M}}^T{\parallel }_{\mathsf{\Omega }}^2\ge ϵ{\parallel \mathit{T}\parallel }_{\mathsf{\Omega }}^2$do 6:     $k\leftarrow k+1$. 7:     Update ${\mathit{B}}^{\left(k\right)}\leftarrow {\mathcal{D}}_{\tau }\left({\mathit{Y}}^{(k-1)}\mathit{M}\right)$, where ${\mathcal{D}}_{\tau }$ is defined in (21). 8:     Update ${\mathit{Y}}^{\left(k\right)}\leftarrow {\mathit{Y}}^{(k-1)}+{\delta }_k{\mathit{P}}_{\mathsf{\Omega }}({\mathit{T}}_{\left(n\right)}-{\mathit{B}}^{\left(k\right)}{\mathit{M}}^T)$. 9: end while 10: $\mathcal{B}\leftarrow {\mathit{B}}^{\left(k\right)}$ folding at mode-n.

4.3. The Rank-Adaptive Tensor Completion Algorithm

This section summarizes our proposed method, namely the rank-adaptive tensor completion algorithm. Based on Tucker decomposition, this algorithm can complete a tensor with only a few observed entries and adaptively estimate its multilinear rank.

Given a tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times \cdots \times I_N}$, only entries with indices in the set $\mathsf{\Omega }$ are observed. Our goal is to predict the unobserved entries based on the premise that $\mathcal{T}$ has a low-rank structure. In this work, we consider the multilinear rank of the tensor. We solve this problem by finding a low multilinear rank tensor $\mathcal{X}$ whose entries on the observation part satisfy ${\parallel \mathcal{X}-\mathcal{T}\parallel }_{\mathsf{\Omega }}^2<ϵ{\parallel \mathcal{T}\parallel }_{\mathsf{\Omega }}^2$, where $ϵ$ is a given tolerance. Here, we restate our problem setting (15)

$$\begin{array}{c}\hfill \underset{\mathcal{G},{\mathit{A}}_1\cdots ,{\mathit{A}}_N}{\min }(R_1,R_2,\cdots ,R_N)s.t.{\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{X}\right)={\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right),\\ \hfill \mathcal{X}=\mathcal{G}{\times }_1{\mathit{A}}_1\cdots {\times }_N{\mathit{A}}_N,{\mathit{A}}_n^T{\mathit{A}}_n=I^{R_n\times R_n},n=1,\cdots ,N.\end{array}$$

Similar to HOOI, we solve (15) by alternative iterations. By fixing all the factor matrices except one, we obtain the rank minimization subproblem (17). Since it is NP-hard, we instead consider its convex relaxation form (18), i.e., minimizing the nuclear norm, which can be handled by the SVT algorithm. The solution can be computed using Algorithm 3. After optimized $\mathcal{B}$ in (18) is obtained, we update $R_n^{\left(k\right)}=\min (R_n^{(k-1)},\mathrm{rank}\left({\mathit{B}}_{\left(n\right)}\right))$ to ensure that $R_n^{\left(k\right)}$ is monotone decreasing. We perform the iterations until some stopping criteria are satisfied, e.g., the estimated rank $(R_1,R_2,\cdots ,R_N)$ remains unchanged, and the factor matrices have slight improvement, or the maximum number of iterations is reached. We present the detailed procedure of our proposed method in Algorithm 4.

Algorithm 4 The rank-adaptive tensor completion based on Tucker decomposition (RATC-TD) algorithm

Input: Set of observed indices $\mathsf{\Omega }$, tensor with observed entries ${\mathcal{P}}_{\mathsf{\Omega }}\left(\mathcal{T}\right)\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ and initial guess of factor matrices ${\mathit{A}}_n^{\left(0\right)}\in {\mathbb{R}}^{I_n\times R_n^{\left(0\right)}}$ for $n=1,2,\cdots ,N$. Output: Estimated rank $(R_1^{\left(k\right)},R_2^{\left(k\right)},\cdots ,R_N^{\left(k\right)})$, core tensor $\mathcal{G}\in {\mathbb{R}}^{R_1^{\left(k\right)}\times \cdots \times R_N^{\left(k\right)}}$ and factor matrices ${\mathit{A}}_n^{\left(k\right)}\in {\mathbb{R}}^{I_n\times R_n^{\left(k\right)}}$ for $n=1,2,\cdots ,N$. 1: $k\leftarrow 0$. 2: Initialize $(R_1^{\left(0\right)},R_2^{\left(0\right)},\cdots ,R_N^{\left(0\right)})\leftarrow (I_1,I_2,\cdots ,I_N)$. 3: while not converge do 4:     for all $n=1,2,\cdots ,N$ do 5:         $k\leftarrow k+1$. 6:         $\mathcal{B}\leftarrow$ solution to (18) using Algorithm 3. 7:         ${\mathit{B}}_{\left(n\right)}\leftarrow$ mode-n unfolding matrix of $\mathcal{B}$. 8:         $\mathit{U},\mathsf{\Sigma },{\mathit{V}}^T\leftarrow$ full-SVD of ${\mathit{B}}_{\left(n\right)}$. 9:         $R_n^{\left(k\right)}\leftarrow \min (R_n^{(k-1)},\mathrm{rank}\left({\mathit{B}}_{\left(n\right)}\right))$. 10:         ${\mathit{A}}_n^{\left(k\right)}\leftarrow {\mathit{U}}_{:,1:R_n^{\left(k\right)}}$. 11:     end for 12:     $\mathcal{G}\leftarrow {\mathsf{\Sigma }}_{1:R_n^{\left(k\right)},1:R_n^{\left(k\right)}}{\mathit{V}}_{:,1:R_n^{\left(k\right)}}^T$ folding at mode-n. 13: end while

5. Numerical Experiments

In this section, numerical experiments are conducted to the effectiveness of our proposed rank-adaptive tensor completion based on Tucker decomposition (RATC-TD) algorithm.

5.1. Test Problem 1: The Recovery of Third-Order Tensors

Many tensor recovery algorithms in the Tucker scheme require a priori given the rank of the tensor, but our method does not require this. In this test problem, synthetic data are considered, and we set the random composition tensor in the same way as [9]. We consider the tensor $\mathcal{T}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$ and obtain tensors by multiplying a randomly generated kernel $\mathcal{G}$ of size $r_1\times r_2\times r_3$ with randomly generated factor matrices ${\mathit{A}}_i$, where ${\mathit{A}}_i\in {\mathbb{R}}^{I_i\times r_i},i=1,2,3$. So, the tensor $\mathcal{T}$ can be represented by $\mathcal{G}{\times }_1{\mathit{A}}_1{\times }_2{\mathit{A}}_2{\times }_3{\mathit{A}}_3$, and the number of entries of $\mathcal{T}$ is denoted by $\left|\mathcal{T}\right|$. By setting tensors in this way, we can obtain random tensors and ensure that the obtained tensors have a low-rank structure, which can be handled well by the classical low-rank tensor recovery methods.

In this experiment, we set the tensor size as $50\times 50\times 50$, the kernel size as $5\times 5\times 5$, and the kernel data are randomly generated from the interval $[0,1]$. The factor matrices are of size $50\times 5$, and the data of factor matrices are randomly generated from $[-0.5,0.5]$. In order to show the robustness and reliability of our proposed algorithm, we add a small perturbation to the synthetic data. We add Gaussian noise with zero mean and standard deviation as $0.1$ times the element mean and take the tensor with Gaussian noise as the ground truth data. As described above, our algorithm does not require the rank in advance, and the algorithm takes the size of the tensor as the initial tensor rank. With the continuous operation of the algorithm, the tucker rank is constantly reduced, and the exact rank can be approximated.

We compare the proposed algorithm with other Tucker-decomposition-based tensor completion algorithms (that require appropriate ranks given a priori), including gHOI [24], M${}^2$SA [11] and Tucker–Wopt [5], and test the sample estimation error $erro{r}_{obs}$ and out-of-sample estimation error $erro{r}_{val}$ with different initial ranks at sampling rates r (the proportion of known data in the total data) of $0.05$, $0.1$ and $0.2$, respectively. The size of the observed index set $\mathsf{\Omega }$ is $r\left|\mathcal{T}\right|$, and each index in $\mathsf{\Omega }$ is randomly generated through the uniform distribution. $\overline{\mathsf{\Omega }}$ denotes the complementary set of $\mathsf{\Omega }$. The specific definitions of two kinds of errors are given by the following formulas:

$$\begin{array}{c}\hfill erro{r}_{obs}=\frac{\parallel {\mathit{P}}_{\mathsf{\Omega }}{(\mathcal{X}-\mathcal{T})\parallel }_F^2}{\parallel {\mathit{P}}_{\mathsf{\Omega }}{\left(\mathcal{T}\right)\parallel }_F^2},\end{array}$$

(32)

$$\begin{array}{c}\hfill erro{r}_{val}=\frac{\parallel {\mathit{P}}_{\overline{\mathsf{\Omega }}}{(\mathcal{X}-\mathcal{T})\parallel }_F^2}{\parallel {\mathit{P}}_{\overline{\mathsf{\Omega }}}{\left(\mathcal{T}\right)\parallel }_F^2},\end{array}$$

(33)

where $\mathcal{X}$ is the result obtained by completion algorithms.

To ensure the generalization of the proposed algorithms, we set the error tolerance $ϵ$ on the observed data to 0.0025, which means we stop optimizing $\mathcal{B}$ with SVT when the relative error is less than $ϵ$ or the maximum number of iterations is reached. When using SVT to optimize $\mathcal{B}$, we set the relative error $erro{r}_{SVT}$ in the following format:

$$\begin{array}{c}\hfill erro{r}_{SVT}=\frac{\parallel {\mathit{P}}_{\mathsf{\Omega }}(\mathit{T}-{\mathit{B}}^{\left(k\right)}{\mathit{M}}^T){\parallel }_F^2}{\parallel {\mathit{P}}_{\mathsf{\Omega }}{\left(\mathit{T}\right)\parallel }_F^2}.\end{array}$$

(34)

For the estimated tensor $\mathcal{X}$ obtained at each iteration in Algorithm 4, the relative error is assessed as (32) and (33), for other methods considered, the relative errors are also assessed using the same form. In addition to the relative error and the maximum number of iterations, Algorithm 4 also terminates if the difference of the relative errors obtained from the algorithm for two consecutive steps is less than a certain threshold $\eta$, i.e.,

$$\begin{array}{c}\hfill \left|\frac{erro{r}_{iter}^{(k+1)}-erro{r}_{iter}^{\left(k\right)}}{erro{r}_{iter}^{\left(k\right)}}\right|<\eta ,\end{array}$$

(35)

where the $erro{r}_{iter}^{\left(k\right)}$ represents the relative error obtained at the kth iteration.

The initial rank is set to $R^{\left(0\right)}=[r,r,r],r=5,10,15$.

Table 1. Relative error (32) comparison with different initial rank on observed data, test problem 1.

Sampling Ratio	Initial Rank	M${}^2$SA	gHOI	Tucker–Wopt	RATC-TD
0.05	5, 5, 5	0.0354	0.0431	0.0354	0.0354
0.05	10, 10, 10	0.0322	0.0404	0.0284	0.0354
0.05	15, 15, 15	0.0270	0.0362	0.1128	0.0354
0.1	5, 5, 5	0.0264	0.0278	0.0264	0.0264
0.1	10, 10, 10	0.0250	0.0267	0.0233	0.0264
0.1	15, 15, 15	0.0229	0.0229	0.0199	0.0264
0.2	5, 5, 5	0.0235	0.0236	0.0235	0.0235
0.2	10, 10, 10	0.0230	0.0229	0.0225	0.0235
0.2	15, 15, 15	0.0223	0.0215	0.0208	0.0235

and

Table 2. Relative error (33) comparison with different initial rank on unknown data, test problem 1.

Sampling Ratio	Initial Rank	M${}^2$SA	gHOI	Tucker–Wopt	RATC-TD
0.05	5, 5, 5	0.0187	0.0518	0.0187	0.0186
0.05	10, 10, 10	0.3438	0.4192	0.0429	0.0186
0.05	15, 15, 15	0.6262	0.6143	1.4444	0.0186
0.1	5, 5, 5	0.0153	0.0231	0.0153	0.0153
0.1	10, 10, 10	0.2300	0.3102	0.0241	0.0153
0.1	15, 15, 15	0.4935	0.4521	0.0610	0.0153
0.2	5, 5, 5	0.0145	0.0162	0.0145	0.0145
0.2	10, 10, 10	0.2480	0.2390	0.0181	0.0145
0.2	15, 15, 15	0.4872	0.2800	0.0222	0.0145

respectively show recovery error (32) and (33) of test problem 1. It can be seen that when the given initial rank does not match the ground truth data, gHOI, M${}^2$SA and Tucker–Wopt have large errors, while our RATC-TD has small errors. We emphasize that our proposed method does not require a given tensor rank; the initial rank of our proposed algorithm is set to the size of the tensor we aim to complete, i.e., $R^{\left(0\right)}=$ size$\left(\mathcal{T}\right)$.

In this experiment, we use our RATC-TD algorithm to estimate the tensor ranks. We next use the M${}^2$SA method to improve the estimated results. Our algorithm can be considered an initial step for other algorithms, and the effects of tensor completion are also sound when only our method is used for optimization. The recovery results without using M${}^2$SA are shown in

Table 3. Numerical results only using our proposed algorithm RATC-TD, test problem 1.

Sampling Ratio	${\mathit{error}}_{\mathit{obs}}$	${\mathit{error}}_{\mathit{val}}$
0.05	0.0500	0.0410
0.1	0.0489	0.0523
0.2	0.0495	0.0443

. As seen from Table 2, since our method was used for optimization in advance, compared with the results obtained by simply using the M${}^2$SA method with a given rank, it gives a better recovery effect. Notably, our proposed algorithm can also stably complete the tensor under a small sampling rate and does not need to give an appropriate rank in advance.

5.2. Test Problem 2: The Recovery of Fourth-Order Tensors

Similarly to test problem 1, following the procedures in [9], we in this test problem generate the object tensor $\mathcal{T}\in {\mathbb{R}}^{30\times 30\times 30\times 30}$ by multiplying a randomly generated kernel $\mathcal{G}\in {\mathbb{R}}^{5\times 5\times 5\times 5}$ with randomly generated factor matrices of size ${\mathbb{R}}^{30\times 5}$. The ground truth data is defined by $\mathcal{T}$ with Gaussian noise added (the same setting as in test problem 1). For gHOI, M${}^2$SA, Tucker–Wopt and our proposed RATC-TD, different initial ranks are tested, and the sampling rate r is set to $0.1$. The corresponding number of entries in the observed data set $\mathsf{\Omega }$ is $r\left|\mathcal{T}\right|$, and the index of each entry in $\mathsf{\Omega }$ is randomly generated through the uniform distribution. The errors for this test problem are shown in

Table 4. The results of tensor completion quality of methods in the case of sampling rate r = 0.1, test problem 2.

Types of Error	Initial Rank	M${}^2$SA	gHOI	Tucker–Wopt	RATC-TD
$erro{r}_{obs}$	5, 5, 5, 5	0.0300	0.0342	0.0300	0.0300
$erro{r}_{val}$	5, 5, 5, 5	0.0147	0.0251	0.0147	0.0147
$erro{r}_{obs}$	15, 15, 15, 15	0.0258	0.0255	0.1552	0.0300
$erro{r}_{val}$	15, 15, 15, 15	0.7130	0.4807	6.6025	0.0147

, where $erro{r}_{obs}$ is the error on the observed data set (32) and $erro{r}_{val}$ is the error on the validation set (33). It can be seen that our RATC-TD method performs well for this test problem.

5.3. Test Problem 3: The Recovery of Real Missing Pictures

In this test problem, real practical data are considered, and the performance of our RATC-TD is compared with that of gHOI and M${}^2$SA. The initial complete images considered are shown in Figure 1. The data format of each image can be regarded as a third-order tensor. Each image here is stored as a tensor with size $256\times 256\times 3$, where the third dimension is 3, representing three color channels of red, green, and blue. The ground truth data $\mathcal{T}$ for this test problem are the images in Figure 1 (details are as follows).

Two ways to construct partial images are considered, and our method is tested to recover the original images using the partial images. First, some black lines are added to the images, which can be considered a kind of structural missing, and the corrupted images are shown in the first column of Figure 2. The black line parts of images correspond to $\overline{\mathsf{\Omega }}$, and the rest correspond to $\mathsf{\Omega }$. Second, after the initial image is converted into the data format of a third-order tensor $\mathcal{T}$, the observed index set $\mathsf{\Omega }$ is constructed with each index randomly generated through the uniform distribution, and the size of $\mathsf{\Omega }$ is $r\left|\mathcal{T}\right|$ (where the sampling rate is $r=0.1$). The images associated with $\mathsf{\Omega }$ are shown in the first column of Figure 3.

From Figure 2, it is clear that the images obtained by our RATC-TD are closer to the ground truth images than the ones obtained by gHOI and M${}^2$SA. The numerical results representing the qualities of recovery are shown in

Table 5. Numerical characterization of the recovery quality of images with black lines, test problem 3.

Girl	PSNR	SSIM	${\mathit{error}}_{\mathit{obs}}$	${\mathit{error}}_{\mathit{val}}$
missing figure	11.4347	0.5044	/	/
gHOI	24.8558	0.9536	0.0236	0.0985
M${}^2$SA	24.8196	0.9532	0.0235	0.1080
RATC-TD	28.9981	0.9833	0.0099	0.0751
Baboon	PSNR	SSIM	${\mathit{error}}_{\mathit{obs}}$	${\mathit{error}}_{\mathit{val}}$
missing figure	8.6372	0.1878	/	/
gHOI	18.3661	0.7477	0.0268	0.1670
M${}^2$SA	17.2923	0.6948	0.0263	0.2341
RATC-TD	20.9857	0.8612	0.0098	0.1004
Boat	PSNR	SSIM	${\mathit{error}}_{\mathit{obs}}$	${\mathit{error}}_{\mathit{val}}$
missing figure	8.5502	0.3229	/	/
gHOI	17.4365	0.8717	0.0282	0.2419
M${}^2$SA	16.9671	0.8556	0.0276	0.2635
RATC-TD	22.0162	0.9524	0.0099	0.0998

. In addition to the errors $erro{r}_{obs}$ (32) on the observation set and the errors $erro{r}_{val}$ (33) on the unobserved data set, we also use SSIM (structure similarity index measure) [29] and PSNR (peak signal-to-noise ratio) parameters to evaluate the effect of our image restoration, which are often used to assess the quality of image restoration. SSIM ranges from 0 to 1; the larger the SSIM is, the smaller the image distortion is. Similarly, the larger the PSNR is, the less distortion there is. PSNR is used to measure the difference between two images. For the restored image $\mathcal{X}$ and the original image $\mathcal{T}$, the PSNR between them is given by the following formula:

$$\begin{array}{c}\hfill \mathrm{PSNR}=10\times lg\left(\frac{{\mathrm{Maxpixel}}^2}{\mathrm{MSE}}\right),\end{array}$$

(36)

where $\mathrm{Maxpixel}$ is the maximum value of the image pixel. In our test problem, $\mathrm{Maxpixel}=255$, and $\mathrm{MSE}$ is the mean square error between $\mathcal{X}$ and $\mathcal{T}$. For the case with sampling rate $r=0.1$, the corresponding restoration results are shown in Figure 3, and the numerical results representing the qualities of recovery are shown in

Table 6. Numerical characterization of the recovery quality of images obtained by random sampling at sampling rate of 0.1, test problem 3.

Girl	PSNR	SSIM	${\mathit{error}}_{\mathit{obs}}$	${\mathit{error}}_{\mathit{val}}$
missing figure	5.5747	0.0329	/	/
gHOI	13.2045	0.6585	0.0309	0.1336
M${}^2$SA	11.9943	0.5977	0.0305	0.1663
RATC-TD	18.9566	0.8049	0.0100	0.0764
Baboon	PSNR	SSIM	${\mathit{error}}_{\mathit{obs}}$	${\mathit{error}}_{\mathit{val}}$
missing figure	5.8165	0.0306	/	/
gHOI	13.6954	0.5862	0.0319	0.1121
M${}^2$SA	12.8840	0.5404	0.0351	0.1441
RATC-TD	17.4101	0.7008	0.0098	0.0968
Boat	PSNR	SSIM	${\mathit{error}}_{\mathit{obs}}$	${\mathit{error}}_{\mathit{val}}$
missing figure	5.6110	0.0346	/	/
gHOI	12.5033	0.7168	0.0312	0.1336
M${}^2$SA	11.6136	0.6549	0.0343	0.1615
RATC-TD	17.4527	0.8735	0.0096	0.0869

. In this case, we only use 10 percent of the ground truth data, but it can be seen that our RATC-TD gives effective recovery results for this test problem.

For the gHOI method and the M${}^2$SA method, which need to be given the rank of the initial tensor in advance, we choose the initial rank as $R^{\left(0\right)}=[20,20,3]$. It is worth noting that for real image data, it is difficult for us to know the exact Tucker rank in advance. When there is no way to know the rank in advance, in order to use gHOI and M${}^2$SA to achieve data completion, the initial rank can only be constantly adjusted through experiments, or the rank is given based on prior experience.

6. Conclusions

In this paper, we propose a novel rank-adaptive tensor completion based on the Tucker decomposition (RATC-TD) method. As the existing tensor completion methods based on Tucker decomposition, such as gHOI, ${\mathrm{M}}^2\mathrm{SA}$ and Tucker–Wopt, typically require the tensor rank given a priori, overestimating or underestimating the tensor ranks can lead to poor results. Inspired by the RaHOOI algorithm, we propose our algorithm based on the HOOI structure. Our proposed algorithm can adaptively estimate the multilinear rank of data in the process of tensor completion. We show the algorithm’s effectiveness through experiments on completing synthetic data and genuine pictures. The results of our algorithm can also provide effective initial data for other tensor completion methods, such as ${\mathrm{M}}^2\mathrm{SA}$. The further work we expect is to extend the algorithm to higher-dimensional problems, which requires us to optimize the algorithm further to reduce the algorithm’s computational time.