# A Proximal Fully Parallel Splitting Method for Stable Principal Component Pursuit

Authors: Hongchun Sun; Jing Liu; Min Sun
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

As a special three-block separable convex programming, the stable principal component pursuit (SPCP) arises in many different disciplines, such as statistical learning, signal processing, and web data ranking. In this paper, we propose a proximal fully parallel splitting method (PFPSM) for solving SPCP, in which the resulting subproblems all admit closed-form solutions and can be solved in distributed manners. Compared with other similar algorithms in the literature, PFPSM attaches a Glowinski relaxation factorη∈3/2,2/3 to the updating formula for its Lagrange multiplier, which can be used to accelerate the convergence of the generated sequence. Under mild conditions, the global convergence of PFPSM is proved. Preliminary computational results show that the proposed algorithm works very well in practice.


---

## Body

## 1. Introduction

Recovering the low-rank and sparse components of a given matrix is a fundamental task in many scientific fields, which can be modeled as the following two-block separable nonconvex and nonsmooth programming:(1)minL,Srank⁡L+τS0,s.t.L+S=D,where D∈Rp×q is a given matrix, rankL is the rank of the matrix L∈Rp×q, S0 denotes the l0-norm of S, which is defined by the number of the nonzeros of its entries, and τ>0 is a trade-off parameter balancing the low rank and sparsity. Recent studies [1–3] indicate that the NP-hard task (1) can be accurately accomplished by solving one of its convex relaxation problems: the principal component pursuit (PCP) problem, in which the rank and the l0-norm are replaced by their tightest convex relaxations, that is, the nuclear norm and the l1-norm, respectively. The PCP problem can be mathematically modeled as the following optimization problem:(2)minL,SL∗+τS1,s.t.L+S=D,where the nuclear norm L∗ is the sum of all singular values of L, which is to induce the low-rank component of the given matrix D, and the l1-norm S1 is the sum of the absolute values of all entries of S, which is to induce the sparse component of the given matrix D.To capture even more applications, Zhou et al. [4] considered that the observation matrix D is corrupted by both small entry-wise noise and gross sparse errors. Then more general case was studied by Tao and Yuan [5] wherein only a fraction of entries of the given matrix D can be observed, and they suggested recovering the low-rank and sparse components of D by solving the following stable principal component pursuit (SPCP) problem:(3)minL,SL∗+τS1,s.t.PΩD-L-S≤δ,where δ>0 is related to the impulsive or Gaussian noise level, · denotes the Frobenius-norm, Ω⊆1,2,…,p×1,2,…,q is the index set of the observable entries of D, and PΩ:Rp×q→Rp×q is the projection operator defined by (4)PΩXij=Xij,ifi,j∈Ω0,otherwise,1≤i≤p,1≤j≤q.In [6], Hou et al. transformed (3) into the following unconstrained problem:(5)minL,SL∗+τS1+12μPΩD-L-S2,where μ>0 is a penalty parameter. Problems (3) and (5) are both convex optimization programming, which can be equivalently reformulated as some semidefinite programming (SDP) problems and then solved by using some state-of-the-art SDP solvers in polynomial time, such as SeDuMi [7] and SDPT3 [8]. However, these solvers are not applicable to high-dimensional (3) and (5), because a linear system needs to be solved (in)exactly at each iteration, which may be more costly and require very large amount of memory in actual implementations [9].By introducing an auxiliary variableU∈Rp×q, the model (5) can be rewritten as the following equivalent form:(6)minL,S,UL∗+τS1+12μPΩU2,s.t.L+S+U=D,which decouples the term PΩD-L-S2 of (5). Recently, there has been an intensive study on the iteration methods for solving (6) [5, 6, 10–12]. Interestingly, these methods usually only exploit the first-order information of the objective function of (6), and most of them are based on the well-knownalternating direction method withmultipliers (ADMM), such as the sequential ADMM [11], the partially parallel ADMMs [5, 6, 10, 12], and the fully parallel ADMMs [13, 14]. As an influential first-order iteration method in optimization, ADMM was independently proposed by Glowinski and Marrocco [15] and Gabay and Mercier [16] about forty years ago, and due to its high efficiency in the numerical simulation, ADMM has been receiving considerable attention in recent years. ADMM can be viewed as an application of the Douglas-Rachford splitting method to the dual of the two-block separable convex programming (2) [17] or a special case of the proximal point method for the general convex programming [18], or a generalization of the classical Uzawa method for solving the saddle-point problems [19]. Chen et al. [20] showed that the sequential ADMM cannot directly extend to three-block separable convex programming by a simple counterexample, and similarly He et al. [13] showed that the full parallel ADMM may be divergent even for two-block separable convex programming. There are three strategies to deal with the issue: (1) the substitution procedure [13, 21]; (2) the technique of regularization [14, 22]; (3) the technique of small step size [23]. For other recent developments of ADMM, including the (sub)linear convergence rate, various acceleration techniques, indefinite regularized matrices, larger step size, and its generalization to solve nonconvex and nonsmooth multiblock separable programming, we refer to [24–28].In this paper, we focus the attention on the fully parallel ADMM, which is based on the substitution procedure and is more suitable for distributed computation. Following the full Jacobian decomposition of the augmented Lagrangian method (FJDALM) in [13], we shall propose a new fully parallel ADMM, termed as the proximal fully parallel splitting method (PFPSM) in the following analysis. Compared to other similar methods in the literature, the new iteration method has three properties. The first one is that all of its subproblems are regularized by some proximal terms, which can enhance its numerical efficiency in practice. The second one is that a Glowinski relaxation factor η∈3/2,2/3 is attached to the updating formula for its Lagrange multiplier, and larger values of η often can accelerate the convergence of the proposed method (as to be reported in Section 4). To the best of our knowledge, this is the first fully parallel ADMM-type method with the Glowinski relaxation factor for solving problem (6). The third one is that it dynamically updates the step size αk.The remainder of this paper is organized as follows. In Section2, we present some necessary notations and useful lemmas. In Section 3, we propose the proximal fully parallel splitting method for solving problem (6) and prove its global convergence step by step. Some numerical results about the proposed method are shown and discussed in Section 4. Finally, in Section 5, some concluding remarks are drawn and several future research directions are discussed also.


## 2. Preliminaries

In this section, we briefly introduce some necessary notations, definitions, and lemmas.For any matrixA=aij∈Rp×q, let A=∑i=1p∑j=1qaij2 denote the Frobenius-norm of A. For any two matrices A,B∈Rp×q, 〈A,B〉≔TraceA⊤B denotes their Frobenius inner product. Then, we have A,B≤AB. If G∈Rp×p is symmetric positive definite matrix, we denote by AG=A,GA=TraceA⊤GA1/2 the G-norm of the matrix A. The effective domain of a function f:X→-∞,+∞ is defined as domf≔x∈X∣fx<+∞. The set of all relative interior points of a given nonempty convex set C is denoted by riC. diagA1,A2,…,Am has Ai as its ith block on the diagonal.A functionf:Rn→R is convex if and only if (7)fαx+1-αy≤αfx+1-αfy,∀x,y∈Rn,α∈0,1.Then, for a convex function f:Rn→R, we have the following basic inequality:(8)fx≥fy+ξ,x-y,∀x,y∈Rn,ξ∈∂fy,where ∂fy=ξ∈Rn:fy¯≥fy+ξ,y¯-y,forally¯∈Rn denotes the subdifferential of f· at the point y. The above definition and property (8) can be easily extended to the matrix-valued function.For the convenience of our analysis, let us relabelL as X1, S as X2, and U as X3 and denote θ1X1=L∗,θ2X2=τS1,θ3X3=1/2μPΩU2. Then θi·i=1,2,3 are all convex function.We make the following standard assumption about problem (6) in this paper.Assumption 1.
The generalized Slater’s condition holds; that is, there is a pointX^1,X^2,X^3∈ridomθ1×domθ2×domθ3∩P, where the set P is defined by (9)P=X1,X2,X3∣X1+X2+X3=D.Under Assumption1, it follows from Theorems 3.22 and 3.23 of [29] that X1∗,X2∗,X3∗ is an optimal solution of problem (6) if and only if there exists a matrix Λ∗∈Rp×q such that X1∗,X2∗,X3∗,Λ∗ is a solution of the following KKT systems:(10)0∈∂θiXi∗-Λ∗,i=1,2,3,X1∗+X2∗+X3∗=D.The following lemma is the basis to reformulate (10) as a mixed variational inequality problem.Lemma 2.
For any two matricesXi∗∈Rp×q,Λ∗∈Rp×q, the relationship 0∈∂θiXi∗-Λ∗ is equivalent to the following mixed variational inequality problem:(11)θiXi-θiXi∗+Xi-Xi∗,-Λ∗≥0,∀Xi∈Rp×q.Proof.
From0∈∂θiXi∗-Λ∗, one has (12)Λ∗∈∂θiXi∗,which together with the convexity of θi· and the subgradient inequality (8) implies (13)θiXi-θiXi∗≥Xi-Xi∗,Λ∗,∀Xi∈Rp×q;that is,(14)θiXi-θiXi∗+Xi-Xi∗,-Λ∗≥0,∀Xi∈Rp×q.
Conversely, fromθiXi-θiXi∗+Xi-Xi∗,-Λ∗≥0,∀Xi∈Rp×q, one has (15)θiXi+Xi,-Λ∗≥θiXi∗+Xi∗,-Λ∗,∀Xi∈Rp×q,which implies(16)Xi∗=argminXi∈Rp×q⁡θiXi+Xi,-Λ∗.From this and Theorem 3.22 of [29], we have 0∈∂θi(Xi∗)-Λ∗. The proof is completed.Remark 3.
Based on (10) and Lemma 2, the matrix X∗=X1∗,X2∗,X3∗ is an optimal solution to problem (6) if and only if there exists a matrix Λ∗∈Rp×q such that(17)θiXi-θiXi∗+Xi-Xi∗,-Λ∗≥0,∀Xi∈Rp×q,i=1,2,3;X1∗+X2∗+X3∗=D.Moreover, any Λ∗∈Rp×q satisfying (17) is an optimal solution to the dual of problem (6).Furthermore, we set(18)X=X1,X2,X3,θX=θ1X1+θ2X2+θ3X3,W=X,Λ.Obviously, (17) can be written as the following mixed variational inequality problem, abbreviated as VIW,F,θ: find a matrix W∗=X1∗,X2∗,X3∗,Λ∗∈W such that(19)θX-θX∗+W-W∗,FW∗≥0,∀W∈W,where W=Rp×q×Rp×q×Rp×q×Rp×q, and(20)FW=-Λ-Λ-ΛX1+X2+X3-D=000-I000-I000-IIII0X1X2X3Λ-000D.The set of the solutions of VIW,F,θ, denoted by W∗, is nonempty by Assumption 1 and Remark 3. Because FW is a linear mapping and its coefficient matrix is skew-symmetric, it satisfies the following nice property:(21)W′-W,FW′-FW=0,∀W′,W∈W.


## 3. The Proximal Fully Parallel Splitting Method

In this section, we shall present the proximal fully parallel splitting method for solving (6) and discuss its convergence results.Firstly, the augmented Lagrangian function of the convex programming (6) can be written as (22)LβX1,X2,X3,Λ=∑i=13θiXi-Λ,X1+X2+X3-D+β2X1+X2+X3-D2,where β>0 is the penalty parameter and Λ∈Rp×q is the Lagrange multiplier for the linear constraint of (6). To simplify our notation in later analysis, we define a block diagonal matrix as follows:(23)G=diag⁡βIp+R1,βIp+R2,βIp+R3,Ipηβ,where Ri∈Rp×pi=1,2,3 are some symmetric positive semidefinite matrices and η∈3/2,2/3 is a predefined parameter. Note that the matrix G is symmetric positive definite.Algorithm 4 (proximal fully parallel splitting method (termed as PFPSM)).
⁡
Step
0.Given the  tolerance ε>0, two constants η∈3/2,2/3,γ∈0,2, and the initial iterate W0=X10,X20,X30,Λ0∈W. Set k≔0.
Step
1 (prediction step). Compute the temporary iterate W~k=X~1k,X~2k,X~3k,Λ~k by(24)X~1k≔arg⁡min⁡LβX1,X2k,X3k,Λk+12X1-X1kR12,X~2k≔arg⁡min⁡LβX1k,X2,X3k,Λk+12X2-X2kR22,X~3k≔arg⁡min⁡LβX1k,X2k,X3,Λk+12X3-X3kR32,Λ~k≔Λk-ηβ∑i=13X~ik-D.
Step
2 (convergence verification).If the stopping criterion(25)maxX1k-X~1k,X2k-X~2k,X3k-X~3k,Λk-Λ~k<εis satisfied, then stop.
Step
3 (correction step).Generate the new iterate by(26)Wk+1=Wk-γαkQ-1GWk-W~k,where Q is any symmetric positive definite matrix and αk is the step size defined by(27)αk=φWk,W~kQ-1GWk-W~kQ2,with(28)φWk,W~k=Wk-W~kG2+1ηΛk-Λ~k,∑i=13Xi-X~i+1-ηη2βΛk-Λ~k2.Set k≔k+1 and go to Step 1.Remark 5.
In PFPSM, we apply the classical proximal point algorithm (PPA) to regularize its subproblems. Compared to the iteration method in [13], the new iteration method PFPSM extends the scope of η from 1 to the interval 3/2,2/3, and larger values of η are often more preferred in practice; see Figure 2. Furthermore, compared to the iteration method in [14], the new iteration method PFPSM removes the restriction imposed on the regularized matrices Rii=1,2,3; see condition (2.10) in [14].The following lemma provides a lower bound of theφWk,W~k, that is, the numerator of the step size αk, and this lower bound plays an important role in the proof of the global convergence of PFPSM.Lemma 6.
LetW~k=(X~1k,X~2k,X~3k,Λ~k) be the matrix generated by PFPSM from the given Wk and φ(Wk,W~k) be defined by (28). Then we have(29)φWk,W~k≥c0Wk-W~kG2,where c0=min2η-3/2η,2-3η/2η.Proof.
By the matrix form of the inequality(30)ab≥-12νa2+1νb2with a=Xik-X~ik,b=Λk-Λ~k and ν=3β, one has(31)Xik-X~ik,Λk-Λ~k≥-123βXik-X~ik2+13βΛk-Λ~k2.It follows from (23) and (28) that (32)φWk,W~k=∑i=13Xik-X~ikβIp+Ri2+1η2βΛk-Λ~k2+1ηΛk-Λ~k,∑i=13Xi-X~i≥∑i=13Xik-X~ikβIp+Ri2+1η2βΛk-Λ~k2-12η3β∑i=13Xik-X~ik2+33βΛk-Λ~k2=2η-3β2η∑i=13Xik-X~ik2+∑i=13Xik-X~ikRi2+2-3η2η2βΛk-Λ~k2≥min⁡2η-32η,2-3η2ηβ∑i=13Xik-X~ik2+∑i=13Xik-X~ikRi2+1ηβΛk-Λ~k2=c0Wk-W~kG2,where the first inequality comes from (31) and the second inequality follows from η∈3/2,2/3. The proof is completed.Remark 7.
It follows from (27) and (29) that(33)αk=φWk,W~kQ-1GWk-W~kQ2≥c0λminGGQ-1G,∀k≥0,where λminG is the minimum eigenvalue of the positive definite matrix G. Therefore the sequence αk generated by PFPSM is bounded away from zero. Specially, if we set Q=G, then (34)αk≥c0Wk-W~kG2Wk-W~kG2=c0,∀k≥0.Now, we shall show that the stopping criterion (25) is reasonable.Lemma 8.
IfXik=X~iki=1,2,3,Λk=Λ~k, then the matrix X1k,X2k,X3k,Λk is a solution of VI(W,F,θ).Proof.
Deriving the first-order optimal condition ofX1 subproblem in (24), we have (35)θ1X1-θ1X~1k+X1-X~1k,-Λk+βX~1k+X2k+X3k-D+R1X~1k-X1k≥0,∀X1∈Rp×q.From the updating formula for Λ in (24), one has (36)Λk=Λ~k+ηβ∑i=13X~ik-D,and, substituting it into the above inequality, we get (37)θ1X1-θ1X~1k+X1-X~1k,-Λ~k+β1-ηX~1k+X~2k+X~3k-D+βX2k-X~2k+βX3k-X~3k+R1X~1k-X1k≥0,∀X1∈Rp×q,which can be reduced to (38)θ1X1-θ1X~1k+X1-X~1k,-Λ~k+1-ηηΛk-Λ~k+β∑i=13Xik-X~ik≥X1-X~1k,βIp+R1X1k-X~1k,∀X1∈Rp×q.Similarly, for the other two variables X2,X3, we have (39)θ2X2-θ2X~2k+X2-X~2k,-Λ~k+1-ηηΛk-Λ~k+β∑i=13Xik-X~ik≥X2-X~2k,βIp+R2X2k-X~2k,∀X2∈Rp×q,θ3X3-θ3X~3k+X3-X~3k,-Λ~k+1-ηηΛk-Λ~k+β∑i=13Xik-X~ik≥X3-X~3k,βIp+R3X3k-X~3k,∀X3∈Rp×q.Furthermore, the updating formula for Λ in (24) can be rewritten as the following variational inequality problem: (40)Λ-Λ~k,∑i=13X~ik-D+1ηβΛ~k-Λk≥0,∀Λ∈Rp×q.Then, adding the above four inequalities, and by manipulation, one has(41)θX-θX~k+W-W~k,FW~k+Ip,Ip,Ip⊤1-ηηΛk-Λ~k+β∑i=13Xik-X~ik≥W-W~k,GWk-W~k,∀W∈W.Therefore, substituting Xik=X~iki=1,2,3,Λk=Λ~k into (41), we get (42)θX-θX~k+W-W~k,FW~k≥0,∀W∈W;that is, (43)θX-θXk+W-Wk,FWk≥0,∀W∈W,which indicates that X1k,X2k,X3k,Λk is a solution of VIW,F,θ. This completes the proof.If the iteration method PFPSM terminates at Step2, then the iterate Wk is an approximate solution of VIW,F,θ. Thus, it is assumed, without loss of generality, that the iteration method PFPSM generates an infinite sequence Wk. Now, we intend to prove that -Q-1GWk-W~k is a descent direction of the merit function 1/2W-W∗Q2 at the point W=Wk, where W∗∈W∗.Lemma 9.
LetWk and W~k be the two sequences generated by PFPSM. Then, one has(44)Wk-W∗,GWk-W~k≥φWk,W~k,∀W∗∈W∗.Proof.
SinceW∗∈W∗⊆W, it follows from (41) that (45)θX∗-θX~k+W∗-W~k,FW~k+Ip,Ip,Ip⊤1-ηηΛk-Λ~k+β∑i=13Xik-X~ik≥W∗-W~k,GWk-W~k.Thus,(46)W~k-W∗,GWk-W~k≥θX~k-θX∗+W~k-W∗,FW~k+Ip,Ip,Ip⊤1-ηηΛk-Λ~k+β∑i=13Xik-X~ik≥W~k-W∗,Ip,Ip,Ip⊤1-ηηΛk-Λ~k+β∑i=13Xik-X~ik=∑i=1X~ik-D,1-ηηΛk-Λ~k+β∑i=13Xik-X~ik=1ηβΛk-Λ~k,1-ηηΛk-Λ~k+β∑i=13Xik-X~ik,where the second inequality follows from W∗∈W∗, and the second equality comes from (24). Therefore, we have (47)Wk-W∗,GWk-W~k≥Wk-W~kG2+1ηβΛk-Λ~k,1-ηηΛk-Λ~k+β∑i=13Xik-X~ik=φWk,W~k,which completes the proof of the lemma.From inequalities (29) and (44), one has (48)QWk-W∗,-Q-1GWk-W~k≤-c0Wk-W~kG2,which indicates that -Q-1GWk-W~k is a descent direction of the merit function 1/2W-W∗Q2 at W=Wk. Therefore, it is natural to generate the new iterate Wk+1 along this direction as the iterative scheme (26). Now, let us investigate how to choose a suitable step size to achieve progress as much as possible. The new iterate with an undeterminate step size α is denoted by Wk+1α, so one has(49)Wk+1α=Wk-αQ-1GWk-W~k.Then, let(50)Θkα≔Wk-W∗Q2-Wk+1α-W∗Q2be the profit function of the kth iteration. However, we cannot maximize Θkα directly because W∗∈W∗ is unknown. Motivated by the strategy in [13], we derive a tight lower bound of Θkα in the next lemma, which is independent of the unknown matrix W∗, to obtain a suboptimal step size.Lemma 10.
For any given matrixWk, let W~k=X~1k,X~2k,X~3k,Λ~k be the matrix generated by (24), φWk,W~k be defined by (28), and Θkα be defined by (50). Then, one has(51)Θkα≥Φkα≔2αφWk,W~k-α2Q-1GWk-W~kQ2.Proof.
The proof is similar to that of Lemma4.6 in [13] and is omitted here.It follows from Lemma9 that Φkα is a lower bound of Θkα for any α>0. Therefore, we can maximize Φkα to accelerate the convergence of PFPSM. Note that Φkα is a quadratic function of α and it reaches its maximum at αk defined by (27). The following theorem shows that the sequence Wk generated by PFPSM is Fejér monotone with respect to the solution set W∗ of VIW,F,θ.Theorem 11.
For anyW∗=X1∗,X2∗,X3∗,Λ∗∈W∗, there exists a positive constant η such that the sequence Wk generated by PFPSM satisfies(52)Wk+1-W∗Q2≤Wk-W∗Q2-τWk-W~kG2,where(53)τ=c02γ2-γλminGG⊤Q-1G.Proof.
It follows from (50) and (51) that(54)Wk+1-W∗G2≤Wk-W∗G2-Φkγαk=Wk-W∗G2-2γαkφWk,W~k+γαk2Q-1GWk-W~kQ2=Wk-W∗G2-γ2-γαkφWk,W~k≤Wk-W∗G2-τWk-W~kG2,where the second equality follows from the definition of αk in (27), and the second inequality comes from (29) and (33). Then, assertion (52) is thus proved.Remark 12.
For some very large-scale (6), computing the step size αk defined by (27) is a fussy work needing much calculating. For these instances, set Q=G in PFPSM, and the step size αk can be fixed as the constant c0 defined in Lemma 6. The new iterate can be generated by(55)Wk+1=Wk-γc0Q-1GWk-W~k.Then, referring to the proof of the above theorem, we have (56)Wk+1-W∗G2≤Wk-W∗G2-Φkγc0=Wk-W∗G2-2γc0φWk,W~k+γc02Wk-W~kG2=Wk-W∗G2-2γc02Wk-W~kG2+γc02Wk-W~kG2≤Wk-W∗G2-γ2-γc02Wk-W~kG2.
Based on the Fejér monotonicity of the sequenceWk generated by PFPSM, the global convergence of PFPSM can be proved by similar arguments to those in [13], which is omitted.Theorem 13.
The sequenceWk generated PFPSM converges to some W∞∈W∗.


## 4. Numerical Results

In this section, we run a series of numerical experiments and apply the proposed PFPSM to solve some concrete examples of SPCP (6). We compare it with the splitting method for separable convex programming (denoted by HTY) in [30], the full Jacobian decomposition of the augmented Lagrangian method (denoted by FJDALM) in [13], and the fully parallel ADMM (denoted by PADMM) in [14]. Through these numerical experiments, we want to illustrate that (1) suitable proximal terms can enhance PFPSM’s numerical efficiency in practice, (2) larger values of the Glowinski relaxation factor η can often accelerate PFPSM’s convergence speed, and 3 compared with the other three ADMM-type methods the dynamically updated step size αk defined in (27) can accelerate the convergence of PFPSM. All codes were written by Matlab R2010a and conducted on a ThinkPad notebook with Pentium(R) Dual-Core CPU T4400@2.2 GHz, 4 GB of memory.Now, let us first introduce two shrinkage operators, which are used to solve the subproblems of the four tested methods. Letc>0 be a constant and E∈Rp×q be a given matrix, and let E=UΣV⊤ be the singular value decomposition (SVD) of E. Then, the solution of the optimization problem(57)minX∈Rp×qcX∗+12X-E2is given by the operator DcE∈Rp×q, which is defined by [31](58)DcE≔Udiag⁡ScΣV⊤.The solution of the optimization problem (59)minS∈Rp×qcS1+12S-E2is given by the operator ScE∈Rp×q, which is defined componentwisely by [32](60)ScEij≔sign⁡Eij·max⁡Eij-c,0,1≤i≤p,1≤j≤q,where sign· is the sign function.In the following, let us discuss how to solve the three subproblems of PFPSM. In fact, these three subproblems are simple enough to have closed-form solutions, which are listed as follows. For simplicity, we setRi=νβIpi=1,2,3 with ν≥0 in PFPSM. Based on the two shrinkage operators (58) and (60), the closed-form solutions of the three subproblems in the prediction step of PFPSM are summarized as follows:(i)
The explicit solutionX~1k is(61)X~1k=argminX1∈Rp×q⁡X1∗+β2X1+X2k+X3k-D-1βΛk2+νβ2X1-X1k2=D1/β1+ν11+νD+1βΛk-X2k-X3k+ν1+νX1k.(ii)
The explicit solutionX~2k is(62)X~2k=argminX2∈Rp×q⁡τX21+β2X1k+X2+X3k-D-1βΛk2+νβ2X2-X2k2=Sτ/β1+ν11+νD+1βΛk-X1k-X3k+ν1+νX2k.(iii)
The subproblem related toX~3 is (63)X~3k=argminX3∈Rp×q⁡12μPΩX32+β2X1k+X2k+X3-D-1βΛk2+νβ2X3-X3k2=argminX3∈Rp×q⁡12μPΩX32+β1+ν2X3-11+νD+1βΛk-X1k-X2k-ν1+νX3k2.Then, the explicit solution X~3k is given as follows:(64)X~3kij=μβ1+ν1+μβ1+νZ^ij,ifi,j∈Ω,Z^ij,otherwise,whereZ^=11+νD+1βΛk-X1k-X2k+ν1+νX3k.


### 4.1. Synthetic Data

In the experiment, following [5], we used the following synthetic data: the true low-rank matrix X1∗=L∗=V1V2⊤ is generated using normally distributed V1∈Rp×r and V2∈Rq×r. Hence, the rank of L∗ is r, and the low-rank ratio rr of L∗ is defined by rr=r/p. The nonzero entries of sparse matrix X2∗ or S∗ are generated uniformly in [-500,500], whose sparsity ratio is denoted by spr, defined by S∗0/pq, where S∗0 represents the number of nonzero entries of S∗. Then, we get the matrix D=L∗+S∗. Here, we consider that only a fraction of entries of D can be observed, and indices of the observed entries are collected in the index set Ω, which is determined randomly with a certain sample ratio (the ratio is denoted by sr=100Ω/pq%, where |Ω| is the cardinality of Ω). Then, we get the matrix PΩ(D), which is further corrupted by adding Gaussian noise with standard deviation σ=10-3. Finally, we get the true observed matrix PΩ(D)+σ*randp,q. Here, we test the following scenarios: p=q∈100,200, sr∈0.9,0.8,0.7, rr∈0.05,0.1, and spr∈0.05,0.1. The parameters τ and μ in (6) are set as τ=1/p,μ=p+8pσ/10, respectively.Throughout, we use the following stopping criterion:(65)maxX~1k-X1k1+X1k,X~2k-X2k1+X2k<10-5,and the initial point is set as X10,X20,X30,Λ0=0,0,0,0. For the penalty parameter β, we take β=0.06Ω/PΩD1 for all tested methods. The parameters in the tested methods are listed as follows:HTY:γ=2.01.FJDALM:γ=1.8.PADMM:Pi=2.4βIp,γ=1.1.PFPSM:Q=G,γ=1.5.Let us first test the sensitivity of PFPSM with respect to the proximal parameterν and the Glowinski relaxation factor η, and the numerical advantage of attaching the parameters ν and η to the iterative scheme of PFPSM is thus clear. We choose ν= 0 : 0.1 : 1 and fix γ=1.5, η=1.15 and p=q=200,rr=spr=0.05,sr=0.8. The maximum number of iterations is set to 500. The numerical results are plotted in Figure 1, in which “RLR” denotes the relative error of the recovered low-rank matrix X1k-X1∗/X1∗ and “RSP” denotes the relative error of the recovered sparse matrix X2k-X2∗/X2∗ when the stopping criterion (65) is satisfied. Figure 1 shows that small values of ν>0 can improve the accuracy of the solutions, while large values of ν>0 can accelerate the convergence speed of PFPSM. Then taking into consideration the numerical performance of other comparing methods, in the following we set ν=0.9. Now let us make the sensitivity analysis of the parameter η. We choose η= 0.86 : 0.02 : 1.15 and fix γ=1.5,ν=0.9 and p=q=100,rr=spr=0.05,sr=0.8. The maximum number of iterations is also set to 500. The numerical results are plotted in Figure 2. The two subplots in Figure 1 indicate that the number of iterations is descent and the elapsed CPU time has a descent tendency as the parameter η increases, which verify the numerical advantage of PFPSM with large values of η. Then, based on this observation, we set η=1.15 in the following experiments.Figure 1
Sensitivity test on the parameterν with p=q=100,rr=spr=0.05,sr=0.8.Figure 2
Sensitivity test on the parameterη with p=q=100,rr=spr=0.05,sr=0.8.Now, we compare PFPSM with PPSA, FJDALM, and PADMM to examine its numerical advantage from aspects of the number of iterations (denoted by “Iter.”), the elapsed CPU time in seconds (denoted by “Time”), the relative error of the recovered low-rank matrix, and the relative error of the recovered sparse matrix when the stopping criterion (65) is satisfied. Some numerical results are listed in Tables 1 and 2, which are the average results over 20 trials for each setting of parameters.Table 1
Numerical results for SPCP withp=q=100 and different rr,spr,sr.
r
r
s
p
r
s
r
Method
Iter.
Time
X
1
k
-
X
1
∗
/
X
1
∗
X
2
k
-
X
2
∗
/
X
2
∗
0.05
0.05
0.9
HTY
38
6.51
3.8250
e
-
004
4.5595
e
-
005
FJDALM
104
15.50
1.5313
e
-
003
5.5291
e
-
005
PADMM
74
6.46
4.3589
e
-
004
4.6012
e
-
005
PFPSM
46
4.35
3.8957
e
-
004
4.5757
e
-
005
0.05
0.05
0.8
HTY
45
7.07
4.7126
e
-
004
4.5788
e
-
005
FJDALM
121
16.41
2.5283
e
-
003
7.3386
e
-
005
PADMM
80
6.37
4.6868
e
-
004
4.1213
e
-
005
PFPSM
51
4.42
4.4273
e
-
004
4.5894
e
-
005
0.05
0.05
0.7
HTY
57
9.24
8.2791
e
-
004
4.9784
e
-
005
FJDALM
160
21.09
4.9488
e
-
003
1.4563
e
-
004
PADMM
88
7.43
5.1563
e
-
004
3.2281
e
-
005
PFPSM
58
5.50
5.2704
e
-
004
3.3050
e
-
005
0.05
0.1
0.9
HTY
47
7.43
1.0740
e
-
003
4.9652
e
-
005
FJDALM
171
25.36
3.6625
e
-
003
8.5186
e
-
005
PADMM
78
7.85
5.4060
e
-
004
3.9414
e
-
005
PFPSM
52
5.36
5.1620
e
-
004
3.9418
e
-
005
0.05
0.1
0.8
HTY
66
10.00
2.0001
e
-
003
6.2475
e
-
005
FJDALM
219
32.29
6.6939
e
-
003
1.5292
e
-
004
PADMM
91
7.82
1.5776
e
-
003
5.3509
e
-
005
PFPSM
63
6.93
1.4964
e
-
003
5.1892
e
-
005
0.05
0.1
0.7
HTY
107
16.97
7.5418
e
-
003
1.8982
e
-
004
FJDALM
299
39.63
1.5251
e
-
002
3.6022
e
-
004
PADMM
123
10.49
8.7675
e
-
003
2.2203
e
-
004
PFPSM
94
9.61
7.7442
e
-
003
1.9931
e
-
004
0.1
0.05
0.9
HTY
36
5.21
3.2385
e
-
004
3.1248
e
-
005
FJDALM
122
17.56
2.7233
e
-
003
1.2950
e
-
004
PADMM
73
5.91
3.5371
e
-
004
3.2098
e
-
005
PFPSM
44
3.79
2.8783
e
-
004
3.1123
e
-
005
0.1
0.05
0.8
HTY
53
7.68
2.4372
e
-
003
1.3577
e
-
004
FJDALM
178
22.78
7.1372
e
-
003
3.6697
e
-
004
PADMM
78
5.91
1.1473
e
-
003
6.8425
e
-
005
PFPSM
51
4.28
1.6505
e
-
003
9.3970
e
-
005
0.1
0.05
0.7
HTY
118
16.91
1.5141
e
-
002
7.8531
e
-
004
FJDALM
292
37.57
2.5434
e
-
002
1.2862
e
-
003
PADMM
140
11.63
1.5502
e
-
002
8.0841
e
-
004
PFPSM
116
9.32
1.4089
e
-
002
7.3816
e
-
004Table 2
Numerical results for SPCP withp=q=200 and different rr,spr,sr.
r
r
s
p
r
s
r
Method
Iter.
Time
X
1
k
-
X
1
∗
/
X
1
∗
X
2
k
-
X
2
∗
/
X
2
∗
0.05
0.05
0.9
HTY
38
5.94
3.8876
e
-
004
4.5610
e
-
005
FJDALM
104
18.69
1.5313
e
-
003
3.6239
e
-
005
PADMM
74
8.42
4.3589
e
-
004
4.6012
e
-
005
PFPSM
44
4.02
3.8750
e
-
004
4.5727
e
-
005
0.05
0.05
0.8
HTY
48
7.53
4.4271
e
-
004
4.5768
e
-
005
FJDALM
121
20.45
2.5283
e
-
003
4.2732
e
-
005
PADMM
80
8.51
4.6868
e
-
004
4.1213
e
-
005
PFPSM
48
3.98
4.3835
e
-
004
4.5828
e
-
005
0.05
0.05
0.7
HTY
59
9.11
1.3515
e
-
003
5.9535
e
-
005
FJDALM
160
27.73
4.9488
e
-
003
5.8452
e
-
005
PADMM
88
10.11
5.1563
e
-
004
3.2281
e
-
005
PFPSM
56
5.23
6.7517
e
-
004
3.4954
e
-
005
0.05
0.1
0.9
HTY
54
8.24
1.2604
e
-
003
5.1543
e
-
005
FJDALM
171
29.41
3.6625
e
-
003
4.3575
e
-
005
PADMM
78
8.94
5.4060
e
-
004
3.9414
e
-
005
PFPSM
49
4.77
6.9882
e
-
004
4.6485
e
-
005
0.05
0.1
0.8
HTY
76
11.54
2.3014
e
-
003
6.7335
e
-
005
FJDALM
219
38.30
6.6939
e
-
003
1.5292
e
-
004
PADMM
91
9.46
1.5776
e
-
003
5.3509
e
-
005
PFPSM
63
6.23
1.6119
e
-
003
5.3545
e
-
005
0.05
0.1
0.7
HTY
121
18.53
8.1084
e
-
003
2.0286
e
-
004
FJDALM
299
47.70
1.5251
e
-
002
4.7190
e
-
005
PADMM
123
13.31
8.7675
e
-
003
2.2203
e
-
004
PFPSM
99
8.51
7.3682
e
-
003
1.9037
e
-
004
0.1
0.05
0.9
HTY
37
5.58
4.7725
e
-
004
4.7778
e
-
005
FJDALM
122
39.53
2.7233
e
-
003
3.2379
e
-
005
PADMM
73
10.73
3.5371
e
-
004
3.2098
e
-
005
PFPSM
44
4.44
2.8783
e
-
004
3.1123
e
-
005
0.1
0.05
0.8
HTY
61
8.88
2.9153
e
-
003
1.6130
e
-
004
FJDALM
178
37.66
7.1372
e
-
003
2.3228
e
-
005
PADMM
78
7.45
1.1473
e
-
003
6.8425
e
-
005
PFPSM
51
5.54
1.6505
e
-
003
9.3970
e
-
005
0.1
0.05
0.7
HTY
133
18.55
1.5750
e
-
002
8.1458
e
-
004
FJDALM
292
38.59
2.5434
e
-
002
1.2862
e
-
003
PADMM
140
11.54
1.5502
e
-
002
8.0841
e
-
004
PFPSM
116
9.91
1.4089
e
-
002
7.3816
e
-
004The numerical results in Tables1 and 2 indicate the following: (1) From the “Iter.” criterion, the tested methods all successfully recover the given matrix D∗; (2) from the “Iter.” and “Time” criteria, our new method PFPSM performs better than the other three tested methods. It is almost the fastest for all the cases, and the reason may be that the suitable choice of the parameter η can accelerate the convergence of PFPSM; (3) compared to HTY, PFPSM is a fully parallel splitting method, which only recorded the maximum time taken by its three subproblems. Thus, we observe in Tables 1 and 2 that PFPSM takes less CPU time than HTY, though it needs to calculate a complicated correction step at each iteration or/and takes more iterations to terminate than HTY. In fact, considering the case p=q=100,rr=spr=0.05,sr=0.9 as an example, we find that the correction step of PFPSM takes about 7.26% of its whole CPU time, while the computation of the recovered sparse matrix and the simple correction step of HTY takes about 13.10% of its whole CPU time. In conclusion, the proposed PFPSM can be used as a surrogate of HTY, FJDALM, and PADMM.In Figures3 and 4, we display the evolution of X1k-X1∗/X1∗ and X2k-X2∗/X2∗ with respect to the iteration counter k, respectively, from which we can see that PFPSM performs better than FJDALM and PADMM and slightly worse than HTY in terms of the accuracy. The reason may be that PFPSM, FJDALM, and PADMM are fully parallel methods, while HTY is a partially parallel method. Therefore, HTY can fully utilize the latest iteration information to speed up its convergence. The relationship is just like that of the Jacobi iteration method and the Gauss-Seidel iteration method for linear system of equations. However, as illustrated in Table 1, the CPU time of PFPSM is almost always smaller than that of HTY, which is the main advantage of the parallel methods over the partially parallel or sequential methods.Figure 3
The relative error of the recovered low-rank matrix withp=q=100,rr=spr=0.05,sr=0.8.Figure 4
The relative error of the recovered sparse matrix withp=q=100,rr=spr=0.05,sr=0.8.


### 4.2. Background Extraction from Surveillance Video with Missing and Noisy Data

In this subsection, we are going to show the efficiency of the proposed method by computing the real-world problem: background extraction from surveillance video with missing and noisy data. In the experiment, following [6], we focus on a well-tested video taken in an airport, which consists of 50 grayscale frames with each frame having 144×176 pixels. We vectorize all frames of the video and get a matrix D∈R25344×50 where each column represents a frame, which is further contaminated with additive Gaussian noise of mean zero and standard deviation σ=10-3. Similar to the last experiment, only a fraction of entries of D can be observed, and indices of the observed entries are also collected in the index set Ω. In the experiment, we set τ=1/p,δ=0.01 and the initial iteration is set to be zero. Other parameters remain the same as those used in the last experiment for all the tested algorithms. The numerical results are given in Table 3, where Obj. represents the value of the objective function of (1) when the tested methods terminate.Table 3
Numerical results for background extraction problem.
s
r
Method
Iter.
Time
rank(X1)
X20
Obj.
0.9
HTY
28
23.0
9
3.25
e
+
005
2.0481
e
+
003
PADMM
45
27.3
7
2.80
e
+
005
1.7668
e
+
003
PFPSM
38
25.5
5
1.22
e
+
005
770.4871
0.8
HTY
28
23.1
7
3.09
e
+
005
1.9492
e
+
003
PADMM
50
28.6
7
2.83
e
+
005
1.7853
e
+
003
PFPSM
41
28.0
4
1.14
e
+
005
723.1297
0.7
HTY
30
23.4
6
3.06
e
+
005
1.9291
e
+
003
PADMM
54
34.5
6
2.81
e
+
005
1.7692
e
+
003
PFPSM
46
34.4
3
1.09
e
+
005
687.3679Table3 shows that although PFPSM implements more iterations than HTY for all test instances, the performance of PFPSM is much better than the other two methods in recovering low-rank and sparse components. In Figure 5, we exhibit a frame in the tested video and the separation results with 30% random missing data. Clearly, all the tested methods can extract the disturbed image quite well.Figure 5
The 10th frame of the clean video and the corresponding corrupted frame with 30% random missing data (a); the extracted background and foreground by HTY (b), PADMM (c), and PFPSM (d).
(a)
(b)
(c)
(d)


## 4.1. Synthetic Data

In the experiment, following [5], we used the following synthetic data: the true low-rank matrix X1∗=L∗=V1V2⊤ is generated using normally distributed V1∈Rp×r and V2∈Rq×r. Hence, the rank of L∗ is r, and the low-rank ratio rr of L∗ is defined by rr=r/p. The nonzero entries of sparse matrix X2∗ or S∗ are generated uniformly in [-500,500], whose sparsity ratio is denoted by spr, defined by S∗0/pq, where S∗0 represents the number of nonzero entries of S∗. Then, we get the matrix D=L∗+S∗. Here, we consider that only a fraction of entries of D can be observed, and indices of the observed entries are collected in the index set Ω, which is determined randomly with a certain sample ratio (the ratio is denoted by sr=100Ω/pq%, where |Ω| is the cardinality of Ω). Then, we get the matrix PΩ(D), which is further corrupted by adding Gaussian noise with standard deviation σ=10-3. Finally, we get the true observed matrix PΩ(D)+σ*randp,q. Here, we test the following scenarios: p=q∈100,200, sr∈0.9,0.8,0.7, rr∈0.05,0.1, and spr∈0.05,0.1. The parameters τ and μ in (6) are set as τ=1/p,μ=p+8pσ/10, respectively.Throughout, we use the following stopping criterion:(65)maxX~1k-X1k1+X1k,X~2k-X2k1+X2k<10-5,and the initial point is set as X10,X20,X30,Λ0=0,0,0,0. For the penalty parameter β, we take β=0.06Ω/PΩD1 for all tested methods. The parameters in the tested methods are listed as follows:HTY:γ=2.01.FJDALM:γ=1.8.PADMM:Pi=2.4βIp,γ=1.1.PFPSM:Q=G,γ=1.5.Let us first test the sensitivity of PFPSM with respect to the proximal parameterν and the Glowinski relaxation factor η, and the numerical advantage of attaching the parameters ν and η to the iterative scheme of PFPSM is thus clear. We choose ν= 0 : 0.1 : 1 and fix γ=1.5, η=1.15 and p=q=200,rr=spr=0.05,sr=0.8. The maximum number of iterations is set to 500. The numerical results are plotted in Figure 1, in which “RLR” denotes the relative error of the recovered low-rank matrix X1k-X1∗/X1∗ and “RSP” denotes the relative error of the recovered sparse matrix X2k-X2∗/X2∗ when the stopping criterion (65) is satisfied. Figure 1 shows that small values of ν>0 can improve the accuracy of the solutions, while large values of ν>0 can accelerate the convergence speed of PFPSM. Then taking into consideration the numerical performance of other comparing methods, in the following we set ν=0.9. Now let us make the sensitivity analysis of the parameter η. We choose η= 0.86 : 0.02 : 1.15 and fix γ=1.5,ν=0.9 and p=q=100,rr=spr=0.05,sr=0.8. The maximum number of iterations is also set to 500. The numerical results are plotted in Figure 2. The two subplots in Figure 1 indicate that the number of iterations is descent and the elapsed CPU time has a descent tendency as the parameter η increases, which verify the numerical advantage of PFPSM with large values of η. Then, based on this observation, we set η=1.15 in the following experiments.Figure 1
Sensitivity test on the parameterν with p=q=100,rr=spr=0.05,sr=0.8.Figure 2
Sensitivity test on the parameterη with p=q=100,rr=spr=0.05,sr=0.8.Now, we compare PFPSM with PPSA, FJDALM, and PADMM to examine its numerical advantage from aspects of the number of iterations (denoted by “Iter.”), the elapsed CPU time in seconds (denoted by “Time”), the relative error of the recovered low-rank matrix, and the relative error of the recovered sparse matrix when the stopping criterion (65) is satisfied. Some numerical results are listed in Tables 1 and 2, which are the average results over 20 trials for each setting of parameters.Table 1
Numerical results for SPCP withp=q=100 and different rr,spr,sr.
r
r
s
p
r
s
r
Method
Iter.
Time
X
1
k
-
X
1
∗
/
X
1
∗
X
2
k
-
X
2
∗
/
X
2
∗
0.05
0.05
0.9
HTY
38
6.51
3.8250
e
-
004
4.5595
e
-
005
FJDALM
104
15.50
1.5313
e
-
003
5.5291
e
-
005
PADMM
74
6.46
4.3589
e
-
004
4.6012
e
-
005
PFPSM
46
4.35
3.8957
e
-
004
4.5757
e
-
005
0.05
0.05
0.8
HTY
45
7.07
4.7126
e
-
004
4.5788
e
-
005
FJDALM
121
16.41
2.5283
e
-
003
7.3386
e
-
005
PADMM
80
6.37
4.6868
e
-
004
4.1213
e
-
005
PFPSM
51
4.42
4.4273
e
-
004
4.5894
e
-
005
0.05
0.05
0.7
HTY
57
9.24
8.2791
e
-
004
4.9784
e
-
005
FJDALM
160
21.09
4.9488
e
-
003
1.4563
e
-
004
PADMM
88
7.43
5.1563
e
-
004
3.2281
e
-
005
PFPSM
58
5.50
5.2704
e
-
004
3.3050
e
-
005
0.05
0.1
0.9
HTY
47
7.43
1.0740
e
-
003
4.9652
e
-
005
FJDALM
171
25.36
3.6625
e
-
003
8.5186
e
-
005
PADMM
78
7.85
5.4060
e
-
004
3.9414
e
-
005
PFPSM
52
5.36
5.1620
e
-
004
3.9418
e
-
005
0.05
0.1
0.8
HTY
66
10.00
2.0001
e
-
003
6.2475
e
-
005
FJDALM
219
32.29
6.6939
e
-
003
1.5292
e
-
004
PADMM
91
7.82
1.5776
e
-
003
5.3509
e
-
005
PFPSM
63
6.93
1.4964
e
-
003
5.1892
e
-
005
0.05
0.1
0.7
HTY
107
16.97
7.5418
e
-
003
1.8982
e
-
004
FJDALM
299
39.63
1.5251
e
-
002
3.6022
e
-
004
PADMM
123
10.49
8.7675
e
-
003
2.2203
e
-
004
PFPSM
94
9.61
7.7442
e
-
003
1.9931
e
-
004
0.1
0.05
0.9
HTY
36
5.21
3.2385
e
-
004
3.1248
e
-
005
FJDALM
122
17.56
2.7233
e
-
003
1.2950
e
-
004
PADMM
73
5.91
3.5371
e
-
004
3.2098
e
-
005
PFPSM
44
3.79
2.8783
e
-
004
3.1123
e
-
005
0.1
0.05
0.8
HTY
53
7.68
2.4372
e
-
003
1.3577
e
-
004
FJDALM
178
22.78
7.1372
e
-
003
3.6697
e
-
004
PADMM
78
5.91
1.1473
e
-
003
6.8425
e
-
005
PFPSM
51
4.28
1.6505
e
-
003
9.3970
e
-
005
0.1
0.05
0.7
HTY
118
16.91
1.5141
e
-
002
7.8531
e
-
004
FJDALM
292
37.57
2.5434
e
-
002
1.2862
e
-
003
PADMM
140
11.63
1.5502
e
-
002
8.0841
e
-
004
PFPSM
116
9.32
1.4089
e
-
002
7.3816
e
-
004Table 2
Numerical results for SPCP withp=q=200 and different rr,spr,sr.
r
r
s
p
r
s
r
Method
Iter.
Time
X
1
k
-
X
1
∗
/
X
1
∗
X
2
k
-
X
2
∗
/
X
2
∗
0.05
0.05
0.9
HTY
38
5.94
3.8876
e
-
004
4.5610
e
-
005
FJDALM
104
18.69
1.5313
e
-
003
3.6239
e
-
005
PADMM
74
8.42
4.3589
e
-
004
4.6012
e
-
005
PFPSM
44
4.02
3.8750
e
-
004
4.5727
e
-
005
0.05
0.05
0.8
HTY
48
7.53
4.4271
e
-
004
4.5768
e
-
005
FJDALM
121
20.45
2.5283
e
-
003
4.2732
e
-
005
PADMM
80
8.51
4.6868
e
-
004
4.1213
e
-
005
PFPSM
48
3.98
4.3835
e
-
004
4.5828
e
-
005
0.05
0.05
0.7
HTY
59
9.11
1.3515
e
-
003
5.9535
e
-
005
FJDALM
160
27.73
4.9488
e
-
003
5.8452
e
-
005
PADMM
88
10.11
5.1563
e
-
004
3.2281
e
-
005
PFPSM
56
5.23
6.7517
e
-
004
3.4954
e
-
005
0.05
0.1
0.9
HTY
54
8.24
1.2604
e
-
003
5.1543
e
-
005
FJDALM
171
29.41
3.6625
e
-
003
4.3575
e
-
005
PADMM
78
8.94
5.4060
e
-
004
3.9414
e
-
005
PFPSM
49
4.77
6.9882
e
-
004
4.6485
e
-
005
0.05
0.1
0.8
HTY
76
11.54
2.3014
e
-
003
6.7335
e
-
005
FJDALM
219
38.30
6.6939
e
-
003
1.5292
e
-
004
PADMM
91
9.46
1.5776
e
-
003
5.3509
e
-
005
PFPSM
63
6.23
1.6119
e
-
003
5.3545
e
-
005
0.05
0.1
0.7
HTY
121
18.53
8.1084
e
-
003
2.0286
e
-
004
FJDALM
299
47.70
1.5251
e
-
002
4.7190
e
-
005
PADMM
123
13.31
8.7675
e
-
003
2.2203
e
-
004
PFPSM
99
8.51
7.3682
e
-
003
1.9037
e
-
004
0.1
0.05
0.9
HTY
37
5.58
4.7725
e
-
004
4.7778
e
-
005
FJDALM
122
39.53
2.7233
e
-
003
3.2379
e
-
005
PADMM
73
10.73
3.5371
e
-
004
3.2098
e
-
005
PFPSM
44
4.44
2.8783
e
-
004
3.1123
e
-
005
0.1
0.05
0.8
HTY
61
8.88
2.9153
e
-
003
1.6130
e
-
004
FJDALM
178
37.66
7.1372
e
-
003
2.3228
e
-
005
PADMM
78
7.45
1.1473
e
-
003
6.8425
e
-
005
PFPSM
51
5.54
1.6505
e
-
003
9.3970
e
-
005
0.1
0.05
0.7
HTY
133
18.55
1.5750
e
-
002
8.1458
e
-
004
FJDALM
292
38.59
2.5434
e
-
002
1.2862
e
-
003
PADMM
140
11.54
1.5502
e
-
002
8.0841
e
-
004
PFPSM
116
9.91
1.4089
e
-
002
7.3816
e
-
004The numerical results in Tables1 and 2 indicate the following: (1) From the “Iter.” criterion, the tested methods all successfully recover the given matrix D∗; (2) from the “Iter.” and “Time” criteria, our new method PFPSM performs better than the other three tested methods. It is almost the fastest for all the cases, and the reason may be that the suitable choice of the parameter η can accelerate the convergence of PFPSM; (3) compared to HTY, PFPSM is a fully parallel splitting method, which only recorded the maximum time taken by its three subproblems. Thus, we observe in Tables 1 and 2 that PFPSM takes less CPU time than HTY, though it needs to calculate a complicated correction step at each iteration or/and takes more iterations to terminate than HTY. In fact, considering the case p=q=100,rr=spr=0.05,sr=0.9 as an example, we find that the correction step of PFPSM takes about 7.26% of its whole CPU time, while the computation of the recovered sparse matrix and the simple correction step of HTY takes about 13.10% of its whole CPU time. In conclusion, the proposed PFPSM can be used as a surrogate of HTY, FJDALM, and PADMM.In Figures3 and 4, we display the evolution of X1k-X1∗/X1∗ and X2k-X2∗/X2∗ with respect to the iteration counter k, respectively, from which we can see that PFPSM performs better than FJDALM and PADMM and slightly worse than HTY in terms of the accuracy. The reason may be that PFPSM, FJDALM, and PADMM are fully parallel methods, while HTY is a partially parallel method. Therefore, HTY can fully utilize the latest iteration information to speed up its convergence. The relationship is just like that of the Jacobi iteration method and the Gauss-Seidel iteration method for linear system of equations. However, as illustrated in Table 1, the CPU time of PFPSM is almost always smaller than that of HTY, which is the main advantage of the parallel methods over the partially parallel or sequential methods.Figure 3
The relative error of the recovered low-rank matrix withp=q=100,rr=spr=0.05,sr=0.8.Figure 4
The relative error of the recovered sparse matrix withp=q=100,rr=spr=0.05,sr=0.8.


## 4.2. Background Extraction from Surveillance Video with Missing and Noisy Data

In this subsection, we are going to show the efficiency of the proposed method by computing the real-world problem: background extraction from surveillance video with missing and noisy data. In the experiment, following [6], we focus on a well-tested video taken in an airport, which consists of 50 grayscale frames with each frame having 144×176 pixels. We vectorize all frames of the video and get a matrix D∈R25344×50 where each column represents a frame, which is further contaminated with additive Gaussian noise of mean zero and standard deviation σ=10-3. Similar to the last experiment, only a fraction of entries of D can be observed, and indices of the observed entries are also collected in the index set Ω. In the experiment, we set τ=1/p,δ=0.01 and the initial iteration is set to be zero. Other parameters remain the same as those used in the last experiment for all the tested algorithms. The numerical results are given in Table 3, where Obj. represents the value of the objective function of (1) when the tested methods terminate.Table 3
Numerical results for background extraction problem.
s
r
Method
Iter.
Time
rank(X1)
X20
Obj.
0.9
HTY
28
23.0
9
3.25
e
+
005
2.0481
e
+
003
PADMM
45
27.3
7
2.80
e
+
005
1.7668
e
+
003
PFPSM
38
25.5
5
1.22
e
+
005
770.4871
0.8
HTY
28
23.1
7
3.09
e
+
005
1.9492
e
+
003
PADMM
50
28.6
7
2.83
e
+
005
1.7853
e
+
003
PFPSM
41
28.0
4
1.14
e
+
005
723.1297
0.7
HTY
30
23.4
6
3.06
e
+
005
1.9291
e
+
003
PADMM
54
34.5
6
2.81
e
+
005
1.7692
e
+
003
PFPSM
46
34.4
3
1.09
e
+
005
687.3679Table3 shows that although PFPSM implements more iterations than HTY for all test instances, the performance of PFPSM is much better than the other two methods in recovering low-rank and sparse components. In Figure 5, we exhibit a frame in the tested video and the separation results with 30% random missing data. Clearly, all the tested methods can extract the disturbed image quite well.Figure 5
The 10th frame of the clean video and the corresponding corrupted frame with 30% random missing data (a); the extracted background and foreground by HTY (b), PADMM (c), and PFPSM (d).
(a)
(b)
(c)
(d)


## 5. Conclusion

In this paper, we have proposed a new fully parallel splitting method for the stable principal component pursuit problem, which inherits all merits of ADMM and parallel method. Compared to other similar methods in the literature, the new method attaches a Glowinski relaxation factorη to the updating formula for the Lagrange multiplier, and large values of η are often beneficial to its convergence in general. Preliminary numerical results indicate that the new method is quite efficient.To ensure the convergence of the proposed iteration method, the scope of the involved parameterη is 3/2,2/3, which is a proper subset of the interval 0,1+5/2, the scope of the Glowinski relaxation factor η of the sequential ADMM-type methods for two-block separable convex programming [25, 33]. Therefore, it is worth extending the scope of the involved parameter η and studying its optimal choice.



---

*Source: 9674528-2017-10-25.xml*