# A General Self-Adaptive Relaxed-PPA Method for Convex Programming with Linear Constraints

Authors: Xiaoling Fu
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We present an efficient
method for solving linearly constrained convex programming. Our algorithmic framework employs an implementable proximal step by a slight relaxation to the subproblem of proximal point algorithm (PPA). In particular, the stepsize choice condition of our algorithm is weaker than some elegant PPA-type methods. This condition is flexible  and effective. Self-adaptive strategies are proposed to  improve the convergence in practice. We theoretically show under mild conditions that our method converges  in a  global sense.  Finally, we discuss applications and perform numerical experiments which confirm the efficiency of the proposed method. Comparisons of our method with some state-of-the-art algorithms are also provided.


---

## Body

## 1.  Introduction

In this paper, we consider the following generic convex programming:(1)min{θ(x)∣Ax=b(orAx≥b),x∈𝒳},
where θ(x):ℜn→ℜ is a convex (not necessary smooth) function, A∈ℜm×n, b∈ℜm, and 𝒳⊂ℜn is a closed convex set. Problem (1) generalizes a wide range of problems that frequently arise in signal and image processing and reconstruction, mechanics, statistics, operations research, and other fields, for example, basis pursuit [1–4], nearest correlation matrix [5–7], matrix completion problem [2, 3, 8–10], and so forth. Before we begin, some assumptions should be presented for problem (1).Assumption 1.
The solution set of (1) is denoted by 𝒳*, and it is assumed to be nonempty.Assumption 2.
The objective function is simple. This means that, for a given constantϱ, the following proximal problem admits a closed-form solution or can be solved efficiently with high precision:
(2)minx∈𝒳θ(x)+ϱ2∥x-a∥2,
where a is any given vector. At first sight, this assumption seems to be quite restrictive, but this is indeed for many problems in practice. For example, nuclear norm function in matrix completion problem, l1-norm function in basis pursuit problem, and so forth.Many fundamental methods have been developed over the past decades to solve problem (1). Proximal point algorithm (PPA) is one of the leading approaches for solving convex optimization problems. It is earlier used for regularized linear equations and has been applied to convex optimization by Martinet [11]. There are some significant theoretical achievements [12–19] in the field of PPA for convex optimization and monotone variational inequalities (VIs). Nowadays, it is still the object of intensive investigation [20] and leads to a variety of primal and dual methods. Common to PPA and its variants is the difficulty of their subproblems; this restricts the practical interest. Augmented Lagrangian method (ALM) [21] is a powerful method for linearly constrained problems. It can be regarded as a variant of PPA applied to the dual problem of (1). However, with the additional regularized term ∥Ax-b∥2, its subproblems require an inverse operator of the form (I+(1/s)ATA)-1 which is hard to implement in some cases. Particularly, ATA is general or large scale, so the computation of inverse operator may fail. Hence, ALM is not sufficiently competitive when the objective function θ(x) is “simple.” Extragradient method (EGM) [22] is a practical method for (1) which employs the information of current iteration. In fact, EGM is an explicit type method and requires two calls to the gradient per iteration; therefore, it might suffer slow convergence. Recently, He and Yuan [23] proposed a contraction method based on PPA (PPA-CM) to solve (1), which is elegant and simple. Inspired by PPA-CM, a Lagrangian PPA-based (LPPA) contraction method is presented in [24] which employs an asymmetrical proximal term [25]. These two PPA-based methods have nice convergence properties that are similar in many ways to PPA, but they both require a quite restrictive condition for convergence and are sensitive to the initial choice of stepsizes.In this paper, we focus on development of PPA-type method for solving (1). Based on LPPA, we propose a general self-adaptive relaxed-PPA method (SRPPA) which is simple yet efficient. Our algorithm capitalizes certain features of PPA, hence, we term it relaxed-PPA. The proposed algorithm has several nice fronts. First, our method is a PPA-type method with asymmetrical linear term, which is clearly a different nature to classical PPA. It relaxes the jointly structure of subproblem to a tractable one. Second, we provide two simple search directions for new iterate. Third, the stepsize choice is flexible, and the condition for convergence guarantee is weaker than both PPA-CM and LPPA. Finally, simple adaptive strategies are employed to choose stepsize, and this appealing aspect is significantly important in practice. We also demonstrate that our method is relevant for various applications whose practical success is made possible by our efficient algorithm.This paper is organized as follows. In Section2, we provide some notations and preliminaries which are useful for subsequent analysis. In Section 3, we review some related works. The general relaxed-PPA and its variant are formally presented in Section 4. Self-adaptive strategies to choose stepsize are also described. Next, in Section 5, the convergence analysis is provided. In Section 6, we present some concrete applications of (1) and elaborate on the implementation of our method; preliminary numerical results are also reported to verify the efficiency of our proposed method. Finally, in Section 7, we conclude the paper with a discussion about the future research directions.


## 2. Preliminaries

In this section, we first establish some important notations used throughout this paper. Then, we describe the variational inequality formulation of (1) which is convenient for the convergence analysis.∂
θ
(
x
) denotes the subdifferential set of the convex function θ(x):
(3)∂θ(x)∶={d∈ℜn∣θ(y)-θ(x)≥dT(y-x),∀y∈ℜn},
and d∈∂θ(x) is called a subgradient of θ(x), see [26]. Let f(x)∈∂θ(x) and f(y)∈∂θ(y), by the convexity of the function θ, we have
(4)(x-y)T(f(x)-f(y))≥0,∀x,y∈ℜn,
which indicates that the mapping f is monotone.Now, we show that (1) can be characterized by a variational inequality; see, for example, [27]. By attaching a Lagrange multiplier vector λ∈Λ to the linear constraint Ax=b(orAx≥b), the Lagrangian function of (1) is
(5)L(x,λ)=θ(x)-λT(Ax-b);
here,
(6)Λ={ℜm,fortheequalityconstraintsAx=b,ℜ+m,fortheinequalityconstraintsAx≥b,
and L(x,λ) is defined on 𝒳×Λ. Then, by the optimality condition, we can easily see that (1) amounts to finding a pair of (x*,λ*) which satisfies
(7)x*∈𝒳,(x-x*)T{f(x*)-ATλ*}≥0,∀x∈𝒳,λ*∈Λ,(λ-λ*)T(Ax*-b)≥0,∀λ∈Λ,
where f(x*)∈∂θ(x*). Denoting
(8)u=(xλ),F(u)=(f(x)-ATλAx-b),Ω=𝒳×Λ,
the system (7) can be characterized by the following variational inequality denoted by VI(Ω,F):
(9)u*∈Ω,(u-u*)TF(u*)≥0,∀u∈Ω.
Recalling the monotonicity of f, it is easy to get that VI(Ω,F) (9) is monotone. Since the solution set of (1) is assumed to be nonempty, the solution set of VI(Ω,F), denoted by Ω*, is also nonempty. Our analysis will be built upon this equivalent VI formulation.


## 3.  The Existing Related Methods

There are basically two lines of research forVI(Ω,F) (9), either deal with it by implicit methods that are in general computationally intractable or concentrate on relaxing it with explicit methods. In this section, we first briefly review the well-known classical PPA and EGM. And then, PPA-CM [23] and LPPA [24] are discussed, which will provide motivation for our general self-adaptive relaxed-PPA.


### 3.1. Classical PPA for the Equivalent Variational Inequality

PPA and its variants are implicit methods which have fast asymptotical convergence rate and produce highly accurate solutions. At each iteration, the subproblem of classical PPA consists of a regularized term, which can be expressed as follows: given any iterateuk=(xk,λk), find u~k=(x~k,λ~k)∈Ω such that
(10)u~k∈Ω,(u-u~k)T{F(u~k)+r(u~k-uk)}≥0,∀u∈Ω.
Then, the update step is taken as follows:
(11)uk+1=uk-γ(uk-u~k),γ∈(0,2).
PPA has a nice convergence property
(12)∥uk+1-u*∥2≤∥uk-u*∥2-γ(2-γ)∥uk-u~k∥2.
Although classical PPA is conceptually appealing, subproblem (10) presents certain computational challenges. More specifically, primal variable x and dual variable λ are tied together, and their subproblems are treated as a joint problem. In most cases, this joint subproblem may be as difficult as the original problem (9). As a result, PPA is “conceptual” rather than implementable. It lacks capability in exploiting potential decomposable/specific structure of subproblem. Variants of classical PPA have been explored in the literature, in order to solve the proximal subproblem (10), inexactly, see, for example, [14, 15, 17, 19]. Unfortunately, inexact variants take expensive computation for obtaining approximative solutions.


### 3.2. The Methods Based on the Simplest Relaxation

To overcome the drawbacks of the classical PPA, it is instinctive to relax subproblem (10) to a solvable one. The most straightforward and simplest relaxation is to replace F(u~k) with F(uk) in the proximal subproblem (10), which amounts to the following subproblem:
(13)u~k∈Ω,(u-u~k)T{F(uk)+r(u~k-uk)}≥0,∀u∈Ω.
The solution of the relaxed problem (13) is given by u~k=PΩ[uk-(1/r)F(uk)]. It is clear that methods with such relaxation are explicit type methods. However, u~k cannot be accepted directly as the new iterate. Using the terminology “predictor-corrector,” such point can be viewed as a predictor. Here, we list two simple methods which employ predictoru~k to obtain corrector as the new iterate. (i)
The extragradient method (EGM) updates the new iterate (corrector) by(14)uk+1=PΩ[uk-1rF(u~k)].(ii)
The projection and contraction methods (PCM) [28–30] perform update as follows:
(15)uk+1=uk-γαk*dkoruk+1=PΩ[uk-γαk*rF(u~k)],
where
(16)dk=(uk-u~k)-1r(F(uk)-F(u~k)),αk*=(uk-u~k)Tdk∥dk∥2.The sequence {uk} generated by the above mentioned EGM or PCM satisfies
(17)∥uk+1-u*∥2≤∥uk-u*∥2-c∥uk-u~k∥2,c>0,
which is similar to PPA. Both EGM and PCM use the simplest relaxation to obtain u~k in kth iteration, hence are computationally practical. These methods have appealing practical aspects; however, such simplest relaxation does not exploit the inner structure of VI(Ω,F) (9). This observation prompts the need for dedicated relaxations.


### 3.3. PPA-Type Contraction Method

The algorithms that are closely related to ours are PPA-CM [23] and LPPA [24]. The PPA-CM obtains the predictor u~k by solving the following subproblem: find (x~k,λ~k)∈Ω such that
(18)u~k∈Ω,(u-u~k)T{F(u~k)+S(u~k-uk)}≥0,∀u∈Ω,
where
(19)S=(rIn-AT-AsIm).
And perform the update
(20)uk+1=uk-γ(uk-u~k),γ∈(0,2).
The framework of LPPA is as follows:
(21)u~k∈Ω,(u-u~k)T{F(u~k)+M(u~k-uk)}≥0,∀u∈Ω,
where
(22)M=(rInAT0sIm).
And the new iterate is defined by
(23)uk+1=uk-γαM(uk-u~k),γ∈(0,2).
Both procedures are simple and can solve subproblem efficiently; but their nice convergence results require a quite restrictive condition, that is; rs>∥ATA∥ in PPA-CM and rs>(1/2)∥ATA∥ in LPPA, respectively. The stepsizes r,s are directly determined by such condition; hence, it is important to estimate ∥ATA∥. Overestimation may lead to poor stepsizes and slow convergence, while underestimation may result in divergence. In addition, they are both quite sensitive to the choice of r,s. To overcome those drawbacks, we propose a general self-adaptive relaxed-PPA, and as mentioned earlier, it can provide improved guarantee for convergence and has potential progress in the choice of stepsize. Furthermore, self-adaptive strategies are designed to improve performance.


## 3.1. Classical PPA for the Equivalent Variational Inequality

PPA and its variants are implicit methods which have fast asymptotical convergence rate and produce highly accurate solutions. At each iteration, the subproblem of classical PPA consists of a regularized term, which can be expressed as follows: given any iterateuk=(xk,λk), find u~k=(x~k,λ~k)∈Ω such that
(10)u~k∈Ω,(u-u~k)T{F(u~k)+r(u~k-uk)}≥0,∀u∈Ω.
Then, the update step is taken as follows:
(11)uk+1=uk-γ(uk-u~k),γ∈(0,2).
PPA has a nice convergence property
(12)∥uk+1-u*∥2≤∥uk-u*∥2-γ(2-γ)∥uk-u~k∥2.
Although classical PPA is conceptually appealing, subproblem (10) presents certain computational challenges. More specifically, primal variable x and dual variable λ are tied together, and their subproblems are treated as a joint problem. In most cases, this joint subproblem may be as difficult as the original problem (9). As a result, PPA is “conceptual” rather than implementable. It lacks capability in exploiting potential decomposable/specific structure of subproblem. Variants of classical PPA have been explored in the literature, in order to solve the proximal subproblem (10), inexactly, see, for example, [14, 15, 17, 19]. Unfortunately, inexact variants take expensive computation for obtaining approximative solutions.


## 3.2. The Methods Based on the Simplest Relaxation

To overcome the drawbacks of the classical PPA, it is instinctive to relax subproblem (10) to a solvable one. The most straightforward and simplest relaxation is to replace F(u~k) with F(uk) in the proximal subproblem (10), which amounts to the following subproblem:
(13)u~k∈Ω,(u-u~k)T{F(uk)+r(u~k-uk)}≥0,∀u∈Ω.
The solution of the relaxed problem (13) is given by u~k=PΩ[uk-(1/r)F(uk)]. It is clear that methods with such relaxation are explicit type methods. However, u~k cannot be accepted directly as the new iterate. Using the terminology “predictor-corrector,” such point can be viewed as a predictor. Here, we list two simple methods which employ predictoru~k to obtain corrector as the new iterate. (i)
The extragradient method (EGM) updates the new iterate (corrector) by(14)uk+1=PΩ[uk-1rF(u~k)].(ii)
The projection and contraction methods (PCM) [28–30] perform update as follows:
(15)uk+1=uk-γαk*dkoruk+1=PΩ[uk-γαk*rF(u~k)],
where
(16)dk=(uk-u~k)-1r(F(uk)-F(u~k)),αk*=(uk-u~k)Tdk∥dk∥2.The sequence {uk} generated by the above mentioned EGM or PCM satisfies
(17)∥uk+1-u*∥2≤∥uk-u*∥2-c∥uk-u~k∥2,c>0,
which is similar to PPA. Both EGM and PCM use the simplest relaxation to obtain u~k in kth iteration, hence are computationally practical. These methods have appealing practical aspects; however, such simplest relaxation does not exploit the inner structure of VI(Ω,F) (9). This observation prompts the need for dedicated relaxations.


## 3.3. PPA-Type Contraction Method

The algorithms that are closely related to ours are PPA-CM [23] and LPPA [24]. The PPA-CM obtains the predictor u~k by solving the following subproblem: find (x~k,λ~k)∈Ω such that
(18)u~k∈Ω,(u-u~k)T{F(u~k)+S(u~k-uk)}≥0,∀u∈Ω,
where
(19)S=(rIn-AT-AsIm).
And perform the update
(20)uk+1=uk-γ(uk-u~k),γ∈(0,2).
The framework of LPPA is as follows:
(21)u~k∈Ω,(u-u~k)T{F(u~k)+M(u~k-uk)}≥0,∀u∈Ω,
where
(22)M=(rInAT0sIm).
And the new iterate is defined by
(23)uk+1=uk-γαM(uk-u~k),γ∈(0,2).
Both procedures are simple and can solve subproblem efficiently; but their nice convergence results require a quite restrictive condition, that is; rs>∥ATA∥ in PPA-CM and rs>(1/2)∥ATA∥ in LPPA, respectively. The stepsizes r,s are directly determined by such condition; hence, it is important to estimate ∥ATA∥. Overestimation may lead to poor stepsizes and slow convergence, while underestimation may result in divergence. In addition, they are both quite sensitive to the choice of r,s. To overcome those drawbacks, we propose a general self-adaptive relaxed-PPA, and as mentioned earlier, it can provide improved guarantee for convergence and has potential progress in the choice of stepsize. Furthermore, self-adaptive strategies are designed to improve performance.


## 4. The General Self-Adaptive Relaxed PPA-Method

In this section, we weave together the ideas of the previous section to present general self-adaptive relaxed-PPA method (SRPPA) which is mostly inspired by LPPA [24]. At first sight, the predictor applied in SRPPA is much the same as LPPA, but the stepsize choice condition for convergence is quite different; moreover, we prove that it is weaker than LPPA. Self-adaptive strategies are elaborately designed to ensure the robustness of our algorithm. Two simple yet efficient constructions of new iterate are also presented which will provide some inspirations for designing various search directions.


### 4.1. General Relaxed-PPA Method

The general primal-dual relaxed-PPA method with implementable structure for (1) is summarized in Algorithm 1. Note that the order of x and λ can be changed to obtain a variant, which is summarized in Algorithm 2. Our relaxed-PPA is intended to blend the implementable properties of EGM (or PCM) with the fast convergence performance of PPA. Now, it is helpful to introduce additional notations that will be used in the rest of this paper. Let G be a positive symmetry definite matrix (we will specify it later),
(24)M=(rInAT0sIm),H=(rIn00sIm),(25)d(uk,u~k)=G-1M(uk-u~k),(26)φ(uk,u~k)=(uk-u~k)TM(uk-u~k).
The relaxed-PPA described here involves two steps. First, we solve the relaxed subproblem (∗), (∗∗) to obtain predictor, which is nice and efficient for the nature of the problem under study. Note that the x-predictor step (∗) involves minimizing θ plus a convex quadratic function, and under Assumption 2, it can be efficiently solved or it admits a closed form solution. And then, λ-predictor step (∗∗) is just a projection onto Λ which is tractable and computationally efficient. It is clear that the prediction step employs a Gauss-Seidel manner to update information efficiently. The correction step ($) only involves matrix-vector multiplication which is very simple and straightforward.Algorithm 1:General primal-dual relaxed-PPA method.
Step  1. Initialization. Let γ∈(0,2) and pick (x0,λ0)∈ℜn×Λ, set k=0.
Step  2. Predictor. Let
x~k=Argmin{θ(x)+r2∥x-[xk+1rATλk]∥2∣x∈𝒳},     (*)
and
λ~k=PΛ[λk-1s(Ax~k-b)].                (**)
Step  3. If the stepsizes r and s are chosen satisfying
φ(uk,u~k)≥14∥d(uk,u~k)∥G2,                  (#)
then go to Step 4. Otherwise, increaser,s and go back to Step 2.
Step  4. Corrector and updating
uk+1=uk-αkG-1M(uk-u~k).                ($)
Step  5. Adjustment
(r,s)={(r2,s2),ifφ(uk,u~k)≥κ*∥d(uk,u~k)∥G2;(r,s),otherwise.Hereκ>4.
Setk:=k+1.Algorithm 2:General dual-primal relaxed-PPA method.
Step  1. Initialization. Let γ∈(0,2) and pick (x0,λ0)∈ℜn×Λ, set k=0.
Step  2. Predictor. Let
λ~k=PΛ[λk-1s(Axk-b)],                 (†)
and
x~k=Argmin{θ(x)+r2∥x-[xk+1rATλ~k]∥2∣x∈𝒳}.        (††)
Step  3. If the stepsize r and s are chosen satisfying
φ(uk,u~k)≥14∥d(uk,u~k)∥G2,                    (‡)
then go to Step 4. Otherwise, increaser,s and go back to Step 2.
Step  4. Corrector and updating
uk+1=uk-αkG-1M(uk-u~k).                (∧)
Step  5. Adjustment
(r,s)={(r2,s2),ifφ(uk,u~k)≥κ*∥d(uk,u~k)∥G2;(r,s),otherwise.Hereκ>4.
Setk:=k+1.Remark 3.
We first make some insight into the correction step in Algorithm1. The obtained  u~k plays no direct role as the new iterate. Instead, we need some kind of  “corrector” defined in ($). Although  matrix G in ($) is just a required positive symmetry definite, our goal here is to fully integrate the information of uk and u~k to construct effective, simple search direction G-1M(uk-u~k) for the corrector. Based on this consideration, we elaborately provide two simple choices of G.Case
1. It is natural to set G=H to induce a simple update form. Then, it is easy to get that
(27)(xk+1λk+1)=(xkλk)-αk(xk-x~k+1rAT(λk-λ~k)λk-λ~k).Case
2. Let G=MH-1MT. This case is a little less intuitive, but it can lead to a simple update form as well as Case 1. The underlying derivation is a little more complicate. Applying G in ($), we get
(28)uk+1=uk-αkM-THM-1M(uk-u~k)=uk-αkM-TH(uk-u~k).RecallingM is a lower triangular matrix, by the fact that its inverse is also a lower triangular, we have
(29)G-1M=M-TH=(In0-AsIm).Plugging the previous relationship to (31), we have
(30)(xk+1λk+1)=(xkλk)-αk(xk-x~kλk-λ~k-1sA(xk-x~k)).In fact, this is a scheme of Gaussian back substitution.Both cases only involve one matrix-vector multiplication which makes the update form simple. And the computational cost is usually inexpensive.Remark 4.
We now study the subproblem described in (∗), (∗∗), and the stepsize choice condition (#). For easy analysis, we characterize (∗), (∗∗) as the following VI formulation.
Find(x~k,λ~k)∈Ω such that
(31)(x-x~kλ-λ~k)T{(f(x~k)-ATλ~kAx~k-b)+(rInAT0sIm)(x~k-xkλ~k-λk)}≥0,∀(xλ)∈Ω,
and its compact form
(32)u~k∈Ω,(u-u~k)T{F(u~k)+M(u~k-uk)}≥0,∀u∈Ω.
We observe that subproblem (32) is similar to (10) in PPA, except for the construction of asymmetry matrix M. As mentioned before, (32) is the same as the prediction subproblem in [24]. Even though they are closely related, the stepsize choice here is quite different. We provide more specific and weaker condition for stepsize r,s. It is clear that condition (#) does not need prior knowledge of matrix A. Furthermore, it only involves matrix-vector multiplication, and so, it is easy to verify, and it is amenable to large-scale A. If r,s fail to meet this convergence condition (#), one should appropriately increase r,s. In the following subsection, we will elaborate on the self-adaptive strategies to increase the stepsizes. At this point, condition (#) may be seen somewhat unmotivated. Some insight into this will be provided later, as we proceed with the convergence analysis. The convergence condition in [24] has a quite different feature: r,s satisfy
(33)rs≥12∥ATA∥.
It is stronger than condition (#). The following lemma is devoted to the proof of this result.Lemma 5.
Let{uk} be the sequence generated by Algorithm 1, H, d(uk,u~k), and φ(uk,u~k)  defined in (24), (25), and (26), respectively. Suppose that condition (33) is satisfied. Then, condition (#) holds immediately.Proof.
Note that(34)φ(uk,u~k)=r∥xk-x~k∥2+s∥λk-λ~k∥2+(λk-λ~k)TA(xk-x~k).
Recall that matrix G described in Algorithm 1 can be designed in two different cases.
Case  1. If  G=H, we immediately have
(35)∥d(uk,u~k)∥G2=∥uk-u~k∥H2+1r∥AT(λk-λ~k)∥2+2(λk-λ~k)TA(xk-x~k).
According to Cauchy-Schwarz inequality, we get(36)φ(uk,u~k)-14∥d(uk,u~k)∥G2=34∥uk-u~k∥H2-14r∥AT(λk-λ~k)∥2+12(λk-λ~k)TA(xk-x~k)≥34∥uk-u~k∥H2-14r∥AT(λk-λ~k)∥2-14(2r∥xk-x~k∥2+12r∥AT(λk-λ~k)∥2)=r4∥xk-x~k∥2+3s4∥λk-λ~k∥2-38r∥AT(λk-λ~k)∥2≥r4∥xk-x~k∥2.
The last inequality follows directly from condition (33).
Case  2. If  G=MH-1MT, by the definition of d(uk,u~k), we obtain
(37)∥d(uk,u~k)∥G2=(uk-u~k)TMTG-1M(uk-u~k)=∥uk-u~k∥H2.
Then, we get(38)φ(uk,u~k)-14∥d(uk,u~k)∥G2=34∥uk-u~k∥H2+(λk-λ~k)TA(xk-x~k)≥34∥uk-u~k∥H2-124r3∥xk-x~k∥2-1234r∥AT(λk-λ~k)∥2=r12∥xk-x~k∥2+3s4∥λk-λ~k∥2-38r∥AT(λk-λ~k)∥2≥r12∥xk-x~k∥2.
The first inequality follows from the Cauchy-Schwarz inequality, and the last one follows directly from condition (33).
Note that, in both cases, we have that condition (#) holds ifrs≥(1/2)∥ATA∥.Condition (33) is not only stronger than Condition (#), but it also requires that matrix M is positive semidefinite, while condition (#) does not. Furthermore, condition (33) may require the explicit expression of A or knowledge of ∥ATA∥. Despite these drawbacks, condition (33) is appealing to the problems in which ∥ATA∥ is known beforehand or easy to compute/obtain. For instance, A is small scale, an identity matrix or a projection operator, and so forth. It is clear that both condition (#) and (33) are more flexible than the one in PPA-CM [23]. The most aggressive condition (#) may lead to further improvement in stepsize choice. Moreover, it is worthwhile to notice that condition (#) is elegantly designed and provides φ(uk,u~k) with favourable property. In fact, for general matrix G, condition (#) also can guarantee convergence.Remark 6.
The update stepsizeαk plays an important role here. In fact, it can be regarded as an optimal stepsize which will be illustrated in the following section.Remark 7.
We should restrict the adjustment in Step 5 of Algorithm1 within a limited number to avoid divergence.In Algorithm1, we carry out the x-predictor before performing λ-predictor. The roles of x and λ are symmetric; hence, sweeping the order will not break the Gauss-Seidel structure. We switch x and λ and obtain a variant of relaxed-PPA with the order of the x-predictor step and λ-predictor step reversed. This variant is illustrated in Algorithm 2. However, there is no a priori information to know which algorithm is superior. Here, we let
(39)M=(rIn0-AsIm).


### 4.2. Adaptive Enhancements

Both PPA-CM and LPPA employ fixed stepsizesr,s. Experiments reveal that they will suffer slow convergence when stepsizes r,s are chosen inappropriately. A natural question is, how to choose the proper initial stepsizes r,s. Here, we propose self-adaptive strategies with the goal of improving the convergence in practice, as well as making performance less dependent on the initial choice of stepsizes. Our strategies dynamically incorporate the information of the current iteration to perform more informative choice of stepsizes for the next iteration [31]. When doing so, the algorithm will be adaptive and free from the initial choice. Denote
(40)(dxkdλk)=((xk-x~k)+1rAT(λk-λ~k)λk-λ~k),
and then, (31) can be rewritten as
(41)(x-x~kλ-λ~k)T{(f(x~k)-ATλ~kAx~k-b)-(rIn00sIm)(dxkdλk)}≥0,∀(xλ)∈Ω.
Under H-norm, the quantity dxk can be viewed as a residual for the dual feasibility condition, and dλk can be viewed as a primal residual. These two residuals converge to zero as relaxed-PPA proceeds. Note that
(42)∥dxk∥r2=r∥xk-x~k∥2+2(xk-x~k)TAT(λk-λ~k)+1r∥AT(λk-λ~k)∥2,∥dλk∥s2=s∥λk-λ~k∥2.
And this implies that small values of s tend to reduce the primal residual but have potential to enlarge violations of dual feasibility and tend to produce larger dual residual. This observation motivates us to balance primal and dual residuals. When condition (#) fails, we increase stepsizes r,s properly according to the  adaptation shown in Algorithm 3.Algorithm 3:Adaptation-I.
If
∥
d
x
k
∥
r
2
≥
τ
1
∥
d
λ
k
∥
s
2
r:=r;s:=s*2;
else if
τ
2
∥
d
x
k
∥
r
2
≤
∥
d
λ
k
∥
s
2
r:=r*2;s:=s;
else
r:=r*1.5;s:=s*1.5.Here,τ1>1, τ2>1. This adaptation strategy increases s when the dual residual ∥dxk∥r2 appears large compared to the primal residual ∥dλk∥s2 and increases r when the dual residual ∥dxk∥r2  seems too small relative to the primal residual  ∥dλk∥s2.As mentioned, condition (33) is stronger than condition (#). If one chooses condition (33), our RPPA also converges. It must have predetermined stepsizes satisfying rs=μ∥ATA∥ (here, μ≥0.5). However, there is no priority knowledge of the choice of individual r or s. Here, we can also adjust r,s automatically when choosing condition (33). Intuitively, we consider expansion of the entire residual under H-norm:
(43)∥dxk∥r2+∥dλk∥s2=r∥xk-x~k∥2+s∥λk-λ~k∥2+1r∥AT(λk-λ~k)∥2+2(xk-x~k)TAT(λk-λ~k);
there are three terms involving r or s, and we intend to balance these terms. For fixed μ, take s=(μ/r)∥AAT∥; then
(44)s∥λk-λ~k∥2=μr∥AAT∥∥λk-λ~k∥2.
Applying (44) into (43), clearly we have
(45)∥dxk∥r2+∥dλk∥s2=r∥xk-x~k∥2+1r(μ∥AAT∥∥λk-λ~k∥2+∥AT(λk-λ~k)∥2)+2(xk-x~k)TAT(λk-λ~k).
Now, we consider adjusting stepsize to balance r∥xk-x~k∥2 and (1/r)(μ∥AAT∥∥λk-λ~k∥2+∥AT(λk-λ~k)∥2) and obtain another adaptation  strategy (see Algorithm 4).Algorithm 4:Adaptation-II.
if
r
∥
x
k
-
x
~
k
∥
2
≥
τ
1
(
s
∥
λ
k
-
λ
~
k
∥
2
+
1
r
∥
A
T
(
λ
k
-
λ
~
k
)
∥
2
)
r:=r2;s=μ∥AAT∥r;
else if
τ
2
r
∥
x
k
-
x
~
k
∥
2
≤
(
s
∥
λ
k
-
λ
~
k
∥
2
+
1
r
∥
A
T
(
λ
k
-
λ
~
k
)
∥
2
)
s:=s2;r=μ∥AAT∥s;
else
r:=r;s:=s.It is worth noting that too many adjustments of stepsizes by Algorithm4 might cause the algorithm to diverge in practice, and we therefore restrict these adaptations  within a limited number of iterations. If one chooses Algorithm 4, there is no need to carry out Step 5 in Algorithm 1 (or Algorithm 2). These techniques embedded into relaxed-PPA automatically choose a “better” stepsize for the next iteration.


## 4.1. General Relaxed-PPA Method

The general primal-dual relaxed-PPA method with implementable structure for (1) is summarized in Algorithm 1. Note that the order of x and λ can be changed to obtain a variant, which is summarized in Algorithm 2. Our relaxed-PPA is intended to blend the implementable properties of EGM (or PCM) with the fast convergence performance of PPA. Now, it is helpful to introduce additional notations that will be used in the rest of this paper. Let G be a positive symmetry definite matrix (we will specify it later),
(24)M=(rInAT0sIm),H=(rIn00sIm),(25)d(uk,u~k)=G-1M(uk-u~k),(26)φ(uk,u~k)=(uk-u~k)TM(uk-u~k).
The relaxed-PPA described here involves two steps. First, we solve the relaxed subproblem (∗), (∗∗) to obtain predictor, which is nice and efficient for the nature of the problem under study. Note that the x-predictor step (∗) involves minimizing θ plus a convex quadratic function, and under Assumption 2, it can be efficiently solved or it admits a closed form solution. And then, λ-predictor step (∗∗) is just a projection onto Λ which is tractable and computationally efficient. It is clear that the prediction step employs a Gauss-Seidel manner to update information efficiently. The correction step ($) only involves matrix-vector multiplication which is very simple and straightforward.Algorithm 1:General primal-dual relaxed-PPA method.
Step  1. Initialization. Let γ∈(0,2) and pick (x0,λ0)∈ℜn×Λ, set k=0.
Step  2. Predictor. Let
x~k=Argmin{θ(x)+r2∥x-[xk+1rATλk]∥2∣x∈𝒳},     (*)
and
λ~k=PΛ[λk-1s(Ax~k-b)].                (**)
Step  3. If the stepsizes r and s are chosen satisfying
φ(uk,u~k)≥14∥d(uk,u~k)∥G2,                  (#)
then go to Step 4. Otherwise, increaser,s and go back to Step 2.
Step  4. Corrector and updating
uk+1=uk-αkG-1M(uk-u~k).                ($)
Step  5. Adjustment
(r,s)={(r2,s2),ifφ(uk,u~k)≥κ*∥d(uk,u~k)∥G2;(r,s),otherwise.Hereκ>4.
Setk:=k+1.Algorithm 2:General dual-primal relaxed-PPA method.
Step  1. Initialization. Let γ∈(0,2) and pick (x0,λ0)∈ℜn×Λ, set k=0.
Step  2. Predictor. Let
λ~k=PΛ[λk-1s(Axk-b)],                 (†)
and
x~k=Argmin{θ(x)+r2∥x-[xk+1rATλ~k]∥2∣x∈𝒳}.        (††)
Step  3. If the stepsize r and s are chosen satisfying
φ(uk,u~k)≥14∥d(uk,u~k)∥G2,                    (‡)
then go to Step 4. Otherwise, increaser,s and go back to Step 2.
Step  4. Corrector and updating
uk+1=uk-αkG-1M(uk-u~k).                (∧)
Step  5. Adjustment
(r,s)={(r2,s2),ifφ(uk,u~k)≥κ*∥d(uk,u~k)∥G2;(r,s),otherwise.Hereκ>4.
Setk:=k+1.Remark 3.
We first make some insight into the correction step in Algorithm1. The obtained  u~k plays no direct role as the new iterate. Instead, we need some kind of  “corrector” defined in ($). Although  matrix G in ($) is just a required positive symmetry definite, our goal here is to fully integrate the information of uk and u~k to construct effective, simple search direction G-1M(uk-u~k) for the corrector. Based on this consideration, we elaborately provide two simple choices of G.Case
1. It is natural to set G=H to induce a simple update form. Then, it is easy to get that
(27)(xk+1λk+1)=(xkλk)-αk(xk-x~k+1rAT(λk-λ~k)λk-λ~k).Case
2. Let G=MH-1MT. This case is a little less intuitive, but it can lead to a simple update form as well as Case 1. The underlying derivation is a little more complicate. Applying G in ($), we get
(28)uk+1=uk-αkM-THM-1M(uk-u~k)=uk-αkM-TH(uk-u~k).RecallingM is a lower triangular matrix, by the fact that its inverse is also a lower triangular, we have
(29)G-1M=M-TH=(In0-AsIm).Plugging the previous relationship to (31), we have
(30)(xk+1λk+1)=(xkλk)-αk(xk-x~kλk-λ~k-1sA(xk-x~k)).In fact, this is a scheme of Gaussian back substitution.Both cases only involve one matrix-vector multiplication which makes the update form simple. And the computational cost is usually inexpensive.Remark 4.
We now study the subproblem described in (∗), (∗∗), and the stepsize choice condition (#). For easy analysis, we characterize (∗), (∗∗) as the following VI formulation.
Find(x~k,λ~k)∈Ω such that
(31)(x-x~kλ-λ~k)T{(f(x~k)-ATλ~kAx~k-b)+(rInAT0sIm)(x~k-xkλ~k-λk)}≥0,∀(xλ)∈Ω,
and its compact form
(32)u~k∈Ω,(u-u~k)T{F(u~k)+M(u~k-uk)}≥0,∀u∈Ω.
We observe that subproblem (32) is similar to (10) in PPA, except for the construction of asymmetry matrix M. As mentioned before, (32) is the same as the prediction subproblem in [24]. Even though they are closely related, the stepsize choice here is quite different. We provide more specific and weaker condition for stepsize r,s. It is clear that condition (#) does not need prior knowledge of matrix A. Furthermore, it only involves matrix-vector multiplication, and so, it is easy to verify, and it is amenable to large-scale A. If r,s fail to meet this convergence condition (#), one should appropriately increase r,s. In the following subsection, we will elaborate on the self-adaptive strategies to increase the stepsizes. At this point, condition (#) may be seen somewhat unmotivated. Some insight into this will be provided later, as we proceed with the convergence analysis. The convergence condition in [24] has a quite different feature: r,s satisfy
(33)rs≥12∥ATA∥.
It is stronger than condition (#). The following lemma is devoted to the proof of this result.Lemma 5.
Let{uk} be the sequence generated by Algorithm 1, H, d(uk,u~k), and φ(uk,u~k)  defined in (24), (25), and (26), respectively. Suppose that condition (33) is satisfied. Then, condition (#) holds immediately.Proof.
Note that(34)φ(uk,u~k)=r∥xk-x~k∥2+s∥λk-λ~k∥2+(λk-λ~k)TA(xk-x~k).
Recall that matrix G described in Algorithm 1 can be designed in two different cases.
Case  1. If  G=H, we immediately have
(35)∥d(uk,u~k)∥G2=∥uk-u~k∥H2+1r∥AT(λk-λ~k)∥2+2(λk-λ~k)TA(xk-x~k).
According to Cauchy-Schwarz inequality, we get(36)φ(uk,u~k)-14∥d(uk,u~k)∥G2=34∥uk-u~k∥H2-14r∥AT(λk-λ~k)∥2+12(λk-λ~k)TA(xk-x~k)≥34∥uk-u~k∥H2-14r∥AT(λk-λ~k)∥2-14(2r∥xk-x~k∥2+12r∥AT(λk-λ~k)∥2)=r4∥xk-x~k∥2+3s4∥λk-λ~k∥2-38r∥AT(λk-λ~k)∥2≥r4∥xk-x~k∥2.
The last inequality follows directly from condition (33).
Case  2. If  G=MH-1MT, by the definition of d(uk,u~k), we obtain
(37)∥d(uk,u~k)∥G2=(uk-u~k)TMTG-1M(uk-u~k)=∥uk-u~k∥H2.
Then, we get(38)φ(uk,u~k)-14∥d(uk,u~k)∥G2=34∥uk-u~k∥H2+(λk-λ~k)TA(xk-x~k)≥34∥uk-u~k∥H2-124r3∥xk-x~k∥2-1234r∥AT(λk-λ~k)∥2=r12∥xk-x~k∥2+3s4∥λk-λ~k∥2-38r∥AT(λk-λ~k)∥2≥r12∥xk-x~k∥2.
The first inequality follows from the Cauchy-Schwarz inequality, and the last one follows directly from condition (33).
Note that, in both cases, we have that condition (#) holds ifrs≥(1/2)∥ATA∥.Condition (33) is not only stronger than Condition (#), but it also requires that matrix M is positive semidefinite, while condition (#) does not. Furthermore, condition (33) may require the explicit expression of A or knowledge of ∥ATA∥. Despite these drawbacks, condition (33) is appealing to the problems in which ∥ATA∥ is known beforehand or easy to compute/obtain. For instance, A is small scale, an identity matrix or a projection operator, and so forth. It is clear that both condition (#) and (33) are more flexible than the one in PPA-CM [23]. The most aggressive condition (#) may lead to further improvement in stepsize choice. Moreover, it is worthwhile to notice that condition (#) is elegantly designed and provides φ(uk,u~k) with favourable property. In fact, for general matrix G, condition (#) also can guarantee convergence.Remark 6.
The update stepsizeαk plays an important role here. In fact, it can be regarded as an optimal stepsize which will be illustrated in the following section.Remark 7.
We should restrict the adjustment in Step 5 of Algorithm1 within a limited number to avoid divergence.In Algorithm1, we carry out the x-predictor before performing λ-predictor. The roles of x and λ are symmetric; hence, sweeping the order will not break the Gauss-Seidel structure. We switch x and λ and obtain a variant of relaxed-PPA with the order of the x-predictor step and λ-predictor step reversed. This variant is illustrated in Algorithm 2. However, there is no a priori information to know which algorithm is superior. Here, we let
(39)M=(rIn0-AsIm).


## 4.2. Adaptive Enhancements

Both PPA-CM and LPPA employ fixed stepsizesr,s. Experiments reveal that they will suffer slow convergence when stepsizes r,s are chosen inappropriately. A natural question is, how to choose the proper initial stepsizes r,s. Here, we propose self-adaptive strategies with the goal of improving the convergence in practice, as well as making performance less dependent on the initial choice of stepsizes. Our strategies dynamically incorporate the information of the current iteration to perform more informative choice of stepsizes for the next iteration [31]. When doing so, the algorithm will be adaptive and free from the initial choice. Denote
(40)(dxkdλk)=((xk-x~k)+1rAT(λk-λ~k)λk-λ~k),
and then, (31) can be rewritten as
(41)(x-x~kλ-λ~k)T{(f(x~k)-ATλ~kAx~k-b)-(rIn00sIm)(dxkdλk)}≥0,∀(xλ)∈Ω.
Under H-norm, the quantity dxk can be viewed as a residual for the dual feasibility condition, and dλk can be viewed as a primal residual. These two residuals converge to zero as relaxed-PPA proceeds. Note that
(42)∥dxk∥r2=r∥xk-x~k∥2+2(xk-x~k)TAT(λk-λ~k)+1r∥AT(λk-λ~k)∥2,∥dλk∥s2=s∥λk-λ~k∥2.
And this implies that small values of s tend to reduce the primal residual but have potential to enlarge violations of dual feasibility and tend to produce larger dual residual. This observation motivates us to balance primal and dual residuals. When condition (#) fails, we increase stepsizes r,s properly according to the  adaptation shown in Algorithm 3.Algorithm 3:Adaptation-I.
If
∥
d
x
k
∥
r
2
≥
τ
1
∥
d
λ
k
∥
s
2
r:=r;s:=s*2;
else if
τ
2
∥
d
x
k
∥
r
2
≤
∥
d
λ
k
∥
s
2
r:=r*2;s:=s;
else
r:=r*1.5;s:=s*1.5.Here,τ1>1, τ2>1. This adaptation strategy increases s when the dual residual ∥dxk∥r2 appears large compared to the primal residual ∥dλk∥s2 and increases r when the dual residual ∥dxk∥r2  seems too small relative to the primal residual  ∥dλk∥s2.As mentioned, condition (33) is stronger than condition (#). If one chooses condition (33), our RPPA also converges. It must have predetermined stepsizes satisfying rs=μ∥ATA∥ (here, μ≥0.5). However, there is no priority knowledge of the choice of individual r or s. Here, we can also adjust r,s automatically when choosing condition (33). Intuitively, we consider expansion of the entire residual under H-norm:
(43)∥dxk∥r2+∥dλk∥s2=r∥xk-x~k∥2+s∥λk-λ~k∥2+1r∥AT(λk-λ~k)∥2+2(xk-x~k)TAT(λk-λ~k);
there are three terms involving r or s, and we intend to balance these terms. For fixed μ, take s=(μ/r)∥AAT∥; then
(44)s∥λk-λ~k∥2=μr∥AAT∥∥λk-λ~k∥2.
Applying (44) into (43), clearly we have
(45)∥dxk∥r2+∥dλk∥s2=r∥xk-x~k∥2+1r(μ∥AAT∥∥λk-λ~k∥2+∥AT(λk-λ~k)∥2)+2(xk-x~k)TAT(λk-λ~k).
Now, we consider adjusting stepsize to balance r∥xk-x~k∥2 and (1/r)(μ∥AAT∥∥λk-λ~k∥2+∥AT(λk-λ~k)∥2) and obtain another adaptation  strategy (see Algorithm 4).Algorithm 4:Adaptation-II.
if
r
∥
x
k
-
x
~
k
∥
2
≥
τ
1
(
s
∥
λ
k
-
λ
~
k
∥
2
+
1
r
∥
A
T
(
λ
k
-
λ
~
k
)
∥
2
)
r:=r2;s=μ∥AAT∥r;
else if
τ
2
r
∥
x
k
-
x
~
k
∥
2
≤
(
s
∥
λ
k
-
λ
~
k
∥
2
+
1
r
∥
A
T
(
λ
k
-
λ
~
k
)
∥
2
)
s:=s2;r=μ∥AAT∥s;
else
r:=r;s:=s.It is worth noting that too many adjustments of stepsizes by Algorithm4 might cause the algorithm to diverge in practice, and we therefore restrict these adaptations  within a limited number of iterations. If one chooses Algorithm 4, there is no need to carry out Step 5 in Algorithm 1 (or Algorithm 2). These techniques embedded into relaxed-PPA automatically choose a “better” stepsize for the next iteration.


## 5. Convergence Analysis

In this section, we analyze convergence of our primal-dual relaxed-PPA. The convergence analysis of dual-primal scheme can follow a similar procedure.Letu*=(x*,λ*) be any solution point, setting u=u* in (32) yields
(46)(u~k-u*)TM(uk-u~k)≥(u~k-u*)TF(u~k).
Since u~k∈Ω, we have (u~k-u*)TF(u*)≥0. Consequently, by using the monotonicity of F, the right hand side of (46) is nonnegative, and thus
(47)(u~k-u*)TM(uk-u~k)≥0.
Now, we are writing the update as
(48)u(α)=uk-αd(uk,u~k)
and
(49)ϑ(α)∶=∥uk-u*∥G2-∥u(α)-u*∥G2,q(α)=2α(uk-u~k)TGd(uk,u~k)-α2∥d(uk,u~k)∥G2.
Here, ϑ(α) can be viewed as a progress function. The following lemma shows that q(α) is a lower bound of ϑ(α).Lemma 8.
Letϑ(α) and q(α) be defined in (49); then one has
(50)ϑ(α)≥q(α).Proof.
Letu* be any solution, from the definition of u(α), we have
(51)ϑ(α)=∥uk-u*∥G2-∥(uk-u*)-αd(uk,u~k)∥G2=2α(uk-u*)TGd(uk,u~k)-α2∥d(uk,u~k)∥G2.
Applying (47) to the first term of (51) gives
(52)(uk-u*)TGd(uk,u~k)=(uk-u~k)TGd(uk,u~k)+(u~k-u*)TGd(uk,u~k)≥(uk-u~k)TGd(uk,u~k).
Substituting (52) into (51), we immediately obtain the assertion.We note thatq(α) is a quadratic function of α and it is natural to maximize q(α) to obtain an “optimal” α:
(53)αk*=(uk-u~k)TGd(uk,u~k)∥d(uk,u~k)∥G2.
We now show that the “optimal” αk* is bounded above from zero in the following Lemma.Lemma 9.
Let sequence{uk} be produced by Algorithm 1, αk*  defined in (53); then, one has
(54)αk*≥14>0.Proof.
Using the definition ofαk* in (53), we have, for all k,
(55)αk*=φ(uk,u~k)∥d(uk,u~k)∥G2≥14.
The inequality follows from  condition (#).Settingα=αk=αk*γ in (50) yields
(56)∥uk+1-u*∥G2≤∥uk-u*∥G2-γ(2-γ)αk*(uk-u~k)TGd(uk,u~k).
Combining Lemmas 8 and 9, we immediately obtain the following convergence theorem.Theorem 10.
Let sequence{uk} be produced by Algorithm 1; then one gets
(57)∥uk+1-u*∥G2≤∥uk-u*∥G2-γ(2-γ)16∥d(uk,u~k)∥G2.Theorem 11.
Let sequence{uk} be generated by Algorithm 1. Then, {uk} converges to some u∞ which is a solution of  VI(Ω,F)  (9).Proof.
First, for eachu∈Ω, we have
(58)(u-u~k)TF(u~k)≥(u-u~k)TM(uk-u~k).
It follows from (57) that {uk} is a bounded sequence and
(59)limk→∞∥uk-u~k∥G=0.
Consequently, {u~k} is also bounded. Since limk→∞∥uk-u~k∥G=0, it follows from (58) that
(60)limk→∞(u-u~k)TF(u~k)≥0,∀u∈Ω.
Because {u~k} is bounded, it has at least a cluster point. Let u∞ be a cluster point of {u~k} and let the subsequence {u~kj} converge to u∞. It follows that
(61)limj→∞(u-u~kj)TF(u~kj)≥0,∀u∈Ω,
and consequently,
(62)(u-u∞)TF(u∞)≥0.∀u∈Ω.
This means that u∞ is a solution of VI(Ω,F). Note that the inequality (57) is true for all solution points of VI(Ω,F), and hence, we have
(63)∥uk+1-u∞∥G2≤∥uk-u∞∥G2,∀k≥0.
Since u~kj→u∞(j→∞) and uk-u~k→0(k→∞), for any given ε>0, there exists an integer l>0 such that
(64)∥u~kl-u∞∥G<ε2,∥ukl-u~kl∥G<ε2.
Therefore, for any k≥kl, it follows from (63) and (64) that
(65)∥uk-u∞∥G≤∥ukl-u∞∥G≤∥ukl-u~kl∥G+∥u~kl-u∞∥G≤ε.
This implies that the sequence {uk} converges to u∞, which is a solution of VI(Ω,F) (9).


## 6.  Applications and Preliminary Numerical Experiments

The general self-adaptive relaxed-PPA (SRPPA) offers a flexible framework for solving many interesting problems. We illustrate our algorithm with different applications: basis pursuit problem, nearest correlation matrix problem. In this section, we describe the results of experiments whose goal is to demonstrate the efficiency of general relaxed-PPA (RPPA) and its self-adaptive version. To that end, we compare RPPA with certain state-of-the-art algorithms on different problems. Our experiments focus on efficiency and speed of convergence and evaluate the methods in terms of their number of iterations and computational times.All the codes were written by Matlab R2009b version, and all the numerical experiments were performed on a Lenovo desktop computer with Intel (R) Core (TM) i5 CPU with 3.2 GHz and 3.5 GB RAM.


### 6.1.  Basis Pursuit Problem

Basis pursuit (BP) finds signal representations in overcomplete dictionaries by equality-constrainedl1 minimization problem. Formally, one solves the problem
(66)min{∥x∥1Ax=b,x∈ℜn}.
And here, ∥·∥1 denotes the l1 norm defined as ∥x∥1:=∑i=1n‍|xi|. BP is a fundamental problem in image processing and modern statistical signal processing, particularly the theory of compressed sensing; see, for example, [1–4] for intensive study. We now discuss our approach to BP problem of over-complete representations. Our experiments in this subsection use synthetic data which were mainly designed to illustrate the nice performance of our RPPA. The synthetic problem that we test here is similar to the one employed in [32]. We generate the data as follows: matrix A is a random m×n matrix, with Gaussian i.i.d. entries of zero mean and variance 1, with m=n/2. xoriginal∈Rn is the original sparse signal, constructed with m/5 nonzero values, randomly from standard normal distribution. We use xoriginal to generate the measurements as b=Axoriginal. It is desirable to use test problems that have a precisely known solution. In fact, when xoriginal is very sparse, it is the solution to the minimization problem (66). Hence, in our synthetic problem, xoriginal is exactly the solution.In our first experiment, we compared general RPPA using two differentG’s mentioned in Section 4.1. For BP problem, we use condition (#) and Algorithm 3. Since A constructed here is a general random matrix, and when A is large scale, ∥ATA∥ might be obtained costly. A simple stopping criterion
(67)err.=∥xk-xoriginal∥≤Tol
was used in this experiment, and the stopping tolerance Tol was set to 10-10. In all the tests, initial stepsizes were set as s=10, r=1, the primal variable x0 was initialized as zeros(n,1), and the dual λ0 was ones(m,1) in Matlab. Table 1 summarizes the performance of general SRPPA. Here, SPDRPPAGi-I(SPDRPPAGi-II) denotes self-adaptive primal-dual RPPA with Algorithm 3 (Algorithm 4), G=H, if i=1, and G=MH-1MT, if i=2. DPRPPA (DPRPPA) denotes dual-primal RPPA version.Table 1
Comparison of behaviours of four SRPPAs.
A
SPDRPPAG1-I
SDPRPPAG1-I
SPDRPPAG2-I
SDPRPPAG2-I
n
It.
CPU
It.
CPU
It.
CPU
It.
CPU
500
360
0.13
241
0.13
513
0.19
875
0.26
1500
379
1.32
259
0.89
387
1.72
594
2.04
2500
501
4.39
343
3.01
429
4.93
648
5.72
3500
399
6.65
325
5.57
663
14.44
1089
18.02
4500
389
12.13
457
14.47
493
20.61
991
31.26
5500
515
24.21
462
21.97
497
30.88
813
38.43
6500
449
25.01
423
23.41
456
33.42
932
51.62
7500
520
45.72
473
41.79
724
63.06
886
77.91Basically, SRPPAs converge very quickly and achieved tight error10-10 in a few hundred iterations. For this experiment, one can see that SDPRPPAG1-I is fastest in all cases. Both SDPRPPAG1-I and SPDRPPAG1-I are Gaussian type methods, with G=H, and they exhibit very similar performance. SDPRPPAG2-I and SPDRPPAG2-I with G=MH-1MT are Gaussian back substitute form methods and perform a little slower than Gaussian type methods. We also plot a figure to graphically illustrate the performance of four SRPPAs. Figure 1 shows the results from the test with n=1000 and n=6000, depicting error versus CPU time. Quality-wise, SPDRPPAG1-I was on par with SDPRPPAG1-I.Comparing SRPPAs applied to BP problem withn=1000 (a) and n=6000 (b). The horizontal axis gives the CPU time; the vertical axis gives the error between the solution and the original.
(a)
(b)In the second experiment, we compare the performance of SPDRPPAG1-I with TFOCS (source code can be found athttp://cvxr.com/tfocs/) [32], ADMM (source code can be found at http://www.stanford.edu/~boyd/papers/admm/) [33], and PPA-CM. To make the comparison independent of the stopping criterion for each algorithm, we first run TFOCS to get its solution xTFOCS and set a benchmark error
(68)benchmarkerr.=∥xTFOCS-xoriginal∥2
and then run other algorithms until they obtain smaller errors than this benchmark. TFOCS was stopped upon
(69)∥xk+1-xk∥max{1,∥xk+1∥}≤Tol.
Since we found that Tol=10-12 is small enough to guarantee very high accuracy, we set Tol=10-12 in TFOCS. The parameters of TFOCS and ADMM were taken with their defaults. To guarantee the convergence, fixed stepsizes r,s were set to s=100, r=1.01*∥AAT∥/s for PPA-CM. In SPDRPPAG1-I, we also choose the same convergence condition (#) and initial step size s=10, r=1 as the previous experiment. We varied the size of A from n=512 (m=n/2) to n=8500. The results of this experiment are summarized in Table 2. There, we report the run time in seconds, the number of iterations, and the error of the recovery solution. In Table 2, “—” means “out of memory.”Table 2
Performance of different iterative methods.
A
TFOCS
ADMM
PPA-CM
SPDRPPAG1-I
n
It./CPU
Err.
It./CPU
Err.
It./CPU
Err.
It./CPU
Err.
512
1681/4.16
7.5
e
-
10
492/0.67
5.0
e
-
10
1027/1.41
6.2
e
-
10
391/0.10
4.0
e
-
10
1024
406/1.68
3.7
e
-
10
274/0.53
3.0
e
-
10
848/0.76
2.3
e
-
10
214/0.28
3.0
e
-
10
2048
507/4.01
4.4
e
-
9
607/3.93
2.9
e
-
9
804/3.65
4.0
e
-
9
239/1.49
3.9
e
-
9
4096
933/19.50
2.6
e
-
9
1070/27.47
2.1
e
-
9
845/15.09
2.3
e
-
9
422/9.84
2.3
e
-
9
4200
461/9.91
2.1
e
-
9
916/25.48
1.8
e
-
9
868/16.27
2.1
e
-
9
391/9.44
1.5
e
-
9
4300
451/9.74
2.2
e
-
9
464/20.37
1.6
e
-
9
884/16.91
2.1
e
-
9
429/10.54
2.1
e
-
9
4600
505/12.37
1.36
e
-
9
2155/56.45
1.2
e
-
9
863/18.95
1.2
e
-
9
425/11.89
1.1
e
-
9
4700
801/20.15
1.3
e
-
11
1517/45.37
8.2
e
-
12
1102/34.27
1.2
e
-
11
425/14.20
1.1
e
-
11
5500
407/13.64
2.4
e
-
6
—/—
—
546/20.13
1.9
e
-
6
308/12.40
2.0
e
-
6
6500
1257/56.11
1.5
e
-
5
—/—
—
508/26.48
1.4
e
-
5
308/17.22
1.3
e
-
5
7500
801/81.42
1.1
e
-
12
—/—
—
1313/813.52
1.1
e
-
12
522/69.52
1.1
e
-
12
8500
842/107.18
2.3
e
-
7
—/—
—
724/100.36
2.2
e
-
7
542/83.07
2.2
e
-
7We observe from Table2 that four algorithms reach high accuracy around 10-9. SPDRPPAG1-I is about two times faster than the first-order method implemented in the TFOCS package, and moreover, it usually outperforms TFOCS in terms of iterations. For medium size problems, SPDRPPAG1-I is clearly faster than ADMM. Even for small size problems, SPDRPPAG1-I shows its superior performance. The main reason lies in that ADMM computed (I+AAT)-1 to solve its subproblem exactly which would take expensive computational cost. Not surprisingly, the general SPDRPPAG1-I performs better than the primary PPA-CM. Here, the total iterations of SPDRPPAG1-I are less than 50% of PPA-CM. As we have mentioned, “optimal” update stepsize αk and more flexible condition for convergence may provide SPDRPPAG1-I improved performance. SPDRPPAG1-I is faster than PPA-CM in terms of CPU times. However, the superiority of CPU time is not so significant as iteration number. For the cases n=4300, it is just about 62% of PPA-CM. This is not particularly surprising; compared to PPA-CM, SPDRPPAG1-I has to take extra computation for convergence condition and “optimal” αk in each iteration.


### 6.2. Nearest Correlation Matrix Problem

The nearest correlation matrix problem is solving the problem(70)min{12∥X-C∥F2∣diag(X)=e,X∈S+n},
where e∈ℜn is the vector whose entries are all 1s, S+n denotes the cone of positive definite symmetric matrices, diag(X) is the vector of diagonal elements of X, and ∥·∥F denotes the matrix Fröbenius norm ∥X∥F=trace(XTX)1/2.Here, we apply PPA-CM, LPPA, and SDPRPPA1-II for solving (70). The standard Matlab Mex interface mexsvd is used to conduct the eigenvalue decomposition. We constructed test data sets and stopping criterion like those of [24]. As mentioned in the prequel, we expect our SRPPA to produce robust performance. To assess the effectiveness of the adaptive strategies proposed in Section 6, we now move on to the description of experiments that focus on the consequences of the initial stepsizes. For investigating, we used dimensions n∈{500,1000} and varied s from 0.05 to 100, and initial points were set to 0 in all cases. Note that A=I; we fixed r=1.01/s for PPA-CM, r=0.65/s for LPPA and chose r=0.65/s as initial start for SDPRPPAG1-II. Since the experiments with other values of n give qualitatively similar results, we therefore do not plot those results to avoid clutter in the figures. The respective numerical results are plotted in Figure 2.CPU times as a function of the initial stepsizes for PPA-CM, LPPA, and SDPRPPAG1-II. The plot on the left is for n=500, while the plot on the right is for n=1000.
(a)
(b)It is clear that, for PPA-CM and LPPA, the convergence performance was a result of the stepsize selection. They are both fairly sensitive to initial stepsizes (or r). The results confirm that, with inappropriate stepsizes, both PPA-CM and LPPA become significantly slow. SDPRPPAG1-II yields significantly robust performance with adaptive strategy. And it is independent of the initial stepsizes and illustrates its superior performance. Furthermore, SDPRPPAG1-II yields competitive results even when PPA-CM and LPPA chose the “good” initial stepsize. This underlies the importance of adaptive strategy in producing good performance. Of course, care should be taken. For instance, the cost of computing optimal stepsize αk here is negligible, compared to the computation of SVD; when they are more costly, general LPPA will be expected to perform  slower than PPA-CM.


## 6.1.  Basis Pursuit Problem

Basis pursuit (BP) finds signal representations in overcomplete dictionaries by equality-constrainedl1 minimization problem. Formally, one solves the problem
(66)min{∥x∥1Ax=b,x∈ℜn}.
And here, ∥·∥1 denotes the l1 norm defined as ∥x∥1:=∑i=1n‍|xi|. BP is a fundamental problem in image processing and modern statistical signal processing, particularly the theory of compressed sensing; see, for example, [1–4] for intensive study. We now discuss our approach to BP problem of over-complete representations. Our experiments in this subsection use synthetic data which were mainly designed to illustrate the nice performance of our RPPA. The synthetic problem that we test here is similar to the one employed in [32]. We generate the data as follows: matrix A is a random m×n matrix, with Gaussian i.i.d. entries of zero mean and variance 1, with m=n/2. xoriginal∈Rn is the original sparse signal, constructed with m/5 nonzero values, randomly from standard normal distribution. We use xoriginal to generate the measurements as b=Axoriginal. It is desirable to use test problems that have a precisely known solution. In fact, when xoriginal is very sparse, it is the solution to the minimization problem (66). Hence, in our synthetic problem, xoriginal is exactly the solution.In our first experiment, we compared general RPPA using two differentG’s mentioned in Section 4.1. For BP problem, we use condition (#) and Algorithm 3. Since A constructed here is a general random matrix, and when A is large scale, ∥ATA∥ might be obtained costly. A simple stopping criterion
(67)err.=∥xk-xoriginal∥≤Tol
was used in this experiment, and the stopping tolerance Tol was set to 10-10. In all the tests, initial stepsizes were set as s=10, r=1, the primal variable x0 was initialized as zeros(n,1), and the dual λ0 was ones(m,1) in Matlab. Table 1 summarizes the performance of general SRPPA. Here, SPDRPPAGi-I(SPDRPPAGi-II) denotes self-adaptive primal-dual RPPA with Algorithm 3 (Algorithm 4), G=H, if i=1, and G=MH-1MT, if i=2. DPRPPA (DPRPPA) denotes dual-primal RPPA version.Table 1
Comparison of behaviours of four SRPPAs.
A
SPDRPPAG1-I
SDPRPPAG1-I
SPDRPPAG2-I
SDPRPPAG2-I
n
It.
CPU
It.
CPU
It.
CPU
It.
CPU
500
360
0.13
241
0.13
513
0.19
875
0.26
1500
379
1.32
259
0.89
387
1.72
594
2.04
2500
501
4.39
343
3.01
429
4.93
648
5.72
3500
399
6.65
325
5.57
663
14.44
1089
18.02
4500
389
12.13
457
14.47
493
20.61
991
31.26
5500
515
24.21
462
21.97
497
30.88
813
38.43
6500
449
25.01
423
23.41
456
33.42
932
51.62
7500
520
45.72
473
41.79
724
63.06
886
77.91Basically, SRPPAs converge very quickly and achieved tight error10-10 in a few hundred iterations. For this experiment, one can see that SDPRPPAG1-I is fastest in all cases. Both SDPRPPAG1-I and SPDRPPAG1-I are Gaussian type methods, with G=H, and they exhibit very similar performance. SDPRPPAG2-I and SPDRPPAG2-I with G=MH-1MT are Gaussian back substitute form methods and perform a little slower than Gaussian type methods. We also plot a figure to graphically illustrate the performance of four SRPPAs. Figure 1 shows the results from the test with n=1000 and n=6000, depicting error versus CPU time. Quality-wise, SPDRPPAG1-I was on par with SDPRPPAG1-I.Comparing SRPPAs applied to BP problem withn=1000 (a) and n=6000 (b). The horizontal axis gives the CPU time; the vertical axis gives the error between the solution and the original.
(a)
(b)In the second experiment, we compare the performance of SPDRPPAG1-I with TFOCS (source code can be found athttp://cvxr.com/tfocs/) [32], ADMM (source code can be found at http://www.stanford.edu/~boyd/papers/admm/) [33], and PPA-CM. To make the comparison independent of the stopping criterion for each algorithm, we first run TFOCS to get its solution xTFOCS and set a benchmark error
(68)benchmarkerr.=∥xTFOCS-xoriginal∥2
and then run other algorithms until they obtain smaller errors than this benchmark. TFOCS was stopped upon
(69)∥xk+1-xk∥max{1,∥xk+1∥}≤Tol.
Since we found that Tol=10-12 is small enough to guarantee very high accuracy, we set Tol=10-12 in TFOCS. The parameters of TFOCS and ADMM were taken with their defaults. To guarantee the convergence, fixed stepsizes r,s were set to s=100, r=1.01*∥AAT∥/s for PPA-CM. In SPDRPPAG1-I, we also choose the same convergence condition (#) and initial step size s=10, r=1 as the previous experiment. We varied the size of A from n=512 (m=n/2) to n=8500. The results of this experiment are summarized in Table 2. There, we report the run time in seconds, the number of iterations, and the error of the recovery solution. In Table 2, “—” means “out of memory.”Table 2
Performance of different iterative methods.
A
TFOCS
ADMM
PPA-CM
SPDRPPAG1-I
n
It./CPU
Err.
It./CPU
Err.
It./CPU
Err.
It./CPU
Err.
512
1681/4.16
7.5
e
-
10
492/0.67
5.0
e
-
10
1027/1.41
6.2
e
-
10
391/0.10
4.0
e
-
10
1024
406/1.68
3.7
e
-
10
274/0.53
3.0
e
-
10
848/0.76
2.3
e
-
10
214/0.28
3.0
e
-
10
2048
507/4.01
4.4
e
-
9
607/3.93
2.9
e
-
9
804/3.65
4.0
e
-
9
239/1.49
3.9
e
-
9
4096
933/19.50
2.6
e
-
9
1070/27.47
2.1
e
-
9
845/15.09
2.3
e
-
9
422/9.84
2.3
e
-
9
4200
461/9.91
2.1
e
-
9
916/25.48
1.8
e
-
9
868/16.27
2.1
e
-
9
391/9.44
1.5
e
-
9
4300
451/9.74
2.2
e
-
9
464/20.37
1.6
e
-
9
884/16.91
2.1
e
-
9
429/10.54
2.1
e
-
9
4600
505/12.37
1.36
e
-
9
2155/56.45
1.2
e
-
9
863/18.95
1.2
e
-
9
425/11.89
1.1
e
-
9
4700
801/20.15
1.3
e
-
11
1517/45.37
8.2
e
-
12
1102/34.27
1.2
e
-
11
425/14.20
1.1
e
-
11
5500
407/13.64
2.4
e
-
6
—/—
—
546/20.13
1.9
e
-
6
308/12.40
2.0
e
-
6
6500
1257/56.11
1.5
e
-
5
—/—
—
508/26.48
1.4
e
-
5
308/17.22
1.3
e
-
5
7500
801/81.42
1.1
e
-
12
—/—
—
1313/813.52
1.1
e
-
12
522/69.52
1.1
e
-
12
8500
842/107.18
2.3
e
-
7
—/—
—
724/100.36
2.2
e
-
7
542/83.07
2.2
e
-
7We observe from Table2 that four algorithms reach high accuracy around 10-9. SPDRPPAG1-I is about two times faster than the first-order method implemented in the TFOCS package, and moreover, it usually outperforms TFOCS in terms of iterations. For medium size problems, SPDRPPAG1-I is clearly faster than ADMM. Even for small size problems, SPDRPPAG1-I shows its superior performance. The main reason lies in that ADMM computed (I+AAT)-1 to solve its subproblem exactly which would take expensive computational cost. Not surprisingly, the general SPDRPPAG1-I performs better than the primary PPA-CM. Here, the total iterations of SPDRPPAG1-I are less than 50% of PPA-CM. As we have mentioned, “optimal” update stepsize αk and more flexible condition for convergence may provide SPDRPPAG1-I improved performance. SPDRPPAG1-I is faster than PPA-CM in terms of CPU times. However, the superiority of CPU time is not so significant as iteration number. For the cases n=4300, it is just about 62% of PPA-CM. This is not particularly surprising; compared to PPA-CM, SPDRPPAG1-I has to take extra computation for convergence condition and “optimal” αk in each iteration.


## 6.2. Nearest Correlation Matrix Problem

The nearest correlation matrix problem is solving the problem(70)min{12∥X-C∥F2∣diag(X)=e,X∈S+n},
where e∈ℜn is the vector whose entries are all 1s, S+n denotes the cone of positive definite symmetric matrices, diag(X) is the vector of diagonal elements of X, and ∥·∥F denotes the matrix Fröbenius norm ∥X∥F=trace(XTX)1/2.Here, we apply PPA-CM, LPPA, and SDPRPPA1-II for solving (70). The standard Matlab Mex interface mexsvd is used to conduct the eigenvalue decomposition. We constructed test data sets and stopping criterion like those of [24]. As mentioned in the prequel, we expect our SRPPA to produce robust performance. To assess the effectiveness of the adaptive strategies proposed in Section 6, we now move on to the description of experiments that focus on the consequences of the initial stepsizes. For investigating, we used dimensions n∈{500,1000} and varied s from 0.05 to 100, and initial points were set to 0 in all cases. Note that A=I; we fixed r=1.01/s for PPA-CM, r=0.65/s for LPPA and chose r=0.65/s as initial start for SDPRPPAG1-II. Since the experiments with other values of n give qualitatively similar results, we therefore do not plot those results to avoid clutter in the figures. The respective numerical results are plotted in Figure 2.CPU times as a function of the initial stepsizes for PPA-CM, LPPA, and SDPRPPAG1-II. The plot on the left is for n=500, while the plot on the right is for n=1000.
(a)
(b)It is clear that, for PPA-CM and LPPA, the convergence performance was a result of the stepsize selection. They are both fairly sensitive to initial stepsizes (or r). The results confirm that, with inappropriate stepsizes, both PPA-CM and LPPA become significantly slow. SDPRPPAG1-II yields significantly robust performance with adaptive strategy. And it is independent of the initial stepsizes and illustrates its superior performance. Furthermore, SDPRPPAG1-II yields competitive results even when PPA-CM and LPPA chose the “good” initial stepsize. This underlies the importance of adaptive strategy in producing good performance. Of course, care should be taken. For instance, the cost of computing optimal stepsize αk here is negligible, compared to the computation of SVD; when they are more costly, general LPPA will be expected to perform  slower than PPA-CM.


## 7. Conclusions

In this paper, we proposed an efficient general self-adaptive relaxed-PPA method for linearly constrained convex programming and provided theoretical convergence analysis for this method. The stepsizes choice condition is flexible and simple. Self-adaptive strategies are provided to make our method more efficient and robust. Experiments of the method in comparison to other state-of-art methods are provided to confirm the efficiency of the proposed method. Our numerical results suggest that SRPPA is effective and simple to implement. There are a few directions for further research, but we list here only two. The first is the question of whether we may modify the algorithm to work with more general constrained convex problems. Second, we aim to provide various relaxations of the subproblem for the practical purpose.



---

*Source: 492305-2013-09-19.xml*