Open AccessArticle

Alternating Direction Method of Multipliers for Generalized Low-Rank Tensor Recovery

Jiarong Shi

Qingyan Yin

Xiuyun Zheng

and

Wei Yang

School of Science, Xi’an University of Architecture and Technology, Xi’an 710055, China

Author to whom correspondence should be addressed.

Algorithms 2016, 9(2), 28; https://doi.org/10.3390/a9020028

Submission received: 23 November 2015 / Revised: 26 March 2016 / Accepted: 13 April 2016 / Published: 19 April 2016

(This article belongs to the Special Issue Manifold Learning and Dimensionality Reduction)

Download

Browse Figures

Figure 1
Comparison of relative errors between Multi-linear Robust Principal Component Analysis (MRPCA) and Generalized Low-Rank Tensor Recovery (GLRTR) with varying λ. "> Figure 2
Relative errors of Generalized Low-Rank Tensor Recovery (GLRTR) with varying τ. "> Figure 3
Relative errors for different a and b. "> Figure 4
Relative errors for different r and p. "> Figure 5
Relative errors versus different Inverse Signal-to-Noise Ratio (ISNR) on sparse noise. "> Figure 6
Relative errors versus different Inverse Signal-to-Noise Ratio (ISNR) on Gaussian noise. "> Figure 7
Background modeling from lobby video. (a) draws the images with missing entries, (b) plots the recovered background, i.e., the low-rank images, (c) shows the recovered foreground, i.e., the sparse noised images, and (d) displays the recovered images. "> Figure 8
Background modeling from bootstrap video. (a) draws the images with missing entries, (b) plots the recovered background, i.e., the low-rank images, (c) shows the recovered foreground, i.e., the sparse noised images, and (d) diplays the recovered images. ">

Versions Notes

Abstract

Low-Rank Tensor Recovery (LRTR), the higher order generalization of Low-Rank Matrix Recovery (LRMR), is especially suitable for analyzing multi-linear data with gross corruptions, outliers and missing values, and it attracts broad attention in the fields of computer vision, machine learning and data mining. This paper considers a generalized model of LRTR and attempts to recover simultaneously the low-rank, the sparse, and the small disturbance components from partial entries of a given data tensor. Specifically, we first describe generalized LRTR as a tensor nuclear norm optimization problem that minimizes a weighted combination of the tensor nuclear norm, the l₁-norm and the Frobenius norm under linear constraints. Then, the technique of Alternating Direction Method of Multipliers (ADMM) is employed to solve the proposed minimization problem. Next, we discuss the weak convergence of the proposed iterative algorithm. Finally, experimental results on synthetic and real-world datasets validate the efficiency and effectiveness of the proposed method.

Keywords:

low-rank tensor recovery; low-rank matrix recovery; nuclear norm minimization; alternating direction method of multipliers

1. Introduction

In the past decade, the low-rank property of some datasets has been explored skillfully to recover both the low-rank and the sparse components or complete the missing entries. The datasets to be analyzed are usually modeled by matrices and the corresponding recovery technique is named as Low-Rank Matrix Recovery (LRMR). LRMR has received a significant amount of attention in some fields of information science such as computer vision, machine learning, pattern recognition, data mining and linear system identification. There are several appealing types of LRMR including Matrix Completion (MC) [1], Robust Principal Component Analysis (RPCA) [2] and Low-Rank Representation (LRR) [3]. Mathematically, we customarily formulate LRMR as matrix nuclear norm minimization problems that can be effectively solved by several scalable methods. Diverse variants of LRMR derive from the aforementioned three models, for instance, MC with noise (or stable MC) [4], stable RPCA [5,6], incomplete RPCA [7], LRR with missing entries [8], and LRMR based on matrix factorization or tri-factorization [9,10].

The fields of image and signal processing usually require processing large amounts of multi-way data, such as video sequences, functional magnetic resonance imaging sequences and direct-sequence code-division multiple access (DS-CDMA) systems [11]. For the datasets with multi-linear structure, traditional matrix-based data analysis is prone to destroy the spatial and temporal structure, and can be affected by the curse of dimensionality. In contrast, tensor-based data representation can avoid or alleviate the above deficiencies to some extent. Furthermore, tensor decompositions can be used to obtain a low-rank approximation of an investigated data tensor. Two of the most popular tensor decompositions are the Tucker model and the PARAFAC (Parallel Factor Analysis) model, and they can be thought of as the higher order generalizations of Singular Value Decomposition (SVD) and Principal Component Analysis (PCA) [12].

Generally speaking, the traditional tensor decomposition models are very sensitive to missing entries, arbitrary outliers and gross corruptions. On the contrary, Low-Rank Tensor Recovery (LRTR) is very effective to recover simultaneously the low-rank component and the sparse noise, or complete the missing entries. As the higher order generalization of LRMR, LRTR is composed mainly of tensor completion [13] and Multi-linear RPCA (MRPCA) [14]. The task of tensor completion is to recover missing or unsampled entries according to the low-rank structure. The alternating least squares method based on the PARAFAC decomposition was originally proposed to complete missing entries [15]. Subsequently, the Tucker decomposition was applied to the problem of tensor completion [16,17,18]. The multi-linear rank is commonly used to describe the low-rank property, although it is hard to be properly estimated. The tensor nuclear norm generalizes the matrix nuclear norm and becomes a new standard for measuring the low-rank structure of a tensor. Hence, LRTR can usually be boiled down to a tensor nuclear norm minimization problem.

Liu et al. [13] established a tensor nuclear norm minimization model for tensor completion and proposed the Alternating Direction Method of Multipliers (ADMM) for efficiently solving this model. Gandy et al. [19] considered the low-n-rank tensor recovery problem with some linear constraints and developed the Douglas–Rachford splitting technique, a dual variant of the ADMM. To reduce the computational complexity of the tensor nuclear norm minimization problems, Liu et al. [20] introduced the matrix factorization idea into the minimization model and obtained a much smaller matrix nuclear norm minimization problem. Tan et al. [21] transformed tensor completion into a linear least squares problem and proposed a nonlinear Gauss–Seidel method to solve the corresponding optimization problem. In [14], Shi et al. extended RPCA to the case of tensors and presented the MRPCA model. MRPCA, also named robust low-rank tensor recovery [22], was described as a tensor nuclear norm minimization which can be solved efficiently by the ADMM.

This paper studies a generalized model of LRTR. In this model, the investigated data tensor is assumed to be the superposition of a low-rank component, a gross sparse tensor and a small dense error tensor. The Generalized Low-Rank Tensor Recovery (GLRTR) aims mainly to recover the low-rank and the sparse components from partially observed entries. For this purpose, we establish a tensor nuclear norm minimization model for GLRTR. In addition, the ADMM method is adopted to solve the proposed convex model.

The rest of this paper is organized as follows. Section 2 provides mathematical notations and introduces preliminaries on tensor algebra. In Section 3, we review related works on LRTR. We present a generalized LRTR model and develop an efficient iterative algorithm for solving the proposed model in Section 4. Section 5 discusses the weak convergence result on the proposed algorithm. We report the experimental results in Section 6. Finally, Section 7 draws the conclusions.

2. Notations and Preliminaries

This section will briefly introduce mathematical notations and review tensor algebra. We denote tensors by boldface Euclid Math One, e.g.,

A

, matrices by boldface capital letters, e.g., A, vectors by boldface letters, e.g., a, and scalars by italic letters, e.g., a or A. A tensor with order N indicates that its entries can be expressed via N indices. For an Nth-order real tensor

A \in ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}}

, its

(i_{1}, i_{2}, \dots, i_{N})

entry is denoted by

a_{i_{1} i_{2} \dots i_{N}}

The n-mode matricization (unfolding) of tensor

A

, denoted by

A_{(n)}

, is an

I_{n} \times \prod_{j \neq n} I_{j}

matrix obtained by rearranging n-mode fibers to be the columns of the resulting matrix. Furthermore, we define an opposite operation “fold” as

{f o l d}_{n} (A_{(n)}) = A

. Given another tensor

ℬ \in ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}}

, the inner product between

A

and

ℬ

is defined as

〈 A, ℬ 〉 = \sum_{i_{1} = 1}^{I_{1}} \sum_{i_{2} = 1}^{I_{2}} \dots \sum_{i_{N} = 1}^{I_{N}} a_{i_{1} i_{2} \dots i_{N}} b_{i_{1} i_{2} \dots i_{N}}

. Then, the Frobenius norm (or l₂-norm) of tensor

A

can be induced by the above inner product, that is,

{‖ A ‖}_{F} = \sqrt{〈 A, A 〉}

. Moreover, we also give other norm definitions of a tensor, which will be needed in the following sections. The l₁-norm of tensor

A

is expressed as

{‖ A ‖}_{1} = \sum_{i_{1} = 1}^{I_{1}} \sum_{i_{2} = 1}^{I_{2}} \dots \sum_{i_{N} = 1}^{I_{N}} | a_{i_{1} i_{2} \dots i_{N}} |

and the tensor nuclear norm is

{‖ A ‖}_{*} = \sum_{n = 1}^{N} w_{n} {‖ A_{(n)} ‖}_{*}

, where

{‖ A_{(n)} ‖}_{*}

, named the nuclear norm of

A_{(n)}

, is the sum of singular values of

A_{(n)}, w_{n} \geq 0

and

\sum_{n = 1}^{N} w_{n} = 1

. Due to the fact that the nuclear norm of one matrix is the tightest convex relaxation of the rank function on the unit ball in the spectral norm, it is easy to verify that the tensor nuclear norm

{‖ A ‖}_{*}

is a convex function with respect to

A

The n-mode product of tensor

A

by a matrix

U \in ℝ^{J_{n} \times I_{n}}

, denoted as

A \times_{n} U

, is an Nth-order tensor with the dimensionality of

I_{1} \times \dots \times I_{n - 1} \times J_{n} \times I_{n - 1} \times \dots \times I_{N}

, and its

(i_{1}, \dots, i_{n - 1}, j, i_{n + 1} \dots, i_{N}) - t h

entry is calculated as

\sum_{i_{n} = 1}^{I_{n}} a_{i_{1} i_{2} \dots i_{N}} u_{j i_{n}}

. The multi-linear rank of tensor

A

is stipulated as a vector

(r a n k (A_{(1)}), r a n k (A_{(2)}), \dots, r a n k (A_{(N)}))

, where

r a n k (\cdot)

indicates the rank function of a matrix. Tensor

A

is low-rank if and only if

r a n k (A_{(n)}) ≪ I_{n}

for some n. Then, the Tucker decomposition of

A

is defined as

A \approx C \times_{1} U^{(1)} \times_{2} U^{(2)} \dots \times_{N} U^{(N)} ≜ C \times_{n = 1}^{N} U^{(n)},

(1)

where the core tensor

C \in ℝ^{J_{1} \times J_{2} \times \dots \times J_{N}}

, mode matrices

U^{(n)} \in ℝ^{J_{n} \times I_{n}}

and

J_{n} \leq I_{n}

. The multi-linear rank of a tensor has superiority over other definitions (e.g., the rank used in PARAFAC decomposition) because it is easy to compute. Refer to the survey [12] for further understanding of tensor algebra.

Within the field of low-rank matrix/tensor recovery, two proximal minimization problems are extensively employed:

min_{X} λ {‖ X ‖}_{1} + \frac{1}{2} {‖ X - A ‖}_{F}^{2},

(2)

min_{X} λ {‖ X ‖}_{*} + \frac{1}{2} {‖ X - A ‖}_{F}^{2},

(3)

where

A

and A are a given data tensor and matrix respectively,

λ \geq 0

is a regularization factor. To address the above two optimization problems, we first define an absolute thresholding operator

S_{λ} (\cdot) : ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}} \to ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}}

as below:

{(S_{λ} (A))}_{i_{1} i_{2} \dots i_{N}} = max (| a_{i_{1} i_{2} \dots i_{N}} | - λ, 0)

. It has been proven that problems (2) and (3) have closed-form solutions denoted by

S_{λ} (A)

[2] and

T_{λ} (A)

[23] respectively, where

T_{λ} (A) = U S_{λ} (Σ) V^{T}

and

U Σ V^{T}

is the singular value decomposition of A.

3. Related Works

Tensor completion and MRPCA are two important and appealing applications of LRTR. In this section, we review the related works on the aforementioned two applications.

As the higher order generalization of matrix completion, tensor completion aims to recover all missing entries with the aid of the low-rank (or approximately low-rank) structure of a data tensor. Although low-rank tensor decompositions are practical in dealing with missing values [15,16,17,18], we have to estimate properly the rank of an incomplete tensor in advance. In the past few years, the matrix nuclear norm minimization model has been extended to the tensor case [13]. Given an incomplete data tensor

D \in ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}}

and a corresponding sampling index set

Ω \subset [I_{1}] \times [I_{2}] \times \dots \times [I_{N}]

, we define a linear operator

P_{Ω} (\cdot) : ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}} \to ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}}

as follows: if

(i_{1}, i_{2}, \dots, i_{N}) \in Ω

{(P_{Ω} (D))}_{i_{1} i_{2} \dots i_{N}} = d_{i_{1} i_{2} \dots i_{N}}

; otherwise,

{(P_{Ω} (D))}_{i_{1} i_{2} \dots i_{N}} = 0

, where

[I_{n}] = {1, 2, \dots, I_{n}}

. Then, the original tensor nuclear norm minimization model for tensor completion [13] is described as:

min_{X} {‖ X ‖}_{*}, s . t . P_{Ω} (X) = P_{Ω} (D) .

(4)

To tackle problem (4), Liu et al. [13] developed three different algorithms and demonstrated experimentally that the ADMM is the most efficient algorithm in obtaining a high accuracy solution. If we further take the dense Gaussian noise into consideration, then the stable version of tensor completion is given as below [19]:

min_{X} {‖ X ‖}_{*}, s . t . {‖ P_{Ω} (X) - P_{Ω} (D) ‖}_{F} \leq η,

(5)

or its corresponding unconstrained formulation:

min_{X} λ {‖ X ‖}_{*} + \frac{1}{2} {‖ P_{Ω} (X) - P_{Ω} (D) ‖}_{F}^{2},

(6)

where the regularized parameter

λ (\geq 0)

is used to balance the low-rankness and the approximate error, and η is a known estimate of the noise level.

In RPCA, a data matrix is decomposed into the sum of a low-rank component and a sparse component, and it is possible to recover simultaneously the two components by principal component pursuit under some suitable assumptions [2]. Shi et al. [14] extended RPCA to the case of tensors and presented the framework of MRPCA, which regards the data tensor

D

as the sum of a low-rank tensor

A

and a sparse noise term

ℰ

. Mathematically, the low-rank and the sparse components can be simultaneously recovered by solving the following tensor nuclear norm minimization problem:

min_{A, ℰ} {‖ A ‖}_{*} + λ {‖ ℰ ‖}_{1}, s . t . D = A + ℰ .

(7)

In [22], MRPCA is also called as robust low-rank tensor recovery.

4. Generalized Low-Rank Tensor Recovery

Both tensor completion and MRPCA do not consider dense Gaussian noise corruptions. In this section, we investigate the model of Generalized Low-Rank Tensor Recovery (GLRTR) and develop a corresponding iterative scheme.

4.1. Model of GLRTR

The datasets contaminated by Gaussian noise are very universal in practical engineering applications. In view of this, we assume the data tensor

D

to be the superposition of the low-rank component

A

, the large sparse corruption

ℰ

and the Gaussian noise

G

. We also consider the case that some entries of

D

are missing. To recover simultaneously the above three terms, we establish a convex GLRTR model as follows:

min_{A, ℰ, G} {‖ A ‖}_{*} + λ {‖ ℰ ‖}_{1} + τ {‖ G ‖}_{F}^{2}, s . t . P_{Ω} (D) = P_{Ω} (A + ℰ + G),

(8)

where the regularization coefficients λ and τ are nonnegative.

If we reinforce the constraints

G = O

(or equivalently

τ \to + \infty

) and

ℰ = O

(or equivalently

λ \to + \infty

), then the GLRTR is transformed into the tensor completion model (4), where the zero tensor

O

has the same dimensionality as

D

. If only the constraint

ℰ = O

is considered, then the GLRTR is equivalent to Equation (6). Furthermore, if we take

G = O

and

Ω = [I_{1}] \times [I_{2}] \times \dots \times [I_{N}]

, then the model of GLRTR becomes the model of MRPCA. In summary, the proposed model is the generalization of the existing LRTR.

For the convenience of using the splitting method, we discard

A

and introduce N + 1 auxiliary tensor variables

ℳ_{1}, ℳ_{2}, \dots, ℳ_{N}, X

, where these auxiliary variables have the same dimensionality as

D

. Let

M_{n (n)}

be the n-mode matricization of

ℳ_{n}

for each

n \in [N]

and

ℳ = {ℳ_{1}, ℳ_{2}, \dots, ℳ_{N}}

. Hence, we have the equivalent formulation of Equation (8):

\begin{array}{l} min_{X, ℰ, G, M} \sum_{N = 1}^{n} w_{n} {‖ M_{n (n)} ‖}_{*} + λ {‖ ℰ ‖}_{1} + τ {‖ G ‖}_{F}^{2}, \\ s . t . X = M_{n} + G + ℰ, P_{Ω} (X) = P_{Ω} (D), n = 1, 2, \dots, N . \end{array}

(9)

As a matter of fact, the discarded low-rank tensor

A

can be represented by

A = \sum_{N = 1}^{n} M_{n} / N

. Now, we explain the low-rankness of

A

from two viewpoints. One is

{‖ A ‖}_{*} \leq \sum_{N = 1}^{n} {‖ M_{n} ‖}_{*} / N \leq max_{1 \leq n \leq N} {‖ M_{n} ‖}_{*}

. The other is that a better solution to Equation (9) will satisfy

ℳ_{1} \approx ℳ_{2} \approx \dots \approx ℳ_{N}

. These viewpoints illustrate that

A

is approximately low-rank along each mode. The aforementioned non-smooth minimization problem is distributed convex. Concretely speaking, the tensor variables can be split into several parts and the objective function is separable across this splitting.

4.2. Optimization Algorithm to GLRTR

As a special splitting method, the ADMM is very efficient to solve a distributed optimization problem with linear equality constraints. It takes the form of the decomposition–coordination procedure and blends the merits of dual decomposition and augmented Lagrangian methods. In this section, we will propose the method of ADMM to solve Equation (9).

We first construct the augmented Lagrangian function of the aforementioned convex optimization problem without considering the constraint

P_{Ω} (X) = P_{Ω} (D)

\begin{array}{l} L_{μ} (X, & ℳ, ℰ, G, Y) = \sum_{n = 1}^{N} w_{n} {‖ M_{n (n)} ‖}_{*} + λ {‖ ℰ ‖}_{1} + τ {‖ G ‖}_{F}^{2} + \\ \sum_{n = 1}^{N} (〈 Y_{n}, X - ℳ_{n} - G - ℰ 〉 + \frac{μ}{2} {‖ X - ℳ_{n} - G - ℰ ‖}_{F}^{2}), \end{array}

(10)

where μ is a positive numerical constant,

Y_{n}

are Lagrangian multiplier tensors and

Y = {Y_{1}, Y_{2}, \dots, Y_{N}}

. There are totally five blocks of variables in

L_{μ} (X, ℳ, ℰ, G, Y)

. The ADMM algorithm updates alternatively each block of variables while letting the other blocks of variables be fixed. If

Y

is given, one block of variables in

{X, ℳ, G, ℰ}

can be updated by minimizing

L_{μ} (X, ℳ, ℰ, G, Y)

with respect to one argument. On the contrary, if

{X, ℳ, G, ℰ}

are given,

Y

is updated by maximizing

L_{μ} (X, ℳ, ℰ, G, Y)

with respect to

Y

. The detailed iterative procedure is outlined as follows:

Computing

X

. Fix

{ℳ, ℰ, G, Y}

and minimize

L_{μ} (X, ℳ, ℰ, G, Y)

with respect to

X

\begin{array}{l} X : & = arg min_{X} L_{μ} (X, ℳ, ℰ, G, Y) \\ = arg min_{X} \sum_{n = 1}^{N} {‖ X - ℳ_{n} - G - ℰ + Y_{n} / μ ‖}_{F}^{2} \\ = \tilde{X}, \end{array}

(11)

where

\tilde{X} = \sum_{n = 1}^{N} (ℳ_{n} - Y_{n} / μ) / N + G + ℰ

. In consideration of the constraint

P_{Ω} (X) = P_{Ω} (D)

, the final update formulation of

X

is revised as

X : = P_{\bar{Ω}} (\tilde{X}) + P_{Ω} (D),

(12)

where

\bar{Ω}

is the complement set of Ω.

Computing

ℳ

. If

ℳ_{n}

is unknown and other variables are fixed, we update

ℳ_{n}

by minimizing

L_{μ} (X, ℳ, ℰ, G, Y)

with respect to

ℳ_{n}

\begin{array}{l} ℳ_{n} : & = arg min_{ℳ_{n}} L_{μ} (X, ℳ, ℰ, G, Y) \\ = arg min_{ℳ_{n}} w_{n} {‖ M_{n (n)} ‖}_{*} + \frac{μ}{2} {‖ (X - G - ℰ + Y_{n} / μ) - ℳ_{n} ‖}_{F}^{2} \\ = arg min_{ℳ_{n}} \frac{w_{n}}{μ} {‖ M_{n (n)} ‖}_{*} + \frac{1}{2} {‖ (Q_{(n)} - M_{n (n)}) - ℳ_{n} ‖}_{F}^{2} \\ = {fold}_{n} (T_{w_{n} / μ} (Q_{(n)})), \end{array}

(13)

where

Q = X - G - ℰ + Y_{n} / μ

and

Q_{(n)}

is the n-mode matricization of

Q

Computing

ℰ

. If

ℰ

is unknown and other variables are fixed, the calculation procedure of

ℰ

is given as follows:

\begin{array}{l} ℰ : & = arg min_{ℰ} L_{μ} (X, ℳ, ℰ, G, Y) \\ = arg min_{ℰ} \frac{λ}{μ N} {‖ ℰ ‖}_{1} + \frac{1}{2} {‖ ℰ - \sum_{n = 1}^{N} (X - ℳ_{n} - G + Y_{n} / μ) / N ‖}_{F}^{2} \\ = S_{(λ / μ N)} (\sum_{n = 1}^{N} (X - ℳ_{n} - G + Y_{n} / μ) / N) . \end{array}

(14)

Computing

G

. The update formulation of

G

is calculated as

\begin{array}{l} G : & = arg min_{G} L_{μ} (X, ℳ, ℰ, G, Y) \\ = arg min_{G} τ {‖ G ‖}_{F}^{2} + \frac{μ N}{2} {‖ G - \sum_{n = 1}^{N} (X - ℳ_{n} - ℰ + Y_{n} / μ) / N ‖}_{F}^{2} \\ = \frac{μ}{2 τ + μ N} \sum_{n = 1}^{N} (X - ℳ_{n} - ℰ + Y_{n} / μ) . \end{array}

(15)

Computing

Y

. If

{X, ℳ, ℰ, G}

are fixed, we can update the Lagrangian multipliers

Y

as follows:

Y_{n} : = Y_{n} + μ (X - ℳ_{n} - G - ℰ), n = 1, 2, \dots, N .

(16)

Once

{X, ℳ, ℰ, G, Y}

are updated, we will increase the value of μ by multiplying it with a constant

ρ (> 1)

during the procedure of iterations. The whole iterative procedure is outlined in Algorithm 1. The stopping condition of Algorithm 1 is set as

max_{1 \leq n \leq N} {‖ X - ℳ_{n} - G - ℰ ‖}_{F}^{2} < ε

or the maximum number of iterations is reached, where ε is a sufficiently small positive number.

Algorithm 1. Solving GLRTR by ADMM

Input: Data tensor

D

, sampling index set Ω, regularization parameters λ and τ.

Initialize:

ℳ, ℰ, G, Y, μ, μ_{max}

and ρ.

Output:

X, ℰ, G

and

ℳ

While not converged do

Update $X$ according to (12).
Update $ℳ_{n}$ according to (13), $n = 1, 2, \dots, N$ .
Update $ℰ$ according to (14).
Update $G$ according to (15).
Update $Y_{n}$ according to (16), $n = 1, 2, \dots, N$ .
Update μ as $μ : = min (ρ μ, μ_{max})$ .

End while

In our implementation,

ℳ, ℰ, G

and

Y

are initialized to zeros tensors. The other parameters are set as follows:

w_{1} = w_{2} = \dots = w_{N} = 1 / N, μ = 10^{- 4}, ρ = 1.1, μ_{max} = 10^{10}

and

ε = 10^{- 8} .

Furthermore, the maximum number of iterations is set to 100.

5. Convergence Analysis

Because the number of block variables in Equation (9) is more than two, it is difficult for us to prove the convergence of Algorithm 1. Nevertheless, the experimental results in the next section demonstrate this algorithm has good convergence behavior. This section will discuss the weak convergence result on our ADMM algorithm.

We consider a special case of Algorithm 1, that is, there is no missing entries. In this case, it holds that

X = D

. Thus, Equation (9) is transformed into

\begin{array}{l} min_{ℳ, ℰ, G} w_{n} {‖ M_{n (n)} ‖}_{*} + λ {‖ ℰ ‖}_{1} + τ {‖ G ‖}_{F}^{2} \\ s . t . D = ℳ_{n} + G + ℰ, n = 1, 2, \dots, N . \end{array}

(17)

We can design a corresponding ADMM algorithm by revising Algorithm 1, namely, the update formulation of

X

is replaced by

X = D

. In fact, the revised algorithm is an inexact version of ADMM. Subsequently, we give the iterative formulations of exact ADMM for Equation (17):

{\begin{cases} ℳ : = arg min_{ℳ} L_{μ} (D, ℳ, ℰ, G, Y), \\ (ℰ, G) : = arg min_{ℰ, G} L_{μ} (D, ℳ, ℰ, G, Y), \\ Y_{n} : = Y_{n} + μ (D - ℳ_{n} - G - ℰ), n = 1, 2, \dots, N . \end{cases}

(18)

In this circumstance, we have the following statement:

Theorem 1.

Let

{ℳ^{(k)}, ℰ^{(k)}, G^{(k)}, Y^{(k)}}

be the sequence produced by (18). If

L_{0} (D, ℳ, ℰ, G, Y)

has a saddle point with respect to

{ℳ, ℰ, G, Y}

, then

{ℳ^{(k)}, ℰ^{(k)}, G^{(k)}}

satisfy that the iterative sequence of objective function in Equation (17) converges to the minimum value.

Proof.

By matricizing each tensor along the one-mode, the constraints of Equation (17) can be rewritten as a linear equation:

(\begin{array}{l} D_{(1)} \\ D_{(1)} \\ ⋮ \\ D_{(1)} \end{array}) = (\begin{array}{l} M_{1 (1)} \\ M_{2 (1)} \\ ⋮ \\ M_{N (1)} \end{array}) + (\begin{array}{l} G_{(1)} \\ G_{(1)} \\ ⋮ \\ G_{(1)} \end{array}) + (\begin{array}{l} E_{(1)} \\ E_{(1)} \\ ⋮ \\ E_{(1)} \end{array}),

(19)

or equivalently,

(\begin{array}{l} I \\ I \\ ⋮ \\ I \end{array}) D_{(1)} = (\begin{array}{l} M_{1 (1)} \\ M_{2 (1)} \\ ⋮ \\ M_{N (1)} \end{array}) + (\begin{array}{l} I \\ I \\ ⋮ \\ I \end{array}) G_{(1)} + (\begin{array}{l} I \\ I \\ ⋮ \\ I \end{array}) E_{(1)},

(20)

where I is an identity matrix of size I₁ × I₁. If we integrate

G_{(1)}

and

E_{(1)}

into one block of variables, then Equation (20) is re-expressed as below:

(\begin{array}{l} I \\ I \\ ⋮ \\ I \end{array}) D_{(1)} = (\begin{array}{l} M_{1 (1)} \\ M_{2 (1)} \\ ⋮ \\ M_{N (1)} \end{array}) + (\begin{matrix} \begin{array}{l} I \\ I \\ ⋮ \\ I \end{array} & \begin{array}{l} I \\ I \\ ⋮ \\ I \end{array} \end{matrix}) (\begin{array}{l} G_{(1)} \\ E_{(1)} \end{array}) .

(21)

We denote

M = {M_{1 (1)}, M_{2 (1)}, \dots, M_{N (1)}}

and partition all variables in Equation (17) into two parts in form of matrices: M and

(E_{(1)}, G_{(1)})

. Essentially,

M, E_{(1)}

and

G_{(1)}

are the matricizations of

ℳ, ℰ

and

G

, respectively. Let

g (M) = \sum_{n = 1}^{N} w_{n} {‖ M_{n (n)} ‖}_{*}, h (E_{(1)}, G_{(1)}) = λ {‖ E_{(1)} ‖}_{1} + τ {‖ G_{(1)} ‖}_{F}^{2}

. Then, the objective function in Equation (17) can be expressed as

g (M) + h (E_{(1)}, G_{(1)})

. It is obvious that

g (M)

and

h (E_{(1)}, G_{(1)})

are two closed, proper and convex functions. According to the basic convergence result given in [24], we have

lim_{k \to + \infty} g (M^{(k)}) + h (E_{(1)}^{(k)}, G_{(1)}^{(k)}) = f^{*}

, where

f^{*}

is the minimum value of

g (M) + h (E_{(1)}, G_{(1)})

under the linear Constraint (21). This ends the proof. □

In the following, we discuss the detailed iterative procedure of

{M, E_{(1)}, G_{(1)}}

{ℳ, ℰ, G}

in Equation (18). It is easy to verify that the update formulation of

ℳ

in Equation (18) is equivalent to:

\begin{array}{l} (M_{1 (1)}, M_{2 (1)}, \dots, M_{N (1)}) : = arg min_{M} L_{μ} (D, ℳ, ℰ, G, Y) \\ = arg min_{(M_{1 (1)}, M_{2 (1)}, \dots, M_{N (1)})} \sum_{n = 1}^{N} (w_{n} {‖ M_{n (n)} ‖}_{*} + \frac{μ}{2} {‖ D - ℳ_{n} - G - ℰ + \frac{1}{μ} Y_{n} ‖}_{F}^{2}) . \end{array}

(22)

The block of variables

M_{n (1)}

can be solved in parallel because the objective function in Equation (22) is separable with respect to

M_{1 (1)}, M_{2 (1)}, \dots, M_{N (1)}

. Hence, the iterative formulation of

ℳ_{n}

is similar to that of Equation (13).

For fixed M and

Y

, we can get the optimal block of variables

(ℰ, G)

(E_{(1)}, G_{(1)})

by minimizing

f (ℰ, G)

, where

f (ℰ, G) = λ {‖ ℰ ‖}_{1} + τ {‖ G ‖}_{F}^{2} + \frac{μ}{2} \sum_{n = 1}^{N} {‖ (D - ℳ_{n} + Y_{n} / μ) - (G + ℰ) ‖}_{F}^{2}

. The partial derivative of

f (ℰ, G)

with respect to

G

\frac{\partial f}{\partial G} = 2 τ G + μ \sum_{n = 1}^{N} ((G + ℰ) - (D - ℳ_{n} + Y_{n} / μ)) .

(23)

By letting

\frac{\partial f}{\partial G} = 0

, we have

G : = δ T - δ N ℰ,

(24)

where

δ = \frac{μ}{2 τ + μ N}, T = \sum_{n = 1}^{N} (D - ℳ_{n} + Y_{n} / μ)

. By substituting Equation (24) into

f (ℰ, G)

, we have the update formulation of

ℰ

\begin{array}{l} ℰ : = arg min_{ℰ} λ {‖ ℰ ‖}_{1} + τ {‖ δ T - δ N ℰ ‖}_{F}^{2} + \frac{μ}{2} \sum_{n = 1}^{N} {‖ (D - ℳ_{n} + Y_{n} / μ) - ℰ - δ T + δ N ℰ ‖}_{F}^{2} \\ = arg min_{ℰ} λ {‖ ℰ ‖}_{1} + τ δ^{2} N^{2} {‖ ℰ - \frac{1}{N} T ‖}_{F}^{2} + \frac{μ {(1 - δ N)}^{2} N}{2} {‖ ℰ - \frac{1}{(1 - δ N) N} \sum_{n = 1}^{N} (D - ℳ_{n} + Y_{n} / μ - δ T) ‖}_{F}^{2} \\ = arg min_{ℰ} λ {‖ ℰ ‖}_{1} + τ δ^{2} N^{2} {‖ \frac{1}{N} T - ℰ ‖}_{F}^{2} + \frac{μ {(1 - δ N)}^{2} N}{2} {‖ ℰ - \frac{1}{N} T ‖}_{F}^{2} \\ = arg min_{ℰ} λ {‖ ℰ ‖}_{1} + \frac{τ δ^{2} N^{2} + μ {(1 - δ N)}^{2} N}{2} {‖ ℰ - \frac{1}{N} T ‖}_{F}^{2} \\ = S_{(λ / (τ δ^{2} N^{2} + μ {(1 - δ N)}^{2} N))} (\frac{1}{N} T) . \end{array}

(25)

A standard extension of the classic ADMM is to use varying parameter μ for each iteration. The goal of this extension is to improve the convergence and reduce the dependence on the initial choice of μ. In the context of the method of multipliers, this approach is proven to be superlinearly convergent if

μ^{(k)} \to + \infty

[25], which inspires us to adopt the non-decreasing sequence

{μ^{(k)}}

. Furthermore, large values of μ result in a large penalty on violations of primal feasibility and are thus inclined to produce small primal residuals. However, it is difficult to prove the convergence of ADMM [24]. A commonly-used choice for

{μ^{(k)}}

μ^{(k + 1)} : = min (ρ μ^{(k)}, μ_{max})

, where

ρ > 1

and

μ_{max}

is the upper limit of μ. Due to the fact that

μ^{(k)}

is fixed after a finite number of iterations, the corresponding ADMM is convergent according to Theorem 1.

6. Experimental Results

We perform experiments on synthetic data and two real-world video sequences, and validate the feasibility and effectiveness of the proposed method. The experimental results of GLRTR are compared with that of MRPCA, where the missing values in MRCPA are replaced by zeros.

6.1. Synthetic Data

In this subsection, we synthesize data tensors with missing entries. First, we generate an Nth-order low-rank tensor as follows:

A = C \times_{n = 1}^{N} U^{(n)},

with the core tensor

C \in ℝ^{J_{1} \times J_{2} \times \dots \times J_{N}}

and mode matrices

U^{(n)} \in ℝ^{I_{n} \times J_{n}} (J_{n} < I_{n})

. The entries of

C

and

U^{(n)}

are independently drawn from the standard normal distribution. Then, we generate randomly a dense noise tensor

G \in ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}}

whose entries also obey the standard normal distribution. Next, we construct a sparse noise tensor

P_{Ω^{'}} (ℰ)

, where

ℰ \in ℝ^{I_{1} \times I_{2} \times \dots \times I_{N}}

is produced by a uniform distribution on the interval

(- a, a)

and the index set

Ω^{'}

is produced by uniformly sampling on

[I_{1}] \times [I_{2}] \times \dots \times [I_{N}]

with probability p′%. Finally, the generation of the sampling index set Ω is similar to Ω′ and the corresponding sampling rate is set to be p%. Therefore, an incomplete data tensor is synthesized as

D_{Ω} = P_{Ω} (A + G + ℰ)

For given data tensor

D_{Ω}

with missing values, its low-rank component recovered by some method is denoted by

\hat{A}

. The Relative Error (RE) is employed to evaluate the recovery performance of the low-rank structure and its definition is given as follows:

RE = {‖ \hat{A} - A ‖}_{F} / {‖ A ‖}_{F}

. Small relative error means good recovery performance. The experiments are carried out on 50 × 50 × 50 tensors and 20 × 20 × 20 × 20 tensors, respectively. Furthermore, we set a = 500 and

J_{1} = J_{2} = \dots = J_{N} = r

For convenience of comparison, we design three groups of experiments. In the first group of experiments, we only consider the case that there are no missing entries, that is,

Ω = [I_{1}] \times [I_{2}] \times \dots \times [I_{N}]

or p = 100. The values of

{‖ ℰ ‖}_{F} / {‖ A ‖}_{F}

and

{‖ G ‖}_{F} / {‖ A ‖}_{F}

are adopted to indicate the Inverse Signal-to-Noise Ratio (ISNR) with respect to the sparse and the Gaussian noise respectively. Three different degrees of sparsity are taken into account, that is, p′ = 5, 10 and 15. In addition, we take

λ = 0.025, τ = 0.02, r \in {3, 5}

for 3rd-order tensors and

λ = 0.03, τ = 0.01, r \in {2, 4}

for 4th-order tensors. For given parameters, we repeat the experiments ten times and report the average results. As a low-rank approximation method for tensors, the Higher-Order SVD (HOSVD) truncation [12] to rank-

(r_{1}, r_{2}, \dots, r_{N})

is not suitable for gross corruptions owing to the fact that its relative error reaches up to 97% or even 100%. Hence, we do not compare this method with our GLRTR in subsequent experiments. The experimental results are shown in Table 1 and Table 2, respectively.

From the above two tables, we have the following observations: (I) Although the values of ISNR on sparse noise are very large, both MRPCA and GLRTR remove efficiently sparse noise to some extent. Meanwhile, large values of ISNR on Gaussian noise are disadvantageous for recovering the low-rank components; (II) GLRTR has better recovery performance than MRPCA. In an average sense, the relative error of GLRTR is 2.68% smaller than that of MRPCA for 50 × 50 × 50 tensors, and 14.17% for 20 × 20 × 20 × 20 tensors; (III) For 3rd-order tensors, GLRTR removes effectively Gaussian noise, and, on average, its relative error is 5.96% smaller than the value of ISNR on Gaussian noise; for 4th-order tensors, GLRTR effectively removes Gaussian noise only in the case r = 2. In summary, GLRTR is more effective than MRPCA in recovering the low-rank components.

The second group of experiments considers four different sampling rates for Ω and one fixed degree of sparsity for Ω′, that is,

p \in {30, 50, 70, 90}

and p′ = 5. We set τ = 0.02 for both 3rd-order and 4th-order tensors, and choose the superior tradeoff parameter λ for each p. The comparisons of experimental results between MRPCA and GLRTR are shown in Table 3 and Table 4, respectively. We can see from these two tables that MRPCA is very sensitive to the sampling rate p%, and it hardly recovers the low-rank components. In contrast, GLRTR achieves better recovery performance for 3rd-order tensors. As for 4th-order tensors, it also has smaller relative error when the sampling rate p% is relatively large. These observations show that GLRTR is more robust to missing values than MRPCA.

We will evaluate the sensitivity of GLRTR to the choice of λ and τ in the last group of experiments. For convenience of designing experiments, we only perform experiments on 50 × 50 × 50 tensors and consider the case that p = 100. The values of λ and τ are set according to the following manner: we vary the value of one parameter while letting the other be fixed. In the first case, the parameter τ is chosen as 0.01. Under this circumstance, the relative errors versus different λ of MRPCA and GLRTR are shown in Figure 1. We take λ = 0.01 in the second case and the relative errors versus different τ of GLRTR are shown in Figure 2.

It can be seen from Figure 1 that the relative errors of MRPCA and GLRTR are about 9.90% and 4.62%, respectively, if 0.02 ≤ λ ≤ 0.07, which means the latter has better recovery performance than the former. Furthermore, their relative errors are relatively stable when λ lies within a certain interval. Figure 2 illustrates that the relative error has the tendency to increase monotonically, and it becomes almost stationary when τ ≥ 1. At this moment, the relative errors lie in the interval (0.037, 0.080). This group of experiments implies that, for our synthetic data, GLRTR is not very sensitive to the choice of λ and τ.

6.2. Influence of Noise and Sampling Rate on the Relative Error

This subsection will evaluate the influence of noise and sampling rate on the relative error. For this purpose, we design four groups of experiments and use the synthetic data generated in the same manner as in the previous subsection. For the sake of convenience, we only carry out experiments on 50 × 50 × 50 tensors.

The first group of experiments aims to investigate the influence of noise on the recovery performance. In the data generation process, we only change the manner for generating

G

, that is, each entry of

G

is drawn independently from the normal distribution with mean zero and standard deviation b. Let

r = 5, p^{'} = 10, p = 100, a = 50 i

and b = 0.2j, where

i, j \in {0, 1, \dots, 10}

. For different combinations of a and b, the relative errors of GLRTR are shown in Figure 3. We can draw two conclusions from this figure. For given b, the relative error is relatively stable with the increasing of a, which means the relative error is not very sensitive to the magnitude of sparse noise. The relative error monotonically increases with the increasing of Gaussian noise level, which validates that large Gaussian noise is disadvantageous for recovering the low-rank component.

Next, we study the influence of sampling rate p% on the relative error for different r. Set

p^{'} = 10, a = 500, b = 1, r = 1 + 2 i, i = 1, 2, \dots, 14

and vary the value of p% from 30 to 100 in steps of size 10. For fixed r and p, we obtain the low-rank component according to GLRTR and then plot the relative errors in Figure 4. From the 3-D colored surface in Figure 4, we can see that both r and p have significant influence on the relative error. This observation indicates that small r or large p is conducive to the recovery of low-rank term.

The third group of experiments will validate the robustness of GLRTR to sparse noise. Concretely speaking, we investigate the recovery performance under different ISNR on sparse noise without consideration of Gaussian noise. Set

r \in {3, 5}, p^{'} = 10, p = 100, b = 0, a = 2^{i} \times 500

i = - 15, - 14, \dots, 10

. The experimental results are shown in Figure 5, where the horizontal and the vertical coordinates represent the ISNR on sparse noise and the relative error, respectively. This figure illustrates that the relative error is less than 4.5% for synthetic 3rd-order tensors, which verifies experimentally that our method is very robust to sparse noise.

In the last group of experiments, we discuss the performance of GLRTR in removing the small dense noise for given large sparse noise. We also propose a combination strategy: GLRTR + HOSVD, that is, GLRTR is followed by HOSVD. The goal of this new method is to improve the denoising performance of GLRTR. Let

r \in {3, 5}, p^{'} = 10, p = 100, a = 500

and

b = 0.1 i, i = 0, 1, \dots, 50

. Different values for b lead to different ISNR on Gaussian noise. We draw four curves to reflect the relationship between the relative error and ISNR on Gaussian noise, as shown in Figure 6, where the black dashed line is a reference line. We have two observations from this figure. When ISNR on Gaussian noise is larger than 3.5%, GLRTR not only successfully separates the sparse noise to some extent but also effectively removes the Gaussian noise. The GLRTR+HOSVD method has better denoising performance than GLRTR in the presence of large Gaussian noise.

6.3. Applications in Background Modeling

In this subsection, we test our method on two real-world surveillance videos for object detection and background subtraction: Lobby and Bootstrap datasets [26]. For convenience of computation, we only consider the first 200 frames for each dataset and transform the color images into the gray-level images. The resolutions of each image in the Lobby and Bootstrap datasets are 128 × 160 and 120 × 160, respectively. We add Gaussian noise with mean zero and standard deviation 5 to each image. Hence, we obtain two data tensors of order 3 and their sizes are 128 × 160 × 200 and 120 × 160 × 200, respectively. For two given tensors, we execute random sampling on them with a probability of 50%.

Considering the fact that MRPCA fails in recovering the low-rank components on the synthetic data with missing values, we only implement the method of GLRTR on the video datasets. Two tradeoff parameters are set as follows: λ = 0.0072 and τ = 0.001. We can obtain the low-rank, the sparse and the completed components from the incomplete data tensors according to the proposed method. Actually, the low-rank terms are the backgrounds and the sparse noise terms correspond to the foregrounds. The experimental results are partially shown in Figure 7 and Figure 8, respectively, where the missing entries in incomplete images are shown in white. From Figure 7 and Figure 8, we can see that GLRTR can recover efficiently the low-rank images and the sparse noise images. Moreover, we observe from the recovered images that a large proportion of missing entries are effectively completed.

To evaluate the completion performance of GLRTR, we define the Relative Approximate Error (RAE) as

RAE = {‖ \hat{D} - D ‖}_{F} / {‖ D ‖}_{F}

, where

D

is the original video tensor without Gaussian noise corruptions and missing entries, and

\hat{D}

is the approximated term of

D

. The RAE of the Lobby dataset is 8.94% and that of the Bootstrap dataset is 20.11%. These results demonstrate GLRTR can complete approximately the missing entries to a certain degree. There are two reasons for that the Bootstrap dataset has relatively large RAE: one is its more complex foreground and the other is that the entries of the foreground can not be recovered when they are missing. In summary, GLRTR is robust to gross corruption, Gaussian noise and missing values.

7. Conclusions

In this paper, we investigate a generalized model of LRTR in which large sparse corruption, missing entries and Gaussian noise are taken into account. For the generalized LRTR, we establish an optimization problem that minimizes the weighted combination of the tensor nuclear norm, the l₁ norm and the Frobenius norm. To address this minimization problem, we present an iterative scheme based on the technique of ADMM. The experimental results on synthetic data and real-world video datasets illustrate that the proposed method is efficient and feasible in recovering the low-rank components and completing missing entries. In the future, we will consider the theoretical conditions for exact recoverability and other scalable algorithms.

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (No. 61403298, No. 11526161 and No. 11401357), and the Natural Science Basic Research Plan in Shaanxi Province of China (No. 2014JQ8323).

Author Contributions

Jiarong Shi constructed the model, developed the algorithm and wrote the manuscript. Qingyan Yin designed the experiments. Xiuyun Zheng and Wei Yang implemented all experiments. All four authors were involved in organizing and refining the manuscript. All authors have read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Candès, E.J.; Recht, B. Exact matrix completion via convex optimization. Found. Comput. Math. 2009, 9, 717–772. [Google Scholar] [CrossRef]
Candès, E.J.; Li, X.; Ma, Y.; Wright, J. Robust principal component analysis? J. ACM. 2011, 58, 37. [Google Scholar] [CrossRef]
Liu, G.; Lin, Z.; Yan, S.; Sun, J.; Yu, Y.; Ma, Y. Robust recovery of subspace structures by low-rank representation. IEEE Trans. Pattern Anal. Mach. Intell. 2013, 35, 171–184. [Google Scholar] [CrossRef] [PubMed]
Candès, E.J.; Plan, Y. Matrix completion with noise. P. IEEE. 2010, 98, 925–936. [Google Scholar] [CrossRef]
Zhou, Z.; Li, X.; Wright, J.; Candès, E.J.; Ma, Y. Stable principal component pursuit. In Proceedings of the 2010 IEEE International Symposium on Information Theory Proceedings (ISIT), Austin, TX, USA, 13–18 June 2010.
Xu, H.; Caramanis, C.; Sanghavi, S. Robust PCA via outlier pursuit. IEEE Trans. Inf. Theory. 2012, 58, 3047–3064. [Google Scholar] [CrossRef]
Shi, J.; Zheng, X.; Yong, L. Incomplete robust principal component analysis. ICIC Express Letters, Part B Appl. 2014, 5, 1531–1538. [Google Scholar]
Shi, J.; Yang, W.; Yong, L.; Zheng, X. Low-rank representation for incomplete data. Math. Probl. Eng. 2014, 10. [Google Scholar] [CrossRef]
Liu, Y.; Jiao, L.C.; Shang, F.; Yin, F.; Liu, F. An efficient matrix bi-factorization alternative optimization method for low-rank matrix recovery and completion. Neural Netw. 2013, 48, 8–18. [Google Scholar] [CrossRef] [PubMed]
Liu, Y.; Jiao, L.C.; Shang, F. A fast tri-factorization method for low-rank matrix recovery and completion. Pattern Recogn. 2013, 46, 163–173. [Google Scholar] [CrossRef]
De Lathauwer, L.; Castaing, J. Tensor-based techniques for the blind separation of DS–CDMA signals. Signal Process. 2007, 87, 322–336. [Google Scholar] [CrossRef]
Kolda, T.G.; Bader, B.W. Tensor decompositions and applications. SIAM Rev. 2009, 51, 455–500. [Google Scholar] [CrossRef]
Liu, J.; Musialski, P.; Wonka, P.; Ye, J. Tensor completion for estimating missing values in visual data. IEEE Trans. Patt. Anal. Mach. Intel. 2013, 35, 208–220. [Google Scholar] [CrossRef] [PubMed]
Shi, J.; Zhou, S.; Zheng, X. Multilinear robust principal component analysis. Acta Electronica Sinica. 2014, 42, 1480–1486. [Google Scholar]
Tomasi, G.; Bro, R. PARAFAC and missing values. Chemom. Intell. Lab. Syst. 2005, 75, 163–180. [Google Scholar] [CrossRef]
Shi, J.; Jiao, L.C.; Shang, F. Tensor completion algorithm and its applications in face recognition. Pattern Recognit. Artif. Intell. 2011, 24, 255–261. [Google Scholar]
Kressner, D.; Steinlechner, M.; Vandereycken, B. Low-rank tensor completion by Riemannian optimization. BIT Numer. Math. 2014, 54, 447–468. [Google Scholar] [CrossRef]
Shi, J.; Yang, W.; Yong, L.; Zheng, X. Low-rank tensor completion via Tucker decompositions. J. Comput. Inf. Syst. 2015, 11, 3759–3768. [Google Scholar]
Gandy, S.; Recht, B.; Yamada, I. Tensor completion and low-n-rank tensor recovery via convex optimization. Inv. Probl. 2011, 27, 19. [Google Scholar] [CrossRef]
Liu, Y.; Shang, F. An efficient matrix factorization method for tensor completion. IEEE Signal Process. Lett. 2013, 20, 307–310. [Google Scholar] [CrossRef]
Tan, H.; Cheng, B.; Wang, W.; Zhang, Y.J.; Ran, B. Tensor completion via a multi-linear low-n-rank factorization model. Neurocomputing. 2014, 133, 161–169. [Google Scholar] [CrossRef]
Goldfarb, D.; Qin, Z. Robust low-rank tensor recovery: Models and algorithms. SIAM J. Matrix Anal. Appl. 2014, 35, 225–253. [Google Scholar] [CrossRef]
Cai, J.F.; Candès, E.J.; Shen, Z. A singular value thresholding algorithm for matrix completion. SIAM J. Optimiz. 2010, 20, 1956–1982. [Google Scholar] [CrossRef]
Boyd, S.; Parikh, N.; Chu, E.; Peleato, B.; Eckstein, J. Distributed optimization and statistical learning via the alternating direction method of multipliers. Found. Trends Mach. Learn. 2011, 3, 1–122. [Google Scholar] [CrossRef]
Rockafellar, R.T. Monotone operators and the proximal point algorithm. SIAM J. Control Optim. 1976, 14, 877–898. [Google Scholar] [CrossRef]
Databases of Lobby and Bootstrap. Available online: http://perception.i2r.a-star.edu.sg/bk_ model/bk_index.html (accessed on 1 November 2015).

Figure 1. Comparison of relative errors between Multi-linear Robust Principal Component Analysis (MRPCA) and Generalized Low-Rank Tensor Recovery (GLRTR) with varying λ.

Figure 2. Relative errors of Generalized Low-Rank Tensor Recovery (GLRTR) with varying τ.

Figure 3. Relative errors for different a and b.

Figure 4. Relative errors for different r and p.

Figure 5. Relative errors versus different Inverse Signal-to-Noise Ratio (ISNR) on sparse noise.

Figure 6. Relative errors versus different Inverse Signal-to-Noise Ratio (ISNR) on Gaussian noise.

Figure 7. Background modeling from lobby video. (a) draws the images with missing entries, (b) plots the recovered background, i.e., the low-rank images, (c) shows the recovered foreground, i.e., the sparse noised images, and (d) displays the recovered images.

Figure 8. Background modeling from bootstrap video. (a) draws the images with missing entries, (b) plots the recovered background, i.e., the low-rank images, (c) shows the recovered foreground, i.e., the sparse noised images, and (d) diplays the recovered images.

Table 1. Comparison of experimental results on 50 × 50 × 50 tensors.

**Table 1.** Comparison of experimental results on 50 × 50 × 50 tensors.
(r, p′)	${‖ ℰ ‖}_{F} / {‖ A ‖}_{F}$ (%)	${‖ G ‖}_{F} / {‖ A ‖}_{F}$ (%)	RE of MRPCA (%)	RE of GLRTR (%)
(3, 5)	1253	18.33	9.71 ± 0.98	6.60 ± 0.69
(3, 10)	1769	19.20	8.04 ± 0.84	7.99 ± 0.84
(3, 15)	2218	20.05	9.84 ± 1.97	8.75 ± 1.87
(5, 5)	569	8.79	7.39 ± 0.79	3.74 ± 0.42
(5, 10)	855	9.40	8.55 ± 0.97	4.47 ± 0.52
(5, 15)	1040	9.16	9.13 ± 0.52	5.05 ± 0.31

Table 2. Comparison of experimental results on 20 × 20 × 20 × 20 tensors.

**Table 2.** Comparison of experimental results on 20 × 20 × 20 × 20 tensors.
(r, p′)	${‖ ℰ ‖}_{F} / {‖ A ‖}_{F}$ (%)	${‖ G ‖}_{F} / {‖ A ‖}_{F}$ (%)	RE of MRPCA (%)	RE of GLRTR (%)
(2, 5)	1722	24.05	29.24 ± 7.95	13.90 ± 4.68
(2, 10)	2823	37.72	48.07 ± 18.47	23.98 ± 9.66
(2, 15)	3052	27.17	44.71 ± 13.66	23.10 ± 8.52
(4, 5)	448	6.97	12.00 ± 1.60	8.96 ± 1.15
(4, 10)	563	6.04	15.12 ± 1.98	6.26 ± 0.87
(4, 15)	758	6.58	24.22 ± 2.61	12.16 ± 2.90

Table 3. Comparison of experimental results on incomplete 50 × 50 × 50 tensors for r = 5.

**Table 3.** Comparison of experimental results on incomplete 50 × 50 × 50 tensors for r = 5.
(p, λ)	(30,0.07)	(50, 0.05)	(70, 0.03)	(90,0.03)
RE of MRPCA (%)	87.81 ± 0.27	77.02 ± 0.46	60.64 ± 1.77	11.49 ± 1.73
RE of GLRTR (%)	11.58 ± 1.72	6.62 ± 0.82	5.51 ± 0.83	4.27 ± 0.47

Table 4. Comparison of experimental results on incomplete 20 × 20 × 20 × 20 tensors for r = 4.

**Table 4.** Comparison of experimental results on incomplete 20 × 20 × 20 × 20 tensors for r = 4.
(p, λ)	(30, 0.03)	(50, 0.035)	(70, 0.03)	(90, 0.03)
RE of MRPCA (%)	88.38 ± 0.25	78.94 ± 0.34	65.06 ± 0.59	40.03 ± 1.34
RE of GLRTR (%)	62.51 ± 1.72	20.09 ± 2.63	9.14 ± 0.12	8.58 ± 1.11

© 2016 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Shi, J.; Yin, Q.; Zheng, X.; Yang, W. Alternating Direction Method of Multipliers for Generalized Low-Rank Tensor Recovery. Algorithms 2016, 9, 28. https://doi.org/10.3390/a9020028

AMA Style

Shi J, Yin Q, Zheng X, Yang W. Alternating Direction Method of Multipliers for Generalized Low-Rank Tensor Recovery. Algorithms. 2016; 9(2):28. https://doi.org/10.3390/a9020028

Chicago/Turabian Style

Shi, Jiarong, Qingyan Yin, Xiuyun Zheng, and Wei Yang. 2016. "Alternating Direction Method of Multipliers for Generalized Low-Rank Tensor Recovery" Algorithms 9, no. 2: 28. https://doi.org/10.3390/a9020028

APA Style

Shi, J., Yin, Q., Zheng, X., & Yang, W. (2016). Alternating Direction Method of Multipliers for Generalized Low-Rank Tensor Recovery. Algorithms, 9(2), 28. https://doi.org/10.3390/a9020028

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Alternating Direction Method of Multipliers for Generalized Low-Rank Tensor Recovery

Abstract

1. Introduction

2. Notations and Preliminaries

3. Related Works

4. Generalized Low-Rank Tensor Recovery

4.1. Model of GLRTR

4.2. Optimization Algorithm to GLRTR

5. Convergence Analysis

6. Experimental Results

6.1. Synthetic Data

6.2. Influence of Noise and Sampling Rate on the Relative Error

6.3. Applications in Background Modeling

7. Conclusions

Acknowledgments

Author Contributions

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI