A Simple but Universal Fully Linearized ADMM Algorithm for Optimization Based Image Reconstruction

Background and Objective: Optimization based image reconstruction algorithm is an advanced algorithm in medical imaging. However, the corresponding solving algorithm is challenging because the optimization model is usually large-scale and non-smooth. This work aims to devise a simple but universal solver for optimization models. Methods: The alternating direction method of multipliers (ADMM) algorithm is a simple and effective solver of the optimization models. However, there always exists a sub-problem that has not closed-form solution. One may use gradient descent algorithm to solve this sub-problem, but the step-size selection via line search is time-consuming. Or, one may use fast Fourier transform (FFT) to get a closed-form solution if the system matrix and the sparse transform matrix are both of special structure. In this work, we propose a simple but universal fully linearized ADMM (FL-ADMM) algorithm that avoids line search to determine step-size and applies to system matrix and sparse transform of any structures. Results: We derive the FL-ADMM algorithm instances for three total variation (TV) models in 2D computed tomography (CT). Further, we validate and evaluate one FL-ADMM algorithm and explore how the two important factors impact convergence rate. Also, we compare this algorithm with the Chambolle-Pock algorithm via real CT phantom reconstructions. These studies show that the FL-ADMM algorithm may accurately solve optimization models in image reconstruction. Conclusion: The FL-ADMM algorithm is a simple, effective, convergent and universal solver of optimization models in image reconstruction. Compared to the existing ADMM algorithms, the new algorithm does not need time-consuming step-size line-search or special demand to system matrix and sparse transform. It is a rapid prototyping tool for optimization based image reconstruction.


Introduction
Image reconstruction is the core technique of medical imaging. There are three types of reconstruction algorithms: analytical algorithm [1], optimization based, iterative algorithm [2,3] and deep learning/ machine learning based algorithm [4]. Among them, the optimization based algorithm has some advantages relative to the two other algorithms. Compared to the analytical algorithm, the optimization based algorithm may incorporating prior information into the optimization model and thus may achieve more accurate reconstruction from incomplete data and/or noisy data. Compared to the deep learning based algorithm, the optimization based algorithm is more stable and is not data dependent [5].
Therefore, the optimization based algorithm is always a hot spot in image reconstruction eld.
The imaging system model of the optimization based algorithm is a linear system of equations. Usually, it is large-scale, ill-posed and often underdetermined. So, it is impossible to solve the linear system by direct inversion. We may continue to construct an optimization model to solve the ill-posed and underdetermined problems by incorporating prior information like sparse prior [6] or low-rank prior [7], etc.
In the optimization models for image reconstruction, total variation (TV) type models have been widely used for they may not only suppress streak artifacts but also suppress noise [8]. In 2006, Sidky et al. proposed a data divergence constrained, TV (DDcTV) minimization model for 2D CT [9]. In 2008, Sidky et al. proposed a DDcTV model for 3D CT and achieved accurate reconstruction [2]. From then on, a variety of TV-type models were proposed in image reconstruction, for example adaptively weighted TV (awTV) [10], anisotropic TV (aTV) [11], high order TV (HOTV) [12,13], non-local TV (NLTV) [14], nuclear TV (nTV) [15], directional TV (dTV) [16], etc. These TV-type models may achieve more accurate reconstructions in their suitable cases. However, the solving problems of these TV-type models are challenging because these models are large-scale and non-smooth. They cannot be solved by simple gradient decent algorithm.
For TV-type models, there are mainly three solving algorithms: adaptive steepest descent-projection onto convex sets (ASD-POCS) algorithm [2,9], Chambolle-Pock (CP) algorithm [17,18], and alternating direction method of multipliers (ADMM) algorithm [19][20][21][22][23][24][25][26]. We have systematically studied ASD-POCS and CP algorithms and applied them to solve TV models for CT and electron paramagnetic resonance imaging (EPRI) [27][28][29][30]. ASD-POCS algorithm is devised according to the physical meaning and the optimization strategy, whereas it is not derived mathematically. Thus, those algorithm parameters in this algorithm should be carefully, empirically chosen to achieve convergence. CP algorithm is obtained by mathematical derivation and may guarantee convergence, but the derivation is not so easy for those who are not familiar with optimization, for example it needs the calculation of convex conjugate [18]. ADMM is another algorithm framework for solving TV-types models. The derivation of algorithm instance according to the ADMM algorithm framework is easier than CP algorithm for the core step is only to split the whole optimization problem into several simple sub-problems by alternating direction technique.
However, for TV-types models in image reconstruction, there is always a sub-problem that has not closedform solution in ADMM algorithm. Thus, this sub-problem is not easy to solve and people began to search techniques to process this issue.
Let's focus on TV minimization in image denoising, restoration and reconstructions and see how the ADMM-like or ADMM algorithm is used to solve the TV models and how the di cult sub-problem is solved. These models have similar formulation to (20) of this paper. To better describe this algorithm development line, we use the symbols in (20) uniformly. Now, let us regard (20) as a general image processing model. If A is an identity matrix, then it is a model for image denoising.
If its function is to blur an image, then it is an image restoration model. Here, A has a special structure for its operation is equivalent to convolution to an image via a blurring mask. In image reconstruction, however, A is just a normal matrix without special structure. In (21), D is the gradient transform matrix, which has special structure for its operation on an image is equivalent to convolution to an image via a nite-difference mask. For the aim of this work is to devise a simple but universal ADMM-type algorithm, we should note the simplicity and universality of each algorithm below.
In 2007 and 2008, Wang et al proposed an alternating minimization algorithm for TV model in image restoration/deblurring [20,31]. In fact, this algorithm is an ADMM-like algorithm. One may regard it as a simpli ed ADMM algorithm. This algorithm introduces an auxiliary variable to replace D u , uses the penalty technique to ensure the equivalence of the substitution, and then use the variable-splitting technique to solve the TV model. After splitting, there are two sub-problems: one is data-delity subproblem and the other one is TV-regularization sub-problem, which may be solved by the commonly used shrinkage operation. However, the data-delity sub-problem has not explicit closed-form solution.
Fortunately, in this sub-problem, A and D are both of special structure for their operations on an image are equivalent to convolution to an image. Thus, this sub-problem may be solved by 2D fast Fourier transforms (FFT). However, we must note that this algorithm cannot be used in image reconstruction for the system matrix, i.e.
A , has not this special structure.
In 2008, Huang et al proposed an alternating minimization algorithm for TV model in image restoration [21]. They introduce a new variable that is equal to the unknown image, then use the penalty technique to formulate the original objective function into the new one. Next, they use alternating minimization to split the problem into two sub-problems. The rst one is the data-delity sub-problem which may be solved by the FFTs. The second one is the TV-regularization sub-problem which may be solved by the Chambolle's projection algorithm. However, we must also note that this solving algorithm still cannot be used in image reconstruction for the system matrix has not special structure, i.e, its operation is not equivalent to a convolution operation.
In 2009, Goldstein et al proposed the split Bregman algorithm for solving TV model in magnetic resonance imaging (MRI) [22]. Similar to the above solving algorithm, the split Bregman algorithm includes two important sub-problems. One may be solved by FFTs, whereas the other one is solved by shrinkage operation. The advantage of split Bregman algorithm over alternating minimization algorithm is that it has not the penalty parameter whose ideal value is in nite. But, this algorithm still cannot be used in CT and EPRI reconstruction for it utilizes the special structure of the system matrix in MRI reconstruction and use FFTs to solve the data-delity sub-problem. Note that the split Bregman algorithm is equivalent to the ADMM algorithm.
In 2010, Yang et al proposed the ADMM algorithm for TV models in MRI reconstruction [25]. It is very similar to the Split Bregman algorithm. Each iteration only involves simple shrinkages and FFTs. Still note that it cannot be used in image reconstruction whose system matrix has not special structure.
Clearly, these ADMM-like or ADMM algorithms all have demand to the system matrix so as to use FFTs, so they are not universal solver. To approach a universal solver, people need to explore new techniques.
In 2010, 2011 and 2013, Li et al proposed their universal ADMM algorithm for TV model in image reconstruction [26,32,33]. They use the nonmonotone line search to decide the step-size of the gradient descent algorithm for solving the data-delity sub-problem. Since this is not a closed-form solution, this inner-iteration should be done many times, for which the time-consuming line search is necessary.
In 2012, Xiao et al proposed the linearized ADMM (L-ADMM) algorithm for TV model in compressed sensing problems [24]. They linearize the data-delity term for solving the corresponding sub-problem and then the FFTs may be used to get a closed-form solution. This algorithm allow the system matrix has a general structure, so it is an universal solver for any type of system matrix. But this algorithm still utilize the special structure of D . If D is another sparse transform that is not equivalent to convolution operation, then the FFTs cannot be used to achieve closed-form solution.
Clearly, linearized ADMM algorithm has potential to approach universal solver for optimization models in image reconstruction.
In 2012, Chan et al proposed an L-ADMM algorithm for constrained linear least-squares problem in image deblurring [34]. They linearized the quadratic regularization term and got a simple closed-form solution for this regularization sub-problem. However, the data-delity sub-problem still uses the special structure of the system matrix whose function is blurring so that the closed-form solution may be achieved by FFTs. This algorithm is not a universal solver for image reconstruction. But we may nd the potential of L-ADMM to achieve closed-form solution.
In 2013, Yang et al proposed an L-ADMM algorithm for nuclear norm minimization [35]. They linearized the data-delity term so as to get a closed-form solution.
In In 2019, Liu et al applied the L-ADMM to a non-convex non-smooth optimization problem and achieved good performance [39].
Clearly, the L-ADMM algorithm may simplify the solving problem of the di cult sub-problems in ADMM algorithm for it may construct closed-form solution. However, we found that these L-ADMM algorithms are not thorough, i.e. people only linearized one quadratic term in the di cult sub-problem. Why do not people linearize all the quadratic terms in the di cult sub-problem? We think it should be deeply investigated. By fully linearization, we expect that a simple but universal solver of closed form may be constructed for the di cult sub-problem in ADMM algorithm instance. 'universal' means that the algorithm does not demand the special structure of the system matrix and the sparse-transform matrix. 'simple' means that the core operations only involve simple matrix-vector multiplications and some simple closed-form operations like shrinkage or projection.
In this work, we proposed a fully linearized ADMM (FL-ADMM) for image reconstruction to simplify the ADMM algorithm by avoiding the search of optimal step-size in gradient descent algorithm and the use of FFT algorithm. The proposed FL-ADMM algorithm framework may be used to derive simple, effective, universal, and convergent algorithm instances for a variety of optimization models in image reconstruction.
To show the potential of the FL-ADMM for prototyping of the optimization models, we derive the FL-ADMM algorithm instances for unconstrained TV (uTV) minimization model, data divergence constrained, TV (DDcTV) minimization model, and TV constrained, data divergence minimization (TVcDM) model for two dimensional (2D) computed tomography (CT).
Also, we validate and evaluate the DDcTV-FL-ADMM algorithm for 2D CT to illuminate that the FL-ADMM algorithm is actually an accurate solver of the DDcTV model. In addition, we explore how the penalty parameter and other algorithm parameters of DDcTV-FL-ADMM impact the convergence rate. Finally, we compare the algorithm with another established universal solver, CP algorithm, to demonstrate its performance.
In Section 2, we give the derivation of the FL-ADMM algorithm instances for three TV models. In Section 3, we perform reconstruction experiments via the proposed DDcTV-FL-ADMM algorithm. We give deep discussions and draw brief conclusions in Sections 4 and 5, respectively.

Preliminary knowledge
In this section, we give some basic optimization knowledge, which will be used in the following parts.

Shrinkage algorithm
(1) One dimensional (1D) shrinkage 1D Shrinkage algorithm or operation may solve the optimization problem shown in Eq. (1). Here, x is a vector of size N, Suppose that Here, S is the 1D shrinkage operator and is the standard sign function, whose value is 1 for positive number, 0 for 0, and − 1 for negative number. is an operator for selecting maximal value.
(2) Two dimensional (2D) shrinkage  Here, x is a 2D-vector-valued vector of size N, The solution of this optimization problem is Here, S 2 is the 2D shrinkage operator and we should note that

Projection algorithm
In optimization theory, projection onto convex sets (POCS) algorithm/operation may solve a special optimization problem as follows.
Here, C is a convex set.
Then, the solution of this optimization problem is Here, P is the POCS operator. The meaning of (6) is to project the point Here, is the There exists accurate algorithm for P , shown in Algorithm 1 [40].

Algorithm 1
Pseudo-codes for 10: returnw The indicator function may be de ned as Thus, a constrained optimization model may be transformed into the unconstrained optimization form. For example, Eq. (5) may be written in the unconstrained form as follows.
In the following sections, we will often use this transformation. Note that, it is just a form-transformation and the optimization meaning of the two forms are completely the same.

Linearization of a quadratic function
In fact, any quadratic function may be linearized at a point according to the Taylor expansion. Here, Here, λ m a x means the largest eigenvalue of a matrix.

ADMM algorithm framework
We consider the structured constrained convex optimization problem Where, are convex but not necessary smooth functions.
The corresponding augmented Lagrange function is The ADMM algorithm framework is 3 ) The sub-problem (14.1) may be further written as The sub-problem (14.2) may be further written as Clearly, the ADMM algorithm may divided the original problems into three simple sub-problems. By this splitting technique, the two convex functions are split and the whole solving problem is potentially simpli ed.
Often, (15) is still hard to solve for usually there is not a simple, close-form solution, for example in case that is also a quadratic function. Thus, people proposed the LADMM algorithm, whose difference from the ADMM algorithm is only linearization to the quadratic function in (15).
According to the linearization method shown in Section 2.1.4, the linearized form of (15) is (17) is potentially be easier to get a close-form solution.

Fully Linearized ADMM (FL-ADMM) algorithm framework
is a simple quadratic function, i.e. the corresponding matrix is an identity matrix or is of special structure, for example it is equivalent to a convolution operation, then (17) may get a close-form solution However, the matrix corresponding to usually is not of special structure in CT and EPRI, so (17) is still di cult to get a close-form solution.
The difference of FL-ADMM from L-ADMM is only that it linearizes any quadratic functions. By use of this fully linearization technique, (17) may get a close-form solution. If , according to (17) and (10), then (17) has the close-form solution below. Where, , and Clearly, FL-ADMM algorithm may make close-form solution and simplify the di cult sub-problem. Otherwise, one must use time-consuming line search to select the step-size for solving (17) via gradient descent or one must use FFT algorithm to solve (17), for which, F and A must be of special structure.
This proposed FL-ADMM algorithm is the core contribution of this work.

FL-ADMM algorithm instances derivation of three TV models
In this section, we will derive three FL-ADMM algorithm instances corresponding to three TV models.
TV models in image reconstruction have been heavily investigated and have achieved accurate reconstructions via sparse-view projections and/or noisy projections. There are three types of TV models: ucTV, DDcTV and TVcDM models. The unconstrained version is often solved by ADMM algorithm, whereas, DDcTV and TVcDM are often solved by ASD-POCS and CP algorithm. Though the derivation below, we will see that FL-ADMM algorithm may solve any type of TV model and each sub-problem is simple to compute.
Without loss of generality, we derive the algorithm instances for 2D CT.

FL-ADMM algorithm for ucTV model (1) The ucTV model
The discrete-to-discrete (D2D) imaging system model of 2D CT may be formulated as where,  The ucTV model may be formulated as where, is the TV norm of image u , and its isotropic form is Here, D is a matrix of size respectively, shown as Here, Thus, the gradient magnitude transform And, the TV norm is (2) The ucTV-FL-ADMM algorithm Now, we have the ucTV model, The FL-ADMM algorithm instance derivation process is as follows.
, then (27) is equivalent to this minimization problem, The corresponding augmented Lagrange function is According to the FL-ADMM algorithm framework, we may derive the algorithm instance.
We perform linearization to each quadratic function in (30), and get Let the gradient of the objective function in (31) be 0, we get Here, According to the 2D shrinkage algorithm in (4), we get , respectively. But, we should note that they are both a 2D vector, so (34) is a 2D shrinkage algorithm.
Now, we get the ucTV-FL-ADMM algorithm. technique. Also, it does not need to use FFT algorithm. The main operations in this algorithm are just simple matrix-vector multiplication and shrinkage algorithm. When we implement this algorithm, the matrices do not need to be explicit form, i.e. we must not store them in the computer memory. We may regard them as some speci c operations, for example, means projection operation, whereas means backprojection operation.
However, we nd that this algorithm converges very slowly, so we propose an improved algorithm with inner iterations as follows.
This model has superior performance than the ucTV model since its model parameter ϵ has clear physical meaning that it embodies the noise level and level of system-inconsistence.
(2) The DDcTV-FL-ADMM algorithm According to (21) and (8), (36) may be written as Let y = D u , and , then (37) is equivalent to this minimization problem, The corresponding augmented Lagrange function is According to the FL-ADMM algorithm framework, we may derive the algorithm instance.
By linearizing the two quadratic functions in (40), letting the gradient of the objective function be 0, we may get Here,  According to the 2D shrinkage algorithm in (4), we get  The TV constrained, data divergence minimization (TVcDM) algorithm is another constrained TV model and has been widely used in CT and EPRI. They are often solved by the CP algorithm. Next, we will see that it may also be solved by FL-ADMM and will see that its derivation is easier than that of CP algorithm for it does not need to calculate the convex conjugate functions.
The TVcDM model may be formulated as (2) The TVcDM-FL-ADMM algorithm By use of the indicator function and the de nition of the TV norm, (47) is equivalent to Then, the corresponding augmented Lagrange function is According to the FL-ADMM algorithm framework, we may derive the corresponding algorithm instance.
It is the same as (30). By full linearization, we may get its close solution.
Here,   Now, we get the TVcDM-FL-ADMM algorithm.  The aim of this work is to design a FL-ADMM algorithm framework and to give prototyping method for optimization models in image reconstruction via the FL-ADMM algorithm framework.

Results
For optimization based image reconstruction, we know that optimization model decides the nal solution, whereas the solving algorithm just impact the convergence rate and path. Since the FL-ADMM algorithm is a solving algorithm, we should validate if it may solve optimization model accurately, and evaluate what may impact its convergence rate.
We have derived three FL-ADMM algorithms for ucTV, DDcTV and TVcDM models. Next, we will validate and evaluate this solving algorithm for DDcTV model in 2D CT, no loss of generality. We

Inverse crime of DDcTV-FL-ADMM algorithm
Inverse crime is a tool to validate the correctness of an inverse problem [41]. If the projection data is complete and exact, then the sign of inverse crime is that the reconstructed image is almost the same with the truth image, i.e. there is not any error except for the computer oating point error.
The imaging con gurations are as follows. The phantom is the Shepp-Logan phantom of size [256,256]. The imaging coordinate system is located at [128,128]. The projection signal at each angle is of size 256 and its coordinate range is [-127, 128]. The virtual detector bin length is 1. The pixel size is also 1. The parallel beam scanning pattern is adopted. The projections are simulated by use of pixel-driven method [42]. We collect 256 projections uniformly distributed in the range of , and the inner iteration number is 50. The gray image has only 256 gray-level. So, if RMSE is less than , the reconstructed image and the truth image will be visually the same, which may be the sign of inverse crime. More strictly, we de ne the sign of inverse crime as At iteration 4570, inverse crime is achieved. The reconstructed images and the corresponding pro les are shown in Fig. 3. Three iteration trends are plotted in Fig. 4.
From Fig. 3, we may see that the reconstructed image is visually the same with the truth image, and that the vertical-center-line pro les of the reconstructed image and the truth image are completely coincident.
This shows that inverse crime is achieved. Now, we may think that the imaging system model, the DDcTV model, the FL-ADMM algorithm and its computer implementation are all correct.
In Fig. 4, we plots three iteration trends to observe the iteration behavior. It may be seen that the RMSE of the reconstructed image and the data error between the guessed data and the truth data may go down and down with the increase of iteration number. At iteration 4570, the RMSE has been less than which means the inverse crime is achieved. We may see that these curves still has a descent trend, showing its good convergence performance. From Fig. 4 (c), we may see that during the later iteration period, the TV value may go down and down with the iteration marching. When the iteration stops, the TV value of the reconstructed image becomes the same with the truth TV. From the image at row 2 and column 5, we may see obvious artifacts. This means that the FL-ADMM algorithm of β value of 100 suffer from the slowest convergence rate. This observation may be more clearly seen in Fig. 6. The order of convergence rate of different β from slow to fast is 100, 0.01, 0.1, 1, and then 10. This indicates that an appropriate β value may achieve fast convergence, whereas too large or too small values always lead to too slow convergence. Even that too huge value may lead to a wrong convergence, i.e. the RMSE will converge to a very large value. However, we want to emphasize that the optimal β value is imaging-condition dependent. In this simulation study, the optimal value is 10. But, in other cases, it may be other value. Usually, one may vary the β value with interval of an order of magnitudes, then select the optimal one to achieve the fastest convergence.

The impact of inner iteration number on convergence rate of the DDcTV-FL-ADMM algorithm
In Algorithm 3-B, i.e. the DDcTV-FL-ADMM algorithm, there is an inner loop process. In this section, we investigate how the inner iteration number impact the convergence rate.
The phantom is still the FORBILD phantom of size [256,256]. The imaging condition is the same with that of Section 3.2. In the DDcTV-FL-ADMM reconstructions, we vary the inner iteration number from 1, 50, 100, 150, to 200 and x the iteration number as 5000, x β as 10, to evaluate their respective convergence rate. Figure 7 shows the reconstructed images and Fig. 8 plots the convergence curves From Fig. 7, we may see that the reconstructed images with inner iteration number of 50, 100,150 and 200 all have higher accuracy. However, if the inner iteration number is 1, i.e. if the FL-ADMM algorithm used is Algorithm 3-A, the reconstructed image suffers from serious artifacts. This means that appropriately selected inner iteration number may achieve faster convergence. The standard FL-ADMM algorithm without inner iteration always leads to too slow convergence rate. This observation may be clearly seen in Fig. 8. Though use of larger inner iteration number may achieve faster convergence rate, the inner loop process will take longer time. Thus, one should select an appropriate inner iteration number to achieve fast convergence and fast inner loop computation. In this case, inner iteration number of 50 is the optimal selection. Also, we want to emphasize that the optimal number of inner iteration is imagingcondition dependent. Usually, one may select some values to perform reconstructions and then select the optimal one.

Comparison with the CP algorithm
CP algorithm has been proposed and widely used in image reconstruction for more than 10 years. We have applied the CP algorithm in EPR imaging and CT. The most important advantage of the CP algorithm is that it may always get closed-form solutions for sub-problems and the nial algorithm instance only involves simple matrix-vector multiplications and some simple operations. Thus, we say the CP algorithm is a fast prototyping tool for optimization based image reconstruction. Similar to CP, the FL-ADMM is also a fast prototyping tool for it may always get closed-form solutions for di cult subproblems and it has not special demand on system matrix and sparse-transform matrix.
Through the studies in subsections 3.1 to 3.3, we have known the correctness of the FL-ADMM algorithm and have realized how the β selection and inner iteration number impact the convergence rate. In this subsection, we investigate the sparse reconstruction capability of the DDcTV-FL-ADMM algorithm by comparison with the DDcTV-CP algorithm whose pseudocode and some explanations are shown in Appendix 1.
The phantom is a real thoracic CT image of size [256,256]. The imaging condition is the same with that of Section 3.2. In the DDcTV-FL-ADMM reconstructions, we x the inner iteration number as 50, x    and PSNR, respectively. RMSE means root mean square error, SSIM means structural similarity index measure, and PSNR means peak signal to noise ratio. From them, we may see that the two algorithms have very close accuracy, which is visually validated by Fig. 9. However, we may also see that the FL-ADMM algorithm is always a little bit better than the CP algorithm. This is because the optimal β selection in FL-ADMM is easier than the optimal selections of

Discussions
In this work, we propose a novel ADMM algoithm, FL-ADMM algorithm, which may be used as a prototyping tool for optimization model in image reconstruction. The key operation is to expand all the quadratic terms so that the corresponding sub-problem may get a simple closed-form solution. Further, we propose the fast FL-ADMM algorithm by use of the inner iteration technique. We have derived three FL-ADMM algorithm instances for three TV models, ucTV, DDcTV, and TVcDM. Further, we validate and evaluate the correctness and sparse recontruction capability of the DDcTV-FL-ADMM algorithm. Also, we analyze how the penalty parameter and the inner iteration number impact the convergence rate. In addition, we compare this algorithm with the SOTA CP algorithm and discuss its potential superioty.
In optimization based image reconstruction, especially in TV-type norm based image reconstruciton, the ADMM algorithm always has a problem that one sub-problem has not simple closed-form solution.
Usually, people has two choices to solve this problem. One is to use gradient desent algorithm to solve this sub-problem. However, the step-size selection is di cult. If one uses the line search method to select the optimal step-size for each iteration, it needs too long time, especially when the imaging model is large scale. The other one is to use FFT technique. In fact, why FFT may be used is because the sparse transform matrix and the system matrix may both be regarded as a convolution operation. So, this method has not universality. Once a sparse transform cannot be regarded as a convolution operation, this method loses e cacy.
These di cuties of ADMM motivated this work whose aim is to devise a method to simplify the implementation of the ADMM algoirthm. Motivated by the linearization technique of the L-ADMM, we propose the FL-ADMM algorithm and propose its acceleration version.
In Section 2, the standard FL-ADMM algorithms are named Algorithm-A, whearas the acclerated FL-ADMM algorithm are named Algoritm-B. For example, for the DDcTV-FL-ADMM algorithm, the standard algorithm is Algorithm 3-A, whereas its improved, accelerated version is Algorithm 3-B.
In fact, the accelerated FL-ADMM algorithm uses a special gradient desenct algorithm. In Algorithm 3-B, Eq. (4.3) is, in fact, a gradient desent equation. But, very importantly, the step-size is , which may be calculated before the whole iteration process. Thus, compared with the ordinary gardient descent algorihtm, this proposed method will be much faster for it doesnot need the time-consuming and complicated step-size search via line search technique.
Compared to the ADMM algorithm using FFT, this proposed algoirthm has universality for it does not need that the sparse transform matrix and the system maxtrix may be both regarded as a convolution operation. Just like the CP algorithm, the FL-ADMM algorithm is not so fast. So, in the future, research on how to speed up this algorithm is necessary. However, this is out of the scope of this work. Just like CP algorithm is a fast prototyping tool for optimization model because of its simplicity and universality, the FL-ADMM algorithm is also a fast prototyping tool. Once one designs a new optimization model for image reconstruction, he/she may derive the FL-ADMM algorithm instance quickly and begin to evaluate the performance of this new model. This is the meaning of fast prototyping tool in optimization based image reconstruction.

Conclusions
Page 31/40 The conclusions of this work may be summed up brie y as the following points.
1. The FL-ADMM algorithm is a universal, simple, effective, and accurate solver of convex optimiztion model in image reconstruction, no matter it is unconstrained or constrained.
2. The FL-ADMM algoirthm improves the traditional ADMM algorithm by avoiding the step-size search for gradient desecnt and the special demands on the sparse transfrom and the system matrix for use of FFT technique. 3. The penalty parameter in this proposed algoritm may impact the convergence rate. Too large or small values both lead to slow convergence. 4. The inner iteration number in this algorithm may impact the convergence rate. One may select the optimal one by running several reconstrutions with different inner iteration number. According to our experience, 30-50 is e cient. But, we note that the optimal inner iteration number is imagingcondition dependent.
5. Compared to the CP algorithm, its algorithm-parameters are easy to tune and thus may achieve higher accuracy.
. Compared to the CP algorithm, it may solve some special optimization model which is di cult for CP because of the use of convex conjugate functions.
In the future, the FL-ADMM algorithm should be focused on its acceleration technique which may borrow the ideas in accelertion of the traditional or linearized ADMM algorithms. Competing Interest The authors declare that they have no known competing nancial interests or personal relationships that could have appeared to in uence the work reported in this paper.

Declarations
Ethics approval This is an observational study. The Shanxi University Research Ethics Committee has con rmed that no ethical approval is required.
Consent to participate Consent to participate is not required for the real CT image comes from public data-set of AAPM.   The reconstructed images of the DDcTV-FL-ADMM algorithm with different β. The number above the images indicate the β value. The images on the second row is the zoomed-in images in the red box shown in the rst row. The display window is [0,1].

Figure 6
Plots of RMSE as function of iteration number. The iteration number is 5000. The two axis are both of logarithm form.

Figure 7
Page 39/40 The reconstructed images of the DDcTV-FL-ADMM algorithm with different inner iteration number. The number above the images indicate the inner iteration number. The images on the second row is the zoomed-in images in the red box shown in the rst row. The display window is [0,1].

Figure 8
Plots of RMSE as function of iteration number. The iteration number is 5000. The two axis are both of logarithm form.