Iterative p-shrinkage thresholding algorithm for low Tucker rank tensor recovery

Abstract

Low-rank tensor recovery, as a higher order extension of low-rank matrix recovery, has generated a great deal of research interests in recent years, such as image inpainting, video inpainting, video decoding, scan completion, and so on. In this paper, we propose an easy-to-implement algorithm based on the framework of alternative direction method, named iterative p-shrinkage thresholding algorithm, for solving the low Tucker rank tensor recovery problem. The performance of the proposed algorithm is investigated on both synthetic and real data. Numerical results on simulation data demonstrate that our algorithm can successfully recover varieties of synthetic low Tucker rank tensors in different sampling ratios with better quality compared to the existing state-of-the-art tensor recovery algorithms. Experiments on real data, including colored image inpainting, MRI image inpainting and hyperspectral image inpainting, further illustrate the effectiveness of the proposed iterative p-shrinkage thresholding algorithm.

Introduction

Nowadays, there are massive multidimensional data to be faced in practical applications, such as color images, hyperspectral images and videos etc. The more complex data we confront, the larger scale problem we need to deal with. Tensors (i.e., multi-way arrays), as a higher order generalization of vectors and matrices, can well characterize the natural structure of multidimensional data in mathematical. In recent years, tensors and their applications have received extensive research and increasing attentions [23]. As a focus point of research, tensor missing values estimation problem aims to recover the original data form its incomplete observations with prior information. It is one among the key tasks in computer vision and graphics applications, such as image inpainting, video decoding, video inpainting, scan completion, appearance acquisition completion, and so on. For image inpainting, when the required image has low-rank structure, the corresponding tensor missing value estimation problem can be viewed as a low-rank tensor recovery problem, i.e., recovering an unknown low-rank tensor from its degraded observations.

Inspired by the development of low-rank matrix recovery research, low-rank tensor recovery problem has obtained increased research on predicting the missing entries from the incomplete or noisy observed data [1], [9], [15], [16], [24], [27], [35], [36], [43]. However, different from the fact that the matrix rank is unique, the definition of tensor rank is not unique [23], which brings some difficulties to the low-rank tensor recovery research. Among the known tensor ranks, Tucker rank and canonical polyadic (CP) rank are two popular ones, which correspond to the Tucker decomposition [37] and CP decomposition [2], [17] of tensor, respectively. In spite of this fact, many scholars have devoted themselves to dealing with low-rank tensor recovery problem because tensors can well exhibit the nature of real world data generally. At present, the related work is mainly based on two aspects, including tensor completion and tensor decomposition approaches.

For the research on tensor completion, one considers tensor completion problem to recover the original tensor from its incomplete observations. In 2009, Liu et al. [26] solved it by applying its convex relaxation model based on Tucker rank, and reported impressive results on visual data recovery, where the popular nuclear norm was given to approximate the ranks of all the unfolding matrices of tensor. Recently, plenty of followed-up studies have been given based on this framework, and the algorithms for solving the proposed nuclear norm relaxation model include: fast composite splitting algorithms [20], Douglas–Rachford splitting method and alternating direction method of multipliers (ADMM) for tensor recovery with exact update rule [16], block coordinate descent (BCD) method [27], fixed point continuation method [46], splitting augmented Lagrange method [38], [45], tensor hard completion method [35], [49], adaptive correction approach [4], etc. Xu et al. [44] proposed an alternating minimization algorithm to recover the low-rank tensor by applying low-rank matrix factorization to the all-mode matricizations of the underlying tensor, as an alternative of the nuclear norm. Furthermore, Ji et al. [21] generalized the work in [44] by adding a total variation regularization, and developed a BCD method to solve it. Jiang et al. [22] designed a block successive upper-bound minimization algorithm to solve tensor completion combined with the framelet regularization based on low-rank matrix factorization. Beyond convex relaxation model, Cao et al. [10] established a concave relaxation model via folded-concave penalty for low Tucker rank tensor recovery problem, and developed a local linear approximation-augmented Lagrange multiplier algorithm to search the local optimum of the proposed model. Li et al. [29] established a non-convex L_p norm relaxation model, and proposed ADMM based on exact and inexact iteratively reweighted algorithms to solve it.

Except for the above research based on Tucker rank, recently there are also some other tensor completion research based on other tensor ranks. Bengua et al. [3] proposed a new approach to tensor completion based on tensor train rank, as well as two optimization formulations, including nuclear norm minimization and parallel matrix factorization for color image and video recovery. Zhang et al. [48] considered the tensor completion with tensor tubal rank, and established the exact recovery condition of it by its convex relaxation minimization. Madathil et al. [32] applied reweighted nuclear norm minimization method to solve tensor completion with twist tensor total variation regularization based on tensor tubal rank. Moreover, Tang et al. [39] proposed a novel tri-clustering tensor completion framework to improve the performance of social image tag refinement, in which they divided the original image-tag-user associated tensor into sub-tensors by tri-clustering method, and well addressed the challenges of large-scale and super-sparse tensor completion. Li et al. [31] proposed a deep collaborative embedding framework for multiple image understanding tasks.

For the research on tensor factorization, one recovers the missing data from the limited observations based on a multilinear non-convex model under the assumption that the rank is given. Acar et al. [1] used CP decomposition (CPD) technique to deal with the missing data, and formulated it as a weighted least squares problem, named CP weighted optimization. Some other related methods were also investigated, such as structured CPD method by using nonlinear least squares [36] and geometric nonlinear conjugate gradient [24]. Zhao et al. [50] considered the low-rank CPD model for tensor completion, and proposed an automatic relevance determination method to solve tuning parameter problem in a Bayesian framework. Geng et al. [15] proposed an expectation-maximization algorithm which used a higher order singular value decomposition (SVD) [12] to characterize the multilinear subspaces as components of low Tucker rank tensor decomposition model. Chen et al. [8] applied tensor decomposition technique for tensor completion by minimizing the nuclear norms of obtained individual factor tensors. Hong et al. [18] considered multi-sided recommendation based on social tensor factorization, and presented an adaptive method to solve it.

To sum up, methods for tensor completion can well exploit the automatic rank optimization without further assumptions, while in contrast, methods using tensor factorization technique are prone to overfitting due to an incorrect tensor rank or prior estimations of latent factors, which may result in severe deterioration of predictive performance. In the following, we pay our attention on tensor completion related research. To describe the low-rank tensor recovery problem more clearly, some notations and models are introduced below.

In this subsection, we give some basic notations and the low-rank tensor recovery models. Throughout the paper, we let $\left[n\right]=\{1,2,\cdots n\}$ and use ${\mathbb{R}}^n,$ ${\mathbb{R}}^{m\times n}$ to denote the space of n dimensional real column vectors and m × n dimensional real matrices, respectively. For $x\in {\mathbb{R}}^n$ and $X\in {\mathbb{R}}^{m\times n},$ we denote $x={(x_1,\cdots x_n)}^T$ and $X=\left(x_{ij}\right),$ i ∈ [m], j ∈ [n]. Let $\text{Diag}\{u_i:i\in [n\left]\right\}$ denote the diagonal matrix whose ith diagonal element is u_i. For a matrix $X\in {\mathbb{R}}^{m\times n},$ we write its SVD as $X=U\Sigma V^T$ with $U\in {\mathbb{R}}^{m\times n},$ $V\in {\mathbb{R}}^{n\times n}$ and $\Sigma =\text{Diag}\left(\sigma \right(X\left)\right)$ with $\sigma \left(X\right)={({\sigma }_1\left(X\right),\cdots ,{\sigma }_n\left(X\right))}^T\in {\mathbb{R}}^n$ being the singular value vector of X. We assume σ₁(X) ≥ ⋅⋅⋅ ≥ σ_n(X) and n ≤ m without lose of generality.

Let $\mathcal{T}:={\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ be the set of tensors with order N and dimension I₁ × I₂ × ⋅⋅⋅I_N. For $\mathcal{X}\in \mathcal{T},$ we denote $\mathcal{X}=\left(x_{i_1{i}_2\cdots i_N}\right),$ where i_n ∈ [I_n], n ∈ [N]. When N ≤ 2, a zeroth order tensor is a scalar, a first order tensor is a vector, and a second order tensor is a matrix. We call a tensor with third order and higher order as a higher-order tensor, which is called tensor for simplicity in this paper. For any n ∈ [N], if fixing all the indexes of $\mathcal{X}$ except for the nth one, we obtain a fiber or mode-n vector of $\mathcal{X},$ and denote it as $x_{i_1\cdots i_{n-1}:i_{n+1}\cdots i_N}$. Obviously, $x_{i_1\cdots i_{n-1}:i_{n+1}\cdots i_N}\in {\mathbb{R}}^{I_n}$ and $\mathcal{X}$ has ∏_j ≠ nI_j mode-n vectors. Moreover, we call the matrix whose columns are mode-n vectors of tensor as mode-n unfolding of $\mathcal{X},$ and denote it as $X_{\left(n\right)}\in {\mathbb{R}}^{I_n\times {\prod }_{k\ne n}{I}_k}$.

For $\mathcal{X}\in \mathcal{T},$ we define its n-rank as the rank of mode-n unfolding X_(n) for any n ∈ [N]. An Nth order tensor has N n-ranks, since it has N mode-n unfoldings. Different from the fact that row rank and column rank of a matrix are equivalent, n-ranks of a tensor are generally different. For $\mathcal{X}\in \mathcal{T},$ we define its Tucker rank as

$${\text{rank}}_T\left(\mathcal{X}\right):=(\text{rank}\left(X_{\left(1\right)}\right),\text{rank}\left(X_{\left(2\right)}\right),\cdots ,\text{rank}\left(X_{\left(N\right)}\right)),$$

which relates to its Tucker decomposition

$$\begin{array}{c}\hfill \mathcal{X}=\mathcal{C}{\times }_1{U}_1{\times }_2{U}_2\cdots {\times }_N{U}_N,\end{array}$$

where $\mathcal{C}\in {\mathbb{R}}^{r_1\times \cdots \times r_N}$ is a core tensor, $U_n\in {\mathbb{R}}^{I_n\times r_n}$ is the left singular matrix of X_(n) for any n ∈ [N]. Here, n-mode product of tensor $\mathcal{S}\in \mathcal{T}$ and matrix $Q\in {\mathbb{R}}^{J\times I_n}$ is defined as $\mathcal{S}{\times }_{n}Q\in {\mathbb{R}}^{I_1\times \cdots \times I_{n-1}\times J\times I_{n+1}\times \cdots \times I_N}$ with

$${\left(\mathcal{S}{\times }_{n}Q\right)}_{i_1\cdots i_{n-1}j{i}_{n+1}\cdots i_N}:=\sum _{i_n=1}^{I_n}{x}_{i_1{i}_2\cdots i_N}{u}_{j{i}_n},j\in \left[J\right],i_k\in \left[I_k\right],k\in \left[N\right].$$

It is known that $\text{rank}\left(X_{\left(n\right)}\right)=r_n$ for n ∈ [N].

Recall that tensor completion problem is given as

$$\begin{array}{c}\hfill \underset{\mathcal{X}\in \mathcal{T}}{min}\mathrm{rank}\left(\mathcal{X}\right)\text{s.t.}{\mathcal{X}}_{\Omega }={\mathcal{M}}_{\Omega },\end{array}$$

where $\mathcal{M}\in \mathcal{T}$ is known, and ${\mathcal{X}}_{\Omega }\in \mathcal{T}$ has the same elements with tensor $\mathcal{M}$ in index set Ω. Also, its nuclear clear norm relaxation model is given as

$$\begin{array}{c}\hfill \underset{\mathcal{X}\in \mathcal{T}}{min}\sum _{i=1}^N{\parallel X_{\left(i\right)}\parallel }_{*}\text{s.t.}{\mathcal{X}}_{\Omega }={\mathcal{M}}_{\Omega },\end{array}$$

where $\parallel X_{\left(i\right)}{\parallel }_{*}:={\sum }_{j=1}^{I_i}{\sigma }_j\left(X_{\left(i\right)}\right)$ is the nuclear norm of X_(i) for each i ∈ [N]. When completion constraint ${\mathcal{X}}_{\Omega }={\mathcal{M}}_{\Omega }$ is extended to general linear constraint $\mathcal{A}\left(\mathcal{X}\right)=b,$ where $\mathcal{A}:\mathcal{T}\to {\mathbb{R}}^s$ is a linear map and $b\in {\mathbb{R}}^s$with $s\le {\prod }_{i=1}^N{I}_i,$ tensor completion turns to a special case of the low-rank tensor recovery. In the following, we consider low Tucker rank tensor recovery as

$$\begin{array}{c}\hfill \underset{\mathcal{X}\in \mathcal{T}}{min}{\mathrm{rank}}_T\left(\mathcal{X}\right)\text{s.t.}\mathcal{A}\left(\mathcal{X}\right)=b.\end{array}$$

Since ${\mathrm{rank}}_T\left(\mathcal{X}\right)$ is an N-way array, problem (2) is a multi-objective programming. Naturally, as discussed in [27], [35], we consider the following low Tucker rank tensor recovery problem:

$$\begin{array}{c}\hfill \underset{\mathcal{X}\in \mathcal{T}}{min}\sum _{i=1}^{\widehat{m}}{\alpha }_i\text{rank}\left(X_{\left(i\right)}\right)\text{s.t.}\mathcal{A}\left(\mathcal{X}\right)=b,\end{array}$$

where α_i > 0, $\widehat{m}\in \left[N\right]$.

Low-rank tensor recovery problem can be seen as a higher order generalization of sparse signal recovery (SSR) problem and low-rank matrix recovery (LMR) problem. In fact, if N takes 1 and 2, tensor reduces to vector and matrix, and problem (2) respectively turns to SSR problem

$$\begin{array}{c}\hfill \underset{x\in {\mathbb{R}}^n}{min}{\parallel x\parallel }_0\text{s.t.}Ax=b,\end{array}$$

and LMR problem

$$\begin{array}{c}\hfill \underset{X\in {\mathbb{R}}^{m\times n}}{\min }\text{rank}\left(X\right)\text{s.t.}{\mathcal{A}}_0\left(X\right)=b_0,\end{array}$$

where ‖x‖₀ means the number of the nonzero elements of x, $A\in {\mathbb{R}}^{s\times n},$ ${\mathcal{A}}_0:{\mathbb{R}}^{m\times n}\to {\mathbb{R}}^{s_0},$ $b\in {\mathbb{R}}^s,$ $b_0\in {\mathbb{R}}^{s_0}$ with s ≤ n and s₀ ≤ mn. Over the last ten years, lots of effective algorithms for solving SSR and LMR problems are proposed, include convex and non-convex relaxation methods, iterative algorithms and greedy type algorithms, etc. In the literature of SSR and LMR, non-convex Schatten-p quasi-norm minimization with 0 < p < 1 has been shown to give better recovery results than convex l₁ or nuclear norm minimization both theoretically and numerically. Especially, in the matrix recovery case, Malek-Mohammadi et al. [34] theoretically proved that the ranks of matrices recovered by Schatten-p quasi-norm minimization are equal or larger than what is guaranteed by nuclear norm minimization, while Liu et al. [25] showed that Schatten-p quasi-norm minimization with $p=2/3$ can guarantee the required matrix recovery in a much weaker condition than nuclear norm minimization. Moreover, Li et al. [30] proposed a robust structured nonnegative matrix factorization learning framework for image representation by leveraging the block-diagonal structure and the l_{2, p}-norm with 0 < p ≤ 1.

Recently, Voronin et al. [40] proposed an effective p-shrinkage thresholding algorithm for solving SSR problem, which is a generalization of the well-known iterative soft thresholding algorithm [11]. The global convergence of the proposed algorithm was proved, and the exact recovery condition for recovering sparse signal was established in [42]. Both theoretical and empirical results of p-shrinkage thresholding algorithm for solving SSR problem greatly support the idea that minimizing penalty function induced by p-shrinkage thresholding operator is an good alternative to classical convex l₁ norm and non-convex Schatten-p quasi-norm minimization relaxation.

Note that data in practical often have some certain spatial structure, and tensors can well represent the data with special structure. Also, the related spatial data recovery problem can be characterized as a low rank tensor recovery problem. A natural question or idea is that how to effectively solve it? The great success of the above p-shrinkage thresholding algorithm for solving SSR problem motives us to consider and apply the similar p-shrinkage strategy and design the related method to solve it. So, we consider to design an iterative p-shrinkage thresholding algorithm for low Tucker rank tensor recovery problem. Meanwhile, due to the developments of ADMM for nonconvex optimization problems [14], [19], [28], [47], we propose an iterative p-shrinkage thresholding algorithm based on ADMM for solving the non-convex relaxation model of problem (3) in this paper.

Moreover, we implement the proposed algorithm both on synthetic data and real data in the experiments, and state the advantages of the proposed algorithm in the following three points. First of all, the proposed algorithm establishes a unified framework of the existing iterative soft and hard thresholding shrinkage algorithms, which shrinks less than the soft thresholding shrinkage for the large elements in the iterations, and overcomes the discontinuity of the hard shrinkage operator. Second, numerical results demonstrate the effectiveness of the proposed algorithm compared to the existing state-of-the-art algorithms, it shows that the proposed algorithm performs well especially in the case of $p=-0.5,$ where the negative value is not involved in the other tensor recovery research. Last, the proposed algorithm is designed based on the popular non-convex ADMM method, which is easy to be implemented.

The rest of this paper is organized as follows. In Section 2, we give some preliminary concepts and results. In Section 3, we propose an iteratively p-shrinkage thresholding algorithm for solving low Tucker rank tensor recovery problem. Numerical experiments for the proposed algorithm based on simulation data and real data recovery are reported in Section 4. In Section 5, we give some conclusions.