University of Southern California Dissertations and Theses

Optimization methods and algorithms for constrained magnetic resonance imaging

(USC Thesis Other)

Optimization methods and algorithms for constrained magnetic resonance imaging

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Optimization Methods and Algorithms for Constrained Magnetic Resonance Imaging by Yunsong Liu A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ELECTRICAL ENGINEERING) December 2022 Copyright 2022 Yunsong Liu Acknowledgments I would like to express my gratitude towards many people who have provided inestimable support for me along with my Ph.D. journey all the way to nish this dissertation. During this journey, there are ups and downs and the Covid pandemic denitely makes things tougher. Without these support, it would be impossible to complete this dissertation. First, I want to greatly appreciate the support and mentorship from my Ph.D. advisor, Professor Justin P. Haldar. It has always been enlightening and fun to discuss research topics and random chat with Justin. I am also grateful to Justin for providing invaluable guidance on being a rigorous researcher, a critical thinker and an ecient communicator. I am grateful to my doctoral committee, Professors Richard M. Leahy and John C. Wood for their valuable support, advices and questions. I am also grateful to my qualifying exam committee, Professors Krishna Nayak and Mahdi Soltanolkotabi for their precious time and insightful questions. I am thankful for my labmates at the Biomedical Imaging Group (BIG) - especially Divya Varadarajan, Daeun Kim, Tae Hyung Kim, Rodrigo Lobos, Jiayang Wang, Chin-Cheng Chan, Hao- Ting Kung, Debdut Mandal, Anand A. Joshi, Takfarinas Medani, Haleh Akrami, Yijun Liu, Clio Gonzlez Zacaras, Wenhui Cui and Omar Zamzam for their help, encouragement and inspirations. I am also grateful to Professor Kawin Setsompop and Dr. Congyu Liao from Stanford University for their great support and expertise in our extensive collaborations in most of the work in this dissertation. The discussion with Kawin and Congyu has added to my Ph.D. journey with much fun and great inspiration. I would like to acknowledge generous nancial support of my graduate studies from the Ming Hsieh Department of Electrical and Computer Engineering and the University of Southern Califor- nia. ii Lastly, I owe a debt of gratitude to my parents and big brother for their unlimited love and support. iii Contents Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Main Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Organization of the Document . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Chapter 2: Phase Constraints Enable Reduced RF Encodings with applications to gSlider-based Diusion MRI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Combining gSlider with Phase Constraints . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2.1 The Standard (Unconstrained) gSlider Model . . . . . . . . . . . . . . . . . . 7 2.2.2 A Theoretical Analysis of Phase Constraints . . . . . . . . . . . . . . . . . . 9 2.2.3 Optimization of G . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.4 Theoretical comparison of dierent RF-encoding schemes . . . . . . . . . . . 15 2.3 Evaluation using Simulated Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Results with 4 RF Encodings . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3.2 Results with 3 RF Encodings . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.1 Practical Implementation Issues . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.2 Beyond gSlider . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 iv Chapter 3: PALMNUT: An Enhanced Proximal Alternating Linearized Mini- mization Algorithm with Application to Separate Regularization of Magnitude and Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Regularized Image Reconstruction with Separate Magnitude and Phase Modeling . . 25 3.2 Alternating Minimization and Proximal Alternating Linearized Minimization Methods 28 3.2.1 Alternating Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2.2 Proximal Alternating Linearized Minimization (PALM) . . . . . . . . . . . . 29 3.3 PALMNUT: A Fast Algorithm for Separate Magnitude and Phase Image Recon- structions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.1 Applying PALM to Eq. (3.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3.1.1 Huber-function regularization of m with Tikhonov regularization of q 37 3.3.1.2 ` 1 regularization of m with Huber-function regularization of q . . . 39 3.3.2 Uncoupled Step Sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.3 Nesterov's Momentum Acceleration and PALMNUT . . . . . . . . . . . . . . 45 3.4 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4.1 Undersampled MRI Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . 48 3.4.2 Regularization-based MRI Denoising . . . . . . . . . . . . . . . . . . . . . . . 50 3.4.3 Phase-Corrected MR Image Combination . . . . . . . . . . . . . . . . . . . . 52 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Chapter 4: An Ecient Algorithm for Spatial-Spectral Partial Volume Compart- ment Mapping with Applications to Multicomponent Diusion and Relaxation MRI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2 Problem Formulation and Existing Methods . . . . . . . . . . . . . . . . . . . . . . . 61 4.2.1 Formal Description of Optimization Problem . . . . . . . . . . . . . . . . . . 61 4.2.2 Classical Voxel-By-Voxel Solutions . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2.3 Spatially-Regularized Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 63 v 4.3 Proposed LADMM-Based Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3.1 Review of LADMM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3.2 Proposed Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.3.3 Solving the f subproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3.4 Solving the z subproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.3.5 Complexity Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3.6 Practical Parameter Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3.6.1 The parameter p . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3.6.2 The penalty parameter . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3.6.3 The rank parameter r . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4.1 Diusion-Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.4.2 IR-MSE Relaxation-Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.4.3 Multi Spin Echo T 2 Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4.4 MR Fingerprinting Relaxation-Relaxation . . . . . . . . . . . . . . . . . . . . 79 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Chapter 5: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Appendix A: Proof of Thm. 3.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 Appendix B: Derivation of Eq. (2.6) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 vi List of Figures 2.1 SNR eciency versus the number of neighborhood voxels (P ) sharing the same phase. Conventional 5 RF gSlider encodings case is chosen as the baseline. Note that the number of subslices to be resolved is S = 5. . . . . . . . . . . . . . . . . . . . . . . . 16 2.2 RF encodings comparison. The optimized RF encodings are obtained from the pro- posed CRLB optimization with parameters P = 25 and = 1:3. . . . . . . . . . . . 17 2.3 Ground truth images used in the simulations. The rst row shows an example slab of magnitude images with 5 sub slices. (f) shows an example of coronal view magnitude image. (g)-(j) show an example of background phase images for each RF encodings where (h)-(j) are obtained by adding a random shift to (g) respectively. . . . . . . . 18 2.4 Reconstructed magnitude images using 4 RF encodings in the simulation without adding noise. (a) is the ground truth. (b) and (c) are minimum norm and smooth phase constrained reconstructions respectively using 4 conventional RF encodings. (d) and (e) are minimum norm and smooth phase constrained reconstructions re- spectively using 4 optimized RF encodings. (f)-(i) are errors of (b)-(e) relative to (a) respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Phase constrained reconstruction in axial view using 4 optimized RF encodings. (a) is the noise contaminated phase in the RF encoded data (d = Gme i +). (b) is the ground truth phase generated by the complex-valued encoding scheme. (c) is the ground truth background phase. (d)-(f) are the reconstructed total, encoding and background phase images respectively. . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Reconstructed magnitude images using 3 RF encodings in the simulation without adding noise. (a) is the ground truth. (b) and (c) are minimum norm and smooth phase constrained reconstructions respectively using 3 conventional RF encodings. (d) and (e) are minimum norm and smooth phase constrained reconstructions re- spectively using 3 optimized RF encodings. (f)-(i) are errors of (b)-(e) relative to (a) respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.7 Reconstructed phase images using 3 optimized RF encodings in the simulation with added Gaussian white noise of standard deviation 0.005. (a) is the noise contami- nated phase in the RF encoded data (d = Gme i +). (b) is the ground truth phase generated by the complex-valued encoding scheme. (c) is the ground truth background phase. (d)-(f) are the reconstructed total, encoding and background phase images respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 vii 2.8 One example of RF pulse design results. The RF pulses are designed by solving the optimization problem in Eq. (2.22) with parameters = 10 6 and = 5 10 6 . Time-bandwidth product of the excitation RF pulses is 12 with time duration of 11 ms. Slice thickness is 1mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.9 SNR eciency comparisons of practical and ideal RF encodings as shown in Fig. 2.8. 24 3.1 The ground truth (a) magnitude and (b) phase images used for the undersampled MRI reconstruction scenario, along with (c) the 8-accelerated k-space sampling mask used to simulate the undersampled acquisition. . . . . . . . . . . . . . . . . . . 48 3.2 Convergence plots and reconstructed images for the undersampled MRI reconstruc- tion scenario. The convergence plots show (a) the cost function value as a function of computation time and (b) the NRMSE value for the complex image as a function of computation time. Also shown are the magnitude and phase images corresponding to (c) zero-lled reconstruction and (d) PALMNUT. . . . . . . . . . . . . . . . . . . 49 3.3 The (a) ground truth and (b) noisy magnitude and phase images used for the regularization-based MRI denoising scenario. . . . . . . . . . . . . . . . . . . . . . . 50 3.4 Convergence plots and reconstructed images for the regularization-based MRI de- noising scenario. The convergence plots show (a) the cost function value as a func- tion of computation time and (b) the NRMSE value for the complex image as a function of computation time. Also shown are the (c) magnitude and phase images corresponding to PALMNUT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.5 (a) The ground truth magnitude and (b) a representative noisy magnitude image, along with (c) the four dierent ground truth phase images used for the scenario with phase-corrected combination of multiple images. . . . . . . . . . . . . . . . . . 52 3.6 Convergence plots and reconstructed images for the phase-corrected MRI image com- bination scenario. The convergence plots show (a) the cost function value as a func- tion of computation time and (b) the NRMSE value as a function of computation time. Also shown is (c) the combined magnitude image obtained from PALMNUT. . 53 3.7 The convergence plots show the cost function value and the NRMSE value as a function of computation time for (a) the undersampled MRI reconstruction scenario and (b) the regularization-based MRI denoising scenario. . . . . . . . . . . . . . . . 55 4.1 Component spectrum locations (a),(f) and spatial maps (b)-(e), (g)-(j) in the DR-CSI mouse spinal cord experiments. Consistent with [1], three components were identied where components 1-3 seem to correspond to white matter structures, gray matter structures and dorsal gray matter dominated structures respectively. The rst and second rows show results with and without spatial regularization respectively where we can see spatial regularization greatly improves noise suppression and results in more continuous spatial maps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 viii 4.2 Convergence curves in the DR-CSI experiments. The proposed LADMM converges 4x-5x faster than the previous ADMM in terms of cost function values. Even more signicant acceleration was observered in terms of dNRMSEs. Also noticeably, the ill-conditioning of the inverse problem in this case causes the dNRMSEs be ever changing without being convergent for both algorithms even after a long running time. 76 4.3 Component spatial maps comparison in the IR-MSE based RR-CSI experiments. No noticeable dierences were observed between the two compared algorithms. As interpreted in [2], component 1 to 6 seem to correspond to white matter structures, gray matter structures with relatively high myelin content, brain structures with high iron content, gray matter without myelin, cerebrospinal uid (CSF) and myelin water respectively. The corresponding spectrum locations of each component are shown in Fig. 4.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4 Component spectrum overlay comparisons in the IR-MSE based RR-CSI experi- ments. C1-C6 stand for components 1 to 6 respectively. No noticeable dierences were observed between the two compared algorithms. The corresponding spatial maps of each component are shown in Fig. 4.3. . . . . . . . . . . . . . . . . . . . . . 78 4.5 Component spectrum overlay (a) and spatial maps (b)-(e) in the multi spin echo based T 2 -relaxation only experiments. Only three components were identied in this experiment which is in stark comparison to the IR-MSE based two dimensional RR-CSI experiment where six components were captured. . . . . . . . . . . . . . . . 78 4.6 Convergence curves of the IR-MSE based RR-CSI experiments. Cost function con- vergence curves show 2-3x faster of speed of the proposed LADMM over the previous ADMM. The acceleration is even more signicant in terms of dNRMSEs. . . . . . . 78 4.7 Convergence curves in the multi spin echo based T 2 -relaxation only experiments. About 2x faster convergence speed was observed for the proposed LADMM over the previous ADMM in terms of both cost function values and dNRMSEs. . . . . . . . . 79 4.8 Convergence curves in the ViSTa MRF based RR-CSI experiments. Faster conver- gence of the proposed LADMM over the previous ADMM were observed for both the cost function values and dNRMSEs. As noticed in the dNRMSEs, the conditioning of the inverse problem is so bad in this case that the dNRMSEs are still changing for both algorithms even after a long running time especially for the previous ADMM. 80 4.9 Component spectrum locations and spatial maps in the ViSTa MRF based RR-CSI experiments. Six anatomically plausible components were identied where compo- nents 1-6 seem to correspond to white matter structures, gray matter structures, CSF, venous, myelin water and scalp respectively. . . . . . . . . . . . . . . . . . . . . 81 4.10 Convergence comparison of the proposed LADMM w/ Eq. (4.26) and Eq. (4.28) in the IR-MSE based RR-CSI experiment. . . . . . . . . . . . . . . . . . . . . . . . . . 82 ix Abstract MRI is a powerful tool in both biomedical research and clinical diagnosis. However, the trade- os among data acquisition time, image signal-to-noise-ratio, and spatial resolution often hinders the applications of MRI in its full potentials. In addition to imaging hardware improvement and advanced pulse sequence development, constrained MRI methods have greatly expanded our tool- boxes in dealing with these trade-os. In constrained MRI, we utilize prior information about the characteristics of the underlying MRI images to perform data acquisition, image reconstruction, and image analysis tasks more eciently. This approach generally requires the use of mathemat- ical optimization techniques, although the optimization problems are often challenging to solve eciently due to their large-scale and non-trivial structure. In this dissertation, three novel contributions in mathematical optimization for constrained MRI were made. First, we proposed a novel method that utilizes phase constraints to accelerate MRI data acquisition based on non-Fourier radiofrequency encoding. While phase constraints are used classically in MRI, we believe that this is the rst time that phase constraints are being applied to enable acceleration along a non-Fourier encoded spatial dimension. We make the novel observation that phase constraints can indeed be successfully used to reduce the number of required non- Fourier encodings, although this requires careful design of the non-Fourier encoding scheme. Results are presented in the context of gSlider, an acquisition method designed for highly-ecient high- resolution diusion MRI. Second, we proposed a new algorithm for the separate regularization of magnitude and phase in MRI reconstruction problems. Our approach is based on a novel application of the proximal alternating linearized minimization algorithm (PALM), and incorporates additional novel features (i.e., Nesterov's momentum and independent selection of the step sizes for each coordinate) to increase convergence speed. Depending on the application, our proposed algorithm x can be hundreds of times faster than existing algorithms for this problem. Finally, we proposed another novel algorithm for spatial-spectral partial volume compartment mapping with applications to multicomponent diusion and relaxation MRI. Our proposed algorithm is based on a novel application of the linearized alternating directions method of multipliers (LADMM) approach that takes advantage of the special structure of the inverse problem, and depending on the dataset, can achieve up to 5-fold acceleration compared to previous algorithms for this problem. xi Chapter 1 Introduction As a non-invasive imaging method that can provide rich tissue contrast information, magnetic resonance imaging (MRI) [3] has been a unique and very important tool in both biomedical research and clinical diagnosis. However, the potentials of MRI in applications are often restricted by the trade-os among data acquisition time, image signal-to-noise ratio (SNR) and resolution [4, 5, 6]. In addition to imaging hardware improvement and advanced pulse sequence development, constrained MRI has been proposed to greatly improve our abilities in dealing with these trade-os by utilizing various prior knowledge about the characteristics of the underlying MR images in data acquisition, image reconstruction and analysis [7, 8]. For instance, smooth phase characteristic has been successfully applied to reduce up to half of the k-space data samples in the well-known partial Fourier technique [7]. In another popular example of compressed sensing MRI [9], the sparsity characteristics of MR images in either the ambient or transform domain are being utilized to both reduce k-space samples and enhance image reconstructions. In constrained MRI, model based methods play an increasingly important role where mathemat- ical optimization has been an essential tool [10]. This results in various challenging optimization problems that could be non-convex and/or large-scale in size. Therefore, how to model and eventu- ally solve constrained MRI problems with optimization methods and algorithms is highly non-trivial and is the focus of this dissertation. More specically, we have proposed three constrained MRI methods for more ecient data acquisition, faster image reconstruction and analysis respectively. 1 Chap.1 Introduction 2 1.1 Main Results For constrained data acquisition, we explored the idea of using smooth phase constraint to reduce samples in the non-Fourier RF-encoded experiments. Although smooth phase has been commonly used to reduce k-space samples, we believe this work is the rst to explore smooth phase in the non-Fourier encoding settings. We have shown from estimation theoretic point view that smooth phase can indeed facilitate data acquisition in non-Fourier encoding settings as long as we use properly designed encoding schemes presented in this work. We also proposed a smooth phase constrained decoding method for image reconstruction with the newly designed encoding schemes. Simulation results veried the proposed method. For constrained image reconstruction, we proposed a fast algorithm named PALMNUT for constrained image reconstruction with separate magnitude and phase regularization. The proposed algorithm is based on PALM and accelerated by incorporating the well-known Nes- terov's momentum technique and a novel uncoupled step size strategy. We show guaranteed cost function convergence results of the proposed algorithm without including Nesterov's mo- mentum. Numerical results show that the proposed algorithms converge signicantly faster than the state of the art algorithms. For constrained image analysis, we proposed a fast algorithm for spatial-spectral partial vol- ume compartment mapping. The proposed algorithm is based on linearized ADMM that can greatly reduce memory consumption and accelerate convergence speed. We tested the algo- rithm in various diusion and relaxation MRI scenarios and show that the proposed algorithm converges drastically faster than the state-of-the-art method. 1.2 Organization of the Document This document is organized as follows. Chapter 2 introduces a novel constrained data acquisition method for accelerating non-Fourier Radiofrequency-encoded data acquisition. Chapter 3 evaluates a new algorithm PALMNUT for constrained image reconstruction. Chapter 4 describes an ecient Chap.1 Introduction 3 algorithm for image analysis with spatial spectral partial volume compartment mapping. Chapter 5 concludes the dissertation. Chapter 2 Phase Constraints Enable Reduced RF Encodings with applications to gSlider-based Diusion MRI 2.1 Introduction 1 Magnetic resonance imaging (MRI) is a powerful and exible tool for biomedical imaging appli- cations, though some of its main limitations include limited signal-to-noise ratio (SNR) and long data acquisition times [49]. Over the decades since the invention of MRI, many dierent com- putational image reconstruction methods have been developed that seek to mitigate these issues by enabling high-quality MRI images to be formed from subsampled and/or low-quality data [50]. One of the earliest computational imaging approaches is based on the observation that, under certain smoothness assumptions about the phase of the MRI image, the Fourier data that is typi- cally acquired in MRI (also called k-space data) will possess certain symmetry characteristics that allow one side of k-space to be predicted from data measurements acquired on the opposite side [50, 51, 52, 53, 54, 7, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64]. This phase-constrained approach is 1 Some of the text and gures in this chapter have been previously published in [48], and are copyright of the ISMRM. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works must be obtained from the ISMRM. 4 Chap.4 Phase Constraints Enable Reduced RF Encodings with applications to gSlider-based Diusion MRI 5 powerful and widely-used because it can reduce traditional data-sampling requirements by nearly a factor of two. Moreover, phase-constraints are complementary to other forms of constraints that are commonly used for computational image reconstruction, allowing for their synergistic combination. The previous literature on phase-constrained MRI reconstruction has generally assumed a Fourier data acquisition model. One reason for this is that phase constraints were originally con- ceived from the symmetry properties of the Fourier transform, and are thus naturally formulated in the context of Fourier data. Another major reason is that the large majority of MRI experi- ments are based on Fourier encoding, and non-Fourier encoding methods have been less thoroughly explored. In this work, we consider the question of whether or not phase constraints might also be benecial in the context of non-Fourier encoded data, with a specic emphasis on non-Fourier encoding methods that perform spatial encoding along certain spatial dimensions using specially- designed spatially-selective radiofrequency (RF) excitation pulses. Non-Fourier RF encoding is well-established in MRI [65], and prominent examples include linescan imaging and Hadamard- encoded MRI [66, 67, 68, 69]; wavelet-encoded MRI [70]; SVD-encoded MRI [71]; and MRI with random encoding [72]. While phase-constrained reconstruction has been combined with non-Fourier RF encoding in previous work [62, 73], the previous work was focused on reducing the number of k-space measurements without changing the number of measured RF encodings. There has also been other work that used q-space regularization constraints to reduce the number of RF encodings needed for each diusion-weighted image in the context of RF-encoded diusion MRI [74, 75, 76], but these methods did not use phase constraints. To the best of our knowledge, the use of phase constraints to reduce the total number of measured RF encodings has never been previously ex- plored. Our personal interest in phase-constrained reconstruction of RF-encoded data is motivated by our recent work on the generalized slice dithered enhanced resolution (gSlider) approach to MRI data acquisition and image reconstruction [77, 78]. gSlider uses Hadamard-like non-Fourier RF encoding methods to gain major SNR eciency advantages and enable ultra-high resolution acqui- sitions in applications like diusion MRI. While this approach enables fast and highly-ecient data Chap.4 Phase Constraints Enable Reduced RF Encodings with applications to gSlider-based Diusion MRI 6 acquisition procedures, experiments could be even faster if the requirement for full RF-encoding could be relaxed by the introduction of phase constraints. This paper is organized as follows. In Sec. 2.2, we introduce the basic gSlider data acquisi- tion model, and perform an estimation-theoretic analysis based on the Cram er-Rao lower bound (CRLB) to evaluate the utility of phase-constraints in this context. Our theoretical analysis reveals that phase-constraints can indeed be used to reduce RF encoding requirements, although the utility of phase constraints will be highly dependent on how the RF encoding scheme is designed. As such, this section also investigates the optimal design of RF encodings. While the theoretical analysis from Sec. 2.2 provides strong theoretical justication for the use of phase constraints to accelerate the RF-encoding dimension in gSlider, our theoretical arguments were based on certain simplifying assumptions that are unlikely to hold in practical MRI applications. As a result, Sec. 2.3 inves- tigates the empirical performance of a phase-regularized reconstruction approach using simulated data with more realistic characteristics. The results of these simulations are consistent with our theoretical analysis, and conrm the potential to use phase constraints to reduce the number of RF encodings in gSlider. The theoretical and empirical results from these two sections (Secs. 2.2 and 2.3) demonstrate the theoretical promise and potential practical value of phase constraints for RF-encoded data, and comprise the bulk of the contribution of this work. Although the practical implementation and prospective evaluation of the proposed approach on real MRI hardware are beyond the scope of this paper, the next section (Sec. 4.5) includes a detailed discussion of next steps that would need to be taken on the path towards practical implementation. This section also includes the discussion of other topics, such as the combination of phase constraints with other types of regularization. Finally, our conclusions are presented in Sec. 4.6. A preliminary account of portions of this work was previouosly presented in a recent conference [48]. Chap.4 Phase Constraints Enable Reduced RF Encodings with applications to gSlider-based Diusion MRI 7 2.2 Combining gSlider with Phase Constraints 2.2.1 The Standard (Unconstrained) gSlider Model The gSlider approach relies on the clever combination of several dierent acquisition and recon- struction approaches, which can make it somewhat complicated to describe. Similar to previous work [78], we will present a simplied description that neglects some of the less-relevant complex- ities, and will allow us to focus narrowly on the aspects of gSlider that are related to our main contributions { readers interested in a more comprehensive description of gSlider are referred to Ref. [77]. We will assume that 3D imaging data is measured with conventional Fourier encoding for the rst two (in-plane) dimensions, with RF encoding applied along the third (slice) dimension. If parallel imaging reconstruction or simultaneous multislice reconstruction methods are needed to reconstruct the in-plane dimensions, we will assume that these methods have already been applied. We will also assume that if multichannel data was acquired, then a phase-preserving coil-combination method was applied to transform the multichannel data into single-channel data (e.g., linear voxel-by-voxel coil-combination using measured sensitivity maps [79]), such that it only remains to reconstruct the RF encoding dimension from complex-valued single-channel data. We assume that we are presented with measured data in the form ofQ dierent complex-valued single-channel RF-encoded 2D imagesd q (r n ), withq = 1;:::;Q, where each 2D image is comprised of N in-plane voxels at positions r n , with n = 1;:::;N. Each of these RF-encoded images is modeled as the linear combination of S dierent unknown subslices m s (r n ), with s = 1;:::;S, and our reconstruction task is to estimate m s (r n ) from the measured data d q (r n ). Specically, we assume a linear response model [65] such that d q (r n ) =e iq (rn) S X s=1 g qs m s (r n ) + q (r n ) (2.1) Chap.4 Phase Constraints Enable Reduced RF Encodings with applications to gSlider-based Diusion MRI 8 for q = 1;:::;Q and n = 1;:::;N, where q (r n ) represents the thermal noise which we assume to be i.i.d. circular complex white Gaussian 2 ; g qs is the complex-valued encoding coecient that the qth RF encoding applies to the sth subslice; and q (r n ) represents an image \background phase" term that is assumed to be shared by all of the S subslices, and represents the combination of all sources contributing to the MR image phase aside from the RF-encoding coecients g qs . Since we have used q (r n ) and the phase of g qs to absorb all sources of phase, we will model the subslice images m s (r n ) as purely real-valued. Note that, due to the random shot-to-shot phase variations associated with diusion MRI, we have assumed that the background phase q (r n ) will be dierent for dierent RF encodings, as in previous gSlider work [77, 78]. In this model, we assume that both the phase values q (r n ) and the subslice values m s (r n ) are unknown and need to be estimated, while g qs is assumed to be known perfectly. For simplicity, we will also sometimes express Eq. (2.1) in an equivalent matrix-vector form as d n = (Gm n )e i n + n (2.2) for n = 1;:::;N, where represents the Hadamard (elementwise) product of two vectors, expo- nentiation of a vector is used to denote elementwise exponentiation of each entry, and d n 2 C Q , m n 2 R S , n 2 R Q , n 2 C Q , and G2 C QS are respectively the vectors and matrices formed from the values of d q (r n ), m s (r n ), q (r n ), q (r n ), and g qs . Combining all N spatial locations together into a single equation can be achieved using d = ((I N G)m)e i +; (2.3) 2 Note that for real data, the actual noise distribution will be dependent on the parallel imaging, simultaneous multi-slice, and coil-combination methods that are applied. While linear methods will preserve the Gaussianity of the original noise distribution, it is likely that some of these reconstruction steps will introduce some degree of spatial noise correlation and spatial heteroskedasticity. As such, our noise model is simpler than real data would be. Nevertheless, we use this simplication because it is helpful for obtaining theoretical characterizations that are easy to interpret, and we believe that the positives outweight the negatives. We note that (1) it would have been straightforward for us to perform theoretical analysis with heteroskedastic noise modeling if so desired, although with substantially increased complexity for displaying and interpreting the results; and (2) we do not expect our main observations about the eciency of using phase constraints to be modied if we had used heteroskedastic correlated noise instead of i.i.d. white noise. Chap.4 Phase Constraints Enable Reduced RF Encodings with applications to gSlider-based Diusion MRI 9 where the vectors d2 C QN , m2 R SN , 2 R QN , and 2 C QN are respectively obtained by concatenating the voxelwise vectors d n , m n , n , and n into long vectors, I N 2 R NN denotes the NN identity matrix, and denotes the standard Kronecker product. Notice that there are S +Q unknown real-valued parameters at each spatial location (the values of m s (r n ) for s = 1;:::;S and q (r n ) for q = 1;:::;Q), while we have Q complex-valued measurements (which can be represented as 2Q real-valued measurements if split into real and imaginary components) in the form of d q (r n ) at each spatial location. As such, a simple degrees- of-freedom argument suggests that, in the absence of additional prior information, we must take more RF encodings than subslices (i.e., QS) in order to have hope in estimating the unknown parameters. The original gSlider implementation [77] is based on this principle, and takes Q =S = 5 with the encoding matrix given by G = 2 6 6 6 6 6 6 6 6 6 6 4 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 7 7 7 7 7 7 7 7 7 7 5 ; (2.4) which is designed so that it is easy to estimate q (r n ) from the measured data d q (r n ). Once the phase is estimated, the data d n can be phase corrected, and the desired images m n can be estimated by simply inverting the matrix G, which is a well-conditioned linear inverse problem [77]. In this work, we will be interested in whether the use of phase constraints will allow us to take Q > < > > : 0; q2V +1; q = 2V: (3.27) The new objective function from Eq. (3.26) now has four terms in it (i.e., J(;), R 1 (), R 2 (), andI V ()), and to apply PALM, it is necessary to associate these with the PALM terms H(;), F (), and G(). The function J(;) is always smooth and involves both m and q, so we will associate it with H(;). The functionI V () is always nonsmooth and only involves q, so we will Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 36 associate it with G(). While these associations are straightforward (we have no other options), we potentially have options for R 1 () and R 2 () depending on their smoothness characteristics. If R 1 () is nonsmooth, then it has to be associated with F (), but if it is smooth (Lipschitz) then we could either choose to associate it with F () or incorporate it into H(;). Similarly, if R 2 () is nonsmooth, then it has to be associated with G(), but if it is smooth (Lipschitz) then we could either choose to associate it with G() or incorporate it into H(;). Although we have dierent options, the remainder of this paper will assume (for simplicity and without loss of generality) that smooth regularization functions will always be incorporated into H(;). This choice leads to function associations that are summarized in Table 3.1. R 1 () is Smooth? R 2 () is Smooth? H(;) F () G() Y Y J(m; q) + R 1 (m) +R 2 (q) 0 I V (q) N Y J(m; q) +R 2 (q) R 1 (m) I V (q) Y N J(m; q) +R 1 (q) 0 R 2 (q) +I V (q) N N J(m; q) R 1 (m) R 2 (q) +I V (q) Table 3.1: Associations between PALM and Eq. (3.26) depending on whether R 1 () and R 2 () are smooth (Lipschitz). With these assignments, the PALM algorithm from Alg. 2 can be directly applied, although it is still necessary to specify the computation of c k , d k , prox F (;), and prox G (;). Although these computations will necessarily depend on the characteristics of R 1 () and R 2 (), we will provide concrete illustrations of these calculations for two typical choices of regularization penalties (we will use these same choices of regularization penalties in the validation study presented later in the paper). Of course, these two illustrations do not encompass every possibility, and interested readers are referred to Refs. [17, 114, 115] for further discussion and examples of computating proximal operators. However, it should be noted that the two illustrations below focus on the case where R 2 () is smooth (the rst two rows of Table 3.1), as we have found that it is frequently nontrivial to derive the prox G (;) operator when G() incorporates both R 2 () andI V () (the last two rows of Table 3.1). Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 37 3.3.1.1 Huber-function regularization of m with Tikhonov regularization of q For our rst illustration, we will consider the case where magnitude regularization takes the form of either R 1 (m) = 1 L X `=1 h ([Bm] ` ) (3.28) or R 1 (m) = 1 L X `=1 h 0 @ v u u t T X t=1 j[B t m] ` j 2 1 A ; (3.29) and phase regularization takes the form of R 2 (q) = 2 2 kCqk 2 2 : (3.30) In these expressions, 1 and 2 are positive scalar regularization parameters that can be respectively adjusted to control the in uence of the magnitude and phase regularization terms, [z] ` is used to denote the`th entry of the vector z, andh () is the Huber function (with parameter> 0) dened as h (t) = 8 > > < > > : 1 2 jtj 2 ; jtj jtj 1 2 ; jtj>: (3.31) The Huber function is a smooth, convex function that is commonly used for both robust statis- tics (to mitigate the eects of outliers) [118] and for edge/discontinuity/sparsity-preserving image regularization [119, 120]. As can be seen from Eq. (3.31), the Huber function is similar to the ` 1 - norm for large values of its argument, but unlike the ` 1 -norm, is also smooth at the origin because it behaves like a squared ` 2 -norm for small values of its argument. The Huber function with a small value of is frequently chosen as a dierentiable surrogate for sparsity-promoting ` 1 -norm regularization, while choosing larger values of can make the Huber function more tolerant to smoothly-varying image regions, more resilient to noise, and easier to characterize theoretically [119, 121, 8]. The regularization in Eq. (3.28) is based on applying the Huber function to a single transfor- mation B2C LN (e.g., a wavelet transform, a nite-dierence transform, etc.) of the magnitude Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 38 vector m. The more general regularization in Eq. (3.29) applies the Huber function to a combination ofT dierent transforms B t 2C LN of m, which can be useful for imposing additional transform- domain structure. For example, combining a horizontal nite-dierence transform with a vertical nite-dierence transform within Eq. (3.29) is a common way to achieve isotropic regularization [122]. In addition, Eq. (3.29) is related to concepts of simultaneous sparsity [123], and our previous work has used regularization penalties of this form to impose the constraint that multi-contrast images of the same anatomy will frequently have correlated edge characteristics [124, 8, 121, 111]. Since Eq. (3.28) is a special case of Eq. (3.29) with T = 1, our description below will assume use of the more general form of Eq. (3.29). The regularization in Eq. (3.30) corresponds to standard quadratic (Tikhonov) regularization. If C is chosen as a spatial nite-dierence operator, this type of regularization can be good at imposing the constraint that the image phase should be spatially smooth without major discontinuities [103, 105, 111, 106, 108]. While not every MRI image will have smooth phase characteristics, most do, and smooth phase is a common constraint within the image reconstruction literature [7, 62]. For this case, both R 1 () and R 2 () are smooth, corresponding to the situation in the rst row of Table 3.1. As such, to implement PALM, we use the assignments: H(;) 1 2 kA (m q) bk 2 2 + 1 L X `=1 h 0 @ v u u t T X t=1 j[B t m] ` j 2 1 A + 2 2 kCqk 2 2 F () 0 G() I V (q): (3.32) The gradients of H(;) needed for PALM are r m H(m; q) =< ( q A H A(m q) A H b + 1 T X t=1 B H t W(m)B t m ) (3.33) and r q H(m; q) = m A H A(m q) A H b + 2 C H Cq; (3.34) Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 39 where W(m) is an LL diagonal matrix depending on m with `th diagonal entry given by [W(m)] `` = 1= max 8 < : ; v u u t T X t=1 j[B t m] ` j 2 9 = ; : (3.35) These gradient expressions give rise to Lipschitz-type upper bounds in the form of Eqs. (3.9) and (3.10), such that the majorization and descent conditions will be satised whenever c k kAk 2 + 1 T X t=1 B H t B t (3.36) and d k kAk 2 k ^ m k k 2 1 + 2 kCk 2 ; (3.37) wherekk 1 denotes the ` 1 -norm. Please refer to [36, 125, 38] for insight on obtaining Lips- chitz constants for this kind of Huber function. The spectral norms in these expressions are not iteration-dependent, and can be precomputed and reused throughout the iterative process. If it is dicult to analytically calculate the spectral norm values, they can also be evaluated using standard computationally-ecient numerical methods like Lanczos iteration [126]. Finally, the proximal operators needed for PALM are given by prox F (w k ;c k ) = w k (3.38) and prox G (z k ;d k ) = z k jz k j ; (3.39) where division is performed elementwise and we choose the convention that 0 0 = 1. 3.3.1.2 ` 1 regularization of m with Huber-function regularization of q For our second illustration, we will consider the case where R 1 (m) = 1 kTmk 1 (3.40) Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 40 and R 2 (q) = 2 L X `=1 h 0 @ v u u t T X t=1 j[B t q] ` j 2 1 A : (3.41) The` 1 -norm penalty with sparsifying transform matrix T from Eq. (3.40) is standard for promoting transform-domain sparsity, and has been previously used to regularize the magnitude vector m in several applications [102, 104, 106, 107, 108, 109, 110]. The characteristics of the Huber function from Eq. (3.41) have been discussed previously. By taking a small value of the parameter , the Huber function can be used as a smooth approximation of the ` 1 -norm in order to enable sparsity- promoting and/or edge-preserving regularization of the phase image, which can useful for some applications with more complicated phase characteristics [109, 106]. As indicated earlier, choosing R 2 () to be a smooth function and including it within H(;) enables simplied computation of prox G (;). In this case, R 1 () is non-smooth while R 2 () is smooth, corresponding to the second row of Table 3.1. As such, to implement PALM, we use the assignments: H(;) 1 2 kA (m q) bk 2 2 + 2 L X `=1 h 0 @ v u u t T X t=1 j[B t q] ` j 2 1 A F () 1 kTmk 1 G() I V (q): (3.42) The gradients of H(;) needed for PALM are r m H(m; q) =< q A H A(m q) A H b (3.43) and r q H(m; q) = m A H A(m q) A H b + 2 T X t=1 B H t W(q)B t q: (3.44) Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 41 These gradient expressions give rise to Lipschitz-type upper bounds in the form of Eqs. (3.9) and (3.10), such that the majorization and descent conditions will be satised whenever c k kAk 2 (3.45) and d k kAk 2 k ^ m k k 2 1 + 2 T X t=1 B H t B t : (3.46) As before, the spectral norms in these expressions are not iteration-dependent, and can be precom- puted and reused throughout the iterative process. Finally, assuming that T is a unitary transform such that T H = T 1 , the proximal operator for F () is given by [115] prox F (w k ;c k ) = T H diag maxfjTw k j 1 =c k ; 0g jTw k j Tw k ; (3.47) where maximization, absolute value, and division operations are performed elementwise. Note that G() is the same as in the previous illustration, and therefore has the same proximal operator (Eq. (3.39)). 3.3.2 Uncoupled Step Sizes Although the PALM algorithm described in the previous section provides a novel and eective approach for solving Eq. (3.3), we observe in this section that it may be very conservative and computationally inecient for PALM to use the same value of d k (and therefore, the same step size 1=d k ) for all elements of the phase vector p. This ineciency stems from the fact that d k is set based on the global Lipschitz constant (eectively, the worst-case rate of change of the gradient along any possible direction), while we have observed that the rate of change of the gradient can be much smaller than the worst-case along specic directions. Concretely, using the global Lipschitz constant means that the step size will depend on the maximum value of m, while we observe that it can be much better for the step size for each coordinate to instead depend on the coordinatewise Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 42 values of m. This observation motivates us to investigate and utilize coordinatewise bounds on the rate of change of the gradient, enabling uncoupled coordinatewise step sizes. Some readers might view our use of uncoupled step sizes as a form of iteration-dependent preconditioning. For the sake of generality, we will rst describe this approach for the general setting of Sec- tion 3.2, where we are given a generic smooth real-valued function H(x; y). The PALM approach utilized majorants of H(x; y) that were derived based on the global (scalar-valued) Lipschitz con- stants L 1 (y) and L 2 (x) of Eqs. (3.9) and (3.10). In this section, we instead make the assumption that a vector L 1 (y)2R N 1 can be found such that, for a xed value of y, we have hr x H(x 1 ; y)r x H(x 2 ; y); x 1 x 2 i p L 1 (y) (x 1 x 2 ) 2 2 (3.48) for8x 1 ; x 2 2 R N 1 , where the square-root operation is applied elementwise. Similarly, we assume that a vector L 2 (x)2R N 2 can be found such that, for a xed value of x, we have hr y H(x; y 1 )r y H(x; y 2 ); y 1 y 2 i p L 2 (x) (y 1 y 2 ) 2 2 (3.49) for8y 1 ; y 2 2R N 2 . It should be noted that if the global Lipschitz continuity conditions of Eqs. (3.9) and (3.10) are known to be satised, then a vector L 1 (y) satisfying Eq. (3.48) can be trivially obtained by setting all of its entries equal to L 1 (y), with an analogous argument holding true for L 2 (x). In particular, the Cauchy-Schwarz inequality combined with Eq. (3.9) implies that hr x H(x 1 ; y)r x H(x 2 ; y); x 1 x 2 ikr x H(x 1 ; y)r x H(x 2 ; y)k 2 kx 1 x 2 k 2 L 1 (y)kx 1 x 2 k 2 2 = p L 1 (y)(x 1 x 2 ) 2 2 : (3.50) However, Eqs. (3.48) and (3.49) are more exible than Eqs. (3.9) and (3.10) because a dierent Lipschitz-like constant can be used for every coordinate, and many of these entries can be potentially much smaller than the global Lipschitz constant (because the function gradient may change much more slowly along these directions). Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 43 From an optimization perspective, Eqs. (3.48) and (3.49) are important because they enable the use of potentially better majorants than were used by PALM, as described by the following theorem. Theorem 3.1. Consider the setting described in Section 3.2, and assume that the smooth real- valued function H(x; y) satises the conditions of Eq. (3.48). Then given a vector c k 2 R N 1 and assuming y is held xed at y = ^ y k1 , the function H(^ x k1 ; ^ y k1 ) +hx ^ x k1 ;r x H(^ x k1 ; ^ y k1 )i + 1 2 k p c k (x ^ x k1 )k 2 2 +F (x) +G(^ y k1 ) (3.51) is a majorant of (x; ^ y k1 ) at the point x = ^ x k1 whenever c k L 1 (^ y k1 ) (elementwise). Simi- larly, given a vector d k 2R N 2 , assumingH(x; y) satises the conditions of Eq. (3.49), and assuming x is held xed at x = ^ x k , the function H(^ x k ; ^ y k1 ) +hy ^ y k1 ;r y H(^ x k ; ^ y k1 )i + 1 2 p d k (y ^ y k1 ) 2 2 +F (^ x k ) +G(y) (3.52) is a majorant of (^ x k ; y) at the point y = ^ y k1 whenever d k L 2 (^ x k ) (elementwise). The proof of this theorem is given in Appendix ??. Although this theorem is stated for real- valued vectors (for consistency with previous descriptions), the same also holds true for complex- valued vectors x and y. Following the same approach as PALM but replacing the original PALM majorants (Eqs. (3.7) and (3.8)) with the new majorants from Thm. 3.1 results in the new PALM algorithm with uncou- pled step sizes given in Alg. 3. This algorithm uses proximal operators with vector-valued c k and d k , which we dene as prox F (w k ; c k ), argmin x2R N 1 F (x) + 1 2 k p c k (x w k )k 2 2 (3.53) Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 44 and prox G (z k ; d k ), argmin y2R N 2 G(y) + 1 2 p d k (y z k ) 2 2 : (3.54) Algorithm 3 PALM for Eq. (3.5) with Uncoupled Step Sizes Initialization: Set k = 1 and initialize ^ x 0 and ^ y 0 . while not converge do Choose c k L 1 (^ y k1 ) (elementwise) w k = ^ x k1 diag(c k ) 1 r x H(^ x k1 ; ^ y k1 ) ^ x k = prox F (w k ; c k ) Choose d k L 2 (^ x k ) (elementwise) z k = ^ y k1 diag(d k ) 1 r y H(^ x k ; ^ y k1 ) ^ y k = prox G (z k ; d k ) k k + 1 end while Just like for the original PALM algorithm, the use of valid majorants in our new algorithm guarantees that it will monotonically decrease the objective function value, and that the objective function value converge if (x; y) is bounded from below. Further, this new algorithm reduces to the original PALM algorithm if the c k and d k are treated as scalars (i.e., the values in dierent entries are always the same) instead of choosing dierent values in each of the entries. We hypothesize that making additional assumptions about the structure of the cost function (e.g., the Kurdyka- Lojasiewicz property that has been used with the PALM algorithm [19]) may enable a proof that this new algorithm has guaranteed convergence to a critical point, although such a proof is beyond the scope of the present work. Now that the general approach has been described, let's consider the application to the specic magnitude and phase optimization problem of interest from Eq. (3.3). Without further assumptions about the structure of A, we do not observe any special coordinatewise structure for the magnitude subproblem. As such, we can simply utilize the global Lipschitz constant in this case, setting c k =c k 1, where 1 is a vector whose entries are all equal to 1 and c k is the value obtained for the PALM algorithm as described in the previous section. Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 45 However, for the phase subproblem with a xed value of m, we observe that hr q J(m; q 1 )r q J(m; q 2 ); q 1 q 2 ikkAkjmj (q 1 q 2 )k 2 2 ; (3.55) where the magnitude operation is applied elementwise to the vector m. Thus, for the rst illustra- tion from the previous subsection (Huber-function regularization of m with Tikhonov regularization of p), our new algorithm can adopt d k kAk 2 j ^ m k j 2 + 2 kCk 2 (3.56) instead of the previous expression from Eq. (3.37). Similarly, for the second illustration from the previous subsection (` 1 -regularization of m with Huber-function regularization of p), our new algorithm can adopt d k kAk 2 j ^ m k j 2 + 2 T X t=1 B H t B t (3.57) instead of the previous expression from Eq. (3.46). Note also that for both of the previous illus- trations, the prox G (;) expressions had no dependence on the actual value of d k , which allows us to simply reuse the same proximal operators for the new algorithm without modication. As such, for these illustrations, the main dierence between the PALM algorithm and our new algorithm is that PALM takes a uniform step size for the phase update that depends on the maximium value of m, while the new algorithm can take larger step sizes for coordinates where the corresponding value of m is small. 3.3.3 Nesterov's Momentum Acceleration and PALMNUT Our proposed PALMNUT (PALM with Nesterov's momentum and Uncoupled sTep sizes) algorithm is obtained by combining Alg. 3 with Nesterov's momentum technique. The basic idea of Nesterov's technique is that, instead of computing the next iterate ^ x k based on values derived from ^ x k1 , it can be better to instead nd the next iterate using values derived from a combination of the previ- ous iterates [43, 44], which can be interpreted as using \momentum" from previous iterations. For Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 46 convex optimization problems, this approach can even result in convergence rates that have optimal order [43, 44]. Of course, the problem of interest in this work is not convex, although it has been shown empirically (oftentimes without rigorous theoretical justication) that Nesterov's momen- tum technique can often substantially accelerate the iterative solution of nonconvex optimization problems. The idea of applying momentum to accelerate the convergence of PALM has been studied in Ref. [127], where the resulting algorithm was called inertial PALM (iPALM). Theoretical conver- gence results for iPALM were proven with restrictive parameter choices (dierent from Nesterov's original parameter choices), although it was also shown empirically that using Nesterov's original parameter choices generally led to much faster convergence, despite the lack of theoretical guaran- tees. Our empirical experience is also consistent with these previous observations, so our proposed PALMNUT algorithm similarly utilizes Nesterov's original parameter choices (following the con- cise form described by Ref. [45]). The nal PALMNUT algorithm, incorporating both uncoupled coordinatewise step sizes and Nesterov's momentum, is given in Alg. 4. Algorithm 4 PALMNUT Initialization: Set k = 1 and initialize ^ x 0 and ^ y 0 . Set u 0 = ^ x 0 and v 0 = ^ y 0 . while not converge do Choose c k L 1 (v k1 ) (elementwise) w k = u k1 diag(c k ) 1 r x H(u k1 ; v k1 ) ^ x k = prox F (w k ; c k ) u k = ^ x k + k1 k+2 (^ x k ^ x k1 ) Choose d k L 2 (u k ) (elementwise) z k = v k1 diag(d k ) 1 r y H(u k ; v k1 ) ^ y k = prox G (z k ; d k ) v k = ^ y k + k1 k+2 (^ y k ^ y k1 ) k k + 1 end while Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 47 3.4 Numerical Experiments In the following subsections, we evaluate PALMNUT in three dierent simulations that are represen- tative of a diverse set of real problems in MRI: sparsity-promoting reconstruction of undersampled k-space data [106, 108, 109], regularization-based denoising of a complex image [111, 121, 124], and using phase correction to enable the combination of multiple images acquired in the presence of experimental phase instabilities [111, 128, 129, 130, 131, 132]. In each of these cases, we compare PALMNUT against AM combined with NCG [103, 105, 106, 108, 111] and AM combined with phase cycling [109]. PALMNUT and AM combined with NCG are directly comparable, because they can both be used to regularize the exponentiated phase e ip , and therefore can both be applied to the exact same optimization problem. As a result, for both of these algorithms, the phase was regularized with R 2 (e ip ) as described previously, and we used the cost function value and the computation time to judge algorithm performance. However, the phase cycling heuristic [109] is intended to be used for the direct regularization of p. For this algorithm, we therefore had to instead use phase regularization of the form R 2 (p). In addition, by its nature, the phase cycling heuristic employs a dierent phase-cycled cost function at each iteration, which makes it dicult to plot a meaningful cost function value. To allow for comparisons with this dierent approach, we therefore also computed a normalized root-mean-squared error (NRMSE) metric for each algorithm. If f represents the ground truth vector and ^ f represents an estimate, the NRMSE of the estimate is given by NRMSE,k ^ f fk 2 =kfk 2 : (3.58) In order to focus on consequential errors, the NRMSE values were always computed after masking out the empty (background) parts of the eld-of-view. This was particularly benecial for AM with phase cycling, which frequently showed higher error levels in the image background compared to the other two approaches. Regularization parameters for each cost function in each scenario were empirically optimized to achieve the smallest possible nal NRMSE values for the two AM algorithms. Our implementation Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 48 (a) Magnitude (b) Phase (c) Sampling Mask Figure 3.1: The ground truth (a) magnitude and (b) phase images used for the undersampled MRI reconstruction scenario, along with (c) the 8-accelerated k-space sampling mask used to simulate the undersampled acquisition. of AM with NCG used the Polak-Ribiere version of NCG [133]. Our implementation of AM with phase cycling was based on code provided by the authors of Ref. [109] (available from https: //mrirecon.github.io/bart/). 3.4.1 Undersampled MRI Reconstruction In the rst set of evaluations, we considered the reconstruction of an MR image from 8-undersampled k-space data. The gold-standard magnitude and phase images, which were obtained from a real fully-sampled in vivo T1-weighted MRI acquisition with 256 256 in-plane matrix size, are shown in Fig. 3.1. This gure also shows the k-space sampling mask (corresponding to 8 undersampling) that we used to simulate an accelerated acquisition. For this simulation, the original fully-sampled k-space data (originally measured with 32 chan- nels) was coil-compressed down to 8 virtual channels to reduce computational complexity, and was then retrospectively undersampled using the aforementioned k-space sampling mask. For recon- struction, the A matrix was chosen according to the standard SENSE model [134], with sensitivity maps estimated using ESPIRiT [135]. The magnitude regularization took the form of an` 1 penalty as given by Eq. (3.40), where, following Ref. [109], the sparsifying transform T was chosen to be the unitary Daubechies-4 wavelet transform. The phase regularization took the form of a Huber- function penalty as given by Eq. (3.41), where the Huber parameter was chosen to be a small Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 49 0 500 1000 1500 2000 2500 3000 3500 4000 Computation time (sec.) 0.5 0.55 0.6 0.65 0.7 0.75 log 10 (Cost function) AM + NCG PALMNUT (a) 0 500 1000 1500 2000 2500 3000 3500 4000 Computation time (sec.) 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 NRMSE AM + NCG PALMNUT AM + Phase Cycling (b) (c) Zero-lled (d) PALMNUT Figure 3.2: Convergence plots and reconstructed images for the undersampled MRI reconstruction scenario. The convergence plots show (a) the cost function value as a function of computation time and (b) the NRMSE value for the complex image as a function of computation time. Also shown are the magnitude and phase images corresponding to (c) zero-lled reconstruction and (d) PALMNUT. Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 50 (a) Ground Truth (b) Noisy Figure 3.3: The (a) ground truth and (b) noisy magnitude and phase images used for the regularization-based MRI denoising scenario. number (i.e., = 0:001) in order to approximate the ` 1 -norm. Following Ref. [109], the transform we used for phase regularization was also a unitary Daubechies-4 wavelet transform. All three algorithms were initialized by applying SENSE-based coil-combination to the multi-channel images obtained by zero-lling the unmeasured data. Convergence results for all three algorithms are shown in Fig. 3.2, along with representative image reconstruction results. As can be seen, the cost function and NRMSE values converge to similar levels for both PALMNUT and AM with NCG, although PALMNUT converged much faster. Although AM with phase cycling converged to a result with a reasonably-good NRMSE, it was substantially worse than PALMNUT in both convergence speed and the nal achieved NRMSE value. 3.4.2 Regularization-based MRI Denoising In the second set of evaluations, we considered regularization-based denoising of a 230180 single- channel T1-weighted MR image, obtained by applying complex coil-combination to an 8-channel dataset and subsequently adding simulated complex Gaussian noise. The ground truth and noisy images are shown in Fig. 3.3. For reconstruction, the A matrix was an identity matrix. Following Refs. [121, 111, 8, 124], the magnitude was regularized using a Huber-function penalty as given by Eq. (3.28), where a Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 51 0 20 40 60 80 100 120 Computation time (sec.) 2 2.5 3 3.5 4 4.5 log 10 (Cost function) AM + NCG PALMNUT (a) 0 20 40 60 80 100 120 Computation time (sec.) 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 NRMSE AM + NCG PALMNUT AM + Phase Cycling (b) (c) PALMNUT Figure 3.4: Convergence plots and reconstructed images for the regularization-based MRI denoising scenario. The convergence plots show (a) the cost function value as a function of computation time and (b) the NRMSE value for the complex image as a function of computation time. Also shown are the (c) magnitude and phase images corresponding to PALMNUT. nite dierence transformation was used to enforce spatial smoothness of the image. Following Ref. [111, 103, 105, 106, 108], the phase was regularized using a Tikhonov penalty as given by Eq. (3.29), also using a nite dierence transformation to enforce spatial smoothness. All algorithms were initialized with the noisy image. Convergence results are shown in Fig. 3.4, along with representative reconstruction results. As can be seen, the results in this case are consistent with the previous case: PALMNUT and AM with NCG converged to similar NRMSE values, although PALMNUT was substantially faster, while AM with phase cycling was reasonably successful yet still substantially worse than the others in both speed and NRMSE. Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 52 (a) Magnitude (b) Noisy (c) Phases Figure 3.5: (a) The ground truth magnitude and (b) a representative noisy magnitude image, along with (c) the four dierent ground truth phase images used for the scenario with phase-corrected combination of multiple images. 3.4.3 Phase-Corrected MR Image Combination In the third experiment, we simulated a scenario that is common in diusion MRI, in which multiple measurements are made of the same image to enable averaging to improve SNR, but the phase varies randomly with each measurement due to experimental instabilities. This case was simulated based on actual diusion MRI magnitude and phase data from Ref. [111]. We simulated a case with four repetitions, based on one ground truth magnitude image and four ground truth phase images, all with matrix size 180 332. The magnitude image was combined with each of the dierent phase images to yield four dierent complex images, and then complex Gaussian noise was added to each result. The ground truth images and a representative noisy image are shown in Fig. 3.5. For reconstruction, we estimated a single shared magnitude image m and four dierent phase images p j corresponding to the four dierent noisy measured images b j ,j = 1; 2; 3; 4, by minimizing the cost function R 1 (m) + 4 X j=1 1 2 kme ip j b j k 2 2 +R 2 (e ip j ); (3.59) where the magnitude and phase regularization penaltiesR 1 () andR 2 () were chosen in exactly the same way as for the denoising scenario from the previous subsection. All optimization algorithms Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 53 0 20 40 60 80 100 120 140 Computation time (sec.) 2.5 3 3.5 4 4.5 log 10 (Cost function) AM + NCG PALMNUT (a) 0 20 40 60 80 100 120 140 Computation time (sec.) 0.06 0.08 0.1 0.12 0.14 0.16 0.18 NRMSE AM + NCG PALMNUT AM + Phase Cycling (b) (c) PALMNUT Figure 3.6: Convergence plots and reconstructed images for the phase-corrected MRI image combi- nation scenario. The convergence plots show (a) the cost function value as a function of computation time and (b) the NRMSE value as a function of computation time. Also shown is (c) the combined magnitude image obtained from PALMNUT. were initialized with the noisy phase images and by taking the average of the noisy magnitude images. The convergence plots shown in Fig. 3.6 show similar characteristics to those observed in the previous scenarios, with PALMNUT having a distinct advantage over the two alternative algo- rithms. 3.5 Discussion The results of the previous section demonstrated that PALMNUT can have major advantages relative to AM methods. However, PALMNUT represents a combination of three dierent ideas (PALM, Nesterov's momentum, and uncoupled stepsizes), and the previous experiments did not investigate which parts of PALMNUT contribute most to its performance. In order to gain more Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 54 insight, we did another set of experiments in which we compared PALMNUT against the origi- nal PALM algorithm (Alg. 2), PALM with Nesterov's momentum but without uncoupled stepsizes (iPALM) 3 , and PALM with uncoupled stepsizes but without Nesterov's momentum (Alg. 3). Re- sults are shown in Fig. 3.7 for both the undersampled MRI reconstruction (Section 3.4.1) and the regularization-based MRI denoising (Section 3.4.2) scenarios described previously. From these plots, we observe that both uncoupled step sizes and Nesterov's momentum independently improve the convergence speed relative to the original PALM method, with Nesterov's momentum contributing a little more than the use of uncoupled step sizes. However, PALMNUT's combination of all of these elements leads to the best overall results (i.e., reaching small cost function values sooner than the other algorithms and achieving the fastest empirical convergence of the cost function value). Interestingly, it was also observed that the NRMSE curves for PALMNUT and iPALM followed similar trends in this set of experiments, and that the NRMSE values appeared to converge sooner than the cost function values did. This is likely due to the well-known fact that NRMSE can be insensitive to some kinds of image features [136], although it should be clear from the cost function behavior that the image estimates are still changing even after the NRMSE curves appear to have converged. In addition, it is important to keep in mind that NRMSE values cannot generally be computed in practical applications, and unlike the cost function behavior, the NRMSE behavior could not be easily used as an algorithm stopping criterion. An interesting phenomenon we observed with PALMNUT is that the magnitude estimate fre- quently converged faster than the phase estimate (results not shown). As a result, it could be potentially benecial to update the phase estimate more frequently than the magnitude estimate, although we believe that such an exploration is beyond the scope of this paper. One of the key ingredients of PALMNUT is the use of coordinatewise step sizes based on coordinatewise Lipschitz-like bounds in the form of Eqs. (3.48) and (3.49). Although these bounds were formulated in the context of PALM, we believe that Lemma 1 (found in the Appendix and used in the proof of Theorem 1) represents a novel coordinatewise majorization relationship, which could also potentially be useful to construct majorants in more general optimization scenarios. 3 iPALM [127] has two momentum step sizes k and k . In our implementation, we set these to be equal with k = k = (k 1)=(k + 2) so that the only dierence between iPALM and PALMNUT is the uncoupled step sizes. Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 55 0 500 1000 1500 2000 2500 3000 Computation time (sec.) 0.5 0.55 0.6 0.65 0.7 0.75 0.8 log 10 (Cost function) PALM iPALM (Nesterov) PALM + Uncoupled step sizes PALMNUT 0 500 1000 1500 2000 2500 3000 Computation time (sec.) 0.07 0.08 0.09 0.1 0.11 0.12 0.13 NRMSE PALM iPALM (Nesterov) PALM + Uncoupled step sizes PALMNUT (a) Undersampled MRI Reconstruction 0 5 10 15 20 25 30 35 40 Computation time (sec.) 2 2.5 3 3.5 4 4.5 log 10 (Cost function) PALM iPALM (Nesterov) PALM + Uncoupled step sizes PALMNUT 0 5 10 15 20 25 30 35 40 Computation time (sec.) 0.11 0.12 0.13 0.14 0.15 0.16 NRMSE PALM iPALM (Nesterov) PALM + Uncoupled step sizes PALMNUT (b) Regularization-based MRI Denoising Figure 3.7: The convergence plots show the cost function value and the NRMSE value as a function of computation time for (a) the undersampled MRI reconstruction scenario and (b) the regularization-based MRI denoising scenario. Finally, we should mention that while PALM combined with uncoupled step sizes has guaranteed monotonic convergence properties, the convergence of the PALMNUT approach (which additionally includes Nesterov's momentum) has not yet been theoretically proven. This represents another potentially interesting topic for future research. Chap.3 PALMNUT: An Enhanced Proximal Alternating Linearized Minimization Algorithm with Application to Separate Regularization of Magnitude and Phase 56 3.6 Conclusion We proposed and evaluated a new algorithm called PALMNUT, which combines the PALM al- gorithm with Nesterov's momentum with uncoupled coordinatewise step sizes derived from coor- dinatewise Lipschitz-like bounds. Although our approach is general and can be applied to other computational imaging scenarios, our evaluation studies focused on MRI scenarios involving sep- arate regularization of the image magnitude and phase. Applying algorithms like PALM and PALMNUT to this problem also required us to reformulate this optimization problem in a novel way that is more compatible with them than the standard formulation. Our empirical results demonstrated that PALMNUT consistently had substantial advantages over previous approaches based on alternating minimization across several dierent MRI scenarios. As a result, we expect that PALMNUT will be useful for these kinds of MRI scenarios, and may also prove useful for more general computational imaging problems with similar optimization structure. Chapter 4 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 4.1 Introduction 1 Partial volume eects occur in virtually all medical imaging experiments, and result from the fact that biological tissues generally contain high-resolution features that are substantially smaller than the practically-achievable spatial resolution. In this context, the signal observed from a macroscopic image voxel will be a mixture of contributions from dierent microscopic sub-voxel tissue compartments. Partial volume mixing is clearly problematic in applications where it is 1 Some of the text and gures in this chapter have been previously published in [137], and are copyright of the ISMRM. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to reuse any copyrighted component of this work in other works must be obtained from the ISMRM. 57 Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 58 important to gain insights into the original (unmixed) high-resolution microscopic features of the tissue. This paper is focused on the problem of unmixing the contributions of dierent microscopic sub-voxel tissue compartments from images acquired with macroscopic resolution. Methods that can successfully solve this problem are highly signicant, because they can provide information about microscopic tissue features that are too small to be directly observed using conventional imaging methods, thereby transcending traditional resolution limits. To solve the partial volume unmixing problem robustly, it is important to acquire imaging data that is sensitive to the dierences between dierent sub-voxel compartments. One common way to achieve this is by acquiring a series of multiple images with various experimental contrast settings, under the assumption that the contrast variations will aect the sub-voxel compartments in distinct ways. For a nite number of components and assuming a linear mixture model, this leads to the data acquisition model m p (r n ) = Q X q=1 b( p ; q )f q (r n ) + p (r n ) (4.1) for p = 1;:::;P and n = 1;:::;N, where P is the total number of measured images in the image sequence; the vector r refers to the spatial location within the image and it is assumed that we observe N voxels at spatial positions r n , n = 1;:::;N; m p (r) is the pth measured image; p (r) is the noise within the pth measured image; p are the contrast-related experimental parameters for the pth measured image; Q is the number of compartments; q are the contrast-related tissue parameters for the qth compartment; f q (r) is a spatial map representing the contribution of the qth sub-voxel tissue compartment to the total signal; and b(; ) is a model of the ideal data corresponding to experimental parameters with tissue parameters . Given a model of this form, the contrast parameters q and spatial maps f q (r) for each sub-voxel compartment can be estimated from the measured data using an appropriate estimation method. Our personal interest in this type of multicomponent mixture modeling stems primarily from our own recent work on multiparametric correlation-spectroscopic imaging in MRI [1, 2] (involving multidimensional diusion, T 1 -relaxation, and/or T 2 -relaxation encoding { see also [138, 139, 140, 141, 142, 143] for related MRI work from other groups). However, multicomponent mixture models Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 59 are also widely used in a range of other experimental paradigms including dynamic PET [144], dynamic contrast-enhanced MRI and CT [145], MRI relaxometry [146], and diusion MRI [147], to name just a few. In many applications, it is common for the compartmental signal model b(; ) to be some form of exponential decay, which emerges as a consequence of the imaging physics. As a result, the unmixing modeling problem reduces to estimating the parameters of a multiexponential decay model. This kind of multiexponential estimation problem (sometimes informally described as a form of \inverse Laplace transform") appears frequently within the physical sciences, and has been studied for hundreds of years by the mathematics and physics communities [148]. Unfortunately, it has also been understood for hundreds of years that multiexponential estimation is a highly ill-posed inverse problem, meaning that there can be many models with very dierent parameters that can all t the measured data similarly well, and also that the solutions to the inverse problem can be highly sensitive to small noise perturbations. Because of this fundamental ill-posedness, it has proven necessary to use specialized methods in order to obtain robust unmixing results. From a data acquisition perspective, it has been demonstrated that multidimensional multiparametric acquisitions (e.g., encodingT 1 -T 2 or diusion- T 2 parameters simultaneously and nonseparably in an MRI context) can lead to inverse problems that are substantially better posed than lower-dimensional single-parameter acquisitions (e.g., that encode T 1 , T 2 , or diusion parameters alone) [2, 149]. From an estimation perspective, it has proven useful to introduce various types of constraints on the solution to the inverse problem. While many dierent types of constraints have been introduced over the decades [150, 151, 152, 153, 154, 1, 138, 140, 139, 141, 143, 2, 142, 155, 156, 157, 158, 159, 160], the work we present in this paper will focus exclusively on methods that assume spatial regularity of the compartmental spatial maps f q (r) [156, 157, 155, 158, 1, 159, 160, 2]. It has been shown theoretically that the use of spatial constraints can dramatically reduce the ill-posedness of the inverse problem in both one-dimensional [155] and multidimensional acquisition [2] scenarios compared to the traditional approach in which mixture modeling is applied independently to each voxel. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 60 While both multidimensional multiparametric acquisitions and spatially-regularized compart- ment estimation can greatly improve the quality of partial volume unmixing, these approaches are both associated with estimation methods that require increased memory and computational com- plexity. Multidimensional multiparametric acquisition naturally leads to increased computational complexity because of the curse of dimensionality. Meanwhile, spatial regularization methods re- quired increased computational complexity because it becomes necessary to estimate the model for multiple voxels simultaneously, instead of solving for each voxel independently. This increased computational complexity can be quite burdensome, and for example, we are aware of only a few examples of methods that estimate all voxels of large-scale spatial maps f q (r),q = 1;:::;Q, at the same time [155, 1, 2, 160], with other spatially-constrained methods opting instead to either only estimate small local image patches or using other heuristics to avoid the computational complexity of the full spatially-regularized inverse problem [156, 157, 158, 159]. The recent methods for solving large-scale spatially-regularized partial volume compartment mapping problems [1, 2, 160] have all made use of algorithms based on the alternating direction method of multipliers (ADMM) [161]. Although there are some dierences in the formulation and algorithmic implementation, all of these ADMM-based approaches end up introducing at least three additional sets of auxiliary variables, which, through clever use of variable-splitting, results in simplied subproblems that can each be solved eciently. While these algorithms produce good results, the ill-posedness of the inverse problem can still mean that they are quite slow to converge. In this work, we introduce a new fast algorithm for solving these kinds of spatially-regularized inverse problems. Our new algorithm is based on two insights: First, that we can save substantial amounts of memory by reducing the number of variable-splittings we introduce; and second, that we can also potentially improve convergence speed by reducing the number of subproblems that need to be solved compared to previous ADMM methods. These insights lead us to propose an algorithm, based on the linearized ADMM (LADMM) [162, 163, 164] (also known as proximal ADMM [165, 166] and generalized ADMM [167]) that only involves a single variable splitting, and results in a smaller number of subproblems that capture more of the structure from the original inverse problem compared to the previous ADMM-based methods, enabling faster convergence with Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 61 much less memory usage. A preliminary account of portions of this work was previously presented in a recent conference [137]. This paper is organized as follows. In Sec. 4.2, we provide a detailed description of the opti- mization problem that we are interested in solving, and review the existing algorithmic approaches for solving it. Then in Sec. 4.3, we review the LADMM approach and describe how we adapt it to the optimization problem at hand. In Sec. 4.4, we compare our new algorithm against existing algorithms in a range of MRI applications, including multicomponent diusion-T 2 relaxation cor- relation spectroscopic imaging (DR-CSI) [1], multicomponent T 1 -T 2 relaxation correlation spectro- scopic imaging (RR-CSI) with both inversion recovery multi spin echo (IR-MSE) [2] and magnetic resonance ngerprinting (MRF) [168] acquisitions, in addition to multicomponent T 2 relaxometry [151, 150, 146]. We conclude the paper with additional discussion and conclusions in Secs. 4.5 and 4.6, respectively. 4.2 Problem Formulation and Existing Methods 4.2.1 Formal Description of Optimization Problem While there are many dierent forms of multicomponent mixture models that can result from Eq. (4.1), our attention in this work is restricted to the class of \spectral" or \continuum" models [152, 142, 2, 169]. In this setting, rather than assuming that the number of compartments Q is a small nite number corresponding to a discrete set of tissue parameters q as in Eq. (4.1), it is instead assumed that the measured data results from a continuous mixture of contributions arising from a potentially-innite set of tissue parameter values. This leads to the integral equation m p (r n ) = Z b( p ; )f( ; r n )d + p (r n ); (4.2) for p = 1;:::;P and n = 1;:::;N, where f( ; r) is a high-dimensional \spectroscopic image" that includes both spatial dimensions r and spectral dimensions . For practical implementation, it Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 62 is common to discretize the spectroscopic image at a nite number of Q spectral positions q , q = 1;:::;Q, which results in m p (r n ) = Q X q=1 w q b( p ; q )f( q ; r n ) + p (r n ); (4.3) for p = 1;:::;P and n = 1;:::;N, where w q are the density normalization coecients (i.e., the numerical quadrature weights) required for accurate approximation of the continuous integral using a nite discrete sum. Assuming that we have measured real-valued magnitude images as is often the case in MRI, Eq. (4.3) can be conveniently expressed in matrix-vector form as m n = Kf n + n n ; (4.4) for n = 1;:::;N, where m n 2R P is the vector of P dierent measured image voxel values m p (r n ) from the nth voxel; n n 2R P is the vector of noise values p (r n ) from the nth voxel; K2R PQ is the \dictionary" describing the ideal imaging physics, with entries w p b( p ; q ); and f n 2R Q is the vector of the Q spectral values f( q ; r n ) from the nth spatial location. To simplify notation, we will also sometimes represent this expression as m = (I N K)f + n; (4.5) where m2R PN is obtained by concatenating the vectors m 1 ;:::; m N into a long vector, and with similar concatenation operations applied to form the spectroscopic image vector f2R QN and the noise vector n2R PN . In this expression I N represents theNN identity matrix, and represents the standard Kronecker product and results in the block-diagonal matrix I N K = 2 6 6 6 6 4 K 0 : : : : : : : : : 0 K 3 7 7 7 7 5 : (4.6) Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 63 4.2.2 Classical Voxel-By-Voxel Solutions Classically, the multicomponent mixture model is estimated by solving a voxel-by-voxel regularized non-negative least squares (NNLS) estimation problem of the form ^ f n = arg min fn2R Q fn0 1 2 km n Kf n k 2 2 +R(f n ); (4.7) for n = 1;:::;N. 2 In Eq. (4.7), the constraint f n 0 should be interpreted elementwise and corresponds to a physically-motivated nonnegativity constraint on the spectroscopic image f( ; r), while R() is an optional voxelwise regularization penalty and is a regularization parameter. Note that ifR() is chosen to be a convex function, then Eq. (4.7) is a convex optimization problem and can be globally optimized using a variety of dierent convex optimization methods [171]. Interestingly, NNLS solutions often naturally produce very sparse solutions even without additional regularization [172]. One of the classic and most commonly-used algorithms to solve NNLS problems in this context is the active set algorithm by Lawson and Hanson[173]. Theoretically, this algorithm converges to a globally optimal solution within a nite number of iterations, where the number of iterations generally scales with the number of unknownsQ [173]. For the voxel-by-voxel estimation approach, this NNLS algorithm can be very fast and ecient, although the lack of spatial regularity constraints can lead to results that are much worse than would have been obtained with spatial regularization. 4.2.3 Spatially-Regularized Solutions In contrast to the voxel-by-voxel approach from the previous subsection, the spatially-regularized approach estimates all voxels simultaneously in a coupled fashion, by solving a problem of the form ^ f = arg min f2R QN f0 1 2 km (I N K)fk 2 2 +R(f); (4.8) 2 Note that the use of least-squares in this expression implicitly corresponds to Gaussian noise-modeling assump- tions, though this same least-squares formulation can be easily adapted to Rician or noncentral chi noise statistics using half-quadratic majorize-minimize techniques [170]. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 64 where the regularization penalty R() is intentionally designed to encourage similarity between the spectra obtained at adjacent voxels. For example, Refs. [1, 2] made use of a Tikhonov regularization penalty of the form R(f) = N X n=1 X m2(n) 1 2 kf n f m k 2 2 , 1 2 kDfk 2 2 ; (4.9) where (n) is the set of indices for all voxel that are adjacent to the nth voxel, and D is a matrix representation of the spatial nite dierence operator. Unless otherwise noted, the remainder of this paper will assume the choice of R() given in Eq. (4.9). In principle, the Lawson-Hanson NNLS algorithm from the previous subsection could also be ap- plied to solve the spatially-regularized problem from Eq. (4.8). However, this requires substantially increasing the size of the NNLS problem from size Q to size QN, which leads to correspondingly- substantial increases in the amount of memory that is needed and decreases in convergence speed. In practice, it has been observed that the Lawson-Hanson algorithm is only worthwhile to use when the number of voxels N is quite small, which has motivated the use of modied problem formulations that only reconstruct small image patches instead of the entire image [157, 158, 159]. While this patch-based spatially-regularized approach can be somewhat computationally eective, the focus on image patches can introduce boundary artifacts at the edges of each patch that would not be observed had the entire image been reconstructed. As noted in the introduction, several large-scale algorithms that can estimate the entire image have also been proposed based on the ADMM algorithm [1, 2, 160]. Below, we present a quick sketch of the algorithm from Ref. [1], which is representative of the algorithms from the other references. This approach uses variable splitting to convert Eq. (4.8) into an equivalent form: ^ f = arg min f min x;y;z 1 2 km (I N K)xk 2 2 +I + (y) + 2 kDzk 2 2 ; (4.10) Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 65 subject to the constraints that f = x, f = y, and f = z, where x, y and z are the new auxiliary variables that have been introduced to simplify subproblems within the alternating ADMM proce- dure. In this expression,I + () has been introduced as the indicator function for the feasible set of the non-negativity constraint [174], such that I + (z) = 8 > < > : 0; z 0 elementwise; +1; otherwise. (4.11) Given this new formulation of the optimization problem, ADMM then alternates between updating the solutions to three separate subproblems, one for each of x, y, and z. For completeness, the full ADMM algorithm from Ref. [1] is reproduced in Algorithm 5. As can be seen, the algorithm needs to keep track of at least 8 vectors that are the same size as the spectroscopic image f, which can be a substantial amount of memory when Q and N are both large. Algorithm 5 Previous ADMM Algorithm [1] 1: Input: ADMM parameters . Problem specication parameters K, m, , D. 2: Initialize: Iteration number k = 0. Primal variables f k , x k , y k , z k . Dual variables d k x , d k y , d k z . 3: M, K T K +I Q 1 . 4: g, (I N K T )m. 5: while not converged do 6: f k+1 = x k + d k x +y k + d k y +z k + d k z =(3) 7: x k+1 = (I N M) g + f k+1 d k x = 8: y k+1 = max 0; f k+1 d k y = (maximization performed elementwise) 9: z k+1 = D T D +I 1 f k+1 d k z 10: d k+1 x = d k x (f k+1 x k+1 ) 11: d k+1 y = d k y (f k+1 y k+1 ) 12: d k+1 z = d k z (f k+1 z k+1 ) 13: k k + 1 14: end while Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 66 4.3 Proposed LADMM-Based Approach 4.3.1 Review of LADMM Our proposed approach is based on the LADMM algorithm [162, 163, 164, 165, 166, 167], and we begin this section with a quick review of LADMM. LADMM is designed to solve optimization problems of the form f^ x; ^ yg = arg min x2R S y2R T f(x) +g(y); s:t: Ax + By = b (4.12) where f : R S ! R[f+1g and g : R T ! R[f+1g are closed proper convex functions. The augmented Lagrangian form of Eq. (4.12) is L (x; y; d) =f(x) +g(y)hAx + By b; di + 2 kAx + By bk 2 2 : (4.13) where d is called a dual variable and > 0 is a penalty parameter that serves as a trade-o between the progress of the primal and dual updates. The conventional ADMM algorithm [161] would approach the optimization problem from Eq. (4.12) by alternatingly updating the primal variables x, y and the dual variable d according to 8 > > > > > > < > > > > > > : y k+1 = arg min y L (x k ; y; d k ) x k+1 = arg min x L (x; y k+1 ; d k ) d k+1 = d k (Ax k+1 + By k+1 b): (4.14) Using a simple completeing-the-square transformation and neglecting constant terms, the x and y subproblems of Eq. (4.14) can be further simplied to 8 > > < > > : y k+1 = arg min y g(y) + 2 By b Ax k + d k = 2 2 x k+1 = arg min x f(x) + 2 Ax b By k+1 + d k = 2 2 : (4.15) Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 67 Indeed, the previous ADMM approach (Algorithm 1) can be obtained by making the associations x2R S ! x T y T z T T 2R 3QN y2R T ! f2R QN A! I 3QN ; B! 2 6 6 6 6 4 I I I 3 7 7 7 7 5 ; b! 02R 3QN f(x) = 1 2 km (I N K)xk 2 2 +I + (y) + 2 kDzk 2 2 g(y)! 0: (4.16) As can be seen, the eciency of ADMM relies on the ease of solving the subproblems in Eq. (4.15), as well as the speed at which the iterations converge. In contrast, LADMM instead solves the subproblems 8 > > > > > > > < > > > > > > > : y k+1 = arg min y L (x k ; y; d k ) + 1 2 (y y k ) T Q(y y k ) x k+1 = arg min x L (x; y k+1 ; d k ) + 1 2 (x x k ) T P(x x k ) d k+1 = d k (Ax k+1 + By k+1 b); (4.17) where the main dierence from ADMM is the introduction of additional quadratic terms involving P and Q matrices which can be chosen to simplify the solution of the subproblems. Detailed discussion of choosing appropriate P and Q matrices is found in Ref. [167]. One of the standard choices in LADMM [162, 163, 164, 165, 166, 167] is to use P = p I S A T A and Q = q I T B T B ; (4.18) with proper choices of parameters p and q . Importantly, these choices cause the x and y subprob- lems from Eq. (4.17) to simplify substantially, in a way that cancels some of the A and B matrices Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 68 that appear in the optimization problem. Specically, again applying a simple completing-the- square transformation and neglecting constant terms, the x and y subproblems from Eq. (4.17) reduce to 8 > > < > > : y k+1 = arg min y g(y) + p 2 ky q k k 2 2 x k+1 = arg min x f(x) + q 2 kx p k k 2 2 ; (4.19) where q k = y k (= q )B T (Ax k + By k b d k =) p k = x k (= p )A T (Ax k + By k+1 b d k =): Interestingly, the two subproblems in Eq. (4.19) are equivalent to computing the proximal operators of the functions g and f. This is benecial because proximal operators are known in closed-form for a variety of interesting functions, including ` 1 -norm and the squared` 2 -norm [175, 176, ?], and enables very simple/fast solution of the subproblems. In our proposed approach, presented in the sequel, we will make dierent choices of the matrices P and Q than those given in Eq. (4.18), although our choice of P is motivated by Eq. (4.18) and the preceding discussion. Convergence of LADMM has been analyzed for dierent scenarios. In [162, 163], it was proven that LADMM will converge to a globally-optimal solution from an arbitrary initialization whenever P and Q are symmetric and positive denite. In [165], it was shown that the same kind of global convergence will be obtained whenever P and Q are symmetric and positive semi-denite. In [29], it was shown that the same kind of global convergence will also be obtained if P and Q are chosen as in Eq. (4.18), but with p > 0:75kA T Ak and q > 0:75kB T Bk (which enables the use of certain indenite P and Qmatrices), wherekk denotes the spectral norm. Several papers have also derived even less-restrictive conditions on P and Q that guarantee global convergence [29,31,32]. The algorithm we propose in the sequel makes use of a positive semi-denite P matrix and an indenite Q matrix that satises the constraints from [164], and thus inherits the global convergence guarantees from previous literature. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 69 4.3.2 Proposed Method Our proposed approach is based on applying LADMM principles to the optimization problem from Eq. (4.8). Instead of using three dierent variable splittings as was done in Eq. (4.10) for the previous ADMM algorithm [1], we propose to use a formulation based on a single variable splitting that is also equivalent to Eq. (4.8): ^ f = arg min f min z 1 2 km (I N K)fk 2 2 +I + (z) + 2 kDzk 2 2 ; (4.20) subject to the constraint that f = z. With this splitting, we can map the LADMM problem described in the previous subsection to this new optimization problem by making the associations: x2R S ! z2R QN ; y2R T ! f2R QN A! I N ; B!I N ; b! 02R QN f(x) =I + (x) +=2kDxk 2 2 ; and g(y)! 1=2km (I N K)yk 2 2 : (4.21) This directly results in the following set of LADMM update equations: f k+1 = arg min f 1 2 km (I N K)fk 2 2 D z k f; d k E + 2 z k f 2 2 + 1 2 (f f k ) T Q(f f k ); (4.22) z k+1 = arg min z I + (z) + 2 kDzk 2 2 D z f k+1 ; d k E + 2 kz f k+1 k 2 2 + 1 2 (z z k ) T P(z z k ); (4.23) and nally d k+1 = d k (z k+1 f k+1 ): (4.24) Our choices of P and Q and the methods we use to solve these subproblems are described in the following subsections. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 70 4.3.3 Solving the f subproblem It is not hard to show that the optimal solution to Eq. (4.22) is given by f k+1 = arg min f 1 2 km (I N K)fk 2 2 + 2 f z k + d k = 2 2 + 1 2 (f f k ) T Q(f f k ): (4.25) Notably, if we choose Q = 0 (which is positive semi-denite), then the resulting problem ends up having a decoupled (separable) block-diagonal structure, which allows the solution to be computed in a simple voxel-by-voxel manner. Compared to the voxel-by-voxel NNLS reconstruction problem from Eq. (4.7), the f subproblem is even simpler because it does not involve any nonnegativity constraints, which allows the optimal solution to be computed in closed form as: f k+1 n = M K T m n +z k n d k n ; (4.26) where z k n 2R Q and d k n 2R Q are respectively the components of z k and d k corresponding to the nth voxel position, and the matrix M2R QQ , K T K +I Q 1 : Note that the matrix M is the same for all voxels n and all iterations k, which allows it to be precomputed once and then reused whenever it is needed [1]. This approach is viable and ecient when Q is small. However, in scenarios where Q is large (e.g., in the multidimensional multiparametric setting [1, 2, 138, 139, 140, 141, 142, 143]), storing and performing matrix-vector multiplications with a QQ matrix can still be relatively computationally burdensome, requiring O(Q 2 ) memory for storage and O(Q 2 ) oating point operations for matrix-vector multiplication. The computational complexity associated with M can be substantially reduced by leveraging the fact that while the matrix K T K is of size QQ, it generally has rank that is smaller than Q. Specically, the matrix will have rank that is no greater than P under the assumption that P < Q (which is generally the case in multidimensional multiparametric scenarios). In practice, it is also common that the dictionary K has columns that are highly-correlated (approximately linearly dependent), allowing it to be accurately approximated by a matrix of rank-r with r < P [152, 177]. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 71 Let K r 2 R PQ be a rank-r approximation of K, where the optimal low-rank approximation is obtained by truncating the singular valued decomposition (SVD) of K. If we assume that the SVD of K r is given by K r = r X i=1 i u i v T i ; (4.27) with singular vectors u i 2 R P and v i 2 R Q and singular values i > 0 for i = 1;:::;r, then it is straightforward to derive that matrix-vector multiplications involving the matrix M r , K T r K r +I Q 1 associated with the rank-r approximation of the dictionary can be calulated simply as M r x = 1 x r X i=1 2 i 2 + 2 i v i v T i x (4.28) for arbitrary vectors x2 R Q . With this approximate representation, storing M r requires only O(rQ) memory (for storing the vectors v i and singular values i ) andO(rQ) oating point opera- tions for computing matrix-vector multiplication, which can be a substantial improvement over the naive O(Q 2 ) approach associated with directly using the matrix M. 4.3.4 Solving the z subproblem It is also not hard to show that the z subproblem from Eq. (4.23) is given by z k+1 = arg min z I + (z) + 2 kz f k+1 d k =k 2 2 + 2 kDzk 2 2 + 1 2 (z z k ) T P(z z k ): (4.29) In this case, the nite dierence operator D causes coupling of the dierent spatial locations, and precludes the use of a voxel-by-voxel solution. In previous ADMM approaches [1, 2, 160], this issue is addressed by introducing more splitting variables, which can increase memory consumption and reduce convergence speed. In our proposed approach, motivated by Eq. (4.18) and the associated discussion from Sec. 4.3.1, we choose the matrix P to simplify the problem so that we can solve for z k+1 non-iteratively. Specically, we choose P = p I QN D T D; (4.30) Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 72 with some parameter p to be determined later. This choice causes full decoupling of every entry of z k+1 . With this choice, if we ignore the nonnegativity constraint on the elements of z k+1 , the optimal solution for the z subproblem becomes z k+1 = 1 p + p z k D T Dz k +f k+1 + d k ; (4.31) which can be computed easily and noniteratively. Because of the simple decoupled structure we obtain from this choice of P, imposing the nonnegativity constraint is as simple as taking the value of z k+1 from Eq. (4.31), while replacing all of the negative entries of z k+1 with 0. The complete description of the proposed algorithm is summarized in Algorithm 6. Algorithm 6 Proposed LADMM 1: Input: LADMM parameters , p , r. Problem specication parameters K, m, , D. 2: Initialize: Iteration number k = 0. Primal variables f k , z k . Dual variables d k . 3: Compute v i and i according to Eq. (4.27) 4: g n , K T m n , for n = 1; 2;:::;N. 5: while not converged do 6: 7: for n = 1; 2;:::;N do 8: c k n = g n +z k n d k n 9: f k+1 n = 1 c k n P r i=1 2 i 2 + 2 i v i v T i c k n 10: end for 11: z k+1 = max 0; pz k D T Dz k +f k+1 +d k p+ 12: d k+1 = d k (z k+1 f k+1 ) 13: k k + 1. 14: end while Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 73 4.3.5 Complexity Considerations Given the choices described in the previous subsections, our proposed LADMM approach uses substantially fewer variables and has substantially fewer subproblems than the previous ADMM approach from Algorithm 5 [1, 2]. In addition, the proposed approach benets from a more ecient matrix inversion for the f subproblem based on the low-rank properties of K, and eliminates the need for matrix inversion for the z subproblem. 4.3.6 Practical Parameter Selection Use of the proposed algorithm requires choosing the values of the parameters p , , and r. 4.3.6.1 The parameter p As described previously, convergence of Algorithm 6 is guaranteed if p is chosen such that p > 0:75kD T Dk [164]. As a result, we use the following setting p = 0:75kD T Dk + (4.32) where is a small positive number (we use = 10 10 ). 4.3.6.2 The penalty parameter Although LADMM will converge to a globally optimal solution regardless of the choice of , the choice of still has a big impact on the convergence speed [161]. Unfortunately, optimal methods for selecting beta are only known for a few special cases [178]. Some adaptive methods for selecting beta have been proposed [179, 180], although these approaches require additional computations that are not attractive for the type of large-scale problem considered in this work. In this work, we use a simple heuristic to select in a context-dependent yet computationally-ecient way. Specically, we rst consider the problem of estimating a small image patch (e.g., 3 3 voxels), and tune beta to maximize convergence speed for the small-scale problem. We then use this same value of beta to solve the full-scale optimization problem. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 74 4.3.6.3 The rank parameter r The choice of the rank parameter r represents a trade-o between computational eciency and the accuracy of the dictionary approximation. For the application examples described in the next sample, although the dictionaries in each case are based on dierent signal models and have dierent values ofQ andP , we observed that the dictionaries for each case were all well-approximated if we choser = 15 (with less than 0:06% approximation error in each case, as measured in the Frobenius norm). 4.4 Results In this section, we evaluate the proposed LADMM algorithm in various CSI scenarios by comparing it to the previous ADMM algorithm proposed in [1]. For convergence speed comparison, we use two metrics: One is the cost function value and the other is the distance from the converged solution (DFCS) dened as DFCS(f k ), kf k f k 2 kf k 2 (4.33) where f k and f are the spectrum variable f at iterationk and nal exit of Algorithm 6 respectively. Due to the high correlation among the dictionary basis functions, we often encounter the cases where very dierent spectrum variables would produce similar cost function values. For this reason, dNRMSE would be a good supplement to evaluate the convergence speed of algorithms. Note that in practice, it makes no sense to use dNRMSE as a stopping criteria as we need to run the algorithm twice. For a fair comparison, we applied the aforementioned small patch based tuning strategy to set the penalty parameters for both of the two compared algorithms. Both algorithms were imple- mented in MATLAB R2022a (The MathWorks, Inc., Natick, Massachusetts, United States.). All the main numerical experiments were run on a AMD EPYC CPU equipped with 32 cores. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 75 4.4.1 Diusion-Relaxation In this subsection, we consider an ex-vivo mouse spinal cord diusion relaxation CSI data that has been reported in [1] (control subject 1). This DR-CSI dataset has 7 b values in the same diusion direction and 4 echo times for a total of N t = 28 number of contrast weightings. We use the same dictionary K as in [1] that is constructed using 70 diusion coecients D (logarithmically distributed from 0:01 to 5 m 2 /ms), and 70 T 2 values (logarithmically distributed from 3 to 300 ms) for a total of N b = 4900 basis functions. The image matrix size is 32 31 and the spatial regularization parameter is set as = 0:001 in this case. As shown in Fig. 4.2, the proposed LADMM converges much faster than the previous ADMM in both cost function values and the dNRMSEs. However, similar to the MR ngerprinting dataset, the ill-conditioning in this case is so bad that the spectrum variables as indicated by dNRMSEs are still changing though the cost functions seem to be convergent very quickly. The component spectrum locations and spatial maps are shown in Fig. 4.1. We also compare the component spectrum and spatial maps results of the spatially regularized reconstruction (second row in Fig. 4.1) generated by the proposed LADMM algorithm and those of NNLS method (rst row in Fig. 4.1). Consistent with [1], three components were identied where components 1-3 seem to correspond to white matter structures, gray matter structures and dorsal gray matter dominated structures respectively. The comparison between the two rows of component spatial maps in Fig. 4.1 show that spatial regularization greatly improves noise suppression and results in more continuous spatial maps. 4.4.2 IR-MSE Relaxation-Relaxation In this subsection, we consider a 2D slice in-vivo human brain RR-CSI data that has been reported in [1] (subject 4) and acquired with an inversion recovery multi spin echo (IR-MSE) sequence. This dataset comes with 7 inversion times and 15 echo times for a total ofN t = 105 contrast weightings. The same dictionary K as in [1] is constructed with 100 T 1 values (ranging from 100 3000 ms spaced logarithmically) and 100 T 2 values (ranging from 2 300 ms spaced logarithmically) for Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 76 NNLS (a) Spectrum (b) Composite (c) C1 (d) C2 (e) C3 LADMM (Spatially Regularized) (f) Spectrum (g) Composite (h) C1 (i) C2 (j) C3 Figure 4.1: Component spectrum locations (a),(f) and spatial maps (b)-(e), (g)-(j) in the DR- CSI mouse spinal cord experiments. Consistent with [1], three components were identied where components 1-3 seem to correspond to white matter structures, gray matter structures and dorsal gray matter dominated structures respectively. The rst and second rows show results with and without spatial regularization respectively where we can see spatial regularization greatly improves noise suppression and results in more continuous spatial maps. (a) (b) Figure 4.2: Convergence curves in the DR-CSI experiments. The proposed LADMM converges 4x-5x faster than the previous ADMM in terms of cost function values. Even more signicant acceleration was observered in terms of dNRMSEs. Also noticeably, the ill-conditioning of the inverse problem in this case causes the dNRMSEs be ever changing without being convergent for both algorithms even after a long running time. a total of N b = 10000 number of basis functions. The image matrix size is 86 66. We set the regularization parameter as = 0:01 according to [2]. As shown in Fig. 4.6, the cost function values of the proposed LADMM converges more than two times faster than that of the previous ADMM. We also notice that due to the ill-conditioning of the problem, the spectrum variables converge much slower than the cost functions for both algorithms. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 77 However, dNRMSEs indicate that the proposed LADMM converges even more signicantly faster than the previous ADMM. Consistent with [2], six anatomically plausible compartments are identied for both compared algorithms. The spectrum locations and spatial maps of the six compartments are shown in Fig. 4.4 and Fig. 4.3 respectively. As can be seen, both the spectra and spatial maps results are almost the same for the two compared algorithms as expected. For more discussion on the interpretation of the six compartments please refer to [2]. Previous ADMM (a) Com- posite (b) Comp. 1 (c) Comp. 2 (d) Comp. 3 (e) Comp. 4 (f) Comp. 5 (g) Comp. 6 Proposed LADMM (h) Com- posite (i) C1 (j) C2 (k) C3 (l) C4 (m) C5 (n) C6 Figure 4.3: Component spatial maps comparison in the IR-MSE based RR-CSI experiments. No noticeable dierences were observed between the two compared algorithms. As interpreted in [2], component 1 to 6 seem to correspond to white matter structures, gray matter structures with relatively high myelin content, brain structures with high iron content, gray matter without myelin, cerebrospinal uid (CSF) and myelin water respectively. The corresponding spectrum locations of each component are shown in Fig. 4.4. 4.4.3 Multi Spin Echo T 2 Relaxation In this subsection, we consider a T 2 relaxation only data extracted from the previous IR-MSE dataset with certain inversion time (TI = 400ms). The dictionary size in this case is reduced to N b = 100. As shown in Fig. 4.7, the proposed LADMM converges at least two times faster than the previous ADMM in terms of both cost function values and dNRMSEs. The compartment spectrum and spatial maps of the proposed LADMM algorithm are shown in Fig. 4.5. We can see that only three compartments are identied which is much less than that in the RR-CSI experiment in the Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 78 (a) Previous ADMM (b) Proposed LADMM Figure 4.4: Component spectrum overlay comparisons in the IR-MSE based RR-CSI experiments. C1-C6 stand for components 1 to 6 respectively. No noticeable dierences were observed between the two compared algorithms. The corresponding spatial maps of each component are shown in Fig. 4.3. (a) Spectrum (b) Composite (c) Comp. 1 (d) Comp. 2 (e) Comp. 3 Figure 4.5: Component spectrum overlay (a) and spatial maps (b)-(e) in the multi spin echo based T 2 -relaxation only experiments. Only three components were identied in this experiment which is in stark comparison to the IR-MSE based two dimensional RR-CSI experiment where six components were captured. (a) (b) Figure 4.6: Convergence curves of the IR-MSE based RR-CSI experiments. Cost function conver- gence curves show 2-3x faster of speed of the proposed LADMM over the previous ADMM. The acceleration is even more signicant in terms of dNRMSEs. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 79 (a) (b) Figure 4.7: Convergence curves in the multi spin echo basedT 2 -relaxation only experiments. About 2x faster convergence speed was observed for the proposed LADMM over the previous ADMM in terms of both cost function values and dNRMSEs. previous subsection. This shows an example of the much improve capability in resolving partial volume eects of the multidimensional multiparametric acquisition method over its one dimensional compartment. 4.4.4 MR Fingerprinting Relaxation-Relaxation In this subsection, we consider a new RR-CSI dataset that was acquired using a modied MR ngerprinting sequence with a ViSTa preparation block to suppress the long T 1 components and thus to emphasize the short T 1 components such as myelin water [181]. A fully sampled (32 spiral interleaves to cover the full k-space) 2D slice in-vivo human brain dataset was acquired with 1 mm in-plane spatial resolution and 10 mm slice thickness in a 3T Siemens Prisma scanner. A total of N t = 540 time points of images with varying ip angles were acquired. The image matrix size is 220 220. To account for theB 1 ( ip angle) inhomogeneity, we acquired aB 1 map and discretized the map using 31 variation values ranging from 0:7 1:3 (baseline is 1 indicating no ip angle variation) with a constant step 0:02. For each of the 31 B 1 variations, we construct a dictionary using 100 T 1 values (ranging from 100 3000 ms spaced logarithmically), 100 T 2 values (ranging from 2300 ms spaced logarithmically) for a total ofN b = 10000 number of basis functions through extended-phase-graph [182] simulations. We thus end up with 31 dictionaries. It is straightforward to incorporate B 1 correction into both algorithms as the dictionary based operations are already spatially separable. Note that the proposed LADMM consumes less than half of the memory Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 80 compared to the previous ADMM, for example in this case, the peak memory consumption are 26GB versus 68GB. The spatial regularization parameter is set as = 0:0001. For this dataset, we were able to identify six anatomically plausible components whose spectrum locations and spatial maps are shown in Fig. 4.9. Although the association of these components to brain structures are beyond the scope of this paper, components 1-3 look like white matter, gray matter and CSF respectively, component 5 looks like myelin water, component 4 could be venous and component 6 may be the scalp. As shown in Fig 4.8, the proposed LADMM converges faster than the previous ADMM in both cost function values and dNRMSEs. However, we do notice that, the ill-conditioning in this case is so bad that even though the cost functions seem to be convergent, the spectrum variables as indicated by dNRMSEs are still changing for both algorithms. (a) (b) Figure 4.8: Convergence curves in the ViSTa MRF based RR-CSI experiments. Faster convergence of the proposed LADMM over the previous ADMM were observed for both the cost function values and dNRMSEs. As noticed in the dNRMSEs, the conditioning of the inverse problem is so bad in this case that the dNRMSEs are still changing for both algorithms even after a long running time especially for the previous ADMM. 4.5 Discussion In the proposed Algorithm 6, besides the LADMM framework that enables greatly reduced number of auxiliary variables for splitting subproblems, another key component in convergence acceleration is due to the SVD based approach in Eq. (4.28). To see this, we compare the convergence speed of two versions of the proposed LADMM: One applied Eq. (4.28) and the other applied the original Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 81 (a) Spectrum (b) Composite (c) C1 (d) C2 (e) C3 (f) C4 (g) C5 (h) C6 Figure 4.9: Component spectrum locations and spatial maps in the ViSTa MRF based RR-CSI experiments. Six anatomically plausible components were identied where components 1-6 seem to correspond to white matter structures, gray matter structures, CSF, venous, myelin water and scalp respectively. matrix inversion based implementation in Eq. (4.26) for the f-subproblem in the IR-MSE based RR-CSI experiment. As shown in Fig. 4.10, the convergence speed of the proposed LADMM is accelerated by 2x-3x by applying Eq. (4.28) instead of Eq. (4.26) for f-subproblem while both versions are faster than the previous ADMM algorithm. Notably, in the proposed LADMM algorithm, the f-subproblem can be implemented voxel-by- voxel independently while the z-subproblem can be implemented spectrum-by-spectrum indepen- dently. This separability in implementations indicate that the proposed LADMM can be further accelerated by parallel computing which is subject to future research. Although in this work we considered the specic Tikhonov based spatial regularization in Eq. (4.9), the LADMM framework and the SVD based approach in Eq. (4.28) can still t for other regularization functions such as ` 1 and total variations. Although, these explorations are beyond the scope of this paper. Chap.5 An Ecient Algorithm for Spatial-Spectral Partial Volume Compartment Mapping with Applications to Multicomponent Diusion and Relaxation MRI 82 Figure 4.10: Convergence comparison of the proposed LADMM w/ Eq. (4.26) and Eq. (4.28) in the IR-MSE based RR-CSI experiment. 4.6 Conclusion We proposed an ecient algorithm to solve the spatial-spectral partial volume compartment map- ping problem. Our proposed algorithm is based on the LADMM framework that only introduces one auxiliary variable for subproblem splittings. We also proposed a low rank SVD based approach for a much faster subproblem solver in the proposed LADMM algorithm. We have shown the signif- icant improve in empirical convergence speed of the proposed LADMM algorithm over the previous ADMM algorithm in various correlation spectroscopic imaging scenarios. Parallel computing could potentially further accelerate the proposed algorithm. Chapter 5 Conclusion The non-invasiveness and rich contrast information of MRI make it a unique and powerful tool in both biomedical research and clinical diagnosis. Although MRI has been existed for nearly ve decades, its full potentials in applications are far from throughly explored largely due to the trade- os among data acquisition time, image SNR and resolution. Complementary to imaging hardware improvement and advanced pulse sequence development, constrained MRI has greatly increased our abilities in dealing with these trade-os by utilizing various prior knowledge of the underlying MRI images to improve data acquisition eciency and image reconstruction/analysis performance. Mathematical optimization has been widely used to model and solve constrained MRI problems. However, the application of mathematical optimization to constrained MRI is often highly non- trivial as MRI dataset often has large size and computational eciency must be considered while designing algorithms for solving such optimization problems. In this dissertation, we have shown three novel methods in constrained MRI with mathematical optimization modeling and algorithms. In Chapter 2, we proposed a phase constrained data acquisition method for non-Fourier RF encoding acquisition. The idea is inspired by the commonly used partial Fourier technique where smooth phase characteristic enables subsampling in k-space due to the Fourier relation. We have shown from estimation theory that smooth phase can also accelerate non-Fourier RF encoding data acquisition with the proposed encoding scheme designs. We also proposed a phase constrained 83 Chap.6 Conclusion 84 decoding method for image reconstruction which we tested in the context of gSlider, a diusion MRI acquisition method, with simulated RF encoded data. In Chapter 3, we proposed a fast and ecient algorithm named PALMNUT for constrained image reconstruction with separate magnitude and phase modeling which has been shown in liter- ature to be more appropriate and/or result in better image qualities in certain scenarios but the resulting optimization problems become much more challenging. Our proposed algorithm is based on a novel application of an existing algorithm called PALM together with two additional features for further acceleration including Nesterov's momentum and a coordinate-wise step size settings strategy. Theoretically, we have shown the convergence guarantee for the proposed algorithm when Nesterov's momentum is not included. Empirically, we compared PALMNUT with two state-of- the-art methods in various MRI image reconstruction tasks and the results show that PLAMNUT can be up-to hundreds times faster than the compared methods. In Chapter 4, we proposed a novel algorithm for a constrained image analysis task of multi compartment mapping with joint spatial-spectral estimation which has been shown in the literature as a challenging ill-posed inverse problem. The proposed algorithm is based on a novel application of the linearized ADMM together with a SVD based ecient implementation of the solution of one of the subproblems. We compared the proposed algorithm with state-of-the-art method in multiple diusion and relaxation MRI settings. Results show the proposed algorithm converges drastically faster than the state-of-the-art. For general optimization in the applied math community, one of the major trends is to de- velop ecient algorithms for large scale problems as data generation and collection become easier nowadays. The MRI application also falls in this trend and directly benets from the development in general mathematical optimization. The algorithms we proposed in this dissertations are good examples of how we can extract nutrients from applied math community. This provides us great opportunities to improve constrained MRI just by `translating' what has been available in applied math community. However, we should be aware of the fact that when researchers in the applied math community develop an algorithm, they often try to make it ts optimization problems that are as general as possible. On the other hand, truly ecient algorithms should try to take advantage Chap.6 Conclusion 85 of the problem specic structure as much as possible. Therefore, general optimization algorithms could be a very good starting point while specic designs are often necessary to develop ecient algorithms for constrained MRI. Although this sounds like an easy formula for algorithm develop- ment in constrained MRI, this process is highly non-trivial and requires deep understanding of both elds. Take an example of the development of PALMNUT, we need to understand MRI physics to rstly be able to come up with an appropriate optimization model and then identies PALM would be a good t for this problem structure. However, as our breakdown experiments show, PALM in itself converges quite slowly. We then observed that the challenging phase subproblem can be much more eciently solved by using the proposed coordinate-wise step sizes strategy that better ts to the special problem structure. With additional momentum acceleration, we then developed PALMNUT which converges up to hudreds times faster than state-of-the-art methods. Appendix A Proof of Thm. 3.1 Theorem 3.1 is an immediate consequence of the following Lemma. Lemma A.1. LetQ(x) :R N !R be smooth, and assume that a non-negative vector L2R N exists such that hrQ(x 1 )rQ(x 2 ); x 1 x 2 i p L (x 1 x 2 ) 2 2 (A.1) for8x 1 ; x 2 2R N . Then for8x; ^ x k 2R N and8c L (elementwise), we have Q(x)Q(^ x k ) +hx ^ x k ;rQ(^ x k )i + 1 2 p c (x ^ x k ) 2 2 : (A.2) Further, the right hand side of Eq. (A.2) can be written compactly as Q(^ x k ) +hx ^ x k ;rQ(^ x k )i + 1 2 p c (x ^ x k ) 2 2 = k + 1 2 p c (x w k ) 2 2 ; (A.3) where w k = ^ x k diag(c) 1 rQ(^ x k ) (A.4) 86 Appendix A 87 and k is a constant that does not depend on the variable x, given by k =Q(^ x k )h^ x k ;rQ(^ x k )i + 1 2 p c ^ x k 2 2 1 2 p c w k 2 2 : (A.5) Proof. Inspired by the proof of Proposition A.24 from [11], dene (t),Q(tx + (1t)^ x k ). Then (0) =Q(^ x k ), (1) =Q(x), and d dt (t) =hrQ(tx + (1t)^ x k ); x ^ x k i: (A.6) We have that Q(x)Q(^ x k ) = Z 1 0 d dt (t)dt = Z 1 0 hrQ(tx + (1t)^ x k ); x ^ x k idt = Z 1 0 hrQ(tx + (1t)^ x k )rQ(^ x k ); x ^ x k idt +hrQ(^ x k ); x ^ x k i = Z 1 0 1 t hrQ(tx + (1t)^ x k )rQ(^ x k );txt^ x k idt +hx ^ x k ;rQ(^ x k )i = Z 1 0 1 t hrQ(tx + (1t)^ x k )rQ(^ x k );tx + (1t)^ x k ^ x k idt +hx ^ x k ;rQ(^ x k )i Z 1 0 1 t p L (tx + (1t)^ x k ^ x k ) 2 2 dt +hx ^ x k ;rQ(^ x k )i = Z 1 0 t p L (x ^ x k ) 2 2 dt +hx ^ x k ;rQ(^ x k )i = 1 2 p L (x ^ x k ) 2 2 +hx ^ x k ;rQ(^ x k )i 1 2 p c (x ^ x k ) 2 2 +hx ^ x k ;rQ(^ x k )i: (A.7) This derivation proves Eq. (A.2), while the simplications leading to Eq. (A.3) come simply from completing the square. Appendix B Derivation of Eq. (2.6) We rst split the real and imaginary parts of the data acquisition model in Eq. (??) that would result in only real valued vectors/matrices 2 6 4 d R d I 3 7 5 = 2 6 4 g R (m;) g I (m;) 3 7 5 + 2 6 4 R I 3 7 5 (B.1) where g R , (I G R )m cos(1 ) (I G I )m sin(1 ) g I , (I G R )m sin(1 ) + (I G I )m cos(1 ): So the Jacobian matrix becomes J g (m;) = 2 6 4 J g R (m) J g R () J g I (m) J g I () 3 7 5 (B.2) Now, let's calculate each of the sub-Jacobian matrix: J g R (m) = diagfcos(1 )g (I G R ) diagfsin(1 )g (I G I ) =I A (B.3) 88 Appendix B 89 where A , diagfcosg G R diagfsing G I , J g R () = 2 6 6 6 6 4 : : : diagfG R m n sin + G I m n cosg : : : 3 7 7 7 7 5 ; (B.4) J g I (m) = diagfsin(1 )g (I G R ) + diagfcos(1 )g (I G I ) =I B (B.5) where B , diagfsing G R + diagfcosg G I , J g R () = 2 6 6 6 6 4 : : : diagfG R m n cos G I m n sing : : : 3 7 7 7 7 5 : (B.6) So, according to Eq. (??) the FIM is I() = J g T J g = 2 6 4 M 11 M 12 M T 12 M 22 3 7 5 (B.7) where M 11 , J g R (m) T J g R (m) + J g I (m) T J g I (m) = I A T A + B T B = G T R G R + G T I G I =<fG H Gg; (B.8) M 22 , J g R () T J g R () + J g I () T J g I () = diag 8 < : Nxy X n=1 jG R m n j 2 +jG I m n j 2 9 = ; = diag 8 < : Nxy X n=1 jGm n j 2 9 = ; (B.9) Appendix B 90 and M 12 = J g R (m) T J g R () + J g I (m) T J g I () = B T n T (B.10) with B n ,A T diagfG R m l sin() + G I m l cos()g + B T diagfG R m l cos() G I m l sin()g = G T I diagfG R m l g G T R diagfG I m l g == G H diagfGm n g : Collecting all the four blocks into the FIM, we obtain Eq. (2.6). Bibliography [1] D. Kim, E. K. Doyle, J. L. Wisnowski, J. H. Kim, and J. P. Haldar, \Diusion-relaxation correlation spectroscopic imaging: A multidimensional approach for probing microstructure," Magn. Reson. Med., vol. 78, pp. 2236{2249, 2017. [2] D. Kim, J. L. Wisnowski, C. T. Nguyen, and J. P. Haldar, \Multidimensional correlation spectroscopic imaging of exponential decays: From theoretical principles to in vivo human applications," NMR Biomed., vol. 33, p. e4244, 2020. [3] P. C. Lauterbur, \Image formation by induced local interactions: Examples employing nuclear magnetic resonance," vol. 242, no. 5394, pp. 190{191, 1973. [4] Z.-P. Liang and P. C. Lauterbur, Principles of Magnetic Resonance Imaging:A Signal Pro- cessing Perspective. Wiley-IEEE Press, 2000. [5] E. M. Haacke, R. W. Brown, M. R. Thompson, and R. Venkatesan, Magnetic Resonance Imaging: Physical Principles and Sequence Design. John Wiley & Sons, Inc., 1999. [6] D. G. Nishimura, Principles of Magnetic Resonance Imaging, 2nd ed. Stanford University, 2016. [7] Z.-P. Liang, F. E. Boada, R. T. Constable, E. M. Haacke, P. C. Lauterbur, and M. R. Smith, \Constrained reconstruction methods in MR imaging," Rev. Magn. Reson. Med., vol. 4, pp. 67{185, 1992. [8] J. P. Haldar, \Constrained imaging: Denoising and sparse sampling," Ph.D. dissertation. [9] M. Lustig, D. Donoho, and J. M. Pauly, \Sparse mri: The application of compressed sensing for rapid mr imaging," Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182{1195, 2007. [10] J. A. Fessler, \Model-based image reconstruction for MRI," IEEE Signal Processing Maga- zine, vol. 27, no. 4, pp. 81{89, 2010. [11] D. P. Bertsekas, Nonlinear Programming, 2nd ed. Athena Scientic, 1999. [12] A. Beck, First-Order Methods in Optimization. SIAM, 2017. [13] F. H. Clark, Optimization and Nonsmooth Analysis. SIAM, 1990. [14] P. Combettes and V. Wajs, \Signal recovery by proximal forward-backward splitting," Mul- tiscale Modeling & Simulation, vol. 4, no. 4, pp. 1168{1200, 2005. 91 Bibliography 92 [15] D. R. Hunter and K. Lange, \A tutorial on MM algorithms," Am. Stat., vol. 58, pp. 30{37, 2004. [16] N. Parikh and S. Boyd, \Proximal algorithms," Found. Trends Optim., vol. 1, no. 3, p. 127239, 2014. [17] P. L. Combettes and J.-C. Pesquet, \Proximal splitting methods in signal processing," in Fixed-Point Algorithms for Inverse Problems in Science and Engineering, H. H. Bauschke, R. S. Burachik, P. L. Combettes, V. Elser, D. R. Luke, and H. Wolkowicz, Eds. New York, NY: Springer New York, 2011, pp. 185{212. [18] H. Attouch, J. Bolte, and B. F. Svaiter, \Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss- Seidel methods," Math. Program., vol. 137, no. 1, pp. 91{129, 2013. [19] J. Bolte, S. Sabach, and M. Teboulle, \Proximal alternating linearized minimization for nonconvex and nonsmooth problems," Math. Program., vol. 146, pp. 459{494, 2014. [20] D. Gabay and B. Mercier, \A dual algorithm for the solution of nonlinear variational problems via nite element approximation," Computers & Mathematics with Applications, vol. 2, no. 1, pp. 17{40, 1976. [21] R. Glowinski, Numerical Methods for Nonlinear Variational Problems. Springer, 1984. [22] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, \Distributed optimization and statistical learning via the alternating direction method of multipliers," Found. Trends Mach. Learn., vol. 3, no. 1, p. 1122, 2011. [23] D. Gabay, \Chapter ix applications of the method of multipliers to variational inequalities," in Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems, ser. Studies in Mathematics and Its Applications, M. Fortin and R. Glowinski, Eds. Elsevier, 1983, vol. 15, pp. 299{331. [24] A. Themelis and P. Patrinos, \Douglas{Rachford splitting and ADMM for nonconvex op- timization: Tight convergence results," SIAM Journal on Optimization, vol. 30, no. 1, pp. 149{181, 2020. [25] R. T. Rockafellar, \Monotone operators and the proximal point algorithm," SIAM Journal on Control and Optimization, vol. 14, no. 5, pp. 877{898, 1976. [26] J. Eckstein and D. P. Bertsekas, \On the DouglasRachford splitting method and the proximal point algorithm for maximal monotone operators," Mathematical Programming, vol. 55, no. 1, pp. 293{318, 1992. [27] M. L. N. Goncalves, J. G. Melo, and R. D. C. Monteiro, \Convergence rate bounds for a proximal ADMM with over-relaxation stepsize parameter for solving nonconvex linearly constrained problems," arXiv e-prints, p. arXiv:1702.01850, 2017. [28] M. Hong, Z.-Q. Luo, and M. Razaviyayn, \Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems," in 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2015, pp. 3836{3840. Bibliography 93 [29] E. Ghadimi, A. Teixeira, I. Shames, and M. Johansson, \Optimal parameter selection for the alternating direction method of multipliers (ADMM): Quadratic problems," IEEE Transac- tions on Automatic Control, vol. 60, no. 3, pp. 644{658, 2015. [30] B. S. He and S. L. Wang, \Alternating direction method with self-adaptive penalty parameters for monotone variational inequalities," Journal of Optimization Theory and Applications, vol. 106, no. 2, p. 337356, 2000. [31] H. Raguet, J. Fadili, and G. Peyr, \A generalized forward-backward splitting," SIAM Journal on Imaging Sciences, vol. 6, no. 3, pp. 1199{1226, 2013. [32] M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, \Fast image recovery using variable splitting and constrained optimization," IEEE Transactions on Image Processing, vol. 19, no. 9, pp. 2345{2356, 2010. [33] D. S. Weller, S. Ramani, and J. A. Fessler, \Augmented Lagrangian with variable splitting for faster non-cartesian` 1 -SPIRiT MR image reconstruction," IEEE Transactions on Medical Imaging, vol. 33, no. 2, pp. 351{361, 2014. [34] C. Chen, B. He, Y. Ye, and X. Yuan, \The direct extension of ADMM for multi-block convex minimization problems is not necessarily convergent," Mathematical Programming, vol. 155, no. 1, p. 5779, 2014. [35] M. Hong and Z.-Q. Luo, \On the linear convergence of the alternating direction method of multipliers," Mathematical Programming, vol. 162, no. 1, p. 165199, 2016. [36] S. Becker, J. Bobin, and E. J. Cand es, \NESTA: A fast and accurate rst-order method for sparse recovery," SIAM J. Imaging Sci., vol. 4, pp. 1{39, 2011. [37] Z. Tan, Y. C. Eldar, A. Beck, and A. Nehorai, \Smoothing and decomposition for analysis sparse recovery," IEEE Transactions on Signal Processing, vol. 62, no. 7, pp. 1762{1774, 2014. [38] Y. Nesterov, \Smooth minimization of non-smooth functions," Math. Program., vol. 103, pp. 127{152, 2005. [39] H. Raguet and L. Landrieu, \Preconditioning of a generalized forward-backward splitting and application to optimization on graphs," SIAM Journal on Imaging Sciences, vol. 8, no. 4, pp. 2706{2739, 2015. [40] T. Pock and A. Chambolle, \Diagonal preconditioning for rst order primal-dual algorithms in convex optimization," in 2011 International Conference on Computer Vision, 2011, pp. 1762{1769. [41] F. Ong, M. Uecker, and M. Lustig, \Accelerating non-Cartesian MRI reconstruction conver- gence using k-space preconditioning," IEEE Transactions on Medical Imaging, vol. 39, no. 5, pp. 1646{1654, 2020. [42] K. Koolstra, J. van Gemert, P. Brnert, A. Webb, and R. Remis, \Accelerating compressed sensing in parallel imaging reconstructions using an ecient circulant preconditioner for Cartesian trajectories," Magnetic Resonance in Medicine, vol. 81, no. 1, pp. 670{685, 2019. Bibliography 94 [43] Y. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course. Boston: Springer, 2004. [44] A. Beck and M. Teboulle, \A fast iterative shrinkage-thresholding algorithm for linear inverse problems," SIAM J. Imaging Sci., vol. 2, pp. 183{202, 2009. [45] W. Su, S. Boyd, and E. Candes, \A dierential equation for modeling Nesterov's accelerated gradient method: Theory and insights," in Proc. NeurIPS, 2014, pp. 2510{2518. [46] T. Goldstein, B. O'Donoghue, S. Setzer, and R. Baraniuk, \Fast alternating direction opti- mization methods," SIAM Journal on Imaging Sciences, vol. 7, no. 3, pp. 1588{1623, 2014. [47] A. Kumar, D. Welti, and R. R. Ernst, \Nmr fourier zeugmatography," Journal of Magnetic Resonance (1969), vol. 18, no. 1, pp. 69{83, 1975. [48] Y. Liu, K. Setsompop, and J. P. Haldar, \Accelerating gSlider-based diusion MRI: Phase constraints enabled reduced RF encoding," in Proc. Int. Soc. Magn. Reson. Med., 2021, p. 1179. [49] Z.-P. Liang and P. C. Lauterbur, Principles of Magnetic Resonance Imaging: A Signal Pro- cessing Perspective. New York: IEEE Press, 2000. [50] J. P. Haldar and Z.-P. Liang, \\Early" constrained reconstruction methods," in Magnetic Res- onance Image Reconstruction: Theory, Methods, and Applications, M. Doneva, M. Akcakaya, and C. Prieto, Eds. Academic Press, 2022. [51] P. Margosian, F. Schmitt, and D. Purdy, \Faster MR imaging: imaging with half the data," Health Care Instrum., vol. 1, pp. 195{197, 1986. [52] D. A. Feinberg, J. D. Hale, J. C. Watts, L. Kaufman, and A. Marks, \Halving MR imaging time by conjugation: Demonstration at 3.5 kG," Radiology, vol. 161, pp. 527{531, 1986. [53] J. J. Cuppen and A. Van Est, \Reducing MR imaging time by one-sided reconstruction," in Topical Conf. Fast MRI Techiques, 1987. [54] D. C. Noll, D. G. Nishimura, and A. Macovski, \Homodyne detection in magnetic resonance imaging," IEEE Trans. Med. Imaging, vol. 10, pp. 154{163, 1991. [55] M. Bydder and M. D. Robson, \Partial Fourier partially parallel imaging," Magn. Reson. Med., vol. 53, pp. 1393{1401, 2005. [56] A. A. Samsonov, E. G. Kholmovski, D. L. Parker, and C. R. Johnson, \POCSENSE: POCS- based reconstruction for sensitivity encoded magnetic resonance imaging," Magn. Reson. Med., vol. 52, pp. 1397{1406, 2004. [57] F. Huang, W. Lin, and Y. Li, \Partial Fourier reconstruction through data tting and con- volution in k-space," Magn. Reson. Med., vol. 62, pp. 1261{1269, 2009. [58] M. Blaimer, M. Gutberlet, P. Kellman, F. A. Breuer, H. K ostler, and M. A. Griswold, \Virtual coil concept for improved parallel MRI employing conjugate symmetric signals," Magn. Reson. Med., vol. 61, pp. 93{102, 2009. Bibliography 95 [59] J. P. Haldar, \Low-rank modeling of localk-space neighborhoods (LORAKS) for constrained MRI," IEEE Trans. Med. Imaging, vol. 33, pp. 668{681, 2014. [60] J. P. Haldar and J. Zhuo, \P-LORAKS: Low-rank modeling of local k-space neighborhoods with parallel imaging data," Magn. Reson. Med., vol. 75, pp. 1499{1514, 2016. [61] T. H. Kim, K. Setsompop, and J. P. Haldar, \LORAKS makes better SENSE: Phase- constrained partial Fourier SENSE reconstruction without phase calibration," Magn. Reson. Med., vol. 77, pp. 2236{2249, 2017. [62] J. P. Haldar and K. Setsompop, \Linear predictability in magnetic resonance imaging recon- struction: Leveraging shift-invariant fourier structure for faster and better imaging," IEEE Signal Process. Mag., vol. 37, pp. 69{82, 2020. [63] F. Ong, J. Cheng, and M. Lustig, \General phase regularized reconstruction using phase cycling," Magn. Reson. Med., vol. 80, pp. 112{125, 2018. [64] Y. Liu and J. P. Haldar, \PALMNUT: An enhanced proximal alternating linearized min- imization algorithm with application to separate regularization of magnitude and phase," IEEE Trans. Comput. Imaging, vol. 7, pp. 530{518, 2021. [65] L. P. Panych, G. P. Zientara, and F. A. Jolesz, \MR image encoding by spatially selective RF excitation: an analysis using linear response models," Int. J. Imag. Syst. Tech., vol. 10, pp. 143{150, 1999. [66] A. A. Maudsley, \Multiple-line-scanning spin density imaging," J. Magn. Reson., vol. 41, pp. 112{126, 1980. [67] C. H. Oh, H. W. Park, and Z. H. Cho, \Line-integral projection reconstruction (LPR) with slice encoding techniques: Multislice regional imaging in NMR tomography," IEEE Trans. Med. Imaging, vol. MI-3, pp. 170{178, 1984. [68] S. P. Souza, J. Szumowski, C. L. Dumoulin, D. P. Plewes, and G. Glover, \SIMA: Simul- taneous multislice acquisition of MR images by Hadamard-encoded excitation," J. Comput. Assist. Tomogr., vol. 12, pp. 1025{1030, 1988. [69] F. A. Breuer, M. Blaimer, R. M. Heidemann, M. F. Mueller, M. A. Griswold, and P. M. Jakob, \Controlled aliasing in parallel imaging results in higher acceleration (CAIPIRINHA) for multi-slice imaging," Magn. Reson. Med., vol. 53, pp. 684{691, 2005. [70] J. B. Weaver, Y. Xu, D. M. Healy, and J. R. Driscoll, \Wavelet-encoded MR imaging," Magn. Reson. Med., vol. 24, pp. 275{287, 1992. [71] G. P. Zientara, L. P. Panych, and F. A. Jolesz, \Near-optimal spatial encoding for dynamically adaptive MRI: Mathematical principles and computational methods," Int. J. Imag. Syst. Tech., vol. 10, pp. 151{165, 1999. [72] J. P. Haldar, D. Hernando, and Z.-P. Liang, \Compressed-sensing MRI with random encod- ing," IEEE Trans. Med. Imaging, vol. 30, pp. 893{903, 2011. Bibliography 96 [73] T. H. Kim and J. P. Haldar, \SMS-LORAKS: Calibrationless simultaneous multislice MRI using low-rank matrix modeling," in Proc. IEEE Int. Symp. Biomed. Imaging, New York City, 2015, pp. 323{326. [74] J. P. Haldar and K. Setsompop, \Fast high-resolution diusion MRI using gSlider-SMS, interlaced subsampling, and SNR-enhancing joint reconstruction," in Proc. Int. Soc. Magn. Reson. Med., Honolulu, 2017, p. 175. [75] G. Ramos-Llorden, L. Ning, C. Liao, R. Mukhometzianov, O. Michailovich, K. Setsompop, and Y. Rathi, \High-delity, accelerated whole-brain submillimeter in vivo diusion MRI using gSlider-spherical ridgelets (gSlider-SR)," Magn. Reson. Med., vol. 84, pp. 1781{1795, 2020. [76] Y. Liu, C. Liao, K. Setsompop, and J. P. Haldar, \An evaluation of regularization strategies for subsampled single-shell diusion MRI," in Proc. IEEE Int. Symp. Biomed. Imaging, 2020, pp. 1437{1440. [77] K. Setsompop, Q. Fan, J. Stockmann, B. Bilgic, S. Huang, S. F. Cauley, A. Nummenmaa, F. Wang, Y. Rathi, T. Witzel, and L. L. Wald, \High-resolution in vivo diusion imaging of the human brain with generalized slice dithered enhanced resolution: Simultaneous multislice (gSlider-SMS)," Magn. Reson. Med., vol. 79, pp. 141{151, 2018. [78] J. P. Haldar, Y. Liu, C. Liao, Q. Fan, and K. Setsompop, \Fast submillimeter diusion MRI using gSlider-SMS and SNR-enhancing joint reconstruction," Magn. Reson. Med., vol. 84, pp. 762{776, 2020. [79] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, \SENSE: Sensitivity encoding for fast MRI," Magn. Reson. Med., vol. 42, pp. 952{962, 1999. [80] S. M. Kay, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory. Up- per Saddle River: Prentice Hall, 1993. [81] F. Pukelsheim, Optimal Design of Experiments. New York: John Wiley & Sons, 1993. [82] S. J. Reeves and Z. Zhe, \Sequential algorithms for observation selection," IEEEE Trans. Signal Process., vol. 47, pp. 123{132, 1999. [83] D. Xu, M. Jacob, and Z.-P. Liang, \Optimal sampling of k-space with Cartesian grids for parallel MR imaging," in Proc. Int. Soc. Magn. Reson. Med., 2005, p. 2450. [84] J. A. Jones, P. Hodgkinson, A. L. Barker, and P. J. Hore, \Optimal sampling strategies for the measurement of spin-spin relaxation times," J. Magn. Reson. B, vol. 113, pp. 25{34, 1996. [85] G. J. Marseille, M. Fuderer, R. de Beer, A. F. Mehlkopf, and D. van Ormondt, \Reduction of MRI scan time through nonuniform sampling and edge-distribution modeling," J. Magn. Reson. B, vol. 103, pp. 292{295, 1994. [86] S. Cavassila, S. Deval, C. Huegen, D. van Ormondt, and D. Graveron-Demilly, \Cram er-Rao bounds: an evaluation tool for quantitation," NMR Biomed., vol. 14, pp. 278{283, 2001. Bibliography 97 [87] A. R. Pineda, S. B. Reeder, Z. Wen, and N. J. Pelc, \Cram er-Rao bounds for three-point decomposition of water and fat," Magn. Reson. Med., vol. 54, pp. 625{635, 2005. [88] D. C. Alexander, \A general framework for experiment design in diusion MRI and its appli- cation in measuring direct tissue-microstructure features," Magn. Reson. Med., vol. 60, pp. 439{448, 2008. [89] J. P. Haldar, D. Hernando, and Z.-P. Liang, \Super-resolution reconstruction of MR image sequences with contrast modeling," in Proc. IEEE Int. Symp. Biomed. Imaging, 2009, pp. 266{269. [90] D. De Naeyer, Y. De Deene, W. P. Ceelen, P. Segers, and P. Verdonck, \Precision analysis of kinetic modelling estimates in dynamic contrast enhanced MRI," vol. 24, pp. 51{66, 2011. [91] B. Zhao, J. P. Haldar, C. Liao, D. Ma, Y. Jiang, M. A. Griswold, K. Setsompop, and L. L. Wald, \Optimal experiment design for magnetic resonance ngerprinting: Cramer-Rao bound meets spin dynamics," IEEE Trans. Med. Imaging, vol. 81, pp. 1620{1633, 2019. [92] C. D. Meyer, Matrix Analysis and Applied Linear Algebra. Philadelphia: SIAM, 2000. [93] M. Gendreau and J.-Y. Potvin, Eds., Handbook of Metaheuristics, 2nd ed. Springer, 2010. [94] Z. Ugray, L. Lasdon, J. Plummer, F. Glover, J. Kelly, and R. Mart , \Scatter search and local NLP solvers: A multistart framework for global optimization," INFORMS Journal on Computing, vol. 19, no. 3, pp. 328{340, 2007. [95] F. Glover, \A template for scatter search and path relinking," in Articial Evolution, J.-K. Hao, E. Lutton, E. Ronald, M. Schoenauer, and D. Snyers, Eds. Springer Berlin Heidelberg, 1998, pp. 1{51. [96] D. G. Nishimura, Principles of Magnetic Resonance Imaging. Stanford University, 1996. [97] S. Conolly, D. Nishimura, and A. Macovski, \Optimal control solutions to the magnetic resonance selective excitation problem," IEEE Transactions on Medical Imaging, vol. 5, no. 2, pp. 106{115, 1986. [98] J. B. Murdoch, A. H. Lent, and M. R. Kritzer, \Computer-optimized narrowband pulses for multislice imaging," Journal of Magnetic Resonance (1969), vol. 74, no. 2, pp. 226{263, 1987. [99] W. A. Grissom, D. Xu, A. B. Kerr, J. A. Fessler, and D. C. Noll, \Fast large-tip-angle multi- dimensional and parallel RF pulse design in MRI," IEEE Transactions on Medical Imaging, vol. 28, no. 10, pp. 1548{1559, 2009. [100] D. O. Brunner and K. P. Pruessmann, \Optimal design of multiple-channel RF pulses under strict power and SAR constraints," Magnetic Resonance in Medicine, vol. 63, no. 5, pp. 1280{1291, 2010. [101] Y. Liu and J. P. Haldar, \NAPALM: An algorithm for MRI reconstruction with separate magnitude and phase regularization," in Proc. Int. Soc. Magn. Reson. Med., 2019, p. 4764. Bibliography 98 [102] M. Cetin and W. C. Karl, \Feature-enhanced synthetic aperture radar image formation based on nonquadratic regularization," IEEE Trans. Image Process., vol. 10, pp. 623{631, 2001. [103] J. A. Fessler and D. C. Noll, \Iterative image reconstruction in MRI with separate magnitude and phase regularization," in Proc. IEEE Int. Symp. Biomed. Imaging, 2004, pp. 209{212. [104] A. Tuysuzoglu, J. M. Kracht, R. O. Cleveland, M. Cetin, and W. C. Karl, \Sparsity driven ultrasound imaging," J. Acoust. Soc. Am., vol. 131, pp. 1271{1281, 2012. [105] J. P. Haldar, Z. Wang, G. Popescu, and Z. P. Liang, \Deconvolved spatial light interference microscopy for live cell imaging," IEEE Trans. Biomed. Eng., vol. 58, pp. 2489{2497, 2011. [106] F. Zhao, D. C. Noll, J. F. Nielsen, and J. A. Fessler, \Separate magnitude and phase regu- larization via compressed sensing," IEEE Trans. Med. Imaging, vol. 31, pp. 1713{1723, 2012. [107] H. E. Guven, A. Gungor, and M. Cetin, \An augmented Lagrangian method for complex- valued compressed SAR imaging," IEEE Trans. Comput. Imaging, vol. 2, pp. 235{250, 2016. [108] M. V. W. Zibetti and A. R. De Pierro, \Improving compressive sensing in MRI with separate magnitude and phase priors," Multidim. Syst. Signal. Process., vol. 28, pp. 1109{1131, 2017. [109] F. Ong, J. Y. Cheng, and M. Lustig, \General phase regularized reconstruction using phase cycling," Magn. Reson. Med., vol. 80, pp. 112{125, 2018. [110] M. Moradikia, S. Samadi, and M. Cetin, \Joint SAR imaging and multi-feature decomposition from 2-D undersampled data via low-rankness plus sparsity priors," IEEE Trans. Comput. Imaging, vol. 5, pp. 1{16, 2018. [111] J. P. Haldar, Y. Liu, C. Liao, Q. Fan, and K. Setsompop, \Fast submillimeter diusion MRI using gSlider-SMS and SNR-enhancing joint reconstruction," Magn. Reson. Med., vol. 84, pp. 762{776, 2020. [112] Y. Liu and J. P. Haldar, \NAPALM: An algorithm for MRI reconstruction with separate magnitude and phase regularization," in Proc. Int. Soc. Magn. Reson. Med., 2019, p. 4764. [113] S. J. Wright, \Coordinate descent algorithms," Math. Program., vol. 151, pp. 3{34, 2015. [114] N. Parikh and S. Boyd, \Proximal algorithms," Found. Trends Optim., vol. 1, pp. 123{231, 2013. [115] A. Beck, First-Order Methods in Optimization. Philadelphia: SIAM, 2017. [116] Z.-P. Liang, \A model-based method for phase unwrapping," IEEE Trans. Med. Imaging, vol. 15, pp. 893{897, 1996. [117] M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, \Fast image recovery using variable splitting and constrained optimization," IEEE Trans. Image Process., vol. 19, pp. 2345{2356, 2010. [118] P. J. Huber, Robust Statistics. New York: John Wiley & Sons, Inc., 1981. Bibliography 99 [119] M. Nikolova and M. K. Ng, \Analysis of half-quadratic minimization methods for signal and image recovery," SIAM J. Sci. Comput., vol. 27, pp. 937{966, 2005. [120] M. J. Black and A. Rangarajan, \On the unication of line processes, outlier rejection, and robust statistics with applications in early vision," Int. J. Comput. Vis., vol. 19, pp. 57{91, 1996. [121] J. P. Haldar, V. J. Wedeen, M. Nezamzadeh, G. Dai, M. W. Weiner, N. Schu, and Z.- P. Liang, \Improved diusion imaging through SNR-enhancing joint reconstruction," Magn. Reson. Med., vol. 69, pp. 277{289, 2013. [122] L. I. Rudin, S. Osher, and E. Fatemi, \Nonlinear total variation based noise removal algo- rithms," Physica D, vol. 60, pp. 259{268, 1992. [123] J. A. Tropp, \Algorithms for simultaneous sparse approximation. Part II: Convex relaxation." Signal Proc., vol. 86, pp. 589{602, 2006. [124] J. P. Haldar and Z.-P. Liang, \Joint reconstruction of noisy high-resolution MR image se- quences," in Proc. IEEE Int. Symp. Biomed. Imaging, 2008, pp. 752{755. [125] P. Weiss, L. Blanc-Feraud, and G. Aubert, \Ecient schemes for total variation minimization under constraints in image processing," SIAM J. Sci. Comput., vol. 31, pp. 2047{2080, 2009. [126] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed. Baltimore: The Johns Hopkins University Press, 2013. [127] T. Pock and S. Sabach, \Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems," SIAM Journal on Imaging Sciences, vol. 9, no. 4, pp. 1756{1787, 2016. [128] M. A. Bernstein, D. M. Thomasson, and W. H. Perman, \Improved detectability in low signal- to-noise ratio magnetic resonance images by means of a phase-corrected real reconstruction," Math. Program., vol. 16, pp. 813{817, 1989. [129] G. C. McKinnon, X. J. Zhou, and N. E. Leeds, \Phase corrected complex averaging for diusion weighted spine imaging," in Proc. Int. Soc. Magn. Reson. Med., 2000, p. 802. [130] C. Liu, R. Bammer, D. h. Kim, and M. E. Moseley, \Self-navigated interleaved spiral (SNAILS): Application to high-resolution diusion tensor imaging," Magn. Reson. Med., vol. 52, pp. 1388{1396, 2004. [131] N.-K. Chen, A. Guidon, H.-C. Chang, and A. W. Song, \A robust multi-shot scan strat- egy for high-resolution diusion weighted MRI enabled by multiplexed sensitivity-encoding (MUSE)," NeuroImage, vol. 72, pp. 41 { 47, 2013. [132] C. Eichner, S. F. Cauley, J. Cohen-Adad, H. E. Moller, R. Turner, K. Setsompop, and L. L. Wald, \Real diusion-weighted MRI enabling true signal averaging and increased diusion contrast," NeuroImage, vol. 122, pp. 373{384, 2015. [133] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C. Cambridge: Cambridge University Press, 1992. Bibliography 100 [134] K. P. Pruessmann, M. Weiger, P. Bornert, and P. Boesiger, \Advances in sensitivity encoding with arbitrary k-space trajectories," Magn. Reson. Med., vol. 46, pp. 638{651, 2001. [135] M. Uecker, P. Lai, M. J. Murphy, P. Virtue, M. Elad, J. M. Pauly, S. S. Vasanawala, and M. Lustig, \ESPIRiT { an eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA," Magn. Reson. Med., vol. 71, pp. 990{1001, 2014. [136] Z. Wang and A. C. Bovik, \Mean squared error: Love it or leave it? a new look at signal delity measures," IEEE Signal Process. Mag., vol. 26, pp. 98{117, 2009. [137] Y. Liu, C. Liao, D. Kim, K. Setsompop, and J. P. Haldar, \Estimating multicomponent 2D relaxation spectra with a ViSTa-MR ngerprinting acquisition," in Proc. Int. Soc. Magn. Reson. Med., London, 2022, p. 4389. [138] D. McGivney, A. Deshmane, Y. Jiang, D. Ma, C. Badve, A. Sloan, V. Gulani, and M. Gris- wold, \Bayesian estimation of multicomponent relaxation parameters in magnetic resonance ngerprinting," Magn. Reson. Med., vol. 80, pp. 159{170, 2018. [139] A. Deshmane, D. F. McGivney, D. Ma, Y. Jiang, C. Badve, V. Gulani, N. Seiberlich, and M. A. Griswold, \Partial volume mapping using magnetic resonance ngerprinting," NMR Biomed., vol. 32, p. e4082, 2019. [140] S. Tang, C. Fernandez-Granda, S. Lannuzel, B. Bernstein, R. Lattanzi, M. Cloos, F. Knoll, and J. Asslander, \Multicompartment magnetic resonance ngerprinting," Inverse Probl., vol. 34, p. 094005, 2018. [141] M. Nagtegaal, P. Koken, T. Amthor, and M. Doneva, \Fast multi-component analysis using a joint sparsity constraint for MR ngerprinting," Magn. Reson. Med., vol. 83, pp. 521{534, 2020. [142] P. J. Slator, M. Palombo, K. Miller, C. F. Westin, F. Laun, D. Kim, J. P. Haldar, D. Ben- jamini, G. Lemberskiy, J. P. de Almeida Martins, and J. Hutter, \Combined diusion- relaxometry microstructure imaging: Current status and future prospects," Magn. Reson. Med., vol. 86, pp. 2987{3011, 2021. [143] D. Benjamini and P. J. Basser, \Multidimensional correlation MRI," NMR Biomed., vol. 33, p. e4226, 2020. [144] F. O'Sullivan, \Imaging radiotracer model parameters in PET: A mixture analysis approach," IEEE Trans. Med. Imaging, vol. 12, pp. 399{412, 1993. [145] M. Ingrisch and S. Sourbron, \Tracer-kinetic modeling of dynamic contrast-enhanced MRI and CT: A primer," J. Pharmacokinet. Pharmacodyn., vol. 40, pp. 281{300, 2013. [146] M. D. Does, \Inferring brain tissue composition and microstructure via MR relaxometry," NeuroImage, vol. 182, pp. 136{148, 2018. [147] D. C. Alexander, T. B. Dyrby, M. Nilsson, and H. Zhang, \Imaging brain microstructure with diusion MRI: Practicality and applications," NMR Biomed., vol. 32, p. e3841, 2019. Bibliography 101 [148] A. A. Istratov and O. F. Vyvenko, \Exponential analysis in physical phenomena," Rev. Sci. Instrum., vol. 70, pp. 1233{1257, 1999. [149] H. Celik, M. Bouhrara, D. A. Reiter, K. W. Fishbein, and R. G. Spencer, \Stabilization of the inverse Laplace transform of multiexponential decay through introduction of a second dimension," J. Magn. Reson., vol. 236, pp. 134{139, 2013. [150] R. M. Kroeker and R. M. Henkelman, \Analysis of biological NMR relaxation data with continuous distributions of relaxation times," J. Magn. Reson., vol. 69, pp. 218{235, 1986. [151] K. P. Whittall and A. L. MacKay, \Quantitative interpretation of NMR relaxation data," J. Magn. Reson., vol. 84, pp. 134{152, 1989. [152] L. Venkataramanan, Y.-Q. Song, and M. D. Hurlimann, \Solving Fredholm integrals of the rst kind with tensor product structure in 2 and 2.5 dimensions," IEEEE Trans. Signal Process., vol. 50, pp. 1017{1026, 2002. [153] R. Bai, A. Cloninger, W. Czaja, and P. J. Basser, \Ecient 2D MRI relaxometry using compressed sensing," J. Magn. Reson., vol. 255, pp. 88{99, 2015. [154] D. Benjamini and P. J. Basser, \Use of marginal distributions constrained optimization (MADCO) for accelerated 2D MRI relaxometry and diusometry," J. Magn. Reson., vol. 271, pp. 40{45, 2016. [155] Y. Lin, J. P. Haldar, Q. Li, P. S. Conti, and R. M. Leahy, \Sparsity constrained mixture modeling for the estimation of kinetic parameters in dynamic PET," IEEE Trans. Med. Imaging, vol. 33, pp. 173{185, 2014. [156] D. Hwang and Y. P. Du, \Improved myelin water quantication using spatially regularized non-negative least squares algorithm," J. Magn. Reson. Imaging, vol. 30, pp. 203{208, 2009. [157] D. Kumar, T. D. Nguyen, S. A. Gauthier, and A. Raj, \Bayesian algorithm using spatial priors for multiexponential T2 relaxometry from multiecho spin echo MRI," Magn. Reson. Med., vol. 68, pp. 1536{1543, 2012. [158] C. Labadie, J.-H. Lee, W. D. Rooney, S. Jarchow, M. Aubert-Frecon, C. S. Springer, Jr, and H. E. Moller, \Myelin water mapping by spatially regularized longitudinal relaxographic imaging at high magnetic elds," Magn. Reson. Med., vol. 71, pp. 375{387, 2014. [159] D. Kumar, H. Hariharan, T. D. Faizy, P. Borchert, S. Siemonsen, J. Fiehler, R. Reddy, and J. Sedlacik, \Using 3D spatial correlations to improve the noise robustness of multi component analysis of 3D multi echo quantitative T2 relaxometry data," NeuroImage, vol. 178, pp. 583{601, 2018. [160] M. Zimmermann, A.-M. Oros-Peusquens, E. Iordanishvili, S. Shin, S. D. Yun, Z. Abbas, and N. J. Shah, \Multi-exponential relaxometry using ` 1 -regularized iterative NNLS (MERLIN) with application to myelin water fraction imaging," IEEE Trans. Med. Imaging, vol. 38, pp. 2676{2686, 2019. Bibliography 102 [161] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, \Distributed optimization and statistical learning via the alternating direction method of multipliers," Found. Trends Mach. Learn., vol. 3, pp. 1{122, 2011. [162] B. He, L.-Z. Liao, D. Han, and H. Yang, \A new inexact alternating directions method for monotone variational inequalities," Math. Program., vol. 92, pp. 103{118, 2002. [163] X. Zhang, M. Burger, and S. Osher, \A unied primal-dual algorithm framework based on Bregman iteration," J. Sci. Comput., vol. 46, pp. 20{46, 2011. [164] B. He, F. Ma, and X. Yuan, \Optimally linearizing the alternating direction method of multipliers for convex programming," Comput. Optim. Appl., vol. 75, pp. 361{388, 2020. [165] M. Fazel, T. K. Pong, D. Sun, and P. Tseng, \Hankel matrix rank minimization with appli- cations to system identication and realization," SIAM J. Matrix Anal. Appl., vol. 34, pp. 946{977, 2013. [166] M. Tao, \Convergence study of indenite proximal ADMM with a relaxation factor," Comput. Optim. Appl., vol. 77, pp. 91{123, 2020. [167] W. Deng and W. Yin, \On the global and linear convergence of the generalized alternating direction method of multipliers," J. Sci. Comput., vol. 66, pp. 889{916, 2016. [168] D. Kim, B. Zhao, L. L. Wald, and J. P. Haldar, \Multidimensional T1 relaxation-T2 relaxation correlation spectroscopic imaging witha magnetic resonance ngerprinting acquisition," in Proc. Int. Soc. Magn. Reson. Med., Montreal, 2019, p. 4991. [169] S. W. Provencher, \A constrained regularization method for inverting data represented by linear algebraic or integral equations," Comput. Phys. Commun., vol. 27, pp. 213{227, 1982. [170] D. Varadarajan and J. P. Haldar, \A majorize-minimize framework for Rician and non-central chi MR images," IEEE Trans. Med. Imaging, vol. 34, pp. 2191{2202, 2015. [171] S. Boyd, Convex Optimization. Cambridge: Cambridge University Press, 2004. [172] M. Slawski and M. Hein, \Non-negative least squares for high-dimensional linear models: Consistency and sparse recovery without regularization," Electronic Journal of Statistics, vol. 7, pp. 3004 { 3056, 2013. [173] C. L. Lawson and R. J. Hanson, Solving Least Squares Problems. Philadelphia: SIAM, 1995. [174] M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, \An augmented Langrangian approach to the constrained optimization formulation of imaging inverse problems," IEEE Trans. Image Process., vol. 20, pp. 681{695, 2011. [175] P. L. Combettes and J.-C. Pesquet, \Proximal splitting methods in signal processing," in Fixed-Point Algorithms for Inverse Problems in Science and Engineering, H. H. Bauschke, R. S. Burachik, P. L. Combettes, V. Elser, D. R. Luke, and H. Wolkowicz, Eds. New York: Springer New York, 2011, pp. 185{212. Appendix B 103 [176] N. Parikh and S. Boyd, \Proximal algorithms," Found. Trends Optim., vol. 1, pp. 127{239, 2014. [177] M. Yang, D. Ma, Y. Jiang, J. Hamilton, N. Seiberlich, M. A. Griswold, and D. McGivney, \Low rank approximation methods for MR ngerprinting with large scale dictionaries," Magn. Reson. Med., vol. 79, pp. 2392{2400, 2018. [178] E. Ghadimi, A. Teixeira, I. Shames, and M. Johansson, \Optimal parameter selection for the alternating direction method of multipliers (admm): Quadratic problems," IEEE Trans. Automat. Contr., vol. 60, pp. 644{658, 2015. [179] B. S. He, H. Yang, and S. L. Wang, \Alternating direction method with self-adaptive penalty parameters for monotone variational inequalities," J. Optim. Theory Appl., vol. 106, pp. 337{356, 2000. [180] Z. Xu, M. Figueiredo, and T. Goldstein, \Adaptive ADMM with Spectral Penalty Parameter Selection," in Proc. Int. Conf. Arti. Intel. Stat., A. Singh and J. Zhu, Eds., vol. 54. PMLR, 2017, pp. 718{727. [181] C. Liao, X. Cao, T. Gong, Z. Wu, Z. Zhou, H. He, J. Zhong, and K. Setsompop, \High- resolution myelin-water fraction (MWF) and T1/T2/proton-density mapping using 3D ViSTa-MR ngerprinting with subspace reconstruction," in Proc. Int. Soc. Magn. Reson. Med., Virtual, 2021, p. 1545. [182] M. Weigel, \Extended phase graphs: Dephasing, RF pulses, and echoes - pure and simple," J. Magn. Reson. Imaging, vol. 41, pp. 266{295, 2015.

Abstract (if available)

Abstract MRI is a powerful tool in both biomedical research and clinical diagnosis. However, the trade-offs among data acquisition time, image signal-to-noise-ratio, and spatial resolution often hinders the applications of MRI in its full potentials. In addition to imaging hardware improvement and advanced pulse sequence development, constrained MRI methods have greatly expanded our toolboxes in dealing with these trade-offs. In constrained MRI, we utilize prior information about the characteristics of the underlying MRI images to perform data acquisition, image reconstruction, and image analysis tasks more efficiently. This approach generally requires the use of mathematical optimization techniques, although the optimization problems are often challenging to solve efficiently due to their large-scale and non-trivial structure.
In this dissertation, three novel contributions in mathematical optimization for constrained MRI were made. First, we proposed a novel method that utilizes phase constraints to accelerate MRI data acquisition based on non-Fourier radiofrequency encoding. While phase constraints are used classically in MRI, we believe that this is the first time that phase constraints are being applied to enable acceleration along a non-Fourier encoded spatial dimension. We make the novel observation that phase constraints can indeed be successfully used to reduce the number of required non-Fourier encodings, although this requires careful design of the non-Fourier encoding scheme. Results are presented in the context of gSlider, an acquisition method designed for highly-efficient high-resolution diffusion MRI. Second, we proposed a new algorithm for the separate regularization of magnitude and phase in MRI reconstruction problems. Our approach is based on a novel application of the proximal alternating linearized minimization algorithm (PALM), and incorporates additional novel features (i.e., Nesterov's momentum and independent selection of the step sizes for each coordinate) to increase convergence speed. Depending on the application, our proposed algorithm can be hundreds of times faster than existing algorithms for this problem. Finally, we proposed another novel algorithm for spatial-spectral partial volume compartment mapping with applications to multicomponent diffusion and relaxation MRI. Our proposed algorithm is based on a novel application of the linearized alternating directions method of multipliers (LADMM) approach that takes advantage of the special structure of the inverse problem, and depending on the dataset, can achieve up to 5-fold acceleration compared to previous algorithms for this problem.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

New theory and methods for accelerated MRI reconstruction

PDF

Shift-invariant autoregressive reconstruction for MRI

PDF

Novel theoretical characterization and optimization of experimental efficiency for diffusion MRI

PDF

Landscape analysis and algorithms for large scale non-convex optimization

PDF

Fast flexible dynamic three-dimensional magnetic resonance imaging

PDF

Correction, coregistration and connectivity analysis of multi-contrast brain MRI

PDF

Architecture design and algorithmic optimizations for accelerating graph analytics on FPGA

PDF

Multi-constrained inversion algorithms for microwave imaging

PDF

High-dimensional magnetic resonance imaging of microstructure

PDF

On the theory and applications of structured bilinear inverse problems to sparse blind deconvolution, active target localization, and delay-Doppler estimation

PDF

Methods for improving reliability and consistency in diffusion MRI analysis

PDF

Improved brain dynamic contrast enhanced MRI using model-based reconstruction

PDF

Solution of inverse scattering problems via hybrid global and local optimization

PDF

Novel optimization tools for structured signals recovery: channels estimation and compressible signal recovery

PDF

Artificial Decision Intelligence: integrating deep learning and combinatorial optimization

PDF

Provable reinforcement learning for constrained and multi-agent control systems

PDF

Thermal analysis and multiobjective optimization for three dimensional integrated circuits

PDF

Learning and control for wireless networks via graph signal processing

PDF

Seeing sleep: real-time MRI methods for the evaluation of sleep apnea

PDF

Train routing and timetabling algorithms for general networks

Asset Metadata

Creator Liu, Yunsong (filename)

Core Title Optimization methods and algorithms for constrained magnetic resonance imaging

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Degree Conferral Date 2022-12

Publication Date 11/07/2022

Defense Date 09/16/2022

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag inverse problem,MRI,OAI-PMH Harvest,optimization,signal processing

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Haldar, Justin P. (committee chair), Leahy, Richard M. (committee member), Wood, John C. (committee member)

Creator Email songyunliu3526@gmail.com,yunsongl@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC112305184

Unique identifier UC112305184

Identifier etd-LiuYunsong-11307.pdf (filename)

Legacy Identifier etd-LiuYunsong-11307

Document Type Dissertation

Format theses (aat)

Rights Liu, Yunsong

Internet Media Type application/pdf

Type texts

Source 20221107-usctheses-batch-990 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu