A parameterized proximal point algorithm for separable convex optimization

TL;DR

This paper introduces a parameterized proximal point algorithm designed for separable convex optimization problems, demonstrating global convergence and superior performance over existing methods in numerical experiments.

Contribution

The paper proposes a novel parameterized proximal point algorithm with proven convergence and improved efficiency for solving separable convex programming problems.

Findings

Global convergence with O(1/t) rate established

Numerical results show better performance than ADMM and relaxed PPA

Effective for sparse optimization in statistical learning

Abstract

In this paper, we develop a parameterized proximal point algorithm (P-PPA) for solving a class of separable convex programming problems subject to linear and convex constraints. The proposed algorithm is provable to be globally convergent with a worst-case O(1/t) convergence rate, wheret denotes the iteration number. By properly choosing the algorithm parameters, numerical experiments on solving a sparse optimization problem arising from statistical learning show that our P-PPA could perform significantly better than other state-of-the-art methods, such as the alternating direction method of multipliers and the relaxed proximal point algorithm.

Equations

$$\begin{array}[]{lll}\min&f(x)+g(y)\\ \textrm{s.t. }&Ax+By=c,\\ &x\in\mathcal{X},y\in\mathcal{Y},\end{array}$$

$$x^{\ast }\in {\Omega }^{\ast }\text{mbox}ifandonlyif{x}^{\ast }\in \Omega ,\phi (x)-\phi (x^{\ast })+⟨x-x^{\ast },\nabla \psi (x^{\ast })⟩\ge 0,\forall x\in \Omega .$$

$$L(x,y,\lambda )=f(x)+g(y)-⟨\lambda ,\tau (Ax+By-c)⟩,$$

$$L(x^{\ast },y^{\ast },\lambda )\le L(x^{\ast },y^{\ast },{\lambda }^{\ast })\le L(x,y,{\lambda }^{\ast }),$$

$$\left\{\begin{array}[]{l}x^{*}=\arg\min\{f(x)-\langle\lambda,\tau Ax\rangle|x\in\mathcal{X}\},\\ y^{*}=\arg\min\{g(y)-\langle\lambda,\tau By\rangle|\ y\in\mathcal{Y}\},\\ \lambda^{*}=\arg\max\{-\left\langle\lambda,\tau(Ax+By-c)\right\rangle|\ \lambda\in\mathcal{R}^{l}\}.\end{array}\right.$$

$$\left\{\begin{array}[]{lll}x^{*}\in\mathcal{X},&f(x)-f(x^{*})+\left\langle x-x^{*},-\tau A^{T}\lambda^{*}\right\rangle\geq 0,&x\in\mathcal{X},\\ y^{*}\in\mathcal{Y},&g(y)-g(y^{*})+\left\langle y-y^{*},-\tau B^{T}\lambda^{*}\right\rangle\geq 0,&y\in\mathcal{Y},\\ \lambda^{*}\in\mathcal{R}^{l},&\left\langle\lambda-\lambda^{*},\tau(Ax^{*}+By^{*}-c)\right\rangle\geq 0,&\lambda\in\mathcal{R}^{l},\end{array}\right.$$

$$\text{VI}(\varphi ,\mathcal{J},\mathcal{M}):\varphi (u)-\varphi (u^{\ast })+⟨w-w^{\ast },\mathcal{J}(w^{\ast })⟩\ge 0,\forall w\in \mathcal{M},$$

$$\varphi (u)=f(x)+g(y),\mathcal{M}=\mathcal{X}\times \mathcal{Y}\times {\mathcal{R}}^l,$$

$$u=\left(\begin{array}[]{c}x\\ y\\ \end{array}\right),\ w=\left(\begin{array}[]{c}x\\ y\\ \lambda\end{array}\right)\ \mbox{and}\ \mathcal{J}(w)=\tau\left(\begin{array}[]{c}-A^{T}\lambda\\ -B^{T}\lambda\\ Ax+By-c\end{array}\right).$$

$$⟨w-\hat{w},\mathcal{J}(w)⟩=⟨w-\hat{w},\mathcal{J}(\hat{w})⟩,\forall w,\hat{w}\in \mathcal{M}.$$

$$\varphi (u)-\varphi (u^{\ast })+⟨w-w^{\ast },\mathcal{J}(w)⟩\ge 0,\forall w\in \mathcal{M}.$$

$$\varphi (u)-\varphi (u^{k+1})+⟨w-w^{k+1},\mathcal{J}(w^{k+1})+G\left(w^{k+1}-w^k\right)⟩\ge 0,\forall w\in \mathcal{M},$$

$$u^{k+1}=\left(\begin{array}[]{c}x^{k+1}\\ y^{k+1}\\ \end{array}\right),\ w^{k+1}=\left(\begin{array}[]{c}x^{k+1}\\ y^{k+1}\\ \lambda^{k+1}\end{array}\right),\ \mathcal{J}(w^{k+1})=\tau\left(\begin{array}[]{c}-A^{T}\lambda^{k+1}\\ -B^{T}\lambda^{k+1}\\ Ax^{k+1}+By^{k+1}-c\end{array}\right).$$

$$G=\left[\begin{array}[]{cc|c}\left(\sigma+\frac{\varepsilon^{2}-1}{s}\right)A^{T}A&&-\varepsilon A^{T}\\ &\left(\rho+\frac{\tau^{2}-1}{s}\right)B^{T}B&-\tau B^{T}\\ \hline\cr-\varepsilon A&-\tau B&s\textbf{I}\end{array}\right],$$

$$s>0,\sigma >\frac{1}{s},(\sigma s-1)(\rho s-1)-{\tau }^2{\epsilon }^2>0,\epsilon \in \mathcal{R}\text{mbox}and\tau \in \mathcal{R}\backslash \{0\}.$$

$$\bar{\sigma }=\sigma +\frac{{\tau }^2-1}{s}\text{mbox}and\bar{\rho }=\rho +\frac{{\tau }^2-1}{s}.$$

$$G=D^T{G}_{0}D,$$

$$D=\left[\begin{array}[]{ccc}A&&\\ &B&\\ &&\textbf{I}\end{array}\right]\ \textrm{and}\quad G_{0}=\left[\begin{array}[]{cc|c}\left(\sigma+\frac{\varepsilon^{2}-1}{s}\right)\textbf{I}&&-\varepsilon\textbf{I}\\ &\left(\rho+\frac{\tau^{2}-1}{s}\right)\textbf{I}&-\tau\textbf{I}\\ \hline\cr-\varepsilon\textbf{I}&-\tau\textbf{I}&s\textbf{I}\end{array}\right].$$

$$\left[\begin{array}[]{ccc}\textbf{I}&&\frac{\varepsilon}{s}\textbf{I}\\ &\textbf{I}&\frac{\tau}{s}\textbf{I}\\ &&\textbf{I}\end{array}\right]\ G_{0}\ \left[\begin{array}[]{ccc}\textbf{I}&&\frac{\varepsilon}{s}\textbf{I}\\ &\textbf{I}&\frac{\tau}{s}\textbf{I}\\ &&\textbf{I}\end{array}\right]^{T}\\ =\left[\begin{array}[]{cc|c}\left(\sigma-\frac{1}{s}\right)\textbf{I}&-\frac{\tau\epsilon}{s}\textbf{I}&\\ -\frac{\tau\epsilon}{s}\textbf{I}&\left(\rho-\frac{1}{s}\right)\textbf{I}&\\ \hline\cr&&s\textbf{I}\end{array}\right]=:\widetilde{G}_{0}.$$

$$⟨\lambda -{\lambda }^{k+1},R_{\lambda }⟩\ge 0,\forall \lambda \in {\mathcal{R}}^l,$$

$$R_{\lambda }$$

$${\lambda }^{k+1}={\lambda }^k-\frac{1}{s}\left[(\tau -\epsilon )A{x}^{k+1}+\epsilon A{x}^k+\tau B{y}^k-\tau c\right].$$

$$x^{k+1}\in \mathcal{X},f(x)-f(x^{k+1})+⟨x-x^{k+1},R_x⟩\ge 0,\forall x\in \mathcal{X},$$

$$R_x$$

$${\bar{\lambda }}^k={\lambda }^k-\frac{\tau +\epsilon }{s}\left(A{x}^k+B{y}^k-c\right).$$

$$x^{k+1}$$

$$y^{k+1}\in \mathcal{Y},g(y)-g(y^{k+1})+⟨y-y^{k+1},R_y⟩\ge 0,\forall y\in \mathcal{Y},$$

$$R_y$$

$${\bar{\lambda }}^{k+\frac{1}{2}}={\lambda }^k-\frac{2}{s}\left[(\tau -\epsilon )A{x}^{k+1}+\epsilon A{x}^k+\tau B{y}^k-\tau c\right].$$

$${\bar{\lambda }}^{k+\frac{1}{2}}$$

Full text

1 Introduction

This article aims to study a novel parameterized proximal point method for solving the following two-block separable convex optimization problem

[TABLE]

where $f(x):\mathcal{R}^{m}\rightarrow\mathcal{R}$ and $\ g(y):\mathcal{R}^{n}\rightarrow\mathcal{R}$ are closed and proper convex functions, but not necessarily smooth; $A\in\mathcal{R}^{l\times m},B\in\mathcal{R}^{l\times n}$ and $c\in\mathcal{R}^{l}$ are given matrices and vectors, respectively; $\mathcal{X}\subset\mathcal{R}^{m}$ and $\mathcal{Y}\subset\mathcal{R}^{n}$ are closed convex sets. Throughout the paper, we assume the solution set of the problem (1) is nonempty and the matrices $A$ and $B$ have full column rank.

It is well-known in the literature that proximal point methods are a class of benchmark methods for solving the problem (1). The Proximal Point Algorithm (PPA) was originally proposed for solving monotone operator inclusion problems Moreau1965 ; Martinet1970 and then became popularized to convex programmings by Rockafellar Rockafellar1976 and EcksteinEckstein1993 . As demonstrated in Rockafellar1976 , the augmented Lagrangian method Powell1969 for solving the problem (1) is actually an application of PPA to its dual problem. And the recently very popular Alternating Direction Method of Multipliers (ADMM) can be also regarded as another special variant of PPA to the dual problem EcksteinBertsekas1992 . Due to its simplicity for implementation, efficiency and strong theoretical background, the PPA has attracted extensive researches in recent years for solving structured convex optimization problems, especially by the optimizers from the areas involving lots of structured data, such as compressed sensing, image processing and machine learning, etc.

There is very rich literature on PPA. Combettes and Pennanen CombettesPennanen2004 showed some conditions for the viability and weak convergence of an inexact Relaxed PPA (R-PPA) for finding a common zero of countably many cohypomonotone operators in Hilbert space. Later, based on the closed-form expressions for the proximity operators CombettesPesquet2007 , Combettes et al.CombettesPustelnik2011 still derived expressions of new proximity operators in product spaces and presented an extension of PPA for solving the multicomponent signal/image processing problems. Recently, PPA was extensively studied for a class of multi-criteria optimization problems with the difference of convex objective functions, whose efficiency was demonstrated by testing a multi-period portfolio minimization problem, see JiGoh2016 for more details. More recently, by using proximal regularization techniques and partially parallel splitting schemes, Wang et al.WangHe2017 developed a proximal partially parallel splitting method for a multi-objective convex minimization problem. For an extensive review on PPA, one may refer GuHeYuan2014 ; HeYuanZhang2013 and the references therein. In what follows, we simply mention a few works closely related to the development of our proposed method. He et al.HeYuanZhang2013 investigated a customized application of the classical PPA for the convex programming with linear constraints, where some image processing problems were tested to show the efficiency of their method. Cai et al.CaiGuHeYuan2013 also proposed a R-PPA for solving (1) and analyzed its global convergence with a worst-case linear convergence rate. Based on the results of CaiGuHeYuan2013 ; GuHeYuan2014 ; HeYuanZhang2013 , Ma and Ni MaNi2016 recently revisited the application of PPA for solving the basis pursuit and matrix completion problem. Our new proposed Parameterized PPA (P-PPA) can be actually regarded as more general extensions of the algorithms developed in HeYuanZhang2013 and MaNi2016 which does not make use of the separable structure of the objective function in (1).

Major contributions of this paper are summarized in the following. Firstly, the proximal matrix in our proposed P-PPA is more general and flexible than those in the previous work HeYuanZhang2013 ; MaNi2016 , due to more induced parameters to take consideration of the problem structure instead of a unique objective function. Secondly, by properly choosing the algorithm parameters, the new P-PPA could significantly outperform some state-of-the-art methods, such as ADMM BoydChu2010 and R-PPA GuHeYuan2014 , for solving the separable convex optimization, especially when the problem size is large and high accurate solutions are required.

The remaining parts are organized as follows. In Section 2, we characterize the solution of the problem (1) as the solution of proper variational inequalities and review the unified framework of PPA. In Section 3, we derive the new P-PPA and discuss its global convergence and worst-case convergence rate in an ergodic sense. At the end of Section 3, we still present a P-PPA with a relaxation step. Some preliminary numerical experiments are performed in Section 4 for comparing our proposed methods with two benchmark methods. We finally conclude the paper in Section 5.

2 Preliminaries

In this section, we first introduce some necessary notations used throughout the paper. Then, we characterize the solution of the problem (1) by the aid of an equivalent variational inequality. Similar approaches have been widely used in the literature, e.g. HeYuanZhang2013 ; HeMaYuan2016 .

For the sake of convenience, let $\mathcal{R}$, $\mathcal{R}^{n}$ and $\mathcal{R}^{n\times m}$ denote the set of real numbers, the set of $n$ dimensional real column vectors and the set of $n\times m$ real matrices, respectively. For any $x,y\in\mathcal{R}^{n}$, $\langle x,y\rangle=x^{T}y$ denotes the standard inner product in $\mathcal{R}^{n}$ and $\|x\|=\sqrt{\langle x,x\rangle}$ is the Euclidean norm. Given any symmetric positive definite matrix $G\in\mathcal{R}^{n\times n}$, the weighted norm $\|x\|_{G}=\sqrt{\langle x,Gx\rangle}$. In addition, we use I and 0 to stand for the identity matrix and the zero vector with proper dimension, respectively.

The following basic lemma given in HeMaYuan2016 will be used as a tool for analyzing the primal-dual solution pair of the problem (1).

Lemma 1

Let $\varphi:\mathcal{R}^{m}\rightarrow\mathcal{R}$ and $\psi:\mathcal{R}^{m}\rightarrow\mathcal{R}$ be two convex functions defined on a closed convex set $\Omega\subset\mathcal{R}^{m}$ and $\psi$ is differentiable. Suppose the solution set $\Omega^{*}=\arg\min\limits_{x\in\Omega}\{\varphi(x)+\psi(x)\}$ is nonempty. Then we have

[TABLE]

Now, given any $\tau\in\mathcal{R}\backslash\{0\}$, that is $\tau\neq 0$, a Lagrangian function of problem (1) can be written as

[TABLE]

where $\lambda\in\mathcal{R}^{l}$ is the Lagrange multiplier. Then, for any primal-dual solution pair $(x^{*},y^{*},\lambda^{*})$ of (1), we have

[TABLE]

which is equivalent to

[TABLE]

By applying Lemma 1, the optimality conditions of the above equations are

[TABLE]

which can be rewritten as a variational inequality (VI)

[TABLE]

where

[TABLE]

[TABLE]

Clearly, the solution set of $\textrm{VI}(\phi,\mathcal{J},\mathcal{M})$, denoted by $\mathcal{M}^{*}$, is nonempty by the assumption of nonempty solution set of the problem (1). Since the affine mapping $\mathcal{J}$ is skew-symmetric, we can obtain

[TABLE]

Hence, the variational inequality (4) is also rewritten as

[TABLE]

When the proximal point algorithms are applied to solve the variational inequality (6) or equivalently (4), they often take the following unified approach: at the $k$-th iteration, find iterate $w^{k+1}$ satisfying

[TABLE]

where the above matrix $G\in\mathcal{R}^{(m+n+l)\times(m+n+l)}$ called proximal matrix is a positive definite, and

[TABLE]

Obviously, different choices of $G$ would result in different proximal point algorithms. We would provide a new choice of $G$ for our proposed P-PPA.

3 Main results

In this section, we first develop the P-PPA for solving (1) in detail and discuss its convergence properties. Then, it is straightforward to extend the method to the case with a relaxation step.

3.1 Development of P-PPA with convergence

Mainly motivated by the proximal matrix of PPA in Eq.(2.5) of HeYuanZhang2013 and Eq.(3.1) of MaNi2016 , we would design the matrix $G$ in (7) having the following structure:

[TABLE]

where $(\sigma,\rho,s,\tau,\varepsilon)$ are parameters satisfying

[TABLE]

For convenience, let us define

[TABLE]

Then, from later analysis we can see that ${1}/{s}$ would play a role of penalty parameter for the equality constraint of (1), while $\bar{\sigma}$ and $\bar{\rho}$ can be regarded as the proximal parameters as those used in the customized PPA HeYuanZhang2013 .

The following lemma ensures that under proper conditions of the parameters, $G$ is a positive definite matrix.

Lemma 2

Suppose that the matrices $A$ and $B$ have full column rank. For any $(\sigma,\rho,s,\tau,\varepsilon)$ satisfying (9), the matrix $G$ defined in (8) is positive definite.

Proof

Clearly, the matrix $G$ is symmetric and can be decomposed into

[TABLE]

where

[TABLE]

By the full column rank assumption of $A$ and $B$, the matrix $G$ is positive definite if and only if $G_{0}$ is positive definite. Noting that

[TABLE]

Hence, $G_{0}$ is positive definite if and only if $\widetilde{G}_{0}$ is positive definite, which is guaranteed if condition (9) holds. Therefore, the proof is completed. $\ \ \ \diamondsuit$

In what follows, we develop our P-PPA in detail. Substituting the matrix $G$ into (7), we have $\lambda^{k+1}\in\mathcal{R}^{l}$ and

[TABLE]

where

[TABLE]

Hence, we have $R_{\lambda}=0$ which leads to

[TABLE]

Meanwhile, by (7) and (12), we also have

[TABLE]

where

[TABLE]

with $\bar{\sigma}$ being defined in (10) and

[TABLE]

Notice that by (9), we get $\bar{\sigma}=(\sigma s+\tau^{2}-1)/s>0$. Hence, by Lemma 1, $x^{k+1}$ is the solution of the following optimization problem

[TABLE]

Similarly, we can also derive from (7) and (12) that

[TABLE]

where

[TABLE]

and

[TABLE]

Furthermore, it follows from (14) and (17) that

[TABLE]

Hence, by Lemma 1, $y^{k+1}$ is the solution of the following optimization problem

[TABLE]

In addition, it follows from (12) and (14) that

[TABLE]

Summarizing all the above discussions, we propose P-PPA as the following algorithm:

Algorithm 1 (P-PPA for solving Problem (1))  1  Choose parameters $(\sigma,\rho,s,\tau,\varepsilon)$ satisfying (9);

2  Initialize $(x^{0},y^{0},\lambda^{0})\in\mathcal{R}^{m}\times\mathcal{R}^{n}\times\mathcal{R}^{l}$ and $r^{0}=Ax^{0}+By^{0}-c$;

3  $\bar{\lambda}^{0}=\lambda^{0}-\frac{\tau+\varepsilon}{s}r^{0}$;

4  For $k=0,1,\cdots,$ **do

**5   $x^{k+1}=\arg\min\limits_{x\in\mathcal{X}}\left\{f(x)+\frac{\bar{\sigma}}{2}\left\|A(x-x^{k})-\frac{\tau}{\bar{\sigma}}\bar{\lambda}^{k}\right\|^{2}\right\}$;

6   $\bar{\lambda}^{k+\frac{1}{2}}=\bar{\lambda}^{k}-\frac{\tau-\varepsilon}{s}\left[A(2x^{k+1}-x^{k})+By^{k}-c\right]$;

7   $y^{k+1}=\arg\min\limits_{y\in\mathcal{Y}}\left\{g(y)+\frac{\bar{\rho}}{2}\left\|B(y-y^{k})-\frac{\tau}{\bar{\rho}}\bar{\lambda}^{k+\frac{1}{2}}\right\|^{2}\right\}$;

8   $r^{k+1}=Ax^{k+1}+By^{k+1}-c$;

9   $\bar{\lambda}^{k+1}=\bar{\lambda}^{k}-\frac{\tau}{s}r^{k+1}-\frac{1}{s}\left[\tau A(x^{k+1}-x^{k})+\varepsilon B(y^{k+1}-y^{k})\right]$.

For the above P-PPA, we have the following remarks.

Remark 1

By taking $\tau=\varepsilon=1,$ the matrix $G$ in (8) would become

[TABLE]

where the parameters $(\sigma,\rho,s)$ satisfy

[TABLE]

In such case, Algorithm 1 would be reduced to

[TABLE]

This algorithm can be considered as a direct extension of the customized PPA HeYuanZhang2013 , which only considers problem with one block structure, to solve the two-block structured problem (1).

Remark 2

The freedom of choosing the parameters $(\sigma,\rho,s,\tau,\varepsilon)$ would allow P-PPA to have more flexibility to select the proximal parameters $\bar{\sigma}$ and $\bar{\rho}$. However, note that they are proportional to the parameters $(\sigma,\tau,1/s)$ and $(\rho,\tau,1/s)$, respectively. As commented in the final part of HeXuYuan2016 , if they are too large, then slow convergence will occur in terms of solving the subproblems, which can significantly affect the overall efficiency of the algorithm.

Next, we discuss the convergence properties of P-PPA in a more general proximal point setting given in (7).

Lemma 3

The sequence $\{w^{k+1}\}$ generated by Algorithm 1 satisfies

[TABLE]

Proof

Setting $w=w^{*}$ in (7), we have

[TABLE]

which, by (4)-(5), implies

[TABLE]

By the above inequality, we obtain

[TABLE]

which immediately gives the inequality (21). $\ \ \ \diamondsuit$

Based on the above lemma, we have the following global convergence theorem.

Theorem 3.1

Suppose that the condition (9) holds and the sequence $\{w^{k+1}\}$ is generated by Algorithm 1. Then, there exists a $w^{\infty}\in\mathcal{M}^{*}$ such that

[TABLE]

Proof

Since condition (9) holds, we have by Lemma 2 that $G$ is positive definite. Then, it follows from (21) that $\{w^{k}\}$ is bounded and

[TABLE]

Let $w^{\infty}$ be any accumulation point of $\{w^{k}\}$. By taking a subsequence of $w^{k}$ in (7) if necessary, it follows from (23) that

[TABLE]

Hence, $w^{\infty}\in\mathcal{M}^{*}$. So, by (21) again, we have

[TABLE]

Then, it follows from $w^{\infty}$ being an accumulation point that (22) holds. $\ \ \ \diamondsuit$

Now, we establish the worst-case $\mathcal{O}(1/t)$ ergodic convergence rate for solving the variational inequality (6). Let

[TABLE]

Theorem 3.2

Suppose that the condition (9) holds and the sequence $\{w^{k+1}\}$ is generated by Algorithm 1. Then, we have

[TABLE]

Proof

By (7) and the property (5), it holds that

[TABLE]

Then, applying the identity

[TABLE]

with

[TABLE]

we can obtain

[TABLE]

which together with (26) imply

[TABLE]

Summing the above inequality over $k=0,1,\cdots,t$, we get

[TABLE]

which by the definition of $\textbf{w}_{t}$ in (24) gives

[TABLE]

By the convexity of the function $\phi(u)$ and the definition of $\textbf{u}_{t}$ in (24), we have

[TABLE]

Then, (25) is true by substituting the above inequality into (27). $\ \ \ \diamondsuit$

3.2 P-PPA with a relaxation step

Applying a relaxation step for PPA is a standard technique to accelerate its convergence GuHeYuan2014 ; Golsshtein1979 ; HeYuanZhang2013 . Combining the PPA approach (7) together with a relaxation step would give the following procedure: at the $k$-th iteration, find $\widetilde{w}^{k+1}$ satisfying

[TABLE]

and then let

[TABLE]

where $\gamma\in(0,2)$ is the relaxation parameter. When $\gamma=1$, the above relaxed PPA will be reduced to the standard PPA. For the relaxed PPA, analogous to (21), it is not difficult (for details, see GuHeYuan2014 ; Golsshtein1979 ) to show

[TABLE]

Then, based on (29), it is straightforward to have global convergence and convergence rate analogous to Theorem 3.1 and Theorem 3.2. More precisely, combining Algorithm 1 with a relaxation step would give the following Relaxed P-PPA (RP-PPA).

Algorithm 2 (RP-PPA for solving Problem (1))  1  Choose parameters $(\sigma,\rho,s,\tau,\varepsilon)$ satisfying (9) and $\gamma\in(0,2)$;

2  Initialize $(x^{0},y^{0},\lambda^{0})\in\mathcal{R}^{m}\times\mathcal{R}^{n}\times\mathcal{R}^{l}$ and $r^{0}=Ax^{0}+By^{0}-c$;

3  $\bar{\lambda}^{0}=\lambda^{0}-\frac{\tau+\varepsilon}{s}r^{0}$;

4  For $k=0,1,\cdots,$ **do

**5   $\widetilde{x}^{k}=\arg\min\limits_{x\in\mathcal{X}}\left\{f(x)+\frac{\bar{\sigma}}{2}\left\|A(x-x^{k})-\frac{\tau}{\bar{\sigma}}\bar{\lambda}^{k}\right\|^{2}\right\}$;

6   $\bar{\lambda}^{k+\frac{1}{2}}=\bar{\lambda}^{k}-\frac{\tau-\varepsilon}{s}\left[A(2\widetilde{x}^{k}-x^{k})+By^{k}-c\right]$;

7   $\widetilde{y}^{k}=\arg\min\limits_{y\in\mathcal{Y}}\left\{g(y)+\frac{\bar{\rho}}{2}\left\|B(y-y^{k})-\frac{\tau}{\bar{\rho}}\bar{\lambda}^{k+\frac{1}{2}}\right\|^{2}\right\}$;

8   $r^{k+1}=A\widetilde{x}^{k}+B\widetilde{y}^{k}-c;\ \Delta x^{k}=\widetilde{x}^{k}-x^{k};\ \Delta y^{k}=\widetilde{y}^{k}-y^{k}$;

9   $\widetilde{\lambda}^{k}=\bar{\lambda}^{k}-\frac{\tau}{s}r^{k+1}-\frac{1}{s}\left[\tau A\Delta x^{k}+\varepsilon B\Delta y^{k}\right]$;

10    Relaxation step: compute

[TABLE]

and

[TABLE]

4 Numerical experiments

In this section, we would perform some numerical experiments for solving the following lasso model problem arising from statistical learning Tibshirani1996 :

[TABLE]

where $\|x\|_{1}=\sum_{i=1}^{n}|x_{i}|$; $\nu>0$ is a scalar regularization parameter; $b\in\mathcal{R}^{l}$ is a response vector; $D\in\mathcal{R}^{l\times n}$ is a design matrix with $l$ and $n$ denoting the number of data points and the number of features, respectively. In typical applications, there are usually many more features than training examples BoydChu2010 , that is $l\ll n$, and the goal is to find a parsimonious model for the data. We refer the readers to ChenDSa2001 ; Tibshirani1996 ; BoydChu2010 ; Hastie2009 for more backgrounds on the lasso model. In the numerical experiments, all algorithms are coded in MATLAB 7.10(R2010a) and run on a PC with Intel Core i5 processor (3.3GHz) with 4 GB memory.

By introducing an auxiliary variable $y\in\mathcal{R}^{n}$, the problem (30) is obviously equivalent to

[TABLE]

which is a special case of (1) with

[TABLE]

Applying Algorithm 1 to solve (31), we have

[TABLE]

for which $x^{k+1}$ can be explicitly obtained by a soft shrinkage operator Donoho Tsaig2008 and $\bar{\sigma}$ is defined in (10). In addition, we can deduce that

[TABLE]

where $\bar{\rho}$ is defined in (10). Notice that the matrix $D^{T}D+\bar{\rho}\textbf{I}$ is positive definite, since $\bar{\rho}=(\rho s+\tau^{2}-1)/s>0$. Though it maybe quite time consuming to compute $(D^{T}D+\bar{\rho}\textbf{I})^{-1}$ and $D^{T}b$ when the problem scale is large, they only need to be computed once before the iteration starts. To save computation, we can actually compute once and cache the Cholesky Factorization of the much smaller matrix $DD^{T}/\bar{\rho}+\textbf{I}$ (note $l\ll n$), which takes about $(l^{2}n+l^{3}/3)$ flops, including the cost of forming $DD^{T}$ and the Cholesky Factorization. Then, all the subsequent $y$-updates can be calculated by the Sherman-Morrison inversion matrix formula together with forward-backward substitutions.

The problem data are generated by the following way. Each entry of the feature matrix $D$ is draw from the standard normal distribution $\mathcal{N}(0,1)$ and then each of its column is normalized. A random sparse vector $x^{true}\in\mathcal{R}^{n}$ with 100 nonzero entries is also drawn from an $\mathcal{N}(0,1)$ distribution. The vector $b$ is computed via $b=Dx^{true}+\epsilon$, where $\epsilon\sim\mathcal{N}(0,10^{-3}\textbf{I})$, and the regularization parameter is set as $\nu=0.12\nu_{max}$ with $\nu_{max}=\|D^{T}b\|_{\infty}$. The following stopping criterion is used for all the comparison algorithms

[TABLE]

where Tol is a given tolerance, $\phi(u^{k})=f(x^{k})+g(y^{k})$, and $\phi^{*}$ is the approximate optimal objective function value obtained by running P-PPA after $2000$ iterations. Then, all the comparison algorithms are set to have the maximum $2000$ number of iterations and they all use the same starting point $(x^{0},y^{0},\lambda^{0})=(\textbf{0},\textbf{0},\textbf{0})$.

4.1 Effects of parameters $(\sigma,\rho,s,\tau,\varepsilon)$

The aim of this subsection is to investigate how the five parameters $(\sigma,\rho,s,\tau,\varepsilon)$ would influence the performance of P-PPA. For this purpose, we first fix the free parameters $(\tau,\varepsilon)$ as $(3,1.5)$ and then change other distinctively constrained parameters $(\sigma,\rho,s)$ to investigate their effects on P-PPA.

Table 1 presents the numerical results of Algorithm 1 (i.e. P-PPA) with different parameters for solving the test problem (31) with dimension $(l,n)=(1800,4000)$. For this set of tests, we fix the tolerance $\textrm{Tol}=1.0\times 10^{-6}$. And in all the numerical Tables, “Iter”, “CPU” and “DRN” denote the iteration numbers, the CPU time in seconds and the dual residual norm $\|y^{k}-y^{k-1}\|$ , respectively. We can observe from Table 1 that:

• •

For the parameters $(\sigma,\rho,s)$, the reported results in each column of IRE, DRN and $\phi(u^{k})$ are nearly the same when fixed any two parameters with one parameter changing.

• •

With the increase of the parameter $\sigma$ or $\rho$, both the iteration number and the CPU time tend to increase(denoted by $\downarrow$);

• •

With the increase of the parameter $s$, both the iteration number and the CPU time decrease firstly and then increase(denoted by ${}^{\uparrow}_{\downarrow}$).

These changing trends identify with Remark 2. Reported results of Table 1 indicate that the choice of the parameters $(\sigma,\rho,s)$ could have a great effect on the performance of P-PPA, the value in the place of the arrow is better than others in each subtable, and it seems that setting $(\sigma,\rho,s)=(0.8,6,3)$ would be a reasonable choice for solving the test problem (31).

Remark 3

Noting from Table 1 that if any four parameters are fixed, then by (9) the remaining one is clearly subjected to a given domain. For instance, in the top subtable of Table 1 we have from (9) that $\sigma>\frac{37.25}{51}\approx 0.73.$ Therefore, we can randomly choose some values in such region to do experiments to find out which one gives relatively better performance.

[TABLE]

Table 1: Results of problem (31) by P-PPA with different parameters $(\sigma,\rho,s)$.

Remark 4

In a similar way as mentioned in Remark 3, we have

[TABLE]

Then, by testing some values of the parameter $\tau(\varepsilon)$ and observing which one performs approximately better, our tuned results are $(\tau,\varepsilon)=(3,1.5).$ Hence, for further comparative experiments with some state-of-the-art methods, we would use the tuned results $(0.8,6,3,3,1.5)$ as the default parameter setting for both P-PPA and RP-PPA.

4.2 Comparative experiments

Now, we would like to compare P-PPA and RP-PPA with other two popular methods for solving the problem (31): ADMM111Available at http://web.stanford.edu/$\sim$boyd/papers/admm/. BoydChu2010 and R-PPA GuHeYuan2014 .

[TABLE]

Table 2: Comparative results of problem (31) with different dimensions222The bold value excluding that of P-PPA is of the smallest in each experiment with respect to (l,n), and the bold value of P-PPA is smaller than that of R-PPA and ADMM..

Table 2 reports the numerical results of all comparison methods for solving the problem (31) with different dimensions $(l,n)$. Actually, ADMM uses the downloaded codes but with penalty parameter 1 and a widely used step-length 1.618 when updating the lagrangian multipliers. For R-PPA, the penalty parameter is set as $10$, which is reasonably good for this method. Both RP-PPA and R-PPA use the same relaxation factor $\gamma=1.2$. The tolerance $\textrm{Tol}=1.0\times 10^{-10}$ is used for all the testing problems in Table 2. The numerical results of solving the test problem with fixed dimensions $(l,n)=(1800,20000)$, but under different accurate tolerances, are presented in Table 3. Moreover, the comparative convergence curves of the objective function $\phi(u^{k})$, the residual error IRE(k) and the dual residual norm $\|y^{k}-y^{k-1}\|$ against the number of iterations are shown in Fig. 1.

[TABLE]

Table 3: Comparative results of problem (31) under different tolerance errors333The bold value of P-PPA is smaller than that of R-PPA and ADMM, and the bold value of RP-PPA is of the smallest..

We can observe from Tables 2-3 that both P-PPA and RP-PPA perform better than ADMM and R-PPA for the relatively large size problems in terms of both the number of iterations and the CPU time. Another outstanding observation is that RP-PPA can clearly shorten the number of iterations and the CPU time of P-PPA. Besides, ADMM performs the worst for large-size problems, while it performs the best for small-size problems (e.g. $n=4000$). Both Table 3 and Fig. 1 illustrate that as the tolerance becomes smaller, P-PPA and RP-PPA could perform significantly better than ADMM and R-PPA. In addition, note from Table 3 that ADMM fails to solve the problem with dimensions $(l,n)=(1800,20000)$ in $2000$ iterations to achieve the accuracy $\textrm{Tol}=1.0\times 10^{-14}$. Reported results show that our proposed methods are efficient when properly choosing the algorithmic parameters.

Remark 5

Although the tuned values of the parameters in the proposed algorithms are not proved theoretically to be the best, the reported numerical results of comparative experiments are sufficient to show that such a choice can make our algorithms outperform the other two algorithms.

5 Conclusion and discussion

By introducing several parameters to the proximal matrix in the framework of the traditional proximal point algorithm, we propose a new Parameterized Proximal Point Algorithm (P-PPA) for solving the separable convex minimization problem. Under certain conditions on these parameters, we show that the P-PPA is globally convergent and would maintain a worst-case $O(1/t)$ ergodic convergence rate. By properly choosing the parameters, the numerical experiments of solving the classical lasso problem in statistical learning indicate our P-PPA and the Relaxed P-PPA (RP-PPA) could perform significantly better than the other two benchmark methods: ADMM and R-PPA, especially for solving large scale problems and high accurate solutions are required.

For the case that the subproblems are not easy to solve, inexact ADMMsHagerZhang2016 are recently developed for solving the general separable convex optimization problems with a linear constraint and with an objective including smooth plus nonsmooth terms, which is particularly useful when the ADMM’s subproblems do not have closed solutions or when the solution of the subproblem is expensive. Also, there are other works in references He2013 ; Solodov2000 which discuss how to solve the subproblems inexactly.

Finally, observe that the P-PPA and RP-PPA can be naturally extended to to solve the problem (1) with inequality constraints or with matrix variables such as

[TABLE]

where both $f$ and $g$ are proper closed convex functions over the matrix variables $X$ and $Y$, $A$ and $B$ are coefficient matrices, $\mathcal{X}$ and $\mathcal{Y}$ are certain closed convex sets. The model problem (32) also arises very often in many important applications in data analysis CandWright2011 ; LiuLiBaiLiu2017 , for example, the robust principal component analysis in image processing, etc. One of the following research tasks could be to extend the P-PPA to solve the multi-block separable convex/nonconvex programming problems.

Acknowledgements

The authors wish to thank the Editor-in-Chief Prof. O.A. Prokopyev and the anonymous referees for providing their valuable suggestions which have significantly improved the quality of the paper.