3616 A Preconditioned ADMM Strategy for Field-Corrected Non-Cartesian MRI Reconstruction

Target Audience: Image reconstruction scientists and users of non-Cartesian k-space trajectories. Purpose: Promoting sparsity during image reconstruction [1] has emerged as a powerful technique for improving signal-to-noise ratio (SNR) and suppressing undersampling artifact in accelerated magnetic resonance imaging (MRI) (see Fig 1). However, sparse reconstruction of non-Cartesian MRI data remains computationally challenging since multiple “gridding” operations must be executed at each iteration of the reconstruction. Recently, Ramani and Fessler [2] proposed an efficient alternating-direction-method-of-multiplier (ADMM) strategy for sparse MRI reconstruction that operates by decomposing the optimization problem into a series of (relatively) easier sub-problems that can be solved serially. For non-Cartesian MRI, the data fidelity sub-problem must also be solved iteratively (e.g., via conjugate gradient (CG) iteration). If off-resonance effects are accounted for in the MRI signal model (e.g., to avoid geometric distortion), standard circulant preconditioners (PC) [2] cannot be used to accelerate this task since the Gram matrix constructed from the sampling operator is not Toeplitz [3]. This challenges the practicality of ADMM for non-Cartesian MRI. In this work, we show that an algebraic reformulation of the ADMM scheme enables the use of simple but effective diagonal PCs for non-Toeplitz models, and demonstrate their practical benefit for undersampled SWIRLS 3D contrast-enhanced MR angiography (CE-MRA) [4]. Methods: Presuming a Dirac delta voxel model, non-Cartesian MRI acquisition can be modeled as , where is the set of measured signals, is a Fourier-based sampling operator that includes off-resonance effects, is the set of target images, and is complex Gaussian noise. Sparse MRI reconstruction comprises solving the following problem: argmin ≜ ‖ ‖ , where ∙ is a sparsity penalty. ADMM performs an alternating minimization of the augmented Lagrangian functional, , , ‖ ‖ ‖ ‖ , where >0 and is a Lagrange multipler vector. The original ADMM scheme [2] for this objective is given in Fig. 2a. For standard sparsity penalties, the second sub-problem can be directly solved (e.g., via softthresholding); and the Lagrange multiplier update is trivial. Although the first sub-problem is linear, it must be solved iteratively due to its size and structure. Since ∗ is not Toeplitz [3], how to effectively PC this system is not obvious. However, note via Woodbury’s identity [5] that ∗ ∗ ∗ ∗ . This allows the first sub-problem to be further decomposed as shown in Fig. 2b. Now, recall that non-iterative “gridding” reconstruction is implemented as ∗ , where ∗ is a positive definite diagonal matrix termed the “density compensation function” (DCF) [6]. This suggests, for the first subproblem in Fig. 2b, use of diagonal PCs of the form . Unlike previous DCF-based PC strategies (e.g., [7]), this approach does not alter our statistical objective and therefore imposes no SNR penalty. To evaluate this strategy, the convergence rate of several algorithms was compared while reconstructing a 3.35x undersampled SWIRLS 3D CE-MRA [4] data set (240x240x240, 8-channel, field map from separate scan). was defined with a joint sparsity penalty [8] ( ‖Ψ ‖ , , Ψ=4-tap 3D Daubechies wavelet [9]) and was implemented as a time-segmented (L=8, Hann window) non-uniform fast Fourier transform (NUFFT) [10] (OF=1.25x, W=5 Kaiser-Bessel kernel). Standard ADMM, Woodbury ADMM w/o and w/ PC, and FISTA [11,12] were implemented in C++/OpenMP and executed with 0.1 for 100 iterations on a dual 8-core 2.6 GHz machine with 128 GB memory. At each ADMM iteration (all cases), the first sub-problem was solved (inexactly) via a single warm-started CG iteration, and was set to yield condition number ∗ 100 [2]. was designed with the method in [4]. Results: Plots of the optimization objective, , versus iteration number and processor (CPU) time are shown in Fig. 3. After termination, all algorithms yielded results visually comparable to the rightmost image in Fig. 1. As expected, standard and Woodbury ADMM w/o PC exhibited similar per-iteration behavior. Although the per-iteration cost of the latter is higher due to more NUFFT operations, this barrier is overcome when DCF-based PCs are used. Woodbury ADMM w/ PC exhibits a 10x improvement in runtime performance over its non-PC analog, 3.5x over standard ADMM, and 2.5x over FISTA. Discussion: The proposed ADMM strategy with DCF-based PC is a promising approach for reducing the computational burden of sparse non-Cartesian MRI reconstruction (including off-resonance correction). Further improvements can potentially be realized by code optimization or use of custom preconditioners. It should also be straightforward to incorporate the proposed methodology into other non-Cartesian MRI reconstruction models that include parallel imaging and/or other regularization strategies (e.g., low-rank). Conclusion: Reformulating the standard ADMM algorithm for non-Cartesian MRI reconstruction with off-resonance correction via Woodbury’s identity enables the use of simple but effective diagonal preconditioners for non-Toeplitz linear sub-problems, and results in improved overall reconstruction efficiency. References: [1] Lustig et al., MRM 58:1182-95, 2007; [2] Ramani and Fessler, IEEE TMI 30:694-706, 2011; [3] Fessler et al., IEEE TSP 53:3393-402, 2005; [4] Shu et al., ISMRM 2011:2656; [5] Golub and Van Loan, Matrix Computations, 1996; [6] Pipe and Menon, MRM 41:179-86, 1999; [7] Pruessmann et al., MRM 46:638-51, 2001; [8] Vasanawala et al., Radiology 256:607-16, 2010; [9] Daubechies and Sweldens, J Fourier Anal Appl 4:247-69, 1998; [10] Sutton et al., IEEE TMI 22:178-88, 2003; [11] Beck and Teboulle, SIAM Imag Sci 2:183-202, 2009; [12] Trzasko et al., RSNA 2012:SSE23-04. FIG 3: Optimization algorithm convergence plots. FIG 1: Thin-slab MIPs of gridding and sparse reconstructions (ADMM-Woodbury-PC) of SWIRLS 3D CE-MRA.