On the preconditioning of the primal form of TFOV-based image deblurring model

To address the staircasing problem in deblurred images generated by a simple total variation (TV) based model, one approach is to use the total fractional-order variation (TFOV) image deblurring model. However, the discretization of the Euler–Lagrange equations for the TFOV-based model results in a nonlinear ill-conditioned system, which adversely influences the performance of computational methods like Krylov subspace algorithms (e.g., Generalized Minimal Residual, Conjugate Gradient). To address this challenge, three novel preconditioned matrices are proposed to improve the conditioning of the primal model when using the conjugate gradient method. These matrices are designed based on circulant approximations of the matrix associated with blurring kernel. Experimental evaluations demonstrate the effectiveness of the proposed preconditioned matrices in enhancing the convergence and accuracy of the conjugate gradient method for solving the primal form of the TFOV-based image deblurring model. The results highlight the importance of appropriate preconditioning strategies in achieving robust and high-quality image deblurring using the TFOV approach.

To address the staircasing problem in deblurred images generated by a simple total variation (TV) based model, one approach is to use the total fractional-order variation (TFOV) image deblurring model.However, the discretization of the Euler-Lagrange equations for the TFOV-based model results in a nonlinear ill-conditioned system, which adversely influences the performance of computational methods like Krylov subspace algorithms (e.g., Generalized Minimal Residual, Conjugate Gradient).To address this challenge, three novel preconditioned matrices are proposed to improve the conditioning of the primal model when using the conjugate gradient method.These matrices are designed based on circulant approximations of the matrix associated with blurring kernel.Experimental evaluations demonstrate the effectiveness of the proposed preconditioned matrices in enhancing the convergence and accuracy of the conjugate gradient method for solving the primal form of the TFOV-based image deblurring model.The results highlight the importance of appropriate preconditioning strategies in achieving robust and high-quality image deblurring using the TFOV approach.
Variational methods are an approach used in image deblurring in the last few decades, to restore sharpness and clarity to blurry images.They involve formulating an optimization problem that aims to find the best estimate of the original sharp image given the observed blurry image.In variational methods, following mathematical model is defined to describe the degradation process that caused the image blur.
where u is an original image, z is a recorded image, ε is the noise function, and If the blurring operator K is given, then the corresponding approach is referred to as non-blind deconvolution [1][2][3] .However, when the blurring operator is unknown, the corresponding approach is referred to as blind deconvolution 4-7 .In this paper, our primary focus is on non-blind deconvolution.Here K represents a first kind Fredholm-integral operator.Therefore, it is compact and the problem (1) becomes ill-posed [8][9][10] .Let be a twodimensional square domain.This model typically involves convolution with a blurring kernel that represents the blurring effect.The goal is to find the original image that, when convolved with the blurring kernel, closely matches the observed blurry image.The following optimization problem is formulated by constructing an objective function that consists of two terms: a data fidelity term (first term) and a regularization term (second term).
The data fidelity term measures the mismatch between the observed blurry image and the estimated sharp image after convolution with the blurring kernel.The regularization term J(u) encourages the restoration algorithm to produce visually desirable solutions by imposing certain constraints or promoting specific image properties.
The total variation of an image measures the amount of variation or changes in pixel intensities across the image.In a blurry image, neighboring pixels tend to have similar intensities, resulting in a low total variation.By promoting images with higher total variation, we can encourage the restoration of sharp edges and fine details.The idea behind TV regularization is to find an image that simultaneously fits the observed blurry image data and has a high total variation.This is achieved by solving an optimization problem that involves minimizing a cost function, which combines a fidelity term that measures the difference between the observed and reconstructed images, and a regularization term that quantifies the total variation.Minimizing the total variation encourages the reconstruction algorithm to produce an image with sharp transitions between regions of different intensities.It helps in removing blur and enhancing edges while preserving important image structures.By incorporating total variation regularization into the deblurring process, it is possible to obtain visually pleasing and sharper images, effectively reducing the impact of blur and noise in the original blurry image.While TV regularization is a widely used method in image deblurring, it does have some drawbacks.One of the major drawbacks is staircase effect.TV regularization can introduce a "staircase effect" artifact, where edges appear as a series of steps rather than smooth transitions.This effect occurs because TV regularization promotes piecewise constant regions, resulting in a blocky appearance around edges instead of accurately representing their continuous nature.
The TFOV regularization offers a powerful regularization approach for image deblurring problems, combining the advantages of edge preservation, flexibility, noise robustness, and reducing staircase effects.Its effectiveness has been demonstrated in numerous studies and has contributed to the advancement of image deblurring techniques.However, discretizing the TFOV-based model's Euler-Lagrange (EL) equations results in a large nonlinear ill-conditioned system.To solve such systems efficiently is quite challenging for numerical methods.Even the powerful numerical algorithms like Krylov subspace methods (Generalized Minimal Residual (GMRES), Conjugate Gradient (CG) etc.) get slow convergence.One of the remedy for slow convergence is to use preconditioning.
Preconditioning is a technique that uses 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 such that this matrix is easy to invert and the preconditioned matrix P −1 A has a good clustering behavior of the eigenvalues.Because rapid convergence is often associated with a clustered spectrum of P −1 A .In the Preconditioning technique, we solve the system P −1 Ax = P −1 b instead of Ax = b because the new system will converge rapidly when we use a suitable preconditioner.To apply the preconditioner matrix P within a Krylov subspace technique, we should calculate a matrix times a vector at each iteration.Hence, evaluating this product must be cheap.In the literature, several preconditioners [26][27][28][29][30][31][32][33][34] are developed for the nonlinear systems.In this study, we consider the following non-linear system of equations This system is derived by discretizing the EL equations associated with the TFOV based image deblurring problem.The coefficient matrix of this system is of size 2N 2 by 2N 2 , where N is the number of pixels.This matrix is non symmetric, ill-conditioned, dense and huge.These properties make the development of an effective computational method more challenging.We know that using direct methods for solving (6) requires O(N 3 ) and hence they are not applicable here.For this system, iterative methods like Krylov subspace (GMRES, CG etc.) methods are applicable.However, their convergence is too slow because they are sensitive to the condition numbers.Hence, the idea of preconditioning is needed to accelerate the convergence of the Krylov subspace methods.In this study, we propose three circulant symmetric positive definite (SPD) preconditioners for system (6).The proposed preconditioners not only increase the convergence rate of numerical method but also contribute in the quality of deblurred images.
The manuscript makes the following contributions: (i) It introduces an efficient and fast algorithm for solving the TFOV-based image deblurring problem.(ii) It introduces three novel preconditioned circulant matrices to address the nonlinearity and ill-conditioned characteristics of the large system in the TFOV model.(iii) It offers an improved treatment for the computationally expensive TFOV regularization functional.
The remaining paper is organized into several sections as follows: section "The TFOV-model" describes the TFOV-based image deblurring model.The Euler-Lagrange equations of the TFOV-model are also discussed in this section.Section "Euler-Lagrange equations" presents the discretization of the TFOV-model and its matrix structure.Section "Discretization of the TFOV-model" introduces the numerical implementation.The proposed circulant preconditioners for the PCG (Preconditioned Conjugate Gradient) method are also explained in the same section.In section "Numerical solution algorithm", the numerical experiments are presented.Section "Conclusion" presents conclusions.

The TFOV-model
www.nature.com/scientificreports/with the α-BV norm ||u|| BV α = ||u|| L 1 + � |∇ α u|dx , where α is the order of the fractional derivatives.TV α is defined by where div α φ = ∂ α φ 1 ∂x + ∂ α φ 2 ∂y .∂ α φ 1 ∂x and ∂ α φ 2 ∂y are the fractional-order derivative along the x and y directions, respectively.The space T denotes where |φ(x)| = � 2 i=1 φ 2 i and C ℓ (�, R 2 ) is the space of α-order continuously differentiable functions.Hence, when the total fractional-order variation (TFOV) model is applied, the previous problem is transformed into the equivalent task of identifying a u ∈ BV α (�) ∩ L 2 (�) that minimizes the functional where TV α β is the modified total fractional-order variation and defined by where 2 and β is employed to make TV α β differentiable at zero.D α x and D α y are the fractional derivatives along the subscript directions.The existence and uniqueness of a minimizer for the problem (7) have been extensively investigated, as discussed in works such as 35,36 .

Fractional-order derivatives
Various definitions of fractional-order derivatives have been proposed in the literature to describe such derivatives 37,38 .In this paper, we will present some of these definitions.For a comprehensive mathematical treatment, a fractional-order derivative is represented as a function operator denoted by D α [a,x] .It is important to note that the order α satisfies the condition 0

Euler-Lagrange equations
For the functional (7) and 1 < α < 2 , the Euler-Lagrange equations are as follows: Here, K * is an adjoint operator and L α (u) is given as Proof From ( 7), we have . By applying the Taylor series, we have where Applying α-order integration by parts 35 , we get Then, it suffices to choose ν ∈ C 1 0 (�, R) and we have Hence, (20) reduces to (16). Case-II: Therefore, the boundary terms in (3) can only become zero if Vol

Discretization of the TFOV-model
To discretize � = (0, 1) × (0, 1) , we divide it into N 2 a uniform grid, where N is a positive integer.Let (x k , y l ) for k, l = 0, 1, . . ., N + 1 within 35,36 .Assuming that u has a homogeneous Dirichlet boundary condition and utilizing the shifted Grünwald approximation approach 39,40 , we consider the following derivative: 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 .By using the homogeneous boundary condition, we can obtain the following expression: Now, we have the following properties: (1) Hence, using the Gershgorin circle theorem, we can obtain B α N which is a symmetric and negative definite Toeplitz matrix.We define the solution matrix at (khx, lhy), where k, l = 1, . . ., N .The ordered solution vector of U is represented as . The discrete version of dif- ferentiation for an α-order derivative is given as Similarly, we have where , and ⊗ denotes the Kro- necker product 38,41 .Now, if we use a finite difference scheme and the discrete fractional derivative shown above, then ( 16) and ( 17) lead to the following primal system Let N F be the number of Fixed Point Iterations and K h be a matrix satisfying where [K h U] ij,lm = h 2 k(x i − x j , y l − y m ) .By utilizing the lexicographical order, the matrix K h is structured as a block Toeplitz with Toeplitz block matrix.The discrete numerical method for the matrix L α h (U m ) is given by: Here, • represents the pointwise multiplication operation, and m corresponds to the m-th Fixed Point Iteration.The matrix U is obtained by reshaping the vector u into an N × N matrix.D 1 (U m ) and D 2 (U m ) are diagonal matrices consisting of the element-wise reciprocals of the non-zero matrices B x α (U m ) and B y α (U m ) , respectively.
Vol:.( 1234567890 where Here, more details about the term L TV h (U m ) can be found in 10 .The structure of matrix B h is as follows: where both G 1 and G 2 are of size N(N − 1) × N 2 , and is a matrix of size (N − 1) × N .The matrix H h is a diagonal matrix, which are calculated by discretizing |∇u| 2 + β 2 .It possesses the following structure: where H x and H y are of sizes (N − 1) × N and N × (N − 1) , respectively.

Numerical solution algorithm
We describe the algorithms for solving the TFOV-based linear system (25).Before delving into the details, we present several essential properties of (25).

The Hessian matrix K
) is extremely large in practical applications.When the value of is small, K * h K h + L α h (U m ) becomes highly ill-conditioned.This is mainly due to the clustering of eigenvalues of K around zero 10 .

K *
h K h is symmetric positive definite.Despite K * h K h is dense, K has the translation invariant property, which allows the application of the FFT method to compute K * h K h u in just O(nlogn) operations 10 .3. In the TFOV model (25), the fractional matrix L α h (U m ) is dense, resulting in an expensive matrix-vector multiplication.On the other hand, in the TV model (28), the non-fractional matrix ) is symmetric positive semidefinite 10 .Consequently, the system (25) is symmetric positive definite.

Preconditioned matrices
The conjugate gradient (CG) algorithm is an appropriate iterative algorithm for solving the system (25).The CG algorithm is an iterative algorithm commonly used to solve systems of linear equations.It is particularly wellsuited for large, sparse, and symmetric positive definite matrices.The CG algorithm iteratively finds the solution by minimizing the residual error in each iteration along conjugate directions.It does not require the explicit storage of the entire matrix, making it memory-efficient for large-scale problems.The CG algorithm achieves convergence in a number of iterations equal to the size of the problem, making it an efficient solver.However, the CG algorithm may converge slowly for ill-conditioned matrices, so preconditioned conjugate gradient (PCG) algorithm are often employed to improve convergence speed.In order to achieve an efficient solution, a preconditioning matrix having SPD property is required.In this context, we introduce two SPD circulant preconditioning matrices.The first one is denoted as P 1 , and it is defined as follows: where γ > 0 , Kh is an approximate circulant form of K h , and I h is an identity matrix.The second precondition- ing matrix P 2 is; where diag(L TV h (U m )) is a diagonal matrix having entries from L TV h (U m ) (29).The third preconditioning matrix P 3 is a product type preconditioner.(25), the inversion of the preconditioner matrices ( P 1 , P 2 , and P 3 ) becomes necessary.The inversion of P 1 and P 2 can be easily performed since their second terms are sparse matrices, and the inversion of K * h Kh requires only O(n log n) floating-point operations using the FFTs, as explained in detail in 10 .The inversion of the preconditioning matrix P 3 involves inverting its middle term, γ diag(L TV h (U m )) + I h .Since this middle term is a sparse matrix, its inversion can be done straightforwardly as well.
For the inversion of ( K * h Kh + γ I h ) 1/2 , we also require only O(n log n) floating-point operations using the FFTs.The PCG method is summarized as follows.
Next, let the eigenvalues of K * h K h and L α h (U m ) be µ K i and µ L α i , respectively, such that µ K i ↓ 0 and µ L α i ↑ ∞ .So the eigenvalues of P −1 1 Ā and P −1 2 Ā are respectively.Here µ L TV i are eigenvalues of Clearly, for γ ≡ , Hence, for γ ≡ , P −1 1 Ā and P −1 2 Ā exhibit better spectrum when compared with the Hessian matrix Ā .Now, let eigenvalues of P −1 3 Ā be Therefore, P −1 3 Ā exhibits better spectrum when compared with Ā. Preconditioning techniques can help mitigate issues like noise amplification and ringing artifacts that are often encountered in deblurring processes.By carefully designing the preconditioner, the optimization process can become more stable and effective, leading to better convergence towards a higher quality deblurred image.The preconditioners we introduce serve a dual purpose: they not only address the ill-posed characteristics of (33) Vol:.( 1234567890) www.nature.com/scientificreports/ the problem but also aid in the retrieval of high-frequency intricacies within the deblurred images.This effect becomes evident through the numerical illustrations we provide.

Numerical examples
Three numerical examples for the TFOV-based image deblurring problem are now presented.Various values of N were employed, resulting in a system with N 2 unknowns.MATLAB program, running on an Intel(R) Core(TM) i7-4510U CPU @ 2.60 GHz, was used for numerical computations and to obtain the numerical results.The evaluation of the deblurred images' quality was conducted using PSNR (Peak Signal to Noise Ratio) and SSIM (Structured Similarity Index Measure).The ke − gen(N, r, σ ) kernel [43][44][45] was employed for numerical calculations.

Initial guess and parameters
The works by 20,35,[46][47][48] delve into the intricate details of automatically selecting parameter , a topic that goes beyond the scope of this current paper.The careful choice of value for holds paramount importance in effectively eliminating both blur and noise.It is advised to avoid extreme values, whether they are exceedingly large or exceptionally small.In the context of the TFOV-based model, the range deemed optimal for lies between 1e−3 and 1e−7.As for the optimal range encompassing α and β , we conducted  www.nature.com/scientificreports/computational experiments employing the cameraman image.Our observations indicate that the most suitable range for α is situated between 1.1 and 1.5, while the optimal realm for β spans from 0.6 to 1.The results of these computations are outlined in Tables 1 and 2. It is evident from our experiments that, for lower values of α , higher values of β are conducive.In our specific experimental setup, we made the deliberate selections of α = 1.4 , = 1e -6, and β = 1 to advance our research.

Example 1
The cameraman image was utilized as an example in this study.It is a complex image consisting of both small-scale texture (peacoat) and large-scale cartoon (face) components.Figure 2 illustrates different aspects of the cameraman image, with each subfigure sized at 512 × 512 .These subfigures include: (a) the exact and (b) blurry images; (c), (d), (e), and (f) the deblurred images using TFOV-based CG, P 1 CG with γ = 1e-3, P 2 CG with γ = 1e-4, and P 3 CG with γ = 1e-5, respectively.Figure 3 illustrates the computation of the relative residual at each iteration with respect to γ .Numerical calculations were performed using the ke − gen(N, 300, 5) kernel.The  5. Additionally, it is noted that both P 1 CG and P 3 CG methods achieve slightly higher PSNR and SSIM compared to P 2 CG.Thus, in this example, the performance of preconditioners P 1 and P 3 proves to be more effective than that of preconditioner P 2 .

Example 2
The Coins image, which consists of both real and synthetic components, was used in this example.
A comparison was made between our TFOV-based algorithm and a TV-based method.Because the TV-based method generates a SPD matrix system, the CG (Conjugate Gradient) method was employed for its solution.www.nature.com/scientificreports/ to 1e−3 for P 1 CG, 1e−5 for P 2 CG, and 1e−4 for P 3 CG.The stopping criterion for the numerical methods was set to a tolerance of tol = 1e-3.

Remark
1.The comparison presented in Table 4 indicates that the TFOV-based methods (CG and PCG) achieve slightly higher PSNR and SSIM metrics compared to the TV-based CG method.This observation is further supported by Fig. 4b-f.Hence, the TFOV-based methods (CG and PCG) generate superior image quality results.2. The effectiveness of preconditioning can be clearly observed from Fig. 5.The number of iterations required by TFOV-based PCG is significantly fewer compared to both TFOV-and TV-based CG methods in order to achieve the desired accuracy of tol = 1e-3.4 shows that the TFOV-based PCG method achieves slightly higher PSNR and SSIM values compared to the regular TFOV-based CG method.However, the PCG method achieves its PSNR and SSIM in significantly fewer iterations.The P 1 CG method requires only 50 iterations, the P 2 CG method requires only 14 iterations, and the P 3 CG method requires only 18 iterations to achieve the desired results.In contrast, the TV-based CG method and the TFOV-based CG method both require more than 50 iterations to reach their respective PSNR and SSIM values.This demonstrates that the TFOV-based PCG algorithm is faster than both the TFOV-based CG method and the TV-based CG algorithm for real and synthetic images.Additionally, the effectiveness of preconditioner P 2 surpasses that of preconditioners P 1 and P 3 (Fig. 5).

Table
Example 3 Here, we have utilized a nontexture Moon image.The various aspects of the Moon image are illustrated in Fig. 6.Each subfigure has a size of 512 × 512 .They represent: (a) the exact, (b) blurry; (c), (d), (e) and (f) deblurred images using TFOV-based CG, P 1 CG, P 2 CG, and P 3 CG, respectively.The numerical calculations were performed using the ke − gen(N, 300, 3) kernel.To facilitate comparison, we considered three different values of N: 128, 256, and 512.The stopping criteria for the numerical methods was set to a tolerance of tol = 1e-4.

Remark
1.The similarity between Fig. 6c-f indicates that all methods produce results of the same quality.
2. Figure 7 clearly shows that for all values of N, the number of iterations required by PCG ( P 1 CG, P 2 CG, and P 3 CG) is significantly fewer compared to TFOV-based CG in order to achieve the desired accuracy of tol = 1e -4. 3. Table 5 demonstrates that the PSNR and SSIM values obtained by the PCG method are almost identical to those achieved by the regular TFOV-based CG method for all values of N.However, the PCG method achieves these PSNR and SSIM values in significantly fewer iterations.For instance, for N = 64 , the P 1 CG method requires only 40 iterations, the P 2 CG method requires only 10 iterations, and the P 3 CG method requires only 19 iterations to attain the desired PSNR and SSIM values.In contrast, the CG algorithm requires over 100 iterations to achieve the same results.Similar observations can be made for other values of N. Thus, the PCG algorithm is faster than TFOV-based CG for nontexture images.Furthermore, the performance of all preconditioners is nearly the same (Fig. 7).

Example 4
In this example, a pair of satellite images from 49 were employed.These images were intentionally subjected to blurring and corruption by Poisson noise, leading to the presence of blurring artifacts.For the blurring process, we utilized a kernel with parameters fspecial( ′ gaussian ′ , 9, sqrt(3)) .The introduction of Poisson noise to the images poses a significant challenge for the majority of deblurring techniques.This type of noise commonly arises in scenarios involving photon counting within various imaging modalities.At the same time, blurring is an inevitable outcome due to the underlying physical principles of the imaging system, which can be conceptualized as the convolution of the image with a point spread function.For the purpose of comparison, we opted to utilize the approach by Chaudhury et al. 49 , known as the non-blind fractional order TV-based algorithm (NFOV).The restored images of the Galaxy can be observed in Fig. 8, each possessing dimensions of 256 × 256 .Similarly, the restored images of Satel are depicted in Fig. 9, each sized at 128 × 128 .For the NFOV method, we configured the parameters as outlined in 49 .The stopping criterion for the computational technique is determined by a tolerance value of tol = 1e-7.Further details regarding this experiment can be found in Table 6.
Remark By referring to Figs. 8, 9, and Table 6, it becomes evident that the outcomes produced by all techniques are nearly indistinguishable.However, our suggested PCG methods yield slightly superior PSNR values while requiring significantly less CPU time.This observation highlights the enhanced effectiveness and swiftness of our proposed PCG methods compared to the NFOV method.

Conclusion
We presented a numerical method (PCG) for solving the primal form of the total fractional order variation (TFOV) based nonlinear image deblurring problem.We introduced three novel circulant preconditioned matrices ( P 1 , P 2 , and P 3 ) and tested them with three examples using PCG.Various types of images (real, complicated, non-texture, synthetic and satellite) were also tested using our new circulant preconditioned matrices.Additionally, we compared the TFOV-based algorithm (CG and PCG) with the TV (total variation) based algorithm and  NFOV method 49 .The convergence rates and residual norms at each iteration were provided for each example.
The numerical tests we conducted highlight the swift convergence achieved by the PCG method when employing the novel circulant preconditioners.Beyond the accelerated speed, the PCG method also exhibits efficacy in addressing the TFOV-based nonlinear image deblurring challenge.For this study, the images under consideration were grayscale; however, we anticipate extending our approach to encompass color images in future research.Additionally, we intend to apply our proposed method to images characterized by varying degrees of blurriness.Moreover, the crucial role in the deblurring process is also attributed to the ideal settings of parameters such as α and β .In the present study, we empirically analyze these parameters, while we intend to establish their precise theoretical constraints in subsequent research.

is the greatest integer function. 1 . 2 . 3 .
Definitions of the left-and right-sided Riemann-Liouville (RL) functional derivatives: and where Definitions of the left-and right-sided Grünwald-Letnikov (GL) functional derivatives: and where Definitions of the left-and right-sided Caputo (C) functional derivatives:

Figure 1 .
Figure 1.Schematic illustration of the blurring kernel.

Figure 5 .
Figure 5.The TV-based CG, TFOV-based CG and PCG convergence at fixed point iteration m = 1 for Example 2. Blue asterisk represents TV-based CG iterations, red asterisk represents TFOV-based CG iterations, yellow circle represents P 1 CG iteration, purple box represents P 2 CG iteration and green line represents P 3 CG iteration.

Figure 7 .
Figure 7.The TFOV-based CG and PCG convergence at fixed point iteration m = 1 for Example 3. Blue asterisk represents TFOV-based CG iterations, red circle represents P 1 CG iteration, yellow box represents P 2 CG iteration and black line represents P 3 CG iteration.

Table 3 .
Comparison of TFOV-based CG and PCG for Example 1.

Table 4 .
Comparison of TV-based CG, TFOV-based CG and PCG for Example 2.

Table 5 .
Comparison of TFOV-based CG and PCG for Example 3.