{'fraction_non_alphanumeric': 0.031439254305124016, 'fraction_numerical': 0.012112275527937227, 'mean_word_length': 4.626666666666667, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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}}
Accelerating optimization-based computed tomography via sparse matrix approximations

Variational formulations of reconstruction in computed tomography have the notable drawback of requiring repeated evaluations of both the forward Radon transform and either its adjoint or an approximate inverse transform which are relatively expensive. We look at two methods for reducing the effect of this resulting computational bottleneck via approximating the transform evaluation with sparse matrix multiplications. The first method is applicable for general iterative optimization algorithms. The second is applicable in error-forgetting algorithms such as split Bregman. We demonstrate these approximations significantly reduce the needed computational time needed for the iterative algorithms needed to solve the reconstruction problem while still providing good reconstructions.


Introduction
In the present work, we focus on computed tomographic reconstruction methods based around solving a minimization problem of the general form (1) min for p = 1, 2 and some linear operator K. For simplicity, we focus on the case of parallel beam tomography. Here, R denotes the Radon transform and M is a linear operator, similar to a "row selector", which restricts to those angles along which measurements are taken in constructing f , which is the given data. Standard methods for tomographic reconstruction typically use filtered back projections [10]. In [2,1], the computational costs of such methods were reduced via efficient computations of the Radon transform and its inverse. However, in the case of situations where the angular measurements are under resolved, filtered back projections lead to high frequency artifacts in the reconstructions. Solving Equation (1), in the specific case where p = 1 and K = ∇, provides improved reconstructions without these artifacts, as demonstrated in [3,9]. Such methods have the drawback of having an increased computational cost. One source of this increase is in the need for repeatedly computing the Radon transform and its adjoint. After describing in some detail this computational bottleneck in section 2, we describe two methods for reducing this cost in section 3. The first relies on the fixed geometries of the tomographic reconstruction problem to replace relatively expensive transforms with sparse matrix multiplications. The second method focuses on taking advantage of the errorforgetting nature of algorithms such as the split Bregman [7] algorithm to reduce the computational costs even further. After describing these methods, we present in section 4 the ability of these two methods to reduce the computational costs of solving the total variation-penalized reconstruction problem.

A computational bottleneck in optimization-based reconstruction
We now turn to specifying the computational bottleneck we are interested in addressing. First, if p = 2, the solution to (1) solves the linear operator equation If p = 1, typically an iterative convex optimization algorithm is needed. Such algorithms often involve solving a series of smooth optimization problems. For instance, the primal-dual method of Chambolle and Pock [5] at each iterate will require us to solve smooth problems of the form which means solving Similarly, the split Bregman method of Goldstein and Osher [7] (used specifically for neutron tomography in [3]) requires repeatedly solving a linear system of the form In contrast to the cases considered in [7] and [5], the solutions to (2),(3), and (4) do not have closed form solutions. This means one must resort to iterative linear solvers. Due to the expense of evaluating R * M * M R, a computational bottleneck arises in the iterative linear solvers. In the case of the primal-dual and the split Bregman methods on large-scale problems, this operator may need to be evaluated hundreds, or even thousands of times, necessitating an efficient implementation of its evaluation.
Let M θi denote the restriction associated with a single direction θ i . Via the Fourier slice theorem, for each angle θ i , the one dimensional Fourier transform of M θi Rµ may be associated with a two dimensional Fourier transform of µ at points on a slice in the Fourier domain [4]. Repeating for all angles used in constructing M in Equation (1), we may rewrite Equation (1) as However, in the discrete setting, F * D µ −f is evaluated on the slices corresponding to the discrete angles for which measurements have been taken, where F D denotes a standard FFT. This results in F * D µ being needed on a non-Cartesian grid; we denote by F NU the associated discrete Fourier transform on nonuniform grids. This means that the split Bregman method, for instance, requires solving the linear system (6) (αF NU F * NU − λK * K µ = αF NUf + λK * v. Similar statements hold for the analogues of Equations (2) and (3). Therefore, this modified minimization problem still poses the same computational difficulties as mentioned above, however now the computational bottleneck lies in the evaluation of F NU F * NU .

Overcoming the computational bottleneck
In this section we look at two methods for reducing the computational bottleneck outlined in section 2. The first method involves replacing F * NU F NU with sparse matrix multiplications and standard fast Fourier transforms. This method is applicable to speeding up the conjugate gradient methods needed in solving Equation (5) via any of the already mentioned iterative optimization methods. The second method is specifically applicable to the split Bregman method. Relying on the error-forgetting properties of the algorithm, we see dramatic reductions on the computational times needed without a loss in the overall performance of the algorithm. In each case, we replace evaluating much of the effort of evaluating F NU F * NU with sparse matrix multiplications. From an application perspective, these matrices depend primarily on the resolution of the instruments, the frequency of angular measurements, and parameters which typically take only a few discrete values. This means for a particular tomographic reconstruction workflow, one may precompute the relevant matrices offline, enabling significant speedup in situations where multiple reconstructions are needed.
3.1. Implementation of convolution fusing. In [8], nonuniform fast Fourier transforms (NUFFT) were developed, providing a computationally efficient method of computing these operators. We briefly review this construction. The one-dimensional, type 1 NUFFT of f (x), is computed by taking where τ is a parameter determining the accuracy desired, and g τ is a periodic heat kernel of the form For notational convenience, we let A be the linear operator defined by A(f )(k) := π τ e k 2 τ f k and let Bf := g τ * f (x).
We may then rewrite F NU F * NU as (8) AF D BB * F * D A (noting that A is self-adjoint). As noted in [8], the main computational task in computing the NUFFT is evaluating In order to do this efficiently, a parameter M sp is chosen in order to achieve a desired level of accuracy in the NUFFT. For m ∈ [0, 2N − 1], we approximate Bf by truncated summations, taking the nearest 2M sp points around x j . Outside of this window, the Gaussians decay sufficiently sharply enough to be ignored.
Defining m j = sup{m : πm/M ≤ x j } and I m : It is noted in [8], drawing from [6], that if M sp = 6, τ = 6/N 2 , the NUFFT is accurate to 6 digits of accuracy and if M sp = 12, τ = 12/N 2 , the NUFFT is accurate to 12 digits of accuracy.
We now turn to B * B, which takes the form This operator is sparse; only grid points within M sp of x j are used, along with appropriate weights taken from the Gaussian sources. In higher dimensions, this sparsity structure holds; for instance in two dimensions, the only grid points used in computing B * Bf (x i,j ) are those points such that |m − i| < M sp and |j − n| < M sp . The structure of B * B is determined by the radial resolution of the data f , the set of angular measurements used, and the value M sp . The first depends on the accuracy of the instrumentation, which typically is fixed. In situations where the set of angular measurements is either fixed, or from one of several canonical sampling regimes, the B * B may be precomputed as a sparse matrix. The resulting sparse matrix fuses the convolutions into a single operation. Thus, the NUFFT's needed in solving Equation (6) may be replaced by 2 FFT's, a single sparse matrix-vector multiplication, and 2 diagonal matrixvector multiplications.
3.2. Acceleration via surrogate matrices. Taking advantage of fusing the convolutions as outline above in section 3.1 means, at each iterate in the split Bregman algorithm, we now solve where now B * B has been replaced with a sparse matrix multiplication. However, the (split) Bregman algorithm has an error forgetting property [7]; that is, one only needs in general to solve Equation (9) approximately in order to still have a convergent algorithm. The benefit for doing exact Bregman updates lies in the ability of Bregman to converge in two exact updates, for sufficiently small initial Bregman distance, as discussed in [12]. However, we note that solving Equation (9) exactly (to machine precision) can be quite expensive. For large-scale imaging problems, the expense of exactly solving Equation (9) may offset the near immediate convergence of Bregman updates for good (with respect to the Bregman distance) initialization. Indeed, in [3], inexact Bregman updates were sufficient for a similar formulation of the tomographic reconstruction problem. In order to further reduce the computational costs associated with solving Equation (9), we look to leverage the error-forgetting nature of the split Bregman algorithm by precomputing a sparse matrix T which approximates AF D (B * B)F * D A. We first consider the cases where u is the simple N × N image given by The result of F * A being applied to this u is constant function. As described above in section 3.1, B * BF * D A(δ i,N/2 δ j,N/2 ) will then be the sum of Gaussians. Therefore, AF D (B * B)F * D A (u) will be a (possibly rescaled) Gaussian centered at (N/2, N/2). This is shown in Figure 1a for N = 256 where we show the central 128 × 128 region when M corresponds to a Shannon sampling frequency in angle. We will use the peakedness seen in that image to construct a T . First, for a chosen We then define as the sum of these translated stencil matrices. This T is used to compute split Bregman updates; instead of Equation (9), we solve (possibly approximately) It is important to recall that for an an image which is a single non-zero pixel near the image's edge, e.g.ũ i,j = δ i,k δ j, for some k < r, AF D (B * B)F * D A ũ will have nontrivial entries near the opposite edge, due to the periodicity of the operators. Our choice of T suppresses these values on the opposite edge. The motivation for this choice is to induce a banded structure in T , as seen in Figures 1b and 1c; in this case r = 3, resulting in a banded matrix with 7 diagonal and off-diagonal entries per row. However, this is reasonable in the case of computed tomography where the sample to be measured is fully contained in the interior of the image. The choice of suppressing values on the opposite edge corresponds to an assumption that the object of interest in the reconstructed image is at least r pixels from each edge of the image and that near the edge, the true µ can be assumed to be 0.

Numerical Implementation
In order to test the efficiency of using the sparse linear matrix representation, we solve Equation (1) for K, p corresponding to anisotropic TV regularization via a split-Bregman algorithm (as discussed for the general problem in [7] and specifically for computed tomography in [11,3]).
Throughout, we use the Shepp-Logan phantom as the reference data and α = λ = 1. All computations reported below were performed in MATLAB on a workstation using Intel 2.20 GHz Xeon Processors with 55 MB caches.
First, in order to compare the cost and accuracy of using the NUFFT as opposed to either the convolution fusing method or the use of sparse surrogates, we look at the time required to compute a fixed number of inexact Bregman iterations. We compute 200 inexact Bregman iterations; each of these iterations involves approximating the solution to Equation (4) via 5 conjugate gradient steps. Therefore, 1000 evaluations of αR * M * M R are needed. We perform this for various image sizes as well as various values of M sp used in computing the NUFFT and its corresponding B. For each value of N x , we choose M that corresponds to using N x π 4 evenly spaced angular measurements. We summarize the results in Tables 1 and 2 as well as Figure 2. First, in Table 1, we report the time required for performing 200 inexact Bregman iterations as well the portion of that time required for performing the conjugate gradient computations. As expected, increasing M sp increases the time required when using the NUFFT and convolution fusing-based optimizations. The sparsity structure of T in the surrogate-based method is dominated by the choice of r (here set to be 1), resulting in the time needed being unchanged for different M sp . Additionally, the time required per evaluation of αR * M * M R using the NUFFT is, for M sp = 2, roughly 3 or 4 times longer than using convolution fusing and 10-20 times longer than using sparse surrogates. At higher values of M sp , the speedup using convolution fusing seems to hold; meanwhile the possible speedup by using sparse surrogates becomes one of several orders. However, the speedup using sparse surrogates is tempered by a significant loss of accuracy in the reconstruction. This can be seen in Table 2, where we report the relative L 1 error of the reconstruction using these 200 split Bregman iterations. The sparse surrogate-based reconstruction results in a relative error that is between 2 and 5 times higher than reconstructions using the other methods. One perhaps surprising result from this test is that the choice of M sp appears to not significantly affect the quality of the reconstruction after a fixed number of iterations. Finally, in Figure 2, we display the relative errors and the L 1 norm of the updates over the iterations. Using either the NUFFT or convolution fusing results in nearly the same iterative reconstruction. Additionally, despite the higher error in the sparse surrogate method, the convergence rates are unchanged. Indeed, the sparse surrogate method leads, after a handful of iterations, to a iteration scheme which is largely monotonically reducing the error of the reconstruction, unlike the other  (6), N x = 512 and 6 digits of accuracy for NUFFT. On the left is the L 1 error at each update; on the right, the L 1 norm of the Bregman updates.
two methods. These observations suggest that for a fixed number of iterations, convolution fusing provides a significant speedup over the NUFFT while not impacting the performance of the image reconstruction procedure. Meanwhile, the use of sparse surrogates decreases the accuracy at any given iterate. This also holds if one increases r: for r = 3, 200 Bregman updates require 136.0654 seconds results in a relative L 1 error of 0.129994 error. As expected, then, increasing r results in more expensive (due to decreasing the sparsity of the resulting T ) albeit more accurate iterations. However, due to the shorter time of each iteration in the sparse surrogate method, there is a possibility that the overall performance of solving (1) can be improved: we may replace relatively accurate iterations with more iterations which, individually, are significantly quicker to perform.  To test this, we look to solve (5) for 512×512 images and K, p again corresponding to anisotropic total variation regularization. We consider M corresponding to taking 512 π 4d angular measurements for d = 1, 2, 4 and T corresponding to r = 3. That is, we consider reconstruction problems where the angular measurements are under-resolved to differing degrees. In [3], the highly under-resolved case was of interest, as angular measurements are relatively expensive in some applications. We will compare the results of using split Bregman to solve (5) with both convolution fusing and with sparse surrogates. That is, we will compare using (9) with using (11) to compute the Bregman update; again, we use at most 5 conjugate gradient steps to compute this update. As a stopping criterion, we terminate either after a maximum of 6000 Bregman updates or if the L 1 norm of a Bregman update is less than 10 −8 of the first Bregman update. In Tables 3 to 5 we list the results of using each method to compute the Bregman update for d = 1, 2, 4, respectively. In each we show the time and number of iterations needed for Bregman updates to fall below a certain threshold, as well as the relative L 1 error of the reconstruction at that iteration. Additionally we show in We see that for any fixed time budget, the sparse surrogate-based method provides improved accuracy in very few iterations. In the case of d = 1, the time needed to get 10 −5 order relative L 1 error via the sparse surrogate method is approximately 10% of the time needed to reach 10 −4     L 1 error via convolution fusing. Similar observations hold for d = 2. However, for d = 4, the sparse surrogate-based method, while taking vastly less time than the convolution fusing-based method, results in a final result which is more inaccurate. This is due to prematurely terminating the reconstruction method, as the optimization routine detects a convergence due to the relatively small size of the Bregman updates. Nonetheless, the time needed to reach 3.4739 × 10 −5 error via convolution fusing is 5104 seconds, versus the 683 seconds needed for sparse surrogates.

Conclusions
Approximating the transform F * F with sparse matrices reduces the computational bottleneck that arises in solving variational tomographic reconstruction problems. Convolution fusing replaces the most computationally expensive portion of the evaluation of this operator with a sparse matrix multiplication, reducing the computation time without any additional loss of accuracy. In cases where an error-forgetting algorithm is to be used, sparse surrogates further reduce the cost of the computation by replacing the entire operator with a sparse matrix; this comes at the cost of a loss of accuracy. However, the reduced time needed per iteration offsets the need for additional iterations, especially in the case of only moderately unresolved angular measurements. In the case of low angular resolution reconstructions where the spares surrogates result in prematurely terminating algorithms, it may prove useful as a warm-start strategy. That is, as we expect sparse surrogate-based iterations may prematurely converge, we use the result as an initial guess for fused convolution-based iterations which are more expensive.