A family of multi-parameterized proximal point algorithms

TL;DR

This paper introduces a new multi-parameterized proximal point algorithm with relaxation for convex optimization with linear constraints, demonstrating improved convergence and performance over existing methods.

Contribution

The paper develops a novel multi-parameterized proximal point algorithm with relaxation, providing convergence analysis and empirical evidence of superior performance.

Findings

Proven global convergence of the algorithm.

Achieved sublinear convergence rate.

Numerical experiments show better performance than existing methods.

Abstract

In this paper, a multi-parameterized proximal point algorithm combining with a relaxation step is developed for solving convex minimization problem subject to linear constraints. We show its global convergence and sublinear convergence rate from the prospective of variational inequality. Preliminary numerical experiments on testing a sparse minimization problem from signal processing indicate that the proposed algorithm performs better than some well-established methods

Equations

$$\begin{array}[]{lll}\min\{f(\mathbf{x})\ |A\mathbf{x}=b,\mathbf{x}\in\mathcal{X}\},\end{array}$$

$$\left\{\begin{array}[]{l}\mathbf{x}^{k+1}=\arg\min\limits_{\mathbf{x}\in\mathcal{X}}f(\mathbf{x})+\frac{\beta}{2}\|A\mathbf{x}-b-\lambda^{k}/\beta\|^{2},\\ \lambda^{k+1}=\lambda^{k}-\beta(A\mathbf{x}^{k+1}-b),\end{array}\right.$$

$$\text{VI}(f,\mathcal{J},\mathcal{M}):f(\mathbf{x})-f({\mathbf{x}}^{\ast })+(w-w^{\ast }{)}^{\mathsf{T}}\mathcal{J}(w^{\ast })\ge 0,\forall w\in \mathcal{M},$$

$$w=\left(\begin{array}[]{c}\mathbf{x}\\ \lambda\\ \end{array}\right),\quad w^{*}=\left(\begin{array}[]{c}\mathbf{x}^{*}\\ \lambda^{*}\\ \end{array}\right),\quad\mathcal{J}(w)=\left(\begin{array}[]{c}-A^{\sf T}\lambda\\ A\mathbf{x}-b\end{array}\right)\quad\textrm{and}\quad\mathcal{M}=\mathcal{X}\times\mathcal{R}^{m}.$$

$$f(\mathbf{x})-f({\mathbf{x}}^{k+1})+(w-w^{k+1}{)}^{\mathsf{T}}\left[\mathcal{J}(w^{k+1})+G(w^{k+1}-w^k)\right]\ge 0,\forall w\in \mathcal{M}.$$

$$G=\left[\begin{array}[]{cc}rI_{n}+\frac{(\theta-1)^{2}-\rho}{s}A^{\sf T}A&(\theta-1)A^{\sf T}\\ (\theta-1)A&sI_{m}\end{array}\right]\in\mathcal{R}^{(n+m)\times(n+m)},$$

$$\rho \in (-\infty ,1],(r,s)\in \left\{(r,s)|r>0,s>0,rs>\Vert A^{\mathsf{T}}A{\Vert }_2\right\}.$$

$$\left\{\begin{array}[]{l}A\mathbf{x}^{k+1}-b+(\theta-1)A(\mathbf{x}^{k+1}-\mathbf{x}^{k})+s(\lambda^{k+1}-\lambda^{k})=0,\\ \mathbf{x}^{k+1}\in\mathcal{X},\quad f(\mathbf{x})-f(\mathbf{x}^{k+1})+(\mathbf{x}-\mathbf{x}^{k+1})^{\sf T}R_{k+1}\geq 0,\quad\forall\mathbf{x}\in\mathcal{X},\end{array}\right.$$

$$R_{k+1}=-A^{\mathsf{T}}{\lambda }^{k+1}+\left[r{I}_n+\frac{(\theta -1{)}^2-\rho }{s}{A}^{\mathsf{T}}A\right]({\mathbf{x}}^{k+1}-{\mathbf{x}}^k)+(\theta -1)A^{\mathsf{T}}({\lambda }^{k+1}-{\lambda }^k).$$

$${\lambda }^{k+1}={\lambda }^k-\frac{1}{s}\left[\theta (A{\mathbf{x}}^{k+1}-b)+(1-\theta )(A{\mathbf{x}}^k-b)\right],$$

$$R_{k+1}$$

$${\mathbf{x}}^{k+1}=\arg \underset{\mathbf{x}\in \mathcal{X}}{\min }\left\{f(\mathbf{x})+\frac{r}{2}{\Vert \mathbf{x}-{\mathbf{x}}^k-\frac{1}{r}{A}^{\mathsf{T}}\left[{\lambda }^k-\frac{2-\theta }{s}(A{\mathbf{x}}^k-b)\right]\Vert }^2+\frac{\rho -1}{2s}{\Vert A(\mathbf{x}-{\mathbf{x}}^k)\Vert }^2\right\}.$$

$$f(\mathbf{x})-f({\tilde{\mathbf{x}}}^k)+(w-{\tilde{w}}^k{)}^{\mathsf{T}}\left[\mathcal{J}({\tilde{w}}^k)+G({\tilde{w}}^k-w^k)\right]\ge 0,\forall w\in \mathcal{M}.$$

$$\frac{\rho -1}{2s}{\Vert A(\mathbf{x}-{\mathbf{x}}^k)\Vert }^2=\frac{\rho -1}{2s}{\Vert (A\mathbf{x}-b)-(A{\mathbf{x}}^k-b)\Vert }^2$$

$$\left\{\begin{array}[]{l}\mathbf{x}^{k+1}=\arg\min\limits_{\mathbf{x}\in\mathcal{X}}\left\{f(\mathbf{x})+\frac{1}{2s}\|A\mathbf{x}-b-s\lambda^{k}\|^{2}+\frac{1}{2}\|\mathbf{x}-\mathbf{x}^{k}\|^{2}_{rI_{n}-\frac{1}{s}A^{\sf T}A}\right\},\\ \lambda^{k+1}=\lambda^{k}-\frac{1}{s}(A\mathbf{x}^{k+1}-b).\end{array}\right.$$

$$\Vert w^{k+1}-w^{\ast }{\Vert }_{\tilde{G}}^2+\frac{1-\rho }{s}\Vert A{\mathbf{x}}^{k+1}-A{\mathbf{x}}^{\ast }{\Vert }^2\le \Vert w^k-w^{\ast }{\Vert }_{\tilde{G}}^2+\frac{1-\rho }{s}\Vert A{\mathbf{x}}^k-A{\mathbf{x}}^{\ast }{\Vert }^2-T_{k+1},\forall w^{\ast }\in {\mathcal{M}}^{\ast },$$

$$\widetilde{G}=\left[\begin{array}[]{cc}rI_{n}+\frac{(\theta-1)^{2}-1}{s}A^{\sf T}A&(\theta-1)A^{\sf T}\\ (\theta-1)A&sI_{m}\end{array}\right].$$

$$(w-\bar{w}{)}^{\mathsf{T}}[\mathcal{J}(w)-\mathcal{J}(\bar{w})]\equiv 0,\text{for any }w,\bar{w}\in \mathcal{M},$$

$${\tilde{w}}^k-w^k=(w^{k+1}-w^k)/\sigma ,$$

$$(w^{\ast }-{\tilde{w}}^k{)}^{\mathsf{T}}G(w^{k+1}-w^k)\ge 0.$$

$$(w^{\ast }-{\tilde{w}}^k{)}^{\mathsf{T}}\tilde{G}(w^{k+1}-w^k)+\frac{1-\rho }{s}(A{\mathbf{x}}^{\ast }-A{\tilde{\mathbf{x}}}^k{)}^{\mathsf{T}}(A{\mathbf{x}}^{k+1}-A{\mathbf{x}}^k)\ge 0.$$

$$2(a-l)\tilde{G}(c-d)=\Vert a-d{\Vert }_{\tilde{G}}^2-\Vert a-c{\Vert }_{\tilde{G}}^2+\Vert c-l{\Vert }_{\tilde{G}}^2-\Vert d-l{\Vert }_{\tilde{G}}^2$$

$$\Vert w^{\ast }-w^{k+1}{\Vert }_{\tilde{G}}^2+\frac{1-\rho }{s}\Vert A{\mathbf{x}}^{\ast }-A{\mathbf{x}}^{k+1}{\Vert }^2\le \Vert w^{\ast }-w^k{\Vert }_{\tilde{G}}^2+\frac{1-\rho }{s}\Vert A{\mathbf{x}}^{\ast }-A{\mathbf{x}}^k{\Vert }^2-T_{k+1},$$

$$T_{k+1}=\Vert w^k-{\tilde{w}}^k{\Vert }_{\tilde{G}}^2-\Vert w^{k+1}-{\tilde{w}}^k{\Vert }_{\tilde{G}}^2+\frac{1-\rho }{s}\left(\Vert A({\mathbf{x}}^k-{\tilde{\mathbf{x}}}^k){\Vert }^2-\Vert A({\mathbf{x}}^{k+1}-{\tilde{\mathbf{x}}}^k){\Vert }^2\right).$$

$$T_{k+1}$$

$$f({\mathbf{x}}_t)-f(\mathbf{x})+({\mathbf{w}}_t-w{)}^{\mathsf{T}}\mathcal{J}(w)\le \frac{1}{2\sigma (t+1)}\left\{\Vert w^0-w{\Vert }_{\tilde{G}}^2+\frac{1-\rho }{s}\Vert A{\mathbf{x}}^0-A\mathbf{x}{\Vert }^2\right\},\forall w\in \mathcal{M}.$$

$${\gamma }_{\xi }=\underset{{\mathbf{x}}^{\ast }\in {\mathcal{X}}^{\ast }}{\inf }\underset{\lambda \in {\Gamma }_{\xi }}{\sup }{\Vert ({\mathbf{x}}^0-{\mathbf{x}}^{\ast };{\lambda }^0-\lambda )\Vert }_{\tilde{G}}^2+\frac{1-\rho }{s}\Vert A{\mathbf{x}}^0-b{\Vert }^2.$$

$$f({\mathbf{x}}_t)-f({\mathbf{x}}^{\ast })+\xi \Vert A{\mathbf{x}}_t-b\Vert \le \frac{{\gamma }_{\xi }}{2\sigma (t+1)},\forall {\mathbf{x}}^{\ast }\in {\mathcal{X}}^{\ast }.$$

$$f({\mathbf{x}}_t)-f({\mathbf{x}}^{\ast })+({\mathbf{w}}_t-w{)}^{\mathsf{T}}\mathcal{J}(w)$$

$$f({\mathbf{x}}_t)-f({\mathbf{x}}^{\ast })+\xi \Vert A{\mathbf{x}}_t-b\Vert$$

Full text

A family of multi-parameterized proximal point algorithms

††thanks: The work was supported by the Natural Science Foundation of China (11801455; 11571178) and the Fundamental Research Funds of China West Normal University (17E084; 18B031).

Jianchao Bai 111Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710129, China (bjc1987@163.com).  Ke Guo 222 (Corresponding author) School of Mathematics and Information, China West Normal University, Nanchong 637002, China (keguo2014@126.com).  Xiaokai Chang 333 School of Science, Lanzhou University of Technology, Lanzhou 730050, China (xkchang@lut.cn).

Abstract

In this paper, a multi-parameterized proximal point algorithm combining with a relaxation step is developed for solving convex minimization problem subject to linear constraints. We show its global convergence and sublinear convergence rate from the prospective of variational inequality. Preliminary numerical experiments on testing a sparse minimization problem from signal processing indicate that the proposed algorithm performs better than some well-established methods.

Keywords: Convex optimization, proximal point algorithm, complexity, signal processing

Mathematics Subject Classification(2010): 65Y20, 90C25, 92C55

1 Introduction

We focus on the following convex minimization problem with linear equality constraints,

[TABLE]

where $f:\mathcal{R}^{n}\rightarrow\mathcal{R}\cup\{+\infty\}$ is a proper closed convex function but possibly nonsmooth, $A\in\mathcal{R}^{m\times n}$ and $b\in\mathcal{R}^{m}$ are given matrix and vector, respectively, $\mathcal{X}\subseteq\mathcal{R}^{n}$ is a closed convex set. Without loss of generality, the solution set of the problem (1) denoted by $\mathcal{X}^{*}$ is assumed to be nonempty.

The augmented Lagrangian method (ALM), independently proposed by Hestenes [8] and Powell [12], is a benchmark method for solving problem (1). Its iteration scheme reads as

[TABLE]

where $\beta,\lambda$ denote the penalty parameter and the Lagrange multiplier w.r.t. the equality constraint, respectively. As analyzed in [13], ALM can be viewed as an application of the well-known proximal point algorithm (PPA) that can date back to the seminal work of Martinet [11] and Rockafellar [14] for the dual problem of (1). Obviously, the efficiency of ALM heavily depends on the solvability of the $\mathbf{x}$-subproblem, that is, whether or not the core $\mathbf{x}$-subproblem has closed-form solution. Unfortunately, in many real applications [2, 4, 9, 10], the coefficient matrix $A$ is not identity matrix (or does not satisfies $AA^{\sf T}=I_{m}$), which makes it difficult even infeasible for solving this subproblem of ALM. To overcome such difficulty, Yang and Yuan [15] proposed a linearized ALM aiming at linearizing the $\mathbf{x}$-subproblem such that its closed-form solution can be easily derived. We refer to the recent progress on this direction [3, 6].

Under basic regularity condition $\mbox{ri~{}dom}(f)\cap X\neq\emptyset$, it is well-known that $x^{*}$ is an optimal solution of (1) if and only if there exists $\lambda^{*}\in\mathcal{R}^{m}$ such that the following variational inequality holds

[TABLE]

where

[TABLE]

From the aforementioned assumption on the solution set of the problem (1), the solution set of (2) denoted by $\mathcal{M}^{*}$ is also nonempty. When PPA is applied to solve the variational inequality $\textrm{VI}(f,\mathcal{J},\mathcal{M})$, it usually reads the unified updating scheme: at the $k$-th iteration, find $w^{k+1}$ satisfying

[TABLE]

We call $G$ the proximal matrix that is usually required to be symmetric positive definite to ensure the convergence of (3). To our knowledge, this idea was initialized by He et al. [HLHY02]. Clearly, different structures of $G$ would result in different versions of PPA. From a computational perspective, our motivation is to design a multi-parameterized PPA for solving problem (1) while maintaining the efficiency as the linearized ALM, although the feasible starting point may be different. Interestingly, many customized proximal matrices shown in [5, 7, 10, 16] turn out to be special cases of our multi-parameterized proximal matrix (See Remark 2.2 for details). In this sense, our proposed algorithm can be viewed as a general customized PPA for solving problem (1). Moreover, we adopt a relaxation strategy to accelerate the convergence.

2 Main Algorithm

In this paper, we design the following multi-parameterized proximal matrix

[TABLE]

where $I_{m}\in\mathcal{R}^{m\times m}$ denotes the identity matrix, $\theta$ is an arbitrary real scalar and

[TABLE]

The notation $\|A^{\sf T}A\|_{2}=\sqrt{\lambda_{\max}(A^{\sf T}A)}$ represents the spectral norm of $A^{\sf T}A.$ It is easy to check that the above matrix $G$ is symmetric positive definite for any parameters $(\rho,r,s)$ satisfying (5).

Now, substituting the matrix $G$ into (3) we have

[TABLE]

with

[TABLE]

By the equation in (6), it can be deduced that

[TABLE]

which further makes (7) become

[TABLE]

Based on the inequality in (6), i.e., the first-order optimality condition of $\mathbf{x}$-subproblem, we obtain

[TABLE]

Then, our relaxed multi-parameterized PPA (RM-PPA) is described as Algorithm 2.1, where we use $\widetilde{w}^{k}:=(\widetilde{\mathbf{x}}^{k},\widetilde{\lambda}^{k})$ to replace the output of (3) with given iterate $(\mathbf{x}^{k},\lambda^{k})$, and we use $(\mathbf{x}^{k+1},\lambda^{k+1})$ to stand for the new iterate after combining a relaxation step. Finally, the inequality (3) becomes

[TABLE]

Algorithm 2.1

*[RM-PPA for solving problem (1)]1 Choose $\sigma\in(0,2),\theta\in\mathcal{R}$ and $(\rho,r,s)$ satisfying (5).

2 Initialize $(\mathbf{x}^{0},\lambda^{0})\in\mathcal{M}.$

3 for $k=1,2,\cdots,$ do

4   $\widetilde{\mathbf{x}}^{k}=\arg\min\limits_{\mathbf{x}\in\mathcal{X}}\left\{f(\mathbf{x})+\frac{r}{2}\left\|\mathbf{x}-\mathbf{x}^{k}-\frac{1}{r}A^{\sf T}\left[\lambda^{k}-\frac{2-\theta}{s}(A\mathbf{x}^{k}-b)\right]\right\|^{2}+\frac{\rho-1}{2s}\left\|A(\mathbf{x}-\mathbf{x}^{k})\right\|^{2}\right\}$.

5   $\widetilde{\lambda}^{k}=\lambda^{k}-\frac{1}{s}\left[\theta(A\widetilde{\mathbf{x}}^{k}-b)+(1-\theta)(A\mathbf{x}^{k}-b)\right]$.

6   \left(\begin{array}[]{c}\mathbf{x}^{k+1}\\ \lambda^{k+1}\\ \end{array}\right)=\left(\begin{array}[]{c}\mathbf{x}^{k}\\ \lambda^{k}\\ \end{array}\right)-\sigma\left(\begin{array}[]{c}\mathbf{x}^{k}-\widetilde{\mathbf{x}}^{k}\\ \lambda^{k}-\widetilde{\lambda}^{k}\\ \end{array}\right).

7 end  *

Remark 2.1

If we set $\rho=1$ in (8), then the $\mathbf{x}$-subproblem amounts to estimating the proximity operator of $f$ when $X=\mathcal{R}^{n}$. The implementation of (8) for such cases is thus extremely simple. Here, we allow $\rho\in(-\infty,1]$ just from the theoretical point of view.

Remark 2.2

Note that ${1}/{s}$ in step 5 actually plays a role of penalty parameter in ALM, while $r$ can be treated as the proximal parameter as used in the customized PPA [7]. The quadratic term

[TABLE]

plays a second penalty role for the equality constraint relating to its $k$-th iteration. By the way of updating $\widetilde{\lambda}^{k}$, it uses the convex combination of the feasibility error at the current iteration and the former iteration when $\theta\in[0,1]$. The parameterized matrix designed in this paper is more general than some in the literature:

• •

If $(\theta,\rho)=(0,1)$, then our matrix $G$ given by (4) will become that in **[7, Eq.(2.5)]**. And in such case, the variable $\lambda$ updates practically in the same way as in ALM, but the core $\mathbf{x}$-subproblem in Algorithm 2.1 is a proximal mapping to have a unique minimum since the subproblem is strongly convex.

• •

If $(\theta,\rho)=(2,1)$, then our parameterized proximal matrix turns to the matrix $Q$ involved in **[5, page 158]**. If $(\theta,\rho)=(\tau+1,1)$, then our matrix $G$ is identical to that in **[10, Eq. (3.1)]** but Algorithm 2.1 uses an additional relaxation step for fast convergence. Moreover, we establish the worst-case $\mathcal{O}(1/t)$ ergodic convergence rate for the objective function value error and the feasibility error.

• •

Regardless of step 6, it is easy to check that Algorithm 2.1 with $\theta=\rho=1$ is a linearization of ALM:

[TABLE]

Specifically, by letting $\beta=1/s$ the scheme (10) is ALM with extra proximal term $\frac{1}{2}\|\mathbf{x}-\mathbf{x}^{k}\|^{2}_{rI_{n}-\frac{1}{s}A^{\sf T}A}$ which eliminates the term $\|A\mathbf{x}\|^{2}$ in the iteration. Algorithm 2.1, in such choice of parameters, is a linearized ALM. Besides, our parameter $\theta$ is more general and flexible than that $\theta\in[-1,1]$ in **[16]**.

3 Convergence Analysis

Before analyzing the global convergence and sublinear convergence rate of Algorithm 2.1, we give a fundamental lemma as the following.

Lemma 3.1

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

[TABLE]

where $T_{k+1}$ is given by (16) and

[TABLE]

**Proof ** According to the inequality (2) and the skew-symmetric property of $\mathcal{J}(w)$, i.e.

[TABLE]

the inequality (9) with setting $w=w^{*}$ gives $(w^{*}-\widetilde{w}^{k})^{\sf T}G(\widetilde{w}^{k}-w^{k})\geq 0.$ Note that the step 6 shows

[TABLE]

so we have

[TABLE]

Since the matrix $G$ can be decomposed as $G=\widetilde{G}+\textrm{Diag}\left((1-\rho)A^{\sf T}A/s,\textbf{0}_{m}\right),$ where $\textbf{0}_{m}$ denotes the zero matrix of size $m\times m$ and $\widetilde{G}$ is given by (11), we thus obtain

[TABLE]

Then, applying the identity

[TABLE]

to the left-hand side of (14), the following inequality holds immediately

[TABLE]

where

[TABLE]

Substituting (13) into the expression of $T_{k+1}$, it can be deduced that

[TABLE]

This completes the whole proof. $\ \ \ \blacksquare$

Lemma 3.1 shows the sequence $\{w^{k}\}$ is contractive under the $\widetilde{G}$-norm w.r.t. the solution set $\mathcal{M}^{*}$, since the matrix $\widetilde{G}$ is positive definite and the term $T_{k+1}\geq 0$. Similar to the convergence proof in e.g. [2] and the proof of Lemma 3.1, the global convergence and sublinear convergence rate of Algorithm 2.1 can be easily established as as below, whose proof is omitted here for the sake of conciseness.

Theorem 3.1

Let $(\rho,r,s)$ satisfy (5) and $\{w_{k}\}$ be generated by Algorithm 2.1. Then,

• •

there exists a $w^{\infty}\in\mathcal{M}^{*}$ such that $\lim_{k\to\infty}w^{k}=w^{\infty}$;

• •

for any $t>0$, let $\mathbf{w}_{t}=\frac{1}{t+1}\sum_{k=0}^{t}\widetilde{w}^{k}$ and $\mathbf{x}_{t}=\frac{1}{t+1}\sum_{k=0}^{t}\widetilde{\mathbf{x}}^{k}$. Then,

[TABLE]

Theorem 3.1 illustrates that Algorithm 2.1 converges globally with a sublinear ergodic convergence rate. Furthermore, we can deduce a compact result as the following corollary by making using of the second result in Theorem 3.1. For any $\xi>0$, let $\Gamma_{\xi}=\{\lambda~{}|~{}\xi\geq\|\lambda\|\}$ and

[TABLE]

Corollary 3.1

Let $\{w_{k}\}$ be generated by Algorithm 2.1. For any $\xi>0$, there exists a $\gamma_{\xi}<\infty$ defined in (17) such that for any $t>0$, we have

[TABLE]

**Proof ** By making use of the identity in (12) and by setting $w=(\mathbf{x}^{*},\lambda)\in\mathcal{X}^{*}\times\mathcal{R}^{m}$ into the second result of Theorem 3.1, we have

[TABLE]

where the second equality and the final inequality use $A\mathbf{x}^{*}=b$. Then, it follows from (18) that

[TABLE]

which, by the definition of $\gamma_{\xi}$ in (17), completes the proof. $\ \ \ \blacksquare$

In a similar analysis to (18) together with (2), we can derive $f(\mathbf{x}_{t})-f(\mathbf{x}^{*})-(\lambda^{*})^{\sf T}(A\mathbf{x}_{t}-b)\geq 0$ showing that $f(\mathbf{x}_{t})-f(\mathbf{x}^{*})\geq-\|\lambda^{*}\|\|A\mathbf{x}_{t}-b\|.$ So, taking $\xi=2\|\lambda^{*}\|+1$ in Corollary 3.1, the following inequality

[TABLE]

holds with $\bar{\gamma}=\gamma_{\xi}<\infty$ given by (17). Rearranging the above inequality, we have

[TABLE]

Hence, we will also have $f(\mathbf{x}_{t})-f(\mathbf{x}^{*})\geq-\|\lambda^{*}\|\|A\mathbf{x}_{t}-b\|\geq-\frac{\bar{\gamma}}{2\sigma(t+1)}$ showing that

[TABLE]

According to (20) and (19), both the objective function value error and the feasibility error at the ergodic iterate $\mathbf{x}_{t}$ will decrease in the order of $\mathcal{O}(1/t)$ as $t$ goes to infinity.

4 Numerical Experiments

In this section, we apply the proposed algorithm to solve the following $l_{1}$-minimization problem from signal processing [9], which aims to reconstruct a length $n$ sparse signal from $m(<n)$ observations:

[TABLE]

Note that this is a special case of (1) with specifications $f=\|\mathbf{x}\|_{1}$ and $\mathcal{X}=\mathcal{R}^{n}$. Applying Algorithm 2.1 to problem (21), we derive444The proximity operator is defined as $\textrm{Prox}_{f,r}(\mathbf{x})=\arg\min\left\{f(\mathbf{x})+\frac{r}{2}\|\mathbf{x}-c\|^{2}~{}|c\in\mathcal{R}^{n}\right\}.$ $\mathbf{x}^{k+1}=\textrm{Prox}_{\|\mathbf{x}\|_{1},r}(\mathbf{x})$ that can be explicitly expressed by the shrinkage operator [4] to be coded by the MATLAB inner function ‘withresh’. Followed by Lemma 3.1, we use the following stopping criterions under given tolerance:

[TABLE]

All of the forthcoming experiments use the same starting points $(\mathbf{x}^{0},\lambda^{0})=(0,0)$ and are tested in MATLAB R2018a (64-bit) on Windows 10 system with an Intel Core i7-8700K CPU (3.70 GHz) and 16GB memory.

Consider an original signal $\mathbf{x}\in\mathcal{R}^{10000}$ containing 180 spikes with amplitude $\pm 1$. The measurement matrix $A\in\mathcal{R}^{3000\times 10000}$ is drawn firstly from the standard norm distribution $\mathcal{N}(0,1)$ and then each of its row is normalized. The observation $b$ is generated by $b=A\mathbf{x}+v$, where $v$ is generated by the Gaussian distribution $\mathcal{N}(0,0.01^{2}I)$ on $\mathcal{R}^{3000}$. With the tuned parameters $(r,s,\rho,\sigma)=(8,1.01\|A^{\sf T}A\|_{2}/r,1,1.4)$, some computational results under different parameter $\theta$ are shown in Table 1 in which we present the number of iterations (IT), the CPU time in seconds (CPU), the final obtained residuals It_err and Eq_err, as well as the recovery error $\textrm{RE}={\|\mathbf{x}^{k}-\mathbf{x}_{\textrm{orig}}\|}/{\|\mathbf{x}_{\textrm{orig}}\|}$. Reported results from Table 1 indicate that the choice of $\theta$ could make a great effect on the performance of our algorithm w.r.t. IT and CPU. And it seems that setting $\theta=0.5$ would be a reasonable choice to save the CPU time and to cost fewer number of iterations. The reconstruction results under $\theta=0.5$ are shown in Fig. 1, from which the solution obtained by our algorithm always has the correct number of pieces and is closer to the original noseless signal.

[TABLE]

Table 1: Results by Algorithm 2.1 with different parameter $\theta$.

Next, we consider the problem (21) with large-scale dimensions $(m,n)=(3000,20000)$. By comparing the proposed Algorithms 2.1 (RM-PPA) with the aforementioned tuned parameters to

• •

Algorithm 2.1 without the relaxation step (“M-PPA”),

• •

The customized PPA (“C-PPA”, [7]) with parameters $(\gamma,r,s)=(1.8,8,1.02\|A^{\sf T}A\|_{2}/r)$,

• •

The parameterized PPA (“P-PPA”, [10]) with parameters $(t,r,s)=(-1,8,1.02\|A^{\sf T}A\|_{2}/r)$,

we show comparative results about the convergence behaviors of the residuals $\textrm{LER(k)}:=\log_{2}(\textrm{Eq\_{err}(k)})$ and $\textrm{LIR(k)}:=\log_{2}(\textrm{It\_{err}(k)})$ in Fig. 2, respectively. The effect on recovering the original signal with different algorithms is shown in Fig. 3. Here, we emphasize that the parameter values in [7, 10] can not terminate the algorithms C-PPA and P-PPA because of the fact $\|A^{\sf T}A\|_{2}=1$ for their examples, so we set $r$ the same value as ours but keep $s$ as the value in their experiments. From Figs. 2-3, we observe that M-PPA is competitive to P-PPA and RM-PPA (that is, Algorithm 2.1) performs better than the rest three algorithms.