1 Introduction

There has been growing interest in the robust principal component analysis (RPCA)^{[1, 2]} recently,namely,recovering a low-rank matrix with an unknown fraction whose entries are arbitrarily corrupted. This problem has wide applications in many areas such as web data ranking and collaborative filtering^[3],computer vision^{[1, 4]},and machine learning^[5].

Mathematically,suppose that we are given a large data matrix D,which is the sum of a low-rank matrix L₀ and a sparse matrix S₀. The aim of the RPCA is to recover the low-rank matrix L₀ from the highly corrupted measure D = L₀ + S₀ exactly and efficiently. Cand`es et al.^[1] and Chandrasekaran et al.^[2] have shown that under some conditions,the RPCA problem can be exactly solved via convex optimization by minimizing a linear combination of the nuclear norm and the l₁norm. That is to solve

where ||·|| denotes the nuclear norm of a matrix,i.e.,the sum of its singular values,||·||₁ denotes the sum of the absolute values of the matrix entries,and λ is a positive weighting parameter.

The convex optimization (1) has received much attention,and many numerical algorithms have been proposed. Liu and Vandenberghe^[6] suggested the interior point method,with which the problem (1) could be easily recast as a semidefinite program and be solved by some in- terior point solvers^[7]. However,since the interior point method depends on second-order in- formation which is difficult to be computed,it is limited to small problems,e.g.,n < 10². Therefore,it cannot handle matrices whose dimension is larger than 102. To overcome the limited scalability of the interior point method,many kinds of first-order methods are devel- oped. Nesterov^[8] presented an accelerated first-order method for optimization,which has been successfully employed for improving the convergence of the iterative thresholding method in l₁ norm minimization^{[9, 10]}. Toh and Yan^[11] designed a proximal gradient algorithm (APG) for the matrix completion inspired by an accelerated version of the iterative thresholding. Ganesh et al.^[12] applied the APG algorithm to low-rank and sparse decomposition,i.e.,the problem (1). Experimental results show that the APG is at least 50 times faster than the iterative thresholding method. Goldfarb et al.^[13] developed the accelerated iterative thresholding and proposed an alternating linearization method and its fast version. They also applied their algo- rithms to the problem (1) and obtained higher accuracy than the APG with the same number of iterations. More recently,Lin et al.^[14] applied the augmented Lagrange multiplier (ALM) for the problem (1),and achieved much higher efficiency,at least five times faster than the APG. Yuan and Yang^[15] proposed the alternating direction method based on the ALM.

However,the above mentioned methods depend strongly on a good parameter-choice rule since they all convert the constraint into the objective function as a penalty term and have to introduce some penalty parameters. A common way is to use the so-called continuation schemes^{[16, 17]} which need to be carefully designed. In this paper,we do not follow the traditional way but transform the problem (1) as an unconstrained version by canceling one variable with the constrained equation and minimize a new objective function for the other variable. That is to minimize

Obviously,the objective function of (2) or (3) is without any introduced parameter,and thus needs no continuation scheme. Based on the Douglas-Rachford splitting method,we design a proximity point algorithm to solve the problem (2) for the exact low-rank matrix recovered from the sparse noise corrupted data matrix D. Although the aim of our algorithm is to recover the low-rank matrix,we shall show that the sparse component S₀ can also be recovered simultaneously by a proximal point sequence.

The remainder of the paper is organized as follows. In Section 2,we introduce the prox- imity operator and the Douglas-Rachford splitting method. In Section 3,we first discuss the uniqueness of the solution to (2) or (3),and then apply the Douglas-Rachford splitting method to the unconstrained problem (2) and design a novel algorithm. Meanwhile,the convergence is proved. In Section 4,we present some implementation details and simulations to illustrate the practical utility of our algorithm. Finally,some discussion is made in Section 5. 2 Method 2.1 Proximity operator

Let H be a real Hilbert space with the inner product <·|·> and a corresponding norm ||·||. Let Γ₀(H) be a class of lower semicontinuous,convex,and proper functions from H to (−∞,+∞]. The proximity operator of a function f ∈ Γ₀(H) is the operator prox_f : H → H,which maps every x ∈ H to the unique minimizer of the function f + 1/2||x − ·||²,i.e.,for all x ∈ H,

The foundational theory of proximity operators was introduced by Moreau^[18] in 1962. Recently, Combettes and Wajs^[19] developed the theory and applied it to signal recovery problems. Here, we provide two useful properties of proximity operators^[20]. One is a subdifferential characteri- zation,and the other is the nonexpansive property. Both of them are essential to establish the convergence of various algorithms based on proximity operators.

Lemma 1 Let f ∈ Γ₀(H). Let x ∈ H and p ∈ H. Then,p = prox_fx if and only if x − p ∈ ∂f(p).

Lemma 2 Let f ∈ Γ₀(H). Then,prox_f is nonexpansive: x,y ∈ H,it holds

For convenience,we introduce the following soft-thresholding operator:

where x ∈ R,and μ > 0. This operator can be extended to vectors and matrices by applying its element-wise. Denote With the soft-thresholding operator,we can provide closed-form formulas for proximity opera- tors of g1 and g2 as follows^[12]: 2.2 Douglas-Rachford splitting method

To solve the following unconstrained convex problem:

where both f₁ and f₂ belong to Γ₀(H). The Douglas-Rachford splitting method^[20] for mini- mizing F(X) is governed by an updating rule as follows: where λ_n ∈ (0,2),r ∈ (0,+∞),and the sequences (A_n)_n∈N and (B_n)_n∈N model tolerances in the implementation of the proximity operators.

Here,we provide a brief summary of notations for describing a convergence result of the above updating rule. The domain of a function f : H → (−∞,+∞] is

Theorem 1^[20] Let fi ∈ Γ₀(H) (i ∈ {1,2}). Let (A_n)_n∈N and (B_n)_n∈N be sequences in H. Let (y_n)_n∈N be a sequence generated by the updating rule (8). Suppose that the following conditions hold: (i) f₁ + f₂ is coercive, (ii) 0 ∈ sri(domf₁ − domf₂), 3 Algorithm 3.1 Proximity point algorithm

The RPCA demonstrates that the solving of the problem (1) can exactly recover the low-rank and sparse components L₀ and S₀. In other words,the solution to the problem (1) is unique. Thereby,we conclude that the solution to (2) or (3),which is equivalent to the problem (1), is also unique. Hence,one can deduce that the first-order optimality condition of (2) or (3) is also sufficient for a point to be a solution. That is the following lemma.

Lemma 3 X is a solution to the problem (2) if and only if it obeys

Y is a solution to the problem (3) if and only if it obeys

Now,we focus on the Hilbert space R^n₁×n₂ with n₁ and n₂ being the dimension of the data matrix D. Let

Before applying the Douglas-Rachford splitting method to the problem (2),we need to compute the proximity operator of f₁(X),i.e.,

Let Z = D − X. Then,we have X = D − Z and

By the closed-form formula for the proximity operator of g₁ = ||X||₁,we deduce that

Applying the Douglas-Rachford splitting method to the problem (2) yields the following updating rule:

where λ_n ∈ (0,2),r ∈ (0,+∞),and the sequences (A_n)_n∈N and (B_n)_n∈N model tolerances in the implementation of the proximity operators. Replacing the proximity operators by their closed-form formulas into the above updating rule leads to our main algorithm of this study, which is described as follows: Algorithm 1 Proximity point algorithm for RPCA Input: D ∈R^n₁×n₂ and λ. (I) Let λ_n ∈ (0,2),r ∈ (0,+∞),Y₀ = D. (II) (U,Σ,V ) := svd(Y_n),P_2,n := USr[Σ]V ^T. (III) X_n := Y_n − 2(Y_n − P_2,n). (IV) Z_n := D − X_n,P_1,n := Srλ[Z_n],P_1,n := D − P_1,n. (V) Y_n+1 := Y_n + λ_n(P_1,n − P_2,n). (VI) n := n + 1. (VII) If the convergence is satisfied,end the computation; otherwise,Y_n+1 := Y₀,and go to Step (II). Output: (P_1,n,P_2,n). 3.2 Convergence analysis

From Theorem 1,we can see that the convergence of the sequence generated by the Douglas- Rachford splitting method is weak. This is because that the discussed space is a general Hilbert space whose dimension may be infinite. But here,our space R^n₁×n₂ is a finite dimensional linear space,and we can deduce a strong convergence result since the weak convergence is equivalent to the strong convergence in any finite dimensional linear space. Furthermore,we shall show that the sparse matrix S₀ can be recovered simultaneously although we do not solve its optimization problem (3).

Theorem 2 Let (A_n)_n∈N and (B_n)_n∈N be sequences in R^n₁×n₂ . Let (Y_n)_n∈N,(P1,n)n∈N, and (P2,n)n∈N be sequences generated by the updating rule (12). We suppose that the sequence (λ_n)_n∈N has a positive lower bound and the following conditions hold:

Proof It is obvious that the sum function f₁ +f₂ is coercive and 0 ∈ sri(domf1 −domf2). Together with the above conditions,the four conditions of Theorem 1 are satisfied. Therefore, we conclude that (Y_n)_n∈N converges weakly to the matrix Y and proxrf2Y is a minimizer of the problem (2). Since the weak convergence turns into a strong convergence in any finite dimensional linear space,(Y_n)_n∈N converges to the matrix Y . By the second condition

and the condition that (λ_n)_n∈N has a positive lower bound,we deduce that both (A_n)_n∈N and (B_n)_n∈N tend to zero as n tends to ∞. Denote By Lemma 2,it holds Then, It follows that (P_1,n)_n∈N converges to the matrix P₁. Since we deduce that D = P₁+P₂ from the fact that the sequence (λ_n)_n∈N has a positive lower bound and sequences (Y_n)_n∈N,(P_1,n)_n∈N,and (P_2,n)_n∈N converge to Y ,P₁,and P₂,respectively. By Lemma 1 and P₂ = prox_rf2Y,we have which follows Y − P₂ ∈ @λrkP1k1 since D = P₁ + P2. Together with Y − P₂ = ∂r||P₂||,we have Canceling one variable by D = P₁ + P₂ yields 4 Numerical result 4.1 Implementation details 4.1.1 Computing SVD

Since we only need to compute few singular values that are larger than a set threshold and their corresponding singular vectors,we use the PROPACK package^[21],which has been employed and recommended by several papers^{[11, 14, 22]}. Before using the PROPACK package, we have to predict the dimension of the principal singular space whose singular values are larger than the set singular. Since the sequences (A_n)_n∈N and (B_n)_n∈N model tolerances in the implementation of the proximity operators and only need to obey the second condition of Theorem 2,the convergence of our algorithm still can be guaranteed by Theorem 2 even if the prediction is less than the practical dimension. This is not the case of other algorithms. We divide our prediction rule into two cases. One case is that we can predict the rank of the low-rank matrix L₀,and the other case is that we have no idea of the rank of the low-rank matrix L₀. Here,we use the same parameter notations in Refs. ^{[11, 14]}. In the first case,we use svpn,which equals to or is somewhat bigger than the rank of the low-rank matrix L₀,as the predicted dimension. In the second case,we use the following prediction rule:

where d = min(n₁,n₂),s is a constant number to avoid our prediction being conservative,svn is the predicted dimension,and svpn is the number of singular values in the svn singular values which are larger than a given threshold. The second rule is similar to the prediction rule^{[11, 14]} which is more complicated and depends on the observation that Lagrangian multiplier matrices are monotonically increasing and become stable at the true rank. Cai et al.^[22] observed the same phenomenon in matrix completion by the singular value thresholding algorithm. However, no person can prove it strictly up to now. Therefore,if the computation of the soft-thresholding operator by their prediction rules is not exact,the convergence of their algorithms needs to be proved by taking the error into account. 4.1.2 Order of minimizing f₁ and f2

As suggested by Lin et al.^[14],we first minimize the function f1,and then minimize the function f₂ in our numerical simulations,considering the huge complexity of the SVD for large dimensional matrices. One can easily check that the order of minimizing does not affect the convergence result. 4.1.3 Setting parameters

Following the theory of RPCA in ^[1],we choose a fixed weighting parameter

for a given data matrix D ∈ R^n₁×n₂ to be decomposed. We set λ_n ≡ 1.5 and r = 20. The stopping criterion is ||D − P_1,n − P_2,n||_F /||D||_F < 10⁻⁷. 4.2 Simulations

In this section,we compare our algorithm,i.e.,the proximity point algorithm (PPA),with the exact augmented Lagrange multiplier (EALM) algorithm and the inexact augmented La- grange multiplier (IALM) algorithm proposed in ^[14]. Both EALM and IALM are faster than the previous state-of-the-art algorithms^[12]. We use the same simulation conditions in ^[14] for impartial comparison,i.e.,we use randomly generated square matrices for our simulations. We create the random low-rank matrices L₀ ∈ R^m×m with the rank r by the following procedures: generate the random matrices U ∈ R^m×r and V ∈ R^m×r with independent and identically dis- tributed Gaussian random variables with zero mean and unit variance,and then set L₀ = UV^T.We generate S₀ as a sparsematrix,whose support is chosen uniformly at random and whose non- zero entries are independent and identically distributed uniformly in the interval [−500,500]. Thus,the ordered pair (L0,S₀) is the true solution. The matrix D = L₀ + S₀ is the input for decomposition with different algorithms and (L,S) denotes the output. We adopt the first pre- diction rule discussed in implementation details,and use relL and rel S to denote the relative errors of the output solutions (L,S),i.e.,

All simulations are performed by MATLAB (version 7.01) on a personal computer equipped with an Intel Pentium-IV 3.06GHZ processor and 512MB of RAM under Windows XP. Brief comparisons of the three algorithms are presented in Tables 1 and 2. As can be seen,for the same triplet {m,rank(L₀),||S₀||₀},the RPCA problem is solved for the same data matrix D by using three different algorithms. Moreover,we can see that the three algorithms all recover the low-rank matrix L₀ and the sparse matrix S₀ with high accuracy. IALM is the fastest one among the three algorithms,and EALM is the slowest one. Since IALM is at least five times faster than APG^{[11, 12]} and IALM is about two times faster than PPA,we conclude that PPA is at least two times faster than APG. Moreover,when the rank of L₀ or kS0k0 increases,the estimate of the number of non-zeros in S₀ (i.e.,||S₀||₀) becomes worse and worse. While errors between the estimated ||S₀||₀ by PPA and IALM and the truth number of non-zeros in S₀ are not beyond two. Although IALM performs better than the others,one has to choose penalty parameters^[14] in each iteration step properly to ensure both optimality and fast convergence.

5 Conclusions

In the paper,we design an algorithm,i.e.,PPA,based on the Douglas-Rachford splitting method for solving the RPCA problem. Compared with the previous state-of-the-art algorithms, our algorithm has at least three advantages. First,the objective function of the problem (2) or (3) is simpler than that of other algorithms,and has only one variable for the low-rank matrix variable L and the sparse matrix variable S. Secondly,it recovers the low-rank and sparse components exactly and simultaneously by solving the problem (2) or (3),and its convergence can be proved strictly. Thirdly,it does not introduce any penalty parameter. In contrast, other algorithms,for example APG and ALM,introduce such penalty parameters for canceling constraint,and hence have to carefully design the continuous schemes to ensure both optimality and fast convergence. Numerical simulations demonstrate the practical utility of our algorithm.