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

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 amaximalmonotone operator and a bounded linear operator. Based on recentwork, in this paperwe 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.


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 where ∈ is the observed noisy image, ‖ ‖ 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 ∈ the vectorized image ∈ × , where = . 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][3][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][11][12], the primal-dual hybrid gradient algorithm [13][14][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 . That is, ‖ ‖ = ( ). 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 Mathematical Problems in Engineering where ∈ , : → (−∞, +∞] is a proper lower semicontinuous convex function and : → is a linear transformation. The convex optimization problem (2) is equivalent to computing the proximity operator of the function ∘ . Recall that the proximity operator of a proper lower semicontinuous convex function : → (−∞, +∞] is defined as 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: where 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: In general, the constrained convex optimization problem (5) can always be reformulated as an unconstrained convex optimization problem. More precisely, where denotes the indicator function of , which is defined by ( ) = 0 if ∈ and ( ) = +∞ 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 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: where is a proper lower semicontinuous convex function. To achieve this goal, let us recall the definition of the resolvent operator. Let be a real Hilbert space and let : → 2 be a maximal monotone operator. The resolvent of is the single-valued operator = ( + ) −1 . Next, let : → (−∞, +∞] be a proper lower semicontinuous convex function and let denote the subdifferential of , and let = . Thus, the resolvent of coincides with the proximity operator as follows: Under certain qualification conditions, the problem considered in (7) can be solved via the resolvent operator * ∘ ∘ + . Overall, this leads to the problem of computing the resolvent operator of + * . More precisely, let ∈ and where : → 2 and : → 2 are two maximal monotone operators and : → is a bounded linear operator from the Hilbert space to the Hilbert space . Following this, let = and = . Hence, the problem (9) reduces to the convex optimization problem (7), because 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 * , where : → is a bounded linear operator with the adjoint * , : → 2 is a maximal monotone operator, and and are two Hilbert spaces. When = , this algorithm utilizes the fixed-point approach to computing the proximity operator ∘ proposed by Micchelli et al. [20]. It is clear that the resolvent of the operator * is a special case of (9) by allowing = 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 + * (9), where : → 2 and : → 2 are two maximal monotone operators and : → is a bounded linear operator. We assume that the operator + * is maximally monotone. This is true, for example, if 0 ∈ ( (dom ) − dom ) [30], where ( ) represents a relative interior of the set . If = 0, then the problem (9) becomes that of computing the resolvent operator of * . Inspired by the work of Moudafi [26] and Combettes et al. [27], we propose a fixedpoint 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.

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 be a real Hilbert space with inner product ⟨⋅, ⋅⟩ and the associated norm ‖ ⋅ ‖. Let denote the identity operator on . We use ⇀ to indicate that the sequence { } converges weakly to and → to indicate that the sequence { } converges strongly to . Let ( , ) denote all bounded linear operators from the Hilbert space to the Hilbert space . Let ∈ ( , ). Then, the adjoint of is the unique operator * ∈ ( , ) such that ⟨ , ⟩ = ⟨ , * ⟩.
Definition 1 (see [28], (maximal monotone operator)). Let : → 2 be a set-valued operator. is said to be monotone if Moreover, is maximally monotone if its graph is not strictly contained in the graph of any other monotone operator.
It is easy to check that firm nonexpansiveness implies nonexpansiveness. If is firmly nonexpansive, then is 1−cocoercive.
Thus, if is − averaged then is also nonexpansive.
The following lemma provides some useful characterizations between an operator and − .

Mathematical Problems in Engineering
Lemma 6 (see [28]). Let be a nonempty subset of and let : → . Then 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 be a nonempty subset of . Furthermore, let 1 : → be 1 − averaged and let 2 : We also make full use of the following lemma.
Lemma 8 (see [28]). Let ∈ and ∈ . Let ∈ and denote the set of real numbers. Then 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.

Computing the Resolvent Operator (9)
Before presenting our main results, we first introduce some notation. For a fixed ∈ , let > 0 define operators : → and : → as follows: and In addition, the following lemma provides a fixed-point characterization of the resolvent operator (9). Lemma 10. Let : → 2 and : → 2 be two maximal monotone operators. Furthermore, let ∈ ( , ) and ∈ . Then, for a given > 0, if and only if is a fixed point of .
Next, we prove the following lemma, which characterizes an important property of the operator .

Lemma 11. Let
: → 2 be a maximal monotone operator, and let ∈ ( , ). For a given ∈ and > 0, define an operator : → : → − ( − * ). Then Proof. (i) Let 1 , 2 ∈ . Then, we have that The first inequality follows from the fact that the resolvent operator is firmly nonexpansive, and the second is trivial.
Lemma 11 shows that, for any On the other hand, − (1/ ) is firmly nonexpansive and is also 1/2−averaged. According to Now, we are ready to present our main results.
Remark 13. We observe that (4 − ‖ ‖ 2 )/2 > 1 for any ∈ (0, 2/‖ ‖ 2 ). Taking = 1 in (20), we obtain a simple iterative algorithm to compute the resolvent operator problem (9). More precisely, for any 0 ∈ the iterative sequences { } and { } are generated as follows: where ∈ (0,2/‖ ‖ 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 With the help of the relation + −1 −1 ∘ −1 = , we obtain from (30) that Let = . Then, the iterative scheme (31) is equivalent to Remark 14. We observe that the resolvent operator of + * in (9) is equivalent to the following monotone inclusion problem: find ∈ , where , , , and 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][35][36][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 find ∈ , By Lemma 11, we know that − ( − * ) is 1/‖ ‖ 2cocoercive. Thus, is a solution of the dual monotone inclusion (34), and = ( − * ) is a solution of the primal monotone inclusion (33). Next, we apply the relaxed forwardbackward splitting algorithm (see, e.g., Theorem 26.14 of [28]) to solve the dual monotone inclusion (34). More precisely, let 0 ∈ and set 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.
Let = 0 in Theorem 12. Then, we obtain the following corollary.
Corollary 15 extends some of the results from [26] in two aspects. (i) The range of { } is expanded and (ii) we obtain strong convergence of the iterative sequence { } (which is not studied in [26]).
Remark 18. The obtained iterative algorithm (38) for computing the resolvent operator + ( ) is different from those proposed by Bauschke and Combettes [21,23]. In fact, similar to (32), the iterative scheme (38) can be rewritten as where = . Bauschke and Combettes [21] proposed a Dykstra-like algorithm to compute the resolvent operator + ( ). Let 0 = , 0 = 0, and 0 = 0. Then, the Dykstra-like algorithm is defined by = ( + ) , In addition, Bauschke and Combettes proved that { } and { } generated by (40) converge strongly to + ( ). Following some simple calculation, the Dykstra-like algorithm (40) can be rewritten as follows: 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 0 ∈ , and define the iterative sequences as follows: where > 0, and { } ⊆ (0, 2) such that ∑ +∞ =0 (2 − ) = +∞ and inf > 0. Next, the sequences { } and { } converge strongly to + ( ). In fact, this iterative algorithm (42) is equivalent to applying the Douglas-Rachford splitting algorithm to the monotone inclusion of Comparing (39), (40), and (42), we find that the obtained iterative sequences generated by all algorithms converge strongly to the resolvent operator + ( ). 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.

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. be a nonzero bounded linear operator such that ∈ core ( (dom ) − dom ), where the core of a subset ⊆ is defined by core = { ∈ | cone ( − ) = } (see, e.g., Definition 6.9 of [28]). Consider the following convex optimization problem: The minimization problem (44) is equivalent to the proximity operator of ( ) + ( − ). As a direct application of the resolvent operator (9), we obtain the following convergence theorem from Theorem 12.

Remark 22.
(1) Comparing (47) with (48), the range of in (47) is clearly larger than that of the iterative sequence (48) introduced by Combettes et al. [27] when the variable step size is constant; i.e., ≡ . In addition, Theorem 20 also recovers the main results of Proposition 28.16 in [28]. However, we require that satisfies ∑ ∞
(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.

Numerical Experiments
In this section, we present numerical experiments to verify the effectiveness of the proposed iterative algorithms for 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 0 and standard deviation 30 is added to the clear images.
We use the signal-to-noise (SNR) and peak-signal-tonoise (PSNR) to evaluate the quality of the restored images. These are defined by and = 20 log 10 255 √ −̃2 , where is the original image,̃is the restored image, and and are the row and column size of the image , respectively. The iterative algorithm is stopped when the following criterion is satisfied: 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: where is the indicator function. It is clear that the optimization problem (52) is a special case of (44). In fact, if ( ) = ( ), ( ) = ‖ ‖ , and = 0, then the proposed iteration scheme (45) can be employed to solve the constrained TV image denoising problem (4). It is well known that ‖ ‖ 2 ≈ 8, where is the usual gradient operator of the total variation (see [10,20]). Let = . Then, (4) (also (52)) is reduced to the wellknown ROF model (1). We select as both the nonnegative and the bounded sets. Moreover, the nonnegative set is defined as { ∈ | ≥ 0}, and the bounded set as { ∈ | 0 ≤ ≤ 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 ( ) = ( ) is the orthogonal projection onto the closed convex set . That is, ( ) = ( ).

Numerical Results and Discussion.
First, we describe the impacts of the parameters for the iterative step size and relaxation variable on the proposed iterative algorithm (45). According to Theorem 20, we can choose ∈ (0, 1/4) and ∈ (0, 2 − 4 ). Figures 2, 3, and 4 demonstrate the performance of the proposed iterative algorithm for solving (52) with different choices of and . The test image was "Barbara," and the stopping criterion was set as 10 −2 .
As shown in Figures 2-4, when the iterative step size is fixed, a larger relaxation parameter 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 ( ). 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 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 = 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 = 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.

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

Data Availability
The image data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare no conflicts of interest.