An Improvement of the Alternating Direction Method of Multipliers to Solve the Convex Optimization Problem

Many problems in the fields of signal and image processing, machine learning, medical image reconstruction, computer vision, and network communications [1,2,3,4,5,6,7,8] can be reduced to solving the following convex optimization problem with linear constraints: where $A\in {ℝ}^{m\times m_1}$ and $B\in {ℝ}^{m\times m_2}$ are the given matrices; $b\in {ℝ}^m$ is the given vector; $f:{ℝ}^{m_1}\to ℝ$ and $g:{ℝ}^{m_2}\to ℝ$ are continuous closed convex functions; and $\mathcal{X}$ and $\mathcal{Y}$ are nonempty closed convex subsets in ${ℝ}^{m_1}$ and ${ℝ}^{m_2}$, respectively. The augmented Lagrangian function of the convex optimization problem (Equation (1)) is the following Equation (2). where $\beta >0$ is the penalty parameter and $\lambda \in {ℝ}^m$ is the Lagrange multiplier vector. Based on the objective function with separable structural properties, the alternating direction method of multipliers (ADMM) has been proposed in the literature [9,10] to solve the problem (Equation (1)) as follows:

The ADMM algorithm can be viewed as the Gauss Seidel form of the augmented Lagrangian multiplier method, and it can be interpreted as the Douglas Rachford splitting algorithm (DRSM) from a dual perspective [11]. In addition, the DRSM can be further interpreted as a proximal point algorithm (PPA) [12]. The ADMM algorithm is an important method for solving separable convex optimization problems. Due to its fast processing speed and good convergence performance, the ADMM algorithm is widely used in statistical learning, machine learning, and other fields. However, the ADMM algorithm has slow convergence speed under the requirement of a large scale and high precision. Therefore, many experts and scholars are working on improving the ADMM algorithm, such as Fortin and Glowinski, who suggested in [13] attaching a relaxation factor to the Lagrange multiplier updating step in Equation (3), which this results in the following scheme: where the parameter γ can be chosen in the interval $\left(0,\frac{1+\sqrt{5}}{2}\right)$, and thus it becomes possible to enlarge the step size for updating the Lagrange multiplier. An advantage of this larger step size in Equation (4) is that it can numerically lead to faster convergence empirically. This scheme (Equation (4)) differs from the original ADMM scheme (Equation (3)) only in the fact that the step size for updating the Lagrange multiplier can be larger than 1. But technically they are actually two distinct families of ADMM algorithms, in which one is derived from the operator splitting framework and the other is derived from Lagrangian splitting. Thus, despite the similarity in notation, the ADMM scheme (Equation (4)) with Fortin et al.’s larger step size and the original ADMM scheme (Equation (3)) are completely different in nature. On the other hand, Glowinski et al. [14] applied the Douglas Rachford splitting method (DRSM) [11] to the dual of Equation (1) and obtained the following scheme:

This scheme (Equation (5)) can be regarded as a symmetric version of the ADMM scheme (Equation (3)) in the sense that the variables $x$ and $y$ are treated equally, each of which is followed consequently by an update of the Lagrange multiplier. However, as shown in [15], the sequence generated by the symmetric ADMM (Equation (5)) is not necessarily strictly contractive with respect to the solution set of Equation (1), while this property can be ensured by the sequence generated by the ADMM (Equation (3)). Because of these deficiencies, Bingsheng He et al. [15] proposed the following scheme: where the parameter α∈ (0, 1) is for shrinking the step sizes in Equation (5). The sequence generated by Equation (6) is strictly contractive with respect to the solution set of Equation (1). Thus, we called it the strictly contractive symmetric version of the ADMM. Limiting α ∈ (0, 1) in the strictly contractive symmetric version of the ADMM (Equation (6)) makes the updating of Lagrange multipliers more conservative and with smaller steps, but smaller steps should be strongly avoided in practical applications. Instead, one wants to seek larger steps where possible to speed up numerical performance. For this purpose, Bingsheng He et al. [16] proposed the following scheme: where $r$ and $s$ are limited to the following region

As for more symmetric alternation method improvements of the ADMM algorithm, we provide the reader with references [17,18,19,20,21,22,23,24,25,26].

In addition, the convergence speed of the ADMM algorithm is directly related to the choice of the penalty parameter $\beta$. Therefore, improvements of the ADMM algorithm with a change to penalty parameter ${\beta }_k$ instead of the fixed penalty parameter $\beta$ are proposed. Such as He, Yang, and Wang [27] proposed an adaptive penalty parameter scheme to solve linearly constrained monotone variational inequalities. In this paper, we, based on the idea of references [16,27], propose an adaptive penalty parameter alternate direction method of multipliers with attach relaxation factors and with the Lagrange multiplier updated twice at each iteration. The iteration format is shown in Algorithm 1.

Algorithm 1: Improved alternate direction multiplier algorithm to solve convex optimization problem (1)
Given error tolerance $\epsilon >0$, controls parameter $\mu \in \left(0,1\right)$, and relaxation factors $(r,s)\in \mathcal{D}$ defined in (Equation (8)). Choose parameter sequence $\left\{{\tau }_k\right\}$, which satisfies ${\tau }_k\ge 0$ and ${\displaystyle \sum _{k=0}^{\infty }{\tau }_k}<+\infty$. Given initial penalty parameter ${\beta }_0>0$ and initial approximation $\left(x_0,y_0,{\lambda }_0\right)\in \left(\mathcal{X},\mathcal{Y},{ℝ}^r\right)$. Set $k=0$. Step 1. Computing
$x_{k+1}=\arg \underset{x\in \mathcal{X}}{\min }\left\{f(x)-{\left({\lambda }_k\right)}^T\left(Ax+B{y}_k-b\right)+\frac{{\beta }_k}{2}{‖Ax+B{y}_k-b‖}_2^2\right\}$,	(9)
${\lambda }_{k+\frac{1}{2}}={\lambda }_k-r{\beta }_k\left(A{x}_{k+1}+B{y}_k-b\right)$,	(10)
$y_{k+1}=\arg \underset{y\in \mathcal{Y}}{\min }\left\{g(y)-{\left({\lambda }_{k+\frac{1}{2}}\right)}^T\left(A{x}_{k+1}+By-b\right)+\frac{{\beta }_k}{2}{‖A{x}_{k+1}+By-b‖}_2^2\right\},$	(11)
${\lambda }_{k+1}={\lambda }_{k+\frac{1}{2}}-s{\beta }_k\left(A{x}_{k+1}+B{y}_{k+1}-b\right)$.	(12)
Step 2. If $\max \left({‖A{x}_{k+1}+B{y}_{k+1}-b‖}_2,{‖B\left(y_k-y_{k+1}\right)‖}_2\right)\le \epsilon$, stop. Otherwise, go to Step 3. Step 3. Updating penalty parameter
${\beta }_{k+1}=\left\{\begin{array}{ll}\left(1+{\tau }_k\right){\beta }_k,& if {‖x_k-P_{\mathcal{X}}\left[x_k-\left(\nabla f\left(x_k\right)-A^T{\lambda }_k\right)\right]‖}_2<\mu {‖A{x}_k+B{y}_k-b‖}_2,\\ {\beta }_k/\left(1+{\tau }_k\right),& if \mu {‖x_k-P_{\mathcal{X}}\left[x_k-\left(\nabla f\left(x_k\right)-A^T{\lambda }_k\right)\right]‖}_2>{‖A{x}_k+B{y}_k-b‖}_2,\\ {\beta }_k,& \mathrm{otherwise}.\end{array}\right.$	(13)
Step 4. Set $k=k+1$, and go to Step 1.

We choose the algorithm described in [28] with the necessary modifications to solve the x-subproblem and y-subproblem in Algorithm 1. The iteration method to solve the x-subproblem can be described as follows in Algorithm 2.

Algorithm 2: Algorithm to solve x-subproblem (9) or y-subproblem (11)

Given constants $\theta \in (0,1),\nu \in (0,1)$, $\mathcal{l}\in (1,\infty )$, error tolerance $\epsilon$, and initial approximation $x_0\in \mathcal{X}$. Set ${\eta }_0=0.9$ and $k=0$. Step 1. Computing ${\overline{x}}_{k+1}=P_{\mathcal{X}}\left[x_k-{\eta }_{k}F\left(x_k\right)\right]$, where $F\left(x\right)=\nabla f\left(x\right)-A^T{\lambda }_k+{\beta }_k{A}^T\left(Ax+B{y}_k-b\right)$. Step 2. If ${\psi }_k={{\eta }_k{‖F\left(x_k\right)-F\left({\overline{x}}_{k+1}\right)‖}_2/‖x_k-{\overline{x}}_{k+1}‖}_2\le \nu$, computing $x_{k+1}=P_{\mathcal{X}}\left[x_k-{\rho }_{k}F\left({\overline{x}}_{k+1}\right)\right],$   where ${\rho }_k=\mathcal{l}{\eta }_k{e}_k^T{d}_k/{‖d_k‖}_2^2$, $e_k=x_k-{\overline{x}}_{k+1}$, $d_k=e_k-{\eta }_k[F(x_k)-F({\overline{x}}_{k+1})]$.   If ${‖x_{k+1}-P_{\mathcal{X}}\left[x_{k+1}-F\left(x_{k+1}\right)\right]‖}_2\le \epsilon$, stop (in this case, $x_{k+1}$ is an approximate solution of the x-subproblem in Algorithm 1). Otherwise, define     ${\eta }_{k+1}\triangleq \left\{\begin{array}{ll}\left(3/2\right){\eta }_k,& if {\psi }_k\le \theta \\ {\eta }_k,& \mathrm{otherwise}\end{array}\right.$ and set $k=k+1$. Go to Step 1.   If all of the above are not true, go to Step 3. Step 3. Define ${\eta }_k\triangleq \left(2/3\right){\eta }_k\ast \min \left\{1,1/{\psi }_k\right\}$ and set ${\overline{x}}_{k+1}=P_{\mathcal{X}}\left[x_k-{\eta }_{k}F\left(x_k\right)\right]$. Go to Step 2.

The rest of this paper is organized as follows: In Section 2, the global convergence of the proposed algorithm is proven. In Section 3, some numerical comparison experiments with existing algorithms are given. Finally, we draw some conclusions.

To conveniently prove the global convergence of Algorithm 1, we first give the following Lemmas 1–3.

Let then $H_k$ is a symmetric semi-definite matrix, and defining we have the following Lemma 7.

To prove the global convergence of Algorithm 1, we divided the domain $\mathcal{D}$ defined in Equation (8) into the following five parts:

Obviously,