# An Implementable First-Order Primal-Dual Algorithm for Structured Convex Optimization

Authors: Feng Ma; Mingfang Ni; Lei Zhu; Zhanke Yu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Many application problems of practical interest can be posed as structured convex optimization models. In this paper, we study a new first-order primaldual algorithm. The method can be easily implementable, provided that the resolvent operators of the component objective functions are simple to evaluate. We show that the proposed method can be interpreted as a proximal point algorithm with a customized metric proximal parameter. Convergence property is established under the analytic contraction framework. Finally, we verify the efficiency of the algorithm by solving the stable principal component pursuit problem.


---

## Body

## 1. Introduction

In this paper, we consider the following separable optimization problem:(1)min⁡x∈𝒳,y∈𝒴⁡f(x)+g(y)s.t.Ax+By=b,
where A∈ℜm×n and B∈ℜm×p are given matrices, b∈ℜm is a given vector, and 𝒳⊂ℜn and 𝒴⊂ℜp are nonempty closed convex sets. f:ℜn→ℜ and g:ℜp→ℜ are convex functions (not necessarily smooth). Throughout this paper, we assume that the solution set of (1) is nonempty. Problem of this type arises in many applications, ranging from machine learning to compressed sensing. We refer to, for example, [1–9], for a few examples of applications.To solve (1), one can use the classical augmented Lagrangian method (ALM). Starting from any initial iterate (x0,y0,λ0), ALM iterates via the following procedure:(2a)(xk+1,yk+1)=argmin⁡x,y⁡ℒ𝒜(x,y;λk),(2b)λk+1=λk-β(Axk+1+Byk+1-b),where the augmented Lagrangian function ℒ𝒜(x,y;λ) associated with the problem (1) is given by
(3)ℒ𝒜(x,y,λ)=f(x)+g(y)-λT(Ax+By-b)+β2∥Ax+By-b∥2,
and λ∈ℜm is the Lagrangian multiplier vector associated with the linear constraint and β>0 is a penalty parameter for the violation of the linear constraint. ALM enjoys very nice convergence and has been shown to be equivalent to a proximal point algorithm applied to the dual of (1) [10]. A noticeable feature of ALM is that it treats (1) as a generic minimization problem and ignores completely the nice separable structure emerging in the objective function. The minimizations of the two functions f and g in (2a) are strongly coupled because of the quadratic term (β/2)∥Ax+By-b∥2. Hence, the implementation of ALM ((2a) and (2b)) can be computationally challenging. To utilize the separable structure of the problem, the well-known alternating direction method (ADM) essentially splits the ALM subproblem into two subproblems with respect to x and y in Gauss-Seidel manner. More specifically, at each iteration, ADM takes the following form:(4a)xk+1=argmin⁡x∈𝒳⁡ℒ𝒜(x,yk;λk),(4b)yk+1=argmin⁡y∈𝒴⁡ℒ𝒜(xk+1,y;λk),(4c)λk+1=λk-β(Axk+1+Byk+1-b).ADM can be interpreted as Douglas-Rachford splitting method applied to the dual problem [11] and proximal point method [12]. We refer to, for example, [13–17] for recent study of ADM and its variants. Comparing to ALM, ADM minimizes (3) with respect to x and y alternatingly at each iteration, rather than with respect to both x and y simultaneously. This splitting procedure makes it possible to exploit the special separable structure of the objective functions. Thus, ADM appears to be a natural fit for solving very large scale distributed machine learning and big-data related optimization problems.WhenA and B in (1) are both identity matrices, ADM can be very efficient. The first two subproblems ((4a) and (4b)) correspond to evaluating the resolvent operators of component functions f and g, respectively. Here, theproximal operator of the function θ:ℜn→ℜ is defined by
(5)u:=proxξθ(a)=argmin⁡⁡{θ(u)+12ξ∥u-a∥2∣u∈ℜn},
where a∈ℜn and ξ>0. In most popular applications of sparse optimization, the proximal operator can be computed exactly and efficiently (e.g., θ(x)=∥x∥1:=∑i=1n‍|xi|orθ(x)=∥x∥2:=∑i=1n‍|xi|2). While when A (resp., B) is not an identity matrix, the resulting ADM subproblems may not be easily solvable, since they involve inverting of A (resp., B), and there is no efficient way of doing so directly. This difficulty could result in inefficiency of the ADM greatly. As a result, first-order algorithms that preserves both the alternating computation feature of ADM and the advantages of not involving the inverse of A and B are highly desirable.In [18, 19], an alternating proximal gradient method (APGM) for solving (1) was given. APGM is based on the framework of ADM, which solves the first two subproblems ((4a) and (4b)) inexactly by taking one proximal gradient step. More specifically, APGM generates the new iterate (xk+1,yk+1,λk+1) via the following procedure:(6a)xk+1=argmin⁡x∈𝒳⁡f(x)+r2∥x-xk-1rAT(λk-β(Axk+Byk-b))∥2,(6b)yk+1=argmin⁡y∈𝒴⁡g(y)+s2∥y-yk-1sBT(λk-β(Axk+1+Byk-b))∥2,(6c)λk+1=λk-β(Axk+1+Byk+1-b),where β>0, r>β∥ATA∥, and s>β∥BTB∥.In this paper, we also focus on the specific scenario of (1), where both A and B are not identity matrices. We propose a new first-order primal-dual algorithm for (1). We show that our algorithm can be viewed as a customized proximal point method with a special proximal regularization parameter, and it is within the framework of the contraction type methods. The proposed algorithm has three main advantages: first, it can be easily implementable, the main computational effort of each iteration is to evaluate the proximal operators of the component objective functions; second, the involved subproblems are solved consecutively in the ADM manner, which makes the algorithm amenable to distributed optimization. Finally, it only uses the first-order information and does not require any matrix inversion or solving linear systems. Hence, the method is well suited for solving large scale distributed optimization problems.The paper is organized as follows. In Section2, we review some preliminaries. In Sections 3 and 4, we present the new method and analyze its convergence. In Section 5 we conduct numerical experiments to compare the proposed method with APGM for solving the stable principal component pursuit problem. We conclude the paper in Section 6.


## 2. Preliminaries

### 2.1. Variational Characterization of (1)

In this section, we reformulate (1) as a variational form, which is useful for succedent algorithmic illustration and convergence analysis.The Lagrangian function associated with (1) is
(7)ℒ(x,y,λ)=f(x)+g(y)-λT(Ax+By-b),
where λ∈ℜm is a Lagrangian multiplier. According to the previous convex assumption of (1), finding optimal solutions of (1) and its dual form is equivalent to finding a saddle point of ℒ. More precisely, let (x*,y*,λ*) be a saddle point of ℒ. We have
(8)ℒλ∈ℜm(x*,y*,λ)≤ℒ(x*,y*,λ*)≤ℒx∈𝒳,y∈𝒴(x,y,λ*).
Then we can directly read off the optimality conditions with variational characterization. More specifically, w*=(x*,y*,λ*)∈Ω is a saddle point of ℒ if and only if it satisfies the following mixed variational inequality (VI):(9a)h(u)-h(u*)+(w′-w*)TF(w*)≥0,∀w′∈Ω,
with
(9b)u:=(xy),w:=(xyλ),h(u):=f(x)+g(y),F(w):=(-ATλ-BTλAx+By-b),(9c)Ω=𝒳×𝒴×ℜm. Hence, a solution of ((9a)–(9c)) yields a solution of (1).Note that the mappingF(·) is said to bemonotonewith respect to Ω if
(10)(F(w)-F(w′))T(w-w′)≥0,∀w,w′∈Ω.
Consequently, it can be easily verified that F(w) is monotone. Under the aforementioned nonempty assumption on the solution set of (1), the solution set of ((9a)–(9c)), denoted by 𝒲*, is also nonempty.


### 2.2. Proximal Point Algorithmic Framework

In this subsection, we review the classical proximal point algorithm (PPA) for solving the VI ((9a)–(9c)).PPA, which was proposed by Martinet [20] and further studied by Rockafellar [10], plays a fundamental rule in optimization. For given iterate wk∈Ω, PPA generates the new iterative wk+1 via the following procedure:
(11)h(u′)-h(uk+1)+(w′-wk+1)T×(F(wk+1)+G(wk+1-wk))≥0,∀w′∈Ω,
where G is a positive definite matrix, playing the role of proximal regularization parameter. A simple choice of G is that G=β·I where β>0 and I is the identity matrix, regularizing the proximal terms wk+1-wk in the uniform way. We refer the reader to example [21–24] for some special choices of G in different scenarios.


## 2.1. Variational Characterization of (1)

In this section, we reformulate (1) as a variational form, which is useful for succedent algorithmic illustration and convergence analysis.The Lagrangian function associated with (1) is
(7)ℒ(x,y,λ)=f(x)+g(y)-λT(Ax+By-b),
where λ∈ℜm is a Lagrangian multiplier. According to the previous convex assumption of (1), finding optimal solutions of (1) and its dual form is equivalent to finding a saddle point of ℒ. More precisely, let (x*,y*,λ*) be a saddle point of ℒ. We have
(8)ℒλ∈ℜm(x*,y*,λ)≤ℒ(x*,y*,λ*)≤ℒx∈𝒳,y∈𝒴(x,y,λ*).
Then we can directly read off the optimality conditions with variational characterization. More specifically, w*=(x*,y*,λ*)∈Ω is a saddle point of ℒ if and only if it satisfies the following mixed variational inequality (VI):(9a)h(u)-h(u*)+(w′-w*)TF(w*)≥0,∀w′∈Ω,
with
(9b)u:=(xy),w:=(xyλ),h(u):=f(x)+g(y),F(w):=(-ATλ-BTλAx+By-b),(9c)Ω=𝒳×𝒴×ℜm. Hence, a solution of ((9a)–(9c)) yields a solution of (1).Note that the mappingF(·) is said to bemonotonewith respect to Ω if
(10)(F(w)-F(w′))T(w-w′)≥0,∀w,w′∈Ω.
Consequently, it can be easily verified that F(w) is monotone. Under the aforementioned nonempty assumption on the solution set of (1), the solution set of ((9a)–(9c)), denoted by 𝒲*, is also nonempty.


## 2.2. Proximal Point Algorithmic Framework

In this subsection, we review the classical proximal point algorithm (PPA) for solving the VI ((9a)–(9c)).PPA, which was proposed by Martinet [20] and further studied by Rockafellar [10], plays a fundamental rule in optimization. For given iterate wk∈Ω, PPA generates the new iterative wk+1 via the following procedure:
(11)h(u′)-h(uk+1)+(w′-wk+1)T×(F(wk+1)+G(wk+1-wk))≥0,∀w′∈Ω,
where G is a positive definite matrix, playing the role of proximal regularization parameter. A simple choice of G is that G=β·I where β>0 and I is the identity matrix, regularizing the proximal terms wk+1-wk in the uniform way. We refer the reader to example [21–24] for some special choices of G in different scenarios.


## 3. The Main Algorithm

### 3.1. Assumption

Before we present our new algorithm, we need to make the following assumption:Assumption 1.
For any givea∈ℜn+p and ξ>0, theproximal operator of h(u) (see (5)) has a closed-form solution or it can be efficiently solved up to a high precision.Whenever the assumption holds, we say that the proximal operator ofh is “easy” to evaluate. Note that h is separable across two variables, that is, h=f+g; according to the definition (5), we have
(12)proxξh(b,c)=(proxξf(b),proxξg(c)),
where b∈ℜn and c∈ℜp. Hence, under the assumption, the proximal operators of the component objective functions f and g are also “easy” to evaluate.


### 3.2. Motivation

The motivation for our algorithm is directly related to the linearized augmented Lagrangian method proposed in [25] and the customized proximal point algorithm proposed in [21] for convex problems with linear constraints.In order to obtain a closed-form solution ofx, we first add a proximal regularization parameter G=(rI0AT000A0(1/β)I) to the variational characterization of the problem ((9a)–(9c)). Then, like the algorithm in [21], we get the following proximal point algorithm:(13a)xk+1=argmin⁡x∈𝒳⁡f(x)+r2∥x-xk-1rATλk∥2,(13b)yk+1=argmin⁡y∈𝒴⁡g(y)+β2∥A(2xk+1-xk)+By-b-λkβ∥2,(13c)λk+1=λk-β(Axk+1+Byk+1-b). However, (13b) is still not implementable. Inspired in [25], in order to alleviate the computation required by y-subproblems, we try to linearize the quadratic term in (13b) by
(14)β2∥A(2xk+1-xk)+By-b-λkβ∥2≈β2∥A(2xk+1-xk)+Byk-b-λkβ∥2+β(y-yk)Tgk+s2∥y-yk∥,
where s>0 is a proximal parameter, and gk=BT(A(2xk+1-xk)+Byk-b-λk/β) is the gradient of (1/2)∥A(2xk+1-xk)+By-b-λk/β∥2 at yk. With a simple manipulation, we get the following approximation to (13b):
(15)yk+1=argmin⁡y∈𝒴⁡g(y)+s2∥y-yk-1sBT(λk-β×(A(2xk+1-xk)+Byk-b))1s∥2.
In this case, the closed-form solutions of the resulting subproblems can be easily obtained.


### 3.3. Description of the Algorithm

In this section, we formally present Algorithm1.Algorithm 1
The proposed algorithm for (1)
Step  0. Input (x0,y0,λ0)∈𝒳×𝒴.
Step  1. Set
xk+1=argmin⁡x∈𝒳f(x)+r2∥x-xk-1rATλk∥2,                 (*)
yk+1=argmin⁡y∈𝒴g(y)+s2∥y-yk-1sBT(λk-β(A(2xk+1-xk)+Byk-b))∥2,
(**)
λk+1=λk-β(A(2xk+1-xk)+Byk+1-b),                  (***)
whereβ>0, and r>β∥ATA∥,s>β∥BTB∥.
Step  2. If a termination criterion is not met, go to Step  1.Remark 2.
From the algorithm, we can see that the minimizations ((*) and (**)) each requires the evaluation of the proximal operators for the component objective functions f and g. It is clear that the implementation of the proposed method is simple under our assumption.


## 3.1. Assumption

Before we present our new algorithm, we need to make the following assumption:Assumption 1.
For any givea∈ℜn+p and ξ>0, theproximal operator of h(u) (see (5)) has a closed-form solution or it can be efficiently solved up to a high precision.Whenever the assumption holds, we say that the proximal operator ofh is “easy” to evaluate. Note that h is separable across two variables, that is, h=f+g; according to the definition (5), we have
(12)proxξh(b,c)=(proxξf(b),proxξg(c)),
where b∈ℜn and c∈ℜp. Hence, under the assumption, the proximal operators of the component objective functions f and g are also “easy” to evaluate.


## 3.2. Motivation

The motivation for our algorithm is directly related to the linearized augmented Lagrangian method proposed in [25] and the customized proximal point algorithm proposed in [21] for convex problems with linear constraints.In order to obtain a closed-form solution ofx, we first add a proximal regularization parameter G=(rI0AT000A0(1/β)I) to the variational characterization of the problem ((9a)–(9c)). Then, like the algorithm in [21], we get the following proximal point algorithm:(13a)xk+1=argmin⁡x∈𝒳⁡f(x)+r2∥x-xk-1rATλk∥2,(13b)yk+1=argmin⁡y∈𝒴⁡g(y)+β2∥A(2xk+1-xk)+By-b-λkβ∥2,(13c)λk+1=λk-β(Axk+1+Byk+1-b). However, (13b) is still not implementable. Inspired in [25], in order to alleviate the computation required by y-subproblems, we try to linearize the quadratic term in (13b) by
(14)β2∥A(2xk+1-xk)+By-b-λkβ∥2≈β2∥A(2xk+1-xk)+Byk-b-λkβ∥2+β(y-yk)Tgk+s2∥y-yk∥,
where s>0 is a proximal parameter, and gk=BT(A(2xk+1-xk)+Byk-b-λk/β) is the gradient of (1/2)∥A(2xk+1-xk)+By-b-λk/β∥2 at yk. With a simple manipulation, we get the following approximation to (13b):
(15)yk+1=argmin⁡y∈𝒴⁡g(y)+s2∥y-yk-1sBT(λk-β×(A(2xk+1-xk)+Byk-b))1s∥2.
In this case, the closed-form solutions of the resulting subproblems can be easily obtained.


## 3.3. Description of the Algorithm

In this section, we formally present Algorithm1.Algorithm 1
The proposed algorithm for (1)
Step  0. Input (x0,y0,λ0)∈𝒳×𝒴.
Step  1. Set
xk+1=argmin⁡x∈𝒳f(x)+r2∥x-xk-1rATλk∥2,                 (*)
yk+1=argmin⁡y∈𝒴g(y)+s2∥y-yk-1sBT(λk-β(A(2xk+1-xk)+Byk-b))∥2,
(**)
λk+1=λk-β(A(2xk+1-xk)+Byk+1-b),                  (***)
whereβ>0, and r>β∥ATA∥,s>β∥BTB∥.
Step  2. If a termination criterion is not met, go to Step  1.Remark 2.
From the algorithm, we can see that the minimizations ((*) and (**)) each requires the evaluation of the proximal operators for the component objective functions f and g. It is clear that the implementation of the proposed method is simple under our assumption.


## 4. Global Convergence

In this section, we show that the proposed method is in some sense equivalent to the proximal point algorithm with a special proximal regularization parameter. Then its convergence can be easily established under the analytic framework of contraction type methods.Lemma 3.
Letwk+1=(xk+1,yk+1,λk+1)T be generated by ((*)–(***)) from a given wk=(xk,yk,λk)T and H=(rI0AT0sI-βBTB0A0(1/β)I). Then
(16)h(u′)-h(uk+1)+(w′-wk+1)T×(F(wk+1)+H(wk+1-wk))≥0,2222222222222222222122222∀w′∈Ω.Proof.
First, by deriving the optimality condition for the subproblem (*), we have
(17)f(x)-f(xk+1)+(x′-xk+1)T(r(xk+1-xk)-ATλk)≥0,222222222222222222222222222222222222222222∀x′∈𝒳.
It can be further rewritten as
(18)f(x)-f(xk+1)+(x′-xk+1)T×(-ATλk+1+r(xk+1-xk)+AT(λk+1-λk))≥0,22222222222222222222222222222222222222∀x′∈𝒳.
Similarly, the optimality condition for the subproblem (**) is
(19)g(y)-g(yk+1)+(y′-yk+1)T(s(yk+1-yk)-BT(λk-β(A(2xk+1-xk)+Byk-b)))≥0,22222222222222222222222222222222222∀y′∈𝒴.
It can be further rewritten as
(20)g(y)-g(yk+1)+(y′-yk+1)T((s-βBTB)(yk+1-yk)-BT(λk-β(A(2xk+1-xk)+Byk+1-b)))≥0,22222222222222222222222222222222222222∀y′∈𝒴.
Substituting (***) into (19), we have
(21)g(y)-g(yk+1)+(y′-yk+1)T×(-BTλk+1+(s-βBTB)(yk+1-yk))≥0,2222222222222222222222222222222222∀y′∈𝒴.
In addition, it follows from (***) that
(22)Axk+1+Byk+1-b+A(xk+1-xk)+1β(λk+1-λk)=0.
Combining (18), (21), and (22) together, we get(23)h(u)-h(uk+1)+(x-xk+1y-yk+1λ-λk+1){(-ATλk+1-BTλk+1Axk+1+Byk+1-b)+(r(xk+1-xk)+AT(λk+1-λk)(s-βBTB)(yk+1-yk)A(xk+1-xk)+1β(λk+1-λk))}≥0,with w=(x,y,λ)T∈Ω. Using the notation of H, the assertion (16) is proved.Note that whenβ>0 and r>β∥ATA∥, s>β∥BTB∥, H is a positive definite matrix. Lemma 3 implies that the proposed algorithm can be viewed as a customized version of the PPA, where H is the metric proximal parameter. Henceforth, we denote the proposed algorithm as CPPA. Now we state and prove the contractive property of our algorithm.Lemma 4.
Suppose the condition(24)β>0,r>β∥ATA∥,s>β∥BTB∥
holds. Then, for any w*=(x*,y*,λ*)T∈Ω*, the sequence wk+1=(xk+1,yk+1,λk+1)T generated by ((*)–(***)) satisfies the following inequality:
(25)∥wk+1-w*∥H2≤∥wk-w*∥H2-∥wk+1-wk∥H2,
where the norm ∥·∥H2 is defined as ∥w∥H2=wTHw.Proof.
Note thatw*∈Ω; it follows from (16) that
(26)h(u*)-h(uk+1)+(w*-wk+1)T×(F(wk+1)+H(wk+1-wk))≥0,
which can be further written as
(27)(w*-wk+1)TH(wk+1-wk)≥h(uk+1)-h(u*)+(wk+1-w*)TF(wk+1).
On the other hand, using the fact that F(·) is a monotone operator, we have
(28)(wk+1-w*)T(F(wk+1)-F(w*))≥0.
Combining (28) and (27), we have
(29)(w*-wk+1)TH(wk+1-wk)≥h(uk+1)-h(u*)+(wk+1-w*)TF(w*).
Since w* is a solution of (9a) and wk+1∈Ω, we have
(30)h(uk+1)-h(u*)+(wk+1-w*)TF(w*)≥0.
Thus, we get
(31)(w*-wk+1)TH(wk+1-wk)≥0.
Since H is positive definite under the requirements of the parameters (24), (31) can be further written as
(32)(w*-wk)TH(wk+1-wk)≥∥wk+1-wk∥H2.
On the other hand,
(33)∥wk+1-w*∥H2=∥wk+1-wk+wk-w*∥H2=∥wk-w*∥H2+2(wk-w*)TH(wk+1-wk)+∥wk+1-wk∥H2.
Inserting (32) into (33), we obtain
(34)∥wk+1-w*∥H2≤∥wk-w*∥H2-∥wk+1-wk∥H2.
This completes the proof.Corollary 5.
Letw* be an arbitrary point in Ω*, and let the sequence {wk} be generated by ((*)–(***)). Then(1)
the sequence{wk} is bounded;(2)
the sequence{∥wk-w*∥H} is nonincreasing;(3)
lim⁡k→∞{∥wk+1-wk∥H}=0.We are now ready to prove the global convergence of the new method.Theorem 6.
Let the sequence{wk} be generated by the proposed algorithm and the condition (24) holds. Then the sequence {wk} converges to an optimal primal-dual solution for (1) from any starting point.Proof.
It follows from Corollary5 that
(35)lim⁡k→∞∥xk+1-xk∥=0,lim⁡k→∞∥yk+1-yk∥=0,lim⁡k→∞∥λk+1-λk∥=0.
Since {wk} is bounded, it has at least one cluster point. Suppose wk has a subsequence {wkj} that converges to w∞={x∞,y∞,λ∞}. It follows from (23) that
(36)lim⁡j→∞h(u)-h(ukj)+(x-xkjy-ykjλ-λkj)T{(-ATλkj-BTλkjAxkj+Bykj-b)}≥0.
Consequently, we have
(37)h(u)-h(u∞)+(x-x∞y-y∞λ-λ∞)T{(-ATλ∞-BTλ∞Ax∞+By∞-b)}≥0,
which implies that w∞∈Ω*.
Becauselim⁡k→∞{∥wk+1-wk∥H}=0, for any given ϵ>0, there exists l0>0 such that
(38)∥wkl+1-wkl∥H<ϵ2,∀l>l0.
Since wkj→w∞, for the ϵ>0 given above, there is an integer kl>l0, such that,
(39)∥wkl-w∞∥H<ϵ2.
Therefore, for any k>kl, it follows from (38) and (39) that
(40)∥wk-w∞∥H≤∥wkl-w∞∥H≤∥wkl+1-wkl∥H+∥wkl+1-w∞∥H<ϵ.
This implies that the sequence {wk} converges to w∞∈Ω*.


## 5. Application to the Stable Principal Component Pursuit Problem

In this section, we apply the proposed method to solve the stable principal component pursuit problem with nonnegative constraints (SPCP). The problem tested is from Example 2 of [19]. Our code was written in Matlab R2009b and all experiments were performed on a laptop with Intel Core 2 Duo @ 2.0 GHz CPU and 2 GB of memory.LetM=L+S+Z be a given observation matrix, where L is a low-rank and non-negative matrix, S is a sparse matrix, Z is a noise matrix. SPCP arising from image processing seeks to recover L and S by solving the following nonsmooth convex optimization problem:
(41)min⁡L,S,Z⁡∥L∥*+ρ∥S∥1+ℐ(∥Z∥F≤σ)+ℐ(L≥0)s.t.L+S+Z=M,
where ∥·∥* is the nuclear norm (defined as the sum of all singular values), ∥·∥1 and ∥·∥F denote the l1 norm and the Frobenius norm of a matrix, respectively, and ℐ(·) is the indicator function for the nonnegative orthant ℜm×n. In the above model, Z denotes the noise matrix and ρ and σ are some fixed parameters.Following the procedure described in [19], by introducing an auxiliary variable K and grouping L and S as one big block [L;S] and grouping Z and K as another big block [Z;K], model (41) can be easily reformulated as
(42)min⁡L,S,Z,K∥L∥*+ρ∥S∥1+ℐ(∥Z∥F≤σ)+ℐ(K≥0)s.t.(III0)(LS)+(I00-I)(ZK)=(M0).
which fits the setting of (1). The proposed algorithm ((*)–(***)) is therefore applicable for (42), and we obtain the following iterative scheme:
(43)Lk+1=argmin⁡⁡∥L∥*+r2∥L-Lk-1r(Λ1k+Λ2k)∥F2,(44)Sk+1=argmin⁡⁡ρ∥S∥1+r2∥S-Sk-1rΛ1k∥F2,(45)Zk+1=argmin⁡⁡ℐ(∥Z∥F≤σ)+s2∥Z-Zk-1s(Λ1k-β(2Lk+1-Lk+2Sk+1-Sk+Zk222222222222222222222222-Kk-M))1s∥F2,(46)Kk+1=argmin⁡⁡ℐ(K≥0)+s2∥K-Kk+1s(Λ2k-β(2Lk+1-Lk-Kk))∥F2,(47)Λ1k+1=Λ1k-β(2Lk+1-Lk+2Sk+1-Sk+Zk+1-M),(48)Λ2k+1=Λ2k-β(2Lk+1-Lk-Kk+1).
The main advantage of CPPA applied to SPCP is that the generated minimizations in ((43)–(46)) are all simple enough to have closed-form solutions. For completeness, we elaborate on the strategy of solving the resulted subproblems at each iteration.(i)
TheL-subproblem (43) amounts to evaluate the proximal operator of the nuclear norm function and is given by the matrix shrinkage operation:
(49)Lk+1:=MatShrink(Lk+1r(Λ1k+Λ2k),1r),
where the matrix shrinkage operator MatShrink(M,ξ)(ξ>0) is defined as
(50)MatShrink(M,ξ):=UDiag(max⁡{σ-ξ,0})VT,
and UDiag(σ)VT is the SVD of matrix M.(ii)
The closed-form solution ofS-subproblem (44) can be given by the l1 shrinkage operation:
(51)Sk+1:=Shrink(Sk+1rΛ1k,ρr),
where the l1 shrinkage operator Shrink(M,ξ) is defined as
(52)[Shrink(M,ξ)]ij:={Mij-ξ,ifMij>ξMij+ξ,ifMij<-ξ0,if|Mij|≤ξ.(iii)
TheZ-subproblem (45) amounts to projecting the matrix Wk:=M+(1/β)Λ1k-(2Lk+1-Lk+2Sk+1-Sk) onto the Euclidean ball ∥Z∥F≤σ, whose closed-form solution is given by
(53)Zk+1:=Wkmax⁡{1,∥Wk∥F/σ}.(iv)
TheK-subproblem (46) amounts to projecting the matrix 2Lk+1-Lk-(1/β)Λ2k onto the the nonnegative orthant, whose closed-form solution is given by:
(54)Kk+1:=max⁡{2Lk+1-Lk-1βΛ2k,0}.For detailed analytical methods of ((43)–(46)), the reader is referred to, for example, [19, 26]. Thus, the simplicity assumption (5) holds for the application (41).For the numerical experiments, we follow [19] to randomly generate the data of (41). We consider the scenario of m=n. For given n, r<n, we generate L*=R1R2T, where R1 and R2 are independent n×r full row rank matrices whose elements are independently and identically (i.i.d.) uniformly distributed in [0,1]. Note that in this experiment, L* is a component-wise nonnegative and low-rank matrix to be recovered. The support of the sparse matrix S* was chosen uniformly and randomly, and the nonzero entries in S* are generated i.i.d. uniformly in [-500,500]. The entries of matrix Z* for noise were generated as i.i.d. Gaussian with standard deviation 10-4; we set M:=L*+S*+Z*.In the following, we compare CPPA with APGM in [18, 19], since they are all designed based on the “simple” assumption and share the same conditions for convergence. For the penalty parameter β, we take β=0.01. For other individual parameters required by these methods, we choose r=2.618β, s=β, and ρ:=1/n. To implement all the compared methods, the initial iterate for both CPPA and APGM is chosen L0=K0=-M, S0=Z0=0, Λ10=Λ20=0. The stopping criterion is set as
(55)resid:=∥L+S+Z-M∥F∥M∥F<ϵr,
where ϵr is the tolerance set as ϵr=10-4. We denoted Rankr:=r/n so that the rank of L* is n*Rankr and Cardr:=cardinality(S*)/(n2) so that the cardinality of S* is n2*Cardr.Some preliminary numerical results are reported in Table1. Since they are synthetic examples with random date, for each scenario we test for 10 times, and the results were averaged over ten runs. Specifically, we reported the number of iteration (Iter.), relative error of the low-rank matrix L(relL), relative error of the sparse matrix S(relS), and CPU times in seconds (CPU(s)), where the relative errors are defined as
(56)relL:=∥L-L*∥F∥L*∥F,relS:=∥S-S*∥F∥S*∥F.Table 1
Numerical results for stable principal component pursuit problem.
n
APGM
CPPA
Iter.
rel
L
rel
S
CPU(s)
Iter.
rel
L
rel
S
CPU(s)
Rank
r
=
0.01
,
Card
r
=
0.01
50
84
9.36
e
-
003
2.67
e
-
005
0.2
88
9.15
e
-
003
2.29
e
-
005
0.2
100
44
8.41
e
-
003
2.47
e
-
005
0.6
52
8.38
e
-
003
1.87
e
-
005
0.6
150
58
5.00
e
-
003
2.16
e
-
005
2.1
57
4.89
e
-
003
1.87
e
-
005
2.0
200
52
4.81
e
-
003
2.01
e
-
005
4.1
45
4.42
e
-
003
2.36
e
-
005
3.4
250
52
3.40
e
-
003
1.97
e
-
005
7.4
46
3.14
e
-
003
2.42
e
-
005
6.4
300
47
3.30
e
-
003
1.90
e
-
005
11.4
42
2.53
e
-
003
3.58
e
-
005
9.8
400
43
2.57
e
-
003
1.83
e
-
005
39.1
42
1.66
e
-
003
4.18
e
-
005
37.7
500
38
2.08
e
-
003
2.77
e
-
005
82.1
42
1.13
e
-
003
4.68
e
-
005
89.2
Rank
r
=
0.02
,
Card
r
=
0.02
50
65
1.38
e
-
002
3.24
e
-
005
0.2
78
1.32
e
-
002
2.68
e
-
005
0.2
100
72
7.05
e
-
003
2.71
e
-
005
1.1
68
6.77
e
-
003
2.27
e
-
005
1.0
150
63
4.88
e
-
003
2.69
e
-
005
2.8
48
4.61
e
-
003
2.42
e
-
005
2.1
200
53
3.66
e
-
003
2.52
e
-
005
7.4
43
3.24
e
-
003
4.50
e
-
005
5.4
250
47
2.98
e
-
003
2.62
e
-
005
7.3
43
2.10
e
-
003
4.51
e
-
005
6.7
300
41
2.55
e
-
003
3.02
e
-
005
10.4
43
1.58
e
-
003
5.38
e
-
005
10.8
400
34
1.70
e
-
003
4.95
e
-
005
32.5
44
1.14
e
-
003
5.74
e
-
005
41.9
500
34
7.67
e
-
004
4.88
e
-
005
75.3
45
9.18
e
-
004
5.73
e
-
005
101.0
Rank
r
=
0.03
,
Card
r
=
0.03
50
98
9.03
e
-
003
3.91
e
-
005
0.4
88
8.86
e
-
003
3.09
e
-
005
0.3
100
78
6.41
e
-
003
3.18
e
-
005
1.3
64
6.13
e
-
003
2.62
e
-
005
1.0
150
56
3.67
e
-
003
3.42
e
-
005
2.5
51
3.41
e
-
003
2.97
e
-
005
2.2
200
45
3.09
e
-
003
3.32
e
-
005
4.2
46
2.54
e
-
003
3.94
e
-
005
4.0
250
38
2.16
e
-
003
4.35
e
-
005
6.1
46
1.67
e
-
003
4.94
e
-
005
7.1
300
35
1.71
e
-
003
4.86
e
-
005
8.7
46
1.32
e
-
003
5.39
e
-
005
11.5
400
34
1.22
e
-
003
8.13
e
-
005
31.9
47
9.56
e
-
004
5.85
e
-
005
45.4
500
34
9.29
e
-
004
8.18
e
-
005
74.5
47
8.43
e
-
004
6.44
e
-
005
106.4From Table1, we observe that both of the methods are efficient. Our method tends to be competitive with APGM and performs reasonably well for moderate dimensions (say, n∈[100–300]). To see the comparison clearly, we focus on the particular case where n=150 and Rankr=0.02, Cardr=0.02; we visualize the iterative processes of different method in Figure 1. More specifically, we plot the evolutions of the relative error relL and relS, with respect to the iterations. According to the curves in Figure 1, the performance of CPPA is slightly worse than that of APGM, when the iterations are small. However, with the iterations going large, CPPA shows a better performance than APGM.(a) Relative errors of low-rank matrixL. (b) Relative errors of sparse matrix S.
(a)
(b)


## 6. Concluding Remarks

We have proposed a new algorithm for solving (1) which admits easy subproblems assuming the proximal mappings of f and g are easy to compute. Our algorithm can be viewed as a customized proximal point method with a special proximal regularization parameter. We established its global convergence under the analytic framework of contraction type methods. The computational results on solving the stable principal component pursuit problem show that our algorithm works reasonably well on large-scale instances.



---

*Source: 396753-2014-03-30.xml*