A Nonconvex Method to Low-rank Matrix Completion

In recent years, the problem of recovering a low-rank matrix from partial entries, known as low-rank matrix completion problem, has attracted much attention in many applications. However, it is a NP-hard problem due to the nonconvexity nature of the matrix rank function. In this paper, a rank Laplace function is studied to recover the low-rank matrices. Firstly, we propose an iterative Laplace thresholding algorithm to solve the regularized Laplace low-rank matrix completion problem. Secondly, some other iterative thresholding algorithms are designed to recover the low-rank matrices. Finally, we provide a series of numerical simulations to test the proposed algorithms on some low-rank matrix completion and image inpainting problems, and the results show that our algorithms perform better than some state-of-art methods in recovering the low-rank matrices.


I. INTRODUCTION
I N recent years, the low-rank matrix completion problem has attracted much attention in many applications such as video denoising [1], [2], image inpainting [3], [4], system identification [5], and so on. In fact, for a low-rank matrix M ∈ R m×n one would like to know as precisely as possible. However, the only information available about M is a sampled set of entries M i, j , (i, j) ∈ Ω, where Ω ⊂ {1, 2, · · · , m} × {1, 2, · · · , n} is the set of indices of samples. The goal of low-rank matrix completion problem is to recover the low-rank matrix from its partial entries M i, j , (i, j) ∈ Ω. In mathematics, the low-rank matrix completion problem can be modeled as the following minimization problem: min X∈R m×n rank(X) s.t. X i, j = M i, j , ∀ (i, j) ∈ Ω. (1) Throughout this paper, without loss of generality, we assume that m ≤ n. For convenience, if we summarize the information available via P Ω (M), where P Ω : R m×n → R m×n is a sampling operator that keeps the entries in Ω unchanged and those outside Ω zeros, i.e., the low-rank matrix completion problem (2) can be reformulated as min X∈R m×n rank(X) s.t. P Ω (X) = P Ω (M).
Unfortunately, due to the nonconvexity nature of the matrix rank function rank(X), problem (3) is NP-hard [6], and can not be solved efficiently by the algorithms. In fact, all known algorithms for exactly solving it are doubly exponential in theory and in practice. To avoid such a difficulty, some researchers introduced the convex surrogate function of rank(X), namely, nuclear-norm X * of X, to relax the rank of matrix X, which leads to the following nuclear-norm lowrank matrix completion problem [6]- [12]: where X * = m i=1 σ i (X) denotes the matrix nuclear norm and σ i (X) is the i-th largest singular value of matrix X ∈ R m×n . As a convex relaxation for the NP-hard problem (3), problem (4) and its regularized version possess many algorithmic advantages. The representative algorithms include semidefinite programming and interior point SDP VOLUME 4, 2016 solver [11], [13], singular value thresholding (SVT) algorithm [7], accelerated proximal gradient (APG) algorithm [14], proximal point algorithms [15], fixed point and Bregman iterative algorithms [8], linearized augmented Lagrangian method and alternating direction method [16], and so on. Although the problem (4) and its regularized version possess many algorithmic advantages, nuclear-norm always performs not very well in recovering the low-rank matrices. In fact, nuclear-norm always needs more observations to recover a real low-rank matrix. Motivated by the good performance of the surrogate functions in compressed sensing, in this paper, we are interested in using the the following Laplace function [17]: to define a rank Laplace function to relax the rank function rank(X) in problem (3), and generating a new low-rank matrix completion model to recover the low-rank matrices. The parameter γ controlling the behavior of the function f γ (t), and the plot of the function f γ (t) is shown in Figure 1.
It is easy to verify that lim Therefore, the function f γ (t) can be used to relax the rank function rank(X). In this paper, the rank Laplace function is defined on the singular values of the matrix X ∈ R m×n , i.e., and it can be considered as an ideal relaxation for the rank function rank(X) because of the relationship In this way, the idea of using the rank Laplace function defined in (6) to relax the rank function rank(X) comes very naturally. Using the rank Laplace function |||X||| γ to relax the rank function rank(X) in problem (3) leads to the following Laplace low-rank matrix completion problem: Usually, solving problem (8) can be transformed into solving its regularized version which is described as the following regularized Laplace low-rank matrix completion problem: where λ > 0 is the regularized parameter. The rest of the paper is organized as follows. In Section II, we present some useful notions and crucial preliminary results that are used in this paper. In Section III, an iterative Laplace thresholding algorithm is proposed to solve problem (9). In Section IV, three types of iterative difference Laplace thresholding algorithms are studied for low-rank matrix completion. In Section V, we provide a series of numerical simulations to test our algorithms on some low-rank matrix completion and image inpainting problems. We conclude this paper in Section VI.

II. NOTIONS AND PRELIMINARY RESULTS
In this section, we present some useful notions and preliminary results that are used in this paper.

A. NOTIONS
The space of m × n real matrices is denoted by R m×n . Given any matrices X, Y ∈ R m×n , X, Y = Tr(Y X) denotes the inner product of matrices X and Y, where Tr(Y X) denotes the trace of the matrix Y X. The Frobenius norm of matrix X ∈ R m×n is denoted by X F = √ Tr(X X). The singular value decomposition (SVD) of matrix X ∈ R m×n is defined as where U is a m × m unitary matrix, V is a n × n unitary matrix, diag({σ i (X)} 1≤i≤m ) ∈ R m×m is a diagonal matrix of the singular values of matrix X arranged in descending order σ 1 (X) ≥ σ 2 (X) ≥ · · · ≥ σ m (X) ≥ 0, and 0 m×(n−m) ∈ R m×(n−m) is a m × (n − m) zero matrix.

B. PRELIMINARY RESULTS
For any t ∈ R, consider the following minimization problem where γ > 0 and λ > 0. For any x ≥ 0, the derivative of θ t (x) is given by Similar to [ [17], Theorem 1], we can get the following key lemma. Theorem 1. Suppose γ > √ λ and x = g λ,γ (t) is the implicit function determined by the equation [θ t (x)] = 0. Let x * be the global optimal solution of problem (10). Then we can get that where t λ,γ = λ γ . Definition 1. (Laplace thresholding function) Suppose t ∈ R, for any λ > 0, the function h λ (t) defined in (12) is called the Laplace thresholding function.
For any λ > 0, we defined the rank Laplace thresholding operator G λ : R m×n → R m×n as Definition 4. (See [18], [19]) T : R m×n → R + is a smooth function of type C 1,1 , i.e., the gradient is Lipschitz continuous, for every X, Z ∈ R m×n , where L(T ) > 0 is the Lipschitz constant of ∇T .
for every X, Z ∈ R m×n .
Moreover, we also need the definition of the singular value shrinkage operator [7], and it simply applies the soft thresholding operator [21], [22] defined on vector to the singular values of a matrix.
then the soft thresholding operator D λ (·) can be expressed as where ν + is the positive part of ν, and ν + = max(ν, 0).

III. ITERATIVE LAPLACE THRESHOLDING ALGORITHM FOR LOW-RANK MATRIX COMPLETION
In this section, the proximal gradient descent algorithm is utilized to solve problem (9) and generate an iterative thresholding algorithm for low-rank matrix completion. Let Then the objective function (9) can be rewritten as Clearly, G(X) is a smooth function of type C 1,1 . Therefore, by Lemma 1, for any Z ∈ R m×n , we have where µ ≥ 1.
In this paper, we use the function to approximate the objective function of problem (9). The function P(X, Z) can also be seen as the linearize of G(X) at Z and add the term λ|||X||| γ , and it can be rewritten as By above approximation, problem (9) can be approximated as the following minimization problem: Ignoring constant terms g(Z) and − 1 2µ ∇g(Z) 2 F , solving problem (24) is equivalent to solving the following minimization problem: where . We evaluate the global optimal solution of problem (24) iteratively by setting Z = X k , which leads to the following iterative scheme: This means that, in each iteration, if the following minimization problem: can be solved, we can easily generate an iterative thresholding algorithm to solve problem (9). According to Theorem 1 and Definition 3, problem (27) permits a thresholding representation theory for its global optimal solution.
Theorem 2. For any γ > λ µ , the global optimal solution X k+1 of the problem (27) can be expressed as where H λ µ (·) is obtained by replacing λ with λ µ in H λ (·). Proof. By von Neumann's trace inequality, we have the inequality Notice that the equality holds if the matrix X admits the singular value decomposition where U B and V B are the left and right orthonormal matrices in the SVD of matrix B µ (X k ). Therefore, solving problem (27) reduces to solve the following minimization problem: The objective function in (30) is separable, solving the minimization problem (30) is equivalent to solving the following n subproblems, for i ∈ {1, 2, · · · , n}, By Theorem 1, for γ > λ µ , the global optimal solution of problem (31) can be expressed as where h λ µ (·) is obtained by replacing λ with λ µ in g λ (·). Therefore, it follows that is the global optimal solution of problem (27).
Initialized X 0 ∈ R m×n , iteration (28) describes an iterative thresholding algorithm to solve problem (9). In this paper, we call this method the iterative Laplace thresholding algorithm (ILTA), and we depict it in Algorithm 1.
2) The sequence {X k } is asymptotically regular, i.e., lim k→∞ X k+1 − X k 2 F = 0. Moreover, the sequence {X k } is bounded and has at least one accumulation point.
Although the ILTA has been generated for low-rank matrix completion in this section, the closed form rank Laplace thresholding operator H λ µ (·) in ILTA is not available. In fact, [θ t (x)] = 0 is a transcendental equation and it has no explicit solution, which means that there is no explicit expression for h λ (·) and ILTA in fact is not used to recover the matrices in practical applications. In the following, we will study some new algorithms for problem (9) to recover the lowrank matrices, which can avoid the confusion caused by the unavailability of rank Laplace thresholding operator H λ µ (·) in ILTA.

IV. ITERATIVE DIFFERENCE LAPLACE THRESHOLDING ALGORITHMS FOR LOW-RANK COMPLETION
In this section, different from the ILTA, three new algorithms for problem (9) are studied to recover the low-rank matrices.

A. ITERATIVE DIFFERENCE LAPLACE THRESHOLDING ALGORITHM
Applying the idea of DCA [24], we transform the objective function where α > 0. In this way, problem (27) can be rewritten as At the current point X k of iteration k, replacing the function N(X) with its affine minorization defined by at the neighborhood of X k , where ∂N(·) represents the subdifferential of N(·) and S k ∈ ∂N(X k ), we can get the following convex minimization problem of the form which generates the algorithm for solving the problem (27), i.e., It is easy to verify that the global optimal solution X k+1 to the problem (35) can be expressed as where T (X k ) = B µ (X k ) + S k and D αλµ (·) is obtained by replacing λ with αλµ in D λ (·).

Remark 1. Let
be the SVD of matrix X k ∈ R m×n . Then, we can set In this paper, we select The more knowledge for the subdifferential of the matrix function can be seen in [25].
For k = 0, 1, 2, · · · , if we threshold the singular values of matrix D αλµ (T (X k )) by the threshold value t λ µ ,γ = λ µγ , we can get the iteration In this paper, we use iteration (37) to recover the low-rank matrices, and we call iteration (37) the iterative difference Laplace thresholding algorithm (IDLTA). The IDLTA for low-rank matrix completion is depicted in Algorithm 2.

B. ADAPTIVE ITERATIVE DIFFERENCE LAPLACE THRESHOLDING ALGORITHM
It is well known that the quality of the solution of a regularized problem depends seriously on the setting of the regularized parameter λ. In general, there is no optimal rule for the selection of proper regularized parameter λ in a regularized problem. In IDLTA, we set in each iteration, where r is the matrix rank of the global optimal solution to the problem (9). When doing so, the IDLTA will be adaptive and free from the choice of the regularized parameter λ in each iteration. In this paper, we call IDLTA with parameter-setting strategy (38) the adaptive iterative difference Laplace thresholding algorithm (AIDLTA), and we depict it in Algorithm 3.

C. FAST ADAPTIVE ITERATIVE DIFFERENCE LAPLACE THRESHOLDING ALGORITHM
In this section, we discuss how to accelerate the proposed AIDLTA for low-rank matrix completion. Inspired by the modified fast iterative shrinkage-thresholding algorithm (FISTA) proposed in [26], which is generated to minimize the problems of the sum of a convex function with Lipschitz continuous gradient and a convex function that its proximity operator is efficiently calculated, we proposed a fast version of AIDLTA for low-rank matrix completion.
When the modified FISTA is applied to solve problem (9) for low-rank matrix completion, it generates the following iterative scheme: where p, q > 0 and ς ∈ (0, 4]. Combining the iterative scheme of (39) into AIDLTA, a fast version of AIDLTA can be given by where p, q > 0, ς ∈ (0, 4] and S k is defined in Remark 1. We call iteration (40) the fast adaptive iterative difference Laplace thresholding algorithm (FAIDLTA), and we depict it in Algorithm 4.

V. NUMERICAL SIMULATIONS
In this section, we present some numerical simulations to test the performance of the proposed algorithms and compare them with some state-of-art methods. Some completion of random matrices and grayscale image inpainting problems are considered in our numerical simulations. In these numerical experiments, we only consider to test the performance of AIDLTA and FAIDLTA, and compared them with TS1-s2 [4], SVP [27] and APG [14].
In all numerical simulations, we set µ = 1.5, α = 0.1 and γ = 2 in AIDLTA, and set µ = 1.5, α = 0.1, γ = 2, p = 0.01, q = 0.2 and ς = 4 in FAIDLTA. The stopping criterion is defined as and the accuracy of the generated solution X opt is measured by the relative error (RE), defined as The numerical simulations are all performed in MATLAB on a personal computer (Intel(R) Core(TM)i5-9500 CPU @ 3.00GHZ, 8.00GB RAM)

A. COMPLETION OF RANDOM MATRICES
In this subsection, we present the numerical simulations to test the performance of AIDLTA and FAIDLTA for the completion of randomly generated matrices of size n × n and rank r, and compared them with TS1-s2, SVP and APG. In numerical simulations, the matrix M ∈ R n×n of rank r is generated as a product M 1 M 2 , where M 1 ∈ R n×r and M 2 ∈ R r×n have independently identically distributed Gaussian entries. The index set Ω of the given s entries is sampled uniformly at random. The sampling ratio is measured by sr = s/n 2 . The freedom ratio is defined as f r = s/(r(2n − r)), which is one quantity that helps to quantify the difficulty of a recovery problem. If f r < 1, there is always an infinite number of matrices of rank r with the observed entries, so it is impossible to recover the original low-rank matrix in this situation.
The numerical results of the algorithms for the completion of random matrices are reported in Table 1 and Table 2. We can see that AIDLTA and FAIDLTA perform better than TS1-s2, SVP and APG, and the convergence speed of FAIDLTA is faster than AIDLTA.

B. IMAGE INPAINTING
In this subsection, we demonstrate the performance of AIDL-TA and FAIDLTA on the image inpainting problems, and compared them with TS1-s2, SVP and APG. These algorithms are tested on two standard 256 × 256 grayscale images (Peppers and Cameraman). We use the SVD to obtain their approximated images with rank r = 30. The original images and their corresponding approximated low-rank images are displayed in Figure 2 and Figure 3, respectively.   The numerical results of the algorithms for image inpainting problems are reported in Table 3. We can see that AIDLTA and FAIDLTA perform better than TS1-s2, SVP and APG in image inpainting problems, and the convergence speed of FAIDLTA is faster than AIDLTA. For sr = 0.40, we display the recovered low-rank Peppers and Cameraman images via the five algorithms in Figure 4 and Figure 5,  respectively.

VI. CONCLUSION
In this paper, we use a rank Laplace function to approximately substitute the rank function in the NP-hard low-rank matrix completion problem, and design some algorithms to recover the low-rank matrices. Numerical results on some low-rank VOLUME 4, 2016