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

Authors: Bao Chen; Yuchao Tang
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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(1)minx∈Rn12x-u2+μxTV,where u∈Rn is the observed noisy image, xTV is the total variation, and μ>0 is the regularization parameter, which balances the data fidelity term and total variation regularization. In the following, we denote by x∈Rn the vectorized image X∈RM×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 φ with a linear transformation B. That is, xTV=φ(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(2)minx∈Rn12x-u2+ϕBx,where u∈Rn, ϕ:Rm→(-∞,+∞] is a proper lower semicontinuous convex function and B:Rn→Rm is a linear transformation. The convex optimization problem (2) is equivalent to computing the proximity operator of the function ϕ∘B. Recall that the proximity operator of a proper lower semicontinuous convex function f:Rn→(-∞,+∞] is defined as(3)proxfu=argminx∈Rn12x-u2+fx.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[0,255] or [0,1], it is natural to incorporate this prior information into the ROF model (1). Thus, we obtain the following constrained ROF model:(4)minx∈C12x-u2+μxTV,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:(5)minx∈C12x-u2+ϕBx.In general, the constrained convex optimization problem (5) can always be reformulated as an unconstrained convex optimization problem. More precisely,(6)minx∈Rn12x-u2+ϕBx+δCx,where δC denotes the indicator function of C, which is defined by δC(x)=0 if x∈C and δC(x)=+∞ 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δ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:(7)minx∈H12x-u2+ϕBx+gx.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→2H be a maximal monotone operator. The resolvent of T is the single-valued operator JT=(I+T)-1. Next, let f:H→(-∞,+∞] be a proper lower semicontinuous convex function and let ∂f denote the subdifferential of f, and let T=∂f. Thus, the resolvent of T coincides with the proximity operator proxf as follows:(8)J∂fu=proxfu.Under certain qualification conditions, the problem considered in (7) can be solved via the resolvent operator JB∗∘∂ϕ∘B+∂g. Overall, this leads to the problem of computing the resolvent operator of S+B∗TB. More precisely, let u∈H and(9)x≔JS+B∗TBu,where S:H→2H and T:G→2G are two maximal monotone operators and B:H→G is a bounded linear operator from the Hilbert space H to the Hilbert space G. Following this, let S=∂g and T=∂ϕ. Hence, the problem (9) reduces to the convex optimization problem (7), because J∂g+B∗∘∂ϕ∘B(u)=proxg+ϕ∘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∗TB, where B:H→G is a bounded linear operator with the adjoint B∗, T:G→2G is a maximal monotone operator, and H and G are two Hilbert spaces. When T=∂ϕ, this algorithm utilizes the fixed-point approach to computing the proximity operator proxϕ∘B proposed by Micchelli et al. [20]. It is clear that the resolvent of the operator B∗TB is a special case of (9) by allowing S=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 operatorS+B∗TB (9), where S:H→2H and T:G→2G are two maximal monotone operators and B:H→G is a bounded linear operator. We assume that the operator S+B∗TB is maximally monotone. This is true, for example, if 0∈ri(B(domS)-domT) [30], where ri(C) represents a relative interior of the set C. If S=0, then the problem (9) becomes that of computing the resolvent operator of B∗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 Section2, 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.


## 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∗TB, where B:H→G is a bounded linear operator with the adjoint B∗, T:G→2G is a maximal monotone operator, and H and G are two Hilbert spaces. When T=∂ϕ, this algorithm utilizes the fixed-point approach to computing the proximity operator proxϕ∘B proposed by Micchelli et al. [20]. It is clear that the resolvent of the operator B∗TB is a special case of (9) by allowing S=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 operatorS+B∗TB (9), where S:H→2H and T:G→2G are two maximal monotone operators and B:H→G is a bounded linear operator. We assume that the operator S+B∗TB is maximally monotone. This is true, for example, if 0∈ri(B(domS)-domT) [30], where ri(C) represents a relative interior of the set C. If S=0, then the problem (9) becomes that of computing the resolvent operator of B∗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 Section2, 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, letH be a real Hilbert space with inner product 〈·,·〉 and the associated norm ·. Let I denote the identity operator on H. We use xk⇀x to indicate that the sequence {xk} converges weakly to x and xk→x to indicate that the sequence {xk} 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∈L(H,G). Then, the adjoint of B is the unique operator B∗∈L(G,H) such that 〈Bx,y〉=〈x,B∗y〉.LetA:H→2H be a set-valued operator. The sets domA={x∈H∣Ax≠∅}, ranA={u∈H∣∃x∈H,u∈Ax}, graA={(x,u)∈H×H∣u∈Ax}, and zerA={x∈H∣0∈Ax} represent the domain, range, graph, and zero set of A, respectively. The inverse A-1 of A is the set-valued operator with the graph {(u,x)∈H×H∣u∈Ax}.We now introduce some definitions and lemmas, most of which can be found in [28, 31].Definition 1 (see [28], (maximal monotone operator)).
LetA:H→2H be a set-valued operator. A is said to be monotone if(10)x-y,u-v≥0,∀x,u∈graA,y,v∈graA.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)).
LetA:H→2H and γ>0. The resolvent of γA is JγA=(I+γA)-1, and the Yosida approximation of γA is γA=1/γ(I-JγA).For anyx∈H, it is easy to derive that γAx∈A(JγAx).Lemma 3 (see [28]).
LetA:H→2H be maximally monotone and γ>0. Then, the following hold:(i)
JγA:H→H, and I-γγA:H→H are firmly nonexpansive and maximally monotone.(ii)
The Yosida approximationγA is γ-cocoercive and maximally monotone.Definition 4 (see [28]).
LetC be a nonempty subset of H and T:C→C. Then(i)
T is nonexpansive if Tx-Ty≤x-y∀x,y∈C.(ii)
T is firmly nonexpansive if Tx-Ty2≤x-y2-(I-T)x-(I-T)y2∀x,y∈C.(iii)
T is α-averaged if there exists a nonexpansive S and α∈(0,1) such that T=(1-α)I+αS, where I denotes the identity operator.(iv)
T is β-cocoercive (or β-inversely strong monotone) for β>0 if 〈x-y,Tx-Ty〉≥βTx-Ty2∀x,y∈C.Remark 5.
(i) An equivalent definition of firm nonexpansiveness is thatTx-Ty2≤〈x-y,Tx-Ty〉∀x,y∈C. It is easy to check that firm nonexpansiveness implies nonexpansiveness. If T is firmly nonexpansive, then T is 1-cocoercive.
(ii) An equivalent definition ofα-averaged is Tx-Ty2≤x-y2-(1-α)/α(I-T)x-(I-T)y2 for all x,y∈C. Thus, if T is α- averaged then T is also nonexpansive.
(iii) LetT be β-cocoercive. Then, for any γ>0γT is β/γ-cocoercive. Furthermore, owing to the Cauchy-Schwartz inequality T is 1/β- Lipschitz continuous; i.e., Tx-Ty≤1/βx-y.The following lemma provides some useful characterizations between an operatorT and I-T.Lemma 6 (see [28]).
LetC be a nonempty subset of H and let T:C→C. Then(i)
T is nonexpansive if and only if I-T is 1/2-cocoercive.(ii)
T is firmly nonexpansive if and only if I-T is firmly nonexpansive.(iii)
T is 1/2v-averaged if and only if I-T is v-cocoercive, where v>1/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]).
LetC be a nonempty subset of H. Furthermore, let T1:C→C be α1- averaged and let T2:C→C be α2- averaged. Then T1T2 is α1+α2-2α1α2/1-α1α2 averaged.We also make full use of the following lemma.Lemma 8 (see [28]).
Letx∈H and y∈H. Let α∈R and R denote the set of real numbers. Then(11)αx+1-αy2+α1-αx-y2=αx2+1-αy2.We end this section by introducing the Krasnoselskii–Mann algorithm. Theorem9 provides a fundamental tool for studying the convergence of many operator splitting methods.Theorem 9 (see [28], (Krasnoselskii–Mann algorithm)).
LetC be a nonempty closed convex subset of H, and let T:C→C be a nonexpansive operator such that Fix(T)≠∅, where Fix(T) denotes the fixed-point set of T. Let λk∈(0,1) be such that ∑k=0∞‍λk(1-λk)=+∞. For any x0∈C, define(12)xk+1=1-λkxk+λkTxk.Then, the following hold:(i)
{xk} is Fejér monotone with respect to Fix(T); i.e., xk+1-p≤xk-p for any p∈Fix(T).(ii)
xk-Txk converges strongly to 0.(iii)
{xk} 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 fixedu∈H, let γ>0 define operators F:G→G and Q:G→G as follows:(13)Fy=y+BJSu-γB∗y,∀y∈G,and(14)Qy=I-J1/γTFy,∀y∈G.In addition, the following lemma provides a fixed-point characterization of the resolvent operator (9).Lemma 10.
LetS:H→2H and T:G→2G be two maximal monotone operators. Furthermore, let B∈L(H,G) and u∈H. Then, for a given γ>0,(15)JS+B∗TBu=JSu-γB∗yif and only if y is a fixed point of Q.Proof.
It follows from the definition of the resolvent operator thatJS+B∗TB(u)=(I+S+B∗TB)-1(u). Thus, we can write(16)JS+B∗TBu∈u-S+B∗TBJS+B∗TBu.Let y satisfy (15). Then we have u-γB∗y∈(I+S)JS+B∗TB(u); that is JS+B∗TB(u)∈u-(γB∗y+SJS+B∗TB(u)). By (16), we have y∈1/γTBJS+B∗TB(u). Thus, y∈1/γTBJS(u-γB∗y). Using the fact that y∈1/γTx⇔x=J1/γT(x+y), we deduce that(17)BJSu-γB∗y=J1/γTBJSu-γB∗y+y,which implies that y is a fixed point of Q.
Conversely, lety be a fixed point of Q. Then, we have that(18)y=I-J1/γTy+BJSu-γB∗y⇔BJSu-γB∗y=J1/γTy+BJSu-γB∗y⇔y+BJSu-γB∗y∈I+1γTBJSu-γB∗y⇔γy∈TBJSu-γB∗y⇒u∈u-γB∗y+B∗TBJSu-γB∗y⇔u∈I+SJSu-γB∗y+B∗TBJSu-γB∗y⇔JS+B∗TBu=JSu-γB∗y.This completes the proof.Next, we prove the following lemma, which characterizes an important property of the operatorF.Lemma 11.
LetS:H→2H be a maximal monotone operator, and let B∈L(H,G). For a given u∈H and γ>0, define an operator R:G→G:y↦-BJS(u-γB∗y). Then(i)
R is 1/γB2-cocoercive.(ii)
For anyγ∈(0,2/B2), F=I-R is γB2/2-averaged.Proof.
(i) Lety1,y2∈G. Then, we have that(19)y1-y2,Ry1-Ry2=y1-y2,-BJSu-γB∗y1+BJSu-γB∗y2=1γu-γB∗y1-u-γB∗y2,JSu-γB∗y1-JSu-γB∗y2≥1γJSu-γB∗y1-JSu-γB∗y22≥1γB2-BJSu-γB∗y1--BJSu-γB∗y22=1γB2Ry1-Ry22.The first inequality follows from the fact that the resolvent operator is firmly nonexpansive, and the second is trivial. Thus, R is 1/γB2-cocoercive.
(ii) For anyγ∈(0,2/B2), we know from Lemma 6(iii) that F=I-R is γB2/2-averaged.Lemma11 shows that, for any γ∈(0,2/B2), F is γB2/2-averaged. On the other hand, I-J1/γT is firmly nonexpansive and is also 1/2-averaged. According to Lemma 7, Q=(I-J1/γT)∘F is 2/(4-γB2)-averaged.Now, we are ready to present our main results.Theorem 12.
LetS:H→2H and T:G→2G be two maximal monotone operators. Furthermore, let B∈L(H,G) and u∈H. For any y0∈G, we define the iterative sequences {xk} and {yk} as follows:(20)xk=JSu-γB∗yk,yk+1=1-λkyk+λkI-J1/γTyk+Bxk,where γ∈(0,2/B2) and λk∈(0,(4-γB2)/2) such that ∑k=0∞‍λk(1-2/(4-γB2)λk)=+∞. Then, the following hold:(i)
{yk} converges weakly to a fixed point of Q.(ii)
Suppose thatinf⁡λk>0. Then, {xk} converges strongly to the resolvent operator (9).Proof.
(i) Because the resolvent operatorJS+B∗TB(u) exists and is unique, the first part of the Lemma 10 implies that Q has a fixed point. Let Q(yk)=(I-J1/γT)(yk+BJS(u-γB∗yk)). Then, the iterative sequence {yk+1} can then be rewritten as follows:(21)yk+1=1-λkyk+λkQyk.Because Q is 2/(4-γB2)-averaged, there exists a nonexpansive mapping U such that Q=(1-2/(4-γB2))I+2/(4-γB2)U and Fix(U)=Fix(Q). Therefore, the iterative sequence (21) is equivalent to(22)yk+1=1-λkyk+λk1-24-γB2yk+24-γB2Uyk=1-24-γB2λkyk+24-γB2λkUyk.Notice that 2/(4-γB2)λk∈(0,1) and ∑k=0∞‍λk(1-2/(4-γB2)λk)=+∞. By Theorem 9, we can conclude the following:(1)
{yk} is Fejér monotone with respect to Fix(U).(2)
yk-Uyk→0 as k→∞.(3)
{yk} converges weakly to a fixed point of U.
Taking into account the fact thatFix(U)=Fix(Q), {yk} also converges weakly to a fixed point of Q, and {yk} is Fejér monotone with respect to Fix(Q). Therefore, limk→∞yk-y exists for any y∈Fix(Q). Further, we have that yk-Qyk=2/(4-γB2)(yk-Uyk)→0, as k→∞.
(ii) Lety∈Fix(Q). By Lemma 11 (ii) and Remark 5 (ii), we then have that(23)Qyk-y2=I-J1/γTyk+BJSu-γB∗yk-I-J1/γTy+BJSu-γB∗y2≤yk+BJSu-γB∗yk-y+BJSu-γB∗y2≤yk-y2-1-αα-BJSu-γB∗yk--BJSu-γB∗y2,where α=γB2/2. From Lemma 8, we have that(24)yk+1-y2=1-λkyk-y+λkQyk-y2=1-λkyk-y2+λkQyk-y2-λk1-λkyk-Qyk2.Substituting (23) into (24), we derive(25)yk+1-y2≤yk-y2-1-ααλkBJSu-γB∗yk-BJSu-γB∗y2-λk1-λkyk-Qyk2,which implies that(26)1-ααλkBJSu-γB∗yk-BJSu-γB∗y2≤yk-y2-yk+1-y2+λkλk-1yk-Qyk2.Because inf⁡λk>0, limk→∞yk-y exists, and limk→∞yk-Qyk=0. Taking k→∞ in the above inequality (26), we obtain(27)limk→∞BJSu-γB∗yk-BJSu-γB∗y=0.
By Lemma10, JS+B∗TB(u)=JS(u-γB∗y), where y∈Fix(Q). Then, we have that(28)xk-JSu-γB∗y2=JSu-γB∗yk-JSu-γB∗y2≤u-γB∗yk-u-γB∗y,JSu-γB∗yk-JSu-γB∗y=γB∗y-B∗yk,JSu-γB∗yk-JSu-γB∗y=γy-yk,BJSu-γB∗yk-BJSu-γB∗y≤γy-ykBJSu-γB∗yk-BJSu-γB∗y.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→∞ on the right of (28), we have that xk-JS(u-γB∗y)→0. Thus, we can conclude that the iterative sequence {xk} converges strongly to the resolvent operator of JS+B∗TB(u). This completes the proof.Remark 13.
We observe that(4-γB2)/2>1 for any γ∈(0,2/B2). Taking λk=1 in (20), we obtain a simple iterative algorithm to compute the resolvent operator problem (9). More precisely, for any y0∈G the iterative sequences {xk} and {yk} are generated as follows:(29)xk=JSu-γB∗yk,yk+1=I-J1/γTyk+Bxk,where γ∈(0,2/B2). 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(30)yk+1=I-J1/γTyk+BJSu-γB∗yk.With the help of the relation JγT+γJγ-1T-1∘γ-1I=I, we obtain from (30) that(31)yk+1=1γJγT-1γyk+BJSu-γB∗yk.Let y¯k=γyk. Then, the iterative scheme (31) is equivalent to(32)y¯k+1=JγT-1y¯k+γBJSu-B∗y¯k.Remark 14.
We observe that the resolvent operator ofS+B∗TB in (9) is equivalent to the following monotone inclusion problem:(33)findx∈H,suchthat0∈x-u+Sx+B∗TBx,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(34)findy∈G,suchthat0∈T-1y-BJSu-B∗y.By Lemma 11, we know that -BJS(u-B∗y) is 1/B2-cocoercive. Thus, y is a solution of the dual monotone inclusion (34), and x=JS(u-B∗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 y0∈G and set(35)yk+1=1-λkyk+λkJγT-1yk+γBJSu-B∗yk,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 splittingLetS=0 in Theorem 12. Then, we obtain the following corollary.Corollary 15.
LetT:G→2G be a maximal monotone operator and B∈L(H,G), and let u∈H. For any y0∈G, we define the iterative sequences {xk} and {yk} as follows:(36)xk=u-γB∗yk,yk+1=1-λkyk+λkI-J1/γTyk+Bxk,where γ and λk are the same as in Theorem 12. Then, the following hold:(i)
{yk} converges weakly to a fixed point of Q, where Q(·)=(I-J1/γT)(Bu+(I-γBB∗)(·)).(ii)
Ifinfλk>0, then {xk} converges strongly to the resolvent of JB∗TB(u) and JB∗TB(u)=u-γB∗y, y∈Fix(Q).Remark 16.
Moudafi [26] proposed the following iterative algorithm to solve the resolvent operator JB∗TB(u). Let y0∈G. Then, the iterative sequence {yk} is defined by(37)yk+1=1-λkyk+λkI-J1/γTBu+I-γBB∗yk,where γ∈(0,2/B2) and λk∈(0,1) such that ∑k=0+∞‍λk(1-λk)=+∞. Moudafi proved that JB∗TB(u)=u-γB∗y∞, where y∞ is the weak limit point of {yk}. Furthermore, he presented several applications of the obtained iterative algorithm.
Corollary15 extends some of the results from [26] in two aspects. (i) The range of {λk} is expanded and (ii) we obtain strong convergence of the iterative sequence {xk} (which is not studied in [26]).LetB=I in Theorem 12. Then, we obtain the following corollary.Corollary 17.
LetS:H→2H and T:G→2G be two maximal monotone operators. In addition, let u∈H. For any y0∈G, we define the iterative sequences {xk} and {yk} as follows:(38)xk=JSu-γyk,yk+1=1-λkyk+λkI-J1/γTyk+xk,where γ∈(0,2) and λk∈(0,(4-γ)/2) such that ∑k=0∞‍λk(1-2/(4-γ)λk)=+∞. Then, the following hold:(i)
{yk} converges weakly to a fixed point of Q, where Q(·)=(I-J1/γT)(·+JS(u-γ·)).(ii)
Suppose thatinfλk>0. Then, {xk} converges strongly to the resolvent operator JS+T(u).Remark 18.
The obtained iterative algorithm (38) for computing the resolvent operator JS+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(39)xk=JSu-y¯k,y¯k+1=1-λky¯k+λkJγT-1y¯k+γxk,where y¯k=γyk.
Bauschke and Combettes [21] proposed a Dykstra-like algorithm to compute the resolvent operator JS+T(u). Let x0=u, p0=0, and q0=0. Then, the Dykstra-like algorithm is defined by(40)yk=JTxk+pk,pk+1=xk+pk-yk,xk+1=JSyk+qk,qk+1=yk+qk-xk+1.In addition, Bauschke and Combettes proved that {xk} and {yk} generated by (40) converge strongly to JS+T(u). Following some simple calculation, the Dykstra-like algorithm (40) can be rewritten as follows:(41)xk+1=JSJTxk+pk+qk,pk+1=JT-1xk+pk,qk+1=JS-1JTxk+pk+qk.
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 z0∈H, and define the iterative sequences as follows:(42)yk=J2γ/2+γS22+γzk+γ2+γu,pk=J2γ/2+γT22+γ2yk-zk+γ2+γu,zk+1=zk+λkpk-yk,where γ>0, and {λk}⊆(0,2) such that ∑k=0+∞‍λk(2-λk)=+∞ and inf⁡λk>0. Next, the sequences {yk} and {pk} converge strongly to JS+T(u). In fact, this iterative algorithm (42) is equivalent to applying the Douglas–Rachford splitting algorithm to the monotone inclusion of(43)0∈12x-12u+Sx+12x-12u+Tx.
Comparing (39), (40), and (42), we find that the obtained iterative sequences generated by all algorithms converge strongly to the resolvent operator JS+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 functionf:H→(-∞,+∞] is called proper if domf={x∈H∣f(x)<+∞}≠∅. Γ0(H) denotes the class of proper lower semicontinuous convex functions from H to (-∞,+∞]. The Fenchel-conjugate of f is defined by f∗:u↦supx∈H〈x,u〉-f(x).Problem 19.
Letr∈G and u∈H, and let f:H→(-∞,+∞] and g:G→(-∞,+∞] be proper lower semicontinuous convex functions. Furthermore, let L:H→G be a nonzero bounded linear operator such that r∈core(L(domf)-domg), where the core of a subset C⊆H is defined by coreC={x∈C∣cone(C-x)=H} (see, e.g., Definition 6.9 of [28]). Consider the following convex optimization problem:(44)minx∈H12x-u2+fx+gLx-r.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 Problem19, let y0∈G and set {xk} and {yk} as(45)xk=proxfu-γL∗yk,yk+1=1-λkyk+λkI-prox1/γgyk+Lxk-r,where γ∈(0,2/L2), and λk∈(0,(4-γL2)/2) such that ∑k=0∞‍λk(1-2/(4-γL2)λk)=+∞ and inf⁡λk>0. Then, the iterative sequence {xk} converges strongly to a solution of (44).Proof.
LetS=∂f, T=∂g(·-r), and B=L in (9). Because f∈Γ0(H) and g∈Γ0(G), S and T are maximal monotone operators. Furthermore, note that JS=(I+∂f)-1=proxf and J1/γT=(I+1/γ∂g(·-r))-1=prox1/γg(·-r)+r. Hence, the iterative sequences (45) follow directly from (20). By Theorem 12, we conclude that {xk} converges strongly to a solution of (44). This completes the proof.Remark 21.
The Moreau identity states that, for anyλ>0, proxλf(u)+λprox1/λf∗(1/λu)=u, where f∗ denotes the Fenchel conjugate of f. It follows from this that the iterative sequences (45) can be rewritten as(46)xk=proxfu-γL∗yk,yk+1=1-λkyk+λkγproxγg∗γyk+γLxk-r.Let y¯k=γyk. Then, (46) is reduced to(47)xk=proxfu-L∗y¯k,y¯k+1=1-λky¯k+λkproxγg∗y¯k+γLxk-r.Combettes et al. [27] proposed the following iterative algorithm to solve the optimization problem (44). For any y0∈G, set {xk} and {yk} as(48)xk=proxfu-L∗yk+bk,yk+1=1-λkyk+λkproxγkg∗yk+γkLxk-r+ak,where λk∈[ϵ,1], γk∈(ϵ,2/L2-ϵ), and ϵ∈(0,min(1,1/L2)). Under the conditions that ∑k=0∞‍ak<+∞ and ∑k=0∞‍bk<+∞, the authors proved that {xk} converges strongly to a solution of (44).Remark 22.
(1) Comparing (47) with (48), the range of λk in (47) is clearly larger than that of the iterative sequence (48) introduced by Combettes et al. [27] when the variable step size γk is constant; i.e., γk≡γ. In addition, Theorem 20 also recovers the main results of Proposition 28.16 in [28]. However, we require that λk satisfies ∑k=0∞‍λk(1-2/(4-γL2)λk)=+∞, while Proposition 28.16 in [28] requires that 0<inf⁡λk≤sup⁡λk<(4-γL2)/2. It is obvious that our conditions on λ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 Figure1). Gaussian noise of mean 0 and standard deviation 30 is added to the clear images.Figure 1
Test images. (a) “Barbara”: pixel range[12,246], (b) “Lena”: pixel range [24,245], (c) “Boat”: pixel range [0,255], and (d) “Goldhill”: pixel range [10,235].
(a) (b) (c) (d)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(49)SNR=20log10f2f-x~2,and(50)PSNR=20log10255MNf-x~2,where f is the original image, 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:(51)yk+1-yk2yk2≤ϵ,where ϵ is a given small constant. We tuned the regularization parameter μ and set it to 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:(52)minx∈Rn12x-f22+μxTV+δCx,where δC is the indicator function. It is clear that the optimization problem (52) is a special case of (44). In fact, if f(x)=δC(x), g(Lx)=xTV, and r=0, then the proposed iteration scheme (45) can be employed to solve the constrained TV image denoising problem (4). It is well known that L2≈8, where L is the usual gradient operator of the total variation (see [10, 20]).LetC=Rn. 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∈Rn∣xi≥0}, and the bounded set C as {x∈Rn∣0≤xi≤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)=δC(x) is the orthogonal projection onto the closed convex set C. That is, proxδC(u)=PC(u).


### 5.1. Numerical Results and Discussion

First, we describe the impacts of the parameters for the iterative step sizeγ and relaxation variable λk on the proposed iterative algorithm (45). According to Theorem 20, we can choose γ∈(0,1/4) and λk∈(0,2-4γ). Figures 2, 3, and 4 demonstrate the performance of the proposed iterative algorithm for solving (52) with different choices of γ and λk. The test image was “Barbara,” and the stopping criterion ϵ was set as 10-2.Figure 2
SNR, PSNR, and objective function values against the iteration number for the ROF model. (a)–(c)γ=1/4, (d)–(f) γ=1/8, and (g)–(i) γ=1/16.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Figure 3
SNR, PSNR, and objective function values against the iteration number for the nonnegative ROF model (N-ROF). (a)–(c)γ=1/4, (d)–(f) γ=1/8, and (g)–(i) γ=1/16.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Figure 4
SNR, PSNR, and objective function against the iteration number for the bounded ROF model (B-ROF). (a)–(c)γ=1/4, (d)–(f) γ=1/8, and (g)–(i) γ=1/16.
(a) (b) (c) (d) (e) (f) (g) (h) (i)As shown in Figures2–4, when the iterative step size γ is fixed, a larger relaxation parameter λ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γ and λk in terms of the SNR, PSNR, and number of iterations (Iter) for the test image “Barbara.”
γ λ k ROF N-ROF B-ROF S N R P S N R I t e r S N R PSNR I t e r S N R PSNR I t e r  1/16 0.2 18.7102 24.5975 41 18.7073 24.5946 41 18.7063 24.5936 41 0.4 18.8872 24.7745 34 18.8925 24.7798 35 18.8924 24.7797 35 0.6 18.9328 24.8200 31 18.9331 24.8204 31 18.9332 24.8205 31 0.8 18.9520 24.8393 29 18.9528 24.8401 29 18.9530 24.8403 29 1 18.9633 24.8506 28 18.9644 24.8516 28 18.9646 24.8519 28 1.2 18.9680 24.8553 26 18.9710 24.8582 27 18.9712 24.8585 27 1.5 18.9745 24.8618 25 18.9758 24.8631 25 18.9761 24.8634 25 1.75 18.9775 24.8648 24 18.9789 24.8661 24 18.9791 24.8664 24  1/8 0.2 18.8906 24.7778 34 18.8897 24.7770 34 18.8896 24.7769 34 0.4 18.9529 24.8401 29 18.9536 24.8409 29 18.9538 24.8411 29 0.6 18.9683 24.8556 26 18.9695 24.8568 26 18.9697 24.8570 26 0.8 18.9753 24.8625 24 18.9766 24.8639 24 18.9769 24.8641 24 1 18.9795 24.8667 23 18.9809 24.8691 23 18.9811 24.8684 23 1.2 18.9816 24.8688 22 18.9831 24.8703 22 18.9833 24.8706 22 1.5 18.9829 24.8702 20 18.9844 24.8717 20 18.9847 24.8720 20  1/4 0.2 18.9543 24.8416 28 18.9547 24.8420 28 18.9548 24.8421 28 0.4 18.9762 24.8635 24 18.9774 24.8647 24 18.9777 24.8649 24 0.6 18.9813 24.8686 21 18.9834 24.8707 22 18.9837 24.8710 22 0.8 18.9836 24.8709 20 18.9851 24.8724 20 18.9854 24.8726 20 1 18.9841 24.8714 19 18.9856 24.8729 19 18.9859 24.8732 19Because 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 Table2 that the iterative step size γ has a significant impact on the performance of the algorithm. Furthermore, the experimental results demonstrate that a larger γ always leads to a faster convergence of the iterative algorithm. Thus, we choose γ=1/4 and λk=1 in the following comparison.Next, we focus on investigating the performances of the constrained and unconstrained ROF models for the test images from Figure1. 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 ϵ=10-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 ϵ = 1 0 - 2 ϵ = 1 0 - 3 ϵ = 1 0 - 4 S N R P S N R I t e r S N R PSNR I t e r S N R PSNR I t e r  “Barbara” ROF 18.9841 24.8714 19 18.9800 24.8672 71 18.9784 24.8656 216 N-ROF 18.9856 24.8729 19 18.9816 24.8688 71 18.9800 24.8672 216 B-ROF 18.9859 24.8732 19 18.9818 24.8691 71 18.9802 24.8675 216  “Lena” ROF 22.4595 28.1409 21 22.4552 28.1366 78 22.4522 28.1336 242 N-ROF 22.4605 28.1419 21 22.4563 28.1377 78 22.4534 28.1347 242 B-ROF 22.4616 28.1430 21 22.4573 28.1387 78 22.4544 28.1358 242  “Boat” ROF 21.6226 26.9652 20 21.6180 26.9606 74 21.6162 26.9588 228 N-ROF 21.6326 26.9752 20 21.6278 26.9704 74 21.6260 26.9686 228 B-ROF 21.6335 26.9761 20 21.6287 26.9713 74 21.6269 26.9695 228  “Goldhill” ROF 20.9503 27.3164 20 20.9353 27.3015 73 20.9328 27.2989 223 N-ROF 20.9527 27.3189 20 20.9379 27.3040 73 20.9353 27.3014 223 B-ROF 20.9577 27.3239 20 20.9429 27.3091 73 20.9403 27.3065 223Figure 5
(a) The noisy image of “Barbara.” (b) Denoised image from the ROF model, (c) denoised image from the N-ROF model, and (d) denoised image from the B-ROF model.
(a) (b) (c) (d)Figure 6
(a) The noisy image of “Lena.” (b) Denoised image from the ROF model, (c) denoised image from the N-ROF model, and (d) denoised image from the B-ROF model.
(a) (b) (c) (d)Figure 7
(a) The noisy image of “Boat.” (b) Denoised image from the ROF model, (c) denoised image from the N-ROF model, and (d) denoised image from the B-ROF model.
(a) (b) (c) (d)Figure 8
(a) The noisy image of “Goldhill.” (b) Denoised image from the ROF model, (c) denoised image from the N-ROF model, and (d) denoised image from the B-ROF model.
(a) (b) (c) (d)


## 5.1. Numerical Results and Discussion

First, we describe the impacts of the parameters for the iterative step sizeγ and relaxation variable λk on the proposed iterative algorithm (45). According to Theorem 20, we can choose γ∈(0,1/4) and λk∈(0,2-4γ). Figures 2, 3, and 4 demonstrate the performance of the proposed iterative algorithm for solving (52) with different choices of γ and λk. The test image was “Barbara,” and the stopping criterion ϵ was set as 10-2.Figure 2
SNR, PSNR, and objective function values against the iteration number for the ROF model. (a)–(c)γ=1/4, (d)–(f) γ=1/8, and (g)–(i) γ=1/16.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Figure 3
SNR, PSNR, and objective function values against the iteration number for the nonnegative ROF model (N-ROF). (a)–(c)γ=1/4, (d)–(f) γ=1/8, and (g)–(i) γ=1/16.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Figure 4
SNR, PSNR, and objective function against the iteration number for the bounded ROF model (B-ROF). (a)–(c)γ=1/4, (d)–(f) γ=1/8, and (g)–(i) γ=1/16.
(a) (b) (c) (d) (e) (f) (g) (h) (i)As shown in Figures2–4, when the iterative step size γ is fixed, a larger relaxation parameter λ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γ and λk in terms of the SNR, PSNR, and number of iterations (Iter) for the test image “Barbara.”
γ λ k ROF N-ROF B-ROF S N R P S N R I t e r S N R PSNR I t e r S N R PSNR I t e r  1/16 0.2 18.7102 24.5975 41 18.7073 24.5946 41 18.7063 24.5936 41 0.4 18.8872 24.7745 34 18.8925 24.7798 35 18.8924 24.7797 35 0.6 18.9328 24.8200 31 18.9331 24.8204 31 18.9332 24.8205 31 0.8 18.9520 24.8393 29 18.9528 24.8401 29 18.9530 24.8403 29 1 18.9633 24.8506 28 18.9644 24.8516 28 18.9646 24.8519 28 1.2 18.9680 24.8553 26 18.9710 24.8582 27 18.9712 24.8585 27 1.5 18.9745 24.8618 25 18.9758 24.8631 25 18.9761 24.8634 25 1.75 18.9775 24.8648 24 18.9789 24.8661 24 18.9791 24.8664 24  1/8 0.2 18.8906 24.7778 34 18.8897 24.7770 34 18.8896 24.7769 34 0.4 18.9529 24.8401 29 18.9536 24.8409 29 18.9538 24.8411 29 0.6 18.9683 24.8556 26 18.9695 24.8568 26 18.9697 24.8570 26 0.8 18.9753 24.8625 24 18.9766 24.8639 24 18.9769 24.8641 24 1 18.9795 24.8667 23 18.9809 24.8691 23 18.9811 24.8684 23 1.2 18.9816 24.8688 22 18.9831 24.8703 22 18.9833 24.8706 22 1.5 18.9829 24.8702 20 18.9844 24.8717 20 18.9847 24.8720 20  1/4 0.2 18.9543 24.8416 28 18.9547 24.8420 28 18.9548 24.8421 28 0.4 18.9762 24.8635 24 18.9774 24.8647 24 18.9777 24.8649 24 0.6 18.9813 24.8686 21 18.9834 24.8707 22 18.9837 24.8710 22 0.8 18.9836 24.8709 20 18.9851 24.8724 20 18.9854 24.8726 20 1 18.9841 24.8714 19 18.9856 24.8729 19 18.9859 24.8732 19Because 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 Table2 that the iterative step size γ has a significant impact on the performance of the algorithm. Furthermore, the experimental results demonstrate that a larger γ always leads to a faster convergence of the iterative algorithm. Thus, we choose γ=1/4 and λk=1 in the following comparison.Next, we focus on investigating the performances of the constrained and unconstrained ROF models for the test images from Figure1. 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 ϵ=10-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 ϵ = 1 0 - 2 ϵ = 1 0 - 3 ϵ = 1 0 - 4 S N R P S N R I t e r S N R PSNR I t e r S N R PSNR I t e r  “Barbara” ROF 18.9841 24.8714 19 18.9800 24.8672 71 18.9784 24.8656 216 N-ROF 18.9856 24.8729 19 18.9816 24.8688 71 18.9800 24.8672 216 B-ROF 18.9859 24.8732 19 18.9818 24.8691 71 18.9802 24.8675 216  “Lena” ROF 22.4595 28.1409 21 22.4552 28.1366 78 22.4522 28.1336 242 N-ROF 22.4605 28.1419 21 22.4563 28.1377 78 22.4534 28.1347 242 B-ROF 22.4616 28.1430 21 22.4573 28.1387 78 22.4544 28.1358 242  “Boat” ROF 21.6226 26.9652 20 21.6180 26.9606 74 21.6162 26.9588 228 N-ROF 21.6326 26.9752 20 21.6278 26.9704 74 21.6260 26.9686 228 B-ROF 21.6335 26.9761 20 21.6287 26.9713 74 21.6269 26.9695 228  “Goldhill” ROF 20.9503 27.3164 20 20.9353 27.3015 73 20.9328 27.2989 223 N-ROF 20.9527 27.3189 20 20.9379 27.3040 73 20.9353 27.3014 223 B-ROF 20.9577 27.3239 20 20.9429 27.3091 73 20.9403 27.3065 223Figure 5
(a) The noisy image of “Barbara.” (b) Denoised image from the ROF model, (c) denoised image from the N-ROF model, and (d) denoised image from the B-ROF model.
(a) (b) (c) (d)Figure 6
(a) The noisy image of “Lena.” (b) Denoised image from the ROF model, (c) denoised image from the N-ROF model, and (d) denoised image from the B-ROF model.
(a) (b) (c) (d)Figure 7
(a) The noisy image of “Boat.” (b) Denoised image from the ROF model, (c) denoised image from the N-ROF model, and (d) denoised image from the B-ROF model.
(a) (b) (c) (d)Figure 8
(a) The noisy image of “Goldhill.” (b) Denoised image from the ROF model, (c) denoised image from the N-ROF model, and (d) denoised image from the B-ROF model.
(a) (b) (c) (d)


## 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.



---

*Source: 7376263-2019-05-05.xml*