Generalized Symmetric ADMM for Separable Convex Optimization

TL;DR

This paper introduces a Generalized Symmetric ADMM that updates dual variables twice with larger stepsizes, improving convergence and flexibility for multi-block convex optimization, especially suitable for big data applications.

Contribution

The paper proposes a novel GS-ADMM algorithm with larger stepsize domain, combining Gauss-Seidel and Jacobi updates, and provides convergence guarantees and practical effectiveness.

Findings

Global convergence with $O(1/t)$ rate established.

Larger stepsize domain compared to existing methods.

Effective in sparse matrix minimization for statistical learning.

Abstract

The Alternating Direction Method of Multipliers (ADMM) has been proved to be effective for solving separable convex optimization subject to linear constraints. In this paper, we propose a Generalized Symmetric ADMM (GS-ADMM), which updates the Lagrange multiplier twice with suitable stepsizes, to solve the multi-block separable convex programming. This GS-ADMM partitions the data into two group variables so that one group consists of $p$ block variables while the other has $q$ block variables, where $p \ge 1$ and $q \ge 1$ are two integers. The two grouped variables are updated in a {\it Gauss-Seidel} scheme, while the variables within each group are updated in a {\it Jacobi} scheme, which would make it very attractive for a big data setting. By adding proper proximal terms to the subproblems, we specify the domain of the stepsizes to guarantee that GS-ADMM is globally convergent with a worst-case $O(1/t)$ ergodic convergence rate. It turns out that our convergence domain of the stepsizes is significantly larger than other convergence domains in the literature. Hence, the GS-ADMM is more flexible and attractive on choosing and using larger stepsizes of the dual variable. Besides, two special cases of GS-ADMM, which allows using zero penalty terms, are also discussed and analyzed. Compared with several state-of-the-art methods, preliminary numerical experiments on solving a sparse matrix minimization problem in the statistical learning show that our proposed method is effective and promising.

Equations

$$\begin{array}[]{lll}\min&\sum\limits_{i=1}^{p}f_{i}(x_{i})+\sum\limits_{j=1}^{q}g_{j}(y_{j})\\ \textrm{s.t. }&\sum\limits_{i=1}^{p}A_{i}x_{i}+\sum\limits_{j=1}^{q}B_{j}y_{j}=c,\\ &x_{i}\in\mathcal{X}_{i},\;i=1,\cdots,p,\\ &y_{j}\in\mathcal{Y}_{j},\;j=1,\cdots,q,\\ \end{array}$$

$${\mathcal{L}}_{\beta }\left(\mathbf{x},\mathbf{y},\lambda \right)=L\left(\mathbf{x},\mathbf{y},\lambda \right)+\frac{\beta }{2}\Vert \mathcal{A}\mathbf{x}+\mathcal{B}\mathbf{y}-c{\Vert }^2,$$

$$L\left(\mathbf{x},\mathbf{y},\lambda \right)=\sum _{i=1}^p{f}_i(x_i)+\sum _{j=1}^q{g}_j(y_j)-⟨\lambda ,\mathcal{A}\mathbf{x}+\mathcal{B}\mathbf{y}-c⟩$$

$$\left\{\begin{array}[]{l}\left({\bf{x}}^{k+1},{\bf{y}}^{k+1}\right)=\arg\min\left\{\mathcal{L}_{\beta}\left({\bf{x}},{\bf{y}},\lambda^{k}\right)|\ {\bf{x}}\in\mathcal{X},{\bf{y}}\in\mathcal{Y}\right\},\\ \lambda^{k+1}=\lambda^{k}-\beta({\cal{A}}{\bf{x}}^{k+1}+{\cal{B}}{\bf{y}}^{k+1}-c).\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\textrm{For}\ i=1,2,\cdots,p,\\ \qquad x_{i}^{k+1}=\arg\min\limits_{x_{i}\in\mathcal{X}_{i}}\mathcal{L}_{\beta}\left(x_{1}^{k+1},\cdots,x_{i-1}^{k+1},x_{i},x_{i+1}^{k},\cdots,x_{p}^{k},{\bf{y}}^{k},\lambda^{k}\right),\\ \textrm{For}\ j=1,2,\cdots,q,\\ \qquad y_{j}^{k+1}=\arg\min\limits_{y_{j}\in\mathcal{Y}_{j}}\mathcal{L}_{\beta}\left({\bf{x}}^{k+1},y_{1}^{k+1},\cdots,y_{j-1}^{k+1},y_{j},y_{j+1}^{k},\cdots,y_{q}^{k},\lambda^{k}\right),\\ \lambda^{k+1}=\lambda^{k}-\beta\left(\mathcal{A}{\bf{x}}^{k+1}+\mathcal{B}{\bf{y}}^{k+1}-c\right).\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\textrm{For}\ i=1,2,\cdots,p,\\ \qquad x_{i}^{k+1}=\arg\min\limits_{x_{i}\in\mathcal{X}_{i}}\mathcal{L}_{\beta}\left(x_{1}^{k},\cdots,x_{i-1}^{k},x_{i},x_{i+1}^{k},\cdots,x_{p}^{k},{\bf{y}}^{k},\lambda^{k}\right),\\ \textrm{For}\ j=1,2,\cdots,q,\\ \qquad y_{j}^{k+1}=\arg\min\limits_{y_{j}\in\mathcal{Y}_{j}}\mathcal{L}_{\beta}\left({\bf{x}}^{k},y_{1}^{k},\cdots,y_{j-1}^{k},y_{j},y_{j+1}^{k},\cdots,y_{q}^{k},\lambda^{k}\right),\\ \lambda^{k+1}=\lambda^{k}-\beta\left(\mathcal{A}{\bf{x}}^{k+1}+\mathcal{B}{\bf{y}}^{k+1}-c\right).\end{array}\right.$$

$$\left\{\begin{array}[]{lll}{\bf{x}}^{k+1}=\arg\min\left\{\mathcal{L}_{\beta}\left({\bf{x}},{\bf{y}}^{k},\lambda^{k}\right)\ |\ {\bf{x}}\in\mathcal{X}\right\},\\ \lambda^{k+\frac{1}{2}}=\lambda^{k}-\tau\beta\left({\cal{A}}{\bf{x}}^{k+1}+{\cal{B}}{\bf{y}}^{k}-c\right),\\ {\bf{y}}^{k+1}=\arg\min\left\{\mathcal{L}_{\beta}({\bf{x}}^{k+1},{\bf{y}},\lambda^{k+\frac{1}{2}})\ |\ {\bf{y}}\in\mathcal{Y}\right\},\\ \lambda^{k+1}=\lambda^{k+\frac{1}{2}}-s\beta\left({\cal{A}}{\bf{x}}^{k+1}+{\cal{B}}{\bf{y}}^{k+1}-c\right),\end{array}\right.$$

$$\mathcal{H}=\left\{(\tau ,s)|s\in (0,(\sqrt{5}+1)/2),\tau +s>0,\tau \in (-1,1),|\tau |<1+s-s^2\right\}$$

$$(\textbf{GS-ADMM})\ \left\{\begin{array}[]{lll}\textrm{For}\ i=1,2,\cdots,p,\\ \quad x_{i}^{k+1}=\arg\min\limits_{x_{i}\in\mathcal{X}_{i}}\mathcal{L}_{\beta}(x_{1}^{k},\cdots,x_{i},\cdots,x_{p}^{k},{\bf{y}}^{k},\lambda^{k})+P_{i}^{k}(x_{i}),\\ \quad\textrm{where }P_{i}^{k}(x_{i})=\frac{\sigma_{1}\beta}{2}\left\|A_{i}(x_{i}-x_{i}^{k})\right\|^{2},\\ \lambda^{k+\frac{1}{2}}=\lambda^{k}-\tau\beta(\mathcal{A}{\bf{x}}^{k+1}+\mathcal{B}{\bf{y}}^{k}-c),\\ \\ \textrm{For}\ j=1,2,\cdots,q,\\ \quad y_{j}^{k+1}=\arg\min\limits_{y_{j}\in\mathcal{Y}_{j}}\mathcal{L}_{\beta}({\bf{x}}^{k+1},y_{1}^{k},\cdots,y_{j},\cdots,y_{q}^{k},\lambda^{k+\frac{1}{2}})+Q_{j}^{k}(y_{j}),\\ \quad\textrm{where }Q_{j}^{k}(y_{j})=\frac{\sigma_{2}\beta}{2}\left\|B_{j}(y_{j}-y_{j}^{k})\right\|^{2},\\ \lambda^{k+1}=\lambda^{k+\frac{1}{2}}-s\beta(\mathcal{A}{\bf{x}}^{k+1}+\mathcal{B}{\bf{y}}^{k+1}-c).\end{array}\right.$$

$$(\tau ,s)\in \mathcal{K}=\left\{(\tau ,s)|\tau +s>0,\tau \le 1,-{\tau }^2-s^2-\tau s+\tau +s+1>0\right\},$$

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

$$L\left({\mathbf{x}}^{\ast },{\mathbf{y}}^{\ast },\lambda \right)\le L\left({\mathbf{x}}^{\ast },{\mathbf{y}}^{\ast },{\lambda }^{\ast }\right)\le L\left(\mathbf{x},\mathbf{y},{\lambda }^{\ast }\right),\forall (\mathbf{x},\mathbf{y},\lambda )\in \mathcal{M}$$

$$\left\{\begin{array}[]{lll}\textrm{For}\ i=1,2,\cdots,p,\\ \qquad x_{i}^{*}=\arg\min\limits_{x_{i}\in\mathcal{X}_{i}}\mathcal{L}_{\beta}\left(x_{1}^{*},\cdots,x_{i-1}^{*},x_{i},x_{i+1}^{*},\cdots,x_{p}^{*},{\bf{y}}^{*},\lambda^{*}\right),\\ \textrm{For}\ j=1,2,\cdots,q,\\ \qquad y_{j}^{*}=\arg\min\limits_{y_{j}\in\mathcal{Y}_{j}}\mathcal{L}_{\beta}\left({\bf{x}}^{*},y_{1}^{*},\cdots,y_{j-1}^{*},y_{j},y_{j+1}^{*},\cdots,y_{q}^{*},\lambda^{*}\right),\\ \lambda^{*}=\arg\max\limits_{\lambda\in\mathcal{R}^{n}}\mathcal{L}_{\beta}\left({\bf{x}}^{*},{\bf{y}}^{*},\lambda\right).\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\textrm{For}\ i=1,2,\cdots,p,\\ \qquad x_{i}^{*}\in\mathcal{X}_{i},\quad f_{i}(x_{i})-f_{i}(x_{i}^{*})+\langle x_{i}-x_{i}^{*},-A_{i}^{\sf T}\lambda^{*}\rangle\geq 0,\forall x_{i}\in\mathcal{X}_{i},\\ \textrm{For}\ j=1,2,\cdots,q,\\ \qquad y_{j}^{*}\in\mathcal{Y}_{j},\quad g_{j}(y_{j})-g_{j}(y_{j}^{*})+\langle y_{j}-y_{j}^{*},-B_{j}^{\sf T}\lambda^{*}\rangle\geq 0,\forall y_{j}\in\mathcal{Y}_{j},\\ \lambda^{*}\in\mathcal{R}^{n},\quad\left\langle\lambda-\lambda^{*},{\cal{A}}{\bf{x}}^{*}+{\cal{B}}{\bf{y}}^{*}-c\right\rangle\geq 0,\forall\lambda\in\mathcal{R}^{n}.\end{array}\right.$$

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

$$h(\mathbf{u})=\sum _{i=1}^p{f}_i(x_i)+\sum _{j=1}^p{g}_j(y_j)$$

$${\bf{u}}=\left(\begin{array}[]{c}{\bf{x}}\\ {\bf{y}}\\ \end{array}\right),\quad{\bf{w}}=\left(\begin{array}[]{c}{\bf{x}}\\ {\bf{y}}\\ \lambda\end{array}\right),\quad\mathcal{J}({\bf{w}})=\left(\begin{array}[]{c}-{\cal{A}}^{\sf T}\lambda\\ -{\cal{B}}^{\sf T}\lambda\\ {\cal{A}}{\bf{x}}+{\cal{B}}{\bf{y}}-c\end{array}\right).$$

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

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

$${\mathcal{M}}^{\ast }=\bigcap _{\mathbf{w}\in \mathcal{M}}\left\{\hat{\mathbf{w}}\in \mathcal{M}|h(\mathbf{u})-h(\hat{\mathbf{u}})+⟨\mathbf{w}-\hat{\mathbf{w}},\mathcal{J}(\mathbf{w})⟩\ge 0\right\}.$$

$$\widetilde{{\bf{x}}}^{k}=\left(\begin{array}[]{c}\widetilde{x}_{1}^{k}\\ \widetilde{x}_{2}^{k}\\ \vdots\\ \widetilde{x}_{p}^{k}\end{array}\right)=\left(\begin{array}[]{c}x_{1}^{k+1}\\ x_{2}^{k+1}\\ \vdots\\ x_{p}^{k+1}\end{array}\right),\quad\widetilde{{\bf{y}}}^{k}=\left(\begin{array}[]{c}\widetilde{y}_{1}^{k}\\ \widetilde{y}_{2}^{k}\\ \vdots\\ \widetilde{y}_{q}^{k}\end{array}\right)=\left(\begin{array}[]{c}y_{1}^{k+1}\\ y_{2}^{k+1}\\ \vdots\\ y_{q}^{k+1}\end{array}\right),$$

$${\tilde{\lambda }}^k={\lambda }^k-\beta \left(\mathcal{A}{\mathbf{x}}^{k+1}+\mathcal{B}{\mathbf{y}}^k-c\right),$$

$$\widetilde{{\bf{u}}}=\left(\begin{array}[]{c}\widetilde{{\bf{x}}}\\ \widetilde{{\bf{y}}}\\ \end{array}\right),\quad\widetilde{{\bf{w}}}^{k}=\left(\begin{array}[]{c}\widetilde{{\bf{x}}}^{k}\\ \widetilde{{\bf{y}}}^{k}\\ \widetilde{\lambda}^{k}\end{array}\right)=\left(\begin{array}[]{c}{\bf{x}}^{k+1}\\ {\bf{y}}^{k+1}\\ \widetilde{\lambda}^{k}\end{array}\right).$$

$$h(\mathbf{u})-h({\tilde{\mathbf{u}}}^k)+⟨\mathbf{w}-{\tilde{\mathbf{w}}}^k,\mathcal{J}({\tilde{\mathbf{w}}}^k)+Q({\tilde{\mathbf{w}}}^k-{\mathbf{w}}^k)⟩\ge 0,\forall \mathbf{w}\in \mathcal{M},$$

$$Q=\left[\begin{array}[]{cc}H_{\mathbf{x}}&\bf{0}\\ \bf{0}&\widetilde{Q}\end{array}\right]$$

$$H_{\mathbf{x}}=\beta\left[\begin{array}[]{ccccccccccccc}\sigma_{1}A_{1}^{\sf T}A_{1}&&&&-A_{1}^{\sf T}A_{2}&&&&\cdots&&&&-A_{1}^{\sf T}A_{p}\\ \\ -A_{2}^{\sf T}A_{1}&&&&\sigma_{1}A_{2}^{\sf T}A_{2}&&&&\cdots&&&&-A_{2}^{\sf T}A_{p}\\ \\ \vdots&&&&\vdots&&&&\ddots&&&&\vdots\\ \\ -A_{p}^{\sf T}A_{1}&&&&-A_{p}^{\sf T}A_{2}&&&&\cdots&&&&\sigma_{1}A_{p}^{\sf T}A_{p}\end{array}\right],$$

$$\widetilde{Q}=\left[\begin{array}[]{ccccccccccccccc|cc}(\sigma_{2}+1)\beta B_{1}^{\sf T}B_{1}&&&&{\bf 0}&&&&\cdots&&&&{\bf 0}&&&&-\tau B_{1}^{\sf T}\\ &&&&&&&&&&&&&&&&\\ {\bf 0}&&&&(\sigma_{2}+1)\beta B_{2}^{\sf T}B_{2}&&&&\cdots&&&&{\bf 0}&&&&-\tau B_{2}^{\sf T}\\ &&&&&&&&&&&&&&&&\\ \vdots&&&&\vdots&&&&\ddots&&&&\vdots&&&&\vdots\\ &&&&&&&&&&&&&&&&\\ {\bf 0}&&&&{\bf 0}&&&&\cdots&&&&(\sigma_{2}+1)\beta B_{q}^{\sf T}B_{q}&&&&-\tau B_{q}^{\sf T}\\ &&&&&&&&&&&&&&&&\\ \hline\cr-B_{1}&&&&-B_{2}&&&&\cdots&&&&-B_{q}&&&&\frac{1}{\beta}I\end{array}\right].$$

$$x_i^{k+1}=\arg \underset{x_i\in {\mathcal{X}}_i}{\min }\left\{f_i(x_i)-⟨{\lambda }^k,A_i{x}_i⟩+\frac{\beta }{2}{\Vert A_i{x}_i-c_{x,i}\Vert }^2+\frac{{\sigma }_1\beta }{2}{\Vert A_i(x_i-x_i^k)\Vert }^2\right\},$$

$$f_i(x_i)-f_i(x_i^{k+1})+⟨A_i(x_i-x_i^{k+1}),-{\lambda }^k+\beta (A_i{x}_i^{k+1}-c_{x,i})+{\sigma }_1\beta A_i(x_i^{k+1}-x_i^k)⟩\ge 0$$

$$f_i(x_i)-f_i({\tilde{x}}_i^k)+⟨A_i(x_i-{\tilde{x}}_i^k),-{\tilde{\lambda }}^k-\beta \sum _{l=1,l\ne i}^p{A}_l({\tilde{x}}_l^k-x_l^k)+{\sigma }_1\beta \sum _{l=1}^p{A}_l({\tilde{x}}_l^k-x_l^k)⟩\ge 0$$

Full text

1 Introduction

We consider the following grouped multi-block separable convex programming problem

[TABLE]

where $f_{i}(x_{i}):\mathcal{R}^{m_{i}}\rightarrow\mathcal{R},\ g_{j}(y_{j}):\mathcal{R}^{d_{j}}\rightarrow\mathcal{R}$ are closed and proper convex functions (possibly nonsmooth); $A_{i}\in\mathcal{R}^{n\times m_{i}},B_{j}\in\mathcal{R}^{n\times d_{j}}$ and $c\in\mathcal{R}^{n}$ are given matrices and vectors, respectively; $\mathcal{X}_{i}\subset\mathcal{R}^{m_{i}}$ and $\mathcal{Y}_{j}\subset\mathcal{R}^{d_{j}}$ are closed convex sets; $p\geq 1$ and $q\geq 1$ are two integers. Throughout this paper, we assume that the solution set of the problem (1) is nonempty and all the matrices $A_{i}$, $i=1,\cdots,p$, and $B_{j}$, $j=1,\cdots,q$, have full column rank. And in the following, we denote $\mathcal{A}=\left(A_{1},\cdots,A_{p}\right),\mathcal{B}=\left(B_{1},\cdots,B_{q}\right)$, ${\bf{x}}=(x_{1}^{\sf T},\cdots,x_{p}^{\sf T})^{\sf T}$, ${\bf{y}}=(y_{1}^{\sf T},\cdots,y_{q}^{\sf T})^{\sf T}$, ${\cal{X}}={\cal{X}}_{1}\times{\cal{X}}_{2}\times\cdots{\cal{X}}_{p}$, ${\cal{Y}}={\cal{Y}}_{1}\times{\cal{Y}}_{2}\times\cdots{\cal{Y}}_{q}$ and $\mathcal{M}=\mathcal{X}\times\mathcal{Y}\times\mathcal{R}^{n}$.

In the last few years, the problem (1) has been extensively investigated due to its wide applications in different fields, such as the sparse inverse covariance estimation problem RothmanBickel2008 in finance and statistics, the model updating problem DongYuTian2015 in the design of vibration structural dynamic system and bridges, the low rank and sparse representations LiuLiBai2017 in image processing and so forth. One standard way to solve the problem (1) is the classical Augmented Lagrangian Method (ALM) Hestenes1969 , which minimizes the following augmented Lagrangian function

[TABLE]

where $\beta>0$ is a penalty parameter for the equality constraint and

[TABLE]

is the Lagrangian function of the problem (1) with the Lagrange multiplier $\lambda\in\mathcal{R}^{n}$. Then, the ALM procedure for solving (1) can be described as follows:

[TABLE]

However, ALM does not make full use of the separable structure of the objective function of (1) and hence, could not take advantage of the special properties of the component objective functions $f_{i}$ and $g_{j}$ in (1). As a result, in many recent real applications involving big data, solving the subproblems of ALM becomes very expensive.

One effective approach to overcome such difficulty is the Alternating Direction Method of Multipliers (ADMM), which was originally proposed in GlowinskiMarrocco1975 and could be regarded as a splitting version of ALM. At each iteration, ADMM first sequentially optimize over one block variable while fixing all the other block variables, and then follows by updating the Lagrange multiplier. A natural extension of ADMM for solving the multi-block case problem (1) takes the following iterations:

[TABLE]

Obviously, the scheme (3) is a serial algorithm which uses the newest information of the variables at each iteration. Although the above scheme was proved to be convergent for the two-block, i.e., $p=q=1$, separable convex minimization (see HeYuan2012 ), as shown in Chen2016 , the direct extension of ADMM (3) for the multi-block case, i.e., $p+q\geq 3$, without proper modifications is not necessarily convergent. Another natural extension of ADMM is to use the Jacobian fashion, where the variables are updated simultaneously after each iteration, that is,

[TABLE]

As shown in HeHouYuan2015 , however, the Jacobian scheme (4) is not necessarily convergent either. To ensure the convergence, He et al. HeTaoYuan2015 proposed a novel ADMM-type splitting method that by adding certain proximal terms, allowed some of the subproblems to be solved in parallel, i.e., in a Jacobian fashion. And in HeTaoYuan2015 , some sparse low-rank models and image painting problems were tested to verify the efficiency of their method.

More recently, a Symmetric ADMM (S-ADMM) was proposed by He et al. HeMaYuan2016 for solving the two-block separable convex minimization, where the algorithm performs the following updating scheme:

[TABLE]

and the stepsizes $(\tau,s)$ were restricted into the domain

[TABLE]

in order to ensure its global convergence. The main improvement of HeMaYuan2016 is that the scheme (5) largely extends the domain of the stepsizes $(\tau,s)$ of other ADMM-type methods HeLiuWangYuan2014 . What’s more, the numerical performance of S-ADMM on solving the widely used basis pursuit model and the total-variational image debarring model significantly outperforms the original ADMM in both the CPU time and the number of iterations. Besides, Gu, et al.GuJiang2015 also studied a semi-proximal-based strictly contractive Peaceman-Rachford splitting method, that is (5) with two additional proximal penalty terms for the $x$ and $y$ update. But their method has a nonsymmetric convergence domain of the stepsize and still focuses on the two-block case problem, which limits its applications for solving large-scale problems with multiple block variables.

Mainly motivated by the work of HeTaoYuan2015 ; HeMaYuan2016 ; GuJiang2015 , we would like to generalize S-ADMM with more wider convergence domain of the stepsizes to tackle the multi-block separable convex programming model (1), which more frequently appears in recent applications involving big data Chandrasekaran2012 ; Ma2017 . Our algorithm framework can be described as follows:

[TABLE]

In the above Generalized Symmetric ADMM (GS-ADMM), $\tau$ and $s$ are two stepsize parameters satisfying

[TABLE]

and $\sigma_{1}\in(p-1,+\infty),\sigma_{2}\in(q-1,+\infty)$ are two proximal parameters111 Note that these two parameters are strictly positive in (7). In Section 4, however, we analyze two special cases of GS-ADMM allowing either $\sigma_{1}$ or $\sigma_{2}$ to be zero. for the regularization terms $P_{i}^{k}(\cdot)$ and $Q_{j}^{k}(\cdot)$. He and YuanHeyuan2015 also investigated the above GS-ADMM (7) but restricted the stepsize $\tau=s\in(0,1)$, which does not exploit the advantages of using flexible stepsizes given in (8) to improve its convergence.

Major contributions of this paper can be summarized as the following four aspects:

• •

Firstly, the new GS-ADMM could deal with the multi-block separable convex programming problem (1), while the original S-ADMM in HeMaYuan2016 only works for the two block case and may not be convenient for solving large-scale problems. In addition, the convergence domain $\mathcal{K}$ for the stepsizes $(\tau,s)$ in (8), shown in Fig. 1, is significantly larger than the domain $\mathcal{H}$ given in (6) and the convergence domain in GuJiang2015 ; Heyuan2015 . For example, the stepsize $s$ can be arbitrarily close to $5/3$ when the stepsize $\tau$ is close to $-1/3$. Moreover, the above domain in (8) is later enlarged to a symmetric domain ${\cal{G}}$ defined in (112), shown in Fig. 2. Numerical experiments in Sec. 5.2.1 also validate that using more flexible and relatively larger stepsizes $(\tau,s)$ can often improve the convergence speed of GS-ADMM. On the other hand, we can see that when $\tau=0$, the stepsize $s$ can be chosen in the interval $(0,(\sqrt{5}+1)/2)$, which was firstly suggested by Fortin and Glowinski in FortinGlowinski1983 ; Glowinski1984 .

• •

Secondly, the global convergence of GS-ADMM as well as its worst-case ${\cal{O}}(1/t)$ ergodic convergence rate are established. What’s more, the total $p+q$ block variables are partitioned into two grouped variables. While a Gauss-Seidel fashion is taken between the two grouped variables, the block variables within each group are updated in a Jacobi scheme. Hence, parallel computing can be implemented for updating the variables within each group, which could be critical in some scenarios for problems involving big data.

• •

Thirdly, we discuss two special cases of GS-ADMM, which is (7) with $p\geq 1,q=1$ and $\sigma_{2}=0$ or with $p=1,q\geq 1$ and $\sigma_{1}=0$. These two special cases of GS-ADMM were not discussed in Heyuan2015 and in fact, to the best of our knowledge, they have not been studied in the literature. We show the convergence domain of the stepsizes $(\tau,s)$ for these two cases is still ${\cal{K}}$ defined in (8) that is larger than $\mathcal{H}$.

• •

Finally, numerical experiments are performed on solving a sparse matrix optimization problem arising from the statistical learning. We have investigated the effects of the stepsizes $(\tau,s)$ and the penal parameter $\beta$ on the performance of GS-ADMM. And our numerical experiments demonstrate that by properly choosing the parameters, GS-ADMM could perform significantly better than other recently quite popular methods developed in BaiLiLi2017 ; HeTaoYuan2015 ; HeXuYuan2016 ; WangSong2017 .

The paper is organized as follows. In Section 2, some preliminaries are given to reformulate the problem (1) into a variational inequality and to interpret the GS-ADMM (7) as a prediction-correction procedure. Section 3 investigates some properties of $\|{\bf{w}}^{k}-{\bf{w}}^{*}\|_{H}^{2}$ and provides a lower bound of $\|{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}\|_{G}^{2}$, where $H$ and $G$ are some particular symmetric matrices. Then, we establish the global convergence of GS-ADMM and show its convergence rate in an ergodic sense. In Section 4, we discuss two special cases of GS-ADMM, in which either the penalty parameters $\sigma_{1}$ or $\sigma_{2}$ is allowed to be zero. Some preliminary numerical experiments are done in Section 5. We finally make some conclusions in Section 6.

1.1 Notation

Denoted by $\mathcal{R},\mathcal{R}^{n},\mathcal{R}^{m\times n}$ be the set of real numbers, the set of $n$ dimensional real column vectors and the set of $m\times n$ real matrices, respectively. For any $x,y\in\mathcal{R}^{n}$, $\langle x,y\rangle=x^{\sf T}y$ denotes their inner product and $\|x\|=\sqrt{\langle x,x\rangle}$ denotes the Euclidean norm of $x$, where the superscript T is the transpose. Given a symmetric matrix $G$, we define $\|x\|_{G}^{2}=x^{\sf T}Gx$. Note that with this convention, $\|x\|_{G}^{2}$ is not necessarily nonnegative unless $G$ is a positive definite matrix ($\succeq\textbf{0}$). For convenience, we use $I$ and $\bf{0}$ to stand respectively for the identity matrix and the zero matrix with proper dimension throughout the context.

2 Preliminaries

In this section, we first use a variational inequality to characterize the solution set of the problem (1). Then, we analyze that GS-ADMM (7) can be treated as a prediction-correction procedure involving a prediction step and a correction step.

2.1 Variational reformulation of (1)

We begin with the following standard lemma whose proof can be found in HeMaYuan2016 ; He2016 .

Lemma 1

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

[TABLE]

It is well-known in optimization that a triple point $({\bf{x}}^{*},{\bf{y}}^{*},\lambda^{*})\in{\cal{M}}$ is called the saddle-point of the Lagrangian function (2) if it satisfies

[TABLE]

which can be also characterized as

[TABLE]

Then, by Lemma 1, the above saddle-point equations can be equivalently expressed as

[TABLE]

Rewriting (9) in a more compact variational inequality (VI) form, we have

[TABLE]

where

[TABLE]

and

[TABLE]

Noticing that the affine mapping ${\cal{J}}$ is skew-symmetric, we immediately get

[TABLE]

Hence, (10) can be also rewritten as

[TABLE]

Because of the nonempty assumption on the solution set of (1), the solution set $\mathcal{M}^{*}$ of the variational inequality $\textrm{VI}(h,\mathcal{J},\mathcal{M})$ is also nonempty and convex, see e.g. Theorem 2.3.5 FacchineiPang2003 for more details. The following theorem given by Theorem 2.1 HeYuan2012 provides a concrete way to characterize the set $\mathcal{M}^{*}$.

Theorem 2.1

The solution set of the variational inequality $\textrm{VI}(h,\mathcal{J},\mathcal{M})$ is convex and can be expressed as

[TABLE]

2.2 A prediction-correction interpretation of GS-ADMM

Following a similar approach in HeMaYuan2016 , we next interpret GS-ADMM as a prediction-correction procedure. First, let

[TABLE]

[TABLE]

and

[TABLE]

Then, by using the above notations, we derive the following basic lemma.

Lemma 2

For the iterates $\widetilde{{\bf{u}}}^{k},\widetilde{{\bf{w}}}^{k}$ defined in (15), we have $\widetilde{{\bf{w}}}^{k}\in\mathcal{M}$ and

[TABLE]

where

[TABLE]

with

[TABLE]

[TABLE]

Proof

Omitting some constants, it is easy to verify that the $x_{i}$-subproblem $(i=1,2,\cdots,p)$ of GS-ADMM can be written as

[TABLE]

where $c_{x,i}=c-\sum\limits_{l=1,l\neq i}^{p}A_{l}x_{l}^{k}-\mathcal{B}{\bf{y}}^{k}$. Hence, by Lemma 1, we have $x_{i}^{k+1}\in\mathcal{X}_{i}$ and

[TABLE]

for any $x_{i}\in\mathcal{X}_{i}$. So, by the definition of (13) and (14), we get

[TABLE]

for any $x_{i}\in\mathcal{X}_{i}$. By the way of generating $\lambda^{k+\frac{1}{2}}$ in (7) and the definition of (14), the following relation holds

[TABLE]

Similarly, the $y_{j}$-subproblem ($j=1,\cdots,q$) of GS-ADMM can be written as

[TABLE]

where $c_{y,j}=c-\mathcal{A}{\bf{x}}^{k+1}-\sum\limits_{l=1,l\neq j}^{q}B_{l}y_{l}^{k}$. Hence, by Lemma 1, we have $y_{j}^{k+1}\in\mathcal{Y}_{j}$ and

[TABLE]

for any $y_{j}\in{\cal{Y}}_{j}$. This inequality, by using (21) and the definition of (13) and (14), can be rewritten as

[TABLE]

for any $y_{j}\in{\cal{Y}}_{j}$. Besides, the equality (14) can be rewritten as

[TABLE]

which is equivalent to

[TABLE]

Then, (16) follows from (20), (22) and (23). $\ \ \ \diamondsuit$

Lemma 3

For the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ generated by GS-ADMM, the following equality holds

[TABLE]

where

[TABLE]

Proof

It follows from the way of generating $\lambda^{k+1}$ in the algorithm and (21) that

[TABLE]

The above equality together with $x_{i}^{k+1}=\widetilde{x}_{l}^{k}$, for $i=1,\cdots,p$, and $y_{j}^{k+1}=\widetilde{y}_{j}^{k}$, for $j=1,\cdots,q$, imply

[TABLE]

which immediately gives (24). $\ \ \ \diamondsuit$

Lemma 2 and Lemma 3 show that our GS-ADMM can be interpreted as a prediction-correction framework, where ${\bf{w}}^{k+1}$ and $\widetilde{{\bf{w}}}^{k}$ are normally called the predictive variable and the correcting variable, respectively.

3 Convergence analysis of GS-ADMM

Compared with (12) and (16), the key to proving the convergence of GS-ADMM is to verify that the extra term in (16) converges to zero, that is,

[TABLE]

In this section, we first investigate some properties of the sequence $\{\|{\bf{w}}^{k}-{\bf{w}}^{*}\|_{H}^{2}\}$. Then, we provide a lower bound of $\|{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}\|_{G}^{2}$. Based on these properties, the global convergence and worst-case $\mathcal{O}(1/t)$ convergence rate of GS-ADMM are established in the end.

3.1 Properties of $\{\left\|{\bf{w}}^{k}-{\bf{w}}^{*}\right\|_{H}^{2}\}$

It follows from (11) and (16) that $\widetilde{{\bf{w}}}^{k}\in\mathcal{M}$ and

[TABLE]

Suppose $\tau+s>0$. Then, the matrix $M$ defined in (25) is nonsingular. So, by (24) and a direct calculation, the right-hand term of (26) is rewritten as

[TABLE]

where

[TABLE]

with $H_{x}$ defined in (18) and

[TABLE]

The following lemma shows that $H$ is a positive definite matrix for proper choice of the parameters $(\sigma_{1},\sigma_{2},\tau,s)$.

Lemma 4

The matrix $H$ defined in (28) is symmetric positive definite if

[TABLE]

Proof

By the block structure of $H$, we only need to show that the blocks $H_{\mathbf{x}}$ and $\widetilde{H}$ in (28) are positive definite if the parameters $(\sigma_{1},\sigma_{2},\tau,s)$ satisfy (30). Note that

[TABLE]

where

[TABLE]

If $\sigma_{1}\in(p-1,+\infty)$, $H_{\mathbf{x},0}$ is positive definite. Then, it follows from (31) that $H_{\mathbf{x}}$ is positive definite if $\sigma_{1}\in(p-1,+\infty)$ and all $A_{i}$, $i=1,\cdots,p$, have full column rank.

Now, note that the matrix $\widetilde{H}$ can be decomposed as

[TABLE]

where

[TABLE]

and

[TABLE]

According to the fact that

[TABLE]

where

[TABLE]

we have $\widetilde{H}_{0}$ is positive definite if and only if $H_{\mathbf{y},0}+(1-\tau)EE^{\sf T}$ is positive definite and $\tau+s>0$. Note that $H_{\mathbf{y},0}$ is positive definite if $\sigma_{2}\in(q-1,+\infty)$, and $(1-\tau)EE^{\sf T}$ is positive semidefinite if $\tau\leq 1$. So, $\widetilde{H}_{0}$ is positive definite if $\sigma_{2}\in(q-1,+\infty)$, $\tau\leq 1$ and $\tau+s>0$. Then, it follows from (33) that $\widetilde{H}$ is positive definite, if $\sigma_{2}\in(q-1,+\infty)$, $\tau\leq 1$, $\tau+s>0$ and all the matrices $B_{j}$, $j=1,\cdots,q$, have full column rank.

Summarizing the above discussions, the matrix $H$ is positive definite if the parameters $(\sigma_{1},\sigma_{2},\tau,s)$ satisfy (30). $\ \ \ \diamondsuit$

Theorem 3.1

The sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ generated by GS-ADMM satisfy

[TABLE]

where

[TABLE]

Proof

By substituting

[TABLE]

into the following identity

[TABLE]

we have

[TABLE]

Now, notice that

[TABLE]

where the second equality holds by (24) and the fourth equality follows from (28). Then, (39) follows from (26)-(27), (41)-(42) and the definition of $G$ in (40). $\ \ \ \diamondsuit$

Theorem 3.2

The sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ generated by GS-ADMM satisfy

[TABLE]

Proof

Setting $\mathbf{w}={\bf{w}}^{*}\in\mathcal{M}^{*}$ in (39) we have

[TABLE]

The above inequality together with (10) reduces to the inequality (43). $\ \ \ \diamondsuit$

It can be observed that if $\left\|{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}\right\|_{G}^{2}$ is positive, then the inequality (43) implies the contractiveness of the sequence $\{{\bf{w}}^{k}-{\bf{w}}^{*}\}$ under the $H$-weighted norm. However, the matrix $G$ defined in (40) is not necessarily positive definite when $\sigma_{1}\in(p-1,+\infty)$, $\sigma_{2}\in(q-1,+\infty)$ and $(\tau,s)\in{\cal{K}}$. Therefore, it is necessary to estimate the lower bound of $\left\|{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}\right\|_{G}^{2}$ for the sake of the convergence analysis.

3.2 Lower bound of $\left\|{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}\right\|_{G}^{2}$

This subsection provides a lower bound of $\left\|{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}\right\|_{G}^{2}$, for $\sigma_{1}\in(p-1,+\infty)$, $\sigma_{2}\in(q-1,+\infty)$ and $(\tau,s)\in{\cal{K}}$, where ${\cal{K}}$ is defined in (8).

By simple calculations, the $G$ given in (40) can be explicitly written as

[TABLE]

where $H_{\mathbf{x}}$ is defined in (18) and

[TABLE]

In addition, we have

[TABLE]

where $\widetilde{D}$ is defined in (34) and

[TABLE]

Now, we present the following lemma which provides a lower bound of $\|{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}\|_{G}$.

Lemma 5

Suppose $\sigma_{1}\in(p-1,+\infty)$ and $\sigma_{2}\in(q-1,+\infty)$. For the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ generated by GS-ADMM, there exists $\xi_{1}>0$ such that

[TABLE]

Proof

First, it is easy to derive that

[TABLE]

where

[TABLE]

with $H_{\mathbf{x},0}$ and $\widetilde{G}_{0}$ defined in (32) and (45), respectively, and

[TABLE]

In the above matrix $\overline{G}_{0}$, $E$ and $H_{\mathbf{y},0}$ are defined in (38). Hence, we have

[TABLE]

Since $\sigma_{1}\in(p-1,+\infty)$ and $\sigma_{2}\in(q-1,+\infty)$, $H_{\mathbf{x},0}$ and $H_{\mathbf{y},0}$ are positive definite. So, there exists a $\xi_{1}>0$ such that

[TABLE]

In view of the definition of $E$ in (38), we have

[TABLE]

Then, the inequality (46) follows from (61) and (3.2). $\ \ \ \diamondsuit$

Lemma 6

Suppose $\tau>-1$. Then, the sequence $\{{\bf{w}}^{k}\}$ generated by GS-ADMM satisfies

[TABLE]

Proof

It follows from the optimality condition of $y_{l}^{k+1}$-subproblem that $y_{l}^{k+1}\in\mathcal{Y}_{l}$ and for any $y_{l}\in{\cal{Y}}_{l}$, we have

[TABLE]

with $c_{y,l}=c-{\cal{A}}{\bf{x}}^{k+1}-\sum\limits_{j=1,j\neq l}^{q}B_{j}y_{j}^{k}$, which implies

[TABLE]

For $l=1,2,\cdots,q$, letting $y_{l}=y_{l}^{k}$ in the above inequality and summing them together, we can deduce that

[TABLE]

where

[TABLE]

and $H_{\mathbf{y},0}$ is defined in (38). Similarly, it follows from the optimality condition of $y_{l}^{k}$-subproblem that

[TABLE]

For $l=1,2,\cdots,q$, letting $y_{l}=y_{l}^{k+1}$ in the above inequality and summing them together, we obtain

[TABLE]

Since $\sigma_{2}\in(q-1,\infty)$ and all $B_{j}$, $j=1,\cdots,q$, have full column rank, we have from (80) that $H_{\mathbf{y}}$ is positive definite. Meanwhile, by the Cauchy-Schwartz inequality, we get

[TABLE]

By adding (3.2) to (3.2) and using (89), we achieve

[TABLE]

From the update of $\lambda^{k+\frac{1}{2}}$, i.e., $\lambda^{k+\frac{1}{2}}=\lambda^{k}-\tau\beta\left(\mathcal{A}{\bf{x}}^{k+1}+\mathcal{B}{\bf{y}}^{k}-c\right)$ and the update of $\lambda^{k}$, i.e., $\lambda^{k}=\lambda^{k-\frac{1}{2}}-s\beta\left(\mathcal{A}{\bf{x}}^{k}+\mathcal{B}{\bf{y}}^{k}-c\right),$ we have

[TABLE]

Substituting the above inequality into the left-term of (3.2), the proof is completed. $\ \ \ \diamondsuit$

Theorem 3.3

Suppose $\sigma_{1}\in(p-1,+\infty)$, $\sigma_{2}\in(q-1,+\infty)$ and $\tau>-1$. For the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ generated by GS-ADMM , there exists $\xi_{1}>0$ such that

[TABLE]

Proof

The inequality (3.3) is directly obtained from (46) and (6). $\ \ \ \diamondsuit$

The following theorem gives another variant of the lower bound of $\|{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}\|_{G}^{2}$, which plays a key role in showing the convergence of GS-ADMM.

Theorem 3.4

Let the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ be generated by GS-ADMM. Then, for any

[TABLE]

where ${\cal{K}}$ is defined in (8), there exist constants $\xi_{i}(i=1,2)>0$ and $\xi_{j}(j=3,4)\geq 0$, such that

[TABLE]

Proof

By the Cauchy-Schwartz inequality, we have

[TABLE]

Since

[TABLE]

we have $\tau>-1$ when $(\tau,s)\in{\cal{K}}$. Then, combining (3.2) with Theorem 3.3, we deduce

[TABLE]

where $\xi_{1}>0$ is the constant in Theorem 3.3. Since $-1<\tau\leq 1$ and $\beta>0$, we have

[TABLE]

In addition, when $(\tau,s)\in{\cal{K}}$, there is

[TABLE]

which, by $\tau>-1$ and $\beta>0$, implies

[TABLE]

Hence, the proof is completed. $\ \ \ \diamondsuit$

3.3 Global convergence

In this subsection, we show the global convergence and the worst-case $O(1/t)$ convergence rate of GS-ADMM. The following corollary is obtained directly from Theorems 3.1-3.2 and Theorem 3.4.

Corollary 1

Let the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ be generated by GS-ADMM. For any $(\sigma_{1},\sigma_{2},\tau,s)$ satisfying (92), there exist constants $\xi_{i}(i=1,2)>0$ and $\xi_{j}(j=3,4)\geq 0$ such that

[TABLE]

and

[TABLE]

Theorem 3.5

Let the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ be generated by GS-ADMM. Then, for any $(\sigma_{1},\sigma_{2},\tau,s)$ satisfying (92), we have

[TABLE]

[TABLE]

and there exists a ${\bf{w}}^{\infty}\in\mathcal{M}^{*}$ such that

[TABLE]

Proof

Summing the inequality (1) over $k=1,2,\cdots,\infty$, we have

[TABLE]

which implies that (100) and (101) hold since $\xi_{1}>0$ and $\xi_{2}>0$.

Because $(\sigma_{1},\sigma_{2},\tau,s)$ satisfy (92), we have by Lemma 4 that $H$ is positive definite. So, it follows from (1) that the sequence $\{{\bf{w}}^{k}\}$ is uniformly bounded. Therefore, there exits a subsequence $\{{\bf{w}}^{k_{j}}\}$ converging to a point ${\bf{w}}^{\infty}=({\bf{x}}^{\infty},{\bf{y}}^{\infty},\lambda^{\infty})\in{\cal{M}}$. In addition, by the definitions of $\widetilde{x}_{k}$, $\widetilde{y}_{k}$ and $\widetilde{\lambda}_{k}$ in (13) and (14), it follows from (100), (101) and the full column rank assumption of all the matrices $A_{i}$ and $B_{j}$ that

[TABLE]

for all $i=1,\cdots,p$ and $j=1,\cdots,q$. So, we have $\lim\limits_{k\to\infty}{\bf{w}}^{k}-\widetilde{{\bf{w}}}^{k}={\bf 0}$. Thus, for any fixed $\mathbf{w}\in{\cal{M}}$, taking $\widetilde{{\bf{w}}}^{k_{j}}$ in (16) and letting $j$ go to $\infty$, we obtain

[TABLE]

Hence, ${\bf{w}}^{\infty}\in{\cal{M}}^{*}$ is a solution point of $\textrm{VI}(h,\mathcal{J},\mathcal{M})$ defined in (12).

Since (1) holds for any ${\bf{w}}^{*}\in{\cal{M}}^{*}$, by (1) and ${\bf{w}}^{\infty}\in{\cal{M}}^{*}$, for all $l\geq k_{j}$, we have

[TABLE]

This together with (101), (103), $\lim\limits_{j\to\infty}{\bf{w}}^{k_{j}}={\bf{w}}^{\infty}$ and the positive definiteness of $H$ illustrate $\lim\limits_{l\to\infty}{\bf{w}}^{l}={\bf{w}}^{\infty}$. Therefore, the whole sequence $\{{\bf{w}}^{k}\}$ converges to the solution ${\bf{w}}^{\infty}\in{\cal{M}}^{*}$. This completes the whole proof. $\ \ \ \diamondsuit$

The above Theorem 3.5 shows the global convergence of our GS-ADMM. Next, we show the ${\cal{O}}(1/t)$ convergence rate for the ergodic iterates

[TABLE]

Theorem 3.6

Let the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ be generated by GS-ADMM. Then, for any $(\sigma_{1},\sigma_{2},\tau,s)$ satisfying (92), there exist $\xi_{j}(j=3,4)\geq 0$ such that

[TABLE]

Proof

For $k=1,\cdots,t$, summing the inequality (1), we have

[TABLE]

Since $(\sigma_{1},\sigma_{2},\tau,s)$ satisfy (92), $H_{\mathbf{y}}$ is positive definite. And by Lemma 4, $H$ is also positive definite. So, it follows from (3.3) that

[TABLE]

By the convexity of $h$ and (105), we have

[TABLE]

Then, (3.6) follows from (3.3). $\ \ \ \diamondsuit$

Remark 1

In the above Theorem 3.5 and Theorem 3.6, we assume the parameters $(\sigma_{1},\sigma_{2},\tau,s)$ satisfy (92). However, because of the symmetric role played by the $x$ and $y$ iterates in the GS-ADMM, substituting the index $k+1$ by $k$ for the $x$ and $\lambda$ iterates, the GS-ADMM algorithm (7) can be clearly written as

[TABLE]

So, by applying Theorem 3.5 and Theorem 3.6 on the algorithm (109), it also converges and has the ${\cal{O}}(1/t)$ convergence rate when $(\sigma_{1},\sigma_{2},\tau,s)$ satisfy

[TABLE]

where

[TABLE]

Hence, the convergence domain ${\cal{K}}$ in Theorem 3.5 and Theorem 3.6 can be easily enlarged to the symmetric domain, shown in Fig. 2,

[TABLE]

Remark 2

Theorem 3.5 implies that we can apply the following easily usable stopping criterion for GS-ADMM:

[TABLE]

where $tol>0$ is a given tolerance error. One the other hand, Theorem 3.6 tells us that for any given compact set $\mathcal{K}\subset\mathcal{M}$, if

[TABLE]

we have

[TABLE]

which shows that GS-ADMM has a worst-case $\mathcal{O}(1/t)$ convergence rate in an ergodic sense.

4 Two special cases of GS-ADMM

Note that in GS-ADMM (7), the two proximal parameters $\sigma_{1}$ and $\sigma_{2}$ are required to be strictly positive for the generalized $p+q$ block separable convex programming. However, when taking $\sigma_{1}=\sigma_{2}=0$ for the two-block problem, i.e., $p=q=1$, GS-ADMM would reduce to the scheme (5), which is globally convergent HeMaYuan2016 . Such observations motivate us to further investigate the following two special cases:

[TABLE]

4.1 Convergence for the case (a)

The corresponding GS-ADMM for the first case (a) can be simplified as follows:

[TABLE]

And, the corresponding matrices $Q,M,H$ and $G$ in this case are the following:

[TABLE]

where $H_{\mathbf{x}}$ is defined in (18) and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

It can be verified that $H$ in (117) is positive definite as long as

[TABLE]

Analogous to (61), we have

[TABLE]

for some $\xi_{1}>0$, since $H_{\mathbf{x},0}$ defined in (32) is positive definite. When $\sigma_{2}=0$, by a slight modification for the proof of Lemma 6, we have the following lemma.

Lemma 7

Suppose $\tau>-1$. The sequence $\{{\bf{w}}^{k}\}$ generated by the algorithm (113) satisfies

[TABLE]

Then, the following two theorems are similar to Theorem 3.4 and Theorem 3.5.

Theorem 4.1

Let the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ be generated by the algorithm (113). For any

[TABLE]

where ${\cal{K}}$ is given in (8), there exist constants $\xi_{i}(i=1,2)>0$ and $\xi_{3}\geq 0$, such that

[TABLE]

Proof

For any $(\tau,s)\in{\cal{K}}$, we have $\tau>-1$. Then, by Lemma 7, the inequality (119) and the Cauchy-Schwartz inequality (3.2), we can deduce that (4.1) holds with $\xi_{1}$ given in (119), $\xi_{2}$ and $\xi_{3}$ given in (97) and (96), respectively. $\ \ \ \diamondsuit$

Theorem 4.2

Let the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ be generated by the algorithm (113). For any

[TABLE]

where

[TABLE]

we have

[TABLE]

and there exists a ${\bf{w}}^{\infty}\in\mathcal{M}^{*}$ such that $\lim\limits_{k\rightarrow\infty}\widetilde{{\bf{w}}}^{k}={\bf{w}}^{\infty}.$

Proof

First, it is clear that Theorem 3.2 still holds, which, combining with Theorem 4.1, gives

[TABLE]

Hence, (127) holds. For $(\sigma_{1},\tau,s)$ satisfying (126), $H$ in (117) is positive definite. So, by (4.1), $\{{\bf{w}}^{k}\}$ is uniformly bounded and therefore, there exits a subsequence $\{{\bf{w}}^{k_{j}}\}$ converging to a point ${\bf{w}}^{\infty}=({\bf{x}}^{\infty},y^{\infty},\lambda^{\infty})\in{\cal{M}}$. So, it follows from (127) and the full column rank assumption of all the matrices $A_{i}$ that

[TABLE]

for all $i=1,\cdots,p$. Hence, by $\lim\limits_{j\to\infty}{\bf{w}}^{k_{j}}={\bf{w}}^{\infty}$ and (127), we have

[TABLE]

and therefore, by the full column rank assumption of $B$ and (127),

[TABLE]

Hence, by (129), we have $\lim\limits_{j\to\infty}{\bf{w}}^{k_{j}}-\widetilde{{\bf{w}}}^{k_{j}}={\bf{0}}$. Thus, by taking $\widetilde{{\bf{w}}}^{k_{j}}$ in (16) and letting $j$ go to $\infty$, the inequality (104) still holds. Then, the rest proof of this theorem follows from the same proof of Theorem 3.5. $\ \ \ \diamondsuit$

Based on Theorem 4.1 and by a similar proof to those of Theorem 3.6, we can also show that the algorithm (113) has the worst-case ${\cal{O}}(1/t)$ convergence rate, which is omitted here for conciseness.

4.2 Convergence for the case (b)

The corresponding GS-ADMM for the case (b) can be simplified as follows

[TABLE]

In this case, the corresponding matrices $Q,M,H$ and $G$ become $\widetilde{Q},\widetilde{M},\widetilde{H}$ and $\widetilde{G}$, which are defined in (19), the lower-right block of $M$ in (25), (29) and (44), respectively.

In what follows, let us define

[TABLE]

Then, by the proof of Theorem 3.4, we can deduce the following theorem.

Theorem 4.3

Let the sequences $\{\mathbf{v}^{k}\}$ and $\{\widetilde{\mathbf{v}}^{k}\}$ be generated by the algorithm (130). For any

[TABLE]

where ${\cal{K}}$ is defined in (8), there exist constants $\xi_{i}(i=1,2)>0$ and $\xi_{j}(j=3,4)\geq 0$ such that

[TABLE]

By slight modifications of the proof of Theorem 3.5 and Theorem 4.2, we have the following global convergence theorem.

Theorem 4.4

Let the sequences $\{{\bf{w}}^{k}\}$ and $\{\widetilde{{\bf{w}}}^{k}\}$ be generated by the algorithm (113). Then, for any

[TABLE]

where ${\cal{K}}$ is defined in (8), we have

[TABLE]

and there exists a ${\bf{w}}^{\infty}\in\mathcal{M}^{*}$ such that $\lim_{k\rightarrow\infty}\widetilde{{\bf{w}}}^{k}={\bf{w}}^{\infty}.$

By a similar proof to that of Theorem 3.6, the algorithm (130) also has the worst-case ${\cal{O}}(1/t)$ convergence rate.

Remark 3

Again, substituting the index $k+1$ by $k$ for the $x$ and $\lambda$ iterates, the algorithm (113) can be also written as

[TABLE]

So, by applying Theorem 4.4 on the above algorithm, we know the algorithm (113) also converges globally when $(\sigma_{1},\tau,s)$ satisfy

[TABLE]

where $\overline{{\cal{K}}}$ is given in (111). Hence, the convergence domain ${\cal{K}}_{1}$ in Theorem 4.2 can be enlarged to the symmetric domain ${\cal{G}}={\cal{K}}_{1}\cup\overline{{\cal{K}}}$ given in (112). By a similar reason, the convergence domain ${\cal{K}}$ in Theorem 4.4 can be enlarged to ${\cal{G}}$ as well.

5 Numerical experiments

In this section, we investigate the performance of the proposed GS-ADMM for solving a class of sparse matrix minimization problems. All the algorithms are coded and simulated in MATLAB 7.10(R2010a) on a PC with Intel Core i5 processor(3.3GHz) with 4 GB memory.

5.1 Test problem

Consider the following Latent Variable Gaussian Graphical Model Selection (LVGGMS) problem arising in the statistical learning Chandrasekaran2012 ; Ma2017 :

[TABLE]

where $C\in\mathcal{R}^{n\times n}$ is the covariance matrix obtained from observation, $\nu$ and $\mu$ are two given positive weight parameters, $tr(L)$ stands for the trace of the matrix $L$ and $\|S\|_{1}=\sum_{ij}|S_{ij}|$. Clearly, by two different ways of partitioning the variables of (131), the GS-ADMM (7) can lead to the following two algorithms:

[TABLE]

[TABLE]

Note that all the subproblems in (132) and (133) have closed formula solutions. Next, we take the scheme (132) for an example to show how to get the explicit solutions of the subproblem. By the first-order optimality condition of the $X$-subproblem in (132), we derive

[TABLE]

which is equivalent to

[TABLE]

Then, from the eigenvalue decomposition

[TABLE]

where $\textrm{Diag}(\rho)$ is a diagonal matrix with $\rho_{i},i=1,\cdots,n$, on the diagonal, we obtain that the solution of the $X$-subproblem in (132) is

[TABLE]

where $\textrm{Diag}(\gamma)$ is the diagonal matrix with diagonal elements

[TABLE]

On the other hand, the $S$-subproblem in (132) is equivalent to

[TABLE]

where $\textrm{Shrink}(\cdot,\cdot)$ is the soft shrinkage operator (see e.g.TaoYuan2014 ). Meanwhile, it is easy to verify that the $L$-subproblem is equivalent to

[TABLE]

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

[TABLE]

5.2 Numerical results

In the following, we investigate the performance of several algorithms for solving the LVGGMS problem, where all the corresponding subproblems can be solved in a similar way as shown in the above analysis. For all algorithms, the maximal number of iterations is set as $1000$, the starting iterative values are set as $(X^{0},S^{0},L^{0},\Lambda^{0})=(\textbf{I},\textbf{2I},\textbf{I},\textbf{0})$, and motivated by Remark 2, the following stopping conditions are used

[TABLE]

together with $\textrm{CER(k)}:=\left\|X^{k}-S^{k}+L^{k}\right\|_{F}\leq 10^{-4}$, where $F^{*}$ is the approximate optimal objective function value obtained by running GS-ADMM (132) after 1000 iterations. In (131), we set $(\nu,\mu)=(0.005,0.05)$ and the given data $C$ is randomly generated by the following MATLAB code with $m=100$, which are downloaded from S. Boyd’s homepage222http://web.stanford.edu/$\sim$boyd/papers/admm/covsel/covsel_example.html:

[TABLE]

5.2.1 Performance of different versions of GS-ADMM

In the following, we denote

[TABLE]

[TABLE]

Table 1: Numerical results of GS-ADMM-I, II, III and IV with different $\beta$.

First, we would like to investigate the performance of the above different versions of GS-ADMM for solving the LVGGMS problem with variance of the penalty parameter $\beta$. The results are reported in Table 1 with $\textrm{TOL}=\textrm{Tol}=1.0\times 10^{-7}$, and $(\tau,s)=(0.8,1.17)$. For GS-ADMM-I and GS-ADMM-II, $(\sigma_{1},\sigma_{2})=(2,3)$. Here, “Iter” and “CPU” denote the iteration number and the CPU time in seconds, and the bold letter indicates the best result of each algorithm. From Table 1, we can observe that:

• •

Both the iteration number and the CPU time of all the algorithms have a similar changing pattern, which decreases originally and then increases along with the decrease of the value of $\beta$.

• •

For the same value of $\beta$, the results of GS-ADMM-III are slightly better than other algorithms in terms of the iteration number, CPU time, and the feasibility errors CER, IER and OER.

• •

GS-ADMM-III with $\beta=0.5$ can terminate after 579 iterations to achieve the tolerance $10^{-7}$, while the other algorithms with $\beta=0.5$ fail to achieve this tolerance within given number of iterations.

In general, the algorithm (132) with $\beta=0.06$ performs better than other cases. Hence, in the following experiments for GS-ADMM, we adapt GS-ADMM-III with default $\beta=0.06$. Also note that $\sigma_{2}=0$, which is not allowed by the algorithms discussed in Heyuan2015 .

[TABLE]

Table 2: Numerical results of GS-ADMM-III with different stepsizes $(\tau,s)$.

Second, we investigate how the stepsizes $(\tau,s)\in{\cal{G}}$ with different values would affect the performance of GS-ADMM-III. Table 2 reports the comparison results with variance of $(\tau,s)$ for $\textrm{TOL}=\textrm{Tol}=1.0\times 10^{-5}$. One obvious observation from Table 2 is that both the iteration number and the CPU time decrease along with the increase of s (or $\tau$) for fixed value of $\tau$ (or s), which indicates that the stepsizes of $(\tau,s)\in{\cal{G}}$ could influence the performance of GS-ADMM significantly. In addition, the results in Table 2 also indicate that using more flexible but with both relatively larger stepsizes $\tau$ and $s$ of the dual variables often gives the best convergence speed. Comparing all the reported results in Table 2, by setting $(\tau,s)=(0.9,1.09)$, GS-ADMM-III gives the relative best performance for solving the problem (131).

5.2.2 Comparison of GS-ADMM with other state-of-the-art algorithms

In this subsection, we would like to carry out some numerical comparison of solving the problem (131) by using GS-ADMM-III and the other four methods:

[TABLE]

Table 3: Comparative results of different algorithms under different tolerances.

• •

The Proximal Jacobian Decomposition of ALM HeXuYuan2016 (denoted by “PJALM”);

• •

The splitting method in HeTaoYuan2015 (denoted by “HTY”);

• •

The generalized parametrized proximal point algorithm BaiLiLi2017 (denoted by “GR-PPA”).

• •

The twisted version of the proximal ADMM WangSong2017 (denoted by “T-ADMM”).

We set $(\tau,s)=(0.9,1.09)$ for GS-ADMM-III and the parameter $\beta=0.05$ for all the comparison algorithms. The default parameter $\mu=2.01$ and $H=\beta{\bf{I}}$ are used for HTY HeTaoYuan2015 . As suggested by the theory and numerical experiments in HeXuYuan2016 , the proximal parameter is set as 2 for PJALM. As shown in BaiLiLi2017 , the relaxation factor of GR-PPA is set as 1.8 and other default parameters are chosen as

[TABLE]

For T-ADMM, the symmetric matrices therein are chosen as $M_{2}=M_{2}=v\textbf{I}$ with $v=\beta$ and the correction factor is set as $a=1.6$ WangSong2017 . The results obtained by the above algorithms under different tolerances are reported in Table 3. With fixed tolerance $\textrm{TOL}=10^{-9}$ and $\textrm{Tol}=10^{-15}$, the convergence behavior of the error measurements IER(k) and OER(k) by the five algorithms using different starting points are shown in Figs. 3-5. From Table 3 and Figures 3-5, we may have the following observation:

• •

Under all different tolerances, GS-ADMM-III performs significantly better than other four algorithms in both the number of iterations and CPU time.

• •

GR-PPA is slightly better than PJALM and HTY, and T-ADMM outperforms PJALM, HTY and GR-PPA.

• •

the convergence curves in Figs. 3-5 illustrate that using different starting points, GS-ADMM-III also converges fastest among the comparing methods.

All these numerical results demonstrate the effectiveness and robustness of GS-ADMM-III, which is perhaps due to the symmetric updating of the Lagrange multipliers and the proper choice of the stepsizes.

6 Conclusion

Since the direct extension of ADMM in a Gauss-Seidel fashion for solving the three-block separable convex optimization problem is not necessarily convergent analyzed by Chen et al. Chen2016 , there has been a constantly increasing interest in developing and improving the theory of the ADMM for solving the multi-block separable convex optimization. In this paper, we propose an algorithm, called GS-ADMM, which could solve the general model (1) by taking advantages of the multi-block structure. In our GS-ADMM, the Gauss-Seidel fashion is taken for updating the two grouped variables, while the block variables within each group are updated in a Jacobi scheme, which would make the algorithm more be attractive and effective for solving big size problems. We provide a new convergence domain for the stepsizes of the dual variables, which is significantly larger than the convergence domains given in the literature. Global convergence as well as the ${\cal{O}}(1/t)$ ergodic convergence rate of the GS-ADMM is established. In addition, two special cases of GS-ADMM, which allows one of the proximal parameters to be zero, are also discussed.

This paper simplifies the analysis in HeMaYuan2016 and provides an easy way to analyze the convergence of the symmetric ADMM. Our preliminary numerical experiments show that with proper choice of parameters, the performance of the GS-ADMM could be very promising. Besides, from the presented convergence analysis, we can see that the theories in the paper can be naturally extended to use more general proximal terms, such as letting $P_{i}^{k}(x_{k}):=\frac{\beta}{2}\|x_{i}-x_{i}^{k}\|_{{\bf{P}}_{i}}$ and $Q_{j}^{k}(y_{j}):=\frac{\beta}{2}\|y_{j}-y_{j}^{k}\|_{{\bf{Q}}_{j}}$ in (7), where ${\bf{P}}_{i}$ and ${\bf{Q}}_{i}$ are matrices such that ${\bf{P}}_{i}\succ(p-1)A_{i}^{\sf T}A_{i}$ and ${\bf{Q}}_{j}\succ(q-1)B_{j}^{\sf T}B_{j}$ for all $i=1,\cdots,p$ and $j=1,\cdots,q$. Finally, the different ways of partitioning the variables of the problem also gives the flexibility of GS-ADMM.

Acknowledgements. The authors would like to thank the anonymous referees for providing very constructive comments. The first author also wish to thank Prof. Defeng Sun at National University of Singapore for his valuable discussions on ADMM and Prof. Pingfan Dai at Xi’an Jiaotong University for discussion on an early version of the paper.