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.

💡 Research Summary

This paper addresses the computational challenge of total‑variation (TV) regularized computed tomography (CT) image reconstruction, especially in three‑dimensional (3‑D) settings where the problem size easily reaches millions of variables. While TV regularization yields high‑quality reconstructions from sparse‑view or limited‑angle data, the associated optimization problem is large, non‑smooth, and often ill‑conditioned, making conventional second‑order methods such as Newton’s method impractical due to prohibitive memory requirements. The authors therefore focus on first‑order gradient‑based algorithms that have modest memory footprints but traditionally converge slowly.

Two accelerated gradient schemes are investigated: (1) GPBB, a gradient‑projection method enhanced with Barzilai‑Borwein (BB) step‑size selection and a non‑monotone line‑search; and (2) UPN (Unknown‑Parameter Nesterov), which adapts Nesterov’s optimal first‑order method to the case where the strong‑convexity constant µ and the Lipschitz constant L are not known a priori. GPBB computes a BB step size at each iteration from the difference of successive iterates and gradients, then reduces a scaling factor β by a factor of 0.95 until a sufficient decrease condition based on the maximum of the last K objective values is satisfied. Although this heuristic dramatically speeds up convergence in practice, the authors note that a rigorous worst‑case rate improvement over the basic gradient method has not been proved.

UPN builds on Nesterov’s method, which for strongly convex problems achieves the optimal rate
 f(x_k) – f* ≤ C·(1 – √(µ/L))^k.
Since µ and L are unavailable in realistic CT problems, the algorithm estimates them on the fly. L is obtained via a back‑tracking line search that repeatedly inflates an estimate until the Lipschitz condition holds. µ is estimated locally between the current iterate x_k and an auxiliary point y_k by finding the largest µ that satisfies the strong‑convexity inequality (4) and then enforcing a non‑increasing sequence of µ values. With these estimates, the algorithm updates the momentum parameter θ and the extrapolation coefficient β according to Nesterov’s recurrence, and constructs the next auxiliary point y_{k+1}. The authors acknowledge that the µ‑estimate may occasionally be too large, potentially breaking convergence guarantees; however, empirical results show that the simple heuristic works well, and a restart mechanism can be added if needed.

A stopping criterion based on the gradient map G_ν(x) = ν