Preconditioning Technique for an Image Deblurring Problem with the Total Fractional-Order Variation Model

: Total fractional-order variation (TFOV) in image deblurring problems can reduce/remove the staircase problems observed with the image deblurring technique by using the standard total variation (TV) model. However, the discretization of the Euler–Lagrange equations associated with the TFOV model generates a saddle point system of equations where the coefﬁcient matrix of this system is dense and ill conditioned (it has a huge condition number). The ill-conditioned property leads to slowing of the convergence of any iterative method, such as Krylov subspace methods. One treatment for the slowness property is to apply the preconditioning technique. In this paper, we propose a block triangular preconditioner because we know that using the exact triangular preconditioner leads to a preconditioned matrix with exactly two distinct eigenvalues. This means that we need at most two iterations to converge to the exact solution. However, we cannot use the exact preconditioner because the Shur complement of our system is of the form S = K ∗ K + λ L α which is a huge and dense matrix. The ﬁrst matrix, K ∗ K , comes from the blurred operator, while the second one is from the TFOV regularization model. To overcome this difﬁculty, we propose two preconditioners based on the circulant and standard TV matrices. In our algorithm, we use the ﬂexible preconditioned GMRES method for the outer iterations, the preconditioned conjugate gradient (PCG) method for the inner iterations, and the ﬁxed point iteration (FPI) method to handle the nonlinearity. Fast convergence was found in the numerical results by using the proposed preconditioners.


Introduction
Although TV regularization is a commonly employed technique in image deblurring problems [1][2][3][4], one significant drawback is the appearance of the "staircase effect", wherein edges are depicted as a sequence of steps rather than smooth transitions.This phenomenon arises because TV regularization encourages the creation of piecewise constant regions, which leads to the appearance of blocks around edges rather than accurately capturing their continuous nature.As a result, researchers are actively investigating and advancing alternative regularization techniques and algorithms to remove or reduce these "staircase effects" and enhance the overall quality of image deblurring methods.An alternative regularization approach is the TFOV model [5][6][7][8].TFOV regularization presents a robust method for enhancing image deblurring, offering a combination of benefits such as edge preservation, flexibility, noise resilience, and reduction in staircase effects.Its effectiveness has been substantiated in numerous studies, significantly contributing to the progression of image deblurring techniques.However, when it comes to discretizing the Euler-Lagrange equations of the TFOV-based model, a substantial nonlinear and ill-conditioned system emerges.Efficiently solving such systems poses a considerable challenge for numerical methods, even with the application of potent numerical algorithms like Krylov subspace methods, such as the generalized minimal residual (GMRES) and conjugate gradient (CG).These methods tend to exhibit slow convergence in this context.One potential remedy to address this slow convergence is the use of preconditioning techniques.Preconditioning is a technique used to transform a linear system of the form Ax = b into another system to improve the spectral properties of the system matrix.A preconditioner is a matrix P that is easy to invert and the preconditioned matrix P −1 A shows good clustering behavior for the eigenvalues.This is because rapid convergence is often associated with a clustered spectrum of P −1 A. In the preconditioning technique, we solve P −1 Ax = P −1 b instead of solving the original Ax = b because the new system P −1 Ax = P −1 b converges rapidly when we use a suitable preconditioner.To apply the preconditioner matrix P within a Krylov subspace method, we need to compute the multiplication of a matrix by a vector of the form z = Pr at each iteration.Hence, evaluating this product must be cheap.In the literature, several preconditioners are developed in [9] for a special linear system, such as block preconditioners and constraint preconditioners.For diagonal preconditioners, we can refer to Silvester and Wathen [10] and Wathen and Silvester [11].For the block triangular preconditioners, we can refer to Bramble and Pasciak [12] and [13][14][15][16], as well as the references therein.For constraint preconditioners, see, for example, Axelsson and Neytcheva [17].Other preconditioners based on Hermitian/skew-Hermitian splitting are studied in [18][19][20].Recently, several new preconditioners for Krylov subspace methods have been introduced.For example, Cao et al. [21] derived two block triangular Schur complement preconditioners from a splitting of the (1, 1)-block of the two-by-two block matrix.Chen and Ma [22] proposed a generalized shift-splitting preconditioner for saddle point problems with a symmetric positive definite (1, 1)-block.Salkuyeh et al. [23] proposed a modified generalized shift-splitting preconditioner where the (1, 1)-block is symmetric positive definite and the (2, 2)-block matrix is symmetric positive semidefinite (not zero).Very recently, block diagonal and block triangular splitting preconditioners were studded by Beik et al. [24], and the authors introduced new variants of the splitting preconditioners and obtained new results for the convergence of the associated stationary iterations and new bounds for the eigenvalues of the corresponding preconditioned matrices.Moreover, they considered inexact versions as preconditioners for flexible Krylov subspace methods.A good survey of preconditioning techniques for general linear systems can be found in [9, 25,26].
In our paper, we consider the following two-by-two block nonlinear system of equations: This system is obtained by discretizing the Euler-Lagrange equations associated with TFOV in image deblurring problems, and the coefficient matrix of this system has a size of 2n by 2n, where n := N 2 and N is the number of pixels in the image.The coefficient matrix of this system is non-symmetric, ill conditioned, dense, and huge.These properties complicate the development of an efficient numerical algorithm.We know that using direct methods for solving Equation (1) requires O(N 3 ) and, hence, they are not applicable here.For this system, iterative methods, like Krylov subspace methods, are applicable.However, their convergence is too slow because they are sensitive to the condition numbers.Hence, preconditioning is needed to accelerate the convergence of the Krylov subspace methods.In this paper, we propose two block triangular preconditioners for Equation (1).In the literature, it has been shown that such preconditioners are among the most effective for solving problems of the saddle point type.Moreover, it is known that using the exact triangular preconditioner leads to a preconditioned matrix with exactly two distinct eigenvalues [25].This means that we need at most two iterations to converge to the exact solution.Since the coefficient matrix A is not symmetric, the suitable outer iterative method is the GMRES method [27], and since the (2, 2)-block in the matrix A is symmetric positive definite, the suitable inner iterative method is the CG method.However, using the GMRES Krylov subspace method as a preconditioner within a different Krylov subspace method (the CG method) may lead to a changing preconditioner.In such cases, the preconditioner matrix changes from step to step.For this reason, we use flexible GMRES (FGMRES) instead of GMRES [27].The flexibility here means that FGMRES is designed to be flexible in terms of the choice of the inner Krylov subspace method and the choice of the preconditioner.This flexibility allows FGMRES to adapt to different Krylov subspace methods.FGMRES can be restarted after a certain number of iterations to control the memory usage and computational cost, especially when solving multiple linear systems with different righthand sides.The main contributions of this work are follows:

•
We propose two block triangular preconditioners and study the bounds of the eigenvalues of the preconditioned matrices.In addition, we demonstrate the effectiveness of our algorithm in the numerical results by starting with the fixed point iteration (FPI) Method as in [28] to linearize the nonlinear primal system . ., then we use the preconditioned conjugate gradient (PCG) method [29] for the inner iterations.After that, we use FGMRES method for the outer iterations.We illustrate the performance of our approach by calculating the peak signal-to-noise ratio (PSNR), CPU-time, residuals and the number of iterations.Finally, we calculate the PSNR for different values of the order of the fractional derivative, α, to show the impact of using the TFOV model.
The remainder of this paper is organized as follows: Section 2 presents the mathematical model of the image deblurring problem, different regularization models, three definitions of the fractional derivative, and the Euler-Lagrange equations associated with the TFOV model.System ( 1) is obtained at the end of this section.Section 3 presents all theoretical contributions of this paper.Section 4 reports some numerical results that show the efficiency of our preconditioners.Section 5 briefly states our conclusions.

Problem Setup
We know that blurring and noise affect the quality of the received image.To deblur an image, we need a mathematical model of how it was blurred.The recorded image z and the true (exact) image u are related by the equation where K denotes the following blurring operator: with translation-invariant kernel, k(x, x ) = k(x − x ), known as the point spread function (PSF).ε is the additive noise function.Ω will denote a square in R 2 on which the image intensity is defined.When K is the identity operator, the problem (2) becomes image de-noising.In this paper, we focus on de-blurring problem.The PSF function must be known.However, if it is unknown, another technique named blind deconvolution can be used [30].The operator K is compact, so the problem (2) is ill-posed [31], and then the resulting matrix systems from the discretization of this problem are highly ill-conditioned.In this case, directly solving this problem is difficult.The most popular approach to obtain a well-posed problem is to add a regularization term.Different regularization terms are used in the literature, for example: 1.
Tikhonov regularization [32] is used to stabilize the problem (2) and also called as penalized least squares.In this case, the problem is then to find a u that minimize the functional with a small positive parameter λ called the regularization parameter that controls the trade-off between the data fitting term (the first term) and the regularization term (the second term).• denotes the norm in L 2 (Ω).The functional J has to be known.Common choices for the functional J are the above functional gives what is known as Tikhonov regularization with the identity, and where | • | denotes Euclidean norm, and ∇ = ∂ ∂x , ∂ ∂y .When u is discontinuous, the functional in (5) often induces either oscillations or ringing.However, in the functional (6), we need to assume that u is a smooth function.Although, this model is easy to use and simple to calculate, it cannot preserve image edges.Hence, both the above choices are unsuitable for image processing applications when we need to recover sharp contrasts modeled by discontinuities in u [28].

2.
Total Variation (TV): One of the most commonly used regularization models is the TV.It was introduced for the first time [33] in edge-preserving image denoising by Rudin, Osher and Fatemi (ROF) and it has improved in recent years for image de-noising, de-blurring, in-painting, blind de-convolution, and processing [1][2][3][4][34][35][36][37][38][39].When using the TV model, the problem is then to find a u that minimizes the functional where Note that, we do not require the continuity of u.Hence, (8) is a good regularization in image processing.However, the Euclidean norm, | ∇u |, is not differentiable at zero.Common modification is to add a small positive parameter β.The resulting is in the modified functional: The well-posedness of the above minimization problem (7) with the functional given in (9) is studied and analyzed in the literature, such as in [1].The success of using TV regularization is that TV gives a balance between the ability to describe piecewise smooth images and the complexity of the resulting algorithms.Moreover, the TV regularization performs very well for removing noise/blur while preserving edges.Despite the good contributions of the TV regularization mentioned above, it favors a piecewise constant solution in the bounded variation (BV) space which often leads to the staircase effect.Thus, stair casing remains one of the drawbacks of the TV regularization.To remove the stair case effects, two modifications to the TV regularization model have been used in the literature.The first approach is to higher the order of the derivatives in the TV regularization term, such as the mean curvature or a nonlinear combination of the first and second derivatives [40][41][42][43][44][45].These modifications remove/reduce the staircase effects and they are effective but they are computationally expensive due to the increasing the order of the derivatives or due to the nonlinearity terms.The second approach is to use the fractional-order derivatives in the TV regularization terms as shown in [46,47].

Fractional-Order Derivative in Image Deblurring
The most important advantage of using fractional differential equations is their nonlocal property.The integer order differential operator is a local operator but the fractional order differential operator is nonlocal.This means that the next state of a system depends not only on its current state but also upon all of its historical states.This is more realistic and it is one reason why fractional calculus has become more and more popular.
In image deblurring problems, the blurring is considered nonlocal in some cases and local in others depending on the cause of the blur.For example, if a body is moving while the background is stationary, then the blur is local and in case the camera is moving then the blur is nonlocal.The blurring operator is a convolution operator that depends on the definition of the kernel functions.In the case of a moving camera, the blurring operator involves each pixel in the image, which means that the blurring process is a nonlocal.The bulrring is nonlocal, in this case, so it is appropriate to choose a regularization operator with the same nonlocal property.This property is available in operators that contain fractional derivatives.Comparative studies have shown that fractional-order differentials are more reliable than the integer-order differentials for enhancing edges in image processing.A similar trend has been observed for the texture and area-preserving properties.Therefore, images processed by fractional differentials are clearer and have higher contrast [48].This approach is widely used in image processing [5][6][7][8]49,50].These works have shown that the fractional-order derivative performs well in achieving a satisfactory compromise between avoiding staircasing and preserving important fine-scale features such as edges and textures.In this paper, we compare the results of the usual TV model with the TFOV model for two image deblurring problems.From Figures 1-6, we can see that the TFOV shows better edge enhancement results than TV in some regions where we can observe that the PSNR at α = 1 is lower than PSNR at α > 1.
Example 1.We used the exact Golden House image plotted in Figure 9, the deblurred Golden House image using the TV-model plotted in Figure 27, and the deblurred Golden House image using the TFOV-model plotted in Figure 29.We took a vertical line almost in the middle of these images and plotted the results of the cross sections in Figure 1.In Figure 1, we highlighted three corners (boxes).From Figures 1-6, we can clearly see that the TFOV-based image deblurring results are smoother than the TV-based image deblurring.Additionally, in each corner TV based image deblurring creates a higher error than TFOV.This shows that TFOV-based image deblurring is better and edge-preserving.The numerical results in the above examples reflect good performance, and motivating us to use the TFOV model in our preconditioners.

The TFOV-Model
Let BV α (Ω) denotes the space of functions of α-bounded variation on Ω defined by The parameter α represents the order of the fractional derivatives, the fractional-order total variation of u, TV α is defined by ∂y denote the fractional-order derivative along the x and y directions respectively.The space T denotes the space of special test functions ) denote the space of α-order continuously differentiable functions.Hence, when the TFOV-model is used, the problem is then to find a u ∈ BV α (Ω) ∩ L 2 (Ω) that minimizes the functional where J α TV β is called the modified total fractional-order variation and defined by where |∇ α u| 2 = (D α x u) 2 + (D α y u) 2 where D α x , D α y are the fractional derivative operators along the x and y directions respectively.Existence and uniqueness of a minimizer to the above problem (10) with the functional (11) are studied and analyzed in the literature [8,51].

Fractional-Order Derivatives
Several definitions have been proposed for fractional-order derivatives [52][53][54].We shall present some of them below.For a systematic presentation of mathematics, a fractionalorder derivative is denoted as function operator D α [a,x] , where a and x are the bounds of the integrals, and α is the order of the fractional derivative such that 0 < n − 1 < α < n where n = [α] + 1 and [•] is the greatest integer function.

1.
Riemann-Liouville (RL) definitions: The left-and right-sided RL derivatives of order α of a function f (x) are given as follows: and where Γ(•) is the gamma function, defined by

Grünwald-Letnikov (GL) definitions:
The left-and right-sided GL derivatives are defined by and where 3.

Caputo (C) definitions:
The left-and right-sided Caputo derivatives are defined by where f (n) denotes the nth-order derivative of function f (x).

Euler-Lagrange Equations
In this subsection, we present the Euler-Lagrange equations associated with the TFOV in image de-blurring problem.
Theorem 1.If α ∈ (1, 2), the Euler-Lagrange equations for the functional given in (10) are: where K * is the adjoint operator of the integral operator K and the nonlinear deferential operator L α (u) is given by: Proof.Let ν ∈ W α 1 (Ω) be a function.Then for u ∈ W α 1 (Ω) ⊂ BV α (Ω), the first order Gateaux derivative of the functional F α (u) of (10) in the direction of ν is where G 1 (u) = 1 2 Ω (Ku − z)dx and G 2 (u) = λJ α TVβ (u).By using the Taylor series in the direction of t, we have where where we know that n = 2 for 1 < α < 2. Case-I: If u(x) Then it suffices to take ν ∈ C 1 0 (Ω, R), this implies Hence (22) reduces to (19).Case-II:

So boundary terms in (23) can only become zero if
This completes the proof.
Note that ( 19) is a nonlinear integro-differential equation of elliptic type.Equation ( 19) can be expressed as a nonlinear first order system [55]: with the dual, or flux, variable We apply the Galerkin method to ( 24)-( 25) together with the midpoint quadrature for the integral term and cell-centered finite difference method for the derivative part.

Discretization of the Fractional Derivative
First, we divide the square domain Ω = (0, 1) × (0, 1) into N 2 equal squares (cells) where N denotes the number of equispaced partitions in the x or y directions.Then, we follow the same discretization in [8,51].We define a spatial partition (x k , y l ) (for all k, l = 0, 1, . . ., N + 1) of image domain Ω.Assume u has a zero Dirichlet boundary condition, we consider the discretization of the α-order fractional derivative at the inner point (x k , y l ) (for all k, l = 0, 1, . . ., N) on Ω along the x-direction by using the shifted Gr ünwald approximation approach [56,57] where f l s = f s,l and ω α j = (−1) j α j j = 0, 1, . . ., N and ω α 0 = 1, ω α j = (1 − 1+α j )ω α j−1 for j > 0. Observe from (27) that the first order estimate of the α-order fractional derivative After incorporating the zero boundary condition in the matrix approximation of fractional derivative, all N equations of fractional derivatives along the x direction in (27) can be written simultaneously in the matrix form From the definition of fractional-order derivative (27), for any 1 < α < 2, the coefficients ω k α have the following properties: (1) Hence by the Gershgorin circle theorem, we can derive that matrix B α N is a symmetric and negative definite Toeplitz matrix (i.e., −B α N is a positive definite Toeplitz matrix).Let U ∈ R N×N denote the solution matrix at all nodes (khx; lhy), k, l = 1, . . ., N corresponding to x-direction and y-direction spatial discretization nodes.Denote by u ∈ R N 2 ×1 , the ordered solution vector of U. The direct and discrete analogue of differentiation of arbitrary α order derivative is Similarly, all values of α-th order y-direction derivative of u(x; y) at these nodes are approximated by where u = u 11 , u 12 , . . . ,u NN and ⊗ denotes the Kronecker product.For more details in the discretization, we refer to [54,58].Now, using the cell center finite difference Method (CCFDM), the fractional discretization shown above, and using the fact that [(−1) n ∇ α •] is the adjoint operator of the operator ∇ α , then ( 24)-( 25) leads to the following system where N F is the number of Fixed-Point Iterations used to linearize the nonlinear term in the square root in (26).The matrix K h is obtained form using the midpoint quadrature for the integral operator as follows: with entries [K h U] ij,lm = h 2 k(x i − x j , y l − y m ).With using the lexicographical order, K h is a block Toeplitz with Toeplitz block (BTTB) matrix.The need for BTTB property will be discussed later in the paper.The discrete scheme of the matrix L α h U is given by: where • is the point wise multiplication and m is the m − th Fixed-Point Iteration.U is an N × N-size reshaped matrix of the vector u and the matrices D 1 (U m ) and D 2 (U m ) are the diagonal of the Hadamard inverses of the non-zero matrices B x α (U m ) and B y α (U m ) respectively.

Difficulties in TFOV-Model Compared to TV-Model
In this subsection, we compare the TFOV-system (1): and the following TV-system: In the TFOV system (1), the fractional matrix L α h is obtained from discretizing a fractional deferential operator and it is dense.The density property leads to an expensive matrixvector multiplication.In this case, the coefficient matrix A α in the system (1) contains three dense submatrices, while in TV system (34), the non-fractional matrix L h is obtained from discretizing a non-fractional deferential operator and it is a sparse matrix, then the coefficient matrix A in the system (34) contains only two dense submatrices.Further, the Schur complement matrix associated with (1) is a sum of two dense matrices while the Schur complement matrix associated with (34) is a sum of one dense matrix and one sparse matrix.

Preconditioning Technique
In the literature, it has been shown that block triangular preconditioners are among the most effective preconditioners for solving saddle point problems.In this paper, we develop two block triangular preconditioners for solving (1).First, we present our main preconditioner matrix [12] and its inverse: where S = K * K + λL α is the Schur complement matrix.We notice that the Schur complement matrix contains the product (K * K) which is not a BTTB matrix.We know that a BTTB matrix-vector product computation cost O(N log N) but using a BCCB extension.Since this extension is not an easy task in some cases, the idea of using a circulant matrix as a preconditioner for a Toeplitz matrix is needed.This idea was first proposed by Strang [59] and Olkin [60] and extended by others to block Toeplitz systems for example Chan et al. [61].Many researchers use Toeplitz preconditioners and block Toeplitz preconditioners for Toeplitz systems.For instance, Chan et al. [62], and Lin and Fu-Rong [63].Band Toeplitz preconditioner and band BTTB preconditioner are proposed by Chan and Raymond [64] and Serra and Stefano [65].In Lin et al.
[66], BTTB preconditioners for BTTB systems are discussed.Several kinds of circulant preconditioners have been proposed to be good preconditioners, see for instance [59,62,[67][68][69].Several kinds of circulant preconditioners have been proposed and proven to be good preconditioners.Therefore, the PCG methods with circulant preconditioners converge very fast when they are used to solve Toeplitz systems.Motivated by these papers, we propose the following two block triangular preconditioners: where S 1 = (I + λL TV ) and S 2 = (C * C + λL TV ).Where I is the denoising operator, the identity matrix is a circulant matrix, and L TV comes from discretaizing the TV model (α = 1) which is a sparse matrix and C is the Strang circulant approximation of the matrix K [59].These circulant approximations are very important to allow us to use the FFT and the convolution theorem.We know that all circulant matrices can be diagonalized by the Fourier matrix, see [70].Also using FFT and the convolution theorem will reduce the cost of the computation from O(N 2 ) into O(N log N).Moreover, all that is needed for computation is the first column or the first row of the circulant matrix, which decreases the amount of required storage.This reduction in the computations and storage leads to efficient solvers for our problem (1).

Preconditioned GMRES Algorithm
In this section, we give a detailed algorithms for using our preconditioner P (P 1 and P 2 ).In Algorithm 1, GMRES method is used to solve the linear system (1).
In Algorithm 1, in Steps 3 and 7, we need to solve a matrix times a vector of the form where S = S 1 or S = S 2 .To do the above multiplications, we use the conjugate gradients method as in Algorithms 2 and 3: 10: 11: end for 12: and e 1 = (1, 0, . . ., 0) T 14: end for 15: Solve for x 2 in the system −S 1 x 2 = b 2 using conjugate gradient method.

Eigenvalues Estimates
In this subsection, we need to study the eigenvalues of the exact preconditioned matrix P −1 A. Since P −1 A and AP −1 are similar matrices, they have the same eigenvalues.Hence we study the eigenvalues of the matrix AP −1 .
Theorem 2. If the linear system (1) is left preconditioned by the matrix P, then the preconditioned matrix is and its minimal polynomial is (ν − 1)(ν + 1) where ν is the eigenvalue of the matrix AP −1 .
Proof.Since AP −1 and P −1 A are similar, it is easy to study the eigenvalues of AP −1 instead of P −1 A. From the form of AP −1 , we notice that the preconditioned matrix has only two distinct eigenvalues ±1 and then we notice that a minimal polynomial of degree at most 2. Hence, when Krylov subspace methods like FGMRES is used, then it converges in 2 iterations or less, in exact arithmetic.This property is of practical use when inexpensive approximations of the Schur complement exist.However, when we approximate the Schur complement matrix S by the matrix S 1 or S 2 , we have the following eigenvalue estimation.
Theorem 3. If the linear system (1) is left preconditioned by the matrix P 1 or P 2 , then the eigenvalues of the preconditioned matrices (39) are described as follows: where σ 1 and σ n are the minimum and the maximum eigenvalues of the matrix (−SS 1 −1 ) or (−SS 2 −1 ).
Example 2. In this example, our aim is to verify that the bounds given in the above theorem are matched.We take N = 16, i.e., n = ( 16) 2 = 256 and we fix α = 1.4,β = 0.1, λ = 0.001.For this task, we use the preconditioner P 1 and we use the test image "Golden House".We notice that the positive eigenvalues are equal to one whereas the negatives are contained in the interval [σ 1 , σ n ] where σ 1 and σ n are defined in the above theorem.In this example σ 1 ∼ = −1.01 and σ n ∼ = −1.The results of this example are plotted and shown in Figures 7 and 8.Moreover, in this experiment, we find that the cond(A) = 3.2915 × 10 4 and cond((P 1 ) −1 A) = 1.6219 which indicate that our preconditioner is effective.From Figures 7 and 8, we notice that the preconditioned matrix has a good clustering behavior of the eigenvalues.The eigenvalues are clustering around 1 and −1.This clustering verifies the above theorem guarantees fast convergence of the FGMRES method.

Numerical Results
In this section, we experimentally study the performance of the FGMRES method with the proposed preconditioners P 1 and P 2 .In the following numerical experiments, we implement Algorithms 1-3, and we take the zero vector to be the initial guess.We stopped the outer iterations (FGMRES) when the residuals satisfies b − Ax k < 10 −7 b where x k = (v k , u k ) is the solution vector in the k − th iteration.We used only one iteration of the Fixed-Point Iteration method to linearize the nonlinear term and then we used the PCG for the inner iterations and it is stopped when the tolerance is 10 −9 .No restarting is used for FGMRES algorithm.For this purpose, two famous 128 × 128 test images, called Retinal Image and Golden House are used in the experiments, as shown in Figures 9 and 10   In order to show the performance of the proposed preconditioners, we need to calculate the PSNR which is commonly used in the signal processing field.It can be calculated by the following formula: where u e and u d are the exact and deblurred images, respectively.Bigger PSNR means better deblurring performance.

The Parameters β and λ Selecting
The value of the parameters β and λ also play a vital role in the performance of the numerical technique used for the image deblurring model.Small values of β affect the convergence rate of the iterations in the numerical technique but do not change the quality of the deblurred images.We have chosen β = 1, β = 0.1 and β = 0.01 which are commonly used in the literature [28].We noticed no significant difference in the results between these values.Regarding the values of the regularization parameters λ, we have chosen λ small enough, 10 −3 , 10 −5 10 −6 and 10 −8 to ensure the best deblurring performance of the corresponding deblurring model.These values are commonly used in the literature [28].
Example 3. In this example, we show the impact of our preconditioners on the convergence speed of the FGMRES algorithm for the fractional-order image deblurring problem.We fix N = 128, β = 0.1, and the regularization parameter λ = 0.001.No restarting is used for the FGMRES algorithm and it is stopped when the tolerance is 10 −7 .We use the test image "Golden House".In each FGMRES iteration, the logarithm of   In Figures 12 and 13, we show the algorithm of the ratio of the current residual norm to the initial residual norm, plotted against the number of FGMRES iteration, for different values of the regularization parameter λ.NP stands for FGMRES without preconditioners, P 1 stands for FGMRES with the preconditioner P 1 and P 2 stands for FGMRES with the preconditioner P 2 .The results in Figures 12 and 13 show that our block triangular preconditioners P 1 and P 2 significantly accelerate the convergence of FGMRES, compared to FGMRES without preconditioners.Additionally, P 1 outperforms P 2 .
Example 4. In this example, we show the effectiveness of our proposed preconditioners in deblurring images.We used two blurred images (of size 128 × 128) shown in Figures 14 and 15.We select the following parameters: α = 1.8, β = 1, and λ = 0.00001.We used our preconditioners P 1 and P 2 to deblur the images and the results are shown in Figures 16-23.From Figures 16-23, the results show that our preconditioners are effective in deblurring images, with significant improvement in the PSNR.For example, the PSNR of deblurred image in Figure 16 is 49.41, compared to the PSNR 22.978 for the blurred image in Figure 23.Example 5.In this example, we compare the total CPU-time (in seconds) required for the convergence of the FGMRES with and without our proposed preconditioners P 1 and P 2 .The results are shown in Table 1 for different N, α, β, λ.From Table 1, the results show that both P 1 and P 2 can significantly reduce the CPU-time required for convergence, compared to FGMRES without preconditioners.For example, for N = 128, α = 1.6, β = 0.01, and λ = 10 −6 , the CPU-time for FGMRES without preconditioning is 76.64 s, while the CPU-time for FGMRES with P 1 is 35.86 s and the CPU-time for FGMRES with P 2 is 38.22.Overall, the results show our proposed preconditioners P 1 and P 2 are effective in accelerating the convergence of FGMRES for the fractional-order image deblurring problem.This can lead to significant reductions in CPU-time, which is important for practical applications.

GMRES versus FGMRES
In this experimental result, we compared the performance of GMRES and FGMRES with our preconditioner P 1 using the following parameters: N = 64, α = 1.4,β = 0.1, and λ = 10 −5 .We used the test image "Golden House".We used both GMRES and FGMRES.In this example, both GMRES and FGMRES were stopped when the tolerance was 10 −7 and no restarting is used.The comparison results are shown in Figure 24, where P 1 GM stand for GMRES with P 1 and P 1 FG stands for FGMRES with P 1 .As shown in the figure, FGMRES is performed slightly better than GMRES.
In the following numerical result, we show the comparison of our TFOV-based algorithm with TV-based algorithm [28].For numerical calculations, we used the motion kernel.For the TV-based method we used β = 1 and λ varying from 10 −2 to 10 −4 , according to [28].The parameters for TFOV-based method are listed in From Table 2, we observed that the PSNR by the TFOV-based PGMRES method is almost the same as that of the ordinary TFOV-based GMRES method, but much higher than that of the TV-based P 1 method for all values of N.However, the P 1 and P 2 methods generate this better PSNR in much fewer iterations.For example, to achieve a better PSNR the P 1 method needs only 18 iterations, and the P 2 method needs only 20 iterations for N = 64.However, the NP method needs 120+ iterations to get the same PSNR.The TV-based P 1 method also takes 120+ iterations to get its lower PSNR.The same is the case for other values of N.This means that the TFOV-based FGMRES method is faster than the TFOV-based GMRES and TV-based P 1 methods.Example 7. In this example, we utilized satellite images used by Chowdhury et al. [71].The images underwent deliberate blurring and were corrupted by Poisson noise, resulting in the presence of blurring artifacts.To achieve the blurring, we applied a kernel with specific parameters, namely, we used the Gaussian build in kernel " f special( gaussian , 9, sqrt(3)) .The introduction of Poisson noise to the image presents a substantial challenge for most deblurring techniques, as this type of noise frequently occurs in scenarios involving photon counting across various imaging methods.Simultaneously, blurring is an inevitable consequence due to the underlying physical principles of the imaging system, which can be thought of as the convolution of the image with a point spread function.For the sake of comparison, we chose to employ the non-blind fractional order TV-based algorithm (NFOV) proposed by Chaudhury et al. [71].The restored satellite images can be seen in Figures 34-38, with each image sized at 128 × 128.We configured the parameters for the NFOV method as specified in the reference by Chowdhury et al. [71].For comparison, we have used two different values of N.These are 64 and 128.Their corresponding blurred PSNR are 20.2985 and 20.4559 respectively.The computational technique's stopping criterion is determined by a tolerance value of tol = 10 −7 .Additional details regarding this experiment can be located in Table 3.
Remark 2. Upon examining  and Table 3, it becomes apparent that the results generated by all methods are virtually indistinguishable.Nevertheless, our proposed methods (GMRES and FGMRES) exhibit slightly higher PSNR values while demanding significantly less CPU time.This observation underscores the improved efficiency and speed of our suggested methods (GMRES and FGMRES) in comparison to the NFOV technique.

Conclusions
In this paper we have proposed two block triangular preconditioners for solving the generalized saddle point system which is derived from discretizing the Euler Lagrange equations associated with the TFOV in image de-blurring based problems.We have investigated the performance of the proposed preconditioners with the FGMRES method.We have tested this method on three types of digital images.We have also compared our algorithm with TV based algorithm.Our experiments show that the block triangular preconditioners are very effective.We have also shown that our technique improves the quality of the reconstruction images via calculation of the PSNR.We showed the performance of both GMRES and FGMRES with our proposed preconditioner and we concluded that FGMRES is slightly better than GMRES.Few iterations and CPU-time are needed to obtain a fast rate of convergence and good de-blurring performance.Circulant approximations are used in the first term of the Shur complement to reduce the cost of the computation from O(N 2 ) into O(N log N) and reduce the storage.The spectrums of the preconditioned matrices are clustered around 1 and −1.
and they are blurred by the motion kernel as shown in Figure11.

Figure 24 .Example 6 .
Figure 24.FGMRES vs. GMRES.Example 6.In this example, we compare our TFOV-based algorithm with the TV-based algorithm on the nontextured peppers image.We use a Gaussian kernel with standard deviation σ = 1.5.The results are shown in Figures 25-30.The size of each subfigure is 256 × 256.The subfigures are as follows: (a) exact image (b) blurry image (c) deblurred image by TV (d) deblurred image by NP (e)deblurred image by P 1 and (f) deblurred image P 2 .For numerical calculations, we used the motion kernel.For the TV-based method we used β = 1 and λ varying from 10 −2 to 10 −4 , according to[28].The parameters for TFOV-based method are listed in Table2.For comparison we used three different values of N: 64, 128 and 256.Their corresponding blurred PSNRs are 20.1827,20.1124 and 20.5531 respectively.For the stopping criteria of the numerical methods, we used tolerance tol = 10 −7 .

Table 1 .
The CPU time comparison of GMRES and FGMRES.

Table 2 .
For comparison we used three different values of N: 64, 128 and 256.Their corresponding blurred PSNRs are 20.1827,20.1124 and 20.5531 respectively.For the stopping criteria of the numerical methods, we used tolerance tol = 10 −7 ., we can clearly see the effectiveness of preconditioning.For all values of N, the number of P 1 and P 2 iterations is much lower than the number of TFOV-based NP and TV-based P 1 iterations to reach the required accuracy tol = 10 −7 .The later fixed-point iterations also have similar results.3.

Table 2 .
The PSNR Comparison of TV, NP, P 1 and P 2 .