Iterative Methods for Computing the Resolvent of the Sum of a Maximal Monotone Operator and Composite Operator with Applications

Abstract

Total variation image denoising models have received considerable attention in the last two decades. To solve constrained total variation image denoising problems, we utilize the computation of a resolvent operator, which consists of a maximal monotone operator and a composite operator. More precisely, the composite operator consists of a maximal monotone operator and a bounded linear operator. Based on recent work, in this paper we propose a fixed-point approach for computing this resolvent operator. Under mild conditions on the iterative parameters, we prove strong convergence of the iterative sequence, which is based on the classical Krasnoselskii–Mann algorithm in general Hilbert spaces. As a direct application, we obtain an effective iterative algorithm for solving the proximity operator of the sum of two convex functions, one of which is the composition of a convex function with a linear transformation. Numerical experiments on image denoising are presented to illustrate the efficiency and effectiveness of the proposed iterative algorithm. In particular, we report the numerical results for the proposed algorithm with different step sizes and relaxation parameters.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In the last two decades, the total variation (TV) image denoising model proposed by Rudin, Osher, and Fatemi [1] has received considerable attention. Often referred to as the ROF model, this takes the form$$\begin{array}{}\text{(1)}& \underset{x\in R^n}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{‖x-u‖}^{\mathrm{2}}+\mu {‖x‖}_{TV},\end{array}$$where $u\in R^n$ is the observed noisy image, ${‖x‖}_{TV}$ is the total variation, and $\mu >\mathrm{0}$ is the regularization parameter, which balances the data fidelity term and total variation regularization. In the following, we denote by $x\in R^n$ the vectorized image $X\in R^{M\times N}$, where $n=MN$. Because TV regularization has the advantage of maintaining image edges when removing noise, it has been extended to many other important image processing problems, including image deblurring [2–4], image inpainting [5], image superresolution [6], image segmentation [7], and image reconstruction [8, 9].

Many efficient iterative algorithms have been proposed to solve the ROF model (1). These include the Chambolle gradient projection algorithm and its variants [10–12], the primal-dual hybrid gradient algorithm [13–15], and the split Bregman algorithm [16, 17]. Isotropic total variation (ITV) and anisotropic total variation (ATV) are the most widely employed methods in the literature. It is worth mentioning that TV includes ITV and ATV, which can both be viewed as compositions of a convex function $\phi$ with a linear transformation $B$. That is, ${‖x‖}_{TV}=\phi (Bx)$. For specific examples, see [18, 19]. Based on this, Micchelli et al. [20] extended the ROF model (1) to the general convex optimization problem$$\begin{array}{}\text{(2)}& \underset{x\in R^n}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{‖x-u‖}^{\mathrm{2}}+\varphi \left(Bx\right),\end{array}$$where $u\in R^n$, $\varphi :R^m\to (-\infty ,+\infty ]$ is a proper lower semicontinuous convex function and $B:R^n\to R^m$ is a linear transformation. The convex optimization problem (2) is equivalent to computing the proximity operator of the function $\varphi \circ B$. Recall that the proximity operator of a proper lower semicontinuous convex function $f:R^n\to (-\infty ,+\infty ]$ is defined as$$\begin{array}{}\text{(3)}& pro{x}_f\left(u\right)=\mathrm{a}\mathrm{r}\mathrm{g}\hspace{0.17em}\underset{x\in R^n}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\frac{\mathrm{1}}{\mathrm{2}}{‖x-u‖}^{\mathrm{2}}+f\left(x\right)\right\}.\end{array}$$Thus, Micchelli et al. [20] proved that the above minimization problem (2) can be solved using a fixed-point equation. The advantage of the fixed-point framework is that it provides a convenient analysis of the convergence of the proposed algorithm and enables the development of efficient numerical algorithms via various fixed-point iterations.

Because the pixel values of grayscale images are generally distributed in $[\mathrm{0,255}]$ or $[\mathrm{0,1}]$, it is natural to incorporate this prior information into the ROF model (1). Thus, we obtain the following constrained ROF model:$$\begin{array}{}\text{(4)}& \underset{x\in C}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{‖x-u‖}^{\mathrm{2}}+\mu {‖x‖}_{TV},\end{array}$$where $C$ is a nonempty closed and convex set, which could be chosen as the pixel range of the images. To solve the constrained ROF model (4), it is sufficient to consider the following constrained convex optimization problem:$$\begin{array}{}\text{(5)}& \underset{x\in C}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{‖x-u‖}^{\mathrm{2}}+\varphi \left(Bx\right).\end{array}$$In general, the constrained convex optimization problem (5) can always be reformulated as an unconstrained convex optimization problem. More precisely,$$\begin{array}{}\text{(6)}& \underset{x\in R^n}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{‖x-u‖}^{\mathrm{2}}+\varphi \left(Bx\right)+{\delta }_C\left(x\right),\end{array}$$where ${\delta }_C$ denotes the indicator function of $C$, which is defined by ${\delta }_C(x)=\mathrm{0}$ if $x\in C$ and ${\delta }_C(x)=+\infty$ otherwise. As the fixed-point algorithm proposed by Micchelli et al. [20] cannot be applied to solve the problem in (6), this motivated us to develop an efficient iterative algorithm for this purpose. Hence, we apply this algorithm to the constrained total variation denoising model (4).

As the indicator function ${\delta }_C$ belongs to the class of proper lower semicontinuous convex functions, we are motivated to solve the following general convex optimization problem in general Hilbert spaces:$$\begin{array}{}\text{(7)}& \underset{x\in H}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{‖x-u‖}^{\mathrm{2}}+\varphi \left(Bx\right)+g\left(x\right).\end{array}$$where $g$ is a proper lower semicontinuous convex function. To achieve this goal, let us recall the definition of the resolvent operator. Let $H$ be a real Hilbert space and let $T:H\to {\mathrm{2}}^H$ be a maximal monotone operator. The resolvent of $T$ is the single-valued operator $J_T=(I+T{)}^{-\mathrm{1}}$. Next, let $f:H\to (-\infty ,+\infty ]$ be a proper lower semicontinuous convex function and let $\partial f$ denote the subdifferential of $f$, and let $T=\partial f$. Thus, the resolvent of $T$ coincides with the proximity operator $pro{x}_f$ as follows:$$\begin{array}{}\text{(8)}& J_{\partial f}\left(u\right)=pro{x}_f\left(u\right).\end{array}$$Under certain qualification conditions, the problem considered in (7) can be solved via the resolvent operator $J_{B^{\ast }\circ \partial \varphi \circ B+\partial g}$. Overall, this leads to the problem of computing the resolvent operator of $S+B^{\ast }TB$. More precisely, let $u\in H$ and$$\begin{array}{}\text{(9)}& x≔J_{S+B^{\ast }TB}\left(u\right),\end{array}$$where $S:H\to {\mathrm{2}}^H$ and $T:G\to {\mathrm{2}}^G$ are two maximal monotone operators and $B:H\to G$ is a bounded linear operator from the Hilbert space $H$ to the Hilbert space $G$. Following this, let $S=\partial g$ and $T=\partial \varphi$. Hence, the problem (9) reduces to the convex optimization problem (7), because $J_{\partial g+B^{\ast }\circ \partial \varphi \circ B}(u)=pro{x}_{g+\varphi \circ B}(u)$.

1.1. Existing Work

Next, let us briefly review some existing work concerning the computation of resolvent operators. Bauschke and Combettes [21] extended the Dykstra algorithm [22] for computing projections onto the intersection of two closed convex sets to compute the resolvent of the sum of two maximal monotone operators. Hence, they obtained an algorithm for finding the proximity operator of the sum of two proper lower semicontinuous convex functions. Combettes [23] proposed two inexact parallel splitting algorithms for computing the resolvent of a weighted sum of maximal monotone operators. The key idea was to reformulate the weighted sum of maximal monotone operators as a sum of two maximal monotone operators in a product space. The two iterative algorithms were based on extensions of the Douglas-Rachford splitting and Dykstra-like methods, respectively. Furthermore, Combettes [23] applied these algorithms when computing the proximity operator of a weighted sum of proper lower semicontinuous convex functions. In more recent work, Artacho and Campoy [24] generalized the averaged alternating modified reflection algorithm [25] to compute the resolvent of the sum of two maximal monotone operators.

In contrast, Moudafi [26] proposed a fixed-point algorithm to compute the resolvent of operator $B^{\ast }TB$, where $B:H\to G$ is a bounded linear operator with the adjoint $B^{\ast }$, $T:G\to {\mathrm{2}}^G$ is a maximal monotone operator, and $H$ and $G$ are two Hilbert spaces. When $T=\partial \varphi$, this algorithm utilizes the fixed-point approach to computing the proximity operator $pro{x}_{\varphi \circ B}$ proposed by Micchelli et al. [20]. It is clear that the resolvent of the operator $B^{\ast }TB$ is a special case of (9) by allowing $S=\mathrm{0}$. In addition, Combettes et al. [27] proposed a dual forward-backward splitting algorithm to solve the convex optimization problem in (7) (see also [28]). The main idea in this case was to first derive the dual problem of (7) and then apply the forward-backward splitting algorithm to solve this problem. Finally, the convergence of the primal iterative sequence was proved using the connection between the primal optimal solution and the dual optimal solution. However, the authors did not apply the obtained iterative algorithm to solve the constrained ROF model (4). Beck and Teboulle [3] solved the constrained ROF model (4) as a subproblem of image deblurring. Chan et al. [29] applied the alternating direction method of multipliers to solve the constrained total variation image deblurring problem. Their numerical results confirmed that the quality of the restored image could be improved by incorporating prior information of the images as a constraint. However, neither study considered the general resolvent operator of (9).

In this study, we focus on computing the resolvent of the operator $S+B^{\ast }TB$ (9), where $S:H\to {\mathrm{2}}^H$ and $T:G\to {\mathrm{2}}^G$ are two maximal monotone operators and $B:H\to G$ is a bounded linear operator. We assume that the operator $S+B^{\ast }TB$ is maximally monotone. This is true, for example, if $\mathrm{0}\in ri(B\mathrm{\hspace{0.17em}\hspace{0.17em}}(\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{\hspace{0.17em}\hspace{0.17em}}S)-\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{\hspace{0.17em}\hspace{0.17em}}T)$ [30], where $ri(C)$ represents a relative interior of the set $C$. If $S=\mathrm{0}$, then the problem (9) becomes that of computing the resolvent operator of $B^{\ast }TB$. Inspired by the work of Moudafi [26] and Combettes et al. [27], we propose a fixed-point approach to computing resolvent operators (9). First, we show that the resolvent operator (9) can be characterized using a fixed-point equation. Subsequently, we propose a fixed-point algorithm and prove the strong convergence of its iterative sequence. Next, we employ the proposed iterative algorithm to solve the convex optimization problem (7) arising in the field of image denoising. Finally, we conduct numerical experiments on image denoising to validate the effectiveness of the proposed algorithm. In particular, we show how the performance of the algorithm is influenced by the selection of the step size and relaxation parameters.

The remainder of this paper is organized as follows. In Section 2, we introduce some notation and present useful definitions and lemmas. In Section 3, we present the main fixed-point algorithm and prove its strong convergence. In Section 4, we employ the obtained iterative algorithm to solve a particular convex optimization problem, which is related to the calculation of the resolvent operator (9). In Section 5, we present some numerical experiments on image denoising to illustrate the performance of our proposed algorithm. Finally, we provide some conclusions and ideas for future work in Section 6.

2. Preliminaries

In this section, we review some basic definitions and lemmas in monotone operator theory and convex analysis, which will be used throughout this paper. First, let $H$ be a real Hilbert space with inner product $⟨·,·⟩$ and the associated norm $‖·‖$. Let $I$ denote the identity operator on $H$. We use $x^k\rightharpoonup x$ to indicate that the sequence $\{x^k\}$ converges weakly to $x$ and $x^k\to x$ to indicate that the sequence $\{x^k\}$ converges strongly to $x$. Let $L(H,G)$ denote all bounded linear operators from the Hilbert space $H$ to the Hilbert space $G$. Let $B\in L(H,G)$. Then, the adjoint of $B$ is the unique operator $B^{\ast }\in L(G,H)$ such that $⟨Bx,y⟩=⟨x,B^{\ast }y⟩$.

Let $A:H\to {\mathrm{2}}^H$ be a set-valued operator. The sets $\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{\hspace{0.17em}\hspace{0.17em}}A=\{x\in H\mid Ax\ne \varnothing \}$, $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{\hspace{0.17em}\hspace{0.17em}}A=\{u\in H\mid \exists x\in H,\hspace{0.17em}\hspace{0.17em}u\in Ax\}$, $gra\mathrm{\hspace{0.17em}\hspace{0.17em}}A=\{(x,u)\in H\times H\mid u\in Ax\}$, and $zer\mathrm{\hspace{0.17em}\hspace{0.17em}}A=\{x\in H\mid \mathrm{0}\in Ax\}$ represent the domain, range, graph, and zero set of $A$, respectively. The inverse $A^{-\mathrm{1}}$ of $A$ is the set-valued operator with the graph $\{(u,x)\in H\times H\mid u\in Ax\}$.

We now introduce some definitions and lemmas, most of which can be found in [28, 31].

Definition 1 (see [28], (maximal monotone operator)).

Let $A:H\to {\mathrm{2}}^H$ be a set-valued operator. $A$ is said to be monotone if$$\begin{array}{}\text{(10)}& ⟨x-y,u-v⟩\ge \mathrm{0},\hspace{1em}\forall \left(x,u\right)\in gra\hspace{0.17em}A,\hspace{0.17em}\hspace{0.17em}\left(y,v\right)\in gra\hspace{0.17em}A.\end{array}$$Moreover, $A$ is maximally monotone if its graph is not strictly contained in the graph of any other monotone operator.

Definition 2 (see [28], (resolvent and Yosida approximation)).

Let $A:H\to {\mathrm{2}}^H$ and $\gamma >\mathrm{0}$. The resolvent of $\gamma A$ is $J_{\gamma A}=(I+\gamma A{)}^{-\mathrm{1}}$, and the Yosida approximation of $\gamma A$ is ${\mathrm{\hspace{0.17em}\hspace{0.17em}}}^{\gamma }A=\left(\mathrm{1}/\gamma \right)(I-J_{\gamma A})$.

For any $x\in H$, it is easy to derive that ${\mathrm{\hspace{0.17em}\hspace{0.17em}}}^{\gamma }Ax\in A(J_{\gamma A}x)$.

Lemma 3 (see [28]).

Let $A:H\to {\mathrm{2}}^H$ be maximally monotone and $\gamma >\mathrm{0}$. Then, the following hold:

(i)

$J_{\gamma A}:H\to H$, and $I-{\gamma }_{\gamma A}:H\to H$ are firmly nonexpansive and maximally monotone.

(ii)

The Yosida approximation ${\mathrm{\hspace{0.17em}\hspace{0.17em}}}^{\gamma }A$ is $\gamma -$cocoercive and maximally monotone.

 

Definition 4 (see [28]).

Let $C$ be a nonempty subset of $H$ and $T:C\to C$. Then

(i)

$T$ is nonexpansive if $‖Tx-Ty‖\le ‖x-y‖\hspace{0.17em}\hspace{0.17em}\forall x,y\in C$.

(ii)

$T$ is firmly nonexpansive if ${‖Tx-Ty‖}^{\mathrm{2}}\le {‖x-y‖}^{\mathrm{2}}-{‖(I-T)x-(I-T)y‖}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\forall x,y\in C$.

(iii)

$T$ is $\alpha -$averaged if there exists a nonexpansive $S$ and $\alpha \in (\mathrm{0,1})$ such that $T=(\mathrm{1}-\alpha )I+\alpha S$, where $I$ denotes the identity operator.

(iv)

$T$ is $\beta -$cocoercive (or $\beta -$inversely strong monotone) for $\beta >\mathrm{0}$ if $⟨x-y,Tx-Ty⟩\ge \beta {‖Tx-Ty‖}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\forall x,y\in C$.

 

Remark 5.

(i) An equivalent definition of firm nonexpansiveness is that ${‖Tx-Ty‖}^{\mathrm{2}}\le ⟨x-y,Tx-Ty⟩\hspace{0.17em}\hspace{0.17em}\forall x,y\in C$. It is easy to check that firm nonexpansiveness implies nonexpansiveness. If $T$ is firmly nonexpansive, then $T$ is $\mathrm{1}-$cocoercive.

(ii) An equivalent definition of $\alpha -$averaged is ${‖Tx-Ty‖}^{\mathrm{2}}\le {‖x-y‖}^{\mathrm{2}}-\left((\mathrm{1}-\alpha )/\alpha \right){‖(I-T)x-(I-T)y‖}^{\mathrm{2}}$ for all $x,y\in C$. Thus, if $T$ is $\alpha -$ averaged then $T$ is also nonexpansive.

(iii) Let $T$ be $\beta -$cocoercive. Then, for any $\gamma >\mathrm{0}\hspace{0.17em}\hspace{0.17em}\gamma T$ is $\beta /\gamma -$cocoercive. Furthermore, owing to the Cauchy-Schwartz inequality $T$ is $\mathrm{1}/\beta -$ Lipschitz continuous; i.e., $‖Tx-Ty‖\le \left(\mathrm{1}/\beta \right)‖x-y‖$.

The following lemma provides some useful characterizations between an operator $T$ and $I-T$.

Lemma 6 (see [28]).

Let $C$ be a nonempty subset of $H$ and let $T:C\to C$. Then

(i)

$T$ is nonexpansive if and only if $I-T$ is $\mathrm{1}/\mathrm{2}-$cocoercive.

(ii)

$T$ is firmly nonexpansive if and only if $I-T$ is firmly nonexpansive.

(iii)

$T$ is $\mathrm{1}/\mathrm{2}v-$averaged if and only if $I-T$ is $v-$cocoercive, where $v>\mathrm{1}/\mathrm{2}$.

 

The following lemma shows that a composition of two averaged operators is also an average. This result first appeared in the work of Ogura and Yamada [32]. Combettes and Yamada [33] subsequently confirmed it with a different proof.

Lemma 7 (see [32]).

Let $C$ be a nonempty subset of $H$. Furthermore, let $T_{\mathrm{1}}:C\to C$ be ${\alpha }_{\mathrm{1}}-$ averaged and let $T_{\mathrm{2}}:C\to C$ be ${\alpha }_{\mathrm{2}}-$ averaged. Then $T_{\mathrm{1}}{T}_{\mathrm{2}}$ is $\left({\alpha }_{\mathrm{1}}+{\alpha }_{\mathrm{2}}-\mathrm{2}{\alpha }_{\mathrm{1}}{\alpha }_{\mathrm{2}}\right)/\left(\mathrm{1}-{\alpha }_{\mathrm{1}}{\alpha }_{\mathrm{2}}\right)$ averaged.

We also make full use of the following lemma.

Lemma 8 (see [28]).

Let $x\in H$ and $y\in H$. Let $\alpha \in R$ and $R$ denote the set of real numbers. Then$$\begin{array}{}\text{(11)}& {‖\alpha x+\left(\mathrm{1}-\alpha \right)y‖}^{\mathrm{2}}+\alpha \left(\mathrm{1}-\alpha \right){‖x-y‖}^{\mathrm{2}}=\alpha {‖x‖}^{\mathrm{2}}+\left(\mathrm{1}-\alpha \right){‖y‖}^{\mathrm{2}}.\end{array}$$

We end this section by introducing the Krasnoselskii–Mann algorithm. Theorem 9 provides a fundamental tool for studying the convergence of many operator splitting methods.

Theorem 9 (see [28], (Krasnoselskii–Mann algorithm)).

Let $C$ be a nonempty closed convex subset of $H$, and let $T:C\to C$ be a nonexpansive operator such that $Fix(T)\ne \varnothing$, where $Fix(T)$ denotes the fixed-point set of $T$. Let ${\lambda }_k\in (\mathrm{0,1})$ be such that ${\sum }_{k=\mathrm{0}}^{\infty }‍{\lambda }_k(\mathrm{1}-{\lambda }_k)=+\infty$. For any $x_{\mathrm{0}}\in C$, define$$\begin{array}{}\text{(12)}& x_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)x_k+{\lambda }_{k}T{x}_k.\end{array}$$Then, the following hold:

(i)

$\{x_k\}$ is Fejér monotone with respect to $Fix(T)$; i.e., $‖x_{k+\mathrm{1}}-p‖\le ‖x_k-p‖$ for any $p\in Fix(T)$.

(ii)

$‖x_k-T{x}_k‖$ converges strongly to $\mathrm{0}$.

(iii)

$\{x_k\}$ converges weakly to a fixed point in $Fix(T)$.

 

3. Computing the Resolvent Operator (9)

Before presenting our main results, we first introduce some notation. For a fixed $u\in H$, let $\gamma >\mathrm{0}$ define operators $F:G\to G$ and $Q:G\to G$ as follows:$$\begin{array}{}\text{(13)}& F\left(y\right)=y+B{J}_S\left(u-\gamma B^{\ast }y\right),\hspace{1em}\forall y\in G,\end{array}$$and$$\begin{array}{}\text{(14)}& Q\left(y\right)=\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)F\left(y\right),\hspace{1em}\forall y\in G.\end{array}$$

In addition, the following lemma provides a fixed-point characterization of the resolvent operator (9).

Lemma 10.

Let $S:H\to {\mathrm{2}}^H$ and $T:G\to {\mathrm{2}}^G$ be two maximal monotone operators. Furthermore, let $B\in L(H,G)$ and $u\in H$. Then, for a given $\gamma >\mathrm{0}$,$$\begin{array}{}\text{(15)}& J_{S+B^{\ast }TB}\left(u\right)=J_S\left(u-\gamma B^{\ast }y\right)\end{array}$$if and only if $y$ is a fixed point of $Q$.

Proof.

It follows from the definition of the resolvent operator that $J_{S+B^{\ast }TB}(u)=(I+S+B^{\ast }TB{)}^{-\mathrm{1}}(u)$. Thus, we can write$$\begin{array}{}\text{(16)}& J_{S+B^{\ast }TB}\left(u\right)\in u-\left(S+B^{\ast }TB\right)J_{S+B^{\ast }TB}\left(u\right).\end{array}$$Let $y$ satisfy (15). Then we have $u-\gamma B^{\ast }y\in (I+S)J_{S+B^{\ast }TB}(u)$; that is $J_{S+B^{\ast }TB}(u)\in u-(\gamma B^{\ast }y+S{J}_{S+B^{\ast }TB}(u))$. By (16), we have $y\in \left(\mathrm{1}/\gamma \right)TB{J}_{S+B^{\ast }TB}(u)$. Thus, $y\in \left(\mathrm{1}/\gamma \right)TB{J}_S(u-\gamma B^{\ast }y)$. Using the fact that $y\in \left(\mathrm{1}/\gamma \right)Tx\iff x=J_{\left(\mathrm{1}/\gamma \right)T}(x+y)$, we deduce that$$\begin{array}{}\text{(17)}& B{J}_S\left(u-\gamma B^{\ast }y\right)=J_{\left(\mathrm{1}/\gamma \right)T}\left(B{J}_S\left(u-\gamma B^{\ast }y\right)+y\right),\end{array}$$which implies that $y$ is a fixed point of $Q$.

Conversely, let $y$ be a fixed point of $Q$. Then, we have that$$\begin{array}{c}\text{(18)}& y=\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(y+B{J}_S\left(u-\gamma B^{\ast }y\right)\right)\iff B{J}_S\left(u-\gamma B^{\ast }y\right)=J_{\left(\mathrm{1}/\gamma \right)T}\left(y+B{J}_S\left(u-\gamma B^{\ast }y\right)\right)\iff \\ y+B{J}_S\left(u-\gamma B^{\ast }y\right)\in \left(I+\frac{\mathrm{1}}{\gamma }T\right)\left(B{J}_S\left(u-\gamma B^{\ast }y\right)\right)\iff \\ \gamma y\in TB{J}_S\left(u-\gamma B^{\ast }y\right)\Rightarrow \\ u\in u-\gamma B^{\ast }y+B^{\ast }TB{J}_S\left(u-\gamma B^{\ast }y\right)\iff \\ u\in \left(I+S\right)J_S\left(u-\gamma B^{\ast }y\right)+B^{\ast }TB{J}_S\left(u-\gamma B^{\ast }y\right)\iff \\ J_{S+B^{\ast }TB}\left(u\right)=J_S\left(u-\gamma B^{\ast }y\right).\end{array}$$This completes the proof.

Next, we prove the following lemma, which characterizes an important property of the operator $F$.

Lemma 11.

Let $S:H\to {\mathrm{2}}^H$ be a maximal monotone operator, and let $B\in L(H,G)$. For a given $u\in H$ and $\gamma >\mathrm{0}$, define an operator $R:G\to G:y\mapsto -B{J}_S(u-\gamma B^{\ast }y)$. Then

(i)

$R$ is $\mathrm{1}/\gamma {‖B‖}^{\mathrm{2}}-$cocoercive.

(ii)

For any $\gamma \in (\mathrm{0},\mathrm{2}/{‖B‖}^{\mathrm{2}})$, $F=I-R$ is $\gamma {‖B‖}^{\mathrm{2}}/\mathrm{2}-$averaged.

 

Proof.

(i) Let $y_{\mathrm{1}},y_{\mathrm{2}}\in G$. Then, we have that$$\begin{array}{}\text{(19)}& ⟨y_{\mathrm{1}}-y_{\mathrm{2}},R\left(y_{\mathrm{1}}\right)-R\left(y_{\mathrm{2}}\right)⟩=⟨y_{\mathrm{1}}-y_{\mathrm{2}},-B{J}_S\left(u-\gamma B^{\ast }{y}_{\mathrm{1}}\right)+B{J}_S\left(u-\gamma B^{\ast }{y}_{\mathrm{2}}\right)⟩=\frac{\mathrm{1}}{\gamma }⟨\left(u-\gamma B^{\ast }{y}_{\mathrm{1}}\right)-\left(u-\gamma B^{\ast }{y}_{\mathrm{2}}\right),J_S\left(u-\gamma B^{\ast }{y}_{\mathrm{1}}\right)-J_S\left(u-\gamma B^{\ast }{y}_{\mathrm{2}}\right)⟩\ge \frac{\mathrm{1}}{\gamma }{‖J_S\left(u-\gamma B^{\ast }{y}_{\mathrm{1}}\right)-J_S\left(u-\gamma B^{\ast }{y}_{\mathrm{2}}\right)‖}^{\mathrm{2}}\ge \frac{\mathrm{1}}{\gamma {‖B‖}^{\mathrm{2}}}{‖\left(-B{J}_S\left(u-\gamma B^{\ast }{y}_{\mathrm{1}}\right)\right)-\left(-B{J}_S\left(u-\gamma B^{\ast }{y}_{\mathrm{2}}\right)\right)‖}^{\mathrm{2}}=\frac{\mathrm{1}}{\gamma {‖B‖}^{\mathrm{2}}}{‖R\left(y_{\mathrm{1}}\right)-R\left(y_{\mathrm{2}}\right)‖}^{\mathrm{2}}.\end{array}$$The first inequality follows from the fact that the resolvent operator is firmly nonexpansive, and the second is trivial. Thus, $R$ is $\mathrm{1}/\gamma {‖B‖}^{\mathrm{2}}-$cocoercive.

(ii) For any $\gamma \in (\mathrm{0},\mathrm{2}/{‖B‖}^{\mathrm{2}})$, we know from Lemma 6(iii) that $F=I-R$ is $\gamma {‖B‖}^{\mathrm{2}}/\mathrm{2}-$averaged.

Lemma 11 shows that, for any $\gamma \in (\mathrm{0},\mathrm{2}/{‖B‖}^{\mathrm{2}})$, $F$ is $\gamma {‖B‖}^{\mathrm{2}}/\mathrm{2}-$averaged. On the other hand, $I-J_{\left(\mathrm{1}/\gamma \right)T}$ is firmly nonexpansive and is also $\mathrm{1}/\mathrm{2}-$averaged. According to Lemma 7, $Q=(I-J_{\left(\mathrm{1}/\gamma \right)T})\circ F$ is $\mathrm{2}/(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})-$averaged.

Now, we are ready to present our main results.

Theorem 12.

Let $S:H\to {\mathrm{2}}^H$ and $T:G\to {\mathrm{2}}^G$ be two maximal monotone operators. Furthermore, let $B\in L(H,G)$ and $u\in H$. For any $y_{\mathrm{0}}\in G$, we define the iterative sequences $\{x_k\}$ and $\{y_k\}$ as follows:$$\begin{array}{c}\text{(20)}& x_k=J_S\left(u-\gamma B^{\ast }{y}_k\right),y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_k\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(y_k+B{x}_k\right),\end{array}$$where $\gamma \in (\mathrm{0,2}/{‖B‖}^{\mathrm{2}})$ and ${\lambda }_k\in (\mathrm{0},(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})/\mathrm{2})$ such that ${\sum }_{k=\mathrm{0}}^{\infty }‍{\lambda }_k(\mathrm{1}-\left(\mathrm{2}/(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})\right){\lambda }_k)=+\infty$. Then, the following hold:

(i)

$\{y_k\}$ converges weakly to a fixed point of $Q$.

(ii)

Suppose that $\inf {\lambda }_k>\mathrm{0}$. Then, $\{x_k\}$ converges strongly to the resolvent operator (9).

 

Proof.

(i) Because the resolvent operator $J_{S+B^{\ast }TB}(u)$ exists and is unique, the first part of the Lemma 10 implies that $Q$ has a fixed point. Let $Q(y_k)=(I-J_{\left(\mathrm{1}/\gamma \right)T})(y_k+B{J}_S(u-\gamma B^{\ast }{y}_k))$. Then, the iterative sequence $\{y_{k+\mathrm{1}}\}$ can then be rewritten as follows:$$\begin{array}{}\text{(21)}& y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_{k}Q{y}_k.\end{array}$$Because $Q$ is $\mathrm{2}/(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})-$averaged, there exists a nonexpansive mapping $U$ such that $Q=(\mathrm{1}-\mathrm{2}/(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}}))I+\left(\mathrm{2}/(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})\right)U$ and $Fix(U)=Fix(Q)$. Therefore, the iterative sequence (21) is equivalent to$$\begin{array}{}\text{(22)}& y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_k\left(\left(\mathrm{1}-\frac{\mathrm{2}}{\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}}}\right)y_k+\frac{\mathrm{2}}{\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}}}U{y}_k\right)=\left(\mathrm{1}-\frac{\mathrm{2}}{\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}}}{\lambda }_k\right)y_k+\frac{\mathrm{2}}{\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}}}{\lambda }_{k}U{y}_k.\end{array}$$Notice that $\left(\mathrm{2}/(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})\right){\lambda }_k\in (\mathrm{0,1})$ and ${\sum }_{k=\mathrm{0}}^{\infty }‍{\lambda }_k(\mathrm{1}-\left(\mathrm{2}/(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})\right){\lambda }_k)=+\infty$. By Theorem 9, we can conclude the following:

(1)

$\{y_k\}$ is Fejér monotone with respect to $Fix(U)$.

(2)

$‖y_k-U{y}_k‖\to \mathrm{0}$ as $k\to \infty$.

(3)

$\{y_k\}$ converges weakly to a fixed point of $U$.

Taking into account the fact that $Fix(U)=Fix(Q)$, $\{y_k\}$ also converges weakly to a fixed point of $Q$, and $\{y_k\}$ is Fejér monotone with respect to $Fix(Q)$. Therefore, ${\mathrm{l}\mathrm{i}\mathrm{m}}_{k\to \infty }‖y_k-y‖$ exists for any $y\in Fix(Q)$. Further, we have that $‖y_k-Q{y}_k‖=‖\left(\mathrm{2}/(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})\right)(y_k-U{y}_k)‖\to \mathrm{0}$, as $k\to \infty$.

(ii) Let $y\in Fix(Q)$. By Lemma 11 (ii) and Remark 5 (ii), we then have that$$\begin{array}{}\text{(23)}& {‖Q{y}_k-y‖}^{\mathrm{2}}={‖\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(y_k+B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)\right)-\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(y+B{J}_S\left(u-\gamma B^{\ast }y\right)\right)‖}^{\mathrm{2}}\le {‖\left(y_k+B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)\right)-\left(y+B{J}_S\left(u-\gamma B^{\ast }y\right)\right)‖}^{\mathrm{2}}\le {‖y_k-y‖}^{\mathrm{2}}-\frac{\mathrm{1}-\alpha }{\alpha }{‖-B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)-\left(-B{J}_S\left(u-\gamma B^{\ast }y\right)\right)‖}^{\mathrm{2}},\end{array}$$where $\alpha =\gamma {‖B‖}^{\mathrm{2}}/\mathrm{2}$. From Lemma 8, we have that$$\begin{array}{}\text{(24)}& {‖y_{k+\mathrm{1}}-y‖}^{\mathrm{2}}={‖\left(\mathrm{1}-{\lambda }_k\right)\left(y_k-y\right)+{\lambda }_k\left(Q{y}_k-y\right)‖}^{\mathrm{2}}=\left(\mathrm{1}-{\lambda }_k\right){‖y_k-y‖}^{\mathrm{2}}+{\lambda }_k{‖Q{y}_k-y‖}^{\mathrm{2}}-{\lambda }_k\left(\mathrm{1}-{\lambda }_k\right){‖y_k-Q{y}_k‖}^{\mathrm{2}}.\end{array}$$Substituting (23) into (24), we derive$$\begin{array}{}\text{(25)}& {‖y_{k+\mathrm{1}}-y‖}^{\mathrm{2}}\le {‖y_k-y‖}^{\mathrm{2}}-\frac{\mathrm{1}-\alpha }{\alpha }{\lambda }_k{‖B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)-B{J}_S\left(u-\gamma B^{\ast }y\right)‖}^{\mathrm{2}}-{\lambda }_k\left(\mathrm{1}-{\lambda }_k\right){‖y_k-Q{y}_k‖}^{\mathrm{2}},\end{array}$$which implies that$$\begin{array}{}\text{(26)}& \frac{\mathrm{1}-\alpha }{\alpha }{\lambda }_k{‖B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)-B{J}_S\left(u-\gamma B^{\ast }y\right)‖}^{\mathrm{2}}\le {‖y_k-y‖}^{\mathrm{2}}-{‖y_{k+\mathrm{1}}-y‖}^{\mathrm{2}}+{\lambda }_k\left({\lambda }_k-\mathrm{1}\right){‖y_k-Q{y}_k‖}^{\mathrm{2}}.\end{array}$$Because $\inf {\lambda }_k>\mathrm{0}$, ${\mathrm{l}\mathrm{i}\mathrm{m}}_{k\to \infty }‖y_k-y‖$ exists, and ${\mathrm{l}\mathrm{i}\mathrm{m}}_{k\to \infty }‖y_k-Q{y}_k‖=\mathrm{0}$. Taking $k\to \infty$ in the above inequality (26), we obtain$$\begin{array}{}\text{(27)}& \underset{k\to \infty }{\mathrm{l}\mathrm{i}\mathrm{m}}‖B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)-B{J}_S\left(u-\gamma B^{\ast }y\right)‖=\mathrm{0}.\end{array}$$

By Lemma 10, $J_{S+B^{\ast }TB}(u)=J_S(u-\gamma B^{\ast }y)$, where $y\in Fix(Q)$. Then, we have that$$\begin{array}{}\text{(28)}& {‖x_k-J_S\left(u-\gamma B^{\ast }y\right)‖}^{\mathrm{2}}={‖J_S\left(u-\gamma B^{\ast }{y}_k\right)-J_S\left(u-\gamma B^{\ast }y\right)‖}^{\mathrm{2}}\le ⟨\left(u-\gamma B^{\ast }{y}_k\right)-\left(u-\gamma B^{\ast }y\right),J_S\left(u-\gamma B^{\ast }{y}_k\right)-J_S\left(u-\gamma B^{\ast }y\right)⟩=\gamma ⟨B^{\ast }y-B^{\ast }{y}_k,J_S\left(u-\gamma B^{\ast }{y}_k\right)-J_S\left(u-\gamma B^{\ast }y\right)⟩=\gamma ⟨y-y_k,B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)-B{J}_S\left(u-\gamma B^{\ast }y\right)⟩\le \gamma ‖y-y_k‖‖B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)-B{J}_S\left(u-\gamma B^{\ast }y\right)‖.\end{array}$$The first inequality follows from the fact that the resolvent operator is firmly nonexpansive, and the final inequality is derived from the Cauchy-Schwarz inequality. Taking $k\to \infty$ on the right of (28), we have that $‖x_k-J_S(u-\gamma B^{\ast }y)‖\to \mathrm{0}$. Thus, we can conclude that the iterative sequence $\{x_k\}$ converges strongly to the resolvent operator of $J_{S+B^{\ast }TB}(u)$. This completes the proof.

Remark 13.

We observe that $(\mathrm{4}-\gamma {‖B‖}^{\mathrm{2}})/\mathrm{2}>\mathrm{1}$ for any $\gamma \in (\mathrm{0},\mathrm{2}/{‖B‖}^{\mathrm{2}})$. Taking ${\lambda }_k=\mathrm{1}$ in (20), we obtain a simple iterative algorithm to compute the resolvent operator problem (9). More precisely, for any $y_{\mathrm{0}}\in G$ the iterative sequences $\{x_k\}$ and $\{y_k\}$ are generated as follows:$$\begin{array}{c}\text{(29)}& x_k=J_S\left(u-\gamma B^{\ast }{y}_k\right),y_{k+\mathrm{1}}=\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(y_k+B{x}_k\right),\end{array}$$where $\gamma \in (\mathrm{0},\mathrm{2}/{‖B‖}^{\mathrm{2}})$. The iterative algorithm (29) is equivalent to the Picard iteration scheme, which can be easily applied in practice. In fact, the iterative scheme (29) can be rewritten as$$\begin{array}{}\text{(30)}& y_{k+\mathrm{1}}=\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(y_k+B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)\right).\end{array}$$With the help of the relation $J_{\gamma T}+\gamma J_{{\gamma }^{-\mathrm{1}}{T}^{-\mathrm{1}}}\circ {\gamma }^{-\mathrm{1}}I=I$, we obtain from (30) that$$\begin{array}{}\text{(31)}& y_{k+\mathrm{1}}=\frac{\mathrm{1}}{\gamma }{J}_{\gamma T^{-\mathrm{1}}}\left(\gamma \left(y_k+B{J}_S\left(u-\gamma B^{\ast }{y}_k\right)\right)\right).\end{array}$$Let ${\overline{y}}_k=\gamma y_k$. Then, the iterative scheme (31) is equivalent to$$\begin{array}{}\text{(32)}& {\overline{y}}_{k+\mathrm{1}}=J_{\gamma T^{-\mathrm{1}}}\left({\overline{y}}_k+\gamma B{J}_S\left(u-B^{\ast }{\overline{y}}_k\right)\right).\end{array}$$

Remark 14.

We observe that the resolvent operator of $S+B^{\ast }TB$ in (9) is equivalent to the following monotone inclusion problem:$$\begin{array}{c}\text{(33)}& \mathrm{f}\mathrm{i}\mathrm{n}\mathrm{d}\hspace{0.17em}\hspace{0.17em}x\in H,\mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\in x-u+Sx+B^{\ast }TBx,\end{array}$$where $u$, $S$, $T$, and $B$ are the same as in (9). This monotone inclusion (33) can be solved using existing methods for more general monotone inclusion problems (such as [34–37]). However, these iterative algorithms are not identical to our proposed algorithm. To illustrate this point, we will explain the proposed algorithm from the perspective of duality. In fact, the dual monotone inclusion of (33) is$$\begin{array}{c}\text{(34)}& \mathrm{f}\mathrm{i}\mathrm{n}\mathrm{d}\hspace{0.17em}\hspace{0.17em}y\in G,\mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{0}\in T^{-\mathrm{1}}y-B{J}_S\left(u-B^{\ast }y\right).\end{array}$$By Lemma 11, we know that $-B{J}_S(u-B^{\ast }y)$ is $\mathrm{1}/{‖B‖}^{\mathrm{2}}$-cocoercive. Thus, $y$ is a solution of the dual monotone inclusion (34), and $x=J_S(u-B^{\ast }y)$ is a solution of the primal monotone inclusion (33). Next, we apply the relaxed forward-backward splitting algorithm (see, e.g., Theorem 26.14 of [28]) to solve the dual monotone inclusion (34). More precisely, let $y_{\mathrm{0}}\in G$ and set$$\begin{array}{}\text{(35)}& y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_k{J}_{\gamma T^{-\mathrm{1}}}\left(y_k+\gamma B{J}_S\left(u-B^{\ast }{y}_k\right)\right),\end{array}$$which is equivalent to the iterative algorithm introduced in Theorem 12. For convenience, we summarize the differences between the proposed algorithm and existing algorithms in Table 1.

Table 1

Comparison of the proposed algorithm with existing algorithms.

Method	Iterative algorithm	Operator splitting algorithms
Primal-dual	[34, 35]	Forward-backward-forward splitting
	[36]	Variable metric forward-backward splitting
	[37]	Forward-backward-half forward splitting
Dual	Proposed (20)	Relaxed forward-backward splitting

Let $S=\mathrm{0}$ in Theorem 12. Then, we obtain the following corollary.

Corollary 15.

Let $T:G\to {\mathrm{2}}^G$ be a maximal monotone operator and $B\in L(H,G)$, and let $u\in H$. For any $y_{\mathrm{0}}\in G$, we define the iterative sequences $\{x_k\}$ and $\{y_k\}$ as follows:$$\begin{array}{c}\text{(36)}& x_k=u-\gamma B^{\ast }{y}_k,y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_k\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(y_k+B{x}_k\right),\end{array}$$where $\gamma$ and ${\lambda }_k$ are the same as in Theorem 12. Then, the following hold:

(i)

$\{y_k\}$ converges weakly to a fixed point of $Q$, where $Q(·)=(I-J_{\left(\mathrm{1}/\gamma \right)T})(Bu+(I-\gamma B{B}^{\ast })(·))$.

(ii)

If $\mathrm{i}\mathrm{n}\mathrm{f}{\lambda }_k>\mathrm{0}$, then $\{x_k\}$ converges strongly to the resolvent of $J_{B^{\ast }TB}(u)$ and $J_{B^{\ast }TB}(u)=u-\gamma B^{\ast }y$, $y\in Fix(Q)$.

 

Remark 16.

Moudafi [26] proposed the following iterative algorithm to solve the resolvent operator $J_{B^{\ast }TB}(u)$. Let $y_{\mathrm{0}}\in G$. Then, the iterative sequence $\{y_k\}$ is defined by$$\begin{array}{}\text{(37)}& y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_k\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(Bu+\left(I-\gamma B{B}^{\ast }{y}_k\right)\right),\end{array}$$where $\gamma \in (\mathrm{0},\mathrm{2}/{‖B‖}^{\mathrm{2}})$ and ${\lambda }_k\in (\mathrm{0,1})$ such that ${\sum }_{k=\mathrm{0}}^{+\infty }‍{\lambda }_k(\mathrm{1}-{\lambda }_k)=+\infty$. Moudafi proved that $J_{B^{\ast }TB}(u)=u-\gamma B^{\ast }{y}^{\infty }$, where $y^{\infty }$ is the weak limit point of $\{y_k\}$. Furthermore, he presented several applications of the obtained iterative algorithm.

Corollary 15 extends some of the results from [26] in two aspects. (i) The range of $\{{\lambda }_k\}$ is expanded and (ii) we obtain strong convergence of the iterative sequence $\{x_k\}$ (which is not studied in [26]).

Let $B=I$ in Theorem 12. Then, we obtain the following corollary.

Corollary 17.

Let $S:H\to {\mathrm{2}}^H$ and $T:G\to {\mathrm{2}}^G$ be two maximal monotone operators. In addition, let $u\in H$. For any $y_{\mathrm{0}}\in G$, we define the iterative sequences $\{x_k\}$ and $\{y_k\}$ as follows:$$\begin{array}{c}\text{(38)}& x_k=J_S\left(u-\gamma y_k\right),y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_k\left(I-J_{\left(\mathrm{1}/\gamma \right)T}\right)\left(y_k+x_k\right),\end{array}$$where $\gamma \in (\mathrm{0,2})$ and ${\lambda }_k\in (\mathrm{0},(\mathrm{4}-\gamma )/\mathrm{2})$ such that ${\sum }_{k=\mathrm{0}}^{\infty }‍{\lambda }_k(\mathrm{1}-\left(\mathrm{2}/(\mathrm{4}-\gamma )\right){\lambda }_k)=+\infty$. Then, the following hold:

(i)

$\{y_k\}$ converges weakly to a fixed point of $Q$, where $Q(·)=(I-J_{\left(\mathrm{1}/\gamma \right)T})(·+J_S(u-\gamma ·))$.

(ii)

Suppose that $\mathrm{i}\mathrm{n}\mathrm{f}{\lambda }_k>\mathrm{0}$. Then, $\{x_k\}$ converges strongly to the resolvent operator $J_{S+T}(u)$.

Remark 18.

The obtained iterative algorithm (38) for computing the resolvent operator $J_{S+T}(u)$ is different from those proposed by Bauschke and Combettes [21, 23]. In fact, similar to (32), the iterative scheme (38) can be rewritten as$$\begin{array}{c}\text{(39)}& x_k=J_S\left(u-{\overline{y}}_k\right),{\overline{y}}_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right){\overline{y}}_k+{\lambda }_k{J}_{\gamma T^{-\mathrm{1}}}\left({\overline{y}}_k+\gamma x_k\right),\end{array}$$where ${\overline{y}}_k=\gamma y_k$.

Bauschke and Combettes [21] proposed a Dykstra-like algorithm to compute the resolvent operator $J_{S+T}(u)$. Let $x_{\mathrm{0}}=u$, $p_{\mathrm{0}}=\mathrm{0}$, and $q_{\mathrm{0}}=\mathrm{0}$. Then, the Dykstra-like algorithm is defined by$$\begin{array}{c}\text{(40)}& y_k=J_T\left(x_k+p_k\right),p_{k+\mathrm{1}}=x_k+p_k-y_k,\\ x_{k+\mathrm{1}}=J_S\left(y_k+q_k\right),\\ q_{k+\mathrm{1}}=y_k+q_k-x_{k+\mathrm{1}}.\end{array}$$In addition, Bauschke and Combettes proved that $\{x_k\}$ and $\{y_k\}$ generated by (40) converge strongly to $J_{S+T}(u)$. Following some simple calculation, the Dykstra-like algorithm (40) can be rewritten as follows:$$\begin{array}{c}\text{(41)}& x_{k+\mathrm{1}}=J_S\left(J_T\left(x_k+p_k\right)+q_k\right),p_{k+\mathrm{1}}=J_{T^{-\mathrm{1}}}\left(x_k+p_k\right),\\ q_{k+\mathrm{1}}=J_{S^{-\mathrm{1}}}\left(J_T\left(x_k+p_k\right)+q_k\right).\end{array}$$

On the other hand, Combettes [23] proposed an inexact Douglas–Rachford splitting algorithm and an inexact Dykstra-like algorithm for computing the resolvent of the sum of a finite family of maximal monotone operators. For the resolvent of the sum of two maximal monotone operators, the inexact Dykstra-like algorithm without errors coincided with the iterative algorithm (40). For simplicity, we have presented the inexact Douglas-Rachford splitting algorithm without errors for computing the resolvent of the sum of two maximal monotone operators. Let $z_{\mathrm{0}}\in H$, and define the iterative sequences as follows:$$\begin{array}{c}\text{(42)}& y_k=J_{\left(\mathrm{2}\gamma /\left(\mathrm{2}+\gamma \right)\right)S}\left(\frac{\mathrm{2}}{\mathrm{2}+\gamma }{z}_k+\frac{\gamma }{\mathrm{2}+\gamma }u\right),p_k=J_{\left(\mathrm{2}\gamma /\left(\mathrm{2}+\gamma \right)\right)T}\left(\frac{\mathrm{2}}{\mathrm{2}+\gamma }\left(\mathrm{2}{y}_k-z_k\right)+\frac{\gamma }{\mathrm{2}+\gamma }u\right),\\ z_{k+\mathrm{1}}=z_k+{\lambda }_k\left(p_k-y_k\right),\end{array}$$where $\gamma >\mathrm{0}$, and $\{{\lambda }_k\}\subseteq (\mathrm{0,2})$ such that ${\sum }_{k=\mathrm{0}}^{+\infty }‍{\lambda }_k(\mathrm{2}-{\lambda }_k)=+\infty$ and $\inf {\lambda }_k>\mathrm{0}$. Next, the sequences $\{y_k\}$ and $\{p_k\}$ converge strongly to $J_{S+T}(u)$. In fact, this iterative algorithm (42) is equivalent to applying the Douglas–Rachford splitting algorithm to the monotone inclusion of$$\begin{array}{}\text{(43)}& \mathrm{0}\in \left(\frac{\mathrm{1}}{\mathrm{2}}x-\frac{\mathrm{1}}{\mathrm{2}}u+Sx\right)+\left(\frac{\mathrm{1}}{\mathrm{2}}x-\frac{\mathrm{1}}{\mathrm{2}}u+Tx\right).\end{array}$$

Comparing (39), (40), and (42), we find that the obtained iterative sequences generated by all algorithms converge strongly to the resolvent operator $J_{S+T}(u)$. The iterative algorithms (39) and (42) do not have any requirements for the initial value, whereas (40) requires the selection of a fixed initial value. Unlike (39) and (42), the Dykstra-like algorithm (40) does not require tuning of the parameters.

4. Application to Convex Optimization Problem

In this section, we apply the obtained results to solve a particular convex optimization problem that has been studied in the literature.

For convenience, we introduce some notation. A function $f:H\to (-\infty ,+\infty ]$ is called proper if $\mathrm{d}\mathrm{o}\mathrm{m}\hspace{0.17em}f=\{x\in H\mid f(x)<+\infty \}\ne \varnothing$. ${\Gamma }_{\mathrm{0}}(H)$ denotes the class of proper lower semicontinuous convex functions from $H$ to $(-\infty ,+\infty ]$. The Fenchel-conjugate of $f$ is defined by $f^{\ast }:u\mapsto {\mathrm{s}\mathrm{u}\mathrm{p}}_{x\in H}⟨x,u⟩-f(x)$.

Problem 19.

Let $r\in G$ and $u\in H$, and let $f:H\to (-\infty ,+\infty ]$ and $g:G\to (-\infty ,+\infty ]$ be proper lower semicontinuous convex functions. Furthermore, let $L:H\to G$ be a nonzero bounded linear operator such that $r\in \mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{\hspace{0.17em}\hspace{0.17em}}(L(\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{\hspace{0.17em}\hspace{0.17em}}f)-\mathrm{d}\mathrm{o}\mathrm{m}\mathrm{\hspace{0.17em}\hspace{0.17em}}g)$, where the core of a subset $C\subseteq H$ is defined by $\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{\hspace{0.17em}\hspace{0.17em}}C=\{x\in C\mid \mathrm{c}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{\hspace{0.17em}\hspace{0.17em}}(C-x)=H\}$ (see, e.g., Definition 6.9 of [28]). Consider the following convex optimization problem:$$\begin{array}{}\text{(44)}& \underset{x\in H}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{‖x-u‖}^{\mathrm{2}}+f\left(x\right)+g\left(Lx-r\right).\end{array}$$The minimization problem (44) is equivalent to the proximity operator of $f(x)+g(Lx-r)$. As a direct application of the resolvent operator (9), we obtain the following convergence theorem from Theorem 12.

Theorem 20.

Under the conditions of Problem 19, let $y_{\mathrm{0}}\in G$ and set $\{x_k\}$ and $\{y_k\}$ as$$\begin{array}{c}\text{(45)}& x_k=pro{x}_f\left(u-\gamma L^{\ast }{y}_k\right),y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_k\left(I-pro{x}_{\left(\mathrm{1}/\gamma \right)g}\right)\left(y_k+L{x}_k-r\right),\end{array}$$where $\gamma \in (\mathrm{0},\mathrm{2}/{‖L‖}^{\mathrm{2}})$, and ${\lambda }_k\in (\mathrm{0},(\mathrm{4}-\gamma {‖L‖}^{\mathrm{2}})/\mathrm{2})$ such that ${\sum }_{k=\mathrm{0}}^{\infty }‍{\lambda }_k(\mathrm{1}-\left(\mathrm{2}/(\mathrm{4}-\gamma {‖L‖}^{\mathrm{2}}\mathrm{\hspace{0.17em}\hspace{0.17em}})\right){\lambda }_k)=+\infty$ and $\inf {\lambda }_k>\mathrm{0}$. Then, the iterative sequence $\{x_k\}$ converges strongly to a solution of (44).

Proof.

Let $S=\partial f$, $T=\partial g(·-r)$, and $B=L$ in (9). Because $f\in {\Gamma }_{\mathrm{0}}(H)$ and $g\in {\Gamma }_{\mathrm{0}}(G)$, $S$ and $T$ are maximal monotone operators. Furthermore, note that $J_S=(I+\partial f{)}^{-\mathrm{1}}=pro{x}_f$ and $J_{\left(\mathrm{1}/\gamma \right)T}=(I+\left(\mathrm{1}/\gamma \right)\partial g(·-r){)}^{-\mathrm{1}}=pro{x}_{\left(\mathrm{1}/\gamma \right)g}(·-r)+r$. Hence, the iterative sequences (45) follow directly from (20). By Theorem 12, we conclude that $\{x_k\}$ converges strongly to a solution of (44). This completes the proof.

Remark 21.

The Moreau identity states that, for any $\lambda >\mathrm{0}$, $pro{x}_{\lambda f}(u)+\lambda pro{x}_{\left(\mathrm{1}/\lambda \right)f^{\ast }}(\left(\mathrm{1}/\lambda \right)u)=u$, where $f^{\ast }$ denotes the Fenchel conjugate of $f$. It follows from this that the iterative sequences (45) can be rewritten as$$\begin{array}{c}\text{(46)}& x_k=pro{x}_f\left(u-\gamma L^{\ast }{y}_k\right),y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+\frac{{\lambda }_k}{\gamma }pro{x}_{\gamma g^{\ast }}\left(\gamma y_k+\gamma \left(L{x}_k-r\right)\right).\end{array}$$Let ${\overline{y}}_k=\gamma y_k$. Then, (46) is reduced to$$\begin{array}{c}\text{(47)}& x_k=pro{x}_f\left(u-L^{\ast }{\overline{y}}_k\right),{\overline{y}}_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right){\overline{y}}_k+{\lambda }_{k}pro{x}_{\gamma g^{\ast }}\left({\overline{y}}_k+\gamma \left(L{x}_k-r\right)\right).\end{array}$$

Combettes et al. [27] proposed the following iterative algorithm to solve the optimization problem (44). For any $y_{\mathrm{0}}\in \mathrm{G}$, set $\{x^k\}$ and $\{y^k\}$ as$$\begin{array}{c}\text{(48)}& x_k=pro{x}_f\left(u-L^{\ast }{y}_k\right)+b_k,y_{k+\mathrm{1}}=\left(\mathrm{1}-{\lambda }_k\right)y_k+{\lambda }_{k}pro{x}_{{\gamma }_k{g}^{\ast }}\left(y_k+{\gamma }_k\left(L{x}_k-r\right)+a_k\right),\end{array}$$where ${\lambda }_k\in [ϵ,\mathrm{1}]$, ${\gamma }_k\in (ϵ,\mathrm{2}/{‖L‖}^{\mathrm{2}}-ϵ)$, and $ϵ\in (\mathrm{0},\mathrm{m}\mathrm{i}\mathrm{n}(\mathrm{1},\mathrm{1}/{‖L‖}^{\mathrm{2}}))$. Under the conditions that ${\sum }_{k=\mathrm{0}}^{\infty }‍‖a_k‖<+\infty$ and ${\sum }_{k=\mathrm{0}}^{\infty }‍‖b_k‖<+\infty$, the authors proved that $\{x_k\}$ converges strongly to a solution of (44).

Remark 22.

(1) Comparing (47) with (48), the range of ${\lambda }_k$ in (47) is clearly larger than that of the iterative sequence (48) introduced by Combettes et al. [27] when the variable step size ${\gamma }_k$ is constant; i.e., ${\gamma }_k\equiv \gamma$. In addition, Theorem 20 also recovers the main results of Proposition 28.16 in [28]. However, we require that ${\lambda }_k$ satisfies ${\sum }_{k=\mathrm{0}}^{\infty }‍{\lambda }_k(\mathrm{1}-\left(\mathrm{2}/(\mathrm{4}-\gamma {‖L‖}^{\mathrm{2}})\right){\lambda }_k)=+\infty$, while Proposition 28.16 in [28] requires that $\mathrm{0}<\inf {\lambda }_k\le \sup {\lambda }_k<(\mathrm{4}-\gamma {‖L‖}^{\mathrm{2}})/\mathrm{2}$. It is obvious that our conditions on ${\lambda }_k$ are weaker than those of Proposition 28.16 in [28].

(2) Although our proposed iterative sequences (45) are error-free, it is not difficult to add error sequences in corresponding locations, as in (48). Because the proof is almost identical to that of Theorem 12, we have omitted it here.

5. Numerical Experiments

In this section, we present numerical experiments to verify the effectiveness of the proposed iterative algorithms for solving the constrained total variation model (4) for image denoising problems. All experiments are conducted on a Lenovo laptop with an Intel 2.3 GHz CPU and 4 GB RAM. We run the program with MATLAB R2014a.

We select “Barbara,” “Lena,” “Boat,” and “Goldhill” as the test images (see Figure 1). Gaussian noise of mean $\mathrm{0}$ and standard deviation $\mathrm{30}$ is added to the clear images.

[figures omitted; refer to PDF]

We use the signal-to-noise (SNR) and peak-signal-to-noise (PSNR) to evaluate the quality of the restored images. These are defined by$$\begin{array}{}\text{(49)}& SNR=\mathrm{20}\hspace{0.17em}\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{10}}\frac{{‖f‖}_{\mathrm{2}}}{{‖f-\tilde{x}‖}_{\mathrm{2}}},\end{array}$$and$$\begin{array}{}\text{(50)}& PSNR=\mathrm{20}\hspace{0.17em}\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{10}}\frac{\mathrm{255}\sqrt{MN}}{{‖f-\tilde{x}‖}_{\mathrm{2}}},\end{array}$$where $f$ is the original image, $\tilde{x}$ is the restored image, and $M$ and $N$ are the row and column size of the image $f$, respectively. The iterative algorithm is stopped when the following criterion is satisfied:$$\begin{array}{}\text{(51)}& \frac{{‖y_{k+\mathrm{1}}-y_k‖}_{\mathrm{2}}}{{‖y_k‖}_{\mathrm{2}}}\le ϵ,\end{array}$$where $ϵ$ is a given small constant. We tuned the regularization parameter $\mu$ and set it to $\mathrm{15}$ for optimal denoised image quality.

We aim mainly to solve the constrained total variation (TV) image denoising problem (4). In particular, we choose the anisotropic total variation as the regularization term during testing. By using the indicator function, the constrained (TV) denoising problem (4) can be reformulated as the following unconstrained optimization problem:$$\begin{array}{}\text{(52)}& \underset{x\in R^n}{\mathrm{m}\mathrm{i}\mathrm{n}}\hspace{1em}\frac{\mathrm{1}}{\mathrm{2}}{‖x-f‖}_{\mathrm{2}}^{\mathrm{2}}+\mu {‖x‖}_{TV}+{\delta }_C\left(x\right),\end{array}$$where ${\delta }_C$ is the indicator function. It is clear that the optimization problem (52) is a special case of (44). In fact, if $f(x)={\delta }_C(x)$, $g(Lx)={‖x‖}_{TV}$, and $r=\mathrm{0}$, then the proposed iteration scheme (45) can be employed to solve the constrained TV image denoising problem (4). It is well known that ${‖L‖}^{\mathrm{2}}\approx \mathrm{8}$, where $L$ is the usual gradient operator of the total variation (see [10, 20]).

Let $C=R^n$. Then, (4) (also (52)) is reduced to the well-known ROF model (1). We select $C$ as both the nonnegative and the bounded sets. Moreover, the nonnegative set $C$ is defined as $\{x\in R^n\mid x_i\ge \mathrm{0}\}$, and the bounded set $C$ as $\{x\in R^n\mid \mathrm{0}\le x_i\le \mathrm{255}\}$. The corresponding minimization problems (52) are referred to as the nonnegative ROF (N-ROF) model and bounded ROF (B-ROF) model, respectively. It is worth mentioning that the proximity operator of the indicator function is the projection operator, which has closed-form solutions based on our choices. Therefore, in (45) the proximity operator of $f(x)={\delta }_C(x)$ is the orthogonal projection onto the closed convex set $C$. That is, $pro{x}_{{\delta }_C}(u)=P_C(u)$.

5.1. Numerical Results and Discussion

First, we describe the impacts of the parameters for the iterative step size $\gamma$ and relaxation variable ${\lambda }_k$ on the proposed iterative algorithm (45). According to Theorem 20, we can choose $\gamma \in (\mathrm{0,1}/\mathrm{4})$ and ${\lambda }_k\in (\mathrm{0,2}-\mathrm{4}\gamma )$. Figures 2, 3, and 4 demonstrate the performance of the proposed iterative algorithm for solving (52) with different choices of $\gamma$ and ${\lambda }_k$. The test image was “Barbara,” and the stopping criterion $ϵ$ was set as $\mathrm{1}{\mathrm{0}}^{-\mathrm{2}}$.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

As shown in Figures 2–4, when the iterative step size $\gamma$ is fixed, a larger relaxation parameter ${\lambda }_k$ leads to a faster convergence of the iterative algorithm. Table 2 presents the corresponding numerical results in terms of the SNR, PSNR, and number of iterations ($Iter$).

Table 2

Numerical results for different choices of $\gamma$ and ${\lambda }_k$ in terms of the SNR, PSNR, and number of iterations ($Iter$) for the test image “Barbara.”

$\gamma$	${\lambda }_k$	ROF			N-ROF			B-ROF
		$SNR$	$PSNR$	$Iter$	$SNR$	$PSNR$	$Iter$	$SNR$	$PSNR$	$Iter$
$\mathrm{1}/\mathrm{16}$	$\mathrm{0.2}$	$\mathrm{18.7102}$	$\mathrm{24.5975}$	$\mathrm{41}$	$\mathrm{18.7073}$	$\mathrm{24.5946}$	$\mathrm{41}$	$\mathrm{18.7063}$	$\mathrm{24.5936}$	$\mathrm{41}$
	$\mathrm{0.4}$	$\mathrm{18.8872}$	$\mathrm{24.7745}$	$\mathrm{34}$	$\mathrm{18.8925}$	$\mathrm{24.7798}$	$\mathrm{35}$	$\mathrm{18.8924}$	$\mathrm{24.7797}$	$\mathrm{35}$
	$\mathrm{0.6}$	$\mathrm{18.9328}$	$\mathrm{24.8200}$	$\mathrm{31}$	$\mathrm{18.9331}$	$\mathrm{24.8204}$	$\mathrm{31}$	$\mathrm{18.9332}$	$\mathrm{24.8205}$	$\mathrm{31}$
	$\mathrm{0.8}$	$\mathrm{18.9520}$	$\mathrm{24.8393}$	$\mathrm{29}$	$\mathrm{18.9528}$	$\mathrm{24.8401}$	$\mathrm{29}$	$\mathrm{18.9530}$	$\mathrm{24.8403}$	$\mathrm{29}$
	$\mathrm{1}$	$\mathrm{18.9633}$	$\mathrm{24.8506}$	$\mathrm{28}$	$\mathrm{18.9644}$	$\mathrm{24.8516}$	$\mathrm{28}$	$\mathrm{18.9646}$	$\mathrm{24.8519}$	$\mathrm{28}$
	$\mathrm{1.2}$	$\mathrm{18.9680}$	$\mathrm{24.8553}$	$\mathrm{26}$	$\mathrm{18.9710}$	$\mathrm{24.8582}$	$\mathrm{27}$	$\mathrm{18.9712}$	$\mathrm{24.8585}$	$\mathrm{27}$
	$\mathrm{1.5}$	$\mathrm{18.9745}$	$\mathrm{24.8618}$	$\mathrm{25}$	$\mathrm{18.9758}$	$\mathrm{24.8631}$	$\mathrm{25}$	$\mathrm{18.9761}$	$\mathrm{24.8634}$	$\mathrm{25}$
	$\mathrm{1.75}$	$\mathrm{18.9775}$	$\mathrm{24.8648}$	$\mathrm{24}$	$\mathrm{18.9789}$	$\mathrm{24.8661}$	$\mathrm{24}$	$\mathrm{18.9791}$	$\mathrm{24.8664}$	$\mathrm{24}$
$\mathrm{1}/\mathrm{8}$	$\mathrm{0.2}$	$\mathrm{18.8906}$	$\mathrm{24.7778}$	$\mathrm{34}$	$\mathrm{18.8897}$	$\mathrm{24.7770}$	$\mathrm{34}$	$\mathrm{18.8896}$	$\mathrm{24.7769}$	$\mathrm{34}$
	$\mathrm{0.4}$	$\mathrm{18.9529}$	$\mathrm{24.8401}$	$\mathrm{29}$	$\mathrm{18.9536}$	$\mathrm{24.8409}$	$\mathrm{29}$	$\mathrm{18.9538}$	$\mathrm{24.8411}$	$\mathrm{29}$
	$\mathrm{0.6}$	$\mathrm{18.9683}$	$\mathrm{24.8556}$	$\mathrm{26}$	$\mathrm{18.9695}$	$\mathrm{24.8568}$	$\mathrm{26}$	$\mathrm{18.9697}$	$\mathrm{24.8570}$	$\mathrm{26}$
	$\mathrm{0.8}$	$\mathrm{18.9753}$	$\mathrm{24.8625}$	$\mathrm{24}$	$\mathrm{18.9766}$	$\mathrm{24.8639}$	$\mathrm{24}$	$\mathrm{18.9769}$	$\mathrm{24.8641}$	$\mathrm{24}$
	$\mathrm{1}$	$\mathrm{18.9795}$	$\mathrm{24.8667}$	$\mathrm{23}$	$\mathrm{18.9809}$	$\mathrm{24.8691}$	$\mathrm{23}$	$\mathrm{18.9811}$	$\mathrm{24.8684}$	$\mathrm{23}$
	$\mathrm{1.2}$	$\mathrm{18.9816}$	$\mathrm{24.8688}$	$\mathrm{22}$	$\mathrm{18.9831}$	$\mathrm{24.8703}$	$\mathrm{22}$	$\mathrm{18.9833}$	$\mathrm{24.8706}$	$\mathrm{22}$
	$\mathrm{1.5}$	$\mathrm{18.9829}$	$\mathrm{24.8702}$	$\mathrm{20}$	$\mathrm{18.9844}$	$\mathrm{24.8717}$	$\mathrm{20}$	$\mathrm{18.9847}$	$\mathrm{24.8720}$	$\mathrm{20}$
$\mathrm{1}/\mathrm{4}$	$\mathrm{0.2}$	$\mathrm{18.9543}$	$\mathrm{24.8416}$	$\mathrm{28}$	$\mathrm{18.9547}$	$\mathrm{24.8420}$	$\mathrm{28}$	$\mathrm{18.9548}$	$\mathrm{24.8421}$	$\mathrm{28}$
	$\mathrm{0.4}$	$\mathrm{18.9762}$	$\mathrm{24.8635}$	$\mathrm{24}$	$\mathrm{18.9774}$	$\mathrm{24.8647}$	$\mathrm{24}$	$\mathrm{18.9777}$	$\mathrm{24.8649}$	$\mathrm{24}$
	$\mathrm{0.6}$	$\mathrm{18.9813}$	$\mathrm{24.8686}$	$\mathrm{21}$	$\mathrm{18.9834}$	$\mathrm{24.8707}$	$\mathrm{22}$	$\mathrm{18.9837}$	$\mathrm{24.8710}$	$\mathrm{22}$
	$\mathrm{0.8}$	$\mathrm{18.9836}$	$\mathrm{24.8709}$	$\mathrm{20}$	$\mathrm{18.9851}$	$\mathrm{24.8724}$	$\mathrm{20}$	$\mathrm{18.9854}$	$\mathrm{24.8726}$	$\mathrm{20}$
	$\mathrm{1}$	$\mathrm{18.9841}$	$\mathrm{24.8714}$	$\mathrm{19}$	$\mathrm{18.9856}$	$\mathrm{24.8729}$	$\mathrm{19}$	$\mathrm{18.9859}$	$\mathrm{24.8732}$	$\mathrm{19}$

Because the prior pixel information of the image is introduced as a constraint, the performance of the constrained ROF model is superior to that of the unconstrained model. Numerical results confirm the advantages of the constrained ROF model. We can see from Table 2 that the iterative step size $\gamma$ has a significant impact on the performance of the algorithm. Furthermore, the experimental results demonstrate that a larger $\gamma$ always leads to a faster convergence of the iterative algorithm. Thus, we choose $\gamma =\mathrm{1}/\mathrm{4}$ and ${\lambda }_k=\mathrm{1}$ in the following comparison.

Next, we focus on investigating the performances of the constrained and unconstrained ROF models for the test images from Figure 1. The numerical results are presented in Table 3. We notice that the SNR and PSNR values slowly decrease with an increasing number of iterations. Because more iterations do not improve the quality of the image, we should stop the iterative algorithm in the early stages. Figures 5–8 present denoised images for $ϵ=\mathrm{1}{\mathrm{0}}^{-\mathrm{2}}$ in Table 3. From the visual point of view, the images restored by the constrained model are closer to the original images than those restored by the unconstrained model. This further confirms the benefits of introducing constraints in the image recovery model.

Table 3

Numerical results in terms of the SNR, PSNR, and number of iterations ($Iter$) on different test images for different tolerance values $ϵ$.

Images	Model	$ϵ=\mathrm{1}{\mathrm{0}}^{-\mathrm{2}}$			$ϵ=\mathrm{1}{\mathrm{0}}^{-\mathrm{3}}$			$ϵ=\mathrm{1}{\mathrm{0}}^{-\mathrm{4}}$
		$SNR$	$PSNR$	$Iter$	$SNR$	$PSNR$	$Iter$	$SNR$	$PSNR$	$Iter$
“Barbara”	ROF	$\mathrm{18.9841}$	$\mathrm{24.8714}$	$\mathrm{19}$	$\mathrm{18.9800}$	$\mathrm{24.8672}$	$\mathrm{71}$	$\mathrm{18.9784}$	$\mathrm{24.8656}$	$\mathrm{216}$
	N-ROF	$\mathrm{18.9856}$	$\mathrm{24.8729}$	$\mathrm{19}$	$\mathrm{18.9816}$	$\mathrm{24.8688}$	$\mathrm{71}$	$\mathrm{18.9800}$	$\mathrm{24.8672}$	$\mathrm{216}$
	B-ROF	$\mathrm{18.9859}$	$\mathrm{24.8732}$	$\mathrm{19}$	$\mathrm{18.9818}$	$\mathrm{24.8691}$	$\mathrm{71}$	$\mathrm{18.9802}$	$\mathrm{24.8675}$	$\mathrm{216}$
“Lena”	ROF	$\mathrm{22.4595}$	$\mathrm{28.1409}$	$\mathrm{21}$	$\mathrm{22.4552}$	$\mathrm{28.1366}$	$\mathrm{78}$	$\mathrm{22.4522}$	$\mathrm{28.1336}$	$\mathrm{242}$
	N-ROF	$\mathrm{22.4605}$	$\mathrm{28.1419}$	$\mathrm{21}$	$\mathrm{22.4563}$	$\mathrm{28.1377}$	$\mathrm{78}$	$\mathrm{22.4534}$	$\mathrm{28.1347}$	$\mathrm{242}$
	B-ROF	$\mathrm{22.4616}$	$\mathrm{28.1430}$	$\mathrm{21}$	$\mathrm{22.4573}$	$\mathrm{28.1387}$	$\mathrm{78}$	$\mathrm{22.4544}$	$\mathrm{28.1358}$	$\mathrm{242}$
“Boat”	ROF	$\mathrm{21.6226}$	$\mathrm{26.9652}$	$\mathrm{20}$	$\mathrm{21.6180}$	$\mathrm{26.9606}$	$\mathrm{74}$	$\mathrm{21.6162}$	$\mathrm{26.9588}$	$\mathrm{228}$
	N-ROF	$\mathrm{21.6326}$	$\mathrm{26.9752}$	$\mathrm{20}$	$\mathrm{21.6278}$	$\mathrm{26.9704}$	$\mathrm{74}$	$\mathrm{21.6260}$	$\mathrm{26.9686}$	$\mathrm{228}$
	B-ROF	$\mathrm{21.6335}$	$\mathrm{26.9761}$	$\mathrm{20}$	$\mathrm{21.6287}$	$\mathrm{26.9713}$	$\mathrm{74}$	$\mathrm{21.6269}$	$\mathrm{26.9695}$	$\mathrm{228}$
“Goldhill”	ROF	$\mathrm{20.9503}$	$\mathrm{27.3164}$	$\mathrm{20}$	$\mathrm{20.9353}$	$\mathrm{27.3015}$	$\mathrm{73}$	$\mathrm{20.9328}$	$\mathrm{27.2989}$	$\mathrm{223}$
	N-ROF	$\mathrm{20.9527}$	$\mathrm{27.3189}$	$\mathrm{20}$	$\mathrm{20.9379}$	$\mathrm{27.3040}$	$\mathrm{73}$	$\mathrm{20.9353}$	$\mathrm{27.3014}$	$\mathrm{223}$
	B-ROF	$\mathrm{20.9577}$	$\mathrm{27.3239}$	$\mathrm{20}$	$\mathrm{20.9429}$	$\mathrm{27.3091}$	$\mathrm{73}$	$\mathrm{20.9403}$	$\mathrm{27.3065}$	$\mathrm{223}$

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

6. Conclusions

The total variation can be viewed as a composition of a convex function with a linear transformation. Thus, Micchelli et al. [20] introduced a fixed-point algorithm based on proximity operators to produce a total variation model for image denoising (1). Inspired by the work of Moudafi [26], we studied the calculation of the resolvent of the sum of a maximal monotone operator and a composite operator (9) to produce a constrained total variation model (4). Subsequently, we proposed a fixed-point algorithm for this resolvent operator. Based on the fixed-point theory of nonexpansive mappings, we proved the strong convergence of the obtained iterative sequence. The advantage of the fixed-point approach is that it provides the potential to develop additional fast iterative algorithms. Numerical simulations on image denoising illustrated the performance of the proposed algorithm. In particular, we found that the step size had a significant impact on the convergence speed of the algorithm. In general, when the iterative step size was fixed, larger relaxation parameters resulted in a faster iterative algorithm convergence. Numerical results also confirmed that the constrained ROF model achieved a superior performance compared with the unconstrained ROF model.

Finally, we wish to note that the constrained TV model (4) can also be derived using other iterative algorithms, such as the primal-dual Chambolle–Pock algorithm [15], the alternating direction method of multipliers [29, 38, 39], and the preconditioned primal-dual algorithm [40, 41]. We have not presented the corresponding numerical results here. Thus, we will further examine the convergence rate of our proposed iterative algorithm and include these comparative results in future work.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

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

References

[1] L. I. Rudin, S. Osher, E. Fatemi, "Nonlinear total variation based noise removal algorithms," Physica D: Nonlinear Phenomena, vol. 60 no. 1–4, pp. 259-268, DOI: 10.1016/0167-2789(92)90242-F, 1992.

[2] Y. Wang, J. Yang, W. Yin, Y. Zhang, "A new alternating minimization algorithm for total variation image reconstruction," SIAM Journal on Imaging Sciences, vol. 1 no. 3, pp. 248-272, DOI: 10.1137/080724265, 2008.

[3] A. Beck, M. Teboulle, "Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems," IEEE Transactions on Image Processing, vol. 18 no. 11, pp. 2419-2434, DOI: 10.1109/TIP.2009.2028250, 2009.

[4] C. Wu, J. Zhang, X.-C. Tai, "Augmented Lagrangian method for total variation restoration with non-quadratic fidelity," Inverse Problems and Imaging, vol. 5 no. 1, pp. 237-261, DOI: 10.3934/ipi.2011.5.237, 2011.

[5] T. F. Chan, J. Shen, "Mathematical models for local nontexture inpaintings," SIAM Journal on Applied Mathematics, vol. 62 no. 3, pp. 1019-1043, DOI: 10.1137/S0036139900368844, 2002.

[6] A. Marquina, S. J. Osher, "Image super-resolution by TV-regularization and Bregman iteration," Journal of Scientific Computing, vol. 37 no. 3, pp. 367-382, DOI: 10.1007/s10915-008-9214-8, 2008.

[7] J. Liu, Y.-B. Ku, S. Leung, "Expectation-maximization algorithm with total variation regularization for vector-valued image segmentation," Journal of Visual Communication and Image Representation, vol. 23 no. 8, pp. 1234-1244, DOI: 10.1016/j.jvcir.2012.09.002, 2012.

[8] E. Y. Sidky, X. Pan, "Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization," Physics in Medicine and Biology, vol. 53 no. 17, pp. 4777-4807, DOI: 10.1088/0031-9155/53/17/021, 2008.

[9] Y. Chen, W. W. Hager, M. Yashtini, X. Ye, H. Zhang, "Bregman operator splitting with variable stepsize for total variation image reconstruction," Computational optimization and applications, vol. 54 no. 2, pp. 317-342, DOI: 10.1007/s10589-012-9519-2, 2013.

[10] A. Chambolle, "An algorithm for total variation minimization and applications," Journal of Mathematical Imaging and Vision, vol. 20 no. 1-2, pp. 89-97, DOI: 10.1023/B:JMIV.0000011320.81911.38, 2004.

[11] G. Yu, L. Qi, Y. Dai, "On nonmonotone Chambolle gradient projection algorithms for total variation image restoration," Journal of Mathematical Imaging and Vision, vol. 35 no. 2, pp. 143-154, DOI: 10.1007/s10851-009-0160-3, 2009.

[12] J.-F. Aujol, "Some first-order algorithms for total variation based image restoration," Journal of Mathematical Imaging and Vision, vol. 34 no. 3, pp. 307-327, DOI: 10.1007/s10851-009-0149-y, 2009.

[13] M. Zhu, T. F. Chan, "An efficient primal-dual hybrid gradient algorithm for total variation image restoration," Technical Report CAM Report, May 2008.

[14] E. Esser, X. Zhang, T. F. Chan, "A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science," SIAM Journal on Imaging Sciences, vol. 3 no. 4, pp. 1015-1046, DOI: 10.1137/09076934X, 2010.

[15] A. Chambolle, T. Pock, "A first-order primal-dual algorithm for convex problems with applications to imaging," Journal of Mathematical Imaging and Vision, vol. 40 no. 1, pp. 120-145, DOI: 10.1007/s10851-010-0251-1, 2011.

[16] T. Goldstein, S. Osher, "The split Bregman method for L 1-regularized problems," SIAM Journal on Imaging Sciences, vol. 2 no. 2, pp. 323-343, DOI: 10.1137/080725891, 2009.

[17] C. Wu, X.-C. Tai, "Augmented Lagrangian method, dual methods, and split Bregman iteration for ROF, vectorial TV, and high order models," SIAM Journal on Imaging Sciences, vol. 3 no. 3, pp. 300-339, DOI: 10.1137/090767558, 2010.

[18] G. Steidl, "A note on the dual treatment of higher-order regularization functionals," Computing: Archives for Scientific Computing, vol. 76 no. 1-2, pp. 135-148, DOI: 10.1007/s00607-005-0129-z, 2006.

[19] P. Chen, J. Huang, X. Zhang, "A primal-dual fixed point algorithm for convex separable minimization with applications to image restoration," Inverse Problems, vol. 29 no. 2,DOI: 10.1088/0266-5611/29/2/025011, 2013.

[20] C. A. Micchelli, L. Shen, Y. Xu, "Proximity algorithms for image models: denoising," Inverse Problems, vol. 27 no. 4,DOI: 10.1088/0266-5611/27/4/045009, 2011.

[21] H. H. Bauschke, P. L. Combettes, "A Dykstra-like algorithm for two monotone operators," Pacific Journal of Optimization, vol. 4 no. 3, pp. 383-391, 2008.

[22] R. L. Dykstra, "An algorithm for restricted least squares regression," Journal of the American Statistical Association, vol. 78 no. 384, pp. 837-842, DOI: 10.1080/01621459.1983.10477029, 1983.

[23] P. L. Combettes, "Iterative construction of the resolvent of a sum of maximal monotone operators," Journal of Convex Analysis, vol. 16 no. 3-4, pp. 727-748, 2009.

[24] F. J. Aragón Artacho, R. Campoy, "Computing the resolvent of the sum of maximally monotone operators with the averaged alternating modified reflections algorithm," Journal of Optimization Theory and Applications,DOI: 10.1007/s10957-019-01481-3, 2019.

[25] F. J. Aragón Artacho, R. Campoy, "A new projection method for finding the closest point in the intersection of convex sets," Computational optimization and applications, vol. 69 no. 1, pp. 99-132, DOI: 10.1007/s10589-017-9942-5, 2018.

[26] A. Moudafi, "Computing the resolvent of composite operators," Cubo: A Mathematical Journal, vol. 16 no. 3, pp. 87-96, DOI: 10.4067/s0719-06462014000300007, 2014.

[27] P. L. Combettes, D. Dung, B. C. Vu, "Dualization of signal recovery problems," Set-Valued and Variational Analysis, vol. 18 no. 3-4, pp. 373-404, DOI: 10.1007/s11228-010-0147-7, 2010.

[28] H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces,DOI: 10.1007/978-3-319-48311-5, 2017.

[29] R. H. Chan, M. Tao, X. Yuan, "Constrained total variation deblurring models and fast algorithms based on alternating direction method of multipliers," SIAM Journal on Imaging Sciences, vol. 6 no. 1, pp. 680-697, DOI: 10.1137/110860185, 2013.

[30] T. Pennanen, "Dualization of generalized equations of maximal monotone type," SIAM Journal on Optimization, vol. 10 no. 3, pp. 809-835, DOI: 10.1137/S1052623498340448, 2000.

[31] C. Byrne, "A unified treatment of some iterative algorithms in signal processing and image reconstruction," Inverse Problems, vol. 20 no. 1, pp. 103-120, DOI: 10.1088/0266-5611/20/1/006, 2004.

[32] N. Ogura, I. Yamada, "Non-strictly convex minimization over the fixed point set of an asymptotically shrinking nonexpansive mapping," Numerical Functional Analysis and Optimization, vol. 23 no. 1-2, pp. 113-137, DOI: 10.1081/NFA-120003674, 2002.

[33] P. L. Combettes, I. Yamada, "Compositions and convex combinations of averaged nonexpansive operators," Journal of Mathematical Analysis and Applications, vol. 425 no. 1, pp. 55-70, DOI: 10.1016/j.jmaa.2014.11.044, 2015.

[34] L. M. Briceño-Arias, P. L. Combettes, "A monotone+skew splitting model for composite monotone inclusions in duality," SIAM Journal on Optimization, vol. 21 no. 4, pp. 1230-1250, DOI: 10.1137/10081602x, 2011.

[35] P. L. Combettes, J.-C. Pesquet, "Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators," Set-Valued and Variational Analysis, vol. 20 no. 2, pp. 307-330, DOI: 10.1007/s11228-011-0191-y, 2012.

[36] B. C. Vũ, "A splitting algorithm for dual monotone inclusions involving cocoercive operators," Advances in Computational Mathematics, vol. 38 no. 3, pp. 667-681, DOI: 10.1007/s10444-011-9254-8, 2013.

[37] L. M. Briceno-Arias, D. Davis, "Forward-backward-half forward algorithm for solving monotone inclusions," SIAM Journal on Optimization, vol. 28 no. 4, pp. 2839-2871, DOI: 10.1137/17M1120099, 2018.

[38] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers," Foundations and Trends in Machine Learning, vol. 3 no. 1,DOI: 10.1561/2200000016, 2010.

[39] M. K. Ng, P. Weiss, X. M. Yuan, "Solving constrained total-variation image restoration and reconstruction problems via alternating direction methods," SIAM Journal on Scientific Computing, vol. 32 no. 5, pp. 2710-2736, DOI: 10.1137/090774823, 2010.

[40] T. Pock, A. Chambolle, "Diagonal preconditioning for first order primal-dual algorithms in convex optimization," Proceedings of the IEEE International Conference on Computer Vision (ICCV '11), pp. 1762-1769, DOI: 10.1109/ICCV.2011.6126441, .

[41] M. Wen, J. Peng, Y. Tang, C. Zhu, S. Yue, "A preconditioning technique for first-order primal-dual splitting method in convex optimization," Mathematical Problems in Engineering, vol. 2017,DOI: 10.1155/2017/3694525, 2017.