A general framework for solving convex optimization problems involving the sum of three convex functions

TL;DR

This paper introduces a unified framework for solving convex optimization problems involving three functions, leveraging operator splitting methods, and demonstrates its effectiveness through applications in image processing and signal recovery.

Contribution

It develops a general iterative algorithm framework that unifies existing methods and introduces a new algorithm for complex convex optimization problems involving three functions.

Findings

Existing algorithms are special cases of the proposed framework.

The new algorithm converges under mild conditions.

Numerical experiments show improved efficiency and effectiveness.

Abstract

In this paper, we consider solving a class of convex optimization problem which minimizes the sum of three convex functions $f(x)+g(x)+h(Bx)$, where $f(x)$ is differentiable with a Lipschitz continuous gradient, $g(x)$ and $h(x)$ have a closed-form expression of their proximity operators and $B$ is a bounded linear operator. This type of optimization problem has wide application in signal recovery and image processing. To make full use of the differentiability function in the optimization problem, we take advantage of two operator splitting methods: the forward-backward splitting method and the three operator splitting method. In the iteration scheme derived from the two operator splitting methods, we need to compute the proximity operator of $g+h \circ B$ and $h \circ B$, respectively. Although these proximity operators do not have a closed-form solution in general, they can be solved very efficiently. We mainly employ two different approaches to solve these proximity operators: one is dual and the other is primal-dual. Following this way, we fortunately find that three existing iterative algorithms including Condat and Vu algorithm, primal-dual fixed point (PDFP) algorithm and primal-dual three operator (PD3O) algorithm are a special case of our proposed iterative algorithms. Moreover, we discover a new kind of iterative algorithm to solve the considered optimization problem, which is not covered by the existing ones. Under mild conditions, we prove the convergence of the proposed iterative algorithms. Numerical experiments applied on fused Lasso problem, constrained total variation regularization in computed tomography (CT) image reconstruction and low-rank total variation image super-resolution problem demonstrate the effectiveness and efficiency of the proposed iterative algorithms.

Tables

Table 1. Table 1: Iterative parameters for the proposed algorithms

Methods	Parameters type I	Parameters type II
Algorithm 3.1 and Algorithm 4.1	${\lambda }_1=1.9/{\lambda }_{max}\left(B{B}^T\right)$	${\lambda }_2=1/{\lambda }_{max}\left(B{B}^T\right)$
Algorithm 3.2 and Algorithm 4.2	${\sigma }_1=1/{\Vert B\Vert }^2,{\tau }_1=1$	${\sigma }_2={\tau }_2=1/\Vert B\Vert$

Table 2. Table 2: Numerical results obtained by Algorithm 3.1 , Algorithm 3.2 , Algorithm 4.1 and Algorithm 4.2 in terms of NMSD, SNR(dB) and the iteration numbers (Iter).

Parameters	Methods	$ϵ={10}^{-4}$				$ϵ={10}^{-8}$
type		$NMSD$	$SNR\left(dB\right)$	$Iter$		$NMSD$	$SNR\left(dB\right)$	$Iter$
$I$	Algorithm 3.1	$0.0063$	$44.2488$	$-$		$0.0063$	$44.2488$	$-$
	Algorithm 3.2	$0.0066$	$43.6446$	$757$		$0.0060$	$44.5044$	$986$
	Algorithm 4.1	$0.0067$	$43.6584$	$-$		$0.0067$	$43.6584$	$-$
	Algorithm 4.2	$0.0066$	$43.6411$	$758$		$0.0060$	$44.5044$	$987$
$II$	Algorithm 3.1	$0.0060$	$44.4103$	$387$		$0.0060$	$44.5044$	$626$
	Algorithm 3.2	$0.0076$	$42.3389$	$1119$		$0.0060$	$44.5044$	$1471$
	Algorithm 4.1	$0.0060$	$44.4048$	$387$		$0.0060$	$44.5044$	$627$
	Algorithm 4.2	$0.0076$	$42.3340$	$1120$		$0.0060$	$44.5044$	$1471$

Table 3. Table 3: Numerical results obtained by Algorithm 3.1 , Algorithm 3.2 , Algorithm 4.1 and Algorithm 4.2 with the choice of parameters type I.

Inner iteration	Methods	$ϵ={10}^{-4}$				$ϵ={10}^{-8}$
numbers		$NMSD$	$SNR\left(dB\right)$	$Iter$		$NMSD$	$SNR\left(dB\right)$	$Iter$
$2$	Algorithm 3.1	$0.0062$	$44.2197$	$385$		$0.0060$	$44.5044$	$500$
	Algorithm 3.2	$0.0063$	$44.0661$	$510$		$0.0060$	$44.5044$	$659$
	Algorithm 4.1	$0.0061$	$44.2254$	$385$		$0.0060$	$44.5044$	$500$
	Algorithm 4.2	$0.0063$	$44.0643$	$510$		$0.0060$	$44.5044$	$659$
$10$	Algorithm 3.1	$0.0062$	$44.1763$	$385$		$0.0060$	$44.5044$	$505$
	Algorithm 3.2	$0.0062$	$44.2187$	$386$		$0.0060$	$44.5044$	$505$
	Algorithm 4.1	$0.0062$	$44.1763$	$385$		$0.0060$	$44.5044$	$505$
	Algorithm 4.2	$0.0062$	$44.2176$	$386$		$0.0060$	$44.5044$	$505$
$20$	Algorithm 3.1	$0.0062$	$44.1723$	$385$		$0.0060$	$44.5044$	$506$
	Algorithm 3.2	$0.0062$	$44.1762$	$385$		$0.0060$	$44.5044$	$505$
	Algorithm 4.1	$0.0062$	$44.1720$	$385$		$0.0060$	$44.5044$	$506$
	Algorithm 4.2	$0.0062$	$44.1752$	$385$		$0.0060$	$44.5044$	$505$

Table 4. Table 4: Numerical results obtained by Algorithm 3.1 , Algorithm 3.2 , Algorithm 4.1 and Algorithm 4.2 with the choice of parameters type II.

Inner iteration	Methods	$ϵ={10}^{-4}$				$ϵ={10}^{-8}$
numbers		$NMSD$	$SNR\left(dB\right)$	$Iter$		$NMSD$	$SNR\left(dB\right)$	$Iter$
$2$	Algorithm 3.1	$0.0061$	$44.3020$	$386$		$0.0060$	$44.5044$	$506$
	Algorithm 3.2	$0.0066$	$43.6520$	$683$		$0.0060$	$44.5044$	$895$
	Algorithm 4.1	$0.0061$	$44.2980$	$386$		$0.0060$	$44.5044$	$500$
	Algorithm 4.2	$0.0066$	$43.6464$	$683$		$0.0060$	$44.5044$	$895$
$10$	Algorithm 3.1	$0.0062$	$44.1838$	$385$		$0.0060$	$44.5044$	$505$
	Algorithm 3.2	$0.0062$	$44.1873$	$392$		$0.0060$	$44.5044$	$514$
	Algorithm 4.1	$0.0062$	$44.1826$	$385$		$0.0060$	$44.5044$	$505$
	Algorithm 4.2	$0.0062$	$44.1864$	$392$		$0.0060$	$44.5044$	$514$
$20$	Algorithm 3.1	$0.0062$	$44.1756$	$385$		$0.0060$	$44.5044$	$506$
	Algorithm 3.2	$0.0062$	$44.2183$	$386$		$0.0060$	$44.5044$	$506$
	Algorithm 4.1	$0.0062$	$44.1747$	$385$		$0.0060$	$44.5044$	$506$
	Algorithm 4.2	$0.0062$	$44.2179$	$386$		$0.0060$	$44.5044$	$506$

Table 5. Table 5: Numerical results obtained by Algorithm 3.1 , Algorithm 3.2 , Algorithm 4.1 and Algorithm 4.2 for solving the constrained TV problem in CT image reconstruction.

Inner iteration	Methods	$ϵ={10}^{-4}$				$ϵ={10}^{-6}$				$ϵ={10}^{-8}$
numbers		$NMSD$	$SNR\left(dB\right)$	$Iter$		$NMSD$	$SNR\left(dB\right)$	$Iter$		$NMSD$	$SNR\left(dB\right)$	$Iter$
$1$	Algorithm 3.1	$0.1756$	$15.1103$	$1036$		$0.0213$	$33.4403$	$9524$		$0.0209$	$33.6169$	$14615$
	Algorithm 3.2	$0.2344$	$12.6014$	$1021$		$0.0219$	$33.2001$	$16799$		$0.0209$	$33.6160$	$27983$
	Algorithm 4.1	$0.1756$	$15.1102$	$1036$		$0.0213$	$33.4404$	$9523$		$0.0209$	$33.6169$	$14615$
	Algorithm 4.2	$0.2344$	$12.6014$	$1021$		$0.0219$	$33.2002$	$16799$		$0.0209$	$33.6160$	$27982$
$2$	Algorithm 3.1	$0.1756$	$15.1071$	$1035$		$0.0213$	$33.4399$	$9535$		$0.0209$	$33.6169$	$14635$
	Algorithm 3.2	$0.2035$	$13.8278$	$1004$		$0.0217$	$33.2875$	$11600$		$0.0209$	$33.6166$	$19156$
	Algorithm 4.1	$0.1756$	$15.1071$	$1035$		$0.0213$	$33.4399$	$9534$		$0.0209$	$33.6169$	$14634$
	Algorithm 4.2	$0.2035$	$13.8278$	$1004$		$0.0217$	$33.2874$	$11599$		$0.0209$	$33.6166$	$19156$
$10$	Algorithm 3.1	$0.1757$	$15.1038$	$1034$		$0.0213$	$33.4396$	$9545$		$0.0209$	$33.6169$	$14652$
	Algorithm 3.2	$0.1758$	$15.0997$	$1034$		$0.0213$	$33.4388$	$9548$		$0.0209$	$33.6169$	$14665$
	Algorithm 4.1	$0.1757$	$15.1038$	$1034$		$0.0213$	$33.4397$	$9545$		$0.0209$	$33.6169$	$14652$
	Algorithm 4.2	$0.1758$	$15.0997$	$1034$		$0.0213$	$33.4388$	$9548$		$0.0209$	$33.6169$	$14665$

Table 6. Table 6: Numerical results of the proposed iterative algorithms for solving the low-rank total variation image super-resolution optimization problem ( 1.7 ).

Inner iteration	Methods	$ϵ={10}^{-6}$					$ϵ={10}^{-8}$
numbers		$NMSD$	$SNR\left(dB\right)$	$SSIM$	$Iter$		$NMSD$	$SNR\left(dB\right)$	$SSIM$	$Iter$
$1$	Algorithm 3.1	$0.0781$	$22.1440$	$0.9796$	$34573$		$0.0797$	$21.9701$	$0.9642$	$-$
	Algorithm 3.2	$0.0890$	$21.0146$	$0.9762$	$19046$		$0.0773$	$22.2393$	$0.9635$	$-$
	Algorithm 4.1	$0.0781$	$22.1440$	$0.9796$	$34573$		$0.0797$	$21.9701$	$0.9642$	$-$
	Algorithm 4.2	$0.0890$	$21.0146$	$0.9762$	$19046$		$0.0773$	$22.2393$	$0.9635$	$-$
$2$	Algorithm 3.1	$0.0781$	$22.1440$	$0.9602$	$34573$		$0.0797$	$21.9701$	$0.9642$	$-$
	Algorithm 3.2	$0.0814$	$21.7906$	$0.9549$	$29252$		$0.0783$	$22.1244$	$0.9650$	$-$
	Algorithm 4.1	$0.0781$	$22.1440$	$0.9602$	$34573$		$0.0797$	$21.9701$	$0.9642$	$-$
	Algorithm 4.2	$0.0814$	$21.7906$	$0.9549$	$29252$		$0.0783$	$22.1244$	$0.9650$	$-$
$10$	Algorithm 3.1	$0.0781$	$22.1440$	$0.9602$	$34573$		$0.0797$	$21.9701$	$0.9642$	$-$
	Algorithm 3.2	$0.0781$	$22.1435$	$0.9602$	$34576$		$0.0797$	$21.9707$	$0.9642$	$-$
	Algorithm 4.1	$0.0781$	$22.1440$	$0.9602$	$34573$		$0.0797$	$21.9701$	$0.9642$	$-$
	Algorithm 4.2	$0.0781$	$22.1435$	$0.9602$	$34576$		$0.0797$	$21.9707$	$0.9642$	$-$

Equations

$$\underset{x\in X}{\min }f(x)+g(x),$$

$$\underset{x\in X}{\min }f(x)+h(Bx),$$

$$\underset{x\in X}{\min }f(x)+g(x)+h(Bx),$$

$$\underset{x\in R^n}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2+{\mu }_1\Vert x{\Vert }_1+{\mu }_2\Vert Dx{\Vert }_1,$$

$$D=\left(\begin{array}[]{ccccc}-1&1&&&\\ &-1&1&&\\ &&\cdots&&\\ &&&-1&1\\ \end{array}\right).$$

$$\underset{x\in C}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\mu \Vert x{\Vert }_{TV},$$

$$\underset{x\in R^n}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\mu \Vert x{\Vert }_{TV}+{\delta }_C(x),$$

$$\underset{X\in R^{m\times n}}{\min }\frac{1}{2}\Vert DSX-T{\Vert }_F^2+{\lambda }_1\Vert X{\Vert }_{\ast }+{\lambda }_2\Vert X{\Vert }_{TV},$$

$$f^{\ast }(u)=\underset{x}{\sup }⟨x,u⟩-f(x).$$

$$prox_{\lambda f}(v)=\arg\min_{x}\big{\{}\frac{1}{2}\|x-v\|^{2}+\lambda f(x)\big{\}}.$$

$$\Vert pro{x}_{\lambda f}(x)-pro{x}_{\lambda f}(y){\Vert }^2+\Vert (x-pro{x}_{\lambda f}(x))-(y-pro{x}_{\lambda f}(y)){\Vert }^2\le \Vert x-y{\Vert }^2,\forall x,y\in X.$$

$$\Vert pro{x}_{\lambda f}(x)-pro{x}_{\lambda f}(y){\Vert }^2\le ⟨x-y,pro{x}_{\lambda f}(x)-pro{x}_{\lambda f}(y)⟩,\forall x,y\in X,$$

$${\tilde{f}}_{\lambda }(x)=\underset{y}{\inf }\{f(y)+\frac{1}{2\lambda }\Vert x-y{\Vert }^2\}.$$

$$\nabla {\tilde{f}}_{\lambda }(x)=\frac{1}{\lambda }(x-pro{x}_{\lambda f}(x)),$$

$$pro{x}_{\lambda f}(u)+\lambda pro{x}_{\frac{1}{\lambda }{f}^{\ast }}(\frac{1}{\lambda }u)=u.$$

$$pro{x}_{(\lambda f{)}^{\ast }}(u)=\lambda pro{x}_{\frac{1}{\lambda }{f}^{\ast }}(\frac{1}{\lambda }u).$$

$$(\lambda f{)}^{\ast }(u)$$

$$=\lambda \underset{x}{\sup }⟨x,\frac{1}{\lambda }u⟩-f(x),$$

$$=\lambda f^{\ast }(\frac{1}{\lambda }u).$$

$$pro{x}_{(\lambda f{)}^{\ast }}(u)$$

$$=\arg \underset{x}{\min }\{\frac{1}{2}\Vert x-u{\Vert }^2+\lambda f^{\ast }(\frac{1}{\lambda }x)\}.$$

$$\underset{x}{\min }\frac{1}{2}\Vert x-u{\Vert }^2+\lambda f^{\ast }(\frac{1}{\lambda }x),$$

$$=$$

$$\underset{x\in X}{\min }\frac{1}{2}\Vert x-u{\Vert }^2+g(x)+h(Bx-r),$$

$$\underset{y\in Y}{\max }$$

$$\underset{x}{\min }{f}_1(x)+f_2(Bx).$$

$$x$$

$$y$$

$$x^{k+1}=pro{x}_{\gamma (g+h\circ B)}(x^k-\gamma \nabla f(x^k))+e_k,$$

$$\underset{x\in X}{\min }\left\{\frac{1}{2}\Vert x-(x^k-\gamma \nabla f(x^k)){\Vert }^2+\gamma g(x)+\gamma h(Bx)\right\},$$

$$\underset{y\in Y}{\max }-\frac{1}{2}\Vert B^{\ast }y-u{\Vert }^2+\gamma {\tilde{g}}_{\gamma }(u-B^{\ast }y)-(\gamma h{)}^{\ast }(y)+\frac{1}{2}\Vert u{\Vert }^2$$

$$=\underset{y\in Y}{\max }-\frac{1}{2}\Vert B^{\ast }y-u{\Vert }^2+\gamma {\tilde{g}}_{\gamma }(u-B^{\ast }y)-\gamma h^{\ast }(\frac{1}{\gamma }y)+\frac{1}{2}\Vert u{\Vert }^2.$$

$$\underset{y^{\prime }\in Y}{\max }{\gamma }^2(-\frac{1}{2}\Vert B^{\ast }{y}^{\prime }-\frac{1}{\gamma }u{\Vert }^2+\frac{1}{\gamma }{\tilde{g}}_{\gamma }(u-\gamma B^{\ast }{y}^{\prime })-\frac{1}{\gamma }{h}^{\ast }(y^{\prime })+\frac{1}{2{\gamma }^2}\Vert u{\Vert }^2).$$

Full text

Keywords: Forward-backward splitting method; Three operator splitting method; Dual; Primal-dual; Total variation.

2010 Mathematics Subject Classification: 90C25; 65K10.

1 Introduction

Many problems in signal and image processing, machine learning and statistical prediction can be transformed into the minimization of the sum of two convex functions, a number of efficient iterative algorithms have been proposed in the last ten decades. Here we briefly review some general iterative algorithms. In particular, the forward-backward splitting algorithm [1] is widely used to solve the following optimization problem,

[TABLE]

where $X$ is a Hilbert space, $f\in\Gamma_{0}(X)$ is differentiable with a Lipschitz continuous gradient, $g\in\Gamma_{0}(X)$ is simple and maybe nonsmooth. Here and in what follows, $\Gamma_{0}(X)$ denotes the class of all the proper lower-semicontinuous (lsc) convex functions from a Hilbert space $X$ to $(-\infty,+\infty]$. The iterative shrinkage-thresholding algorithm (ISTA) [2] could be viewed as a special case of the forward-backward splitting algorithm applied to (1.1) when $f(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}$ and $g(x)=\lambda\|x\|_{1}$, where $A\in R^{m\times n}$, $b\in R^{n}$ and $\lambda>0$. As a generalization of the ISTA, Beck and Teboulle [3] proposed a fast iterative shrinkage-thresholding algorithm (FISTA) on the basis of the Nesterov’s accelerated gradient method [4]. The $\ell_{1}-$norm promotes sparsity for one-dimensional signal and plays an important role in the compressive sensing theory [5]. For two dimensional signal, such as image, the total variation (TV) regularization [6] has been widely used for image restoration and image reconstruction since its ability to maintain sharp edges of piecewise constant images(see, for example, [7, 8, 9, 10, 11, 12]). Beck and Teboulle [13] applied the FISTA to solve the total variation image deblurring problem. Since the proximity operator of the total variation has no closed-form solution, they applied the dual gradient method [14] to solve this sub-problem. It is observed that the total variation can be represented by a combination of a convex function with a first-order difference matrix, the following general optimization problem was gained much attention,

[TABLE]

where $X$ and $Y$ are two Hilbert spaces, $f\in\Gamma_{0}(X)$ and $h\in\Gamma_{0}(Y)$, $B:X\rightarrow Y$ is a linear transform. In the literature, the proximity operator of $f$ and $h$ are always assumed to be easily computed. Chambolle and Pock [15] proposed a primal-dual proximity algorithm to solve (1.2) and its’ dual problem in the sense of Fenchel conjugate. The primal-dual proximity algorithm belongs to the primal-dual method, it updates both a primal and dual variable in each iteration. Some other primal-dual methods to solve (1.2) can be found in [16, 17]. It was pointed out in [18] that the primal-dual proximity algorithm can be interpreted as a proximal point algorithm (PPA). Another important method to solve (1.2) is the alternating direction of multiplier method (ADMM) [19]. The ADMM was shown to be equivalent to the Douglas-Rachford splitting method [20] applied to the dual problem of (1.2) in [21]. Compare the primal-dual proximity algorithm with ADMM, the former only needs the matrix-vector multiplications of operator $B$ and its adjoint operator, while ADMM needs to solve a non-trivial sub-problem. If $f$ is differentiable with a Lipschitz continuous gradient, Chen et al. [22] proposed a primal-dual fixed point algorithm based on the proximity operator (PDFP2O) to solve (1.2). When $f(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}$, the PDFP2O reduces to the generalized soft-thresholding algorithm proposed by Loris and Verhoeven [23]. The PDFP2O could be viewed as a combination of the forward-backward splitting algorithm [1] with the fixed point algorithm [24] for computing the proximity operator of $h\circ B$, which is also a primal-dual method. We refer the interested reader to [25] for a comprehensive review of the recent development of the primal-dual related algorithms.

In recent years, the minimization of the sum of three convex functions was received much attention that takes the form of

[TABLE]

where $X$ and $Y$ are two Hilbert spaces, $f\in\Gamma_{0}(X)$ is differentiable with a Lipschitz continuous gradient, $g\in\Gamma_{0}(X)$ and $h\in\Gamma_{0}(Y)$ maybe nonsmooth, $B:X\rightarrow Y$ is a bounded linear operator. In the following, we always assume that the proximal function associated with $g$ and $h$ are easy to be computed. It is obvious that (1.1) and (1.2) are a special case of (1.3). Although the iterative algorithms mentioned before can be exploited to solve the considered three sum of convex functions problem (1.3), the corresponding iterative algorithm would either introduce more dual variables or doesn’t make use of the differentiable property of the function $f(x)$. To overcome these drawbacks, there are several efficient iterative algorithms have been proposed to solve the optimization problem (1.3). Condat [26] proposed a primal-dual splitting algorithm to solve (1.3). Vu [27] proposed an iterative algorithm to solve a zero of a monotone operator inclusion problem, which was based on the primal-dual splitting framework proposed in [28]. It’s observed that Vu’s algorithm coincides with Condat’s algorithm when applied to solve the particular optimization problem (1.3). In the following, we use Condat-Vu algorithm to refer the iterative algorithm proposed by Condat [26] and Vu [27] independently. The Condat and Vu algorithm includes several existing iterative algorithms as its special case. When $f=0$, it becomes Chambolle-Pock algorithm [15, 29], also known as the primal-dual hybrid gradient method (PDHG) [30, 16], and proximal forward-backward splitting algorithm when $h=0$ [1]. It has two iterative parameters restricted by the operator norm $B$ and the Lipschitz constant of $\nabla f$ together. Li and Zhang [31] proved the convergence of a general primal-dual splitting algorithm, which included Condat-Vu algorithm as a special case. They proved the iteration schemes have $O(1/k)$ convergence rate both in the ergodic sense and in the sense of partial primal-dual gap. They also introduced a quasi-Newton and an overrelaxation strategies for accelerating the corresponding algorithms. Wen et al. [32] proposed a self-adaptive primal-dual splitting algorithm, which didn’t require to know the operator norm $B$. Latafat and Patrinos [33] proposed an asymmetric forward-backward-adjoint (AFBA) splitting algorithm, which included the Condat and Vu algorithm. As a generalization of the primal-dual fixed point algorithm based on proximity operator [22, 23] and the preconditioned alternating projection algorithm (PAPA) [34], Chen et al. [35] proposed a so-called primal-dual fixed point (PDFP) algorithm to solve the optimization problem (1.3). They proved the convergence of the PDFP based on the fixed point theory and also obtain the convergence rate of the iteration scheme under suitable conditions. Very recently, Yan [36] proposed a new primal-dual algorithm for solving the optimization problem (1.3), namely primal-dual three operator (PD3O). It reduces to the primal-dual proximity algorithm of Chambolle-Pock algorithm [15, 29] when $f=0$ and PDFP2O [22] when $g=0$. In addition, it recovers the three-operator splitting scheme developed by Davis and Yin [37] when $B$ is the identity operator.

As far as we know, Condat-Vu algorithm [26, 27], PDFP [35] and PD3O [36] are three main type of methods that designed to solve the sum of three convex functions (1.3). However, there are no results to show the relation of these iterative algorithms with other iterative algorithms. In this paper, we will present a unified framework to derive efficient iterative algorithms for solving the optimization problem (1.3), which are based on the forward-backward splitting method [1] and the three operator splitting method [37], respectively. Our main works are divided into two parts: (1) Based on the forward-backward splitting method, we need to compute the proximity operator $g+h\circ B$, which has no closed-form solution in general. To tackle this difficulty, we exploit two approaches: one is dual and the other is primal-dual. In the dual approach, we derive the dual of the minimization problem related to the proximity operator $g+h\circ B$. Thanks to the simple form of function $g(x)$, we show that the corresponding dual problem is the sum of two convex functions with one is differentiable with a Lipschitz continuous gradient and the other is the Fenchel conjugate of $h$. Then, the dual problem is solved once again by the forward-backward splitting algorithm and the primal optimal solution is obtained via the dual optimal solution. In the primal-dual approach, we solve the proximity minimization of $g(x)+h(Bx)$ based on the primal-dual proximity algorithm. We also point out the connection of our proposed iterative algorithms with PDFP [35] and Condat-Vu algorithm [26, 27], respectively. (2) We employ the three operator splitting method to solve the considered optimization problem (1.3). In the three operator iteration scheme, we need to compute the proximity operator of $h\circ B$. At this point, we also use the dual and primal-dual method to solve the proximity operator of $h\circ B$. The connection of the PD3O [36] with our proposed iterative algorithm is presented. We show that PDFP [35] and PD3O [36] are identical. We also obtain a new iterative algorithm to solve the optimization problem (1.3) from the primal-dual approach. To demonstrate the efficiency and effectiveness of the proposed iterative algorithms, we focus on applying them to solve the fused Lasso problem (1.4), the constrained total variation regularization problem (1.5) and the low-rank total variation image super-resolution problem (1.7). We present a detailed discussion on how to select the parameters for the proposed iterative algorithms based on the numerical results.

Fused Lasso problem. The Fused Lasso problem was first introduced by Tibshirani et al. [38] which is an extension of the Lasso problem [39]. It is dedicated to finding a vector with sparsity in both of the coefficients and their successive differences. The optimization problem is

[TABLE]

where $\mu_{1}>0$ and $\mu_{2}>0$ are regularization parameters, $A\in R^{m\times n}$, $b\in R^{n}$ and the difference matrix $D_{(n-1)\times n}$ is defined by

[TABLE]

The fused Lasso problem was applied in many fields, see for example [40, 41, 42, 43]. When $A=I$, the fused Lasso problem reduces to proximity operator of $\mu_{1}\|x\|_{1}+\mu_{2}\|Dx\|_{1}$, which was first studied by Frideman et al. [44]. Liu et al. [41] simplified the method of Frideman et al. [44]. Recently, Shi et al. [45] considered the proximity operator of $g(x)+\mu_{2}\|Dx\|_{1}$, where $g\in\Gamma_{0}(X)$, which generalized the results of Frideman et al. [44] and Liu et al. [41]. Our considered problem (1.3) is more general than the results of [44, 41, 45].

Constrained total variation regularization problem. The total variation (TV) [6] regularization is widely used in the image restoration and image reconstruction problems [46, 47]. To recover an image $x$ is often reduces to solve the following constrained TV regularization problem

[TABLE]

Or equally

[TABLE]

where $A\in R^{m\times n}$ is the system matrix, $b\in R^{m}$ is the collected data which is corrupted by noise, and $\mu>0$ is the regularization parameter which gives the balance between the data error term and the regularization term; $C$ is an nonempty closed convex set which is set according to the prior information on the image, such as nonnegative or bound range. $\delta_{C}$ denotes the indicator function of the closed convex set $C$ such that $\delta_{C}(x)=0$ if $x\in C$ and $+\infty$ otherwise. We note that the total variation $\|x\|_{TV}$ can be represented by a combination of convex function with a discrete gradient operator (See Definition 5.1).

Low-rank total variation image super-resolution problem. Besides the TV regularization, low-rank regularization is widely used in image restoration for recovering missing values in an image. A low-rank total variation (LRTV) image super-resolution optimization problem takes the form of

[TABLE]

where $\lambda_{1}>0$ and $\lambda_{2}>0$ are two regularization parameters, $T$ denotes the observed lower-resolution image, $D$ is a down-sampling operator, $S$ is a blurring operator and $X$ is the high-resolution image that we want to recover. Here, $\|\cdot\|_{F}$ represents the Frobenius-norm of matrix and $\|X\|_{*}$ denotes the nuclear norm of matrix $X$. The LRTV image super-resolution optimization problem (1.7) was first appeared in [48] and further studied in [49]. It is a generalization of image super-resolution by TV-regularization [50].

The rest of the paper is organized as follows. In Section 2, we provide some basic definitions and lemmas from convex analysis. In Section 3, we present the forward-backward splitting method combined with the dual and primal-dual approach to solve the considered optimization problem (1.3). We prove the convergence of the proposed iterative algorithms in finite-dimensional Hilbert spaces. We also point out the relation of some existing iterative algorithms with our iterative algorithms. In Section 4, we use the three operator splitting method together with the dual and primal-dual methods to solve (1.3). We show the convergence of the related iterative algorithms. The relationship between our iterative algorithms and the existing ones will also be included. In Section 5, we apply the proposed iterative algorithms to solve the fused Lasso problem, the constrained total variation regularization problem and the low-rank total variation image super-resolution problem. Finally, we give some conclusions.

2 Preliminaries

Throughout the paper, let $X$ be a Hilbert space equipped with the inner product $\langle\cdot,\cdot\rangle$ and associated norm $\|\cdot\|$. Let $I$ be the identity operator on $X$. The operator norm of $B$ is denoted by $\|B\|=\sup_{x\in X}\frac{\|Bx\|}{\|x\|}$. The adjoint of an operator $B:X\rightarrow Y$ by $B^{*}$ such that $\langle x,By\rangle=\langle B^{*}x,y\rangle$. We will denote $C\subseteq X$ as a nonempty, closed and convex set, the interior of $C$ will be denoted int $C$ and its relative interior ri $C$. Denote $\Gamma_{0}(X)$ as the set of all proper lower-semicontinuous (lsc) convex functions from $X$ to $(-\infty,+\infty]$.

A real-valued function $f:X\rightarrow(-\infty,+\infty]$ is coercive, if $\lim_{\|x\|\rightarrow+\infty}f(x)=+\infty$. The domain of $f$ is defined by dom $f=\{x\in X:f(x)<+\infty\}$ and $f$ is proper if dom $f\neq\emptyset$. We say that a real-valued function $f$ is lower semi-continuous (lsc) if $\lim\inf_{x\rightarrow x_{0}}f(x)\geq f(x_{0})$. The subdifferential of $f$ is the set-valued operator $\partial f:X\rightarrow 2^{X}:x\mapsto\{u\in X|f(y)\geq f(x)+\langle u,y-x\rangle,\forall y\in X\}$.

The following fact will be required.

Lemma 2.1**.**

([51]) Let $g\in\Gamma_{0}(X)$, $h\in\Gamma_{0}(Y)$ and $B:X\rightarrow Y$ is a bounded linear operator such that $0\in\emph{int}(B(\emph{dom}g)-\emph{dom}h)$, then $\partial(g+h\circ B)=\partial g+B^{*}\circ\partial h\circ B$.

We will use the concept of Fenchel Conjugate or Convex Conjugate to derive the dual of the considered optimization problem, which we recall below.

Definition 2.1**.**

(Fenchel Conjugate) Let $f:X\rightarrow(-\infty,\infty]$, the Fenchel Conjugate of $f$ denoted by $f^{*}$ is defined by

[TABLE]

When $f$ is a proper lsc convex function, then $f^{**}=f$. For any $(x,u)\in X\times X$, $u\in\partial f(x)\Leftrightarrow x\in\partial f^{*}(u)$.

1962, Moreau [52] introduced the notion of proximity operator, which plays a very important role in designing proximal algorithms to solve convex optimization problems. A comprehensive review of the proximal algorithms can be found in [53].

Definition 2.2**.**

(Proximity operator) Let $f:X\rightarrow R\cup\{+\infty\}$ be a proper lsc convex function. For $\lambda>0$, the proximity operator $prox_{\lambda f}:X\rightarrow X$ of $\lambda f$ is defined by

[TABLE]

The proximity operator has many important properties. For instance, $prox_{\lambda f}$ is firmly nonexpansive, i.e.,

[TABLE]

Equivalent,

[TABLE]

which means that the proximity operator is nonexpansive due to the Cauchy-Schwarz inequality. The proximity operator is indeed an extension of the orthogonal projection operator. In fact, let $f(x)=\delta_{C}(x)$, then the proximity operator $prox_{\lambda\delta_{C}}(v)=P_{C}(v)$. The proximity operator of $\ell_{1}$-norm is the soft-thresholding operator which is commonly used in the sparse related optimization problems. For other interesting functions, we refer the readers to [54] for the explicit form of proximity operators.

Definition 2.3**.**

(Moreau envelope) Let $f\in\Gamma_{0}(X)$ and $\lambda>0$, the Moreau envelope or Moreau-Yosida regularization is given by

[TABLE]

The next lemma shows that the Moreau-Yosida regularization is differentiable on $X$.

Lemma 2.2**.**

([54]) Let $\lambda>0$ and $f\in\Gamma_{0}(X)$, then the function $f_{\lambda}(X):X\rightarrow R$ is differentiable and its gradient

[TABLE]

with $1/\lambda$-Lipschitz continuous.

The Moreau equality presents the computation of proximity operator a proper lsc convex function $f$ from its convex conjugate $f^{*}$ or converse.

Lemma 2.3**.**

([54]) For any $\lambda>0$ and a vector $u$, the Moreau equality takes the form of

[TABLE]

By the definition of convex conjugate and proximity operator, it is easy to prove the following result.

Lemma 2.4**.**

Let $\lambda>0$ and $f\in\Gamma_{0}(X)$, for any $u\in X$, then

[TABLE]

Proof.

By the definition of convex conjugate, for any $u\in X$, we have

[TABLE]

It follows from the above relation and also together with the definition of proximity operator, we get

[TABLE]

Let $x^{\prime}=\frac{1}{\lambda}x$, then

[TABLE]

Therefore, $prox_{(\lambda f)^{*}}(u)=\lambda prox_{\frac{1}{\lambda}f^{*}}(\frac{1}{\lambda}u)$. ∎

We shall make full use of the following lemmas to derive our iterative algorithm. The lemma was proved in [55]. See also [54].

Lemma 2.5**.**

([55]) Let $g\in\Gamma_{0}(X)$ and $h\in\Gamma_{0}(Y)$. Let $u\in X$ and $r\in Y$. Let $B:X\rightarrow Y$ be a bounded linear operator satisfying the following condition:$r\in\emph{int}(B(\emph{dom}\ g)-\emph{dom}\ h)$. Consider the minimization problem of

[TABLE]

Then the following hold:

(i)* The dual problem of (2.10) is*

[TABLE]

(ii)* Let $y^{*}$ is an optimal solution of the dual problem (2.11) and $x^{*}=prox_{g}(u-B^{*}y^{*})$. Then, $x^{*}$ is an optimal solution of the primal problem (2.10).*

(iii)* The optimal value of the primal problem (2.10) is equal to the optimal value of the dual problem (2.11).*

The following lemma can be found in [31].

Lemma 2.6**.**

([31]) Let $f_{1}\in\Gamma_{0}(X)$ and $f_{2}\in\Gamma_{0}(Y)$. Let $B:X\rightarrow Y$ be a bounded linear operator such that $0\in\emph{int}(B(\emph{dom}\ f_{1})-\emph{dom}\ f_{2})$. Consider the following general optimization problem

[TABLE]

Let $x$ is a solution of (2.12), then for any $\sigma>0,\tau>0$, there exists a vector $y\in Y$ such that

[TABLE]

Conversely, if there exists $\sigma>0,\tau>0$, $x\in X$ and $y\in Y$ satisfying (2.12), then $x$ is a solution of (2.12).

3 A forward-backward splitting method to solve (1.3)

In this section, we will propose efficient iterative algorithms to solve the optimization problem (1.3). We make the following two assumptions throughout the paper: (1) The minimizers of the optimization problem (1.3) is always exists; (2) $0\in\emph{int}(B(\emph{dom}\ g)-\emph{dom}\ h)$. As we have mentioned in the introduction, since the function $f(x)$ in the optimization problem (1.3) is smooth, we can employ the forward-backward splitting algorithm with errors [1] and obtain the iteration scheme as follows. For any $x^{0}\in X$,

[TABLE]

where $\gamma\in(0,2/L)$, $L$ is the Lipschitz constant of $\nabla f$, $e_{k}$ denotes the error between $x^{k+1}$ and the proximity operator $prox_{\gamma(g+h\circ B)}(x^{k}-\gamma\nabla f(x^{k}))$. The following convergence theorem was proved in [1].

Theorem 3.1**.**

([1]) Let $\gamma\in(0,2/L)$, $L$ is the Lipschitz constant of $\nabla f$. For any $x^{0}\in X$, the iterative sequence $\{x^{k}\}$ is defined by (3.1). We assume that $\sum_{n=0}^{\infty}\|e_{k}\|<+\infty$. Then the iterative sequence $\{x^{k}\}$ converges weakly to a solution of the optimization problem (1.3).

Since the proximity operator of the function $\gamma(g+h\circ B)$ has no closed-form solution, so we cannot obtain the exact value of $x^{k+1}$ with the error vector $e_{k}=0$. However, if we can obtain a sufficient approximation of the proximity operator so that the error vector $e_{k}$ satisfies the requirement of Theorem 3.1, the obtained iterative sequences will still converge to the solution of the optimization problem (1.3). To get the updated sequences $\{x^{k+1}\}$, the key problem is to efficient compute the proximity operator of $\gamma(g+h\circ B)$, i.e., $prox_{\gamma(g+h\circ B)}$. In the following, we consider two efficient approaches to solving the proximity operator of function $\gamma(g+h\circ B)$ at the point $x^{k}-\gamma\nabla f(x^{k})$: one is from dual and the other is from primal-dual. Both of the two methods can produce an approximation solution of the proximity operator $prox_{\gamma(g+h\circ B)}(x^{k}-\gamma\nabla f(x^{k}))$ and the error will satisfy the requirement of Theorem 3.1.

3.1 Dual approach

First, we consider the dual approach. Recall that the proximity operator of $prox_{\gamma(g+h\circ B)}(x^{k}-\gamma\nabla f(x^{k}))$ is the minimizer of the following optimization problem

[TABLE]

which is a special case of (2.10). In Lemma 2.5, let $r=0$, $u=x^{k}-\gamma\nabla f(x^{k})$. Define $g:=\gamma g$ and $h:=\gamma h$. Then we obtain the dual formulation of the minimization problem (3.2) is

[TABLE]

The equality in (3.1) is due to the fact that $(\gamma h)^{*}(y)=\gamma h^{*}(\frac{1}{\gamma}y)$. Let $y$ is a solution of the dual problem (3.1), by Lemma 2.5 (ii), we know that $x=prox_{\gamma g}(u-B^{*}y)$ is the optimal solution of (3.2). In (3.1), let $y^{\prime}=\frac{1}{\gamma}y$, then the dual problem (3.1) reduces to

[TABLE]

The corresponding optimal solution $x$ of (3.2) becomes $x=prox_{\gamma g}(u-\gamma B^{*}y^{\prime})$. Let $F(y^{\prime})=\frac{1}{2}\|B^{*}y^{\prime}-\frac{1}{\gamma}u\|^{2}-\frac{1}{\gamma}\widetilde{g}_{\gamma}(u-\gamma B^{*}y^{\prime})$, then the optimal solution of (3.4) is equal to the minimizer of the following minimization problem

[TABLE]

By Lemma 2.2, we have $\nabla F(y^{\prime})=-\frac{1}{\gamma}Bprox_{\gamma g}(u-\gamma B^{*}y^{\prime})$ and $\nabla F(y^{\prime})$ is Lipschitz continuous with Lipschitz constant $\lambda_{max}(BB^{*})$, where $\lambda_{max}(BB^{*})$ denotes the spectral radius of operator $BB^{*}$. In fact, since the proximity operator is nonexpansive, then for any $y_{1},y_{2}\in Y$, we have

[TABLE]

Therefore, we obtain the following iteration scheme to solve (3.5), which is based on the forward-backward splitting algorithm. For any $y^{0}\in Y$, choose $0<\lambda<2/\lambda_{max}(B^{*}B)$,

[TABLE]

In conclusion, we propose the following iterative algorithm to solve the optimization problem (1.3). The iterative algorithm includes an outer iteration step and an inner iteration step, respectively.

We prove the convergence of Algorithm 3.1 in finite-dimensional Hilbert spaces.

Theorem 3.2**.**

Let $\gamma\in(0,2/L)$ and $\lambda\in(0,2/\lambda_{max}(BB^{*}))$. For any $x^{0}\in X$ and $y^{0}\in Y$, the iterative sequence $\{x^{k}\}$ and $\{y^{k}\}$ are generated by Algorithm 3.1. Then the iterative sequence $\{x^{k}\}$ converges to a solution of the optimization problem (1.3).

Proof.

Let $y$ is the optimal solution of the minimization problem (3.5), it follows from Lemma 2.5 that $x=prox_{\gamma g}(u^{k}-\gamma B^{*}y)$ is the optimal solution of (3.2). According to the classical convergence of the forward-backward splitting algorithm, we have $\|y^{k}-y\|\rightarrow 0$ as $k\rightarrow\infty$. Therefore, for a given constant $1/(\gamma\|B\|k^{2})>0$, there exists an integer $j_{k}$, when $J_{k}\geq j_{k}$, we have $\|y^{J_{k}}-y\|\leq 1/(\gamma\|B\|k^{2})$. Notice that $x^{k+1}=prox_{\gamma g}(u^{k}-\gamma B^{*}y^{J_{k}})$, we get

[TABLE]

The first inequality comes from the fact that the proximity operator is nonexpansive and the second inequality is due to the estimation between $y^{J_{k}}$ and $y$. Then $\sum_{n=0}^{\infty}\|e_{k}\|<+\infty$. By Theorem 3.1, we can conclude that the iterative sequence $\{x^{k}\}$ converges to a solution of the optimization problem (1.3). This completes the proof.

∎

3.2 Primal-dual approach

In this part, we present a primal-dual method to solve the minimization problem (3.2). We will see that the primal-dual method can also obtain an accurate approximation to the optimal solution of the minimization problem (3.2). We employ the primal-dual proximity algorithm developed in [15] to solve (3.2) and obtain the following iteration schemes. For any $\overline{x}^{0}\in X$ and $y^{0}\in Y$, choose $\tau>0$ and $\sigma>0$ satisfying $\tau\sigma\|B\|^{2}<1$, let $j_{k}=0,1,2,\cdots$, do

[TABLE]

After simple calculation, we obtain

[TABLE]

Therefore, the detailed iterative algorithm is summarized in Algorithm 3.2.

We prove the convergence of Algorithm 3.2 in finite-dimensional Hilbert spaces.

Theorem 3.3**.**

Let $\gamma\in(0,2/L)$. Let $\sigma>0$ and $\tau>0$ satisfy the condition $\tau\sigma<\frac{1}{\|B\|^{2}}$. For any $x^{0},\overline{x}^{0}\in X$ and $y^{0}\in Y$, the iterative sequences $\{x^{k}\}$, $\{\overline{x}^{k}\}$ and $\{y^{k}\}$ are generated by Algorithm 3.2. Then the iterative sequence $\{x^{k}\}$ converges to a solution of the optimization problem (1.3).

Proof.

In Lemma 2.6, let $f_{1}(x)=\frac{1}{2}\|x-u^{k}\|^{2}+\gamma g(x)$, $f_{2}(x)=\gamma h(x)$. Let $x$ is an optimal solution of (3.2), by Lemma 2.6, for any $\tau>0$ and $\sigma>0$, there exists a vector $y$ such that

[TABLE]

It follows from the definition of the proximity operator, the equation (3.13) reduces to

[TABLE]

Due to Theorem 1 of [15], we know that $\{\overline{x}^{k}\}$ and $\{y^{k}\}$ converge to $x$ and $y$, respectively. Therefore, for a given constant $\frac{1+\tau}{k^{2}}>0$, there exists an integer $j_{k}$, when $J_{k}\geq j_{k}$, we have $\|\overline{x}^{J_{k}}-x\|+\tau\|B\|\|y^{J_{k}}-y\|\leq(1+\tau)\frac{1}{k^{2}}$. Then, we obtain

[TABLE]

Therefore, $\sum_{n=0}^{\infty}\|e^{k}\|<+\infty$. By Theorem 3.1, we can conclude that the iterative sequence $\{x^{k}\}$ converges to a solution of the optimization problem (1.3). This completes the proof.

∎

3.3 Connections to other existing iterative algorithms

In this subsection, we present the connections of the proposed iterative algorithm to some existing iterative algorithms.

In Algorithm 3.1, let $j_{k}=k$ and the number of inner iteration equals to one, then it is reduced to the PDFP [35],

[TABLE]

With the help of Moreau equality (Lemma 2.3), the updated sequence $y^{k+1}$ is equal to

[TABLE]

Let $\overline{y}^{k}=\frac{\gamma}{\lambda}y^{k}$, then the PDFP iteration scheme (3.16) can be rewritten as

[TABLE]

Chen et al. [35] proved the convergence of (3.16) under the conditions that $0<\lambda<1/\lambda_{max}(BB^{*})$ and $0<\gamma<2/L$. Our proposed Algorithm 3.1 provides a more wide selection of the iterative parameter $\lambda$ than the PDFP [35]. It’s observed that the PDFP [35] coincides with the dual forward-backward algorithm proposed in Combettes et al. [55] when $f(x)=\frac{1}{2}\|x-u\|_{2}^{2}$ and the iterative parameter $\gamma=1$.

In the next, we show the connection between Algorithm 3.2 and the Condat-Vu algorithm. In Algorithm 3.2, let $\overline{x}^{0}=x^{0}$, $j_{k}=k$ and fix the number of inner iteration with one, then the iteration scheme of Algorithm 3.2 reduces to

[TABLE]

Let $\overline{y}^{k}=\frac{1}{\gamma}y^{k}$ and after simple calculation, then the iteration scheme (3.19) can be represented as

[TABLE]

Let $\sigma^{\prime}=\frac{\sigma}{\gamma}$ and $\tau^{\prime}=\frac{\tau\gamma}{1+\tau}$, then the above iteration scheme (3.20) recovers the primal-dual splitting algorithm proposed in [26],

[TABLE]

Condat [26] proved the iteration scheme (3.21) converges to a solution of the optimization problem (1.3) under the condition $\frac{1}{\tau^{\prime}}-\sigma^{\prime}\|B\|^{2}>L/2$. This condition is satisfied with the choice of iterative parameters in Algorithm 3.2. In fact, since $\sigma=\gamma\sigma^{\prime}$ and $\tau=\frac{\tau^{\prime}}{\gamma-\tau^{\prime}}$, it follows from the conditions of $\sigma\tau\|B\|^{2}<1$ and $0<\gamma<2/L$, we have

[TABLE]

In the above inequality, all the iterative parameters are mixed. While in our proposed Algorithm 3.2, we have more freedom to choose the iterative parameters than the Condat-Vu algorithm. We note that Yan [36] presented a formulation of Condat-Vu algorithm as follows,

[TABLE]

Yan [36] pointed out the iteration scheme (3.23) may diverges with the condition of $0<\gamma<2/L$ and $0<\sigma<1/\|B\|^{2}$. Under the same condition of the iterative parameters $\gamma$ and $\sigma$, we give another formulation of the Condat-Vu algorithm below, which is converged. In fact, let $\tau=1$ in the iteration scheme (3.20), it reduces to

[TABLE]

The difference between (3.23) and (3.24) is the iterative parameter $\gamma$ involved in the calculation of $x^{k+1}$. Under the choice of $\gamma$ belongs to $(0,2/L)$, the iteration scheme (3.24) maintains converge while the iteration scheme (3.23) may diverges as stated by Yan [36].

4 A three operator splitting method to solve (1.3)

Recently, Yan [36] proposed a primal-dual three operator (PD3O) splitting algorithm to solve the optimization problem (1.3). In particular, when $B=I$, the PD3O [36] coincides with the three operator splitting algorithm proposed by Davis and Yin [37]. The three operator splitting method is a generalization of many well-known operators splitting methods, such as Forward-Backward splitting method [56, 57], Douglas-Rachford splitting method [58, 59], Peaceman-Rachford splitting method [60, 61] and Forward-Douglas-Rachford splitting method [62]. In this section, we will apply the three operator splitting algorithm to solve the optimization problem (1.3). The iteration scheme of the three operator splitting algorithm with errors when applying to the optimization problem (1.3) takes the form of, for any $z^{0}\in X$,

[TABLE]

where $\gamma\in(0,2/L)$, $e_{k}$ is an error vector. We take the following convergence results of the three operator splitting algorithm (4.1) from [37].

Theorem 4.1**.**

([37]) Let $\gamma\in(0,2/L)$, $L$ is the Lipschitz constant of $\nabla f$. For any $x^{0}\in X$, the iterative sequence $\{x^{k}\},\{s^{k}\}$ and $\{z^{k}\}$ are defined by (4.1). We assume that $\sum_{n=0}^{\infty}\|e_{k}\|<+\infty$. Then the iterative sequence $\{x^{k}\}$ and $\{s^{k}\}$ converge weakly to a solution of the optimization problem (1.3).

The key implementation of the above iteration scheme (4.1) is to compute proximity operator of function $\gamma h(Bx)$. Although it still doesn’t have a closed-form solution, as we have done in the last section, we can still obtain a precision solution from two approaches: one is dual and the other is primal-dual. In the following, we will present the details. In particular, we will show that the PD3O [36] is actually a special case of our proposed iterative algorithm.

4.1 Dual approach

In this part, we present how to get the updated iterative sequences $s^{k}$ in (4.1) via the dual approach. For convenience, let $u^{k}=2x^{k}-z^{k}-\gamma\nabla f(x^{k})$, we have

[TABLE]

Let $g:=0$ and $h:=\gamma h$ in Lemma 2.5, it follows from the derivation of (3.1) and (3.4), we obtain the dual problem of the minimization problem (4.2) is

[TABLE]

and the primal optimal solution of (4.2) $v^{*}=u^{k}-\gamma B^{*}y^{*}$, where $y^{*}$ is the dual optimal solution of (4.3). The dual optimization problem (4.3) can be rewritten as follows,

[TABLE]

Both of the optimization problem (4.3) and (4.4) have the same optimal solution $y^{*}$. The corresponding minimization problem (4.4) can be solved by the forward-backward splitting algorithm and the iteration scheme is presented below. For any $y^{0}\in Y$, choose $0<\lambda<2/\lambda_{max}(BB^{*})$,

[TABLE]

Let $s^{k}=u^{k}-\gamma B^{*}y^{J_{k}}$, where $y^{J_{k}}$ is the limit point of the iteration scheme (4.5). For the updated sequence $z^{k+1}$ in (4.1), we have

[TABLE]

In conclusion, we obtain the following iterative algorithm to solve the optimization problem (1.3), which is based on the three operator splitting scheme (4.1).

Similar to Theorem 3.2, we can prove the following convergence theorem of Algorithm 4.1 in finite-dimensional Hilbert spaces.

Theorem 4.2**.**

Let $\gamma\in(0,2/L)$ and $\lambda\in(0,2/\lambda_{max}(BB^{*}))$. For any $x^{0}\in X$ and $y^{0}\in Y$, the iterative sequences $\{x^{k}\}$ and $\{s^{k}\}$ are generated by Algorithm 4.1. Then the iterative sequences $\{x^{k}\}$ and $\{s^{k}\}$ converge to a solution of the optimization problem (1.3).

Proof.

Let $y$ is an optimal solution of (4.3), by Lemma 2.5, $v=u^{k}-\gamma B^{*}y$ is the optimal solution of (4.2), i.e., $prox_{\gamma(h\circ B)}(u^{k})=u^{k}-\gamma B^{*}y$. Then it follows from the classical convergence of the forward-backward splitting algorithm, the iterative sequence $y^{k}\rightarrow y$ as $k\rightarrow\infty$. For a given constant $\frac{1}{\gamma\|B\|k^{2}}>0$, there exists an integer $j_{k}>0$, when $J_{k}\geq j_{k}$, we have $\|y^{J_{k}}-y\|\leq\frac{1}{\gamma\|B\|k^{2}}$. Noticing that $s^{k}=u^{k}-\gamma B^{*}y^{J_{k}}$, we have

[TABLE]

Therefore, $\sum_{k=0}^{\infty}\|e_{k}\|<+\infty$. By Theorem 4.1, we can conclude that the iterative sequences $\{x^{k}\}$ and $\{s^{k}\}$ converge to a solution of the optimization problem (1.3).

∎

4.2 Primal-dual approach

In this part, we employ the primal-dual proximity algorithm to solve the optimization problem (4.2). Given arbitrary $v^{0}\in X$ and $y^{0}\in Y$, for $j_{k}=0,1,2,\cdots$, the iteration scheme is defined by

[TABLE]

where $\tau>0$ and $\sigma>0$ satisfy $\tau\sigma\|B\|^{2}<1$. By the definition of proximity operator and Lemma 2.4, the iteration scheme (4.11) and (4.12) can be simplified as

[TABLE]

In conclusion, we obtain the following iterative algorithm to solve the optimization problem (1.3).

We can also prove the convergence of Algorithm 4.2 in finite-dimensional Hilbert spaces. The proof method is similar to Theorem 3.3.

Theorem 4.3**.**

Let $\gamma\in(0,2/L)$. Let $\sigma>0$ and $\tau>0$ satisfy the condition that $\tau\sigma<\frac{1}{\|B\|^{2}}$. For any $z^{0},v^{0}\in X$ and $y^{0}\in Y$, the iterative sequences $\{x^{k}\}$, $\{z^{k}\}$, $\{v^{k}\}$ and $\{y^{k}\}$ are generated by Algorithm 4.2. Then the iterative sequences $\{x^{k}\}$ and $\{v^{k}\}$ converge to a solution of the optimization problem (1.3).

Proof.

In Lemma 2.6, define $f_{1}:=\frac{1}{2}\|\cdot-u^{k}\|^{2}$ and $f_{2}:=\gamma h$. Let $v$ is the optimal solution of (4.2), then by Lemma 2.6, for any $\sigma>0$ and $\tau>0$, there exists a vector $y$ such that

[TABLE]

It follows from the convergence of the primal-dual proximity algorithm (Theorem 1 of [15]), we have $\{v^{k}\}$ and $\{y^{k}\}$ converge to $v$ and $y$, respectively. Then for a given constant $\frac{1+\tau}{k^{2}}>0$, there exits an integer $j_{k}$, when $J_{k}\geq j_{k}$, we have $\|v^{J_{k}}-v\|+\tau\|B\|\|y^{J_{k}}-y\|\leq\frac{1+\tau}{k^{2}}$. Therefore, we get

[TABLE]

Then $\sum_{n=0}^{\infty}\|e_{k}\|<+\infty$. By Theorem 4.1, we can conclude that the iterative sequences $\{x^{k}\}$ and $\{v^{k}\}$ converge to a solution of the optimization problem (1.3). This completes the proof.

∎

4.3 Connection to existing iterative algorithms

In this part, we show that the proposed Algorithm 4.1 recovers the PD3O [36] proposed by Yan [36]. In fact, let $j_{k}=k$ and set the number of inner iteration equals to one in Algorithm 4.1, then we have the iteration scheme of

[TABLE]

Yan [36] proved the convergence of the PD3O (4.17) under the condition $0<\gamma<2/L$ and $0<\lambda<1/\lambda_{max}(BB^{*})$. Algorithm 4.1 provides a larger range of acceptable parameters $\lambda$ than the PD3O (4.17). Yan [36] showed that the PD3O (4.17) is equivalent to the three operator splitting algorithm [37] for solving the optimization problem (1.3) when $B=I$. In fact, let $\lambda=1$ and $B=I$ in the PD3O (4.17), then it reduces to

[TABLE]

which is exactly the three operator splitting algorithm proposed in [37]. In the next, we show that the PDFP (3.16) is equivalent to the PD3O (4.17). Let $z^{k}=x^{k}-\gamma\nabla f(x^{k})-\gamma B^{*}y^{k}$ and $x^{k}=prox_{\gamma g}(z^{k})$ in the iteration scheme (3.16), then we have

[TABLE]

For the iterative sequence $\{y^{k+1}\}$ in (4.19), we obtain

[TABLE]

Further, let $z^{k+1}=x^{k}-\gamma\nabla f(x^{k})-\gamma B^{*}y^{k+1}$, it follows from $x^{k}=prox_{\gamma g}(z^{k})$ that $x^{k+1}=prox_{\gamma g}(z^{k+1})=prox_{\gamma g}(x^{k}-\gamma\nabla f(x^{k})-\gamma B^{*}y^{k+1})$, which is the same as the updated iterative sequence $\{x^{k+1}\}$ in the iteration scheme of PDFP (3.16). This confirms that the PDFP [35] is equivalent to the PD3O [36]. The PD3O [36] only contains the computation of proximity operator of function $g$ for one time, while the PDFP [35] needs to compute it two times.

Let $j_{k}=k$ and the number of inner iteration equals to one, then the iteration scheme of Algorithm 4.2 is reduced to

[TABLE]

The iteration scheme (4.21) is different from PDFP [35], Condat-Vu algorithm [26, 27] and PD3O [36]. To the best of our knowledge, there is no existing iterative algorithm, which is equivalent to the iterative algorithm (4.21). Besides the before mentioned three iterative algorithms, the iteration scheme (4.21) could be viewed as the fourth type of iterative algorithm to solve the optimization problem (1.3).

5 Numerical results

In this section, we study the performance of the proposed iterative algorithms to solve the fused Lasso problem (1.4), the constrained total variation regularization problem (1.5) and the low-rank total variation image super-resolution problem (1.7). All the experiments are performed in a standard Lenovo Laptop with Intel(R) Core(TM) i7-4712MQ CPU 2.3GHz and 4GB RAM under MATLAB (2013a) software.

(1) Parameters setting. Generally speaking, a large selection of iterative parameter $\gamma$ will reduce to fast convergence of the iterative algorithms. So we fix $\gamma=1.9/L$ for all the proposed iterative algorithms, where $L$ is the Lipschitz constant of $\nabla f$. According to the numerical results of Yan [36], the choice of the iterative parameter $\lambda$ has little influence on the convergence speed. Since our proposed iterative algorithms provide a larger range of acceptable parameters of $\lambda$ than the PDFP [35] and PD3O [36], so we provide two choices of the iterative parameter $\lambda$ in Algorithm 3.1 and Algorithm 4.1. We also give two different choice of $\sigma$ and $\tau$ for Algorithm 3.2 and Algorithm 4.2. The detailed parameters selection are summarized in Table 1.

(2) Performance evaluations. The Signal-to-Noise (SNR) and Normalized Mean Square Distance (NMSD) are used to measure the quality of the reconstructed signal or image, where

[TABLE]

and

[TABLE]

where $\overline{x}$ is the average value of the ideal $x$ and $x_{r}$ is the reconstructed signal or image.

(3) Stopping criteria. We set the relative error between two successive iterative sequences is less than a prescribed tolerance value as the stopping criteria. That is

[TABLE]

where $\epsilon$ is a given small number. If the recontruction $x^{k}$ is an image, then the 2-norm is replaced by Frobenius norm, respectively.

5.1 Fused Lasso problem

To solve the Fused Lasso problem (1.4), let’s define $f(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}$, $g(x)=\mu_{1}\|x\|_{1}$ and $h(Bx)=\mu_{2}\|Dx\|_{1}$ in the optimization problem (1.3), then we can apply the proposed iterative algorithms. We follow Ye and Xiao [42]’s method to generate synthetic dataset. The true coefficient $x=(x_{1},x_{2},\cdots,x_{n})\in R^{n}$ is generated according to

[TABLE]

$A\in R^{m\times n}$ is a random matrix whose elements follow the standard Gaussian distribution and $b=Ax+e$, where $e$ is additive Gaussian noise with mean [math] and variance $0.1$. Here, we set $m=100$ and $n=200$. It’s known that the Lispchitz constant of $\nabla f$ is $\|A\|^{2}$ and the eigenvalues of $DD^{T}$ are $2-2cos(i\pi/n)$ [35], $i=1,2,\cdots,n-1$. So we take $\lambda_{max}(DD^{T})=4$ and $\|D\|=2$. The regularization parameters are set as $\mu_{1}=0.2$ and $\mu_{2}=0.8$.

In the first experiment, we set the number of inner iteration equals to one for all the proposed iterative algorithms, including Algorithm 3.1, Algorithm 3.2, Algorithm 4.1 and Algorithm 4.2. In this case, Algorithm 3.1 and Algorithm 4.1 reduce to the corresponding FDFP [35] and PD3O [36], respectively. Algorithm 3.2 reduces to the Condat-Vu algorithm [26, 27]. We test the performance of these iterative algorithms with given parameters according to Table 1. The numerical results are reported in Table 2. The symbol ${}^{\prime}-^{\prime}$ in Table 2 means it exceeds the maximum $5000$ iteration numbers. We can see from Table 2 that Algorithm 3.2 and Algorithm 4.2 perform nearly the same. Algorithm 3.1 and Algorithm 4.1 take the parameters type I, which do not converge within the required maximum iteration numbers. Under the choice of parameters type II, we find that all these iterative algorithms converge to the same solution when they converge. However, Algorithm 3.2 and Algorithm 4.2 require more iteration numbers than Algorithm 3.1 and Algorithm 4.1. The proposed Algorithm 3.1 and Algorithm 4.1 perform more robustness in the parameters type II than in the parameters type I. For Algorithm 3.2 and Algorithm 4.2, the number of iterations of parameters type I is less than the number of iterations of parameters type II. But this difference is not too much.

Further, we plot the objective function values and SNR values versus the number of iterations in Figure 1 and Figure 2, respectively. Figure 3 shows the recovered signal and the true sparse signal using Algorithm 3.1, Algorithm 3.2, Algorithm 4.1 and Algorithm 4.2 which taking parameters type I.

In the second experiment, we demonstrate how the performance of the proposed iterative algorithms is influenced by the number of inner iterations. The results are reported in Table 3 and Table 4, respectively. We can see from Table 2 that the Algorithm 3.1 and Algorithm 4.1 do not converge within the given maximum iteration numbers when the parameters type I is selected. In Table 3, when the number of inner iterations exceeds one, we find that the Algorithm 3.1 and Algorithm 4.1 converge. It can be seen from Table 3 and Table 4 that by increasing the number of inner iterations, the number of outer iterations required by all the proposed iterative algorithms is reduced, but this decreasing trend will stop over a certain value.

Take inner iteration numbers ten and parameters type I, Figure 4 shows the objective function values and the SNR values versus the number of iterations. We can see from Figure 4 that the convergence of the four iterative algorithms is nearly the same.

5.2 Constrained total variation regularization problem

In this subsection, we employ the proposed iterative algorithms to solve the constrained total variation regularization problem (1.5) arising in computed tomography (CT) image reconstruction. In the following, let’s recall the discrete definition of the total variation.

Definition 5.1**.**

([63, 24]) Let $x\in R^{n}$ be a $\sqrt{n}\times\sqrt{n}$ image. The norm $\|\cdot\|_{2}$ and $\|\cdot\|_{1}$ denote the usual 2-norm and 1-norm of vectors, respectively. First, we define a first order difference matrix $B$ as follows,

[TABLE]

and the matrix $D$ is defined by

[TABLE]

where $I$ denotes the identity matrix and $\otimes$ denotes the Kronecker inner product.

(i)  The isotropic total variation (ITV) is defined by $\|x\|_{ITV}=\|Dx\|_{2,1}$, where $\|y\|_{2,1}=\sum_{i=1}^{n}\sqrt{y_{i}^{2}+y_{n+i}^{2}}$, $y\in R^{2n\times 1}$;

(ii) The anisotropic total variation (ATV) is defined by $\|x\|_{ATV}=\|Dx\|_{1}$.

In this case, let $f(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}$, $g(x)=\delta_{C}(x)$ ($\delta_{C}(x)$ is the indicator function of the closed convex set $C$), and $h(Bx)=\mu\|x\|_{TV}$. The closed convex set $C$ is set as nonnegative set, i.e., $C=\{x\in R^{n}|x_{i}\geq 0\}$. We use the standard Shepp-Logan phantom as the reconstructed image (See Figure 5). The phantom is scanned by fan-beam way with 20 views distributed randomly from 0 to 360 and 320 rays in each view. So the size of the system matrix $A$ is $6400\times 65536$. The simulated projection data is generated by AIRtools [64]. The data vector $b$ is assumed to be corrupted by random Gaussian noise with zero mean and 0.01 variance. For the iterative parameters, we let $\lambda=1/\lambda_{max}(DD^{T})$ for Algorithm 3.1 and Algorithm 4.1, and $\sigma=1/\lambda_{max}(DD^{T}),\tau=1$ for Algorithm 3.2 and Algorithm 4.2. It is known that the discrete gradient matrix $D$ (5.2) has $\lambda_{max}(DD^{T})=8$. The Lipschitz constant of $\nabla f(x)$ is estimated via the power iteration method. The regularization parameter $\mu=0.5$ is used. The numerical results are reported in Table 5.

We can see from Table 5 that the performance of Algorithm 3.1 is the same as Algorithm 4.1. This result confirms that the Algorithm 3.1 and Algorithm 4.1 are equivalent. Similarly, Algorithm 3.2 and Algorithm 4.2 are the same. When the number of inner iterations is one, the primal-dual based Algorithm 3.2 and Algorithm 4.2 require more iterations than the Algorithm 3.1 and Algorithm 4.1, which are proposed by the dual method. By increasing the number of inner iterations, the number of outer iterations required to reach the same accuracy by Algorithm 3.1, Algorithm 3.2, Algorithm 4.1 and Algorithm 4.2 is getting closer. When the number of inner iterations is ten, we can see that the four iterative algorithms converge to the optimal solution at the same rate. The corresponding objective function values are shown in the Figure 6 left, and the SNR values versus the number of iterations is plotted in the Figure 6 right. The reconstructed images of the four iterative algorithms are shown in Figure 7.

5.3 Low-rank total variation image super-resolution problem

High-resolution (HR) images are always needed in various medical imaging diagnoses, such as CT, MRI, and PET et al. However, in practice, due to the limitations of the imaging acquisition systems or to fast the medical imaging reconstruction time, lower-resolution (LR) images are obtained sometimes. There is much demand to generate high-resolution images from lower-resolution images. Super-resolution (SR) image reconstruction is the technique that can achieve this goal. In this subsection, we apply the proposed iterative algorithms to solve the low rank total variation (LRTV) image super-resolution optimization problem (1.7). In the LRTV optimization problem (1.7), there is two regularization terms, one is the nuclear norm $\|X\|_{*}$, which is a convex relaxation of lower rank constraint. The other is the total variation $\|X\|_{TV}$, which can be represented by a combination of a convex function with a linear operator. It is easy to notice that the LRTV optimization problem (1.7) is a special case of the optimization problem (1.3) by letting $f(X)=\frac{1}{2}\|DSX-T\|_{F}^{2}$, $g(X)=\lambda_{1}\|X\|_{*}$ and $h(BX)=\lambda_{2}\|X\|_{TV}$. It is worth mentioning that the proximity operator of nuclear norm has a closed-form solution according to [65].

In the numerical experiment, we use the same data-set as [49]. In detail, a representative 2D slice from T1 MR phantom in Brainweb111http://www.bic.mni.mcgill.ca/brainweb is selected, which has a size of $217\times 181$ with a resolution of 1mm (Figure 8).

The blurring operator $S$ is implemented using a Gaussian Kernel with a standard deviation of 1 pixel. The blurred image is then down-sampled by averaging every 4 pixels. The upsampled operator is implemented by nearest-neighbor interpolation. Besides SNR and NMSD, we also use Structural Similarity Index (SSIM) [66] to evaluate the quality of reconstruction images. The SSIM is defined as,

[TABLE]

where $\mu_{f}$ and $\mu_{g}$ are the mean values respectively in the original HR image $f$ and recovered image $g$, $\sigma_{f}^{2}$ and $\sigma_{g}^{2}$ are the variances, $\sigma_{fg}$ is the covariance of two images, $c_{1}=(k_{1}L)^{2}$ and $c_{2}=(k_{2}L)^{2}$ with $k_{1}=0.01$ and $k_{2}=0.03$, and $L$ is the dynamic range of pixel values. SSIM ranges from 0 to 1, and 1 means perfect recovery.

We test the performance of the proposed iterative algorithms with a choice of regularization parameters $\lambda_{1}=0.01$ and $\lambda_{2}=0.01$ according to [49]. In this numerical tests, we don’t know exactly the operator norm of the down-sampling operator $D$ and the blurring operator $S$, so we tune the parameter $\gamma$ as far as possible to ensure that the proposed iterative algorithms converge. We set $\gamma=0.1$. The other iterative parameters are chosen the same as Subsection 5.2. The numerical results are summarized in Table 6. They symbol ${}^{\prime}-^{\prime}$ means that the maximum iteration number $1\times 10^{5}$ exceeds. We can observe that the proposed iterative algorithms perform the same in terms of NMSD, SNR and SSIM when the stopping criterion $\epsilon=10^{-8}$ and the number of inner iteration number equals to ten. The corresponding SNR and objective function values versus iteration numbers for all the proposed algorithms are plotted in Figure 9. We can see from Figure 9 that the SNR is not monotonically increasing, but first increases and then decreases. For all the proposed algorithms, the maximum SNR is the same and equals to $22.2399(dB)$. The number of iterations used is $48671$, $48718$, $48671$ and $48718$, respectively. This result is consistent with our previous numerical results. In order to have visual inspection, Figure 10 shows the reconstructed images for the proposed iterative algorithms.

6 Conclusions

In this paper, we study a class of convex optimization problem (1.3), which minimizes the sum of three convex functions. We make full use of the gradient of the differentiable convex function of the objective function. To derive effective iterative algorithms for solving the considered optimization problem (1.3), we employ two monotone operator splitting methods: the forward-backward splitting method and the three operator splitting method. In both of the iteration schemes, we are required to compute proximity operator of $g+h\circ B$ and $h\circ B$, respectively. Although these proximity operators have no closed-form solution, it can be solved effectively from dual and primal-dual approach. It is interesting that we find three existing iterative algorithms for solving the optimization problem (1.3) are a special case of our proposed iterative algorithms, respectively. We prove the convergence of the proposed iterative algorithm in finite-dimensional Hilbert spaces. We also provide a decompose way to select iterative parameters of the Condat-Vu algorithm. The equivalence between the PDFP [35] and the PD3O [36] is presented. Numerical experiments on the fused Lasso problem (1.4), the constrained total variation regularization problem (1.5) and the low-rank total variation image super-resolution problem (1.7) demonstrate the effectiveness of our proposed iterative algorithms. For the number of inner iterations, our numerical results confirm it could be fixed with a small number to keep the convergence of the iterative algorithms. In the future work, we will consider to improve the proposed iterative algorithms by using line-search method without the requirement of the Lipschitz constant of $\nabla f$ and also the operator norm of $B$.

Competing interests

The authors declare that they have no competing interests.

ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundations of China (11401293, 11661056, 11771198), the Natural Science Foundations of Jiangxi Province (20151BAB211010), the China Postdoctoral Science Foundation (2015M571989) and the Jiangxi Province Postdoctoral Science Foundation (2015KY51).