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

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.

1. Introduction

In this paper, we consider the following separable optimization problem: min_{𝑥∈𝒳,𝑦∈𝒴}𝑓(𝑥)+𝑔𝑦s.t.𝐴𝑥+𝐵𝑦=𝑏, (1) where 𝐴∈ℜ^𝑚×𝑛 and 𝐵∈ℜ^𝑚×𝑝 are given matrices, 𝑏∈ℜ^𝑚 is a given vector, and 𝒳⊂ℜ^𝑛 and 𝒴⊂ℜ^𝑝 are nonempty closed convex sets. 𝑓∶ℜ^𝑛→ℜ and 𝑔∶ℜ^𝑝→ℜ 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 (𝑥⁰,𝑦⁰,𝜆⁰), ALM iterates via the following procedure:𝑥^𝑘+1,𝑦^𝑘+1=argmin_𝑥,𝑦ℒ_𝒜𝑥,𝑦;𝜆^𝑘,𝜆^𝑘+1=𝜆^𝑘−𝛽𝐴𝑥^𝑘+1+𝐵𝑦^𝑘+1−𝑏, (2a)(2b)where the augmented Lagrangian function ℒ_𝒜(𝑥,𝑦;𝜆) associated with the problem (1) is given by ℒ_𝒜𝑥,𝑦,𝜆=𝑓(𝑥)+𝑔𝑦−𝜆^𝑇+𝛽𝐴𝑥+𝐵𝑦−𝑏2𝐴𝑥+𝐵𝑦−𝑏², (3) and 𝜆∈ℜ^𝑚 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 𝑓 and 𝑔 in (2a) are strongly coupled because of the quadratic term (𝛽/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 𝑥 and 𝑦 in Gauss-Seidel manner. More specifically, at each iteration, ADM takes the following form:𝑥^𝑘+1=argmin_𝑥∈𝒳ℒ_𝒜𝑥,𝑦^𝑘;𝜆^𝑘,𝑦^𝑘+1=argmin_𝑦∈𝒴ℒ_𝒜𝑥^𝑘+1,𝑦;𝜆^𝑘,𝜆^𝑘+1=𝜆^𝑘−𝛽𝐴𝑥^𝑘+1+𝐵𝑦^𝑘+1−𝑏. (4a)(4b)(4c)

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 𝑥 and 𝑦 alternatingly at each iteration, rather than with respect to both 𝑥 and 𝑦 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.

When 𝐴 and 𝐵 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 𝑓 and 𝑔, respectively. Here, the proximal operator of the function 𝜃∶ℜ^𝑛→ℜ is defined by 𝑢∶=𝐩𝐫𝐨𝐱_𝜉𝜃(𝑎)𝜃(𝑢)+1=argmin‖‖2𝜉𝑢−𝑎²∣𝑢∈ℜ^𝑛, (5) where 𝑎∈ℜ^𝑛 and 𝜉>0. In most popular applications of sparse optimization, the proximal operator can be computed exactly and efficiently (e.g., 𝜃(𝑥)=‖𝑥‖₁∶=∑^𝑛_𝑖=1|𝑥_𝑖|or𝜃(𝑥)=‖𝑥‖₂∶=∑^𝑛_𝑖=1|𝑥_𝑖|²). While when 𝐴 (resp., 𝐵) is not an identity matrix, the resulting ADM subproblems may not be easily solvable, since they involve inverting of 𝐴 (resp., 𝐵), 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 𝐴 and 𝐵 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 (𝑥^𝑘+1,𝑦^𝑘+1,𝜆^𝑘+1) via the following procedure:𝑥^𝑘+1=argmin_𝑥∈𝒳𝑓(𝑥)+𝑟2𝑥−𝑥^𝑘−1𝑟𝐴^𝑇𝜆^𝑘−𝛽𝐴𝑥^𝑘+𝐵𝑦^𝑘−𝑏²,𝑦^𝑘+1=argmin_𝑦∈𝒴𝑔𝑦+𝑠2𝑦−𝑦^𝑘−1𝑠𝐵^𝑇𝜆^𝑘−𝛽𝐴𝑥^𝑘+1+𝐵𝑦^𝑘−𝑏²,𝜆^𝑘+1=𝜆^𝑘−𝛽𝐴𝑥^𝑘+1+𝐵𝑦^𝑘+1−𝑏, (6a)(6b)(6c)where 𝛽>0, 𝑟>𝛽‖𝐴^𝑇𝐴‖, and 𝑠>𝛽‖𝐵^𝑇𝐵‖.

In this paper, we also focus on the specific scenario of (1), where both 𝐴 and 𝐵 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 Section 2, 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) where 𝜆∈ℜ^𝑚 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 (𝑥^∗,𝑦^∗,𝜆^∗) be a saddle point of ℒ. We have ℒ_𝜆∈ℜ^𝑚𝑥^∗,𝑦^∗,𝜆≤ℒ𝑥^∗,𝑦^∗,𝜆^∗≤ℒ_{𝑥∈𝒳,𝑦∈𝒴}𝑥,𝑦,𝜆^∗. (8) Then we can directly read off the optimality conditions with variational characterization. More specifically, 𝑤^∗=(𝑥^∗,𝑦^∗,𝜆^∗)∈Ω is a saddle point of ℒ if and only if it satisfies the following mixed variational inequality (VI):ℎ(𝑢)−ℎ𝑢^∗+𝑤^′−𝑤^∗𝑇𝐹𝑤^∗≥0,∀𝑤^′∈Ω, (9a) with 𝑥𝑦𝑥𝑦𝜆(𝑢)(𝑥)𝑦,𝐹(𝑤)𝑢∶=,𝑤∶=,ℎ∶=𝑓+𝑔∶=−𝐴^𝑇𝜆−𝐵^𝑇𝜆,𝐴𝑥+𝐵𝑦−𝑏Ω=𝒳×𝒴×ℜ^𝑚. (9b)(9c) Hence, a solution of ((9a)–(9c)) yields a solution of (1).

Note that the mapping 𝐹(⋅) is said to be monotone with respect to Ω if 𝐹(𝑤)−𝐹𝑤^′𝑇𝑤−𝑤^′≥0,∀𝑤,𝑤^′∈Ω. (10) Consequently, it can be easily verified that 𝐹(𝑤) 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 𝑤^𝑘∈Ω, PPA generates the new iterative 𝑤^𝑘+1 via the following procedure: ℎ𝑢^′−ℎ𝑢^𝑘+1+𝑤^′−𝑤^𝑘+1𝑇×𝐹𝑤^𝑘+1+𝐺𝑤^𝑘+1−𝑤^𝑘≥0,∀𝑤^′∈Ω, (11) where 𝐺 is a positive definite matrix, playing the role of proximal regularization parameter. A simple choice of 𝐺 is that 𝐺=𝛽⋅𝐼 where 𝛽>0 and 𝐼 is the identity matrix, regularizing the proximal terms 𝑤^𝑘+1−𝑤^𝑘 in the uniform way. We refer the reader to example [21–24] for some special choices of 𝐺 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 give 𝑎∈ℜ^𝑛+𝑝 and 𝜉>0, the proximal operator of ℎ(𝑢) (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 of ℎ is “easy” to evaluate. Note that ℎ is separable across two variables, that is, ℎ=𝑓+𝑔; according to the definition (5), we have 𝐩𝐫𝐨𝐱_𝜉ℎ(𝑏,𝑐)=𝐩𝐫𝐨𝐱_𝜉𝑓(𝑏),𝐩𝐫𝐨𝐱_𝜉𝑔(𝑐), (12) where 𝑏∈ℜ^𝑛 and 𝑐∈ℜ^𝑝. Hence, under the assumption, the proximal operators of the component objective functions 𝑓 and 𝑔 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 of 𝑥, we first add a proximal regularization parameter 𝐺=^{𝑟𝐼0𝐴𝑇000}_{𝐴0(1/𝛽)𝐼} to the variational characterization of the problem ((9a)–(9c)). Then, like the algorithm in [21], we get the following proximal point algorithm:𝑥^𝑘+1=argmin_𝑥∈𝒳𝑓(𝑥)+𝑟2𝑥−𝑥^𝑘−1𝑟𝐴^𝑇𝜆^𝑘2,𝑦^𝑘+1=argmin_𝑦∈𝒴𝑔𝑦+𝛽2𝐴2𝑥^𝑘+1−𝑥^𝑘𝜆+𝐵𝑦−𝑏−^𝑘𝛽²,𝜆^𝑘+1=𝜆^𝑘−𝛽𝐴𝑥^𝑘+1+𝐵𝑦^𝑘+1−𝑏. (13a)(13b)(13c) However, (13b) is still not implementable. Inspired in [25], in order to alleviate the computation required by 𝑦-subproblems, we try to linearize the quadratic term in (13b) by 𝛽2𝐴2𝑥^𝑘+1−𝑥^𝑘𝜆+𝐵𝑦−𝑏−^𝑘𝛽²≈𝛽2𝐴2𝑥^𝑘+1−𝑥^𝑘+𝐵𝑦^𝑘𝜆−𝑏−^𝑘𝛽²+𝛽𝑦−𝑦^𝑘𝑇𝑔^𝑘+𝑠2𝑦−𝑦^𝑘, (14) where 𝑠>0 is a proximal parameter, and 𝑔^𝑘=𝐵^𝑇(𝐴(2𝑥^𝑘+1−𝑥^𝑘)+𝐵𝑦^𝑘−𝑏−𝜆^𝑘/𝛽) is the gradient of (1/2)‖𝐴(2𝑥^𝑘+1−𝑥^𝑘)+𝐵𝑦−𝑏−𝜆^𝑘/𝛽‖² at 𝑦^𝑘. With a simple manipulation, we get the following approximation to (13b): 𝑦^𝑘+1=argmin_𝑦∈𝒴𝑔𝑦+𝑠2𝑦−𝑦^𝑘−1𝑠𝐵^𝑇𝜆^𝑘×𝐴−𝛽2𝑥^𝑘+1−𝑥^𝑘+𝐵𝑦^𝑘−𝑏². (15) 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 Algorithm 1.

The proposed algorithm for (1) Step  0. Input 𝑥0,𝑦0,𝜆0∈𝒳×𝒴. Step  1. Set      𝑥𝑘+1=argmin𝑥∈𝒳𝑓(𝑥)+𝑟2𝑥−𝑥𝑘−1𝑟𝐴𝑇𝜆𝑘2,                 (∗)      𝑦𝑘+1=argmin𝑦∈𝒴𝑔𝑦+𝑠2𝑦−𝑦𝑘−1𝑠𝐵𝑇𝜆𝑘−𝛽𝐴2𝑥𝑘+1−𝑥𝑘+𝐵𝑦𝑘−𝑏2,    (∗∗)      𝜆𝑘+1=𝜆𝑘−𝛽𝐴2𝑥𝑘+1−𝑥𝑘+𝐵𝑦𝑘+1−𝑏,                  (∗∗∗) where 𝛽>0, and 𝑟>𝛽𝐴𝑇𝐴,𝑠>𝛽𝐵𝑇𝐵. Step  2. If a termination criterion is not met, go to Step  1.

Algorithm 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 𝑓 and 𝑔. 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. Let 𝑤^𝑘+1=(𝑥^𝑘+1,𝑦^𝑘+1,𝜆^𝑘+1)^𝑇 be generated by ((∗)–(∗∗∗)) from a given 𝑤^𝑘=(𝑥^𝑘,𝑦^𝑘,𝜆^𝑘)^𝑇 and 𝐻=^{𝑟𝐼0𝐴𝑇0𝑠𝐼−𝛽𝐵𝑇𝐵0}_{𝐴0(1/𝛽)𝐼}. Then ℎ𝑢^′−ℎ𝑢^𝑘+1+𝑤^′−𝑤^𝑘+1𝑇×𝐹𝑤^𝑘+1+𝐻𝑤^𝑘+1−𝑤^𝑘≥0,∀𝑤^′∈Ω. (16)

Proof. First, by deriving the optimality condition for the subproblem (∗), we have 𝑓(𝑥)−𝑓𝑥^𝑘+1+𝑥^′−𝑥^𝑘+1𝑇𝑟𝑥^𝑘+1−𝑥^𝑘−𝐴^𝑇𝜆^𝑘≥0,∀𝑥^′∈𝒳. (17) It can be further rewritten as 𝑓(𝑥)−𝑓𝑥^𝑘+1+𝑥^′−𝑥^𝑘+1𝑇×−𝐴^𝑇𝜆^𝑘+1+𝑟𝑥^𝑘+1−𝑥^𝑘+𝐴^𝑇𝜆^𝑘+1−𝜆^𝑘≥0,∀𝑥^′∈𝒳. (18) Similarly, the optimality condition for the subproblem (∗∗) is 𝑔𝑦−𝑔𝑦^𝑘+1+𝑦^′−𝑦^𝑘+1𝑇𝑠𝑦^𝑘+1−𝑦^𝑘−𝐵^𝑇𝜆^𝑘−𝛽𝐴2𝑥^𝑘+1−𝑥^𝑘+𝐵𝑦^𝑘−𝑏≥0,∀𝑦^′∈𝒴. (19) It can be further rewritten as 𝑔𝑦−𝑔𝑦^𝑘+1+𝑦^′−𝑦^𝑘+1𝑇𝑠−𝛽𝐵^𝑇𝐵𝑦^𝑘+1−𝑦^𝑘−𝐵^𝑇𝜆^𝑘−𝛽𝐴2𝑥^𝑘+1−𝑥^𝑘+𝐵𝑦^𝑘+1−𝑏≥0,∀𝑦^′∈𝒴. (20) Substituting (∗∗∗) into (19), we have 𝑔𝑦−𝑔𝑦^𝑘+1+𝑦^′−𝑦^𝑘+1𝑇×−𝐵^𝑇𝜆^𝑘+1+𝑠−𝛽𝐵^𝑇𝐵𝑦^𝑘+1−𝑦^𝑘≥0,∀𝑦^′∈𝒴. (21) In addition, it follows from (∗∗∗) that 𝐴𝑥^𝑘+1+𝐵𝑦^𝑘+1−𝑏+𝐴𝑥^𝑘+1−𝑥^𝑘+1𝛽𝜆^𝑘+1−𝜆^𝑘=0. (22) Combining (18), (21), and (22) together, we getℎ(𝑢)−ℎ𝑢^𝑘+1+𝑥−𝑥^𝑘+1𝑦−𝑦^𝑘+1𝜆−𝜆^𝑘+1−𝐴^𝑇𝜆^𝑘+1−𝐵^𝑇𝜆^𝑘+1𝐴𝑥^𝑘+1+𝐵𝑦^𝑘+1+𝑟𝑥−𝑏^𝑘+1−𝑥^𝑘+𝐴^𝑇𝜆^𝑘+1−𝜆^𝑘𝑠−𝛽𝐵^𝑇𝐵𝑦^𝑘+1−𝑦^𝑘𝐴𝑥^𝑘+1−𝑥^𝑘+1𝛽𝜆^𝑘+1−𝜆^𝑘≥0, (23)with 𝑤=(𝑥,𝑦,𝜆)^𝑇∈Ω. Using the notation of 𝐻, the assertion (16) is proved.

Note that when 𝛽>0 and 𝑟>𝛽‖𝐴^𝑇𝐴‖, 𝑠>𝛽‖𝐵^𝑇𝐵‖, 𝐻 is a positive definite matrix. Lemma 3 implies that the proposed algorithm can be viewed as a customized version of the PPA, where 𝐻 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 𝛽>0,𝑟>𝛽𝐴^𝑇𝐴,𝑠>𝛽𝐵^𝑇𝐵 (24) holds. Then, for any 𝑤^∗=(𝑥^∗,𝑦^∗,𝜆^∗)^𝑇∈Ω^∗, the sequence 𝑤^𝑘+1=(𝑥^𝑘+1,𝑦^𝑘+1,𝜆^𝑘+1)^𝑇 generated by ((∗)–(∗∗∗)) satisfies the following inequality: 𝑤^𝑘+1−𝑤^∗2_𝐻≤𝑤^𝑘−𝑤^∗2_𝐻−𝑤^𝑘+1−𝑤^𝑘2_𝐻, (25) where the norm ‖⋅‖²_𝐻 is defined as ‖𝑤‖²_𝐻=𝑤^𝑇𝐻𝑤.

Proof. Note that 𝑤^∗∈Ω; it follows from (16) that ℎ𝑢^∗−ℎ𝑢^𝑘+1+𝑤^∗−𝑤^𝑘+1𝑇×𝐹𝑤^𝑘+1+𝐻𝑤^𝑘+1−𝑤^𝑘≥0, (26) which can be further written as 𝑤^∗−𝑤^𝑘+1𝑇𝐻𝑤^𝑘+1−𝑤^𝑘≥ℎ𝑢^𝑘+1−ℎ𝑢^∗+𝑤^𝑘+1−𝑤^∗𝑇𝐹𝑤^𝑘+1. (27) On the other hand, using the fact that 𝐹(⋅) is a monotone operator, we have 𝑤^𝑘+1−𝑤^∗𝑇𝐹𝑤^𝑘+1−𝐹𝑤^∗≥0. (28) Combining (28) and (27), we have 𝑤^∗−𝑤^𝑘+1𝑇𝐻𝑤^𝑘+1−𝑤^𝑘≥ℎ𝑢^𝑘+1−ℎ𝑢^∗+𝑤^𝑘+1−𝑤^∗𝑇𝐹𝑤^∗. (29) Since 𝑤^∗ is a solution of (9a) and 𝑤^𝑘+1∈Ω, we have ℎ𝑢^𝑘+1−ℎ𝑢^∗+𝑤^𝑘+1−𝑤^∗𝑇𝐹𝑤^∗≥0. (30) Thus, we get 𝑤^∗−𝑤^𝑘+1𝑇𝐻𝑤^𝑘+1−𝑤^𝑘≥0. (31) Since 𝐻 is positive definite under the requirements of the parameters (24), (31) can be further written as 𝑤^∗−𝑤^𝑘𝑇𝐻𝑤^𝑘+1−𝑤^𝑘≥𝑤^𝑘+1−𝑤^𝑘2_𝐻. (32) On the other hand, 𝑤^𝑘+1−𝑤^∗2_𝐻=𝑤^𝑘+1−𝑤^𝑘+𝑤^𝑘−𝑤^∗2_𝐻=𝑤^𝑘−𝑤^∗2_𝐻+2𝑤^𝑘−𝑤^∗𝑇𝐻𝑤^𝑘+1−𝑤^𝑘+𝑤^𝑘+1−𝑤^𝑘2_𝐻. (33) Inserting (32) into (33), we obtain 𝑤^𝑘+1−𝑤^∗2_𝐻≤𝑤^𝑘−𝑤^∗2_𝐻−𝑤^𝑘+1−𝑤^𝑘2_𝐻. (34) This completes the proof.

Corollary 5. Let 𝑤^∗ be an arbitrary point in Ω^∗, and let the sequence {𝑤^𝑘} be generated by ((∗)–(∗∗∗)). Then(1)the sequence {𝑤^𝑘} is bounded;(2)the sequence {‖𝑤^𝑘−𝑤^∗‖_𝐻} is nonincreasing;(3)lim_𝑘→∞{‖𝑤^𝑘+1−𝑤^𝑘‖_𝐻}=0.

We are now ready to prove the global convergence of the new method.

Theorem 6. Let the sequence {𝑤^𝑘} be generated by the proposed algorithm and the condition (24) holds. Then the sequence {𝑤^𝑘} converges to an optimal primal-dual solution for (1) from any starting point.

Proof. It follows from Corollary 5 that lim_𝑘→∞𝑥^𝑘+1−𝑥^𝑘=0,lim_𝑘→∞𝑦^𝑘+1−𝑦^𝑘=0,lim_𝑘→∞𝜆^𝑘+1−𝜆^𝑘=0. (35) Since {𝑤^𝑘} is bounded, it has at least one cluster point. Suppose 𝑤^𝑘 has a subsequence {𝑤^𝑘𝑗} that converges to 𝑤^∞={𝑥^∞,𝑦^∞,𝜆^∞}. It follows from (23) that lim_𝑗→∞ℎ(𝑢)−ℎ𝑢^𝑘𝑗+𝑥−𝑥^𝑘𝑗𝑦−𝑦^𝑘𝑗𝜆−𝜆^𝑘𝑗𝑇−𝐴^𝑇𝜆^𝑘𝑗−𝐵^𝑇𝜆^𝑘𝑗𝐴𝑥^𝑘𝑗+𝐵𝑦^𝑘𝑗−𝑏≥0. (36) Consequently, we have ℎ(𝑢)−ℎ𝑢^∞+𝑥−𝑥^∞𝑦−𝑦^∞𝜆−𝜆^∞𝑇−𝐴^𝑇𝜆^∞−𝐵^𝑇𝜆^∞𝐴𝑥^∞+𝐵𝑦^∞−𝑏≥0, (37) which implies that 𝑤^∞∈Ω^∗.
Because lim_𝑘→∞{‖𝑤^𝑘+1−𝑤^𝑘‖_𝐻}=0, for any given 𝜖>0, there exists 𝑙₀>0 such that 𝑤^𝑘𝑙+1−𝑤^𝑘𝑙_𝐻<𝜖2,∀𝑙>𝑙₀. (38) Since 𝑤^𝑘𝑗→𝑤^∞, for the 𝜖>0 given above, there is an integer 𝑘_𝑙>𝑙₀, such that, 𝑤^𝑘𝑙−𝑤^∞_𝐻<𝜖2. (39) Therefore, for any 𝑘>𝑘_𝑙, it follows from (38) and (39) that 𝑤^𝑘−𝑤^∞_𝐻≤𝑤^𝑘𝑙−𝑤^∞_𝐻≤𝑤^𝑘𝑙+1−𝑤^𝑘𝑙_𝐻+𝑤^𝑘𝑙+1−𝑤^∞_𝐻<𝜖. (40) This implies that the sequence {𝑤^𝑘} converges to 𝑤^∞∈Ω^∗.

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.

Let 𝑀=𝐿+𝑆+𝑍 be a given observation matrix, where 𝐿 is a low-rank and non-negative matrix, 𝑆 is a sparse matrix, 𝑍 is a noise matrix. SPCP arising from image processing seeks to recover 𝐿 and 𝑆 by solving the following nonsmooth convex optimization problem: min_{𝐿,𝑆,𝑍}‖𝐿‖_∗+𝜌‖𝑆‖₁+ℐ‖𝑍‖_𝐹≤𝜎+ℐ(𝐿≥0)s.t.𝐿+𝑆+𝑍=𝑀, (41) where ‖⋅‖_∗ is the nuclear norm (defined as the sum of all singular values), ‖⋅‖₁ and ‖⋅‖_𝐹 denote the 𝑙₁ norm and the Frobenius norm of a matrix, respectively, and ℐ(⋅) is the indicator function for the nonnegative orthant ℜ^𝑚×𝑛. In the above model, 𝑍 denotes the noise matrix and 𝜌 and 𝜎 are some fixed parameters.

Following the procedure described in [19], by introducing an auxiliary variable 𝐾 and grouping 𝐿 and 𝑆 as one big block [𝐿;𝑆] and grouping 𝑍 and 𝐾 as another big block [𝑍;𝐾], model (41) can be easily reformulated as min_{𝐿,𝑆,𝑍,𝐾}‖𝐿‖_∗+𝜌‖𝑆‖₁+ℐ‖𝑍‖_𝐹()𝐿𝑆+𝑍𝐾=𝑀0.≤𝜎+ℐ𝐾≥0s.t.𝐼𝐼𝐼0𝐼00−𝐼 (42) which fits the setting of (1). The proposed algorithm ((∗)–(∗∗∗)) is therefore applicable for (42), and we obtain the following iterative scheme: 𝐿^𝑘+1=argmin‖𝐿‖_∗+𝑟2𝐿−𝐿^𝑘−1𝑟Λ^𝑘₁+Λ^𝑘₂²_𝐹,𝑆^𝑘+1=argmin𝜌‖𝑆‖₁+𝑟2𝑆−𝑆^𝑘−1𝑟Λ^𝑘₁²_𝐹,𝑍^𝑘+1=argminℐ‖𝑍‖_𝐹+𝑠≤𝜎2𝑍−𝑍^𝑘−1𝑠Λ^𝑘₁−𝛽2𝐿^𝑘+1−𝐿^𝑘+2𝑆^𝑘+1−𝑆^𝑘+𝑍^𝑘−𝐾^𝑘−𝑀²_𝐹,𝐾^𝑘+1()+𝑠=argminℐ𝐾≥02𝐾−𝐾^𝑘+1𝑠Λ^𝑘₂−𝛽2𝐿^𝑘+1−𝐿^𝑘−𝐾^𝑘2_𝐹,Λ₁^𝑘+1=Λ^𝑘₁−𝛽2𝐿^𝑘+1−𝐿^𝑘+2𝑆^𝑘+1−𝑆^𝑘+𝑍^𝑘+1,Λ−𝑀_2𝑘+1=Λ^𝑘₂−𝛽2𝐿^𝑘+1−𝐿^𝑘−𝐾^𝑘+1. (43)(44)(45)(46)(47)(48) 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)The 𝐿-subproblem (43) amounts to evaluate the proximal operator of the nuclear norm function and is given by the matrix shrinkage operation: 𝐿^𝑘+1∶=MatShrink𝐿^𝑘+1𝑟Λ^𝑘₁+Λ^𝑘₂,1𝑟, (49) where the matrix shrinkage operator MatShrink(𝑀,𝜉)(𝜉>0) is defined as MatShrink(𝑀,𝜉)∶=𝑈Diag(max{𝜎−𝜉,0})𝑉^𝑇, (50) and 𝑈Diag(𝜎)𝑉^𝑇 is the SVD of matrix 𝑀.(ii)The closed-form solution of 𝑆-subproblem (44) can be given by the 𝑙₁ shrinkage operation: 𝑆^𝑘+1∶=Shrink𝑆^𝑘+1𝑟Λ^𝑘₁,𝜌𝑟, (51) where the 𝑙₁ shrinkage operator Shrink(𝑀,𝜉) is defined as Shrink(𝑀,𝜉)_𝑖𝑗𝑀∶=^𝑖𝑗−𝜉,if𝑀_𝑖𝑗𝑀>𝜉_𝑖𝑗+𝜉,if𝑀_𝑖𝑗𝑀<−𝜉0,if_𝑖𝑗≤𝜉. (52)(iii)The 𝑍-subproblem (45) amounts to projecting the matrix 𝑊^𝑘∶=𝑀+(1/𝛽)Λ^𝑘₁−(2𝐿^𝑘+1−𝐿^𝑘+2𝑆^𝑘+1−𝑆^𝑘) onto the Euclidean ball ‖𝑍‖_𝐹≤𝜎, whose closed-form solution is given by 𝑍^𝑘+1𝑊∶=^𝑘max1,𝑊^𝑘_𝐹./𝜎 (53)(iv)The 𝐾-subproblem (46) amounts to projecting the matrix 2𝐿^𝑘+1−𝐿^𝑘−(1/𝛽)Λ^𝑘₂ onto the the nonnegative orthant, whose closed-form solution is given by: 𝐾^𝑘+1∶=max2𝐿^𝑘+1−𝐿^𝑘−1𝛽Λ^𝑘₂,0. (54)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 𝑚=𝑛. For given 𝑛, 𝑟<𝑛, we generate 𝐿^∗=𝑅₁𝑅^𝑇₂, where 𝑅₁ and 𝑅₂ are independent 𝑛×𝑟 full row rank matrices whose elements are independently and identically (i.i.d.) uniformly distributed in [0,1]. Note that in this experiment, 𝐿^∗ is a component-wise nonnegative and low-rank matrix to be recovered. The support of the sparse matrix 𝑆^∗ was chosen uniformly and randomly, and the nonzero entries in 𝑆^∗ are generated i.i.d. uniformly in [−500,500]. The entries of matrix 𝑍^∗ for noise were generated as i.i.d. Gaussian with standard deviation 10⁻⁴; we set 𝑀∶=𝐿^∗+𝑆^∗+𝑍^∗.

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 𝑟=2.618𝛽, 𝑠=𝛽, and 𝜌∶=1/√𝑛. To implement all the compared methods, the initial iterate for both CPPA and APGM is chosen 𝐿⁰=𝐾⁰=−𝑀, 𝑆⁰=𝑍⁰=0, Λ⁰₁=Λ⁰₂=0. The stopping criterion is set as ‖‖resid∶=𝐿+𝑆+𝑍−𝑀_𝐹‖𝑀‖_𝐹<𝜖_𝑟, (55) where 𝜖_𝑟 is the tolerance set as 𝜖_𝑟=10⁻⁴. We denoted Rank_𝑟∶=𝑟/𝑛 so that the rank of 𝐿^∗ is 𝑛∗Rank_𝑟 and Card_𝑟∶=cardinality(𝑆^∗)/(𝑛²) so that the cardinality of 𝑆^∗ is 𝑛²∗Card_𝑟.

Some preliminary numerical results are reported in Table 1. 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 𝐿(rel_𝐿), relative error of the sparse matrix 𝑆(rel_𝑆), and CPU times in seconds (CPU(s)), where the relative errors are defined as rel_𝐿∶=𝐿−𝐿^∗_𝐹‖𝐿^∗‖_𝐹,rel_𝑆∶=𝑆−𝑆^∗_𝐹𝑆^∗_𝐹. (56)

Table 1 

Numerical results for stable principal component pursuit problem.

From Table 1, 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, 𝑛∈[100–300]). To see the comparison clearly, we focus on the particular case where 𝑛=150 and Rank_𝑟=0.02, Card_𝑟=0.02; we visualize the iterative processes of different method in Figure 1. More specifically, we plot the evolutions of the relative error rel_𝐿 and rel_𝑆, 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)

(b)

Figure 1 

(a) Relative errors of low-rank matrix 𝐿. (b) Relative errors of sparse matrix 𝑆.

6. Concluding Remarks

We have proposed a new algorithm for solving (1) which admits easy subproblems assuming the proximal mappings of 𝑓 and 𝑔 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The work is supported in part by the Natural Science Foundation of China, Grant NSFC-70971136.