An improved fast iterative shrinkage thresholding algorithm with an error for image deblurring problem

In this paper, we introduce a new iterative forward-backward splitting method with an error for solving the variational inclusion problem of the sum of two monotone operators in real Hilbert spaces. We suggest and analyze this method under some mild appropriate conditions imposed on the parameters such that another strong convergence theorem for these problem is obtained. We also apply our main result to improve the fast iterative shrinkage thresholding algorithm (IFISTA) with an error for solving the image deblurring problem. Finally, we provide numerical experiments to illustrate the convergence behavior and show the effectiveness of the sequence constructed by the inertial technique to the fast processing with high performance and the fast convergence with good performance of IFISTA.

To solve the variational inclusion problem (1.1) via fixed point theory, we define the mapping J A,B r : H → D(B) as follows: where J B r = (I + rB) -1 is the resolvent operator of B for r > 0. For x ∈ H, we see that which shows that the fixed point set of J A,B r coincides with the solutions set of (A + B) -1 (0).This suggests the following iteration process: x 1 ∈ C and where {r n } ⊂ (0, ∞) and D(B) ⊂ C.This method is called a forward-backward splitting algorithm [16,17].
In applications, we always let A = ∇F and B = ∂G such that F : H → R is a convex and differentiable function and G : H → R ∪ {+∞} is a convex and lower semi-continuous function, where ∇F is the gradient of F with L-Lipschitz continuous and ∂G is the subdifferential of G which defined by Then problem (1.1) is reduced to the following convex minimization problem: Recall that the proximity operator prox G of G is defined for all x ∈ H as follows: For x ∈ H and r > 0, we see that Therefore, Many researchers have proposed and analyzed the iterative shrinkage thresholding algorithms for solving the convex minimization problem (1.2) under a few specific conditions as follows.
In the weak convergence theorems, Lions and Mercier [16] first introduced forwardbackward splitting (FBS) algorithm: where x 1 ∈ H and {λ n } ⊂ (0, 2/L).Later, Moudafi and Oliny [18] introduced the iterative forward-backward splitting (IFBS) algorithm: In our research, we focus attention on the inertial parameter θ n which controls the momentum of x nx n-1 in the fast iterative shrinkage thresholding algorithm (FISTA) of Beck and Teboulle [19] as follows: where y 1 = x 0 ∈ H and t 1 = 1.In FISTA, we observe that y n is known before x n , where the sequence {x n } converges weakly to the solution of the convex minimization problem (1.2).Recently, Hanjing and Suantai [20] introduced the forward-backward modified Walgorithm (FBMWA) as follows: In the same way, Padcharoen and Kuman [21] introduced the forward-backward modified MM-algorithm (FBMMMA) as follows: Other weak convergence theorems of all those algorithms are obtained.
In the strong convergence theorems, Verma and Shukla [22] introduced a new accelerated proximal gradient algorithm (NAGA) as follows: where x 0 , x 1 ∈ H, {α n }, {θ n } ⊂ (0, 1], and {λ n } ⊂ (0, 2/L).They proved that the sequence {x n } of NAGA converges strongly under the condition x n -x n-1 θ n → 0 as n → ∞.How to choose the parameter θ n ?We leave it for the reader to verify.In their proof, we observe that NAGA still holds under conditions α n → 0 and θ n α n x nx n-1 → 0 as n → ∞, and the parameter θ n can be chosen as where {ω n } is a positive sequence such that ω n = o(α n ).Cholamjiak et al. [23] introduced the strong convergence theorem for the inclusion problem (SCTIP) by letting S = I, A = ∇F, and B = ∂G as follows: , where x 0 , x 1 ∈ C and f is a contraction of C into itself, and {α n }, {β n } ⊂ (0, 1), {λ n } ⊂ (0, 2/L), and {θ n } ⊂ [0, θ ] such that θ ∈ [0, 1).They proved that the sequence {x n } of SCTIP converges strongly under the following conditions: Moreover, many researchers have proposed and analyzed the iterative forward-backward scheme with a variable step size, which does not depend on the Lipschitz constant of the operator A = ∇F (see also [24,25]).
In our research, we consider the forward-backward splitting method with an error as follows: x 1 ∈ C and where {λ n } ⊂ (0, ∞), {e n } ⊂ H, D(B) ⊂ C, and J B λ n = (I + λ n B) -1 .We introduce a new iterative forward-backward splitting method with an error for solving the variational inclusion problem (1.1) as follows: where x 0 , x 1 ∈ C and f is a contraction of C into itself, and {α n } ⊂ (0, 1), {λ n } ⊂ (0, 2/L), {e n } ⊂ H, and {θ n } ⊂ [0, θ ] such that θ ∈ [0, 1).Moreover, it can be applied to improve the fast iterative shrinkage thresholding algorithm (IFISTA) with an error for solving the convex minimization problem (1.2) by letting A = ∇F and B = ∂G as follows: , which obtains a self-adaptive scheme with fast convergence properties under some mild conditions when compared to the existing algorithms in the literature.The outline of our research is as follows: in Sect.2, we give some well-known definitions and lemmas which are used in Sect. 3 to prove the strong convergence theorem of IFISTA for solving the variational inclusion problem (1.1), and we also apply its result in Sect. 4 for solving the image deblurring problem, which is a special case of convex minimization problem (1.2); and in Sect.5, we provide numerical experiments to illustrate the fast processing with high performance and the fast convergence with good performance of IFISTA by the inertial technique.

Preliminaries
Let C be a nonempty closed convex subset of a real Hilbert space H.We will use the notation: → to denote the strong convergence, to denote the weak convergence, to denote the weak limit set of {x n }, and Fix(T) = {x : x = Tx} to denote the fixed point set of the mapping T.
Recall that the metric projection P C : H → C is defined as follows: for each x ∈ H, P C x is the unique point in C satisfying Let B be a mapping of H into 2 H .The domain and the range of B are denoted by and v ∈ By.A monotone operator B on H is said to be maximal if its graph is not strictly contained in the graph of any other monotone operator on H.For a maximal monotone operator B on H and r > 0, we define the single-valued resolvent operator J B r : r is firmly nonexpansive and Fix(J B r ) = B -1 (0).We collect together some known lemmas which are the main tools in proving our result.

Lemma 2.1 ([26]) Let C be a nonempty closed convex subset of a real Hilbert space H.
Then:

Lemma 2.2 ([27]
) Let H and K be two real Hilbert spaces, and let T : K → K be a firmly nonexpansive mapping such that

Lemma 2.3 ([27]) Let H be a real Hilbert space and T : H → H be an operator. The following statements are equivalent:
(i) T is firmly nonexpansive, (ii) Tx -Ty 2 ≤ xy, Tx -Ty , ∀x, y ∈ H, (iii) I -T is firmly nonexpansive.

Lemma 2.4 ([28]
) Let C be a nonempty closed convex subset of a real Hilbert space H. Let the mapping A : C → H be α-inverse strongly monotone and r > 0 be a constant.Then we have for all x, y ∈ C. In particular, if 0 < r ≤ 2α, then I -rA is nonexpansive.
where {δ n } is a sequence in (0, 1) and {b n } is a real sequence.Assume that ∞ n=0 c n < ∞.Then the following results hold: Lemma 2.7 ([31]) Assume that {s n } is a sequence of nonnegative real numbers such that where {γ n } is a sequence in (0, 1), {η n } is a sequence of nonnegative real numbers, and {δ n }, {ρ n } are real sequences such that

Theorem 3.1 Let C be a nonempty closed convex subset of a real Hilbert space H. Let A be an α-inverse strongly monotone mapping of H into itself and B be a maximal monotone operator on H such that the domain of B is included in C, and assume that
be the resolvent of B for λ > 0 and f be a k-contraction mapping of C into itself.Let x 0 , x 1 ∈ C and {x n } ⊂ C be a sequence generated by for all n ∈ N, where {α n } ⊂ (0, 1), {λ n } ⊂ (0, 2α), {e n } ⊂ H, and {θ n } ⊂ [0, θ ] such that θ ∈ [0, 1) satisfy the following conditions: Then the sequence {x n } converges strongly to a point x * ∈ (A + B) -1 (0) where x * = P (A+B) -1 (0) f (x * ).
Proof Picking z ∈ (A + B) -1 (0) and fixing n ∈ N, it follows that z = J B λ n (zλ n Az).Firstly, we will show that {x n }, {y n }, and {z n } are bounded.Since therefore, by nonexpansiveness of J B λ n and Iλ n A, we have It follows by the same arguments again that So, by condition (C4) and putting ) ≥ 0 in Lemma 2.6 (i), we conclude that the sequence { x nz } is bounded.That is the sequence {x n } is bounded, and so is {z n }.Moreover, by condition (C4), ∞ n=1 e n < ∞ implies lim n→∞ e n = 0, that is, lim n→∞ e n = 0, it follows that the sequence {y n } is also bounded.
Since P (A+B) -1 (0) f is k-contraction on C, by Banach's contraction principle there exists a unique element x * ∈ C such that x * = P (A+B) -1 (0) f (x * ), that is, x * ∈ (A + B) -1 (0), it follows that x * = J B λ n (x *λ n Ax * ).Now, we will show that x n → x * as n → ∞.On the other hand, we have This implies that It follows by (3.1), Lemma 2.4, and the firm nonexpansiveness of J B λ n that We also have This implies that Hence, by (3.1), (3.2), (3.3), the nonexpansiveness of J B λ n , and Iλ n A, we obtain It follows that Therefore, using conditions (C1), (C3), and (C4), we can check that all those sequences satisfy conditions (i) and (ii) in Lemma 2.7.To complete the proof, we verify that condition (iii) in Lemma 2.7 is satisfied.Let lim i→∞ η n i = 0.Then, by condition (C2), we have and It follows by conditions (C2), (C4) and (3.4) that As {x n } is bounded, so is {x n i }, there exists a subsequence {x n i j } of {x n i } which converges weakly to x ∈ C. Without loss of generality, we can assume that x n i x as i → ∞.On the other hand, by conditions (C1) and (C4), we have It follows that z n i x as i → ∞.Hence, by (3.5) and the demiclosedness at zero in Lemma 2.5, we obtain x ∈ Fix(J B then, by (3.5) and conditions (C1) and (C4), we obtain

It implies that y n i
x as i → ∞.Therefore, by Lemma 2.1(iii), we obtain It follows by conditions (C1), (C3), and (C4) that lim sup i→∞ δ n i ≤ 0. So, by Lemma 2.7, we conclude that x n → x * as n → ∞.This completes the proof.
Remark 3.2 Indeed, the parameter θ n can be chosen as follows: where N ∈ N and {ω n } is a positive sequence such that ω n = o(α n ).
We focus on the image restoration using the fixed point optimization algorithm in Theorem 4.1.The image deblurring problem is in the form where A ∈ R m×n represents a known blurring operator (which is called the point spread function: PSF), b ∈ R m is a known observed blurred (and additive noisy) image, ε ∈ R m is an unknown additive white Gaussian noise, and x ∈ R n is an unknown signal/image to be restored (estimated).Both b and x are formed by stacking the columns of their corresponding two-dimensional images.
In order to solve problem (4.1), we introduce the least absolute shrinkage and selection operator (LASSO) of Tibshirani [34] for solving the following minimization problem: where λ > 0 is a regularization parameter and W : R n → R n represents the orthogonal or tight frame wavelet synthesis, which is a special case of the convex minimization problem (1.2) when such that A T stands the transpose of A, and A is the largest singular value of A (i.e., the square root of the largest eigenvalue of the matrix A T A) or the spectral norm A 2 .
In this image deblurring case using Theorem 4.1, if the blurring operator A is smaller than the observed blurred image b and the restored image x, then it is changed by padPSF in MATLAB to embed its array to the matrix A big ∈ R m×n , and followed by a transformation to the signal matrix A sig ∈ R m×n for calculating the matrix A eig = (a eig ij ) ∈ R m×n of eigenvalues of the signal matrix A sig using the discrete fast Fourier transformation (FFT) or the discrete cosine transformation (DCT).That is, We set m = n and the process of gradient ∇F always maps the signal x to 2 times of the signal A T (Axb), where x, A T , A, and b are in the form of the signal transformation FFT or DCT.That is, , where A T eig = A -1 eig such that A T eig and A -1 eig stand for the transpose and the inverse signal transform of the eigenvalues matrix A eig , respectively.
By [35] and the reference therein, for all u = (u 1 , u 2 , . . ., u m ) T ∈ R m and for each n ∈ N, we have Algorithm 1 Improved fast iterative shrinkage thresholding algorithm (IFISTA).

procedure IFISTA
Choose the initials x 0 , x 1 ∈ R m arbitrarily.Set M is the maximum loops to stop and tol is a prescribed tolerance value.Set the operator A and the mapping f in backing tracks.n ← 0 repeat n ← n + 1 Update the parameters α n , λ n , and θ n , and the error data e n ∈ R m .
x n 2 ≤ tol) 0} for all i = 1, 2, . . ., m.When the process of prox λ n G has been finished, it returns to W -1 (prox λ n G (u)), where W -1 stands for the inverse wavelet synthesis such that W -1 (•) = W T (•) before continuing other processes.That is, In the next section, we present IFISTA in Algorithm 1 to the improved fast iterative shrinkage thresholding algorithm [19] in the same programming techniques [36].

Applications and numerical examples
In this section, we illustrate the performance of IFISTA compared with IFBS, FISTA, FBMWA, FBMMMA, NAGA, and SCTIP for solving the image deblurring problem (4.1) through LASSO problem (4.2) with λ = 10 -4 .We implemented them in MATLAB R2019a to solve and run on personal laptop Intel(R) Core(TM) i5-8250U CPU @1.80 GHz 8 GB RAM.We use the quality measures (it is better if it is larger value) of the image restoration as follows.
Let x, x n ∈ R M×N represent the original image and the estimated image at n th iteration(s), respectively.xx n 2 2 .
(2) For looking at how signal peak is, the measure is the peak signal-to-noise ratio (PSNR) of the images x and x n , which is defined (measured in decibels: dB) by where MAX is the maximum possible pixel value of the m-unit class (m-bit) image such that MAX = 2 m -1 (for instance, MAX = 255 for 8-bits image and MAX = 65535 for 16-bits image), and MSE(x, x n ) is the mean squared error (MSE) of the images x and x n , which is defined by MSE(x, 2 such that the images x and x n are c-multichannel image (for instance, c = 1 for gray or monochrome image, c = 3 for RGB color image, and c = 4 for CMYK color image).
Similarly, this measure is the improvement in signal-to-noise ratio (ISNR) of the images x, x n , and b where the image b ∈ R M×N represents the observed blurred (and additive noisy) c-multichannel image, which is defined (measured in decibels: dB) by .
For comparison, we consider the standard test images downloaded from [37] for Cameraman, Woman, Pirate, and Living room, with each monochrome images consisting of 512 × 512 pixels, which represent the original images x ∈ R 512×512 , and we converted them   to the double class type by im2double(imread('image_name')) in MATLAB, and also its 2D three-stage Haar wavelet transform Wx ∈ R 512×512 as in Fig. 1.The original images went through a Gaussian blur of size 9 × 9 and standard deviation 4 as a point spread function (PSF) which represents the blurring operator A by fspecial or psfGauss in MATLAB, and went through imfilter in MATLAB (computed by mirrorreflecting as the array across-the array border or symmetric) and followed by an additive zero-mean white Gaussian noise with standard deviation 10 -3 , which represents the observed blurred and noisy image b ∈ R 512×512 as in Fig. 2. The PSF A was changed by padPSF in MATLAB to embed its array to the matrix A big ∈ R 512×512 , and it transformed to a signal matrix A sig ∈ R 512×512 for calculating the matrix A eig = (a eig ij ) ∈ R 512×512 of eigenvalues of the signal matrix A sig using the discrete cosine transformation (DCT).That is, In compared algorithms, all parameters have been set to their high performance.For each n ∈ N, we set for FBMWA, NAGA, SCTIP, IFISTA, 1 2 (1 -10 -6 n+1 ) for FBMMMA, for SCTIP, 10 -6 n+1 for FBMWA, for FBMMMA, and by [35], we introduced the best choice types of testing the parameter λ n for the fast convergence as in Table 1 (see also, Tables 1 for NAGA (B2 type), give an evaluation order for those algorithms as in Table 4 as follows: CPU ≤ CPU, SNR ≥ SNR, ISNR ≥ ISNR, and n ≤ n, respectively.We see that all evaluations of IFISTA are satisfied, while only CPU times, SNR, and ISNR evaluations of both IFBS and FISTA are satisfied, where SNR and ISNR of IFBS are greater than FISTA, and then, we conclude that IFISTA, IFBS, and FISTA are in the 1 st , 2 nd , and 3 rd place, respectively, of the top three fast processing with high performance for compared image deblurring as Fig. 5, Fig. 6, Fig. 7, and Fig. 8. From results of each algorithms at first 1000 th iterations as in Table 3, we see that the quantities of SNR and ISNR of all algorithms are vastly different.We give an evaluation order for those algorithms as in Table 4 as follows: CPU ≤ CPU, SNR ≥ SNR, ISNR ≥ ISNR, n ≤ n, and tol ≤ tol, respectively.We see that all evaluations of IFBS, NAGA, and IFISTA     gramming techniques [36] and setting all parameters to their high performance, we obtain the following results.2. For looking at the fast convergence with good performance for compared image deblurring, NAGA, IFBS, and IFISTA are in the 1 st , 2 nd , and 3 rd place, respectively, and all are better than FISTA, FBMWA, FBMMMA, and SCTIP.

Lemma 2 . 5 (
[29] (Demiclosedness principle)) Let C be a nonempty closed convex subset of a real Hilbert space H, and let S : C → C be a nonexpansive mapping with Fix(S) = ∅.If the sequence {x n } ⊂ C converges weakly to x and the sequence {(I -S)x n } converges strongly to y, then (I -S)x = y; in particular, if y = 0, then x ∈ Fix(S).

Lemma 2 . 6 (
[30]) Let {a n } and {c n } be sequences of nonnegative real numbers such that

Figure 1
Figure 1 Original images and their 2D three-stage Haar wavelet transform

Figure 2
Figure 2 Observed blurred and noisy images

Figure 3 Figure 4
Figure 3 The SNR convergence behavior of image deblurring

Figure 5
Figure 5 Top three fast processing with high performance for deblurring of Cameraman

Figure 6
Figure 6 Top three fast processing with high performance for deblurring of Woman

Figure 7
Figure 7 Top three fast processing with high performance for deblurring of Pirate

Figure 8 of 25 Figure 9
Figure 8 Top three fast processing with high performance for deblurring of Living room

Figure 10
Figure 10 Top three fast convergence with good performance for deblurring of Woman

Figure 11
Figure 11 Top three fast convergence with good performance for deblurring of Pirate

Figure 12
Figure 12 Top three fast convergence with good performance for deblurring of Living room

Table 1
The best choice types of testing the parameter λ n for the fast convergence