General parameterized proximal point algorithm with applications in statistical learning

TL;DR

This paper introduces a flexible, parameterized proximal point algorithm with a relaxation step for multi-block convex optimization, demonstrating improved efficiency and convergence properties in statistical learning applications.

Contribution

It develops a general parameterized PPA with relaxation for multi-block convex problems, establishing convergence and efficiency improvements over existing methods.

Findings

Proposed algorithm has a global convergence guarantee.

Achieves an $ ext{O}(1/t)$ convergence rate.

Numerical experiments show superior efficiency in sparse matrix minimization.

Abstract

In the literature, there are a few researches to design some parameters in the Proximal Point Algorithm (PPA), especially for the multi-objective convex optimizations. Introducing some parameters to PPA can make it more flexible and attractive. Mainly motivated by our recent work (Bai et al., A parameterized proximal point algorithm for separable convex optimization, Optim. Lett. (2017) doi: 10.1007/s11590-017-1195-9), in this paper we develop a general parameterized PPA with a relaxation step for solving the multi-block separable structured convex programming. By making use of the variational inequality and some mathematical identities, the global convergence and the worst-case $\mathcal{O}(1/t)$ convergence rate of the proposed algorithm are established. Preliminary numerical experiments on solving a sparse matrix minimization problem from statistical learning validate that our algorithm is more efficient than several state-of-the-art algorithms.

Tables

Parameter	GR-PPA						PJALM						HTY
${ϵ}_1({ϵ}_2={10}^{-10}$  $/{ϵ}_3={10}^{-6})$	IT	CPU	IER	OER	CER		IT	CPU	IER	OER	CER		IT	CPU	IER	OER	CER
${10}^{-4}$	110	3.15	2.55e-7	2.89e-11	2.45e-7		132	3.55	1.28e-7	9.50e-11	5.59e-8		128	3.36	1.51e-9	9.26e-11	1.73e-8
${10}^{-7}$	116	3.28	9.39e-8	7.73e-11	3.09e-7		136	3.67	9.53e-8	6.88e-11	4.17e-8		128	3.36	1.51e-9	9.26e-11	1.73e-8
${10}^{-12}$	215	6.08	9.25e-13	4.44e-16	5.49e-12		293	7.63	9.96e-13	8.87e-16	4.35e-13		234	6.21	9.56e-13	2.44e-15	1.18e-11
${ϵ}_2({ϵ}_1={10}^{-6}$  $/{ϵ}_3={10}^{-6})$
${10}^{-4}$	105	2.72	4.17e-7	1.46e-10	6.58e-7		107	3.08	9.34e-7	9.54e-10	4.06e-7		79	2.03	1.03e-7	1.35e-8	9.20e-7
${10}^{-11}$	124	3.54	5.11e-8	3.31e-12	5.86e-8		161	4.24	1.53e-8	9.76e-12	6.70e-9		151	4.01	2.83e-10	9.40e-12	3.41e-9
${10}^{-15}$	214	6.15	1.36e-12	4.44e-16	1.71e-12		279	7.49	2.74e-12	2.22e-16	1.21e-12		243	6.48	5.20e-13	8.87e-16	6.42e-12
${ϵ}_3({ϵ}_1={10}^{-6}$  $/{ϵ}_2={10}^{-8})$
${10}^{-8}$	141	4.09	6.05e-9	1.60e-12	7.80e-9		156	5.12	2.20e-8	1.43e-11	9.65e-9		136	3.89	8.32e-10	4.17e-11	9.76e-9
${10}^{-11}$	200	6.44	7.24e-12	2.44e-15	6.41e-12		251	7.19	2.12e-11	1.26e-14	9.33e-12		237	6.96	7.83e-13	3.10e-15	9.66e-12
${10}^{-12}$	225	7.31	1.50e-12	6.66e-16	9.87e-13		282	8.65	2.22e-12	3.11e-15	6.70e-13		271	8.10	9.94e-14	2.22e-16	9.59e-13

Equations

$$\begin{array}[]{lll}\min&\sum\limits_{i=1}^{p}f_{i}(x_{i})\\ \textrm{s.t. }&\sum\limits_{i=1}^{p}A_{i}x_{i}=b,\ x_{i}\in\mathcal{X}_{i},\end{array}$$

$$A=[A_1,\cdots ,A_p],x=(x_1^{\mathsf{T}},\cdots ,x_p^{\mathsf{T}}{)}^{\mathsf{T}}$$

$$L(x_1,\cdots ,x_p,\lambda )=\sum _{i=1}^p{f}_i(x_i)-\tau ⟨\lambda ,\sum _{i=1}^p{A}_i{x}_i-b⟩,$$

$$L(x_1^{\ast },\cdots ,x_p^{\ast },\lambda )\le L(x_1^{\ast },\cdots ,x_p^{\ast },{\lambda }^{\ast })\le L(x_1,\cdots ,x_p,{\lambda }^{\ast }),$$

$$\left\{\begin{array}[]{lll}x_{1}^{*}=\arg\min\limits_{x_{1}\in\mathcal{X}_{1}}\left\{f_{1}(x_{1})-\tau\left\langle\lambda^{*},A_{1}x_{1}\right\rangle\right\},\\ \ \vdots\\ x_{p}^{*}=\arg\min\limits_{x_{p}\in\mathcal{X}_{p}}\left\{f_{p}(x_{p})-\tau\left\langle\lambda^{*},A_{p}x_{p}\right\rangle\right\},\\ \lambda^{*}=\arg\max\limits_{\lambda\in\mathbb{R}^{m}}-\tau\left\langle\lambda,\sum\limits_{i=1}^{p}A_{i}x_{i}^{*}-b\right\rangle,\end{array}\right.$$

$$\left\{\begin{array}[]{lll}x_{1}^{*}\in\mathcal{X}_{1},&f_{1}(x_{1})-f_{1}(x_{1}^{*})+\left\langle x_{1}-x_{1}^{*},-\tau A_{1}^{\sf T}\lambda^{*}\right\rangle\geq 0,&\forall\ x_{1}\in\mathcal{X}_{1},\\ \quad\vdots&&\\ x_{p}^{*}\in\mathcal{X}_{p},&f_{p}(x_{p})-f_{p}(x_{p}^{*})+\left\langle x_{p}-x_{p}^{*},-\tau A_{p}^{\sf T}\lambda^{*}\right\rangle\geq 0,&\forall\ x_{p}\in\mathcal{X}_{p},\\ \lambda^{*}\in\mathbb{R}^{m},&\left\langle\lambda-\lambda^{*},\tau(\sum\limits_{i=1}^{p}A_{i}x_{i}^{*}-b)\right\rangle\geq 0,&\forall\ \lambda\in\mathbb{R}^{m}.\end{array}\right.$$

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

$$\varphi (u)=\sum _{i=1}^p{f}_i(x_i),\Omega ={\mathcal{X}}_1\times \cdots \times {\mathcal{X}}_p\times {\mathbb{R}}^m,$$

$$u=\left(\begin{array}[]{c}x_{1}\\ x_{2}\\ \vdots\\ x_{p}\\ \end{array}\right),\ w=\left(\begin{array}[]{c}x_{1}\\ \vdots\\ x_{p}\\ \lambda\end{array}\right)\ \mbox{and}\ \mathcal{J}(w)=\tau\left(\begin{array}[]{c}-A_{1}^{\sf T}\lambda\\ \vdots\\ -A_{p}^{\sf T}\lambda\\ \sum\limits_{i=1}^{p}A_{i}x_{i}-b\end{array}\right).$$

$$u^{k}=\left(\begin{array}[]{c}x_{1}^{k}\\ x_{2}^{k}\\ \vdots\\ x_{p}^{k}\\ \end{array}\right),\ w^{k}=\left(\begin{array}[]{c}x_{1}^{k}\\ \vdots\\ x_{p}^{k}\\ \lambda^{k}\end{array}\right)\ \mbox{and}\ \mathcal{J}(w^{k})=\tau\left(\begin{array}[]{c}-A_{1}^{\sf T}\lambda^{k}\\ \vdots\\ -A_{p}^{\sf T}\lambda^{k}\\ \sum\limits_{i=1}^{p}A_{i}x_{i}^{k}-b\end{array}\right).$$

$$⟨w^{k+1}-w,\mathcal{J}(w^{k+1})⟩=⟨w^{k+1}-w,\mathcal{J}(w)⟩,\forall w,w^{k+1}\in \Omega .$$

$$w^{k+1}\in \Omega ,\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 \Omega ,$$

$$G=\left[\begin{array}[]{cccc|c}\left(\sigma_{1}+\frac{\varepsilon^{2}-1}{s}\right)A_{1}^{\sf T}A_{1}&&&&-\varepsilon A_{1}^{\sf T}\\ &\left(\sigma_{2}+\frac{\tau^{2}-1}{s}\right)A_{2}^{\sf T}A_{2}&&&-\tau A_{2}^{\sf T}\\ &&\ddots&&\vdots\\ &&&\left(\sigma_{p}+\frac{\tau^{2}-1}{s}\right)A_{p}^{\sf T}A_{p}&-\tau A_{p}^{\sf T}\\ \hline\cr-\varepsilon A_{1}&-\tau A_{2}&\cdots&-\tau A_{p}&s\textbf{I}\end{array}\right],$$

$$\mathcal{K}=\left\{{\sigma }_1>\frac{1+(p-1)\tau |\epsilon |}{s},{\sigma }_i>\frac{1+(p-2){\tau }^2+\tau |\epsilon |}{s},s>0|\epsilon \in \mathbb{R},\tau \in {\mathbb{R}}^{+},i=2,\cdots ,p\right\}.$$

$${\bar{\sigma }}_i:={\sigma }_i+\frac{{\tau }^2-1}{s},\forall i=1,\cdots ,p.$$

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

$$G_{0}=\left[\begin{array}[]{cccc|c}\left(\sigma_{1}+\frac{\varepsilon^{2}-1}{s}\right)\textbf{I}&&&&-\varepsilon\textbf{I}\\ &\left(\sigma_{2}+\frac{\tau^{2}-1}{s}\right)\textbf{I}&&&-\tau\textbf{I}\\ &&\ddots&&\vdots\\ &&&\left(\sigma_{p}+\frac{\tau^{2}-1}{s}\right)\textbf{I}&-\tau\textbf{I}\\ \hline\cr-\varepsilon\textbf{I}&-\tau\textbf{I}&\cdots&-\tau\textbf{I}&s\textbf{I}\end{array}\right].$$

$$\begin{array}[]{l}\left[\begin{array}[]{ccccc}\textbf{I}&&&&\frac{\varepsilon}{s}\textbf{I}\\ &\textbf{I}&&&\frac{\tau}{s}\textbf{I}\\ &&\ddots&&\vdots\\ &&&\textbf{I}&\frac{\tau}{s}\textbf{I}\\ &&&&\textbf{I}\end{array}\right]G_{0}\left[\begin{array}[]{ccccc}\textbf{I}&&&&\frac{\varepsilon}{s}\textbf{I}\\ &\textbf{I}&&&\frac{\tau}{s}\textbf{I}\\ &&\ddots&&\vdots\\ &&&\textbf{I}&\frac{\tau}{s}\textbf{I}\\ &&&&\textbf{I}\end{array}\right]^{\sf T}\\ \\ =\left[\begin{array}[]{ccccc|c}\left(\sigma_{1}-\frac{1}{s}\right)\textbf{I}&-\frac{\varepsilon\tau}{s}\textbf{I}&-\frac{\varepsilon\tau}{s}\textbf{I}&\cdots&-\frac{\varepsilon\tau}{s}\textbf{I}&\\ -\frac{\tau\varepsilon}{s}\textbf{I}&\left(\sigma_{2}-\frac{1}{s}\right)\textbf{I}&-\frac{\tau^{2}}{s}\textbf{I}&\cdots&-\frac{\tau^{2}}{s}\textbf{I}&\\ \vdots&\vdots&\vdots&\ddots&\vdots&\\ -\frac{\tau\varepsilon}{s}\textbf{I}&-\frac{\tau^{2}}{s}\textbf{I}&-\frac{\tau^{2}}{s}\textbf{I}&\cdots&\left(\sigma_{p}-\frac{1}{s}\right)\textbf{I}&\\ \hline\cr&&&&&s\textbf{I}\end{array}\right]=:\widetilde{G}_{0},\end{array}$$

$$G=P^{\mathsf{T}}{\bar{G}}_{0}P$$

$$\bar{G}_{0}=\left[\begin{array}[]{cccc|c}\nu\textbf{I}&-\textbf{I}&\cdots&-\textbf{I}&\textbf{0}\\ -\textbf{I}&\nu\textbf{I}&\cdots&-\textbf{I}&\textbf{0}\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ -\textbf{I}&-\textbf{I}&\cdots&\nu\textbf{I}&\textbf{0}\\ \hline\cr\textbf{0}&\textbf{0}&\cdots&\textbf{0}&\textbf{I}\end{array}\right],$$

$$(s,\tau ,\epsilon ):=(1,1,1)\&{\sigma }_i:=\nu +1,\forall i=1,2,\cdots ,p,$$

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

$$0=R_{\lambda }:=\tau \left(\sum _{i=1}^p{A}_i{x}_i^{k+1}-b\right)-\epsilon A_1(x_1^{k+1}-x_1^k)-\tau \sum _{i=2}^p{A}_i(x_i^{k+1}-x_i^k)+s\left({\lambda }^{k+1}-{\lambda }^k\right).$$

$${\lambda }^{k+1}={\lambda }^k-\frac{1}{s}\left[(\tau -\epsilon )A_1{x}_1^{k+1}+\epsilon A_1{x}_1^k+\tau \sum _{i=2}^p{A}_i{x}_i^k-b\tau \right].$$

$$x_1^{k+1}\in {\mathcal{X}}_1,f_1(x_1)-f_1(x_1^{k+1})+⟨x_1-x_1^{k+1},R_{x_1}⟩\ge 0,\forall x_1\in {\mathcal{X}}_1,$$

$$\begin{array}[]{lll}R_{x_{1}}&=&-\tau A_{1}^{\sf T}\lambda^{k+1}-\varepsilon A_{1}^{\sf T}(\lambda^{k+1}-\lambda^{k})+\left(\sigma_{1}+\frac{\varepsilon^{2}-1}{s}\right)A_{1}^{\sf T}A_{1}(x_{1}^{k+1}-x_{1}^{k})\\ \\ &=&\varepsilon A_{1}^{\sf T}\lambda^{k}-(\tau+\varepsilon)A_{1}^{\sf T}\left\{\lambda^{k}-\frac{1}{s}\left[(\tau-\varepsilon)A_{1}x_{1}^{k+1}+\varepsilon A_{1}x_{1}^{k}+\tau\sum\limits_{i=2}^{p}A_{i}x_{i}^{k}-b\tau\right]\right\}\\ \\ &&+\left(\sigma_{1}+\frac{\varepsilon^{2}-1}{s}\right)A_{1}^{\sf T}A_{1}(x_{1}^{k+1}-x_{1}^{k})\\ \\ &=&-\tau A_{1}^{\sf T}\lambda^{k}+\frac{\tau^{2}-\varepsilon^{2}}{s}A_{1}^{\sf T}A_{1}x_{1}^{k+1}+\frac{(\tau+\varepsilon)\tau}{s}A_{1}^{\sf T}\left(\sum\limits_{i=2}^{p}A_{i}x_{i}^{k}-b\right)+\left(\sigma_{1}+\frac{\varepsilon^{2}-1}{s}\right)A_{1}^{\sf T}A_{1}x_{1}^{k+1}\\ \\ &&-\left(\sigma_{1}+\frac{\varepsilon^{2}-1}{s}\right)A_{1}^{\sf T}A_{1}x_{1}^{k}+\frac{(\tau+\varepsilon)\varepsilon}{s}A_{1}^{\sf T}A_{1}x_{1}^{k}\\ \\ &=&-\tau A_{1}^{\sf T}\lambda^{k}+\left(\sigma_{1}+\frac{\tau^{2}-1}{s}\right)A_{1}^{\sf T}A_{1}x_{1}^{k+1}+\frac{(\tau+\varepsilon)\tau}{s}A_{1}^{\sf T}\left(\sum\limits_{i=2}^{p}A_{i}x_{i}^{k}-b\right)\\ \\ &&-\left[\sigma_{1}+\frac{\tau^{2}-1-\tau(\tau+\varepsilon)}{s}\right]A_{1}^{\sf T}A_{1}x_{1}^{k}\\ \\ &=&-\tau A_{1}^{\sf T}\bar{\lambda}^{k}+\bar{\sigma}_{1}A_{1}^{\sf T}A_{1}(x_{1}^{k+1}-x_{1}^{k})\end{array}$$

$${\bar{\lambda }}^k={\lambda }^k-\frac{\tau +\epsilon }{s}\left(\sum _{i=1}^p{A}_i{x}_i^k-b\right).$$

$$x_1^{k+1}$$

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

$$x_i^{k+1}\in {\mathcal{X}}_i,f_i(x_i)-f_i(x_i^{k+1})+⟨x_i-x_i^{k+1},R_{x_i}⟩\ge 0,\forall x_i\in {\mathcal{X}}_i,$$

Full text

General parameterized proximal point algorithm with applications in statistical learning

††thanks: The work was supported by the National Science Foundation of China under grant 11671318 and the Natural Science Foundation of Fujian province under grant 2016J01028.

Jianchao Baia111Corresponding author. E-mail address: bjc1987@163.com.  Jicheng Lia  Pingfan Daia,b  Jiaofen Lic

*a**School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, P.R. China

bDepartment of Information Engineering, Sanming University, Sanming 365004, P.R. China

cCollege of Mathematics and Computational Science, Guilin University of

Electronic Technology, Guilin 541004, P.R. China*

Abstract In the literature, there are a few researches to design some parameters in the Proximal Point Algorithm (PPA), especially for the multi-objective convex optimizations. Introducing some parameters to PPA can make it more flexible and attractive. Mainly motivated by our recent work (Bai et al., A parameterized proximal point algorithm for separable convex optimization, Optim. Lett. (2017) doi: 10.1007/s11590-017-1195-9), in this paper we develop a general parameterized PPA with a relaxation step for solving the multi-block separable structured convex programming. By making use of the variational inequality and some mathematical identities, the global convergence and the worst-case $\mathcal{O}(1/t)$ convergence rate of the proposed algorithm are established. Preliminary numerical experiments on solving a sparse matrix minimization problem from statistical learning validate that our algorithm is more efficient than several state-of-the-art algorithms. Keywords: Structured convex programming; Proximal point algorithm; Relaxation step; Complexity; Statistical learning

AMS subject classifications(2010): 65C60; 65Y20; 90C25

1 Introduction

Throughout this paper, let $\mathbb{R}(\mathbb{R}^{+}),\mathbb{R}^{m},\mathbb{R}^{m\times n}$ be the set of real (positive) numbers, the set of $m$ dimensional real column vectors and the set of $m\times n$ dimensional real matrices, respectively. The symbol $\|z\|_{2}$ denotes the Euclidean norm of $z\in\mathbb{R}^{m}$, which is defined by $\|z\|_{2}=\sqrt{\langle z,z\rangle}$ with the standard inner product $\langle\cdot,\cdot\rangle$. For any symmetric positive definite matrix $G\in\mathbb{R}^{m\times m}$, $\|z\|_{G}=\sqrt{\langle z,Gz\rangle}$ represents the weighted $G$-norm of $z$. We also use T and I to stand for the transpose of a vector/matrix and the identity matrix with proper dimensions, respectively.

Consider the following multi-block separable structured convex optimization

[TABLE]

where $p>1$ is a positive integer, $f_{i}(\cdot):\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}$ are closed convex functions (possibly nonsmooth); $A_{i}\in\mathbb{R}^{m\times n_{i}},b\in\mathbb{R}^{m}$ are respectively given matrices and vector; all structured sets $\mathcal{X}_{i}\subset\mathbb{R}^{n_{i}}(i=1,\cdots,p)$ are closed and convex. Throughout the discussions, we make the following assumption:

Assumption 1.1

The solution set of the problem (1.1) is nonempty and all the matrices $A_{i}(i=1,\cdots,p)$ have full column rank.

Note that the first part of Assumption 1.1 is basic and the second part is necessary. For instance, when dealing with a minimization problem subject to the linearly equality constraint $Ax=b$ involving large-size coefficient matrix and variable, we can split

[TABLE]

to reduce the dimensions of $A$ and $x$, then in such case each $A_{i}$ has full column rank. In fact, many practical application problems do contain at least two different variables (or can be transformed into an equivalent problem with at least two variables), e.g. the total-variational image deblurring problems [11, 20], the transformed joint sparse recovery problem [12], the sparse inverse covariance estimation problem [6], the low-rank and sparse problem [13, 19] and so forth.

The Proximal Point Algorithm (PPA), which was originally proposed to tackle the monotone operator inclusion problems [14, 15], is regarded as a powerful benchmark method for solving the convex problems like (1.1). As verified by Rockafellar [18], the well-known Augmented Lagrangian Method (ALM) for (1.1) with $p=2$ was actually an application of PPA to its dual problem. Moreover, the classical Alternating Direction Method of Multipliers (ADMM) can be also treated as a special variant of PPA to the dual problem [5]. In the last several years, there has been a constantly increasing interest in developing the theories and applications of PPA. For example, He et al.[8] showed a customized application of the classical PPA to the model (1.1) with $p=1$, where some image processing problems were carried out to show the efficiency of the method therein. Ma and Ni [16] proposed a parameterized PPA for (1.1) with $p=1$, where both the basis pursuit problem and the matrix completion problem were tested in experiments to exam the numerical performance of their algorithm. Recently, Cai et al.[3] designed a PPA with a relaxation step for the model (1.1) with $p=2$, whose global convergence and the worst-case sub-linear convergence rate were analyzed in detail. More recently, by introducing some parameters to the metric proximal matrix, an extended parameterized PPA based on [16] was developed for the two block separable convex programming [2], whose effectiveness and robustness was demonstrated by testing a sparse vector optimization problem in the statistical learning compared with two popular algorithms.

To the best of our knowledge, there are a few researches on the parameterized PPA for solving the multi-block model (1.1) with at least three variables. Based on such observation and motivated by our recent published work [2], the aims of this article are to design a general parameterized PPA with a relaxation step (GR-PPA) for solving (1.1) and to test some practical examples having more than two variables to investigate the performance of our GR-PPA. In the remaining parts, Section 2 shows the details of constructing the proposed GR-PPA and analyzing its convergence theories. In Section 3, some numerical experiments are carried out in terms of different parameters, tolerances and initial values. Finally, we conclude the paper in Section 4.

2 Main results

In this section, we first construct a parameterized proximal matrix to design a novel GR-PPA for (1.1), where its convergence is analyzed in detail afterwards. The whole convergence analysis is based on the variational inequality and uses some special techniques to simplify its proof.

2.1 Formation of GR-PPA

For any $\tau\in\mathbb{R}^{+}$, the Lagrangian function of (1.1) is constructed as

[TABLE]

where $\lambda\in\mathbb{R}^{m}$ denotes the Lagrange multiplier with respect to the equality constraint. Let $(x_{1}^{*},\cdots,x_{p}^{*},\lambda^{*})$ be the saddle-point belonging to the solution set $\Omega^{*}$ of (1.1). Then, the following basic inequalities

[TABLE]

imply

[TABLE]

whose optimality conditions are derived as follows

[TABLE]

It is not hard to reformulate the above optimality conditions into a variational inequality:

[TABLE]

where

[TABLE]

[TABLE]

In the following, we also denote

[TABLE]

Clearly, we can get from the skew-symmetric property of the mapping $\mathcal{J}(w)$ that

[TABLE]

It is well-known that the standard PPA with given iterate $w^{k}$ reads the following unified framework:

[TABLE]

where $G$ is a symmetric positive definite matrix (called the proximal matrix). Followed by such framework, a general parameterized proximal matrix will be designed for the new GR-PPA. Concisely, let the matrix $G$ be of following block form

[TABLE]

where $(\sigma_{1},\cdots,\sigma_{p},s)$ are parameters restricted into the domain

[TABLE]

For the sake of analysis convenience, here and next we denote

[TABLE]

Lemma 2.1

Let $\mathcal{K}$ be defined in (2.7) and the matrices $A_{i}(i=1,2,\cdots,p)$ have full column rank. Then, the matrix $G$ in (2.6) is symmetric positive definite for any $(\sigma_{1},\cdots,\sigma_{p},s)\in\mathcal{K}$.

Proof Notice that the matrix $G$ is symmetric and can be decomposed as

[TABLE]

where $D={\rm Diag}(A_{1},\cdots,A_{p},\textbf{I})$ and

[TABLE]

Since all the matrices $A_{i}$ are assumed to have full column rank, therefore, the matrix $G$ is positive definite if and only if $G_{0}$ is positive definite. Using the identity

[TABLE]

it is easy to demonstrate that $G_{0}$ is positive definite if and only if $\widetilde{G}_{0}$ is positive definite, which is guaranteed by the region $\mathcal{K}$ defined in (2.7) because in such case both the upper-left and lower-right block matrices of $\widetilde{G}_{0}$ are positive definite. $\ \ \square$

Remark 2.1

According to the equation (3.7) in [9], the proximal matrix of their proposed algorithm can be decomposed as

[TABLE]

with $P={\rm Diag}(\sqrt{\beta}A_{1},\cdots,\sqrt{\beta}A_{p},\mathbf{I}/\sqrt{\beta})$ and

[TABLE]

in which $\beta>0$ is a penalty parameter with respect to the equality constraint of (1.1) and $\nu\geq m-1$ denotes the proximal parameter. Compared (2.10) to (2.11), if we set

[TABLE]

then it is clear that $\widetilde{G}_{0}=\bar{G}_{0},$ which implies that the proximal matrix of [9] in essence is a special case of (2.6). That is, the method in [9] is a special case of our GR-PPA without relaxation step. A similar way can be used to analyze the proximal matrix in equation (8.6) in [7]. Such relations show that our parameterized proximal matrix is more flexible and general than some in the past.

In what follows, the concrete iterative principles of our GR-PPA are analyzed one by one. Substituting the matrix $G$ into (2.5), we obtain

[TABLE]

which is equivalent to

[TABLE]

Then, it can be derived from the above equality that

[TABLE]

Combining (2.5)-(2.6) and (2.12), we have

[TABLE]

where

[TABLE]

with $\bar{\sigma}_{1}$ defined in (2.8) is positive by (2.7) and

[TABLE]

By (2.13) together with $R_{x_{1}}$, obviously, $x_{1}^{k+1}$ is the solution of the following problem

[TABLE]

Here we get from (2.12) and (2.14) that

[TABLE]

Analogously, for $i=2,3,\cdots,p,$ it follows from (2.5)-(2.6) and (2.12) that

[TABLE]

where

[TABLE]

and

[TABLE]

Therefore, $x_{i}^{k+1}$ is the solution of the following problem

[TABLE]

Consequently, the algorithmic framework of GR-PPA is described in Algorithm 2.1, where we use $\widetilde{w}^{k}=(\widetilde{x}_{1}^{k},\cdots,\widetilde{x}_{p}^{k},\widetilde{\lambda}^{k})$ to replace the original output of the $x_{i}$-subproblems in (2.15), (2.19) and Lagrange multipliers with given iterate $w^{k}=(x_{1}^{k},\cdots,x_{p}^{k},\bar{\lambda}^{k})$, and we use $w^{k+1}=(x_{1}^{k+1},\cdots,x_{p}^{k+1},\bar{\lambda}^{k+1})$ to stand for the new iterate after adding a relaxation step.

Remark 2.2

Noticing that the steps 4-9 in Algorithm 2.1 are actually the PPA updates whose proximal term is implicitly simplified into the quadratic term of each subproblem. Algorithm 2.1 is an extension of our method [2] for the two-block separable convex problem, but the domain of the parameters restricted into (2.7) is not a direct extension of the past. Besides, Algorithm 2.1 adopts a serial idea between $x_{1}$-subproblem and other subproblems, while the parallel idea is used among the $x_{i}$-subproblems ($i=2,\cdots,p$). From the relaxation step of Algorithm 2.1, we immediately have the following certain relationship

[TABLE]

Algorithm 2.1 (GR-PPA for solving Problem (1.1))  1   Choose parameters $(\sigma_{1},\cdots,\sigma_{p},s)\in\mathcal{K},\gamma\in(0,2)$ and initialize $(x_{1}^{0},\cdots,x_{p}^{0},\lambda^{0})\in\Omega$;

2    $r^{0}=\sum\limits_{i=1}^{p}A_{i}x_{i}^{0}-b;\ \bar{\lambda}^{0}=\lambda^{0}-\frac{\tau+\varepsilon}{s}r^{0}$ by (2.14);

3   For $k=0,1,\cdots,$ if not converge, **do

**4   $\widetilde{x}_{1}^{k}=\arg\min\limits_{x_{1}\in\mathcal{X}_{1}}\left\{f_{1}(x_{1})+\frac{\bar{\sigma}_{1}}{2}\left\|A_{1}(x_{1}-x_{1}^{k})-\frac{\tau}{\bar{\sigma}_{1}}\bar{\lambda}^{k}\right\|^{2}\right\}$;

5   $r^{k}=\sum\limits_{i=1}^{p}A_{i}x_{i}^{k}-b;\ \Delta x_{1}^{k}=\widetilde{x}_{1}^{k}-x_{1}^{k}$;

6   $\bar{\lambda}^{k+\frac{1}{2}}=\bar{\lambda}^{k}-\frac{\tau-\varepsilon}{s}\left(2A_{1}\Delta x_{1}^{k}+r^{k}\right)$;

7   Update the $x_{i}$-subproblem for $i=2,\cdots,p$ by (2.19):

[TABLE]

8   $\Delta x_{i}^{k}=\widetilde{x}_{i}^{k}-x_{i}^{k},\ \forall\ i=2,\cdots,p$;

9   $\widetilde{\lambda}^{k}=\bar{\lambda}^{k}-\frac{\tau+\varepsilon}{s}\sum\limits_{i=1}^{p}A_{i}\Delta x_{i}^{k}-\frac{1}{s}\left[(\tau-\varepsilon)A_{1}\Delta x_{1}^{k}+r^{k}\tau\right]$ by (2.16);

10    Relaxation step:

[TABLE]

2.2 Convergence analysis of GR-PPA

This subsection analyzes the global convergence and the ergodic convergence rate of Algorithm 2.1. First of all, we present an important lemma described as follows.

Lemma 2.2

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

[TABLE]

Proof By Algorithm 2.1 and (2.5), we have

[TABLE]

which, by using (2.2) and (2.4) with setting $w=w^{*}$, leads to

[TABLE]

Then, it follows from (2.20) and (2.23) that

[TABLE]

which immediately completes the proof. $\ \ \square$

Lemma 2.2 shows that the sequence $\{w^{k+1}-w^{*}\}$ is strictly contractive under the weighted $G$-norm. Moreover, the following global convergence theorem holds.

Theorem 2.1

Let the parameters $(\sigma_{1},\cdots,\sigma_{p},s)\in\mathcal{K}$ and the sequence $\{w^{k+1}\}$ be generated by Algorithm 2.1. Then, under Assumption 1.1 there exists a point $w^{\infty}\in\Omega^{*}$ such that

[TABLE]

Proof See the proof of Theorem 1 [2]. $\ \ \square$

In order to establish the convergence rate of Algorithm 2.1 in an ergodic sense, we first need to characterize the solution set of $\textrm{VI}(\phi,\mathcal{J},\Omega)$, which had been given by e.g. [11] in the following:

Theorem 2.2

The solution set of $\textrm{VI}(\phi,\mathcal{J},\Omega)$ in (2.2) is convex and can be characterized as

[TABLE]

Theorem 2.3

For any $(\sigma_{1},\cdots,\sigma_{p},s)\in\mathcal{K}$, let

[TABLE]

where $\{\widetilde{w}^{k}\}$ is the iterative sequence of Algorithm 2.1. Then, under Assumption 1.1 we have

[TABLE]

Proof Combining (2.22) with (2.4), we get

[TABLE]

Meanwhile, by (2.20) and the following identity

[TABLE]

with substitutions $(\widetilde{w}^{k},w,w^{k+1},w^{k})=(a,b,c,d)$, it follows that

[TABLE]

where the first two terms

[TABLE]

Combining (2.26)-(2.28), we deduce

[TABLE]

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

[TABLE]

which by the definitions of $\textbf{w}_{t}$ and $\textbf{u}_{t}$ results in

[TABLE]

Because of the convexity of the function $\phi(u)$ (since all $f_{i}$ are assumed to be convex), the following inequality holds

[TABLE]

which, by substituting it into (2.29), completes the whole proof. $\ \ \square$

Remark 2.3

Theorem 2.3 illustrates the worst-case $\mathcal{O}(1/t)$ convergence rate of Algorithm 2.1 in an ergodic sense. By the region $\gamma\in(0,2)$ in Algorithm 2.1 and the inequality (2.25), one may choose a larger value $\gamma$ approximating to 2 so that the value of right-hand of (2.25) is much smaller.

3 Numerical experiments

In this section, we investigate the numerical performance of our proposed GR-PPA for solving a popular sparse matrix optimization problem. All experiments are simulated in MATLAB 7.14 (R2012a) on a lenovo-PC with Intel Core i5 processor (2.70GHz) and 4 GB memory.

3.1 Test problem

The Latent Variable Gaussian Graphical Model Selection (LVGGMS) problem [4, 17] arising in statistical learning is of the following form:

[TABLE]

where $C\in\mathbb{R}^{n\times n}$ is a given covariance matrix obtained from the sample variables, $\nu$ and $\mu$ are given positive weighting factors, $tr(L)$ denotes the trace of the matrix $L$, $\|S\|_{1}=\sum_{i,j=1}^{n}|S_{ij}|$ stands for the $l_{1}$-norm of the matrix $S$ and $S_{ij}$ means its $ij$-th entry.

Clearly, the LVGGMS problem (3.1) can be regarded as a special case of (1.1). And applying Algorithm 2.1 it is easy to write the three corresponding subproblems as the following

[TABLE]

Observe that the above subproblems have closed formula solutions. According to the first-order optimality condition of the $X$-subproblem in (3.2), we derive

[TABLE]

Then, by using the eigenvalue decomposition

[TABLE]

where $\textrm{Diag}(\rho)$ is a diagonal matrix with diagonal entries $\rho_{i}(i=1,\cdots,n)$, we get its explicit solution

[TABLE]

in which $\textrm{Diag}(\gamma)$ is the diagonal matrix with diagonal entries

[TABLE]

Applying the soft shrinkage operator $\textrm{Shrink}(\cdot,\cdot)$, see e.g.[19], the solution of the $S$-subproblem is

[TABLE]

Besides, it is obvious that the $L$-subproblem in (3.2) is equivalent to

[TABLE]

where $\max\{\widetilde{\rho},\textbf{0}\}$ is taken component-wise and $V\textrm{Diag}(\widetilde{\rho})V^{\sf T}$ is the eigenvalue decomposition of the matrix

[TABLE]

3.2 Numerical results

In the following experiments, the parameters $(\nu,\mu)=(0.005,0.05)$ in (3.1) come from [17] and the matrix $C$ is randomly generated by the MATLAB codes of Boyd’s homepage222http://web.stanford.edu/$\sim$boyd/papers/admm/covsel/covsel_example.html. with $n=100$. The maximal number of iterations is set as $1000$ and the following termination criterions are simultaneously used for all algorithms:

[TABLE]

where $(\epsilon_{1},\epsilon_{2},\epsilon_{3})$ are given tolerances, $(X^{k},S^{k},L^{k})$ are the $k$-th iterative values and $F^{*}$ is the approximate optimal objective function value by running Algorithm 2.1 after 1000 iterations. The penalty parameter and the relaxation factor in all involved algorithms are set as $(\beta,\gamma)=(0.05,1.8).$ We choose the feasible points $(X^{0},S^{0},L^{0},\lambda^{0})=(\textbf{I},\textbf{4I},3\textbf{I},\textbf{0})$ as the initial starting iterates.

First of all, for the two-block classical lasso problem, as shown in [2] that Algorithm 2.1 could perform significantly better than ADMM and the customized relaxed PPA [7] when properly choosing the parameters, especially for solving large-size problems and high accurate solutions are required. Followed by the similar way of adjusting the parameters [2], we first would like to investigate the effects of the parameters $(\sigma_{1},\sigma_{2},\sigma_{3},s)$ on Algorithm 2.1 for solving the three-block case (3.1). For this purpose, given $(\varepsilon,\tau)$ we randomly choose some values of one parameter and fix the remaining as known values, then it is clear from (2.7) that

[TABLE]

Here we can note that $\sigma_{1}$ and $\sigma_{i}(i=2,3)$ play the same role if taking $\varepsilon=\tau$, since in such case their regions are the same. Therefore, in the following we would take $\varepsilon=\tau=\frac{\sqrt{5}-1}{2}$ (the golden section ratio) to adjust the relative better value of e.g. $\sigma_{1}$ for an example. We also investigate the effect of the parameter $s$ on Algorithm 2.1.

Under the tolerance $\epsilon_{i}=10^{-8}\ (i=1,2,3)$, experimental results of adjusting the parameters $\sigma_{1}$ and $s$ are respectively reported in Tables 1-2, where the parameter in Table 1 is restricted by $\sigma_{1}>\frac{2+(\sqrt{5}-1)^{2}}{20}\approx 0.1764$ and the parameter in Table 2 is restricted by $s>\frac{2+(\sqrt{5}-1)^{2}}{0.178}\approx 9.9097$. The notations “IT”, “CPU” denote the number of iterations and the CPU time in seconds, respectively. We can observe from Tables 1-2 that:

• •

For the parameters $(\sigma_{1},\sigma_{2},\sigma_{3},s)$, the experimental results in each column of IER and CER are nearly the same when fixing any three parameters with one parameter changing;

• •

With the increase of the parameter $\sigma_{1}$ ($s$), both the iteration number and the CPU time tend to increase, which verifies the comments in the final part of [9]. That is, slow convergence would occur in terms of solving the proximal subproblem if the proximal parameter is set too large, which may significantly affect the whole computational efficiency of the algorithm;

• •

Both Tables 1-2 show that $(\sigma_{1},\sigma_{2},\sigma_{3},s)=\left(0.178,0.178,0.178,10\right)$ would be a relative reasonable choice for tackling the test problem (3.1) and we would use them as the default parameter setting for GR-PPA.

[TABLE]

Table 1: Results of GR-PPA with fixed parameters $(\sigma_{2},\sigma_{3},s)=(0.2,0.2,10)$.

[TABLE]

Table 2: Results of GR-PPA with fixed parameters $(\sigma_{1},\sigma_{2},\sigma_{3})=(0.178,0.178,0.178)$.

Secondly, we would like to use our GR-PPA (Algorithm 2.1) with the tuned parameters to compare with two state-of-the-art algorithms: the Proximal Jacobian decomposition of the ALM [9] (“PJALM”) and the splitting algorithm of the ALM [10] (“HTY”). The default parameter $\mu=2.01$ and $H=\beta\mathbf{I}$ are used for HTY [10], and the proximal parameter of PJALM [9] is fixed as 2 as suggested by the theory therein. Experimental results of these algorithms under different tolerances are reported in Table 3. Furthermore, with fixed tolerances $(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(10^{-12},10^{-14},10^{-12})$, convergence curves of the residuals IER and OER against the number of iterations by different algorithms with different starting iterates are respectively depicted in Figs. 1-2.

[FIGURE:]

It is clear from Table 3 that under higher tolerances our GR-PPA performs significantly better than the others in terms of the number of iterations and the CPU time, although HTY could perform slightly better than GR-PPA under some lower tolerances. Besides, the comparative convergence curves depicted in Figs. 1-2 show that GR-PPA converges faster than the rest two algorithms for different starting feasible points, which illustrates that GR-PPA could performs well for any initial iterative values. Reported results of Table 1 and convergence curves in Figs. 1-2 demonstrate the effectiveness and robustness of the proposed algorithm.

4 Conclusion remarks

In this paper, we develop a relaxed parameterized PPA for solving a class of separable convex minimization problems. The global convergence and a worst-case sub-linear convergence rate of the algorithm are established. Numerical experiments on testing a sparse matrix minimization problem in statistical learning also verify that our proposed algorithm outperforms two popular algorithms when properly choosing the relaxation factor and the parameters in the proximal matrix.

Note that, in Sec.3, we only take a three-block problem of (1.1) for an example to investigate the performance of our algorithm, since its performance for the two-block lasso problem had been verified in [2] compared with the classical ADMM and the relaxed PPA [7]. For the problem with more than three variables, the parameters in GR-PPA need to adjust afresh via experiments. From the framework of PPA in (2.5), one may have other choices to construct a novel parameterized proximal matrix to develop the corresponding PPA only if it is symmetric positive definite. Besides, the proposed algorithm is applicable to the separable convex optimization problem where the coefficient matrices in the linear constraint have full column rank. Hence, whether there exists a PPA for the non-separable case with nonlinear constraints? In addition, we have tried our best to analyze the algorithmic convergence rate in the non-ergodic sense but without any results. These questions need further investigations in the future.

AcknowledgementsThe authors wish to thank the Editor-in-Chief Prof. Choi-Hong Lai and the anonymous referees for providing their valuable suggestions, which have significantly improved the quality of our paper.