{'fraction_non_alphanumeric': 0.035984414393765755, 'fraction_numerical': 0.008174803269921307, 'mean_word_length': 4.696257615317667, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerated gradient methods for total-variation-based CT image reconstruction

Total-variation (TV)-based Computed Tomography (CT) image reconstruction has shown experimentally to be capable of producing accurate reconstructions from sparse-view data. In particular TV-based reconstruction is very well suited for images with piecewise nearly constant regions. Computationally, however, TV-based reconstruction is much more demanding, especially for 3D imaging, and the reconstruction from clinical data sets is far from being close to real-time. This is undesirable from a clinical perspective, and thus there is an incentive to accelerate the solution of the underlying optimization problem. The TV reconstruction can in principle be found by any optimization method, but in practice the large-scale systems arising in CT image reconstruction preclude the use of memory-demanding methods such as Newton's method. The simple gradient method has much lower memory requirements, but exhibits slow convergence. In the present work we consider the use of two accelerated gradient-based methods, GPBB and UPN, for reducing the number of gradient method iterations needed to achieve a high-accuracy TV solution in CT image reconstruction. The former incorporates several heuristics from the optimization literature such as Barzilai-Borwein (BB) step size selection and nonmonotone line search. The latter uses a cleverly chosen sequence of auxiliary points to achieve a better convergence rate. The methods are memory efficient and equipped with a stopping criterion to ensure that the TV reconstruction has indeed been found. An implementation of the methods (in C with interface to Matlab) is available for download from http://www2.imm.dtu.dk/~pch/TVReg/. We compare the proposed methods with the standard gradient method, applied to a 3D test problem with synthetic few-view data. We find experimentally that for realistic parameters the proposed methods significantly outperform the gradient method.


I. INTRODUCTION
Algorithm development for image reconstruction from incomplete data has experienced renewed interest in the past years. Incomplete data arises for instance in case of a small number of views, and the development of algorithms for incomplete data thus has the potential for a reduction in imaging time and the delivered dosage. Total-varation (TV)-based image reconstruction is a promising direction, as experiments have documented the potential for accurate image reconstruction under conditions such as few-view and limited-angle data, see e.g. [1].
However, it also known that it is difficult to design fast algorithms for obtaining exact TV reconstructions due to nonlinearity and nonsmoothness of the underlying optimization problem. Many different approaches have been developed, such as time marching [2], fixed-point iteration [3], and various minimization-based methods such as sub-gradient methods [4], second-order cone programming (SOCP) [5], and dualitybased methods [6], [7] -but for large-scale applications such as CT image reconstruction the computational burden is still unacceptable. As a consequence heuristic and much faster techniques such as the one in [1] for approximating the TV solution have been developed. In such inaccurate, but efficient, TV-minimization solvers the resulting image depends on several algorithm parameters, which introduces an unavoidable variability. In contrast, for the accurate TV algorithms considered here, the resulting image can be considered dependent only on the parameters of the optimization problem.
In this work we present two accelerated gradient-based optimization methods that are capable of computing the TV reconstruction of 3D volumes to within an accuracy specified by the user.

A. Total-variation-based image reconstruction
In this work we consider total-variation (TV)-based image reconstruction for computed tomography. The 3D reconstruction is represented by the vector x * which is the solution to the minimization problem Here, x is the unknown image, Q is the set of feasible x, A is the system matrix, b is the projection data stacked into a column vector, and α is the TV regularization parameter specifying the relative weighting between the fidelity term and the TV term. x TV is the discrete total-variation of x, where N is the number of voxels and D j is the forward difference approximation to the gradient at voxel j.

B. Smooth and strongly convex functions
We recall that a continuously differentiable function f is convex if for all x, y ∈ Q. A stronger notion of convexity is strong convexity: f is said to be strongly convex with strong convexity parameter µ if there exists a µ ≥ 0 such that for all x, y ∈ Q. Furthermore, f has Lipschitz continuous for all x, y ∈ Q. The ratio µ/L is important for the convergence rate of gradient methods we will consider. The problem (1) can be shown [8] to be strongly convex and have Lipschitz continuous gradient in the case where A specifies a full-rank overdetermined linear system. In the rank deficient or underdetermined case, which occurs for instance for few-view data, the strong convexity assumption is violated. However, as we shall see, this turns out not to pose a problem for the gradient methods we consider.

A. Gradient projection methods
The optimization problem (1) can, in principle, be solved by use of a simple gradient projection (GP) method (6) where P Q denotes projection onto the set Q of feasible x, and θ k is the step size at the kth step. The worst-case convergence rate of GP with µ > 0 and constant step size is where For large-scale imaging modalities, such as CT, this slow convergence renders the simple gradient method impractical. On the other hand the simplicity and the low memory requirements of the gradient method remain attractive. Various modifications have been suggested in the optimization literature. For instance, a significant acceleration is often observed empirically if the gradient method is equipped with a Barzilai-Borwein (BB) step size strategy and a nonmonotone line search [10], [11], [12], [13], [14], see Algorithm 1: GPBB for a pseudo-code. Empirically we have found K = 2 and σ = 0.1 to be satisfactory parameter choices. However, it remains unproven that GPBB achieves a better worst-case convergence rate than (7).

B. Nesterov's optimal method
Nesterov [15] proposed a gradient-based method that for given µ > 0 achieves the convergence rate where C N is a constant, and he proved the method to be optimal, i.e., that no gradient-based method can achieve better worst-case convergence rate on the class of strongly convex problems. Comparing (7) and (8), we see how the ratio µ/L affects the predicted worst-case convergence rates: When µ/L decreases, both convergence rates become slower, but less in (8) due to the square root. We therefore expect Nesterov's method to show better convergence for smaller µ/L. Small µ/L arise for instance when the number of views is small, see [8].
Nesterov's method requires that both µ and L are given by the user, and in order for the method to be convergent µ must be chosen sufficiently small and L sufficiently large. For real world applications such as CT, µ and L are seldom known, which makes the method impractical. Taking overly conservative estimates can depreciate the better convergence rate (8); hence, accurate estimates of µ and L are important.

C. Estimating µ and L
A sufficiently large L can be chosen using back-tracking line search [16], [17], see Algorithm 2: BT for pseudo-code. Essentially, an estimateL of L is increased by multiplication with a constant ρ L > 1 until (5) is satisfied.
Accurately estimating µ, such that (4) is satisfied globally, is more difficult. Here, we propose a simple and computationally inexpensive heuristic: In the kth iteration choose an estimate µ k as the largest value of µ that satisfies (4) between x (k) and y (k) , and make the µ k -sequence non-increasing: Algorithm 2: BT input : y,L output: x,L 1L ←L ; We call the Nesterov method equipped with estimation of µ and L Unknown Parameter Nesterov (UPN) and pseudo-code is given in Algorithm 3: UPN.
Unfortunately, convergence of UPN is not guaranteed, since the estimate (9) can be too large. However, we have found empirically that an estimate sufficient for convergence is typically effectively determined by (9).
It is possible to ensure convergence by introducing a restarting procedure [8] at the price of lowering the convergence rate bound and thereby losing optimality of the method. However, we have found empirically that the restarting procedure is seldom needed, and for realistic parameters the simple heuristic (9) is sufficient.

D. Stopping criterion
For an unconstrained convex optimization problem such as (1) the norm of the gradient is a measure of closeness to the minimizer through the first-order optimality conditions [18]. For a constrained convex optimization problem it is possible to express a similar optimality condition, namely in terms of the gradient map defined by where ν > 0 is a scalar. A point x is optimal if and only if G ν (x ) = 0 for any ν > 0 [17]. We can use this to design a stopping criterion: Stop the algorithm after iteration k if G ν (x (k) ) 2 /N ≤ , where is a user-specified tolerance.
For an under-determined problem, e.g. in the few-view case, the objective function in (1) is nearly flat at the minimizer, which makes it difficult to determine when a sufficiently accurate reconstruction has been found. Here, the gradient map provides a simple, yet sensitive, stopping criterion.

A. Simulation setup
At this point we emphasize that our objective is to obtain an accurate TV reconstruction while reducing the required number of gradient method iterations. We include two reconstructions merely to demonstrate that the methods indeed are successful in solving (1), thereby reconstructing the desired image. In [8] dependence of the convergence with respect to parameter variation is explored.
To demonstrate and compare the convergence of GP, GPBB and UPN we set up a simple test problem. As test image x true we use the threedimensional FORBILD head phantom discretized into 64 3 voxels. We simulate a parallel beam geometry with view directions evenly distributed over the unit sphere. Projections are computed as the forward mapping of the discretized image subject to additive Gaussian white noise e of relative magnitude e 2 / Ax true 2 = 0.01, i.e., We enforce nonnegativity by taking Q = R 64 3 + . We consider two reconstructions: A "many-view" using 55 views and a "few-view" using only 19 views of size 91 2 pixels. In the latter case A has less rows than columns, which can be shown [8] to lead to violation of the assumption on strong convexity. The iterations are continued to the tolerance = 10 −8 is met. Fig. 1 shows the middle (33rd) axial voxel slice through the original 3D volume together with many-view and few-view UPN reconstructions using α = 0.01. Both reconstructions reproduce the orignal features accurately, except for two small features are missing in the few-view reconstruction. Fig. 2 shows the convergence of the three methods in terms of objective value φ (k) relative to the true minimal objective value φ as function of the iterations k. As φ is unknown, we have approximated it by computing the UPN solution for an two orders of magnitude lower than the value used in the iterations.

B. Simulation results
In both cases we see that UPN converges to a satisfactory accuracy within 2000 iterations, whereas GP does not, and GPBB only does in the former case. In the many-view case UPN and GPBB both produce a significant (and comparable) acceleration over GP. In the few-view case, we also observe acceleration for both, but UPN stands out with much faster convergence. This is in accordance with the expectation stated in Section III-B.
The adequacy of the stopping criterion is evaluated by a simple visual comparison of the few-view simulation gradient map norm decay (Fig. 2 right) and the objective decay (Fig. 2 center). Apart from the erratic decay for GPBB (which is caused by highly irregular step length selection) there is a pronounced correspondence, and we therefore consider the stopping criterion effective.
Although UPN was designed for strongly convex problems, we conclude that the method also works in the non-strongly convex case of having few-view data -in fact, from the preliminary results the non-strongly convex case is where UPN shows its biggest potential by exhibiting a much faster convergence than GP and GPBB.

V. CONCLUSION
We have described the gradient-based optimization methods GPBB and UPN and their worst-case convergence rates. Our simulations show that both algorithms are able to significantly accelerate high-accuracy TV-based CT image reconstruction compared to a simple gradient method. In particular UPN shows much faster convergence when applied to few-view data. The software (implementation in C with an interface to Matlab) is available from http://www2.imm.dtu.dk/ ∼ pch/ TVReg/.