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

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.

1. Introduction

In this paper, we consider the following generic convex programming: min{𝜃(𝑥)∣𝐴𝑥=𝑏(or𝐴𝑥≥𝑏),𝑥∈𝒳}, (1) where 𝜃(𝑥)∶ℜ^𝑛→ℜ is a convex (not necessary smooth) function, 𝐴∈ℜ^𝑚×𝑛, 𝑏∈ℜ^𝑚, and 𝒳⊂ℜ^𝑛 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: min_𝑥∈𝒳𝜃(𝑥)+𝜚2‖𝑥−𝑎‖², (2) where 𝑎 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, 𝑙₁-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 ‖𝐴𝑥−𝑏‖², its subproblems require an inverse operator of the form (𝐼+(1/𝑠)𝐴^𝑇𝐴)⁻¹ which is hard to implement in some cases. Particularly, 𝐴^𝑇𝐴 is general or large scale, so the computation of inverse operator may fail. Hence, ALM is not sufficiently competitive when the objective function 𝜃(𝑥) 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 Section 2, 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.

𝜕𝜃(𝑥) denotes the subdifferential set of the convex function 𝜃(𝑥): 𝜕𝜃(𝑥)∶=𝑑∈ℜ^𝑛∣𝜃𝑦−𝜃(𝑥)≥𝑑^𝑇𝑦−𝑥,∀𝑦∈ℜ^𝑛, (3) and 𝑑∈𝜕𝜃(𝑥) is called a subgradient of 𝜃(𝑥), see [26]. Let 𝑓(𝑥)∈𝜕𝜃(𝑥) and 𝑓(𝑦)∈𝜕𝜃(𝑦), by the convexity of the function 𝜃, we have 𝑥−𝑦^𝑇𝑓(𝑥)−𝑓𝑦≥0,∀𝑥,𝑦∈ℜ^𝑛, (4) which indicates that the mapping 𝑓 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 𝐴𝑥=𝑏(or𝐴𝑥≥𝑏), the Lagrangian function of (1) is 𝐿(𝑥,𝜆)=𝜃(𝑥)−𝜆^𝑇(𝐴𝑥−𝑏); (5) here, ℜΛ=^𝑚ℜ,fortheequalityconstraints𝐴𝑥=𝑏,_𝑚+,fortheinequalityconstraints𝐴𝑥≥𝑏, (6) and 𝐿(𝑥,𝜆) is defined on 𝒳×Λ. Then, by the optimality condition, we can easily see that (1) amounts to finding a pair of (𝑥^∗,𝜆^∗) which satisfies 𝑥^∗∈𝒳,𝑥−𝑥^∗𝑇𝑓𝑥^∗−𝐴^𝑇𝜆^∗𝜆≥0,∀𝑥∈𝒳,_∗∈Λ,𝜆−𝜆^∗𝑇𝐴𝑥^∗−𝑏≥0,∀𝜆∈Λ, (7) where 𝑓(𝑥^∗)∈𝜕𝜃(𝑥^∗). Denoting 𝑥𝜆(𝑢)=𝑓(𝑥)𝑢=,𝐹−𝐴^𝑇𝜆𝐴𝑥−𝑏,Ω=𝒳×Λ, (8) the system (7) can be characterized by the following variational inequality denoted by VI(Ω,𝐹): 𝑢^∗∈Ω,𝑢−𝑢^∗𝑇𝐹𝑢^∗≥0,∀𝑢∈Ω. (9) Recalling the monotonicity of 𝑓, it is easy to get that VI(Ω,𝐹) (9) is monotone. Since the solution set of (1) is assumed to be nonempty, the solution set of VI(Ω,𝐹), 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 for VI(Ω,𝐹) (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 iterate 𝑢^𝑘=(𝑥^𝑘,𝜆^𝑘), find 𝑢^𝑘=(𝑥^𝑘,𝜆^𝑘)∈Ω such that 𝑢^𝑘∈Ω,𝑢−𝑢^𝑘𝑇𝐹𝑢^𝑘+𝑟𝑢^𝑘−𝑢^𝑘≥0,∀𝑢∈Ω. (10) Then, the update step is taken as follows: 𝑢^𝑘+1=𝑢^𝑘−𝛾𝑢^𝑘−𝑢^𝑘,𝛾∈(0,2). (11) PPA has a nice convergence property 𝑢^𝑘+1−𝑢^∗2≤𝑢^𝑘−𝑢^∗2−𝛾2−𝛾𝑢^𝑘−𝑢^𝑘2. (12) Although classical PPA is conceptually appealing, subproblem (10) presents certain computational challenges. More specifically, primal variable 𝑥 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 𝐹(𝑢^𝑘) with 𝐹(𝑢^𝑘) in the proximal subproblem (10), which amounts to the following subproblem: 𝑢^𝑘∈Ω,𝑢−𝑢^𝑘𝑇𝐹𝑢^𝑘+𝑟𝑢^𝑘−𝑢^𝑘≥0,∀𝑢∈Ω. (13) The solution of the relaxed problem (13) is given by 𝑢^𝑘=𝑃_Ω[𝑢^𝑘−(1/𝑟)𝐹(𝑢^𝑘)]. It is clear that methods with such relaxation are explicit type methods. However, 𝑢^𝑘 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 predictor 𝑢^𝑘 to obtain corrector as the new iterate. (i)The extragradient method (EGM) updates the new iterate (corrector) by 𝑢^𝑘+1=𝑃_Ω𝑢^𝑘−1𝑟𝐹𝑢^𝑘. (14)(ii)The projection and contraction methods (PCM) [28–30] perform update as follows: 𝑢^𝑘+1=𝑢^𝑘−𝛾𝛼^∗_𝑘𝑑^𝑘or𝑢^𝑘+1=𝑃_Ω𝑢^𝑘−𝛾𝛼^∗_𝑘𝑟𝐹𝑢^𝑘, (15) where 𝑑^𝑘=𝑢^𝑘−𝑢^𝑘−1𝑟𝐹𝑢^𝑘−𝐹𝑢^𝑘,𝛼^∗_𝑘=𝑢^𝑘−𝑢^𝑘𝑇𝑑^𝑘𝑑^𝑘2. (16)The sequence {𝑢^𝑘} generated by the above mentioned EGM or PCM satisfies 𝑢^𝑘+1−𝑢^∗2≤𝑢^𝑘−𝑢^∗2−𝑐𝑢^𝑘−𝑢^𝑘2,𝑐>0, (17) which is similar to PPA. Both EGM and PCM use the simplest relaxation to obtain 𝑢^𝑘 in 𝑘th iteration, hence are computationally practical. These methods have appealing practical aspects; however, such simplest relaxation does not exploit the inner structure of VI(Ω,𝐹) (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 𝑢^𝑘 by solving the following subproblem: find (𝑥^𝑘,𝜆^𝑘)∈Ω such that 𝑢^𝑘∈Ω,𝑢−𝑢^𝑘𝑇𝐹𝑢^𝑘+𝑆𝑢^𝑘−𝑢^𝑘≥0,∀𝑢∈Ω, (18) where 𝑆=𝑟𝐼_𝑛−𝐴^𝑇−𝐴𝑠𝐼_𝑚. (19) And perform the update 𝑢^𝑘+1=𝑢^𝑘−𝛾𝑢^𝑘−𝑢^𝑘,𝛾∈(0,2). (20) The framework of LPPA is as follows: 𝑢^𝑘∈Ω,𝑢−𝑢^𝑘𝑇𝐹𝑢^𝑘+𝑀𝑢^𝑘−𝑢^𝑘≥0,∀𝑢∈Ω, (21) where 𝑀=𝑟𝐼_𝑛𝐴^𝑇0𝑠𝐼_𝑚. (22) And the new iterate is defined by 𝑢^𝑘+1=𝑢^𝑘−𝛾𝛼𝑀𝑢^𝑘−𝑢^𝑘,𝛾∈(0,2). (23) Both procedures are simple and can solve subproblem efficiently; but their nice convergence results require a quite restrictive condition, that is; 𝑟𝑠>‖𝐴^𝑇𝐴‖ in PPA-CM and 𝑟𝑠>(1/2)‖𝐴^𝑇𝐴‖ in LPPA, respectively. The stepsizes 𝑟,𝑠 are directly determined by such condition; hence, it is important to estimate ‖𝐴^𝑇𝐴‖. 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 𝑟,𝑠. 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 𝑥 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 𝐺 be a positive symmetry definite matrix (we will specify it later), 𝑀=𝑟𝐼_𝑛𝐴^𝑇0𝑠𝐼_𝑚,𝐻=𝑟𝐼_𝑛00𝑠𝐼_𝑚𝑑𝑢_𝑘,,𝑢^𝑘=𝐺⁻¹𝑀𝑢^𝑘−𝑢^𝑘,𝜑𝑢^𝑘,𝑢^𝑘=𝑢^𝑘−𝑢^𝑘𝑇𝑀𝑢^𝑘−𝑢^𝑘. (24)(25)(26) 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 𝑥-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.

Step  1. Initialization. Let 𝛾∈(0,2) and pick (𝑥0,𝜆0)∈ℜ𝑛×Λ, set 𝑘=0. Step  2. Predictor. Let      𝑥𝑘𝑟=Argmin𝜃(𝑥)+2𝑥−𝑥𝑘+1𝑟𝐴𝑇𝜆𝑘2∣𝑥∈𝒳,     (*) and          𝜆𝑘=𝑃Λ𝜆𝑘−1𝑠(𝐴𝑥𝑘−𝑏).                 (**) Step  3. If the stepsizes 𝑟 and 𝑠 are chosen satisfying          𝜑(𝑢𝑘,𝑢𝑘1)≥4𝑑(𝑢𝑘,𝑢𝑘)2𝐺,                  (#) then go to Step 4. Otherwise, increase 𝑟,𝑠 and go back to Step 2. Step  4. Corrector and updating          𝑢𝑘+1=𝑢𝑘−𝛼𝑘𝐺−1𝑀(𝑢𝑘−𝑢𝑘).                ($) Step  5. Adjustment   𝑟(𝑟,𝑠)=2,𝑠2,if𝜑(𝑢𝑘,𝑢𝑘)≥𝜅∗𝑑(𝑢𝑘,𝑢𝑘)2𝐺;(𝑟,𝑠),otherwise.Here𝜅>4. Set 𝑘∶=𝑘+1.

Algorithm 1 

General primal-dual relaxed-PPA method.

Step  1. Initialization. Let 𝛾∈(0,2) and pick (𝑥0,𝜆0)∈ℜ𝑛×Λ, set 𝑘=0. Step  2. Predictor. Let           𝜆𝑘=𝑃Λ𝜆𝑘−1𝑠(𝐴𝑥𝑘−𝑏),                 (†) and      𝑥𝑘𝑟=Argmin𝜃(𝑥)+2𝑥−𝑥𝑘+1𝑟𝐴𝑇𝜆𝑘2∣𝑥∈𝒳.        (††) Step  3. If the stepsize 𝑟 and 𝑠 are chosen satisfying           𝜑𝑢𝑘,𝑢𝑘≥14𝑑𝑢𝑘,𝑢𝑘2𝐺,                    (‡) then go to Step 4. Otherwise, increase 𝑟,𝑠 and go back to Step 2. Step  4. Corrector and updating          𝑢𝑘+1=𝑢𝑘−𝛼𝑘𝐺−1𝑀(𝑢𝑘−𝑢𝑘).                    (∧) Step  5. Adjustment   𝑟(𝑟,𝑠)=2,𝑠2,if𝜑𝑢𝑘,𝑢𝑘≥𝜅∗𝑑𝑢𝑘,𝑢𝑘2𝐺;(𝑟,𝑠),otherwise.Here𝜅>4. Set 𝑘∶=𝑘+1.

Algorithm 2 

General dual-primal relaxed-PPA method.

Remark 3. We first make some insight into the correction step in Algorithm 1. The obtained 𝑢^𝑘 plays no direct role as the new iterate. Instead, we need some kind of “corrector” defined in ($). Although matrix 𝐺 in ($) is just a required positive symmetry definite, our goal here is to fully integrate the information of 𝑢^𝑘 and 𝑢^𝑘 to construct effective, simple search direction 𝐺⁻¹𝑀(𝑢^𝑘−𝑢^𝑘) for the corrector. Based on this consideration, we elaborately provide two simple choices of 𝐺.

Case 1. It is natural to set 𝐺=𝐻 to induce a simple update form. Then, it is easy to get that 𝑥^𝑘+1𝜆^𝑘+1=𝑥^𝑘𝜆^𝑘−𝛼_𝑘𝑥^𝑘−𝑥^𝑘+1𝑟𝐴^𝑇𝜆^𝑘−𝜆^𝑘𝜆^𝑘−𝜆^𝑘. (27)

Case 2. Let 𝐺=𝑀𝐻⁻¹𝑀^𝑇. 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 𝐺 in ($), we get 𝑢^𝑘+1=𝑢^𝑘−𝛼_𝑘𝑀^−𝑇𝐻𝑀⁻¹𝑀𝑢^𝑘−𝑢^𝑘=𝑢^𝑘−𝛼_𝑘𝑀^−𝑇𝐻𝑢^𝑘−𝑢^𝑘. (28)

Recalling 𝑀 is a lower triangular matrix, by the fact that its inverse is also a lower triangular, we have 𝐺⁻¹𝑀=𝑀^−𝑇𝐼𝐻=^𝑛0−𝐴𝑠𝐼_𝑚. (29)

Plugging the previous relationship to (31), we have 𝑥^𝑘+1𝜆^𝑘+1=𝑥^𝑘𝜆^𝑘−𝛼_𝑘𝑥^𝑘−𝑥^𝑘𝜆^𝑘−𝜆^𝑘−1𝑠𝐴𝑥^𝑘−𝑥^𝑘. (30)

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 (𝑥^𝑘,𝜆^𝑘)∈Ω such that 𝑥−𝑥^𝑘𝜆−𝜆^𝑘𝑇𝑓𝑥^𝑘−𝐴^𝑇𝜆^𝑘𝐴𝑥^𝑘+−𝑏𝑟𝐼_𝑛𝐴^𝑇0𝑠𝐼_𝑚𝑥^𝑘−𝑥^𝑘𝜆^𝑘−𝜆^𝑘𝑥𝜆≥0,∀∈Ω, (31) and its compact form 𝑢^𝑘∈Ω,𝑢−𝑢^𝑘𝑇𝐹𝑢^𝑘+𝑀𝑢^𝑘−𝑢^𝑘≥0,∀𝑢∈Ω. (32) We observe that subproblem (32) is similar to (10) in PPA, except for the construction of asymmetry matrix 𝑀. 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 𝑟,𝑠. It is clear that condition (#) does not need prior knowledge of matrix 𝐴. Furthermore, it only involves matrix-vector multiplication, and so, it is easy to verify, and it is amenable to large-scale 𝐴. If 𝑟,𝑠 fail to meet this convergence condition (#), one should appropriately increase 𝑟,𝑠. 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: 𝑟,𝑠 satisfy 1𝑟𝑠≥2𝐴^𝑇𝐴. (33) It is stronger than condition (#). The following lemma is devoted to the proof of this result.

Lemma 5. Let {𝑢^𝑘} be the sequence generated by Algorithm 1, 𝐻, 𝑑(𝑢^𝑘,𝑢^𝑘), and 𝜑(𝑢^𝑘,𝑢^𝑘) defined in (24), (25), and (26), respectively. Suppose that condition (33) is satisfied. Then, condition (#) holds immediately.

Proof. Note that 𝜑𝑢^𝑘,𝑢^𝑘=𝑟𝑥^𝑘−𝑥^𝑘2+𝑠𝜆^𝑘−𝜆^𝑘2+𝜆^𝑘−𝜆^𝑘𝑇𝐴𝑥^𝑘−𝑥^𝑘. (34) Recall that matrix 𝐺 described in Algorithm 1 can be designed in two different cases.
Case  1. If 𝐺=𝐻, we immediately have 𝑑𝑢^𝑘,𝑢^𝑘2_𝐺=𝑢^𝑘−𝑢^𝑘2_𝐻+1𝑟𝐴^𝑇(𝜆^𝑘−𝜆^𝑘)²+2𝜆^𝑘−𝜆^𝑘𝑇𝐴𝑥^𝑘−𝑥^𝑘. (35)
According to Cauchy-Schwarz inequality, we get 𝜑𝑢^𝑘,𝑢^𝑘−14𝑑𝑢^𝑘,𝑢^𝑘2_𝐺=34𝑢^𝑘−𝑢^𝑘2_𝐻−1𝐴4𝑟^𝑇𝜆^𝑘−𝜆^𝑘2+12𝜆^𝑘−𝜆^𝑘𝑇𝐴𝑥^𝑘−𝑥^𝑘≥34𝑢^𝑘−𝑢^𝑘2_𝐻−1𝐴4𝑟^𝑇(𝜆^𝑘−𝜆^𝑘)²−142𝑟𝑥^𝑘−𝑥^𝑘2+1𝐴2𝑟^𝑇𝜆^𝑘−𝜆^𝑘2=𝑟4𝑥^𝑘−𝑥^𝑘2+3𝑠4𝜆^𝑘−𝜆^𝑘2−3𝐴8𝑟^𝑇(𝜆^𝑘−𝜆^𝑘)²≥𝑟4𝑥^𝑘−𝑥^𝑘2. (36)
The last inequality follows directly from condition (33).
Case  2. If 𝐺=𝑀𝐻⁻¹𝑀^𝑇, by the definition of 𝑑(𝑢^𝑘,𝑢^𝑘), we obtain 𝑑𝑢^𝑘,𝑢^𝑘2_𝐺=𝑢^𝑘−𝑢^𝑘𝑇𝑀^𝑇𝐺⁻¹𝑀𝑢^𝑘−𝑢^𝑘=𝑢^𝑘−𝑢^𝑘2_𝐻. (37)
Then, we get 𝜑𝑢^𝑘,𝑢^𝑘−14𝑑𝑢^𝑘,𝑢^𝑘2_𝐺=34𝑢^𝑘−𝑢^𝑘2_𝐻+𝜆^𝑘−𝜆^𝑘𝑇𝐴𝑥^𝑘−𝑥^𝑘≥34𝑢^𝑘−𝑢^𝑘2_𝐻−124𝑟3𝑥^𝑘−𝑥^𝑘2−123𝐴4𝑟^𝑇𝜆^𝑘−𝜆^𝑘2=𝑟𝑥12^𝑘−𝑥^𝑘2+3𝑠4𝜆^𝑘−𝜆^𝑘2−3𝐴8𝑟^𝑇(𝜆^𝑘−𝜆^𝑘)²≥𝑟𝑥12^𝑘−𝑥^𝑘2. (38)
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 if 𝑟𝑠≥(1/2)‖𝐴^𝑇𝐴‖.

Condition (33) is not only stronger than Condition (#), but it also requires that matrix 𝑀 is positive semidefinite, while condition (#) does not. Furthermore, condition (33) may require the explicit expression of 𝐴 or knowledge of ‖𝐴^𝑇𝐴‖. Despite these drawbacks, condition (33) is appealing to the problems in which ‖𝐴^𝑇𝐴‖ is known beforehand or easy to compute/obtain. For instance, 𝐴 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 𝜑(𝑢^𝑘,𝑢^𝑘) with favourable property. In fact, for general matrix 𝐺, condition (#) also can guarantee convergence.

Remark 6. The update stepsize 𝛼_𝑘 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 Algorithm 1 within a limited number to avoid divergence.

In Algorithm 1, we carry out the 𝑥-predictor before performing 𝜆-predictor. The roles of 𝑥 and 𝜆 are symmetric; hence, sweeping the order will not break the Gauss-Seidel structure. We switch 𝑥 and 𝜆 and obtain a variant of relaxed-PPA with the order of the 𝑥-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 𝑀=𝑟𝐼_𝑛0−𝐴𝑠𝐼_𝑚. (39)

4.2. Adaptive Enhancements

Both PPA-CM and LPPA employ fixed stepsizes 𝑟,𝑠. Experiments reveal that they will suffer slow convergence when stepsizes 𝑟,𝑠 are chosen inappropriately. A natural question is, how to choose the proper initial stepsizes 𝑟,𝑠. 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 𝑑^𝑘_𝑥𝑑^𝑘_𝜆=𝑥^𝑘−𝑥^𝑘+1𝑟𝐴^𝑇𝜆^𝑘−𝜆^𝑘𝜆^𝑘−𝜆^𝑘, (40) and then, (31) can be rewritten as 𝑥−𝑥^𝑘𝜆−𝜆^𝑘𝑇𝑓𝑥^𝑘−𝐴^𝑇𝜆^𝑘𝐴𝑥^𝑘−−𝑏𝑟𝐼_𝑛00𝑠𝐼_𝑚𝑑^𝑘_𝑥𝑑^𝑘_𝜆∀𝑥𝜆≥0,∈Ω. (41) Under 𝐻-norm, the quantity 𝑑^𝑘_𝑥 can be viewed as a residual for the dual feasibility condition, and 𝑑^𝑘_𝜆 can be viewed as a primal residual. These two residuals converge to zero as relaxed-PPA proceeds. Note that 𝑑^𝑘_𝑥²_𝑟=𝑟𝑥^𝑘−𝑥^𝑘2+2𝑥^𝑘−𝑥^𝑘𝑇𝐴^𝑇𝜆^𝑘−𝜆^𝑘+1𝑟𝐴^𝑇𝜆^𝑘−𝜆^𝑘2,𝑑^𝑘_𝜆²_𝑠=𝑠𝜆^𝑘−𝜆^𝑘2. (42) And this implies that small values of 𝑠 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 𝑟,𝑠 properly according to the adaptation shown in Algorithm 3.

If𝑑𝑘𝑥2𝑟≥𝜏1𝑑𝑘𝜆2𝑠  𝑟∶=𝑟;𝑠∶=𝑠∗2; elseif𝜏2𝑑𝑘𝑥2𝑟≤𝑑𝑘𝜆2𝑠  𝑟∶=𝑟∗2;𝑠∶=𝑠; else  𝑟∶=𝑟∗1.5;𝑠∶=𝑠∗1.5.

Algorithm 3 

Adaptation-I.

Here, 𝜏₁>1, 𝜏₂>1. This adaptation strategy increases 𝑠 when the dual residual ‖𝑑^𝑘_𝑥‖²_𝑟 appears large compared to the primal residual ‖𝑑^𝑘_𝜆‖²_𝑠 and increases 𝑟 when the dual residual ‖𝑑^𝑘_𝑥‖²_𝑟 seems too small relative to the primal residual ‖𝑑^𝑘_𝜆‖²_𝑠.

As mentioned, condition (33) is stronger than condition (#). If one chooses condition (33), our RPPA also converges. It must have predetermined stepsizes satisfying 𝑟𝑠=𝜇‖𝐴^𝑇𝐴‖ (here, 𝜇≥0.5). However, there is no priority knowledge of the choice of individual 𝑟 or 𝑠. Here, we can also adjust 𝑟,𝑠 automatically when choosing condition (33). Intuitively, we consider expansion of the entire residual under 𝐻-norm: 𝑑^𝑘_𝑥²_𝑟+𝑑^𝑘_𝜆²_𝑠=𝑟𝑥^𝑘−𝑥^𝑘2+𝑠𝜆^𝑘−𝜆^𝑘2+1𝑟𝐴^𝑇𝜆^𝑘−𝜆^𝑘2+2𝑥^𝑘−𝑥^𝑘𝑇𝐴^𝑇𝜆^𝑘−𝜆^𝑘; (43) there are three terms involving 𝑟 or 𝑠, and we intend to balance these terms. For fixed 𝜇, take 𝑠=(𝜇/𝑟)‖𝐴𝐴^𝑇‖; then 𝑠𝜆^𝑘−𝜆^𝑘2=𝜇𝑟𝐴𝐴^𝑇𝜆^𝑘−𝜆^𝑘2. (44) Applying (44) into (43), clearly we have 𝑑^𝑘_𝑥²_𝑟+𝑑^𝑘_𝜆²_𝑠=𝑟𝑥^𝑘−𝑥^𝑘2+1𝑟𝜇𝐴𝐴^𝑇𝜆^𝑘−𝜆^𝑘2+𝐴^𝑇𝜆^𝑘−𝜆^𝑘2+2𝑥^𝑘−𝑥^𝑘𝑇𝐴^𝑇𝜆^𝑘−𝜆^𝑘. (45) Now, we consider adjusting stepsize to balance 𝑟‖𝑥^𝑘−𝑥^𝑘‖² and (1/𝑟)(𝜇‖𝐴𝐴^𝑇‖‖𝜆^𝑘−𝜆^𝑘‖²+‖𝐴^𝑇(𝜆^𝑘−𝜆^𝑘)‖²) and obtain another adaptation strategy (see Algorithm 4).

if𝑟𝑥𝑘−𝑥𝑘2≥𝜏1𝑠𝜆𝑘−𝜆𝑘2+1𝑟𝐴𝑇(𝜆𝑘−𝜆𝑘)2  𝑟𝑟∶=2𝜇;𝑠=𝐴𝐴𝑇𝑟; elseif𝜏2𝑟𝑥𝑘−𝑥𝑘2≤𝑠𝜆𝑘−𝜆𝑘2+1𝑟𝐴𝑇(𝜆𝑘−𝜆𝑘)2  𝑠𝑠∶=2𝜇;𝑟=𝐴𝐴𝑇𝑠; else  𝑟∶=𝑟;𝑠∶=𝑠.

Algorithm 4 

Adaptation-II.

It is worth noting that too many adjustments of stepsizes by Algorithm 4 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.

Let 𝑢^∗=(𝑥^∗,𝜆^∗) be any solution point, setting 𝑢=𝑢^∗ in (32) yields 𝑢^𝑘−𝑢^∗𝑇𝑀𝑢^𝑘−𝑢^𝑘≥𝑢^𝑘−𝑢^∗𝑇𝐹𝑢^𝑘. (46) Since 𝑢^𝑘∈Ω, we have (𝑢^𝑘−𝑢^∗)^𝑇𝐹(𝑢^∗)≥0. Consequently, by using the monotonicity of 𝐹, the right hand side of (46) is nonnegative, and thus 𝑢^𝑘−𝑢^∗𝑇𝑀𝑢^𝑘−𝑢^𝑘≥0. (47) Now, we are writing the update as 𝑢(𝛼)=𝑢^𝑘−𝛼𝑑𝑢^𝑘,𝑢^𝑘 (48) and 𝜗(𝛼)∶=𝑢^𝑘−𝑢^∗2_𝐺−𝑢(𝛼)−𝑢^∗2_𝐺,𝑞(𝛼)=2𝛼𝑢^𝑘−𝑢^𝑘𝑇𝐺𝑑𝑢^𝑘,𝑢^𝑘−𝛼²𝑑𝑢^𝑘,𝑢^𝑘2_𝐺. (49) Here, 𝜗(𝛼) can be viewed as a progress function. The following lemma shows that 𝑞(𝛼) is a lower bound of 𝜗(𝛼).

Lemma 8. Let 𝜗(𝛼) and 𝑞(𝛼) be defined in (49); then one has 𝜗(𝛼)≥𝑞(𝛼). (50)

Proof. Let 𝑢^∗ be any solution, from the definition of 𝑢(𝛼), we have 𝜗(𝛼)=𝑢^𝑘−𝑢^∗2_𝐺−𝑢^𝑘−𝑢^∗−𝛼𝑑𝑢^𝑘,𝑢^𝑘2_𝐺=2𝛼𝑢^𝑘−𝑢^∗𝑇𝐺𝑑𝑢^𝑘,𝑢^𝑘−𝛼²𝑑𝑢^𝑘,𝑢^𝑘2_𝐺. (51) Applying (47) to the first term of (51) gives 𝑢^𝑘−𝑢^∗𝑇𝐺𝑑𝑢^𝑘,𝑢^𝑘=𝑢^𝑘−𝑢^𝑘𝑇𝐺𝑑𝑢^𝑘,𝑢^𝑘+𝑢^𝑘−𝑢^∗𝑇𝐺𝑑𝑢^𝑘,𝑢^𝑘≥𝑢^𝑘−𝑢^𝑘𝑇𝐺𝑑𝑢^𝑘,𝑢^𝑘. (52) Substituting (52) into (51), we immediately obtain the assertion.

We note that 𝑞(𝛼) is a quadratic function of 𝛼 and it is natural to maximize 𝑞(𝛼) to obtain an “optimal” 𝛼: 𝛼^∗_𝑘=𝑢^𝑘−𝑢^𝑘𝑇𝐺𝑑𝑢^𝑘,𝑢^𝑘𝑑𝑢^𝑘,𝑢^𝑘2_𝐺. (53) We now show that the “optimal” 𝛼^∗_𝑘 is bounded above from zero in the following Lemma.

Lemma 9. Let sequence {𝑢^𝑘} be produced by Algorithm 1, 𝛼^∗_𝑘 defined in (53); then, one has 𝛼^∗_𝑘≥14>0. (54)

Proof. Using the definition of 𝛼^∗_𝑘 in (53), we have, for all 𝑘, 𝛼^∗_𝑘=𝜑𝑢^𝑘,𝑢^𝑘𝑑𝑢^𝑘,𝑢^𝑘2_𝐺≥14. (55) The inequality follows from condition (#).

Setting 𝛼=𝛼_𝑘=𝛼^∗_𝑘𝛾 in (50) yields 𝑢^𝑘+1−𝑢^∗2_𝐺≤𝑢^𝑘−𝑢^∗2_𝐺−𝛾2−𝛾𝛼^∗_𝑘𝑢^𝑘−𝑢^𝑘𝑇𝐺𝑑𝑢^𝑘,𝑢^𝑘. (56) Combining Lemmas 8 and 9, we immediately obtain the following convergence theorem.

Theorem 10. Let sequence {𝑢^𝑘} be produced by Algorithm 1; then one gets 𝑢^𝑘+1−𝑢^∗2_𝐺≤𝑢^𝑘−𝑢^∗2_𝐺−𝛾2−𝛾𝑑𝑢16^𝑘,𝑢^𝑘2_𝐺. (57)

Theorem 11. Let sequence {𝑢^𝑘} be generated by Algorithm 1. Then, {𝑢^𝑘} converges to some 𝑢^∞ which is a solution of  𝑉𝐼(Ω,𝐹) (9).

Proof. First, for each 𝑢∈Ω, we have 𝑢−𝑢^𝑘𝑇𝐹𝑢^𝑘≥𝑢−𝑢^𝑘𝑇𝑀𝑢^𝑘−𝑢^𝑘. (58) It follows from (57) that {𝑢^𝑘} is a bounded sequence and lim_𝑘→∞𝑢^𝑘−𝑢^𝑘_𝐺=0. (59) Consequently, {𝑢^𝑘} is also bounded. Since lim_𝑘→∞‖𝑢^𝑘−𝑢^𝑘‖_𝐺=0, it follows from (58) that lim_𝑘→∞𝑢−𝑢^𝑘𝑇𝐹𝑢^𝑘≥0,∀𝑢∈Ω. (60) Because {𝑢^𝑘} is bounded, it has at least a cluster point. Let 𝑢^∞ be a cluster point of {𝑢^𝑘} and let the subsequence {𝑢^𝑘𝑗} converge to 𝑢^∞. It follows that lim_𝑗→∞𝑢−𝑢^𝑘𝑗𝑇𝐹𝑢^𝑘𝑗≥0,∀𝑢∈Ω, (61) and consequently, 𝑢−𝑢^∞𝑇𝐹𝑢^∞≥0.∀𝑢∈Ω. (62) This means that 𝑢^∞ is a solution of VI(Ω,𝐹). Note that the inequality (57) is true for all solution points of VI(Ω,𝐹), and hence, we have 𝑢^𝑘+1−𝑢^∞2_𝐺≤𝑢^𝑘−𝑢^∞2_𝐺,∀𝑘≥0. (63) Since 𝑢^𝑘𝑗→𝑢^∞(𝑗→∞) and 𝑢^𝑘−𝑢^𝑘→0(𝑘→∞), for any given 𝜀>0, there exists an integer 𝑙>0 such that 𝑢^𝑘𝑙−𝑢^∞_𝐺<𝜀2,𝑢^𝑘𝑙−𝑢^𝑘𝑙_𝐺<𝜀2. (64) Therefore, for any 𝑘≥𝑘_𝑙, it follows from (63) and (64) that 𝑢^𝑘−𝑢^∞_𝐺≤𝑢^𝑘𝑙−𝑢^∞_𝐺≤𝑢^𝑘𝑙−𝑢^𝑘𝑙_𝐺+𝑢^𝑘𝑙−𝑢^∞_𝐺≤𝜀. (65) This implies that the sequence {𝑢^𝑘} converges to 𝑢^∞, which is a solution of VI(Ω,𝐹) (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-constrained 𝑙₁ minimization problem. Formally, one solves the problem min‖𝑥‖₁𝐴𝑥=𝑏,𝑥∈ℜ^𝑛. (66) And here, ‖⋅‖₁ denotes the 𝑙₁ norm defined as ‖𝑥‖₁∶=∑^𝑛_𝑖=1|𝑥_𝑖|. 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 𝐴 is a random 𝑚×𝑛 matrix, with Gaussian i.i.d. entries of zero mean and variance 1, with 𝑚=𝑛/2. 𝑥_original∈𝑅^𝑛 is the original sparse signal, constructed with 𝑚/5 nonzero values, randomly from standard normal distribution. We use 𝑥_original to generate the measurements as 𝑏=𝐴𝑥_original. It is desirable to use test problems that have a precisely known solution. In fact, when 𝑥_original is very sparse, it is the solution to the minimization problem (66). Hence, in our synthetic problem, 𝑥_original is exactly the solution.

In our first experiment, we compared general RPPA using two different 𝐺’s mentioned in Section 4.1. For BP problem, we use condition (#) and Algorithm 3. Since 𝐴 constructed here is a general random matrix, and when 𝐴 is large scale, ‖𝐴^𝑇𝐴‖ might be obtained costly. A simple stopping criterion err.=𝑥^𝑘−𝑥_original≤Tol (67) was used in this experiment, and the stopping tolerance Tol was set to 10⁻¹⁰. In all the tests, initial stepsizes were set as 𝑠=10, 𝑟=1, the primal variable 𝑥⁰ was initialized as zeros(𝑛,1), and the dual 𝜆⁰ was ones(𝑚,1) in Matlab. Table 1 summarizes the performance of general SRPPA. Here, SPDRPPAG𝑖-I(SPDRPPAG𝑖-II) denotes self-adaptive primal-dual RPPA with Algorithm 3 (Algorithm 4), 𝐺=𝐻, if 𝑖=1, and 𝐺=𝑀𝐻⁻¹𝑀^𝑇, if 𝑖=2. DPRPPA (DPRPPA) denotes dual-primal RPPA version.

Table 1 

Comparison of behaviours of four SRPPAs.

Basically, SRPPAs converge very quickly and achieved tight error 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 𝐺=𝐻, and they exhibit very similar performance. SDPRPPAG2-I and SPDRPPAG2-I with 𝐺=𝑀𝐻⁻¹𝑀^𝑇 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 𝑛=1000 and 𝑛=6000, depicting error versus CPU time. Quality-wise, SPDRPPAG1-I was on par with SDPRPPAG1-I.

(a)

(b)

Figure 1 

Comparing SRPPAs applied to BP problem with 𝑛=1000 (a) and 𝑛=6000 (b). The horizontal axis gives the CPU time; the vertical axis gives the error between the solution and the original.

In the second experiment, we compare the performance of SPDRPPAG1-I with TFOCS (source code can be found at http://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 𝑥_TFOCS and set a benchmark error benchmarkerr.=𝑥_TFOCS−𝑥_original2 (68) and then run other algorithms until they obtain smaller errors than this benchmark. TFOCS was stopped upon 𝑥^𝑘+1−𝑥^𝑘max1,𝑥^𝑘+1≤Tol. (69) Since we found that Tol=10⁻¹² is small enough to guarantee very high accuracy, we set Tol=10⁻¹² in TFOCS. The parameters of TFOCS and ADMM were taken with their defaults. To guarantee the convergence, fixed stepsizes 𝑟,𝑠 were set to 𝑠=100, 𝑟=1.01∗‖𝐴𝐴^𝑇‖/𝑠 for PPA-CM. In SPDRPPAG1-I, we also choose the same convergence condition (#) and initial step size 𝑠=10, 𝑟=1 as the previous experiment. We varied the size of 𝐴 from 𝑛=512 (𝑚=𝑛/2) to 𝑛=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.

We observe from Table 2 that four algorithms reach high accuracy around 10⁻⁹. 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 (𝐼+𝐴𝐴^𝑇)⁻¹ 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 𝛼_𝑘 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 𝑛=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” 𝛼_𝑘 in each iteration.

6.2. Nearest Correlation Matrix Problem

The nearest correlation matrix problem is solving the problem 1min2‖𝑋−𝐶‖²_𝐹∣diag(𝑋)=𝑒,𝑋∈𝑆^𝑛₊, (70) where 𝑒∈ℜ^𝑛 is the vector whose entries are all 1𝑠, 𝑆^𝑛₊ denotes the cone of positive definite symmetric matrices, diag(𝑋) is the vector of diagonal elements of 𝑋, and ‖⋅‖_𝐹 denotes the matrix Fröbenius norm ‖𝑋‖_𝐹=trace(𝑋^𝑇𝑋)^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 𝑛∈{500,1000} and varied 𝑠 from 0.05 to 100, and initial points were set to 0 in all cases. Note that 𝐴=𝐼; we fixed 𝑟=1.01/𝑠 for PPA-CM, 𝑟=0.65/𝑠 for LPPA and chose 𝑟=0.65/𝑠 as initial start for SDPRPPAG1-II. Since the experiments with other values of 𝑛 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.

(a)

(b)

Figure 2 

CPU times as a function of the initial stepsize 𝑠 for PPA-CM, LPPA, and SDPRPPAG1-II. The plot on the left is for 𝑛=500, while the plot on the right is for 𝑛=1000.

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 stepsize 𝑠 (or 𝑟). 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 𝛼_𝑘 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.

Acknowledgments

The author is grateful to Caihua Chen, Wenxing Zhang, and Xiangfeng Wang for interesting discussions on nearest correlation matrix problem. Xiaoling Fu was supported by the NSFC Grant 70901018.