# Convergence Analysis of a Modified Forward-Backward Splitting Algorithm for Minimization and Application to Image Recovery

Authors: Kunrada Kankam; Watcharaporn Cholamjiak; Prasit Cholamjiak
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

Many applications in applied sciences and engineering can be considered as the convex minimization problem with the sum of two functions. One of the most popular techniques to solve this problem is the forward-backward algorithm. In this work, we aim to present a new version of splitting algorithms by adapting with Tseng’s extragradient method and using the linesearch technique with inertial conditions. We obtain its convergence result under mild assumptions. Moreover, as applications, we provide numerical experiments to solve image recovery problem. We also compare our algorithm and demonstrate the efficiency to some known algorithms.


---

## Body

## 1. Introduction and Preliminaries

In various fields of applied sciences and engineering such as signal recovery, image restoration and machine learning [1–9] can be formulated as convex minimization problem (CMP) in the term of sum of nonsmooth and smooth functions. Let H be a real Hilbert space. CMP is modeled as follows:
(1)minfk+gk:k∈H,where f:H⟶−∞,+∞ and g:H⟶ℝ are two proper lower semicontinuous convex functions such that f is differentiable on H. For any λ>0, it is known that k∗ is an optimal solution to (1) if
(2)0∈λ∇fk∗+λ∂gk∗⇔I−λ∇fk∗∈I+λ∂gk∗⇔k∗=I+λ∂g−1k∗−λ∇fk∗,where ∇f is the gradient of f is linear function, which is defined by
(3)∇fk,y=f′k,y,where the derivative of f at k in the direction y is f′k,y=fk+ty−fk/t and ∂g· is the classical subdifferential of g which is given by
(4)∂g·=z∈H:gu−g·−z,u−·≥0,∀u∈ℝn.It is known that∂g is maximal monotone and if g is differentiable, then ∂g is the gradient of g denoted by ∇g. This leads to the classical forward-backward splitting algorithm (FBS) [10, 11] which is defined by k0∈H and
(5)kn+1=proxλngkn−λn∇fkn,where λn>0 and proxλng=I+λn∂g−1 is the proximal operator. On the one hand, (5) includes the gradient algorithm kn+1=kn−λn∇fkn, where λn>0 and f is a Lipschitz continuous gradient. Moreover, (5) includes the proximal point algorithm kn+1=proxλngkn, where λn>0 and g is a nondifferentiable function. We know that the proximal operator is single-valued and is characterized by
(6)k−proxλgkλ∈∂gproxλgk,for all k∈H and λ>0. The iteration (5) has been attracted extensively by many researchers. See, for example, [12–18]. One popular method for solving (1) is the modified forward-backward splitting method (MFBS) or Tseng’s extragradient method [19]; MFBS is generated by k0∈H and
(7)kn+1=proxαngkn−αn∇fkn−αn∇fproxαngkn−αn∇fkn−∇fkn,where αn⊂0,+∞ is a real sequence. The convergence rate is well known for the speed of O1/n. Later, various schemes were proposed to improve the convergence and accelerate the method. Among them, Lorenz and Pock [20] have improved the convergence speed of FBS from the standard O1/n to O1/n2.Recently, Beck and Teboulle [21] introduced a fast iterative shrinkage-thresholding algorithm (FISTA-BT) by the following scheme.Algorithm 1.
FISTA-BT algorithm.
Initialization:t0=1 and α=1/L.
Iterative step: letk0=k1∈H and calculate kn+1 as follows:
Step 1. Compute the inertial step:(8)xn=kn+θnkn−kn−1,
wheretn=1+1+4tn−12/2 and θn=tn−1−1/tn.
Step 2. Compute thekn+1 step:
(9)kn+1=proxαgxn−α∇fxn.
Setn≔n+1 and return to Step 1.Without the Lipschitz condition on the gradient of functions, Cruz and Nghia [22] proposed a new version of the forward-backward method (FISTA-CN) based on the linesearch rule.Algorithm 2.
FISTA-CN algorithm.
Initialization:t0=1, σ>0, θ∈0,1, and δ∈0,1/2.
Iterative step: letk0=k1∈H and calculate kn+1 as follows:
Step 1. Compute the inertial step:(10)xn=kn+θnkn−kn−1,yn=PΩxn,
wheretn=1+1+4tn−12/2 and θn=tn−1−1/tn.
Step 2. Compute thekn+1 step:
(11)kn+1=proxαngyn−αn∇fyn,
whereαn=σθmn and mn is the smallest number such that
(12)σθmn∇fproxαngyn−αn∇fyn−∇fyn≤δproxαngyn−αn∇fyn−yn.
Stop criteria ifkn+1=yn, then stop.
Ifkn+1≠yn, then set n≔n+1 and return to Step 1.In 2017, Verma and Shukla [23] introduced the new accelerated proximal gradient algorithm (NAGA) which is generated by the following.Algorithm 3.
NAGA algorithm.
Iterative step: letk0=k1∈H and calculate kn+1 as follows:
Step 1. Compute the inertial step:(13)xn=kn+θnkn−kn−1.
Step 2. Compute(14)yn=1−αnxn+αnproxαngxn−αn∇fxn,kn+1=proxαngyn−αn∇fyn,
whereαn∈0,2/L. Set n≔n+1 and return to Step 1.This work presents a new splitting method called a new modified forward-backward splitting algorithm (NMFBS) for convex minimization problems. Our results extend and improve the corresponding results of Tseng [19] and Cruz and Nghia [22]. The step size defined in this work does not require the Lipschitz condition of the gradient functions. Finally, we also present the numerical experiments of our algorithm for solving image recovery problems and show the comparison of our proposed method to FISTA-BT [21], FISTA-CN [22], and NAGA [23].


## 2. Main Theorem

We assume thatf:H⟶ℝ∪+∞ and g:H⟶ℝ∪+∞ are proper, lower semicontinuous, and convex functions; f is uniformly continuous on bounded sets; and ∇f is bounded on bounded sets. The following is our algorithm.Algorithm 4.
The new modified forward-backward splitting algorithm (NMFBS)
Initialization: givenσ>0, θ>0, δ∈0,1, and θ1>0.
Iterative step: letk0=k1∈H and calculate kn+1 as follows:
Step 1. Compute the inertial step:(15)xn=kn+θnkn−kn−1.
Step 2. Compute:(16)yn=proxαngxn−αn∇fxn+αn∇fxn−∇fproxαngxn−αn∇fxn,
whereαn=σθmn and mn is the smallest number such that
(17)αn2∇fxn−∇fproxαngxn−αn∇fxn2+∇fyn−∇fproxαngyn−αn∇fyn2≤δ2xn−proxαngxn−αn∇fxn2+yn−proxαngyn−αn∇fyn2.
Step 3. Compute thekn+1 step:
(18)kn+1=proxαngyn−αn∇fyn+αn∇fyn−∇fproxαngyn−αn∇fyn.
Setn≔n+1 and return to Step 1.Following the proof as in [24], we can show the following lemma.Lemma 1.
The linesearch (17) has a finite step.Theorem 2.
Suppose thatαn≥α for some α>0, θn≥0, and ∑n=1∞θn<+∞. Then, kn generated by Algorithm 4 converges weakly to a minimizer of f+g.Proof.
Letk∗∈argminf+g, and set pn=proxαngxn−αn∇fxn. Then, we obtain
(19)yn=pn+αn∇fxn−∇fpn.
Moreover, we have(20)xn−pn−αn∇fxn∈αn∂gpn.
Using (19), we see that
(21)αn∇fxn=yn−pn+αn∇fpn.
Combining (20) and (21), we have
(22)xn−yn−αn∇fpn∈αn∂gpn.
Now, setrn=proxαngyn−αn∇fyn. Then, we obtain
(23)kn+1=rn+αn∇fyn−∇frn.
Also, we have(24)yn−kn+1−αn∇frn∈αn∂grn.
Sincek∗∈argminf+g, we obtain −αn∇fk∗∈αn∂k∗. Thus, by (22), (24), and the monotonicity of ∂g, we have
(25)xn−yn−αn∇fpn−∇fk∗,pn−k∗≥0,yn−kn+1−αn∇frn−∇fk∗,rn−k∗≥0.
So, we havexn−yn,pn−k∗≥0 and yn−kn+1,rn−k∗≥0 by the monotonicity of ∇f. Thus, we have
(26)xn−yn,pn−yn+xn−yn,yn−k∗≥0,(27)yn−kn+1,rn−kn+1+yn−kn+1,kn+1−k∗≥0.
We note thatx±y2=x2±2x,y+y2 for all x,y∈H. Using (26), we have
(28)2xn−yn,pn−yn+xn−yn,yn−k∗=xn−k∗2−yn−k∗2+yn−pn2−xn−pn2.
Using (27), we have
(29)2yn−kn+1,rn−kn+1+yn−kn+1,kn+1−k∗=yn−k∗2−kn+1−k∗2+kn+1−rn2−yn−rn2.
From (26), (27), (28), and (29), we obtain
(30)kn+1−k∗2≤xn−k∗2+yn−pn2−xn−pn2+kn+1−rn2−yn−rn2.
Using (17), (19), (23), and (30), we obtain
(31)kn+1−k∗2≤xn−k∗2+pn+αn∇fxn−∇fpn−pn2−xn−pn2+rn+αn∇fyn−∇frn−rn2−yn−rn2=xn−k∗2+αn2∇fxn−∇fpn2+∇fyn−∇frn2−xn−pn2−yn−rn2≤xn−k∗2+δ2xn−pn2+yn−rn2−xn−pn2−yn−rn2=xn−k∗2−1−δ2xn−pn2−1−δ2yn−rn2.
Next, we will show thatlimn⟶∞kn−k∗ exists. From (31), we see that
(32)kn+1−k∗≤xn−k∗=kn−k∗+θnkn−kn−1≤kn−k∗+θnkn−k∗+kn−1−k∗=1+θnkn−k∗+θnkn−1−k∗.
By Lemma 5 in [1], we have
(33)kn+1−k∗≤K·∏j=1n1+2θj,where K=maxk1−k∗,k2−k∗. Since ∑n=1∞θn<+∞, we have kn which is bounded. Thus, ∑n=1∞θnkn−kn−1<+∞. By Lemma 1 in [25] and (32), we have limn⟶∞kn−k∗ that exists. From (31), we see that
(34)kn+1−k∗2≤xn−k∗2−1−δ2xn−pn2−1−δ2yn−rn2=kn+θnkn−kn−1−k∗2−1−δ2xn−pn2−1−δ2yn−rn2=kn−k∗2+2θnkn−k∗kn−kn−1+θn2kn−kn−12−1−δ2xn−pn2−1−δ2yn−rn2.
Notinglimn⟶∞θnkn−kn−1=0, limn⟶∞kn−k∗ exists and δ∈0,1, we have
(35)limn⟶∞xn−pn=0,(36)limn⟶∞yn−rn=0.
Since∇f is uniformly continuous on bounded sets, we have
(37)limn⟶∞∇fxn−∇fpn=0,(38)limn⟶∞∇fyn−∇frn=0.
By definition ofxn, it is easy to see that limn⟶∞xn−kn=0. Then,
(39)pn−kn≤xn−pn+xn−kn⟶0asn⟶∞.
From (35), (36), (37), and (39), we obtain
(40)rn−kn≤rn−yn+yn−pn+pn−kn=rn−yn+pn+αn∇fxn−∇fpn−pn+pn−kn=rn−yn+αn∇fxn−∇fpn+pn−kn⟶0asn⟶∞.
By the boundedness ofkn, we assume that k¯ is a weak limit point of kn; i.e., there is a subsequence kni of kn such that kni⇀k¯. Since limi⟶∞rni−kni=0, we also obtain rni⇀k¯ as i⟶∞. Using (6), we obtain
(41)yni−rniαni−∇fyni=yni−proxαnigyni−αni∇fyni−αni∇fyniαni∈∂gyni−αni∇fyni.
It follows that(42)yni−rniαni+∇frni−∇fyni∈∇frni+∂grni⊆∂f+grni.
By passingi⟶∞ and using (36) and (38), we have 0∈∇fk¯+∂gk¯ by Fact 2.2 in [22]. Hence, by Theorem 5.5 in [26], we can conclude that kn converges weakly to a point in argminf+g. We thus complete the proof.Remark 3.
The condition thatαn≥α>0 for some α can be dropedd in case ∇f is Lipschitz continuous on H since it is bounded below from 0 (see Proposition 4.4,4.11 [22]).Remark 4.
In the main theorem, we use the linesearch technique to calculate our step size at each iteration unlike the result of [4, 5, 16, 17]. It is worth mentioning here also that choice of the step size in our algorithm does not depend on the Lipschitz condition of the gradient function. Our proposed algorithms can be applied in image recovery which are more applicable than those of [4, 5, 16, 17].


## 3. Numerical Experiments

Medical imaging plays a crucial role in modern medicine and image data which are found in various clinical specialties, for routine diagnostics in X-ray imaging, monitoring intraoperative progress during surgical procedures and guidance and diagnosis in ailing. In practice, the degradations are unavoidable because the medical imaging systems limit the intensity of the incident radiation to protect the patient’s health. So how to improve image quality is a good choice for medical analysis. Image processing mainly consists of image deblurring, image denoising, and image inpainting which is a branch that usually can be employed optimization techniques to solve it.The image restoration problem can be explained as follows:(43)b=Ak+w,where b∈ℝm×1 is the observed image, A∈ℝm×n is the blurring matrix, k∈ℝn×1 is an original image, and w is additive noise. To solve problem (43), we aim to approximate the original image by transforming (43) to the following LASSO problem [27]:
(44)mink12b−Ak22+λk1,where ·1 is ℓ1-norm. In general, (44) can be formulated in a general form by estimating the minimizer of sum of two functions when fk=1/2b−Ak22 and gk=λk1. We next present our algorithm (NMFBS) for LASSO problem with λ=10−7 and also compare its efficiency with FISTA-BT [21], FISTA-CN [22], and NAGA [23]. All computational experiments were written in Matlab 2020b and performed on a 64-bit MacBook Pro Chip Apple M1 and 8 GB of RAM.Letk be the original images size 448×2993 and 386×608×3, respectively. These are shown in Figure 1. To measure the quality of restored images, we use the peak signal-to-noise ratio (PSNR) in decibel (dB) [28] and the structural similarity index metric (SSIM) [29]. The iteration numbers for all algorithms is 1200th.Figure 1
(a) The original image size448×2993 and (b) 386×608×3, respectively.
(a)(b)All parameters are chosen as in Table1. The initial points k0=k1 are vectors of ones with the size of original images for all algorithms. The blurred images are shown in Figures 2–4. The parameter θn of FISTA-BT, FISTA-CN, and NAGA is defined as in Algorithm 1.Table 1
Chosen parameters of each algorithm.
AlgorithmsParametersγθδσt1FISTA-BT1/A———1FISTA-CN—0.90.40.21NAGA1/A———1NMFBS—0.90.90.21Figure 2
The blurred RGB images of Fig(a) and Fig(b) for out-of-focus blur matrices with radiusr=7.
(a)(b)Figure 3
The blurred RGB image of Fig(a) and Fig(b) for Gaussian blur of the filter size7×7 with standard deviation σ=15.
(a)(b)Figure 4
The blurred RGB image of Fig(a) and Fig(b) for motion blur specified with the motion length of 15 pixels and motion orientationθ=15.
(a)(b)The numerical results are reported in Table2 and Figures 5–8.Table 2
The results of deblurred images for each algorithm.
Blur typesMethodsFig(a)Fig(b)PSNRSSIMPSNRSSIMOut-of-focus blur (disk) with radiusr=7 (see Figure 2)FISTA-BT34.16710.895144.60860.9855FISTA-CN36.39110.932745.89920.9888NAGA35.61760.922245.48530.9878NMFBS40.19680.972847.38280.9916Gaussian blur of the filter size5×5 with standard deviation σ=5 (see Figure 3)FISTA-BT34.84320.963145.85460.9906FISTA-CN36.51380.973547.03150.9926NAGA35.95270.970346.65340.9920NMFBS38.98910.983648.48200.9946Motion blur specified with the motion length of45 pixels and motion orientation θ=45 (see Figure 4)FISTA-BT40.54800.954848.60890.9940FISTA-CN44.03010.978250.02190.9955NAGA42.75660.970949.64780.9952NMFBS50.54680.994951.37040.9967Figure 5
The blurred image of Fig(a) by out of focussing and the restored images for FISTA-BT (PSNR: 34.1671, SSIM: 0.89951), FISTA-CN (PSNR: 36.3911, SSIM: 0.9327), NAGA (PSNR: 35.6176 SSIM: 0.9222), and NMFBS (PSNR: 40.1968, SSIM: 0.9728), respectively.Figure 6
Graphs of PSNR and SSIM for Fig(a) by out of focussing, respectively.Figure 7
The blurred image of Fig(b) by motion blurring and the restored images for FISTA-BT (PSNR: 48.6089, SSIM: 0.9940), FISTA-CN (PSNR: 50.0219, SSIM: 0.9955), NAGA (PSNR: 49.6478 SSIM: 0.9952), and NMFBS (PSNR: 51.3704, SSIM: 0.9967), respectively.Figure 8
Graphs of PSNR and SSIM for Fig(b) by motion blurring, respectively.From Table2, we see that numerical experiments of NMFBS are better than those of FISTA-BT, FISTA-CN, and NAGA in terms of PSNR and SSIM for all blur types.We next provide some experiments of the recovered images for two cases to illustrate the convergence behavior of all algorithms in comparison. We plot the number of iterations versus PSNR and SSIM in Figures6 and 8.


## 4. Conclusion

We have introduced the modified forward-backward algorithm for solving the convex minimization problem of the sum of two functions in a real Hilbert space. The proposed algorithm does not need to compute the Lipschitz constant of the gradient of functions. We have proved that the sequence generated by the algorithm weakly converges to a minimizer under some mild conditions. Our result can be applied effectively to solve image recovery as shown in numerical experiments. The comparative experiments showed that the proposed algorithm has a better efficiency than FISTA-BT [21], FISTA-CN [22], and NAGA [23] in terms of PSNR and SSIM for all blur types.



---

*Source: 3455998-2022-10-06.xml*