Jeroen van Gemert^{1}, Kirsten Koolstra^{2}, Peter BĂ¶rnert^{3}, Andrew Webb^{2}, and Rob Remis^{1}

Reconstruction methods in parallel imaging and compressed sensing problems are generally very time consuming, especially for a large number of coil elements. In this work, the image is reconstructed using the Split Bregman algorithm (SB). We present an efficient and effective preconditioner that reduces the number of iterations in the linear least squares step of SB by almost a factor of 5 as alternative to extra variable splitting. The designed preconditioner works for Cartesian sampling schemes and for different coil configurations. It has negligible initialization time and leads to an overall speedup factor of 2.5.

For each SB iteration a linear least squares problem of the form $$$\text{A}\textbf{x}=\textbf{b}$$$ has to be solved, where $$$\textbf{x}$$$ is the true image, $$$\textbf{b}$$$ is the weighted undersampled k-space data, and

$$\text{A}=\mu\sum_{i=1}^{\text{N}_\text{c}}\left(\text{R}\text{F}\text{S}_i\right)^H\text{R}\text{F}\text{S}_i +\lambda\left(\text{D}_x^H\text{D}_x+\text{D}_y^H\text{D}_y\right)+\gamma\text{W}^H\text{W}\qquad(1)$$

is the system matrix that remains constant throughout the SB iterations. Here, the $$$\text{S}_i $$$ are diagonal matrices representing complex coil sensitivity maps for each of the $$$\text{N}_\text{c}$$$ channels. Furthermore, $$$\text{F}$$$ is the discrete two-dimensional Fourier transform matrix, and the undersampling pattern is specified by the binary diagonal sampling matrix $$$\text{R}$$$. The regularization parameters $$$\lambda$$$ and $$$\gamma$$$ in (1) are promoting sparsity via the total variation constraint (employing differentiation matrices $$$\text{D}_x$$$ and $$$\text{D}_x$$$) and the unitary wavelet transform $$$\text{W}$$$, respectively.

The second and third term in (1) are block circulant with circulant blocks (BCCB) matrices and are diagonalized by the Fourier matrix $$$\text{F}$$$. Writing

$$\text{A}=\text{F}^H\text{F}\text{A}\text{F}^H\text{F}=\text{F}^H\text{K}\text{F}$$

we have

$$\text{K}=\mu\underbrace{\sum_{i=1}^{\text{N}_\text{c}}\text{F} \text{S}_i^H\text{F}^H\text{R}^H\text{R}\text{F}\text{S}_i \text{F}^H}_{\text{K}_c} + \lambda \underbrace{ \vphantom{\sum_{i=1}^{\text{N}\text{c}}}\text{F} \left( \text{D}_x^H\text{D}_x + \text{D}_y^H\text{D}_y \right)\text{F}^H}_{\text{K}_d}+ \gamma \underbrace{\vphantom{\sum_{i=1}^{\text{N}\text{c}}}\text{I}}_{\text{K}_w}$$

where $$$\text{K}_d$$$ and $$$\text{K}_w$$$ are diagonal, but $$$\text{K}_c$$$ is not. Most of the
energy in $$$\text{K}_c$$$ is located
around its main diagonal, however. Hence, we approximate the matrix $$$\text{K}_c$$$ by its diagonal to arrive
at a BCCB approximation of $$$\text{A}$$$ as
preconditioner for solving the linear system. The required steps can be efficiently
performed using FFTs.

*MR Data Acquisition*: A Philips Ingenia
3T dual transmit MR system was used to acquire *in vivo* data. A 12-element posterior receiver array and 15-channel
head coil were used for reception in spine and brain, respectively, whereas the
body coil was used for RF transmission. T_{1}-weighted images were
acquired using a standard clinical TSE sequence and IR-TSE sequence for spine
and brain, respectively.

*Reconstruction*: The SB
algorithm was implemented in MATLAB on a standard Windows 64-bit PC with 8 GB
of internal memory. Unprocessed k-space data was stored per channel and used to
construct cropped complex coil sensitivity maps^{14}. Both the individual coil images and the coil
sensitivity maps were normalized. Reconstruction for the spine data set was
performed without and with coil compression^{15,16}. A random line pattern and a fully random pattern, both with variable density,
were studied for undersampling factors four (*R*=4) and eight (*R*=8).

Figure 1 shows the TSE spine images for
Cartesian undersampling factors *R*=4 and *R*=8. The full system matrix is
constructed for illustration purposes only and compared with its circulant approximation
in Fig. 2.

Figure 3a shows the number of iterations required for the preconditioned conjugate-gradient (PCG) solver to converge in each SB iteration without preconditioner, with the Jacobi preconditioner (the diagonal of matrix $$$\text{A}$$$), and with the circulant preconditioner proposed in this abstract. The reduced number of iterations leads to a shortened reconstruction time, plotted in Fig. 4a for the total PCG part in the SB algorithm (4.65x), and in Fig. 4b for the entire reconstruction algorithm (2.5x). For all matrix sizes the initialization time of the preconditioner is less than 2% of the image reconstruction time. Similar results are obtained for the 15-channel head coil as shown in Fig. 5.

We would like to thank Ad Moerland from Philips Healthcare Best (The Netherlands) and Mariya Doneva from Philips Research Hamburg (Germany) for helpful discussions on reconstruction.

1. Pruessmann K et al. SENSE: Sensitivity Encoding for Fast MRI. Magnetic Resonance in Medicine, 1999; 42(5):952-962.

2. Griswold M et al. Generalized Autocalibrating Partially Parallel Acquisitions (GRAPPA). Magnetic Resonance in Medicine. 2002; 47(6):1202-1210.

3. Blaimer M et al. SMASH, SENSE, PILS, GRAPPA: How to Choose the Optimal Method. Topics in Magnetic Resonance Imaging. 2004; 15(4):223-236.

4. Deshmane A et al. Parallel MR imaging. Journal of Magnetic Resonance Imaging. 2012; 36(1):55-72

5. Donoho D et al. Compressed Sensing. IEEE Transactions on Information Theory. 2006; 52(4):1289-1306.

6. Lustig M et al. Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging. Magnetic Resonance in Medicine. 2007; 58(6):1182-1195

7. Goldstein T et al. The split Bregman Method for L1-Regularized Problems. SIAM Journal on Imaging Sciences. 2009; 2(2):323-343.

8. Ramani A et al. Parallel MR Image Reconstruction using Augmented Lagrangian Methods. IEEE Trans Med Imaging. 2011;30(3):694-706.

9. Weller D et al. Augmented Lagrangian with Variable Splitting for Faster Non-Cartesian L1-SPIRiT MR Image Reconstruction. IEEE Trans Med Imaging. 2014;33(2):351-361.

10. Saad Y. Iterative methods for sparse linear systems. Society for Industrial and Applied Mathematics. 2003; 2.

11. Chen C et al. Real Time Dynamic MRI with Dynamic Total Variation. International Conference on Medical Image Computing and Computer-Assisted Intervention. 2014;138-145

12. Li R et al. Fast Preconditioning for Accelerated Multi-contrast MRI Reconstruction. International Conference on Medical Image Computing and Computer-Assisted Intervention. 2015,700-707.

13. Xu Z et al. Efficient Preconditioning in Joint Total Variation Regularized Parallel MRI Reconstruction. International Conference on Medical Image Computing and Computer- Assisted Intervention. 2015;563-570.

14. Uecker M et al. ESPIRiT - An Eigenvalue Approach to Autocalibrating Parallel MRI: Where SENSE Meets GRAPPA. Magnetic Resonance in Medicine. 2014; 71(3):990-1001.

15. Huang F et al. A Software channel compression technique for faster reconstruction with many channels. Magnetic Resonance Imaging. 2008; 26(1):133-141.

16. Zhang T et al. Coil Compression for Accelerated Imaging with Cartesian Sampling. Magnetic Resonance in Medicine. 2013; 69(2):571-582.

Figure 1: Reconstruction results for a T_{1}-weighted spine scan using different Cartesian
undersampling factors. (a) shows the fully sampled scan as a reference, whereas
(b) and (c) depict the reconstruction results for undersampling factors four
(*R*=4) and eight (*R*=8), respectively. The absolute difference is shown in (d)
and (e) for *R*=4 and *R*=8, respectively. The reconstruction matrix is 512x512.
Regularization parameters were set to µ =1·10^{-3}, λ=4·10^{-3}, and γ=1·10^{-3}. Images were
acquired using a TSE sequence with FOV=340x340 mm^{2}; in-plane
resolution 0.66x0.66 mm^{2}; slice thickness = 4 mm, TE/ TR/TSE factor
= 8 ms/ 648 ms/ 8.

Figure 2: System matrix and its circulant approximation. The first two columns and the last two columns
show the system matrix elements for random Phase-Encoding (PE) (a) and fully random (b) undersampling, respectively. The top
row depicts the elementwise magnitude and phase of the true system matrix $$$\text{A}$$$,
whereas the bottom row depicts the elementwise circulant approximation. The zeros in $$$\text{A}$$$ are not present in the circulant approximation, since the circulant property is enforced. The newly introduced entries do not add relevant information as the image vector on which the system matrix acts contains zero signal in the corresponding region.

Figure 3: Number of PCG iterations required in each Split Bregman
iteration. The circulant preconditioner reduces the number of iterations
considerably compared with the non-preconditioned case. The Jacobi
preconditioner hardly reduces the number of PCG iterations. (a) depicts the number
of PCG iterations for Set 1: µ = 1·10^{-3}, λ=4·10^{-3}, γ=1·10^{-3},
whereas (b) depicts the number of PCG iterations for Set 1, Set 2: µ =1·10^{-2},
λ=4·10^{-3}, γ=1·10^{-3}, and Set 3: µ=1·10^{-3}, λ=4·10^{-3},
γ= 4·10^{-3}. The preconditioner shows the largest speed up factor when
the regularization parameters are well-balanced.

Figure 4: Computation time for 20 Split Bregman iterations
and different problem sizes. (a) Using the preconditioner, the total
computation time for the PCG part in 20 Bregman iterations is reduced by a
factor larger than 4.5 for all studied problem sizes. (b) The computation time
for 20 Split Bregman iterations of the complete algorithm (including Split
Bregman updating steps) is reduced approximately by a factor of 2.5 for the
considered problem sizes.

Figure 5: Reconstruction results for a T_{1}-weighted
brain scan. (a) shows a fully sampled scan as a reference and (b) the reconstruction
results for an undersampling factor *R*=4. The absolute difference is shown in
(c). The reconstruction matrix has dimensions 256x256 and regularization
parameters where chosen as µ = 1·10^{-3}, λ=4·10^{-3}, γ=2·10^{-3}.
The convergence results for the PCG part with and without preconditioner are
plotted in (d). Images were acquired using an IR TSE sequence with FOV =
230x230 mm^{2}; in-plane resolution 0.90x0.90 mm^{2}; slice
thickness = 4 mm; TE/TR/TSE factor = 20ms/2000ms/8; refocusing angle = 120°; IR
delay: 800ms.