Orthogonal tubal rank-1 tensor pursuit for tensor completion

Accepted Manuscript

Orthogonal Tubal Rank-1 Tensor Pursuit for Tensor Completion Weize Sun, Lei Huang, H.C. So, Jiajia Wang PII: DOI: Reference:

S0165-1684(18)30381-5 https://doi.org/10.1016/j.sigpro.2018.11.015 SIGPRO 6989

To appear in:

Signal Processing

Received date: Revised date: Accepted date:

14 November 2017 18 October 2018 20 November 2018

Please cite this article as: Weize Sun, Lei Huang, H.C. So, Jiajia Wang, Orthogonal Tubal Rank-1 Tensor Pursuit for Tensor Completion, Signal Processing (2018), doi: https://doi.org/10.1016/j.sigpro.2018.11.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A novel tensor completion method based on tensor tubal rank minimiza-

CR IP T

tion and orthogonal matching pursuit is proposed. • The convergence under the noise-free case is proved.

• The algorithm is modified to overcome the tubal impulsive noises by employing tensor tubal lp-norm optimization with 1 < p < 2.

• Experimental results show that the proposed methods outperform the

AC

CE

PT

ED

M

putational complexity.

AN US

state-of-the art algorithms in terms of RMSE, PSNR, SSIM and/or com-

1

ACCEPTED MANUSCRIPT

Orthogonal Tubal Rank-1 Tensor Pursuit for Tensor Completion

a College

CR IP T

Weize Suna , Lei Huanga , H.C. Sob,∗, Jiajia Wanga of Information Engineering, Shenzhen University, Shenzhen, Guangdong, China. of Electronic Engineering, City University of Hong Kong, Hong Kong

b Department

Abstract

AN US

This work addresses the issue of tensor completion. The properties of the tensor tubal rank are firstly discussed. It is shown that the tensor tubal rank has similar properties like that of matrix rank derived from SVD. The completion algorithm for the case that the measurements are noise-free or corrupted by Gaussian noise is then proposed based on an orthogonal pursuit on tubal rank-1 tensors. The philosophy behind the devised approach is to relax the problem of tensor tubal

M

rank minimization into tensor Frobenius-norm optimization with a constraint on the maximum number of orthogonal tensors. An iterative procedure which

ED

calculates one orthogonal tensor at each iterative step is then suggested, and the local convergence under the noise-free case is also proved. Furthermore, the proposed method is generalized to the situation where the observations are

PT

corrupted by impulsive noise in a tubal form. To tackle the impulsive noise, we formulate the problem of tensor completion as minimization of tensor tubal `p -norm with 1 < p < 2. An iteratively reweighted procedure is employed to

CE

compute the orthogonal tensors. The algorithms are compared with the stateof-the-art approaches using both synthetic data and real data sets.

AC

Keywords: Tensor completion, tensor SVD, tensor pursuit, rank-one tensor, tensor tubal rank.

∗ Corresponding

author work described in this paper was supported by the National Natural Science Foundation of China (NSFC) under Grant 61501300 and U1501253. 1 The

Preprint submitted to Elsevier

November 26, 2018

ACCEPTED MANUSCRIPT

1. Introduction Recovery of multi-dimensional data under limited number of measurements, which is referred to as matrix completion or tensor completion, is an important

CR IP T

problem which has drawn increasing attention in recent years. This problem is widely seen in various fields such as recommendation systems [1], dimensionality reduction [2], speech signal processing [3][4], MIMO radar [5][6], data mining [7][8], multi-class learning [9], and computer vision [10][11].

The idea of recovering unmeasured entries of a matrix basically relies on the

framework of finding a low-rank matrix that can model the original data, as

AN US

turns out to be the task of minimizing the matrix rank or the `0 -norm of the singular values of the matrix in the general case. However, it is a non-convex

NP-hard problem. One conventional solution is to minimize the matrix nuclear norm or the `1 -norm of the singular values instead of matrix rank, therefore relaxing the non-convex problem to a tractable and convex one [12][13][14]. Based

M

on this concept, a number of algorithms, such as Singular Value Thresholding (SVT) [15], Grassman Manifold based method [16], Singular Value Projection [17] and Fixed Point Continuation with Approximate SVD (FPCA) [18], have

ED

been proposed.

For tensors which are regarded as a generalization of vectors and matrices,

PT

although the process of data recovery is more natural as it employs the intrinsically multi-dimensional structure, the problem is more complicated to solve because the lack of a unique definition of tensor rank. One type of methods is to

CE

use the minimum number of rank-1 tensors from the CANDECOMP/PARAFAC (CP) decomposition [19][20][21]. However, the decomposition is very costly, and sometimes fails to give reasonable result because of the ill-posedness property

AC

[22]. Other researchers minimize the Tucker rank as the tensor rank from the Tucker decomposition to achieve good reconstruction results, but these algorithms require the knowledge of tensor rank [23] or a rank estimate [5] thereby being only valid for special applications. Moreover, some of the tensor completion methods based on the CP or Tucker rank minimization rely on the usage

3

ACCEPTED MANUSCRIPT

of optimization toolbox [24] thus leading to a very high complexity [25]. Other researchers, on the other hand, have proposed to solve it as a low-n-rank tensor recovery problem, which first unfolds an R-D tensor to R matrices in R ways

CR IP T

and then minimizes the sum of the ranks or nuclear norms of all the matrices for data recovery [26]. However, it may fail to exploit the tensor structure with the unfolding operations and thereby leads to a suboptimal procedure [27].

Furthermore, when the tensor is very large, these methods have to perform Singular Value Decomposition (SVD) on several large matrices at each iteration, resulting in considerably high computational complexity [28].

AN US

Recently, a new type of tensor decomposition based on tensor product is

proposed [29]. This factorization method first defines the tensor tubal rank as the tensor rank and then decomposes the tensor into several sub-tensors using the idea of unraveling the multi-dimensional structure by constructing group rings along the tensor tubes. The advantage of such an approach is that the resulting algebra and analysis are very close to those of matrix algebra and

M

analysis [30], and it is shown to be suitable for tensor completion [11]. However, similar to the matrix nuclear norm based methods [15], these algorithms require

ED

a selection of several user-defined parameters. With inappropriate parameters, a severe downgrade of performance will be resulted. Moreover, most existing methods cannot deal with the situation when there are outlier measurements

PT

[31][32]. In this work, we propose a novel tensor completion method based on tensor tubal rank [29][30] minimization as well as orthogonal matching pursuit

CE

[33]. A modification is then applied for tensor completion under tubal impulsive noises by employing tensor tubal `p -norm optimization on the residual between

AC

the measurements and recovered data with p < 2. The rest of the paper is organized as follows. In Section 2, the notations

and required definitions are first presented. In Section 3, our tensor recovery approaches for both random sampling and tubal sampling cases are developed, and the convergence is also proved. In Section 4, a modified algorithm is devised to tackle the situation when the measured entries are corrupted by tubal outliers. Finally, numerical results are provided in Section 5, and conclusions are drawn 4

ACCEPTED MANUSCRIPT

in Section 6.

2. Notations, Definitions and Preliminaries

CR IP T

Before going to the main result, we present the notations and preliminaries

in this paper. We mainly follow the definitions in [29][30], and some notations are newly defined in this section.

Scalars, vectors, matrices and tensors are denoted by italic, bold lower-case, bold upper-case and bold calligraphic symbols, respectively. The transpose and conjugate transpose of a vector or a matrix are written as T and H , and the i × i

AN US

identity matrix is symbolized as Ii . To refer to the (m1 , m2 , . . . , mR ) entry of an

A}(m1 ,m2 ,...,mR ) . The element-wise R-D tensor A ∈ CM1 ×M2 ×···×MR , we use {A multiplication operator between two tensors is defined as , while the symbol t represents the concatenation operator where A = A 1 tr A 2 is obtained by stacking A 2 to the end of A 1 along the rth dimension. C

M

For a compact notation, the m3 th frontal slice of a 3-D tensor A ∈

M1 ×M2 ×M3

is defined as A {m3 } = A (:, :, m3 ) for m3 = 1, 2, · · · , M3 , that

is, A = A {1} t3 A {2} t3 · · · t3 A {M3 } , and an M1 × M1 × M3 identity tensor is

ED

{1} denoted by I M1 M1 M3 with I M1 M1 M3 = IM1 and all the other entries being 0.

Furthermore, the tensor conjugate transpose of a 3-D tensor A ∈ CM1 ×M2 ×M3

PT

is an M2 × M1 × M3 tensor A H which is obtained by taking the conjugate transpose of each of the frontal slices and then reversing the order from the second through the M3 th conjugate transposed frontal slices.

CE

The Fast Fourier Transform (FFT) along the third dimension of all the tubes

e = fft(A A, [ ], 3). Therefore, in the same way, A can of a tensor A is denoted as A

be computed by applying inverse FFT (IFFT) along the third dimension of all

AC

e . We can then define the t-product [29] between two tensors: the tubes of A

Definition 2.1. t-product: The t-product ∗ between two tensors A ∈

RM1 ×M4 ×M3 and B ∈ RM4 ×M2 ×M3 is defined as

C = A ∗ B ∈ RM1 ×M2 ×M3

5

(1)

ACCEPTED MANUSCRIPT

where {m3 } e {m3 } e {m3 }B Ce =A

(2)

CR IP T

for m3 = 1, 2, · · · , M3 , that is, the actual computation of C can be divided into

3 steps: FFT along the third dimension of A and B , matrix multiplications, and

IFFT along the third dimension of Ce. Now the orthogonal tensor is defined as

follows.

Definition 2.2. Orthogonal tensor: A tensor Q ∈ RM1 ×M1 ×M3 is called

AN US

orthogonal if Q H ∗ Q = Q ∗ Q H = I M1 M1 M3 .

(3)

Furthermore, a tensor Q ∈ RM1 ×M2 ×M3 , with M1 ≤ M2 , is called left orthogonal if Q H ∗ Q = I M1 M1 M3 .

Definition 2.3. Block diagonal form: Let A denote the block diagonal

e: matrix of tensor A

M

h {1} e ) = blockdiag( A e A , blockdiag(A

ED

e {2} A

···

i e {M3 } ) ∈ CM1 M3 ×M2 M3 A

(4)

According to this definition, we can see that the block diagonal matrix of A H

PT

is equal to the conjugate transpose of the block diagonal matrix A [30]:

CE

Definition 2.4.

H

AH = A .

(5)

Tensor Singular Value Decomposition (t-SVD):

Analogous to the SVD of a 2-D matrix, the t-SVD of a 3-D tensor X ∈

AC

RM1 ×M2 ×M3 is written as [29]

X = U ∗ S ∗ VH =

Ms X

m=1

Um ∗ Sm ∗ VH m

(6)

where Ms = min{M1 , M2 }, U = U 1 t2 U 2 t2 · · · t2 U Ms ∈ RM1 ×Ms ×M3 , V =

V 1 t2 V 2 t2 · · · t2 V Ms ∈ RM2 ×Ms ×M3 , U m = U (:, m, :) and V m = V (:, m, :) for m = 1, 2, . . . , Ms . The S ∈ RM1 ×M2 ×M3 is defined as a rectangular f-diagonal 6

ACCEPTED MANUSCRIPT

tensor with all frontal slices being diagonal matrices. Note that U and V are

orthogonal tensors. Here we name S m ∈ R1×1×M3 as the mth singular value

V m ) as the left (right) singular vector slices. Actually, tube of tensor X and U m (V

CR IP T

this decomposition can be obtained by applying matrix SVD on all the frontal e to get Ue , Se and Ve and then performing IFFT on them [29]. slices of X

Based on t-SVD, one can now define the tensor ranks [29] as follows.

Definition 2.5. Tensor multi-rank and tensor tubal rank (TTR):

The tensor multi-rank of a tensor X ∈ RM1 ×M2 ×M3 is a vector r ∈ ZM3 ×1 with e . The TTR the m3 th element indicating the rank of the m3 th frontal slice of X

AN US

X ), is defined as the number of non-zero singular value of X , written as TTR(X tubes of S , which is r = max{r}. If a third-order tensor X is of full TTR, it means that r = min{M1 , M2 }.

Definition 2.6. Inner product of tensors: If A and B are third-order tensors of dimensions M1 × M2 × M3 , the inner product between A and B is defined as:

A}(m1 ,m2 ,m3 ) {B B }(m1 ,m2 ,m3 ) {A

(7)

m1 =1 m2 =1 m3 =1

ED

and

M1 X M2 X M3 X

M

< A , B >,

< A , B >=

1 1 H B A) = trace(B < A, B > . M3 M3

(8)

PT

Note that according to the conjugate symmetric property of FFT and (4), it can be easily verified.

CE

Definition 2.7. Tensor Frobenius-norm: The Frobenius-norm of a tenqP PM2 PM3 M1 A||F = A}2(m1 ,m2 ,m3 ) , has the folsor A , defined as ||A m1 =1 m2 =1 m3 =1 {A

AC

lowing relation according to the above definition A||F =< A , A >1/2 = √ ||A

1 e ||F = √ 1 ||A A||F . ||A M3 M3

(9)

Remark 2.1. According to these definitions, by defining the normalized

A} = A /||A A||F we get tensor of A ∈ RM1 ×M2 ×M3 as normalize{A e} = normalize{A

p

7

M3

e /||A e ||F . A

(10)

ACCEPTED MANUSCRIPT

Definition 2.8. Tensor mode-2 multiplication: The tensor mode-2 e ∈ CM1 ×M4 ×M3 and B e ∈ CM4 ×M2 ×M3 is multiplication between two tensors A

e ?B e ∈ CM1 ×M2 ×M3 , where Ce defined as Ce = A

{m3 }

{m3 }

e =A

e {m3 } , which in turn B

e and B e. A

The novel tubal `p -norm is defined as follows.

CR IP T

means that each frontal slice of Ce is a product of corresponding frontal slices of Definition 2.9. Tensor tubal `p -norm: The tensor tubal `p -norm of a

tensor A is computed by first calculating an l2 -norm of all its tubes and then computing the `p -norm of the resulting matrix: M1 X M2 X

AN US

A||t,p = ||A||t,p = ( ||A

{A}p(m1 ,m2 ) )1/p

(11)

m1 =1 m2 =1

A(m1 , m2 , :)||2 or A , ||A A||2,: . Here, A (m1 , m2 , :) is the where {A}(m1 ,m2 ) = ||A (m1 , m2 )th tube of A . Note that (11) can be taken as a tensor form of the `2,p -norm [34].

M

Remark 2.2. According to Definition 2.9, we have

ED

A||t,p = √ ||A

1 e ||t,p . ||A M3

(12)

Remark 2.3. The tensor tubal `2 -norm of a tensor is its Frobenius-norm.

PT

3. Main Result

In this section, we derive the tensor recovery algorithm using 3-D data as an

CE

example. Note that for higher-dimensional data, we can recursively apply the Fourier transform over successive dimensions higher than three [35] and then

AC

unfold the whole tensor into a 3-D one. 3.1. Rank-1 Tensor Pursuit with Random Sampling We first formally formulate the tensor completion problem with random

sampling, which means that any of the measurements in a tensor Y = S + N ∈

RM1 ×M2 ×M3 is sampled independently with a probability q1 ∈ (0, 1), while S

and N refer to the signal and noise terms. The q1 is also named as observation 8

ACCEPTED MANUSCRIPT

percentage. Our objective is to recover S from Y Ω , namely, the observed entries of Y where Ω is an M1 × M2 × M3 tensor with elements being 1 or 0 indicating whether the corresponding entries are observed or not, respectively. Generally

CR IP T

speaking, this problem can be formulated as X) min rank(X

(13)

X

s.t. X Ω = S Ω .

(14)

The tensor rank can be described using the CP [21] or Tucker [23] criterion, or P3 relaxed to a summation of several matrix ranks as i=1 αi rank(Xi ) where Xi

AN US

is the i-way unfolding of the tensor X [26][27]. In this work, we assume that the tensor S has a limited tubal rank r < min{M1 , M2 }. Then it is easy to verify that this tensor can be written as a linear combination of several tubal rank-1 tensors:

M

S = M (l) (θ (l) ) =

X

l∈{L}

θl M l

(15)

Ml : l ∈ {L}} is the set of all the tubal rank-1 tensors with unit where {M

ED

tensor Frobenius-norm, and θl is the magnitude of the l-th rank-1 tensor. The low tubal rank tensor completion problem can then be transformed into the

min ||θ||0 θ

M(θ))Ω = Y Ω = S Ω . s.t. (M

(16) (17)

CE

PT

following constrained minimization:

M(θ))Ω − Y Ω || ≤  where  is a small For the noisy scenarios, (17) becomes ||(M value. It is worth noting that by relaxing the `0 -norm to `1 -norm, the resultant

AC

problem is identical to minimizing the tensor nuclear norm (TNN) [30]. We can then reformulate the problem as Y Ω − (M M(θ))Ω ||2F min ||Y

(18)

s.t. ||θ||0 ≤ R

(19)

θ

9

ACCEPTED MANUSCRIPT

where R indicates the maximum number of tubal rank. According to (9), (18) is equivalent to Y Ω − (M M(θ))Ω ||2F min ||Y

(20)

CR IP T

θ

s.t. ||θ||0 ≤ R

(21)

which can be solved by the orthogonal matching pursuit. In the kth iteration, the tubal rank-1 basis tensors M 1 , M 2 , · · · , M k−1 and their weights θ (k−1) have

already been computed in the former (k−1) iterations. The recovered tensor and

X (k−1) )Ω , the residual term are X (k−1) = M (k−1) (θ (k−1) ) and R (k) = Y Ω − (X

AN US

respectively. When R (k) = 0, which is a zero matrix, the problem is solved and

the rank of the recovered tensor is at most k − 1. On the other hand, for a nonzero R (k) , the next M k can be solved by

Mk ||F = 1} max{< M k , R (k) >: M k = U k ∗ S k ∗ V H k , ||M Mk

(22)

M

where U k , S k , V k are tensors of the dimensions M1 × 1 × M3 , 1 × 1 × M3 , M2 × 1 × M3 , respectively. According to (4), the optimization of the block

ED

diagonal term is substituted by a series of optimization of each block: gk max {< M

gk {m3 } M

{m3 }

g (k) ,R

{m3 }

gk >: M

{m3 }

= Ufk

{m3 }

Sfk

{m3 }

H

(Vfk ){m3 } } (23)

Ufk

PT

Mk ||F = 1. Then we obtain the for m3 = 1, 2, · · · , M3 under the constraint ||M {m3 }

H

gk and (Vfk ){m3 } as the top left and right singular vectors of M

[33] which can be solved using the power method [36] [37], and Sfk

CE

gk corresponding maximum singular value. Once all the M

{m3 }

{m3 }

{m3 }

is the

are computed,

AC

the normalization process of M k can be carried out according to (10) or gk ← M

p

gk /||Sfk ||F M3 M

(24)

and the norm of M k is thus restricted to 1.

Then the new θ (k) is determined by solving the least squares (LS) problem: ˙ (k) θ (k) ||2 min ||y˙ − M θ (k)

10

(25)

ACCEPTED MANUSCRIPT

h i ˙ (k) = m Y Ω }, M Mi )Ω } for ˙ i = vec{(M where y˙ = vec{Y ˙ 1, m ˙ 2. · · · , m ˙ k and m i = 1, 2, · · · , k.

The above two steps can run iteratively until a stopping criterion is satisfied,

CR IP T

and this Tubal Rank-1 Tensor Pursuit (TR1TP) algorithm is summarized in

Algorithm 1. Furthermore, we prove that the algorithm converges and achieves a linear convergence rate in Appendix A, which satisfies:

R(k+1) ||22 ≤ (1 − 1/ min{M1 , M2 })||R Rs (k)||22 ≤ (1 − 1/ min{M1 , M2 })k ||Y Y Ω ||22 . ||R

(26)

Step 1: Initialize R (1) = Y Ω .

AN US

Algorithm 1: Tubal Rank-1 Tensor Pursuit Algorithm for Tensor Completion

g (k) . gi , i = 1, 2, · · · , k − 1, calculate its Fourier domain R Step 2: Given R (k) and M for m3 = 1, 2, · · · , M3

Step 3: Find a pair of top left and right singular vectors (Ufk gk set M

{m3 }

end

= Ufk

{m3 }

M

corresponding singular value Sfk

{m3 }

Sfk

{m3 }

Vfk

{m3 }

g (k) of R

{m3 }

{m3 }

, Vfk

{m3 }

) and its

using the power method, and

.

ED

gk according to (24). Step 4: Normalize M

Step 5: Compute θ (k) as the LS solution of (25).

PT

X (k) )Ω . Step 6: Assign X (k) = M (k) (θ (k) ) and update R (k+1) = Y Ω − (X Step 7: Repeat Steps 2 to 6 until the maximum iteration number R, or

CE

X (k) − X (k−1) ||2 /||X X (k−1) ||2 < , is reached, and output X . ||X Note that this algorithm has only one user-defined parameter, which is the

maximum number of iterations or the tubal rank of the tensor, denoted as R.

AC

However, when the number of iterations is large, the LS solution of (25) needs to compute the inverse of a large matrix, making it computationally unattractive. To overcome this problem, an economic updating scheme which updates α = h i α1 , α2 instead of θ in Algorithm 1 is suggested, and we solve ˙ k α2 ||2 min ||y˙ − x˙ (k−1) α1 − m αk

11

(27)

ACCEPTED MANUSCRIPT

X (k−1) )Ω } and X (k) = α1X (k−1) + instead of (25) by assigning x˙ (k−1) = vec{(X α2M k . As the convergence analysis of the economic updating scheme is similar

CR IP T

to that of Algorithm 1, it is ignored herein. 3.2. Rank-1 Tensor Pursuit with Random Tubal Sampling

Sometimes the received data are sampled in tubal format [30] [38], which means that any one of M1 M2 tubes along the third dimension of a tensor Y ∈

RM1 ×M2 ×M3 is recorded with a probability q1 independently. In this case, if the Y }(m1 ,m2 ,m3 ) is recorded, then all the entries along the (m1 , m2 ) tube entry {Y

AN US

are available. Therefore, it is called ‘random tubal sampling’, and a sampling matrix Ω of dimensions M1 × M2 with elements being 1 or 0 is used to refer to the measured and non-measured tubes of Y instead of Ω in (17). However, for consistency, we follow the notation Y Ω as the observed entries tubes of Y like that in previous sections.

It is easy to see that the Fourier transform operation along the third dimen-

M

sion does not change the structure of the observed entries, therefore we can use e Ω to represent the observed entries of Y e and rewrite (20) as the notation Y

ED

e Ω − (M f (θ))Ω ||2 min ||Y F θ

(28)

which indicates that in this case, we only need to perform the FFT once and

PT

then finish all the computations in the Fourier domain. The modified method, named as the TR1TP for random Tubal sampling (TR1TPT), is shown below

CE

as Algorithm 2. Following the same techniques in Appendix A, Algorithm 2 can also be proved to achieve a linear convergence rate. For compactness of

AC

presentation, it is not shown here.

4. Rank-1 Tensor Pursuit with Random Tubal Outlier In this section, we focus on the problem that the observed data are corrupted

by tubal impulsive noise, that is, some tubes of a tensor Y ∈ RM1 ×M2 ×M3 are

corrupted by impulsive noise while the other tubes are clean or only added by

12

ACCEPTED MANUSCRIPT

Algorithm 2: TR1TPT Algorithm for Tensor Completion. (1)

e Step 2: Initialize R

for m3 = 1, 2, · · · , M3

eΩ . =Y

CR IP T

e. Step 1: Calculate Fourier transform of Y as Y

e (k) and M gi , i = 1, 2, · · · , k − 1, find a pair of top left Step 3: Given R and right singular vectors (Ufk

singular value Sfk

gk set M

{m3 }

end

{m3 }

= Ufk

{m3 }

{m3 }

, Vfk

{m3 }

) and its corresponding

e (k−1) ){m3 } using the power method, and of (R

Sfk

{m3 }

(Vfk

{m3 }

)H .

(k)

e (k) = M f Step 6: Assign X

AN US

gk according to (24). Step 4: Normalize M h i ˙ †˙ ˙ f f ˙ ˙ ˙ e Step 5: Compute θ (k) as θ (k) = M y , where M = e e e m , m . · · · , m k k 1 2 k .

e (k+1) = Y e Ω − (X e (k) )Ω . (θ (k) ) and update R

Step 7: Repeat Steps 3 to 6 until the maximum iteration number R, or e (k) − X e (k−1) ||2 /||X e (k−1) ||2 < , is reached, and output the IFFT ||X

M

e. of X

additive Gaussian noise. This phenomenon has been well elaborated in [32][39].

ED

For ease of presentation, we assume that the observed tensor is sampled by

PT

random tubal sampling, and define the additive noise as E 1 + q2E 2 E = (1 − q2 )E

(29)

where the entries in E 1 are zero or follows a zero mean white Gaussian distri-

CE

bution, and E 2 are outlier tubes. The q2 ∈ (0, 1) in (29) is the probability that the (m1 , m2 ) tube of E is sampled from E 1 and E 2 with probabilities of 1 − q2

AC

and q2 , respectively. The observed data become S + E )Ω Y Ω = (S

(30)

where S has a low tensor tubal rank, and the target is to recover a low-rank tensor X → S from Y Ω . When E 1 = 0, (30) can be formulated as a minimization

13

ACCEPTED MANUSCRIPT

of tubal `0 -norm Y Ω − X Ω ||t,0 min ||Y

(31)

s.t. ||θ||0 ≤ R,

X = M (θ)

CR IP T

θ

where the tubal `0 -norm is used to indicate the location of outliers. However, this is an NP hard problem and is not suitable for the more general case that E 1 6= 0 . Therefore, we relax (31) to tubal `p -norm for p ∈ (1, 2): Y Ω − (M M(θ))Ω ||t,p min f = ||Y θ

(32)

AN US

which can also be solved by utilizing the matching pursuit type algorithm. Note

that when the data are sampled by random tubal sampling and corrupted by random tubal outliers, the completion can be carried out in Fourier domain as in Section 3.2. Here we follow the procedure of economic pursuit method. The equations of data recovery under random sampling and/or using the full weighting vector θ are similar and thus not shown here.

M

e (k−1) In the kth iteration, the recovered tensor and the residual term are X

e (k) = Y eΩ − X e (k−1) , respectively. Subsequently, we relax the Fourier and R Ω

ED

formulation of (32) to a two-step problem that first solves min

R(k) − (U Uk ∗ VH f = ||R k )Ω ||t,p

PT

U k ∗V VH M k =U k

(33)

where U k and V k are tensors of the dimensions M1 × 1 × M3 and 1 × M2 × M3 .

CE

According to (12), (33) is equivalent to min

g H ) f k =(U e k ?V M k Ω

g H ) || e (k) − (U ek ? V f = ||R k Ω t,p

(34)

AC

and then we can solve the problem using an alternately block updating apg H under a given U e k and then Ue k+1 with a fixed proach, that is, updating V k

g H iteratively. Then we can compute V k min

α k =[α1 ,α2 ]

e Ω − α1 (X e (k−1) )Ω − α2 (M f k )Ω ||t,p ||Y

f k according to (24). after normalizing M

14

(35)

ACCEPTED MANUSCRIPT

Assume that in the i-th iteration, the (Ue k )(i) is already computed, and the g g H are both of tubal rankH . As long as U e k and V target is now to determine V k k

g H = H t H t ···t H e = Z 1 t2 Z 2 t2 · · · t2 Z M and V 1, we can write R k 1 2 2 2 2 M2 2

Z m2 − G m2 ? H m2 ||t,p min ||Z

H m2

CR IP T

e where Z m2 ∈ RM1 ×1×M3 and H m2 ∈ R1×1×M3 are the m2 -th slice or tube of R g H , respectively, and change (34) into solving and V k

(36)

e k Ω m and Ω = Ω 1 t2 Ω 2 t2 · · · t2 Ω M with Ω m being the where G m2 = U 2 2 2

same size of Ue k with values 0 or 1. Noting that in (34), the missing tubes of

AN US

Z m2 and G m2 are the same, and we can use the tubes that are only observed, thus reducing the data size and computational complexity.

W ||2,: ∈ RM1 ×1 , we realize By defining W = Z m2 − G m2 ? H m2 and W = ||W that the minimization problem becomes

H m2

M1 X

(W)p(m1 ) =

m1 =1

2 (W)p−2 (m1 ) (W)(m1 )

m1 =1

W||22

D W ||22 = ||D D Z m2 − (D D G m2 ) ? H m2 ||2t,2 = ||D (37)

ED

= ||D

M1 X

M

min ||W||pp =

where D {m3 } = D = W(p−2)/2 for m3 = 1, 2, · · · , M3 . According to (4), (9)

PT

and defining Z D = D Z m2 as well as G D = D G m2 , (36) becomes Z D ) − blockdiag(G GD ) min ||blockdiag(Z

Hm2

H m2 )||22 blockdiag(H

(38)

CE

which is a matrix based optimization problem that can be solved by LS fit. {m } {m3 } † {m3 } GD Therefore the m3 -th slice of H m2 is calculated as H m2 3 = (G ) Z D . The g H can now be computed using such as the Tubal Iteratively Reweighted Least V k

AC

g H and D iteratively until convergence. Squares (TIRLS) method by updating V k (i) (i+1) g g H H can be retrieved following the similar Once (V k ) is obtained, (U k ) procedure of (34) - (38), and terminates when a stopping criterion is reached.

We are now at a position to solve (35). To this end, by defining W α =

{m3 } (p−2)/2 e Ω − α1 (X e (k−1) )Ω − α2 (M f k )Ω , Wα = ||W W α ||2,: and D α Y = Dα = Wα

15

ACCEPTED MANUSCRIPT

for m3 = 1, 2, · · · , M3 , (35) becomes min

α (k) =[α1 ,α2 ]

e )Ω − α1 (D e k−1 )Ω − α2 (D e k )Ω ||2 Dα Y Dα X Dα X ||(D

(39)

CR IP T

e )Ω , (D e k−1 )Ω Dα Y Dα X and then α (k) is calculated using LS after vectorizing (D e k )Ω . Dα X and (D

α(k) )ini = [1, 1], it is updated iteratively with With an initialization of (α

{m } D α 3 until a stopping criterion is satisfied. The proposed TIRLS based eco-

nomic TR1TPT algorithm (TIRLS-TR1TPT-Econ) is summarized as Algorithm

AN US

3.

Algorithm 3: TIRLS-TR1TPT-Econ Algorithm for Tensor Completion e and set R e (1) = Y eΩ . Step 1: Calculate Fourier transform of Y as Y (k) (k−1) g H following e e Step 2: Given R and X , calculate Ufk and V k g H . gk = Ufk ? V (34) - (38) and set M k

gk according to (24). Step 3: Normalize M e Step 5: Assign X

(k)

M

Step 4: Compute α (k) according to (39). e = α1X

(k−1)

f k and update R e + α2M

(k+1)

e Ω − (X e =Y

(k)

)Ω .

Step 6: Repeat Steps 2 to 4 until the maximum iteration number R, or

PT

e. of X

ED

e (k) − X e (k−1) ||2 /||X e (k−1) ||2 < , is reached, and output the IFFT ||X

CE

5. Numerical Results In this section, we show the numerical results, based on the synthetic data as

well as real image and video data, of the proposed algorithms and state-of-the-

AC

art algorithms. The first subsection considers the standard tensor completion problem while the second one tackles the problem in the presence of tubal outliers. All the experiments are repeated 100 times to obtain the average performance.

16

ACCEPTED MANUSCRIPT

5.1. Tensor Completion under Random Tubal Sampling and Random Sampling In the first test, we study the convergence of the proposed ‘TR1TP’ and ‘TR1TPT’ methods as well as their economic variants, named as ‘TR1TP-Econ’

CR IP T

and ‘TR1TPT-Econ’. The element-wise Root Mean Square Error (RMSE), deX −S S ||F /||S S ||F , where X and S are the recovered and original tensors, fined as ||X is used to evaluate the performance of data recovery. The dimensions of the synthetic data tensor S is 100 × 100 × 30 of tubal rank 5, and its elements are

sampled under the rules of random tubal sampling with two different sampling

rates q1 = 40% and 80%. The empirical RMSEs versus the iteration number

AN US

R in noise-free and Gaussian noise situations are plotted in Figures 1 and 2, respectively. Note that, for the noisy case, the observed entry tubes are corrupted by an additive white Gaussian noise (AWGN) and the signal-to-noise ratio (SNR) is 30dB. It is shown that for all the methods, as the number of

iterations increases, the recovery performance improves. When R is larger than a certain threshold R0 , such as R > R0 = 30 at q1 = 40% under noise-free

M

case, the proposed algorithms converge, and have good performance for a wide range of R. Therefore, we can simply set R to a large number to get acceptable

ED

results and avoid performance downgrade. For the noisy scenario, after several iterations, the performance of the proposed algorithms degrades slightly as the number of iterations increases, showing that in this case, the new recovered

PT

results from new iterations are actually noise. It is worth noting that for both the noise-free and noisy cases, as the sampling rate increases, the number of

CE

iterations required to converge increases, indicating that R should be chosen in proportional to q1 . Furthermore, due to the de-noising effect of the FFT operation along the 3rd dimension, the recovery performance of ‘TR1TPT’ and

AC

‘TR1TPT-Econ’ methods are slightly better than ‘TR1TPT’ and ‘TR1TPTEcon’, showing the attractiveness of the former tensor completion technique under the tubal sampling model. In the second test, the proposed methods are compared with the Tensor

Nuclear Norm algorithm based on t-SVD (‘t-SVD-TNN’) [30], the economic Rank-1 Tensor Pursuit (‘R1TP-Econ’) method [40] as well as the Fixed Point 17

ACCEPTED MANUSCRIPT

-1

0 TR1TP TR1TP-Econ TR1TPT TR1TPT-Econ

-5

-4

RMSE (dB)

RMSE (dB)

-3

-5 -6

CR IP T

-2

-10

-15

-7 -20 -8

0

20

40

Number of iterations

-25

AN US

-9

60

0

20

40

60

Number of iterations

Figure 1: Average RMSE versus iteration number in noise-free case. Left: q1 = 40%; Right: q1 = 80%.

M

Low-Rank Tensor Recovery (‘FP-LRTC’) approach [26]. The random sampling model is employed thus only the ‘TR1TP’ and ‘TR1TP-Econ’ approaches are applicable. We change the parameter M2 from 80 to 180 and set q1 = 90%, with

ED

the other parameters being the same as those in Figure 1. All the algorithms will terminate when the normalized difference between adjacent iterations is small

PT

X (k) − X (k−1) ||2 /||X X (k−1) ||2 <  where  = 10−6 , or the N th enough, that is, ||X iteration is reached. According to Figure 1, we set N = R = 60 for the proposed and ‘R1TP-Econ’ approaches. For the other algorithms, N = 100 is used. For

CE

the t-SVD-TNN approach, the shrinkage parameter is set to be τ = 0.025P ,

Y Ω ||2 is the `2 -norm of the observed entries, and decreases by a where P = ||Y

AC

factor 0.95 after each iteration. The results of RMSE and computational time in seconds are shown in Figures 3 and 4. In this test, the ‘t-SVD-TNN’ algorithm gives the best performance for small data size, but it is much more complicated than the proposed ‘TR1TP-Econ’ one. The two proposed algorithms give the best performance when the size of tensor is sufficiently large, verifying the efficiency of the economic process according to Figure 4. On the other hand, the

18

ACCEPTED MANUSCRIPT

-1

0 TR1TP TR1TP-Econ TR1TPT TR1TPT-Econ

-2

-4

-5 -6

-8 -10 -12

-7

-14

-8

-16

0

10

20

30

Number of iterations

-18

AN US

-9

CR IP T

-6

-4

RMSE (dB)

RMSE (dB)

-3

-2

40

0

10

20

30

40

Number of iterations

Figure 2: Average RMSE versus iteration number in Gaussian noise. Left: q1 = 40%; Right: q1 = 80%.

to the modeling error.

M

‘R1TP-Econ’ and ‘FP-LRTC’ methods provide unsatisfied results, partly due

The effect of Gaussian noise on the recovery performance is now examined.

ED

The empirical RMSEs of various techniques versus SNR are shown in Figure 5. We employ the same parameters as those in Figure 2 except q1 = 70%

PT

and R = 25. It is observed that the ‘t-SVD-TNN’ and proposed methods give similar recovery performance, and they can retrieve the data tensor correctly when the SNR is sufficiently large. The average computational times

CE

of the ‘TR1TP’, ‘TR1TP-Econ’, ‘t-SVD-TNN’, ‘R1TP-Econ’ and ‘FP-LRTC’ are 5.1025s, 1.4534s, 11.7652s, 1.4645s and 5.5055s, respectively, indicating that

AC

the ‘TR1TP-Econ’ is the most computationally attractive one among all the methods. To evaluate the applicability of real-world problems, we test our methods

and other state-of-the-art algorithms on color image inpainting, that is, to fill in the missing pixel values of a partially observed image. We employ the peak SNR (PSNR), structural similarity index (SSIM) [41] and computational time to

19

ACCEPTED MANUSCRIPT

-5

-10

TR1TP TR1TP-Econ t-SVD-TNN R1TP-Econ FP-LRCT

-20

-25

-30

-35

90

100

110

120

130

140

150

160

170

180

AN US

-40 80

CR IP T

RMSE (dB)

-15

M2

Figure 3: Average RMSE versus M2 .

evaluate these approaches, and use the Berkeley Segmentation database [42] in

M

this experiment. This database contains a wide variety of natural scenes and has a total number of 100 test images with dimensions 321×481×3. For each image,

ED

we discard the last column and row to get the data tensor Y 1 ∈ [0, 255]320×480×3 , reshape it to a higher-order tensor of size 16 × 20 × 20 × 24 × 3, and then unfold the 3rd and 4th dimensions into the 1st and 2nd ones to form a new tensor Y 2

PT

of the same dimensions as Y 1 . This Y 2 can be regarded as a ’detail-background’ representation of the original image Y 1 [43], and the final inpainting results of all the algorithms are the average of recovered values of both tensors.

CE

The image recovery performance in noise-free environment is studied at first

and the random tubal sampling model is used. The proposed ‘R1TPT-Econ’

AC

algorithm is employed and compared with ‘t-SVD-TNN’, ‘R1TP-Econ’ and ‘FPLRTC’ methods at 20% observed tubes of all 100 images. The maximum number of iterations is N = R = 100 for all the algorithms, and the other parameters are the same as those in the previous test. The PSNR results of all the images are shown in Figure 6. The average PSNRs of ‘R1TPT-Econ’, ‘t-SVD-TNN’ and ‘R1TP-Econ’ methods are 22.0965dB, 22.1316dB and 21.5329dB, and the 20

ACCEPTED MANUSCRIPT

45

35

Time (s)

30 25 20 15 10 5

90

100

110

120

130

140

150

160

170

180

AN US

0 80

CR IP T

TR1TP TR1TP-Econ t-SVD-TNN R1TP-Econ FP-LRCT

40

M2

Figure 4: Average computational time versus M2 .

corresponding average SSIM values are 0.8991, 0.8970 and 0.8848, verifying that

M

the ‘t-SVD-TNN’ and proposed methods give the best image recovery performance. The ‘FP-LRTC’ approach, on the other hand, fails to give reasonable

ED

results.

In the second experiment of image inpainting, we test the recovery performance in the presence of noise. A small disturbance or noise N is added to

PT

the observed entries and N follows an independent and identically uniform distribution between [−25, 25]. Note that the values of noise-free tensor S vary from 0 to 255. The parameters are the same as those in previous test except

CE

the sampling rate q1 varies from 5% to 35%. The average PSNR and SSIM results of all the images are shown in Figure 7, and one typical results of 5 im-

AC

ages, namely, ‘Doll’, ‘Elephant’, ‘Dock’, ‘Rocks’ and ‘Landscape’, at q1 = 15% are shown in Figure 8. The average computational times of ‘TR1TPT-Econ’, ‘t-SVD-TNN’, ‘R1TP-Econ’ and ‘FP-LRTC’ are 43.97s, 86.47s, 18.17s, 0.9618s and 505.21s, respectively. It is shown that for small observation percentage, the ‘R1TPT-Econ’ approach gives the best performance with moderate complexities comparing to the state-of-the-art ones, indicating the robustness of the 21

ACCEPTED MANUSCRIPT

0 -2

CR IP T

RMSE (dB)

-4 -6 -8 -10

TR1TP TR1TP-Econ t-SVD-TNN R1TP-Econ FP-LRCT

-12

-16

0

5

10

AN US

-14

15

20

25

30

SNR (dB)

Figure 5: Average RMSE versus SNR

M

proposed methods against this severe environment.

5.2. Tensor Completion under Random Tubal Outliers In this section, the proposed ‘TIRLS-TR1TPT-Econ’ approach is compared

ED

with the state-of-the-art algorithms under the situation where the observed random tubal sampled data are corrupted by tubal impulsive noise. The methods

PT

to be compared with are the ‘IRLS-t-SVD’ [44], ‘IRLS-PARAFAC’ [45], and ‘t-SVD-TNN’ based on `2 -norm minimization. At first, we compare the completion performances of various methods under

CE

different outliers levels for synthetic data, which is a 3-D tensor of dimensions 100 × 100 × 20 and tubal rank 5. Only 60% of the tubal entries are observed,

AC

and 10% of the observed entries are corrupted by tubal outliers. The outliers

E 2 in (29) are generated using a zero mean Gaussian distribution as well as

S ||22 /||E E 2 ||22 ranging uniform distribution, with a signal-to-outlier ratio (SOR) ||S from −10dB to 10dB. For all the IRLS based methods, the p in `p -norm minimization is set to 1.5. The maximum number of iterations is N = R = 20 for the ‘IRLS-t-SVD’ and proposed methods and N = 100 for the others. As

22

ACCEPTED MANUSCRIPT

TR1TP-Econ t-SVD-TNN R1TP-Econ FP-LRCT

30

20 15 10 5 0

0

5

10

15

20

25

30

35

Image ID 30

PSNR(dB)

25 20 15 10

0 50

55

60

65

40

45

AN US

5

CR IP T

PSNR(dB)

25

70

75

80

85

90

95

50

100

Image ID

Figure 6: PSNR results of all 100 images.

M

the ‘IRLS-t-SVD’ approach is a hard threshold method which requires prior knowledge of the tubal rank of the data, the real tubal rank 5 is chosen. The remaining parameters are the same as those in Figure 5. The mean and stan-

ED

dard deviation (std) of RMSE results are shown in Table 1. It is shown that the proposed approach gives the best completion performance. When the SOR is large, the l2 -norm based ‘t-SVD-TNN’ method can achieve an acceptable recov-

PT

ery result, but it fails when the SOR is small. The PARAFAC-type approach, however, fails due to the modeling error.

CE

The images in Berkeley Segmentation database are now employed for image inpainting in the presence of outliers. The images are sampled tubal-wised at q1 = 40% and 60%, and q2 = 10% of the measured tubes are corrupted by salt-

AC

and-pepper noises. The PSNR and SSIM results of the proposed algorithm with p = 1.7 under different numbers of iterations or tensor components R are shown

in Figure 9 to illustrate the convergence of the methodology. Similar to the cases in Figures 1 and 2, it is observed that for a larger observation percentage q1 , a larger R is required. Nevertheless, a smaller R is suggested for tensor

23

ACCEPTED MANUSCRIPT

24 0.9 22

0.85

CR IP T

0.8

20

18

SSIM

PSNR (dB)

0.75

16

0.7 0.65 0.6

14

TR1TPT-Econ t-SVD-TNN R1TP-Econ FP-LRCT

0.1

0.2

0.5 0.45 0.4

AN US

12

10

0.55

0.3

observation percentage

0.1

0.2

0.3

observation percentage

Figure 7: PSNR and SSIM versus observation percentage in the presence of noise.

completion under impulsive noise comparing to the noise-free case under the

M

similar parameter setting. In the remaining experiments, we set R = 15. Then the ‘t-SVD-TNN’ and ‘IRLS-t-SVD’ methods are included for comparison and

ED

for the IRLS based methods, p = 1.7 is used. The maximum number of iterations for the ‘t-SVD-TNN’ and ‘IRLS-t-SVD’ methods are 100 and 50. Furthermore, for ‘IRLS-t-SVD’ approach, the cases of tubal rank varying from 5 to 15 are

PT

tested and only the best result is used. The average PSNR and SSIM results of all the images are shown in Figure 10, and the average computational times of the ‘TIRLS-TR1TPT-Econ’, ‘t-SVD-TNN’ and ‘IRLS-t-SVD’ are 49.56s, 40.93s

CE

and 543.95s, respectively. It is seen that the ‘TIRLS-R1TPT-Econ’ approach can give the best recovery performance for both q2 = 10% and 20%, and the

AC

computational complexity is much smaller than the ‘IRLS-t-SVD’ approach, indicating the efficiency of the proposed method against this severe environment. The ‘t-SVD-TNN’ approach, on the other hand, fails to recover the images correctly. Finally, the video inpainting problem is tackled and the algorithms are evaluated on the widely used YUV Video Sequences [46]. In the experiment, we test 24

ACCEPTED MANUSCRIPT

TR1TPT-Econ

t-SVD-TNN

R1TP-Econ

FP-LRCT

Landscape

AN US

Rocks

Dock

CR IP T

Elephant

Doll

original image observed image

M

Figure 8: Typical recovery result of 5 images.

ED

the methods using three videos, namely, ‘Akiyo’, ‘Container’ and ‘Hall’, and the dimensions of the videos are 144×176×3×30. As the data are of 4-D, in order to apply the t-SVD based and proposed methods, we simply fold the 4th dimension

PT

into the second one. Furthermore, similar to the image cases, the 4-D video is reshaped to a higher-order tensor of dimensions 12 × 12 × 11 × 16 × 3 × 5 × 6, and

CE

then unfold the 3rd and 6th dimensions into the 1st one, and the 4th and 7th dimensions into the 2nd one to form a new tensor of dimensions 660 × 1152 × 3, and the final results of all the algorithms are the average of recovered values

AC

of both tensors. The parameters are the same as those in previous one except q1 = 50% and q2 = 10%. The average PSNR and SSIM results of all the 3 videos are given in Table 2, and one typical recovered frame of the ‘TIRLSTR1TPT-Econ’, ‘t-SVD-TNN’ and ‘IRLS-t-SVD’ approaches are displayed in Figure 11 at p = 1.5. We observe that the ‘TIRLS-TR1TPT-Econ’ scheme also

25

ACCEPTED MANUSCRIPT

Table 1: Completion results of different methods (RMSE (100 times mean ± std), -- means fail to give reasonable result). Uniform noise as outliers

CR IP T

SOR(dB) Algorithm

10

5

0

-5

-10

TIRLS-TR1TPT-Econ

4.58±0.28

5.72±0.54

7.07±1.11

8.22±1.47

10.31± 2.21

IRLS-t-SVD

11.71±1.51

13.36±1.73

15.70±1.78

16.84±2.01

18.37±2.21

IRLS-PARAFAC

98.57±1.22

98.67±1.17

98.44±1.48

98.93±1.09

98.78±1.21

t-SVD-TNN

17.61±0.27

32.03±0.45

56.52±0.79

--

--

Gaussian noise as outliers

10

5

0

-5

-10

TIRLS-TR1TPT-Econ

3.53± 0.20

4.34± 0.29

5.40± 0.44

6.82± 1.34

8.14± 1.94

IRLS-t-SVD

10.53±1.83

11.85±1.69

13.07± 1.89

14.96± 1.61

16.93±1.97

IRLS-PARAFAC

98.47±1.41

98.24±1.54

98.38±1.40

98.84±0.98

98.76±1.17

t-SVD-TNN

8.39±0.25

15.56±0.26

28.43±0.44

50.58±0.65

--

M

Algorithm

AN US

SOR(dB)

has good recovery performance and low complexity for a wide range of p, while

ED

the ‘IRLS-t-SVD’ method fails at p = 1.1. Generally speaking, the proposed ‘TIRLS-TR1TPT-Econ’ approach gives the best recovery and is the most computationally attractive one among all the algorithms, confirming its robustness

PT

and high efficiency against random tubal outliers.

CE

6. Conclusion

In this paper, we devise several data recovery methods for the tensor comple-

AC

tion problem using the idea of tensor tubal rank and greedy orthogonal tensor pursuit. The TR1TP and TR1TPT methods are first developed under the situation when the measured entries are randomly sampled or tubal randomly sampled. Using the idea of unraveling the multi-dimensional structure by performing FFT along the tensor tubes, the recovered orthogonal multi-dimensional tensors can be retrieved. Additionally, the TR1TP approach is proved to be con-

26

0.94

21.5

0.92

21

0.9

20.5

0.88

20

0.86

19.5 19

CR IP T

22

SSIM

PSNR (dB)

ACCEPTED MANUSCRIPT

0.84 0.82

TIRLS-TR1TPT-Econ (q1 = 60%)

18.5

0.8

TIRLS-TR1TPT-Econ (q1 = 40%)

18

0

5

10

15

iteration number

0.76

AN US

17.5

0.78

20

0

5

10

15

20

iteration number

Figure 9: PSNR and SSIM versus number of iterations.

verged linearly. Furthermore, when there are tubal outliers in the measurements,

M

the algorithm is modified by combining the idea of tubal `p -norm optimization. Experimental results show that the proposed methods outperform the state-

complexity.

ED

of-the art algorithms in terms of RMSE, PSNR, SSIM and/or computational

PT

7. Appendix A

We show that the proposed TR1PT algorithm has a linear convergence rate.

CE

Before proving the convergence, some useful properties will be presented first. Mi )Ω >= 0 and < R (k+1) , M i >= 0 for i = Property A.1. < R (k+1) , (M

AC

1, 2, · · · , k.

Proof. Note that θ (k) is the optimal solution of problem (25). As a result, we

M(i) (θ (i) ))Ω , (M Mi )Ω >= 0 for i = 1, 2, · · · , k, where (M M(i) (θ (i) ))Ω have < Y Ω −(M

X (k) )Ω and X (k) = M (k) (θ k ), is defined in (15). As long as R (k+1) = Y Ω − (X

Mi )Ω >= 0. Furthermore, the we reach to the conclusion that < R (k+1) , (M

unobserved entries in R(k+1) are all set to 0 according to the definition, leading 27

ACCEPTED MANUSCRIPT

22

0.95

21

0.9

20 0.85

CR IP T

0.8

18

SSIM

PSNR (dB)

19

17

TIRLS-TR1TPT-Econ (q =10%) 2

0.75

t-SVD-TNN (q2 =10%)

IRLS-t-SVD (q2 =10%)

16

TIRLS-TR1TPT-Econ (q2 =20%)

0.7

t-SVD-TNN (q2 =20%)

15

IRLS-t-SVD (q2 =20%)

0.65

13 0.2

0.3

0.4

0.5

observation percentage

AN US

14

0.6 0.2

0.6

0.3

0.4

0.5

0.6

observation percentage

Figure 10: PSNR and SSIM versus observation percentage under salt-and-pepper noise.

Mi ) >= 0 for i = 1, 2, · · · , k. to < R (k+1) , (M

M

This property indicates that R (k+1) is perpendicular to the previously re-

covered M i .

ED

˙ (k) in (25) has a full column rank. Property A.2. For R (k) 6= 0, M

Proof. According to the definition, when R (k) 6= 0, we have R (i) 6= 0 for all i = 1, 2, · · · , k − 1, otherwise the stopping criterion of the proposed method

PT

˙ (1) 6= 0. Hence, the conclusion is violated. Since R (1) = Y Ω 6= 0, we have M ˙ (k−1) , k ≥ 2, holds for k = 1. We now assume that the property holds for M

CE

˙ (k) . In particular, as M ˙ (k−1) has full column and show that it also holds for M

AC

˙ (k) does not have full column rank, then rank, suppose that for contradiction, M h iT ˙ (k−1)α or ˙k=M there exists an α = α1 α2 · · · αk−1 that satisfies m Mk )Ω = (M

i=k−1 X i=1

Mi )Ω . αi (M

According to Property A.1., we have < R (k+1) , M k

R(k)

R

, M k >= 0, that is, <

R (k) , M

k

(40) >= 0 and <

>= 0 in (22), indicating that the maximum

singular value of all the frontal slices of R(k) is 0, leading to R(k) = R(k) = 0. 28

ACCEPTED MANUSCRIPT

Table 2: PSNR(dB), SSIM and computational complexity (seconds) of different methods under different p. p

Algorithm

1.1

1.3

1.5

1.7

TIRLS-TR1TPT-Econ

20.44

21.86

22.05

IRLS-t-SVD

11.82

17.92

19.98

t-SVD-TNN

16.77

SSIM

p 1.1

1.3

1.5

TIRLS-TR1TPT-Econ

0.900

0.923

0.927

IRLS-t-SVD

0.505

0.815

0.875

t-SVD-TNN

20.21

19.83

19.50

18.67

1.7

1.9

0.899

0.875

0.863

0.838

AN US

Algorithm

1.9

CR IP T

PSNR(dB)

0.771

Complexity(seconds)

p

1.1

1.3

1.5

1.7

1.9

TIRLS-TR1TPT-Econ

129.2

129.8

129.5

127.6

129.6

IRLS-t-SVD

950.8

1358.9

1359.3

1391.1

1245.3

t-SVD-TNN

M

Algorithm

395.9

ED

˙ (k) has a full column This contradicts the fact that R (k) 6= 0. Therefore, M rank.

R(k+1) || and Now we can build the relation between two successive residuals ||R

PT

(k) R(k) || by defining θk−1 ||R , namely, the value of θ(k) at the (k − 1)-th iteration,

CE

as 0. Setting θ (k) = θ (k−1) + ι(k) , it follows from (18) that R(k) − ι(k) = arg min ||R

ι

i=1

(41)

AC

ι

k X ˙ (k) ι(k) ||2 Mi )Ω ||22 = arg min ||r(k) − M {ι(k) }(i) (M 2

˙ (k) is defined in (25). Since M ˙ (k) has full R(k) )} and M where r(k) = vec{(R

column rank according to Property A.2, we have the solution ˙ (k) )T M ˙ (k) )−1 (M ˙ (k) )T r(k) . ι(k) = ((M

29

(42)

ACCEPTED MANUSCRIPT

Observed frame

TIRLS-TR1TPT-Econ

IRLS-t-SVD

t-SVD-TNN

AN US

Hall

CR IP T

Container

Akiyo

Original frame

Figure 11: Typical recovered frame of 3 videos.

Defining

M(k) (ι(k) ))Ω = L (k) = (M

k X i=1

Mi )Ω {ι(k) }(i) (M

(43)

M

we have X k = X k−1 + L (k) and R (k+1) = R (k) − L (k) .

R(k+1) ||22 = ||R R(k) ||22 − ||L L((k)) ||22 and ||L L(k) ||22 ≥< ||R

R (k) , M k >.

ED

Property A.3.

According to Property A.1 and (43), it is clear that < R (k+1) , L (k) >= 0.

Hence, we obtain

PT

R(k+1) ||22 = ||R R(k) − L (k) ||22 ||R

AC

CE

R(k) ||22 − 2 < R (k) , L (k) > +||L L(k) ||22 = ||R R(k) ||22 − 2 < R (k+1) + L (k) , L (k) > +||L L(k) ||22 = ||R R(k) ||22 − ||L L(k) ||22 . = ||R

(44)

L(k) ||22 is bounded as follows. For any R (k) 6= 0, according to Property The ||L

A.1, we have h ˙ (k) )T r(k) = 0 (M

0

···

0

< R (k) , M k >

iT

.

(45)

From the definition, M i for i = 1, 2, · · · , k is normalized to 1. As long as the

˙ (k) , denoted as m Mk )Ω }, contains only part of ˙ k = vec{(M k-th column of M 30

ACCEPTED MANUSCRIPT

˙ (k) = QU be the QR ˙ k || ≤ 1. Let M elements of M k , it has the property ||m

˙ k where QT Q = IM , MΩ = M1 M2 M3 q1 , and U ∈ Rk×k is factorization of M Ω an upper triangular matrix. According to [33], the (k, k) element of U, denoted

CR IP T

as Ukk , has the property ˙ k || ≤ 1. |Ukk | ≤ ||m

(46)

Combining (42) and (46), we have

˙ (k) )T r(k) )k =< R (k) , M k > /|Ukk |2 ≥< R (k) , M k > . {ι(k) }(k) = ((UT U)−1 (M

(47)

L(k) ||22 =< R (k) − R (k+1) , ||L

AN US

L(k) ||22 can be written as Furthermore, ||L

k X Mi )Ω > {ι(k) }(i) (M i=1

Mk )Ω > + < R (k) , =< R (k) , ιk,k (M (k)

}(k) < R

≥< R (k) , M k > .

, Mk >

i=1

Mi )Ω > − < R (k+1) , {ι(k) }(i) (M

M

= {ι

(k)

k−1 X

k X Mi )Ω > {ι(k) }(i) (M i=1

(48)

ED

Now the convergence of the proposed algorithm is investigated. According to the definition, we have

M3 {m3 } {m } 1 1 X g (k) gk 3 , R < R (k) , M k >= M3 M3 1

PT

< R (k) , M k > =

CE

=

M3 {m } {m3 } p {m3 } 1 X g H ){m3 } , R fk 3 > (Sfk M3 /||Sfk ||F ) < Ufk (V k M3 1

(49)

{m }

AC

3 {m3 } {m3 } g (k) where Sfk is the dominant singular value of R defined in (23), Ufk g H ){m3 } are the corresponding singular vectors, and ||S fk || is defined in and (V k

(24). It follows that Sfk

{m3 }

{m3 }

{m } g H ){m3 } , R fk 3 >≥ (V k

= < Ufk r {m } fk 3 ||2 /min{M1 , M2 } ≥ ||R 2

31

r

fk ||R

{m3 }

fk ||22 /rank(R

{m3 }

) (50)

ACCEPTED MANUSCRIPT

and ||Sfk ||22 =

M3 X

(Sfk

{m3 }

m3 =1

)2 ≥

{m }

fk 3 ||2 fk ||2 ||R ||R 2 2 = . min{M , M } min{M , M2 } 1 2 1 =1

M3 X

m3

(51)

CR IP T

Combining the equations (49)-(51) yields g (k) ||2 /(M min{M , M }) = ||R R(k) ||22 /min{M1 , M2 } (52) < R (k) , M k >≥ ||R 3 1 2 2

which in turn leads to

R(k+1) ||22 = ||R R(k) ||22 − ||L L(k) ||22 ≤ ||R R(k) ||22 − < R (k) , M k > ||R Rk ||22 . ≤ (1 − 1/ min{M1 , M2 })||R

(53)

AN US

Hence, the TR1TP algorithm converges linearly and is not dependent on

the length of the third dimension M3 . Furthermore, by using the fact that R(1) ||2 = ||Y Y Ω ||2 , we get (26). ||R References

M

References

[1] N. Boumal and P. A. Absil, “RTRMC: A Riemannian trust-region method

ED

for low-rank matrix completion,” Proceedings of International Conference on Neural Information Processing Systems 24 (NIPS 2011), pp.406-414, Dec. 2011, Granada, Spain.

PT

[2] N. Linial, E. London and Y. Rabinovich, “The geometry of graphs and some of its algorithmic applications,” Combinatorica, vol.15, no.2, pp.215-

CE

245, Jun. 1995.

[3] Y. Pang, S. Wang and Y. Yuan, “Learning regularized LDA by cluster-

AC

ing,” IEEE Transactions on Neural Networks and Learning Systems, vol.25, no.12, pp.2191-2201, Apr. 2014.

[4] M.J. Taghizadeh, R. Parhizkar, P.N. Garner, H. Bourlard and A. Asaei, “Ad hoc microphone array calibration: Euclidean distance matrix completion algorithm and theoretical guarantees,” Signal Processing, vol.107, pp.123-140, Feb. 2015. 32

ACCEPTED MANUSCRIPT

[5] L.T. Huang, A.L.F.De Almeida and H.C. So, “Target estimation in bistatic MIMO radar via tensor completion,” Signal Processing, vol.120, pp.654659, Mar. 2016.

CR IP T

[6] X. Hu, N. Tong, J. Wang, S. Ding and X. Zhao, “Matrix completion-based MIMO radar imaging with sparse planar array,” Signal Processing, vol.131, pp.49-57, Feb. 2017.

[7] T. G. Kolda and J. Sun, “Scalable tensor decompositions for multi-aspect data mining,” 2008 Eighth IEEE International Conference on Data Mining,

AN US

pp.363-372, Dec. 2008, Pisa, Italy.

[8] J. Sun, S. Papadimitriou, C.-Y. Lin, N. Cao, S. Liu and W. Qian, “MultiVis: Content-based social network exploration through multi-way visual analysis,” SIAM International Conference on Data Mining (SDM2009), pp.1063-1074, Sparks, Nevada, USA, Apr. - May, 2009.

M

[9] G. Obozinski, B. Taskar and M. I. Jordan, “Joint covariate selection and joint subspace selection for multiple classification problems,” Statistics

ED

Computing, vol.20, no.2, pp.231-252, 2010. [10] J. Liu, P. Musialski, P. Wonka and J. Ye, “Tensor completion for estimating missing values in visual data,” IEEE Transactions on Pattern Analysis and

PT

Machine Intelligence, vol.35, no.1, pp.208-220, 2013. [11] Z. Zhang, G. Ely, S. Aeron, N. Hao and M. Kilmer, “Novel methods for mul-

CE

tilinear data completion and de-noising based on tensor-SVD,” IEEE Conference on Computer Vision and Pattern Recognition (CVPR)), pp.3842-

AC

3849, Jun. 2014.

[12] M. Fazel, H. Hindi and S. Boyd, “A rank minimization heuristic with application to minimum order system approximation,” Proceedings of the 2001 American Control Conference, vol.6, no.2, pp. 4734-4739, Arlington, USA, 25-27 Jun. 2001.

33

ACCEPTED MANUSCRIPT

[13] B. Recht, M. Fazel and P.A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” SIAM Review, vol.52, pp.471-501, 2010.

CR IP T

[14] E.J. Cand` es and B. Recht, “Exact matrix completion via convex optimiza-

tion,” Foundations of Computational Mathematics, vol.9, no.6, pp.717, 2009.

[15] J.F. Cai, E.J. Cand` es and Z. Shen, “A singular value thresholding algo-

rithm for matrix completion,” SIAM Journal on Optimization, vol.20, no.4,

AN US

pp.1956-1982, Jan. 2010.

[16] R.H. Keshavan and S. Oh, “A gradient descent algorithm on the Grassman manifold for matrix completion,” arXiv preprint arXiv:0910.5260, 2009. [17] P. Jain, R. Meka and I.S. Dhillon, “Guaranteed rank minimization via singular value projection,” Proceedings of International Conference on Neural

Vancouver, Canada.

M

Information Processing Systems 23 (NIPS 2010), pp. 937-945, Dec. 2010,

ED

[18] S.Q. Ma, D. Goldfarb and L. F. Chen, “Fixed point and Bregman iterative methods for matrix rank minimization,” Mathematical Programming,

PT

vol.128, pp.321-353, 2011.

[19] T.G. Kolda and B.W. Bader, “Tensor decompositions and applications,”

CE

SIAM Review, vol.51, no.3, pp.455-500, Sep. 2009. [20] C. Mu, B. Huang, J. Wright and D. Goldfarb, “Square deal: Lower bounds

AC

and improved relaxations for tensor recovery,” Proceedings of the 31st International Conference on Machine Learning (ICML), wol.32, no.2, pp.7381, 22-24 Jun. 2014, Beijin, China.

[21] P. Jain and S. Oh, “Provable tensor factorization with missing data,” Proceedings of International Conference on Neural Information Processing Systems 27 (NIPS 2014), pp.1431-1439, Dec. 2014, Montr´ eal, Canada.

34

ACCEPTED MANUSCRIPT

[22] V. De Silva and L.H. Lim, “Tensor rank and the ill-posedness of the best low-rank approximation problem,” SIAM Journal on Matrix Analysis Applications, vol.30, no.3, pp.1084-1127, 2008.

CR IP T

[23] M. Filipovi´ c and A. Juki´ c, “Tucker factorization with missing data with application to low-n-rank tensor completion,” Multidimensional Systems and Signal Processing, vol.26, no.2, pp.677-692, Jul. 2015.

[24] T.N. Thieu, H.J. Yang, T.D. Vu and S.H. Kim, “Recovering incomplete data using Tucker model for tensor with low-n-rank,” International Journal

AN US

of Contents, vol.12, no.23 pp.22-28, Sep. 2016.

[25] L. Yang, J. Fang, H. Li and B. Zeng, “An iterative reweighted method for Tucker decomposition of incomplete tensors,” IEEE Transactions on Signal Processing, vol.64, no.18, pp.4817-4829, Sep. 2016.

[26] L. Yang, Z.H. Huang and X. Shi, “A fixed point iterative method for low-

M

n-rank tensor pursuit,” IEEE Transactions on Signal Processing, vol.61, no.11, pp.2952-2962, Jun. 2013.

ED

[27] Y. Ming and C. H. Zhang, “On tensor completion via nuclear norm minimization,” Foundations of Computational Mathematics, vol.16, no.4,

PT

pp.1031-1068, Aug. 2016.

[28] Y. Yang, Y. Feng and J. A. K. Suykens, “A rank-one tensor updating algorithm for tensor completion,” IEEE Signal Processing Letters, vol.22,

CE

no.10, pp.1633-1637, Oct. 2015.

AC

[29] M. Kilmer, K. Braman, N. Hao and R. Hoover, “Third-order tensors as operators on matrices: A theoretical and computational framework with applications in imaging,” SIAM Journal on Matrix Analysis and Applications, vol.34, no.1, pp.148-172, Jan. 2013.

[30] Z. Zhang and S. Aeron, “Exact tensor completion using t-SVD,” IEEE Transactions on Signal Processing, vol.65, no.6, pp.1511-1526, Dec. 2016.

35

ACCEPTED MANUSCRIPT

[31] M. Yuan and Y. Lin, “Model selection and estimation in regression with grouped variables,” Journal of the Royal Statistical Society, vol.68, no.1, pp.49-67, 2006.

CR IP T

[32] W. Deng, W. Yin and Y. Zhang, “Group sparse optimization by alternating direction method,” Proceedings of SPIE, Wavelets and Sparsity XV, vol.8858, 26 Sep. 2013.

[33] Z. Wang, M.J. Lai, Z. Lu, W. Fan, H. Davulcu and J.P. Ye, “Orthogonal

rank-one matrix pursuit for low rank matrix completion,” SIAM Journal

AN US

on Scientific Computing, vol.37, no.1, 2014.

[34] X. Liao, H. Li and L. Carin, “Generalized alternating projection for weighted-`2,1 minimization with applications to model-based compressive sensing,” SIAM Journal on Imaging Sciences, vol.7, no. 2, pp.797-823, 2014. [35] C.D. Martin, R. Shafer, and B. LaRue, “An order-p tensor factorization

M

with applications in imaging,” SIAM Journal on Scientific Computing, vol.35, no.1, pp.474-490, 2013.

ED

[36] M. Jaggi and M. Sulovsk, “A simple algorithm for nuclear norm regularized problems,” International Conference on Machine Learning (ICML), pp.471478, Jul. 2010, Haifa, Israel.

PT

[37] X. Zhang, Y. Yu and D. Schuurmans, “Accelerated training for matrixnorm regularization: A boosting approach,” Proceedings of International

CE

Conference on Neural Information Processing Systems 25 (NIPS 2012), pp.2915-2923, 3-6 Dec. 2012, Lake Tahoe, Nevada, USA.

AC

[38] X.Y. Liu, S. Aeron, V. Aggarwal, X. Wang and M.Y. Wu, “Fine-grained indoor localization with adaptively sampled RF fingerprints,” Hungarian Journal of Industry and Chemistry, vol.58, no.3, 2015.

[39] Y. Pang, X. Li and Y. Yuan, “Robust tensor analysis with L1-norm,” IEEE Transactions on Circuits and Systems for Video Technology, vol.20, no.2, pp.172-178, Feb. 2010. 36

ACCEPTED MANUSCRIPT

[40] Y. Yang, S. Mehrkanoon and J.A.K. Suykens, “Higher order matching pursuit for low rank tensor learning,” arXiv:1503.02216 [stat.ML], Mar. 2015. [41] Z. Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli, “Image quality as-

CR IP T

sessment: From error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600-612, Apr. 2004.

[42] D. Martin, C. Fowlkes, D. Tal, and J. Malik, “A database of human seg-

mented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,” in Proceedings of IEEE In-

AN US

ternational Conference of Computer Vision (ICCV), vol. 2, no. 11, pp. 416-423, 2001, Barcelona, Spain.

[43] W. Sun, Y. Chen, L. Huang and H.C. So, “Tensor completion via generalized tensor tubal rank minimization using general unfolding,” IEEE Signal Processing Letters, Vol.25, No.6, pp.868-872, Jun. 2018.

M

[44] W. Sun, X. Lin, H.C. So and L. Huang, “Iteratively reweighted tensor SVD for robust multi-dimensional harmonic retrieval,” Proceedings of the Inter-

ED

national Conference on Acoustics, Speech, and Signal Processing (ICASSP 2016), pp. 4312-3512, Mar. 2016, Shanghai, China. [45] X. Lin, L. Huang and W. Sun, “Lp-PARAFAC for joint DoD and DoA

PT

estimation in bistatic MIMO radar,” Proceedings of 2016 CIE International Conference on Radar, pp.1-5, Oct. 2016, Guangzhou, China.

AC

CE

[46] http://trace.eas.asu.edu/yuv/.

37