1 Introduction
Magnetic resonance imaging (MRI) is a medical imaging approach that uses magnetic fields to create detailed anatomical images. Although MRI provides excellent softtissue contrast without the use of ionizing radiation, it takes a long time to fully sample the measurement space. Thus, it is common to take relatively few measurements and apply sophisticated postprocessing to reconstruct an accurate image. Although our paper focuses on MRI, the methods we propose apply to any application where the goal is to recover a signal from undersampled Fourier measurements.
The measurements collected in coil MRI, known as “kspace” measurements, can be modeled as
(1) 
where
is a vectorized version of the
pixel image we wish to recover,is a unitary 2D discrete Fourier transform (DFT),
is a sampling mask formed fromrows of the identity matrix
, is the th coilsensitivity map, and is additive white Gaussian noise (AWGN) with precision(i.e., variance
). In the special case of singlecoil MRI, and , the allones vector. In this paper, we use the variabledensity sampling masks like that shown in Fig. 1. In MRI, the pixelstomeasurement ratio, , is known as the “acceleration rate.” When , one cannot uniquely determine from due to the nullspace of , and so prior information about is needed for its recovery.To recover images from undersampled MRI measurements, many methods have been proposed. Some are based on iterative optimization [1, 2], as described in the sequel. More recently, deep neural networks (DNNs) that directly map MRI measurements to an image have been proposed, e.g., [3, 4]. Although such DNNs work well, training them requires huge fullysampled kspace datasets (which may be difficult or impossible to obtain) and changes in the acquisition parameters (e.g., sampling mask) from training to testing can degrade performance [5].
In this work, we focus on the “plug and play” approach that iteratively calls a DNN for image denoising, which brings several advantages. First, DNN denoisers can be trained using image patches, implying the need for relatively few images and no kspace data. Second, the denoiser is trained independently of the acquisition parameters, so that it generalizes to any acquisition scenario. Our approach is based on the generalized expectation consistent (GEC) approximation algorithm from [6], which lives in the family of approximate message passing (AMP) algorithms like [7, 8].
2 Background
2.1 Compressedsensingbased methods
The conventional approach to MRI image recovery [1, 2] is to pose and solve an optimization problem of the form
(2) 
where promotes measurement fidelity and is an imagebased regularizer. Typical choices are
(3) 
for (1) and with a suitable transform (e.g., wavelet or totalvariation) and carefully chosen . Such encourage sparsity in the transform coefficients .
Many algorithms have been proposed to solve (2) with convex and [2]. For example, the alternating directions method of multipliers (ADMM) [9] iterates
(4a)  
(4b)  
(4c) 
where the proximal operator is defined as
(5) 
In (4), is a tunable stepsize that affects the speed of ADMM’s convergence but not its fixed point.
2.2 Plugandplay methods
The prox operator (5
) can interpreted as a denoiser, in particular, the maximum a posteriori (MAP) estimator of
with prior from an observation with precision AWGN . Leveraging this fact, Bouman et al. [10] proposed to replace ADMM line (4b) with a call to a highperformance image denoiser like BM3D [11] or DnCNN [12], giving rise to “plugandplay” (PnP) ADMM. PnP extensions of other algorithms, such as primaldual splitting (PDS) [13] and proximal gradient (PG) [14], have also been proposed. As shown in the recent overview paper, PnP methods have been shown to significantly outperform compressedsensingbased approaches in MRI [5]. Note, however, that when (4b) is replaced with a denoising step of the form “,” the stepsize does affect the fixedpoint [5] and thus must be tuned.Although PnP algorithms work well in MRI, there is room for improvement. For example, while image denoisers are typically designed/trained to remove the effects of AWGN, PnP algorithms do not provide the denoiser with an AWGNcorrupted input at each iteration. Rather, the denoiser’s input error has iterationdependent statistics that are difficult to analyze or predict.
2.3 Approximate message passing
In (2), if we interpret as a loglikelihood and as a log prior, then can be interpreted as the MAP estimate of from . However, because image recovery results are often judged by meansquared error (MSE), one may be more interested in the minimum MSE (MMSE) estimate of from . Interestingly, both MMSE and MAP estimation are facilitated by approximate message passing (AMP) methods like [7, 8].
For example, the AMP algorithm from [7] iterates
(6a)  
(6b)  
(6c) 
over , starting from , where is a Lipschitz denoising function, is the trace of the Jacobian of at , and . When configured for MAP estimation, AMP uses the MAP denoiser . When configured for MMSE estimation, AMP instead uses the MMSE denoiser for with .
Importantly, when the forward operator is large and i.i.d. subGaussian, and when is Lipschitz, the macroscopic behavior of AMP is rigorously characterized by a scalar stateevolution [15, 16]. When is the MMSE denoiser and the stateevolution has a unique fixed point, AMP provably converges to the MMSEoptimal estimate [15, 16]. For images, the MMSE denoiser can be approximated by BM3D or a DNN, as proposed in [17], leading to “denoisingAMP” (DAMP). There, the traceJacobian in (6a) is approximated using the MonteCarlo approach [18]
(7) 
with random and small .
More recently, the vector AMP (VAMP) algorithm [8] was proposed, with similar properties as AMP (e.g., rigorous state evolution and provable MMSE estimation) but applicability to a wider class of random matrices: right orthogonally invariant (ROI) ones. Inspired by DAMP, a denoising VAMP (DVAMP) was proposed in [19] and analyzed in [20].
2.4 AMP for MRI
Neither AMP nor VAMP works as intended in MRI because in (1) lacks sufficient randomness. In fact, these algorithms tend to diverge in MRI if applied without modification.
The failure of AMP and VAMP can be understood from their error recursions. For AMP, the error recursion is [21]
(8)  
(9) 
It is important to keep in mind that images have much more energy at low Fourier frequencies than at high ones. The same tends to be true of the output error of an image denoiser. Even so, if was large and i.i.d. (with zero mean and elementwise variance ), then the term in (8) would randomize such that the denoiser input error vector looks like AWGN. In MRI, however, both and have Fourier structure, this randomization does not happen, and AMP behaves unpredictably. A similar behavior plagues VAMP.
Several MRIspecific variations of AMP and VAMP have been proposed to counter these deficiencies. For example, [22] proposed DAMP with a very small , which helps the algorithm converge, but at the cost of degrading its fixed points. [23] proposed a damped DVAMP that, combined with a novel initialization, showed improved performance and runtime over PnPADMM for MRI.
Several other VAMPbased algorithms for MRI have been designed to recover the wavelet coefficients of the image rather than the image itself. The motivation is that, in this case, is a FourierWavelet matrix, which is approximately block diagonal [24], where the blocks correspond to the wavelet subbands. With an appropriate modification of VAMP, the subband error vectors can be made to behave more like AWGN, albeit with different variances. The first incarnation of this idea appeared in [25], where a fixed bandwise normalization procedure was used. Later, for singlecoil MRI with variabledensity sampling masks, a “variable density AMP” (VDAMP) algorithm with bandspecific adaptive wavelet thresholding was proposed in [26], which was able to successfully predict the noise variance in each subband at each iteration. More recently, the DVDAMP algorithm [27] extended VDAMP to DNN denoising in each subband.
Although DVDAMP is the stateoftheart AMP algorithm for MRI, it is based on a nonstandard modification of VAMP with degraded fixed points, which makes early stopping critical for good performance. Also, it is not clear how to extend DVDAMP to multicoil MRI. These issues motivate our approach, which is described next.
3 Proposed Approach
3.1 Denoising GEC
Our approach uses the GEC framework from [6], which is summarized in Alg. 1. When solving a convex optimization problem of the form (2), the functions in Alg. 1 take the form
(10) 
for the generalized proximal operator
(11) 
where and creates a diagonal matrix from its vector argument. Note that if , then . Furthermore, if the vectors were held fixed over the iterations and took the form , then Alg. 1 reduces to a variant of ADMM (4) with two dual updates: (4c) and a similar step between (4a) and (4b). So, GEC can be interpreted as an ADMMlike algorithm with two adaptive vectorvalued stepsizes, and .
In lines 6 and 11, GEC averages the diagonal of the Jacobian separately over coefficient subsets using :
(12) 
In (12), is the size of the th subset and is the th diagonal subblock of the matrix input . When , GEC reduces to VAMP.
We focus on the quadratic loss (3), which yields
(13) 
For , we propose to “plug in” a DNN denoiser. For both and , we approximate the term in (12) using
(14) 
where the th coefficient subset in is i.i.d. unitvariance Gaussian and the others are zero. Inspired by DAMP and DVAMP, we call this approach “denoising GEC” (DGEC).
3.2 DGEC for Image Recovery
Like [25, 26, 27], we recover the wavelet coefficients rather than the image pixels . For orthogonal wavelet transform , we have and , so that we can rewrite (1) as
(15) 
To apply DGEC to recovery, we choose as in (13), but with in place of , and for the diagonalization subsets we choose the subbands of a depth 2D wavelet transform. As in [27], we perform denoising in the wavelet domain using a denoiser that can exploit knowledge of the noise variance in each wavelet subband, as provided by the precision vector .
4 Numerical Experiments
We now compare the proposed DGEC algorithm to the existing DVDAMP [27] and PnPPDS [13] algorithms. Based on the extensive experiments in [27], DVDAMP is stateoftheart among PnP algorithms. PnPPDS is a useful baseline, since it has the same fixed points as PnPADMM and PnPPG.
Denoisers: For a fair comparison to DVDAMP [27], we use the DNN denoiser proposed in [27], which is modification of DnCNN [12]
that accepts the noise standard deviation (SD) in each wavelet subband. The denoiser was trained using noise that was white in each subband but with SD that varies across subbands. In particular, 5 copies of the denoiser were trained using subband noise SDs uniformly distributed in the ranges 010, 1020, 2050, 50120, and 120500, respectively. (Pixel values ranged from 0255.) The DNNs were trained using patches from 70 MRI images of the Stanford 2D FSE dataset available at
http://mridata.org. For PnPPDS, we used a standard DnCNN denoiser trained on the same data with white noise of SD uniformly distributed in 2050, as in
[27]. Because we used realvalued images, the denoisers use only the real part of the input and generate a realvalued output.Test data: For evaluation, we used the ten 352352 MRI images in Fig. 2, which were not in the training dataset. The measurements were constructed using (1) with complex AWGN whose variance was adjusted to give a premasking SNR of 40 dB. For the multicoil experiments, we used coil sensitivities simulated using the BiotSavart law, while in the singlecoil case, we used .
Algorithm parameters: For DGEC and DVDAMP, we used a 2D Haar wavelet transform with levels, giving wavelet subbands. DGEC used the autotuning scheme from [28] and the damping scheme from [23] with parameter . DVDAMP code was obtained from the authors and run under default settings, which are detailed in [27]. PnPPDS was run for 200 iterations using the stepsize that maximized PSNR on the training set.
coil  coils  

method  PSNR  SSIM  PSNR  SSIM  PSNR  SSIM  PSNR  SSIM 
PnPPDS  40.66  0.968  37.38  0.951  34.71  0.935  33.09  0.917 
DVDAMP  42.36  0.972  35.92  0.918  n/a  n/a  n/a  n/a 
DGEC  42.97  0.977  37.65  0.946  45.18  0.993  41.13  0.982 
Singlecoil results: Table 1 shows that DGEC outperformed DVDAMP in all singlecoil experiments and outperformed PnPPDS in all but SSIM at . Figure 5 shows an example of the wavelet coefficients input to DGEC’s denoiser at the 10th iteration, and their error relative to the true coefficients. Figure 3 shows the evolution of the standard deviation at the input to DGEC’s denoiser in each subband; there is a good agreement between true and predicted values. Figure 6 suggests that the subband errors are Gaussian. Figure 4 shows image recoveries and error maps for one test image at .
Multicoil results: Table 1 shows DGEC significantly outperforming PnPPDS in PSNR and SSIM in the 4coil case. DVDAMP does not support multicoil recovery and thus is not shown.
5 Conclusion
We designed a GECbased plugandplay algorithm for MRI that renders the subband errors white and Gaussian with predictable variance, and used it with a denoiser trained to handle subband errors that are white and Gaussian with known variance. Experiments show good performance relative to previous approaches in single and multicoil settings.
References
 [1] M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, ‘‘Compressed sensing MRI,’’ IEEE Signal Process. Mag., vol. 25, no. 2, pp. 7282, Mar. 2008.
 [2] J. A. Fessler, ‘‘Optimization methods for magnetic resonance image reconstruction,’’ IEEE Signal Process. Mag., vol. 37, no. 1, pp. 3340, 2020.

[3]
K. H. Jin, M. T. McCann, E. Froustey, and M. Unser, ‘‘Deep convolutional neural network for inverse problems in imaging,’’
IEEE Trans. Image Process., vol. 26, no. 9, pp. 45094522, 2017.  [4] K. Hammernik, T. Klatzer, E. Kobler, M. P. Rec ht, D. K. Sodickson, T. Pock, and F. Knoll, ‘‘Learning a variational network for reconstruction of accelerated MRI data,’’ Magnetic Resonance Med., vol. 79, no. 6, pp. 30553071, 2018.
 [5] R. Ahmad, C. A. Bouman, G. T. Buzzard, S. Chan, S. Liu, E. T. Reehorst, and P. Schniter, ‘‘Plug and play methods for magnetic resonance imaging,’’ IEEE Signal Process. Mag., vol. 37, no. 1, pp. 105116, 2020.
 [6] A. K. Fletcher, M. SahraeeArdakan, S. Rangan, and P. Schniter, ‘‘Expectation consistent approximate inference: Generalizations and convergence,’’ in Proc. IEEE Int. Symp. Inform. Thy., 2016, pp. 190194.
 [7] D. L. Donoho, A. Maleki, and A. Montanari, ‘‘Message passing algorithms for compressed sensing,’’ Proc. Nat. Acad. Sci., vol. 106, no. 45, pp. 18 91418 919, Nov. 2009.
 [8] S. Rangan, P. Schniter, and A. K. Fletcher, ‘‘Vector approximate message passing,’’ IEEE Trans. Inform. Theory, pp. 66646684, 2019.
 [9] 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, pp. 1122, 2011.
 [10] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, ‘‘Plugandplay priors for model based reconstruction,’’ in Proc. IEEE Global Conf. Signal Info. Process., 2013, pp. 945948.
 [11] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, ‘‘Image denoising by sparse 3D transformdomain collaborative filtering,’’ IEEE Trans. Image Process., vol. 16, no. 8, pp. 20802095, 2007.
 [12] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang, ‘‘Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,’’ IEEE Trans. Image Process., vol. 26, no. 7, pp. 31423155, 2017.
 [13] S. Ono, ‘‘Primaldual plugandplay image restoration,’’ IEEE Signal Process. Lett., vol. 24, no. 8, pp. 11081112, 2017.
 [14] U. Kamilov, H. Mansour, and B. Wohlberg, ‘‘A plugandplay priors approach for solving nonlinear imaging inverse problems,’’ IEEE Signal Process. Lett., vol. 24, no. 12, pp. 18721876, May 2017.
 [15] R. Berthier, A. Montanari, and P.M. Nguyen, ‘‘State evolution for approximate message passing with nonseparable functions,’’ Inform. Inference, 2019.
 [16] M. Bayati and A. Montanari, ‘‘The dynamics of message passing on dense graphs, with applications to compressed sensing,’’ IEEE Trans. Inform. Theory, vol. 57, no. 2, pp. 764785, Feb. 2011.
 [17] C. A. Metzler, A. Maleki, and R. G. Baraniuk, ‘‘BM3DAMP: A new image recovery algorithm based on BM3D denoising,’’ in Proc. IEEE Int. Conf. Image Process., 2015, pp. 31163120.
 [18] S. Ramani, T. Blu, and M. Unser, ‘‘MonteCarlo SURE: A blackbox optimization of regularization parameters for general denoising algorithms,’’ IEEE Trans. Image Process., vol. 17, no. 9, pp. 15401554, 2008.
 [19] P. Schniter, S. Rangan, and A. K. Fletcher, ‘‘Denoisingbased vector approximate message passing,’’ in Proc. Intl. Biomed. Astronom. Signal Process. (BASP) Frontiers Workshop, 2017, p. 77.
 [20] A. K. Fletcher, P. Pandit, S. Rangan, S. Sarkar, and P. Schniter, ‘‘Plugin estimation in highdimensional linear inverse problems: A rigorous analysis,’’ in Proc. Neural Inform. Process. Syst. Conf., 2018, pp. 74407449.
 [21] P. Schniter, ‘‘A simple derivation of AMP and its state evolution via firstorder cancellation,’’ IEEE Trans. Signal Process., vol. 68, pp. 42834292, 2020.
 [22] E. M. Eksioglu and A. K. Tanc, ‘‘Denoising AMP for MRI reconstruction: BM3DAMPMRI,’’ SIAM J. Imag. Sci., vol. 11, no. 3, pp. 20902109, 2018.
 [23] S. Sarkar, R. Ahmad, and P. Schniter, ‘‘MRI image recovery using damped denoising vector AMP,’’ in Proc. IEEE Int. Conf. Acoust. Speech & Signal Process., 2021, pp. 81088112.
 [24] B. Adcock, A. C. Hansen, C. Poon, and B. Roman, ‘‘Breaking the coherence barrier: A new theory for compressed sensing,’’ Forum of Mathematics, Sigma, vol. 5, no. E4, doi:10.1017/fms.2016.32.
 [25] P. Schniter, S. Rangan, and A. K. Fletcher, ‘‘Plugandplay image recovery using vector AMP,’’ presented at the Intl. Biomedical and Astronomical Signal Processing (BASP) Frontiers Workshop, VillarssurOllon, Switzerland, Jan. 2017. [Online]. Available: http://www2.ece.ohiostate.edu/~schniter/pdf/basp17_poster.pdf
 [26] C. Millard, A. T. Hess, B. Mailhé, and J. Tanner, ‘‘Approximate message passing with a colored aliasing model for variable density Fourier sampled images,’’ arXiv:2003.02701, 2020.
 [27] C. A. Metzler and G. Wetzstein, ‘‘DVDAMP: Denoisingbased approximate message passing for compressive MRI,’’ in Proc. IEEE Int. Conf. Acoust. Speech & Signal Process., 2021, pp. 14101414.
 [28] A. K. Fletcher, M. SahraeeArdakan, S. Rangan, and P. Schniter, ‘‘Rigorous dynamics and consistent estimation in arbitrarily conditioned linear systems,’’ in Proc. Neural Inform. Process. Syst. Conf., 2017, pp. 25422551.