1 Introduction
The main disadvantages of Magnetic Resonance Imaging (MRI) are its long scan times and, in consequence, its sensitivity to motion. Exploiting the complementary information from multiple receive coils, parallel imaging is able to recover images from undersampled kspace data and to accelerate the measurement [1, 2, 3, 4, 5, 6, 7]. Because parallel magnetic resonance imaging can be used to accelerate basically any imaging sequence it has many important applications. Parallel imaging brought a fundamental shift in image reconstruction: Image reconstruction changed from a simple direct Fourier transform to the solution of an illconditioned inverse problem. This work gives an overview of image reconstruction from the perspective of inverse problems. After introducing basic concepts such as regularization, discretization, and iterative reconstruction, advanced topics are discussed including algorithms for autocalibration, the connection to approximation theory, and the combination with compressed sensing.
2 Parallel Imaging as Inverse Problem
2.1 Forward Model
The signal from a receivecoil array with channels is given by [8]:
(1) 
The complexvalued magnetization image represents the state of the transverse magnetization of the excited spins in the fieldofview (FOV) at the time of image acquisition. In MRI, typically a volume () or a thin slice () of proton spins is excited using a resonant radiofrequency pulse. The image (or volume) is modulated by the complexvalued sensitivities of all receive coils (Fig. 1). The kspace signals are then given by the Fourier transform of the coil images sampled at discrete time points along a given kspace trajectory . Equation 1 neglects relaxation and offresonance effects during the acquisition, which is possible if all samples at time points are acquired in a small window around the echo time (TE) after excitation. Relaxation and offresonance effects from the excitation of the spins until the echo time are incorporated into the image and define the image contrast. Because in most cases it is not possible to acquire all data in this short acquisition window, samples have to be acquired in a repeated series of identical experiments, which always restore the magnetization image to exactly the same state.^{1}^{1}1Singleshot EchoPlanar Imaging (EPI) and spiral imaging sequences acquire all data in a single acquisition. These sequences are fast, but image quality is compromised by blurring, distortions, and phase cancellation, due to relaxation and offresonance effects. Parallel imaging can be used to shorten the acquisition window and reduce these artifacts. This requirement to repeat the basic experiment many times is the reason for the long scans times in MRI. The goal of parallel imaging is to reduce the amount of data required to reconstruct images by optimally exploiting the complementary information from multiple receive coils. Although there are fundamental limits to the encoding power of the receive sensitivities, it has the potential to accelerate MRI by a factor of about four in each spatial dimension [9].
2.2 Image Reconstruction
If the coil sensitivities are known (e.g. from a prescan), image reconstruction for parallel imaging can be formulated as a linear inverse problem with discrete data [10]. Mathematically, the forward problem is given by an operator which maps the magnetization image of excited spins in a FOV to the sample values . This operator can be thought of as the composition of a physical model for hypothetical multidimensional continuous kspace signals and a sampling operator (see Fig. 3). The operator is given by
(2) 
with the encoding functions defined as and a scalar product defined on (antilinear in the second argument). The sampling operator evaluates the ideal kspace signals at the sample locations in kspace. It is assumed that the sampling process corresponds to the pointevaluation of ideal kspace functions, i.e. . The sample values are corrupted by additive white Gaussian noise. Although the noise is typically correlated between receive channels, this correlation can be removed with a whitening step.
A variational solution to the inverse problem can be defined as the minimizer of a functional, i.e.
(3) 
The functional is composed of a leastsquares data fidelity term (which alternatively may also include weighting or use a robust norm [11, 12]) and an additional regularization term . Discretized versions of this minimization problem are the basis of SMASH and SENSE parallel imaging methods [5, 6, 7, 13, 14]. For parallel MRI, this formulation has two advantages: First, arbitrary Cartesian or nonCartesian sampling schemes can be used [14]. Second, the regularization term can be used to introduce prior knowledge about the solution.^{2}^{2}2A third advantage  which conceptually goes beyond parallel imaging  is the possibility to extend the forward model to include further physical effects in modelbased reconstruction, e.g. field maps [15, 16], motioninduced phase maps [17, 18], motion [19, 20], relaxation maps [21, 22, 23], or diffusion models [24], etc. Figure 4 shows images of a human brain recovered from undersampled data by numerical optimization of Equation 3.
2.3 Regularization
Illconditioning causes noise amplification during image reconstruction, which initially limited the application of parallel imaging to only moderate acceleration. This limitation can be overcome by incorporating prior knowledge about the image using regularization methods [25, 13, 26]. In the simplest case, regularization may consist of a basic quadratic penalty in the framework of a linear reconstruction, or make use of much more sophisticated techniques which exploit the structure of images but demand a nonlinear reconstruction. For a leastsquares problem with quadratic regularization, i.e. with a positive definite operator , the solution of Eq. 3 is explicitly given by the formula
(4) 
In the limit this solution is called the best approximate and is given by the MoorePenrose pseudoinverse
. It has a statistical interpretation  assuming white Gaussian noise  as the best unbiased estimate for the image. Regularization can be interpreted as prior knowledge and the optimizer as a maximum a posteriori (MAP) estimate of the image. Although regularization leads to a fundamental tradeoff between bias and noise  which has to be chosen carefully for optimal image quality  it makes the use of higher acceleration possible. An optimal estimate in terms of mean squared error can only be obtained with regularization.
For optimal results, the prior knowledge should include as much specific knowledge about the image as possible. For example, regularization can exploit smoothness in the time domain [27], or exploit that changes relative to a fixed reference image can be assumed to be small. The later is used successfully for realtime MRI [28] or dynamic contrast enhanced (DCE) MRI [29]. While regularization is simple to implement and already a clear improvement compared to unregularized parallel imaging, much better results can be be obtained when using more advanced techniques such as wavelet regularization, i.e. , total variation, or other edgepreserving penalties [30, 31, 32, 33].
2.4 Discretization
Numerical reconstruction methods make use of discretization, i.e. the unknown image is expanded into a sum of basis functions:
(5) 
Although the choice of this basis has subtle implications for results and interpretation, this topic has not drawn much attention.^{3}^{3}3Discussions in terms of “ideal voxel functions” (or “target voxel shapes”) can be found in earlier works [7, 13]. Most imagedomain formulations based on SENSE use a grid of Dirac pulses to represent the image, because multiplication with the coil sensitivities in the forward model is then simply a pointwise multiplication. In contrast, kspace methods such as SMASH use a finite Fourier basis.
Discretization has a regularizing effect, i.e. the discretized problem might have better condition than the continuous problem. In parallel imaging, this effect can often be seen in the area outside of the sampled kspace region. Extrapolation to these area causes high noise amplification [34]
. A discretization scheme which excludes these degrees of freedom will be less affected by noise. On the other hand, a small basis leads to discretization errors because the solution can not be represented accurately. Both problems can be avoided by using fine discretization with a large number of basis function and explicit regularization to control noise amplification
[35].It should be noted that for parallel imaging the ideal continuous solution of Eq. 3 can usually be computed almost perfectly [36]
. Coil sensitivities are very smooth and can be approximated with a small number of Fourier coefficients (on an oversampled FOV). The forward operator can then be understood as a convolution of the Fourier series of the image with a short filter. Because the acquired kspace data consists of a finite number of samples, also only a finite number of loworder Fourier coefficients from the infinite number of coefficients in the Fourier series of the image actually appear in the result of this convolution. For quadratic regularization, a minimumnorm solution is obtained when the infinite number of remaining higherorder coefficients are set to zero. An implementation of the forward operator requires an aperiodic convolution which can be implemented efficiently using a fast Fourier transform (FFT) algorithm. In practice this differs from a conventional SENSE implementation only by using zero padding and in the exact interpretation of the recovered coefficients.
For nonquadratic regularization, discretization errors may also arise in the implementation the regularization terms. In general, oversampling can be used to reduce these errors. The combination of nonquadratic regularization and oversampling can also avoid artifacts caused by truncation of the signal in the Fourier domain (Gibbs ringing) [37]. Finally, an important aspect related to discretization is a common error called an “inverse crime” [38]: When testing a reconstruction algorithms with simulated data, computing this data using the same discretization scheme as used for the reconstruction can result in highly misleading results. One possibility to avoid this error is the use of analytic phantoms [39].
2.5 Numerical Optimization
For regular sampling schemes and quadratic regularization, a solution can be computed directly with matrix inversion, because the equations decouple into small systems [7]. Although very efficient, this approach is not very flexible. Matrixfree iterative methods can be used instead to efficiently compute the solution for arbitrary sampling schemes [14]
. Matrixfree methods are build from procedural implementations of the matrixvector products
and (or ). For Cartesian sampling, these operations can be implemented using pointwise multiplications and FFT algorithms. For nonCartesian sampling, efficient nonuniform fast Fourier transform (nuFFT) algorithms have been developed to estimate samples at arbitrary kspace locations [40, 41, 42]. Even more efficient algorithms can be designed when considering the combined operator which appears in the gradient of the leastsquares data fidelity. For example, the effect of sampling in the Fourier domain can be computed exactly as a convolution with a truncated pointspread function with the use of two zeropadded FFTs [43]. Overlapadd and overlapsave convolution algorithms can be used to exploit the compact representation of the coil sensitivities in the Fourier domain [44].For quadratic regularization, an efficient iterative algorithm is the conjugate gradient method applied to the normal equations [14]:
(6) 
It should be noted that the use of a density compensation as known from the direct gridding algorithm is neither required nor recommended.^{4}^{4}4In (noniterative) gridding, the density compensation is a diagonal matrix which approximates the inverse of . Combined with the adjoint it yields an approximation of the pseudoinverse, i.e. . Including a density compensation into an iterative optimization method produces solutions different from the optimal leastsquares solution [14]. I.e., naively using instead of as is sometimes suggested to improve the condition yields a different optimization problem:
For regularization, the simplest (and slowest) reconstruction algorithm is iterative softthresholding [45]:
(7)  
(8) 
The first equation is a gradient descent step and the second update uses softthresholding in a transform basis , e.g. a discrete wavelet transform. This scheme converges slowly, but can be accelerated with the addition of a ravine step as in FISTA [46, 47]. Especially when using multiple convex penalties , a very flexible approach is an extension of the Alternating Direction Method of Multipliers (ADMM) [48, 49, 50] that can solve optimization problems of the form
(9) 
This approach is very flexible and has many advantages from an implementation point of view, because it splits the optimization into independent subproblems. Many different kinds of regularization terms can easily be integrated if respective proximal operators of the form
(10) 
are available in a computationally efficient form. For example, the proximal operator for the data fidelity term is the regularized leastsquare inverse which can be computed efficiently with the methods of conjugate gradients. The proximal operator for regularization can be evaluated simply using softthesholding. Efficient implementations of many advanced algorithms which make use of parallel programming can be found Berkeley Advanced Reconstruction Toolbox (BART) [51].
3 Autocalibration
To obtain optimal results in parallel MRI, accurate and uptodate information about the sensitivities of all receive coils is required. While approximate coil sensitivities can be computed from the geometry of the receive coils using the Bios Savart law, exact sensitivities depend on the loading of the coils and need to be determined with high accuracy during the actual measurement. A prescan can provide accurate calibration information, but this requires that experimental conditions stay exactly the same for the duration of the whole examination. Because this is not always guaranteed, autocalibration methods have been developed which perform calibration using a small amount of additional data acquired during each individual scan [52, 53, 54, 55]. Because this reduces overall acceleration, optimal calibration from a minimum amount of data is desired. Two advanced techniques are described in the following: Joint estimation techniques simultaneously estimate image content and coil sensitivities from all data, which minimizes the amount of additional calibration data required. Subspace methods do not directly estimate sensitivities, but learn a signal subspace from calibration data. These algorithms can adapt to experimental conditions that violate the sensitivitybased signal model formulated in Equation 1. For this reason, they are more robust to certain kinds of errors.
3.1 Nonlinear Inverse Reconstruction
Starting with the signal equation (Eq. 1), but now considering both image and coil sensitivities as unknowns, one obtains a nonlinear inverse problem related to blind multichannel deconvolution. Modelling the coil sensitivities as smooth functions in a Sobolev space , the nonlinear version of the forward operator can be written as:
(11) 
Many autocalibrating parallel imaging methods reduce the reconstruction problem to a linear problem by first estimating the sensitivities from a subset of the data and then solving for the image using these fixed estimates using a conventional linear reconstruction. Because this is suboptimal, improved algorithms have been developed, which solve the nonlinear inverse problem [56, 57, 58].
In nonlinear inversion [58], a regularized solution is defined as the solution of the minimization problem
(12) 
Here, a smoothness penalty restricts the solution to smooth coil sensitivities. The Iterative Regularized GaussNewton Method (IRGNM) [59] is used to iteratively update an estimate of the solution based on a linearization of the original problem:
(13) 
Here, is the Frechét derivative of at the current estimate and the regularization parameters are reduced in each iteration step. The smoothness penalty can be chosen as (with some constants ). This penalty can be transformed into a norm by expressing the sensitivities using Fourier coefficients rescaled with a positive definite diagonal matrix, which avoids bad conditioning of the of the reconstruction problem. For quadratic regularization of the image, i.e. , the algorithm then has the explicit update rule
(14) 
Nonlinear reconstruction methods can be applied to nonCartesian sampling [60, 61, 62] and extended to include nonlinear penalties [32, 63].
One limitation of nonlinear methods is that they may need an initial guess close to solution to converge to the correct global minimum. While it is usually sufficient to set the image to a constant value and the coil sensitivities to zero, in some cases a guess closer to the true solution is required. In this case, any direct estimation method can be used to estimate a set of approximate coil sensitivities, which can then be used to initialize the nonlinear method.
3.2 Calibration Matrix
The calibration matrix is a fundamental tool which can be used to formulate many autocalibration methods. Reconstruction kernels in GRAPPA [53] and SPIRiT [64] are nullspace vectors of this matrix [65]. The calibration matrix is a multidimensional multichannel Casorati matrix constructed from fullysampled patches in a calibration area in the center of kspace (see Fig. 6). It is related to the trajectory matrix from singular spectrum analysis (SSA) [66], and also to the lag crosscovariance matrix, which can be estimated as (with a normalization constant). Because coil sensitivities are very smooth, multichannel signals have correlations in small local kspace patches. This implies that the calibration matrix (and the lag crosscovariance matrix) are lowrank, i.e. have a small signal space and large null space (Fig. 7). If the calibration matrix is constructed from an incomplete kspace with missing samples, structured lowrank matrix completion can be used to recover a completed matrix, which is the basis of a calibrationless parallel imaging technique known as SAKE [67].
3.3 ESPIRiT
Coilbycoil reconstruction was originally proposed because combination of all channels in SMASHbased parallel imaging sometimes caused phase cancellation [68, 53]. In combination with autocalibration coilbycoil reconstruction has a very advantageous side effect: The reconstruction becomes robust against certain kinds of inconsistencies  in particular reconstruction in a tight FOV [69]. The fundamental reason is that the coilbycoil reconstruction operator does not enforce the strict signal model of sensitivitybased reconstruction schemes formulated in Eq. 1, but represents a convex relaxation of this model.
ESPIRiT is a new reconstruction algorithm which exploits this. Because of shiftinvariance the nullspace condition should be true for all patches in an ideal multichannel kspace . Let be an operator which extracts a patch around a given kspace position , a leastsquares version of this condition is then given by
(15) 
This can be further transformed to a convolutiontype coilbycoil operator , which reproduces ideal kspace signals:
(16)  
(17)  
(18) 
Here, is the size of a single patch. Transforming into the image domain yields an operator which operates pointwise. Because reproduces ideal kspace signals, the imagedomain version reproduces the vector of coil images at each point , i.e. . In other words, everywhere where the image is nonzero the vector of sensitivities is a pointwise eigenvector to the eigenvalue of the operator . The eigenvector and eigenvalue maps from a pointwise eigendecomposition are shown in Figure 7. Together, these steps form a computational method to extract accurate coil sensitivities from the nullspace of the calibration matrix.
If the kspace is corrupted and does not fit the ideal model, multiple sets of sensitivities can appear in affected image regions as multiple eigenvectors to eigenvalue one. An extended forward model can take this additional information into account. Respective methods offer robustness to certain kinds of errors similar to autocalibrating coilbycoil methods such as GRAPPA [65].
4 Sampling and Reconstruction in kSpace
While the formulation of parallel imaging as an inverse problem is a powerful conceptual framework, additional theoretical tools are required to understand and evaluate different sampling schemes in kspace. For this purpose, a formulation of parallel imaging as approximation in a Reproducing Kernel Hilbert Space (RKHS) has recently been developed [34]. A RKHS is a Hilbert space of functions with continuous (bounded) pointevaluation functionals. This condition guarantees that sampling is compatible with the norm of the Hilbert space. This formulation of parallel imaging yields a unified framework to formulate and analyze sampling and reconstruction in kspace.
By identifying the norm of the signals in kspace with the norm of the corresponding images in , it can be shown that the ideal kspace signals are an RKHS with the matrixvalued kernel
(19) 
This kernel captures all local correlations in kspace induced by the sensitivities and exploited in parallel imaging algorithms for recovery of missing kspace samples. Given this kernel, a standard formula from approximation theory can be applied to obtain interpolation coefficients
for all channels and known kspace samples for interpolation to arbitrary kspace positions:(20) 
With these interpolation coefficients, unknown values in kspace can then be recovered from the acquired samples with the interpolation formula
(21) 
When no regularization is used in the computation of coefficients the recovered ideal kspace corresponds to the best approximate solution defined before. The interpolation formulas used in GRAPPA and SPIRiT and similar methods are local variants of this formula with an empirical estimate of the ideal kernel (Eq. 19). In addition to this interpolation formula, the link to approximation theory yields new insights into sampling in kspace. In particular, a pointwise error bound in kspace can be derived [70]:
(22) 
The power function is computed from the kernel and the interpolation functions and depends only on the coil sensitivities and the sample locations. It can be used to analyze the properties of different sampling schemes for parallel imaging independent from any actual imaging data. Figure 8 shows the Power function for two different sampling schemes computed for a particular set of coil sensitivities.
5 Compressed Sensing Parallel Imaging
Compressed sensing is based on the idea that randomized undersampling schemes produce incoherent noiselike artifacts in a transform domain which can then be suppressed using denoising to iteratively recover the original signal [71, 72]. It exploits the compressibility of the image information, i.e. a sparse representation in a transform domain, to make the sparse signal coefficients stand out from the incoherent noise. Nonlinear regularization terms can then be used to suppress the incoherent artifacts and recover a sparse representation of the image from the undersampled data. Because MRI acquires data in the Fourier domain and is flexible enough to use almost arbitrary sampling schemes, this idea can be applied directly [31, 73].
Parallel imaging can be synergistically combined with compressed sensing [31, 33, 64, 74]. This combination leads to exactly the same optimization problems already considered for parallel imaging alone, but requires incoherent sampling schemes suitable for compressed sensing. The most important schemes in practical use are variabledensity Poissondisc sampling and radial trajectories. Poissondisc sampling guarantees that samples are not too close together. This would waste sampling time, because kspace positions which are close are highly correlated and can already be recovered using parallel imaging. Variabledensity schemes have several advantages: They equalize the power spectrum of the missing samples, provide graceful degradation in case full recovery is not possible, and can be used for autocalibrating parallel imaging when the kspace center is fully sampled. Methods which combine parallel imaging and compressed sensing represent the state of the art in image reconstruction, as demonstrated by their use in demanding applications such as in pediatric imaging without sedation [75, 76, 12]. Figure 9 shows an image from a pediatric patient reconstructed with parallel imaging compressed sensing at an acceleration factor of about seven.
6 Conclusion
Image reconstruction for parallel imaging can be formulated as an inverse problem. Based on this formulation, advanced iterative algorithms can be developed which (1) make use of optimal (Cartesian or nonCartesian) sampling schemes, and (2) extend parallel imaging with advanced nonlinear regularization terms. These ideas are combined in recent methods for compressed sensing parallel imaging, which currently represent the state of the art in image reconstruction.
References
 [1] M. Hutchinson and U. Raft. Fast MRI data acquisition using multiple detectors. Magn Reson Med, 6:87–91, 1988.
 [2] J. R. Kelton, R. L Magin, and S. M. Wright. An algorithm for rapid image acquisition using multiple receiver coils. In Proc. SMRM, 8th Annual Meeting, page 1172, Amsterdam, 1989.
 [3] D. Kwiat, S. Einav, and G. Navon G. A decoupled coil detector array for fast image acquisition in magnetic resonance imaging. Med Phys, 18:251–265, 1991.
 [4] J. W. Carlson and T. Minemura. Imaging time reduction through multiple receiver coil data acquisition and image reconstruction. Magn Reson Med, 29:681–688, 1993.
 [5] J. B. Ra and C. Y. Rim. Fast imaging using subencoding data sets from multiple detectors. Magn Reson Med, 30:142–145, 1993.
 [6] D. K. Sodickson and W. J. Manning. Simultaneous acquisition of spatial harmonics (SMASH): Fast imaging with radiofrequency coil arrays. Magn Reson Med, 38:591–603, 1997.
 [7] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger. SENSE: Sensitivity encoding for fast MRI. Magn Reson Med, 42:952–962, 1999.
 [8] P. B. Roemer, W. A. Edelstein, C. E. Hayes, S. P. Souza, and O. M. Mueller. The NMR phased array. Magn Reson Med, 16:192–225, 1990.
 [9] F. Wiesinger, P. Boesiger, and K. P. Pruessmann. Electrodynamics and ultimate SNR in parallel MR imaging. Magn Reson Med, 52:376–390, 2004.
 [10] M. Bertero, C. De Mol, and E. R. Pike. Linear inverse problems with discrete data. i. general formulation and singular system analysis. Inverse Problems, 1:301–330, 1985.
 [11] K. M. Johnson, W. F. Block S. B. Reeder, and A. Samsonov. Improved least squares mr image reconstruction using estimates of kspace data consistency. Magn Reson Med, 67:1600–1608, 2012.
 [12] J. Y. Cheng, T. Zhang, N. Ruangwattanapaisarn, M. T. Alley, M. Uecker, J. Pauly, M. Lustig, and S. S. Vasanawala. Freebreathing pediatric MRI with nonrigid motion correction and acceleration. J Magn Reson Imaging, 2014. DOI: 10.1002/jmri.24785.
 [13] D. K. Sodickson and C. A. McKenzie. A generalized approach to parallel magnetic resonance imaging. Med Phys, 28:1629–1643, 2001.
 [14] K. P. Pruessmann, M. Weiger, P. Börnert, and P. Boesiger. Advances in sensitivity encoding with arbitrary kspace trajectories. Magn Reson Med, 46:638–651, 2001.
 [15] B. P. Sutton, D. C. Noll, and J. A. Fessler. Fast, iterative image reconstruction for MRI in the presence of field inhomogeneities. IEEE Trans Med Imaging, 22:178–188, 2003.
 [16] B. P. Sutton, D. C. Noll, and J. A. Fessler. Dynamic field map estimation using a spiralin/spiralout acquisition. Magn Reson Med, 51:1194–1204, 2004.
 [17] C. Liu, M. E. Moseley, and R. Bammer. Simultaneous phase correction and SENSE reconstruction for navigated multishot DWI with noncartesian kspace sampling. Magn Reson Med, 54:1412–1422, 2005.
 [18] M. Uecker, A. Karaus, and J. Frahm. Inverse reconstruction method for segmented multishot diffusionweighted MRI with multiple coils. Magn Reson Med, 62:1342–1348, 2009.
 [19] F. Odille, P.A. Vuissoz, P.Y. Marie, and J. Felblinger. Generalized reconstruction by inversion of coupled systems (GRICS) applied to freebreathing MRI. Magn Reson Med, 60:146–157, 2008.
 [20] J. Y. Cheng, M. T. Alley, C. H. Cunningham, S. S. Vasanawala, J. M. Pauly, and M. Lustig. Nonrigid motion correction in 3D using autofocusing with localized linear translations. Magn Reson Med, 68:1785–1797, 2012.
 [21] V. T. Olafsson, D. C. Noll, and J. A. Fessler. Fast joint reconstruction of dynamic r and field maps in functional MRI. IEEE Trans Med Imaging, 27:1177–1188, 2008.
 [22] C. Graff, Z. Li, A. Bilgin, M. I. Altbach, A. F. Gmitro, and E. W. Clarkson. Iterative T2 estimation from highly undersampled radial fast spinecho data. In Proceedings of the 14th ISMRM Annual Meeting, page 925, Seattle, 2006.
 [23] K. T. Block, M. Uecker, and J. Frahm. Modelbased iterative reconstruction for radial fast spinecho MRI. IEEE Trans Med Imaging, 28:1759–1769, 2009.

[24]
C. L. Welsh, E. V. R. DiBella, G. Adluru, and E. W. Hsu.
Modelbased reconstruction of undersampled diffusion tensor kspace data.
Magn Reson Med, 70:429–440, 2013.  [25] K. King and L. Angelos. SENSE image quality improvement using matrix regularization. In Proceedings of the 9th Annual Meeting of ISMRM, page 1771, Glasgow, 2001.
 [26] F.H. Lin, K. K. Kwong, J. W. Belliveau, and L. L. Wald. Parallel imaging reconstruction using automatic regularization. Magn Reson Med, 51:559–567, 2004.
 [27] G. Adluru, S. P. Awate, T. Tasdizen, R. T. Wihtaker, and E. V. R. DiBella. Temporally constrained reconstruction of dynamic cardiac perfusion MRI. Magn Reson Med, 57:1027–1036, 2007.
 [28] M. Uecker, S. Zhang, D. Voit, K.D. Merboldt, and J. Frahm. Real time MRI: Recent advances using radial FLASH. Imaging in Medicine, 4:461–476, 2012.
 [29] B. Xu, P. Spincemaille, G. Chen, M. Agrawal, T. D. Nguyen, M. R. Prince, and Y. Wang. Fast 3d contrast enhanced MRI of the liver using temporal resolution acceleration with constrained evolution reconstruction. Magn Reson Med, 69:370–381, 2013.
 [30] A. Raj, G. Singh, R. Zabih, B. Kressler, Y. Wang, N. Schuff, and M. Weiner. Bayesian parallel imaging with edgepreserving priors. Magn Reson Med, 57:8–21, 2007.
 [31] K. T. Block, M. Uecker, and J. Frahm. Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint. Magn Reson Med, 57:1086–1098, 2007.
 [32] M. Uecker, K. T. Block, and J. Frahm. Nonlinear inversion with l1wavelet regularization  application to autocalibrated parallel imaging. In Proceedings of the 16th Annual Meeting of the ISMRM, page 1479, Toronto, 2008.
 [33] B. Liu, K. King, M. Steckner, J. Xie, J. Sheng, and L. Ying. Regularized sensitivity encoding (SENSE) reconstruction using Bregman iterations. Magn Reson Med, 61:145–152, 2009.
 [34] V. Athalye, M. Lustig, and M. Uecker. Parallel magnetic resonance imaging as approximation in a reproducing kernel Hilbert space. Inverse Problems, 31:045008, 2015.
 [35] J. Tsao, J. Sánchez, P. Boesiger, and K. P. Pruessmann. Minimumnorm reconstruction for optimal spatial response in highresolution SENSE imaging. In Proceedings of the 11th ISMRM Annual Meeting, page 0014, Toronto, 2003.
 [36] M. Uecker. Nonlinear Reconstruction Methods for Parallel Magnetic Resonance Imaging. PhD thesis, GeorgAugustUniversität Göttingen, 2009.
 [37] K. T. Block, M. Uecker, and J. Frahm. Suppression of MRI truncation artifacts using total variation constrained data extrapolation. International Journal of Biomedical Imaging, 2008:184123, 2008.
 [38] D. Colton and R. Kress. Inverse Acoustic and Electromagnetic Scattering Theory. Springer, Berlin, 1992.
 [39] M. GuerquinKern, L. Lejeune, K. P. Pruessmann, and M. Unser. Realistic analytical phantoms for parallel magnetic resonance imaging. IEEE Trans Med Imaging, 31:626–636, 2012.
 [40] J. D. O’Sullivan. A fast sinc function gridding algorithm for Fourier inversion in computer tomography. IEEE Trans Med Imaging, 4:200–207, 1985.
 [41] J. I. Jackson, C.H. Meyer, D.G. Nishimura, and A. Macovski. Selection of a convolution function for Fourier inversion using gridding. IEEE Trans Med Imaging, 3:473–478, 1991.
 [42] P. J. Beatty, D. G. Nishimura, and J. M. Pauly. Rapid gridding reconstruction with a minimal oversampling ratio. IEEE Trans Med Imaging, 24:799–808, 2005.
 [43] F. Wajer and K. P. Pruessmann. Major speedup of reconstruction for sensitivity encoding with arbitrary trajectories. In Proceedings of the 9th Annual Meeting of the ISMRM, page 767, Glasgow, 2001.
 [44] M. Uecker and M. Lustig. Memorysaving iterative reconstruction on overlapping blocks of kspace. In Proceedings of the 21th Annual Meeting of the ISMRM, page 2645, Salt Lake City, 2013.
 [45] I. Daubechies, M. Defrise, and C. De Mol. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Comm Pure Appl Math, 57:1413–1457, 2004.
 [46] Y. Nesterov. A method of solving a convex programming problem with convergence rate o (1/k2). Soviet Mathematics Doklady, 27:372–376, 1983.
 [47] A. Beck and M. Teboulle. A fast iterative shrinkagethresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences 2.1, pages 183–202, 2009.

[48]
S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.
Distributed optimization and statistical learning via the alternating
direction method of multipliers.
Foundations and Trends in Machine Learning
, 3:1–122, 2011.  [49] M. V. Afonso, J. M. BioucasDias, and M. A. Figueiredo. An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems. IEEE Trans Image Process, 20:681–95, 2011.
 [50] S. Ramani and J. A. Fessler. Parallel MR reconstruction using augmented Lagrangian methods. IEEE Trans Med Imaging, 30:694–706, 2011.
 [51] BART Developers. Berkeley Advanced Reconstruction Toolbox (BART). http://mikgroup.github.io/bart/.
 [52] P. M. Jakob, M. A. Griswold, R. R. Edelman, and D. K. Sodickson. AUTOSMASH: A selfcalibrating technique for SMASH imaging. Magnetic Resonance Materials in Physics, Biology and Medicine, 7:42–54, 1998.
 [53] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and A. Haase. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn Reson Med, 47:1202–1210, 2002.
 [54] C. A. McKenzie, E. N. Yeh, M. A. Ohliger, M. D. Price, and D. K. Sodickson. Selfcalibrating parallel imaging with automatic coil sensitivity extraction. Magn Reson Med, 47:529–538, 2002.
 [55] M. A. Griswold, D. O. Walsh, R. M Heidemann, A. Haase, and P. M. Jakob. The use of an adaptive reconstruction for array coil sensitivity mapping and intensity normalization. In Proceedings of the 10th ISMRM Annual Meeting, page 2410, Honolulu, 2002.
 [56] F. Bauer and S. Kannengiesser. An alternative approach to the image reconstruction for parallel data acquisition in MRI. Math Methods Appl Sci, 30:1437–1451, 2007.
 [57] L. Ying and J. Sheng. Joint image reconstruction and sensitivity estimation in SENSE (JSENSE). Magn Reson Med, 57:1196–1202, 2007.
 [58] M. Uecker, T. Hohage, K. T. Block, and J. Frahm. Image reconstruction by regularized nonlinear inversion – joint estimation of coil sensitivities and image content. Magn Reson Med, 60:674–682, 2008.
 [59] AB. Bakushinsky. Iterative methods for nonlinear operator equations without regularity. new approach. In Dokl. Russian Acad. Sci, volume 330, pages 282–284, 1993.
 [60] F. Knoll, C. Clason, M. Uecker, and R. Stollberger. Improved reconstruction in noncartesian parallel imaging by regularized nonlinear inversion. In Proceedings of the 17th ISMRM Annual Meeting, page 2721, Honolulu, 2009.
 [61] J. Sheng, E. Wiener, B. Liu, F. Boada, and L. Ying. Improved selfcalibrated spiral parallel imaging using JSENSE. Medical Engineering & Physics, 31:510–514, 2009.
 [62] M. Uecker, S. Zhang, and J. Frahm. Nonlinear inverse reconstruction for real‐time MRI of the human heart using undersampled radial FLASH. Magn Reson Med, 63:1456–1462, 2010.
 [63] F. Knoll, C. Clason, K. Bredies, M. Uecker, and R. Stollberger. Parallel imaging with nonlinear reconstruction using variational penalties. Magn Reson Med, 67:34–41, 2012.
 [64] M. Lustig and J. M. Pauly. SPIRiT: Iterative selfconsistent parallel imaging reconstruction from arbitrary kspace. Magn Reson Med, 64:457–471, 2010.
 [65] 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, 71:990–1001, 2014.
 [66] N. Golyandina, A. Korobeynikov, A. Shlemov, and K. Usevich. Multivariate and 2D extensions of singular spectrum analysis with the rssa package. Journal of Statistical Software, 2015. (in press), arXiv:1309.5050 [stat.ME].
 [67] P. J. Shin, P. E. Z. Larson, M. A. Ohliger, M. Elad, J. M. Pauly, D. B. Vigneron, and M. Lustig. Calibrationless parallel imaging reconstruction based on structured lowrank matrix completion. Magn Reson Med, 72:959–970, 2014.
 [68] C. A. McKenzie, M. A. Ohliger, E. N. Yeh, M. D. Price, and D. K. Sodickson. Coilbycoil image reconstruction with SMASH. Magn Reson Med, 46:619–623, 2001.
 [69] M. A. Griswold, S. Kannengiesser, R. M. Heidemann, J. Wang, and P. M. Jakob. Fieldofview limitations in parallel imaging. Magn Reson Med, 52:1118–1126, 2004.

[70]
Z. Wu and R. Schaback.
Local error estimates for radial basis function interpolation of scattered data.
PIMA J Numer Anal, 13:13–27, 1993.  [71] D. L. Donoho. Compressed sensing. IEEE Trans Inform Theory, 52:1289–1306, 2006.
 [72] E. J. Candès, J. Romberg, and T. Tao. Robust uncertainty principle: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inform Theory, 52:489–509, 2006.
 [73] M. Lustig, D. L. Donoho, and J. M. Pauly. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magn Reson Med, 58:1182–1195, 2007.
 [74] R. Otazo, K. Dim, L. Axel, and D. K. Sodickson. Combination of compressed sensing and parallel imaging for highly accelerated firstpass cardiac perfusion MRI. Magn Reson Med, 64:767–776, 2010.
 [75] S. S. Vasanawala, M. T. Alley, B. A. Hargreaves, R. A. Barth, J. M. Pauly, and M. Lustig. Improved pediatric MR imaging with compressed sensing. Radiology, 256:607–616, 2010.
 [76] T. Zhang, J. Y. Cheng, A. G. Potnick, R. A. Barth, M. T. Alley, M. Uecker, M. Lustig, J. M. Pauly, and S. S. Vasanawala. Fast pediatric 3D free breathing abdominal dynamic contrast enhanced MRI with a high spatiotemporal resolution. J Magn Reson Imaging, 41:460–473, 2015.