Plug and play methods for magnetic resonance imaging (long version)

TL;DR

This paper explores plug-and-play algorithms for MRI image reconstruction, emphasizing their iterative denoising approach, equilibrium interpretation, and demonstrating their effectiveness in recovering images from undersampled data.

Contribution

It introduces a unifying framework for PnP methods in MRI, including convergence analysis and practical examples of image recovery from undersampled datasets.

Findings

PnP methods effectively reconstruct MRI images from undersampled data

Convergence of PnP algorithms can be analyzed via equilibrium equations

Illustrative examples demonstrate practical success in MRI recovery

Abstract

Magnetic Resonance Imaging (MRI) is a non-invasive diagnostic tool that provides excellent soft-tissue contrast without the use of ionizing radiation. Compared to other clinical imaging modalities (e.g., CT or ultrasound), however, the data acquisition process for MRI is inherently slow, which motivates undersampling and thus drives the need for accurate, efficient reconstruction methods from undersampled datasets. In this article, we describe the use of "plug-and-play" (PnP) algorithms for MRI image recovery. We first describe the linearly approximated inverse problem encountered in MRI. Then we review several PnP methods, where the unifying commonality is to iteratively call a denoising subroutine as one step of a larger optimization-inspired algorithm. Next, we describe how the result of the PnP method can be interpreted as a solution to an equilibrium equation, allowing convergence analysis from the equilibrium perspective. Finally, we present illustrative examples of PnP methods applied to MRI image recovery.

Figures

Equations

$${\mathbf{y}}_i^{(t)}={\mathbf{P}}^{(t)}\mathbf{F}{\mathbf{S}}_i{\mathbf{x}}^{(t)}+{\mathbf{w}}_i^{(t)},$$

$$\mathbf{y}=\mathbf{Ax}+\mathbf{w},\mathbf{w}\sim \mathcal{N}(0,{\sigma }^2\mathbf{I}),$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}}{\mathrm{argmin}}\left\{\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{A}\mathbf{x}{\Vert }_2^2+\varphi (\mathbf{x})\right\},$$

$$\displaystyle\widehat{\boldsymbol{x}}_{\text{\sf map}}=\arg\min_{\boldsymbol{x}}\big{\{}-\ln p(\boldsymbol{y}|\boldsymbol{x})-\ln p(\boldsymbol{x})\big{\}}.$$

$$\mathbf{z}=\mathbf{x}+\mathbf{w},\mathbf{w}\sim \mathcal{N}(0,{\sigma }^2\mathbf{I}).$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}\in {\mathbb{C}}^N}{\mathrm{argmin}}\underset{\mathbf{v}\in {\mathbb{C}}^N}{\min }\left\{\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{A}\mathbf{x}{\Vert }_2^2+\varphi (\mathbf{v})\right\}\text{mbox}subjectto\mathbf{x}=\mathbf{v}.$$

$$\underset{\mathbf{x},\mathbf{v}}{\min }\underset{\mathbf{\lambda }}{\max }\left\{\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }_2^2+\varphi (\mathbf{v})+\mathrm{Re}\{{\mathbf{\lambda }}^{\text{H}}(\mathbf{x}-\mathbf{v})\}+\frac{1}{2\eta }\Vert \mathbf{x}-\mathbf{v}{\Vert }^2\right\},$$

$$\underset{\mathbf{x},\mathbf{v}}{\min }\underset{\mathbf{u}}{\max }\left\{\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }_2^2+\varphi (\mathbf{v})+\frac{1}{2\eta }\Vert \mathbf{x}-\mathbf{v}+\mathbf{u}{\Vert }^2-\frac{1}{2\eta }\Vert \mathbf{u}{\Vert }^2\right\}.$$

$${\mathbf{x}}_k$$

$${\mathbf{v}}_k$$

$${\mathbf{u}}_k$$

$${\mathrm{prox}}_{\varphi }(\mathbf{z};\eta )$$

$$\mathbf{h}(\mathbf{z};\eta )$$

$$={\left({\mathbf{A}}^{\text{H}}\mathbf{A}+\frac{{\sigma }^2}{\eta }\mathbf{I}\right)}^{-1}\left({\mathbf{A}}^{\text{H}}\mathbf{y}+\frac{{\sigma }^2}{\eta }\mathbf{z}\right).$$

$${\mathbf{x}}_k$$

$${\mathbf{v}}_k$$

$${\mathbf{u}}_k$$

$${\mathbf{z}}_k$$

$${\mathbf{x}}_k$$

$${\mathbf{s}}_k$$

$${\mathbf{x}}_k$$

$${\mathbf{v}}_k$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}}{\mathrm{argmin}}\left\{\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }_2^2+\lambda \Vert \mathbf{\Psi }\mathbf{x}{\Vert }_1\right\},$$

$${\mathrm{prox}}_{\varphi }(\mathbf{z};\eta )=\underset{\mathbf{x}}{\mathrm{argmin}}\left\{\lambda \Vert \mathbf{\Psi }\mathbf{x}{\Vert }_1+\frac{1}{2\eta }\Vert \mathbf{x}-\mathbf{z}{\Vert }^2\right\}.$$

$${\mathbf{f}}_{\text{tdt}}(\mathbf{z};\tau )\triangleq {\mathbf{\Psi }}^{\text{H}}\text{soft-thresh}(\mathbf{\Psi }\mathbf{z};\tau )$$

$$\hat{\mathbf{x}}={\mathbf{\Psi }}^{\text{H}}\hat{\mathbf{\alpha }}\text{for}\hat{\mathbf{\alpha }}=\underset{\mathbf{\alpha }\in \mathrm{range}(\mathbf{\Psi })}{\mathrm{argmin}}\left\{\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-{\mathbf{A}\mathbf{\Psi }}^{\text{H}}\mathbf{\alpha }{\Vert }_2^2+\lambda \Vert \mathbf{\alpha }{\Vert }_1\right\},$$

$$\hat{\mathbf{x}}={\mathbf{\Psi }}^{\text{H}}\hat{\mathbf{\alpha }}\text{for}\hat{\mathbf{\alpha }}=\underset{\mathbf{\alpha }}{\mathrm{argmin}}\left\{\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-{\mathbf{A}\mathbf{\Psi }}^{\text{H}}\mathbf{\alpha }{\Vert }_2^2+\frac{\beta }{2}\Vert {\mathbf{P}}_{\mathbf{\Psi }}^{\perp }\mathbf{\alpha }{\Vert }_2^2+\lambda \Vert \mathbf{\alpha }{\Vert }_1\right\}$$

$$0=\frac{1}{{\sigma }^2}{\mathbf{A}}^{\text{T}}(\mathbf{A}\hat{\mathbf{x}}-\mathbf{y})+\frac{1}{\eta }(\hat{\mathbf{x}}-\mathbf{f}(\hat{\mathbf{x}})),$$

$$\displaystyle\widehat{\boldsymbol{x}}=\operatorname*{arg\,min}_{\boldsymbol{x}}\left\{\frac{1}{2\sigma^{2}}\|\boldsymbol{y}-\boldsymbol{Ax}\|^{2}+\phi_{\text{\sf red}}(\boldsymbol{x})\right\}\text{~{}with~{}}\phi_{\text{\sf red}}(\boldsymbol{x})\triangleq\frac{1}{2\eta}\boldsymbol{x}^{\text{{T}}}\big{(}\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})\big{)},$$

$$\widehat{p_{\text{x}}}(\mathbf{x})$$

$$\widetilde{p_{\text{x}}}(\mathbf{x};\eta )$$

$$\hat{\mathbf{x}}=\arg \underset{\mathbf{x}}{\min }\left\{\frac{1}{2{\sigma }^2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }^2-\ln \widetilde{p_{\text{x}}}(\mathbf{x};\eta )\right\}.$$

$$0=\frac{1}{{\sigma }^2}{\mathbf{A}}^{\text{T}}(\mathbf{A}\hat{\mathbf{x}}-\mathbf{y})-\nabla \ln \widetilde{p_{\text{x}}}(\hat{\mathbf{x}};\eta ).$$

$$\nabla \ln \widetilde{p_{\text{x}}}(\mathbf{z};\eta )$$

$$\displaystyle\boldsymbol{0}=\frac{1}{\sigma^{2}}\boldsymbol{A}^{\text{{T}}}(\boldsymbol{A}\widehat{\boldsymbol{x}}-\boldsymbol{y})+\frac{1}{\eta}\big{(}\widehat{\boldsymbol{x}}-\boldsymbol{f}_{\text{\sf mmse}}(\widehat{\boldsymbol{x}};\eta)\big{)},$$

Full text

I Introduction

Magnetic Resonance Imaging (MRI) uses radiofrequency waves to non-invasively evaluate the structure, function, and morphology of soft tissues. MRI has become an indispensable imaging tool for diagnosing and evaluating a host of conditions and diseases. Compared to other clinical imaging modalities (e.g., CT or ultrasound), however, MRI suffers from slow data acquisition. A typical clinical MRI exam consists of multiple scans and can take more than an hour to complete. For each scan, the patient may be asked to stay still for several minutes, with slight motion potentially resulting in image artifacts. Furthermore, dynamic applications demand collecting a series of images in quick succession. Due to the limited time window in many dynamic applications (e.g., contrast enhanced MR angiography), it is not feasible to collect fully sampled datasets. For these reasons, MRI data is often undersampled. Consequently, computationally efficient methods to recover high-quality images from undersampled MRI data have been actively researched for the last two decades.

The combination of parallel (i.e., multi-coil) imaging and compressive sensing (CS) has been shown to benefit a wide range of MRI applications [1] [2], including dynamic applications, and has been included in the default image processing frameworks offered by several major MRI vendors. More recently, learning-based techniques (e.g.,  [3, 4][5, 6, 7, 8, 9] [10, 11, 12, 13, 14]) have been shown to outperform CS methods. Some of these techniques learn the entire end-to-end mapping from undersampled k-space or aliased images to recovered images (e.g.,  [4][7]). Considering that the forward model in MRI changes from one dataset to the next, such end-to-end methods must be either trained over a large and diverse data corpus or limited to a specific application. Other methods train scan-specific convolutional neural networks (CNN) on a fully-sampled region of k-space and then use it to interpolate missing k-space samples [8]. These methods do not require separate training data but demand a fully sampled k-space region. Due to the large number of unknown CNN parameters, such methods require a fully sampled region that is larger than that typically acquired in parallel imaging, limiting the acceleration that can be achieved. Other supervised learning methods are inspired by classic variational optimization methods and iterate between data-consistency enforcement and a trained CNN, which acts as a regularizer [6] [11]. Such methods require a large number of fully sampled, multi-coil k-space datasets, which may be difficult to obtain in many applications. Also, since CNN training occurs in the presence of dataset-specific forward models, generalization from training to test scenarios remains an open question [9]. Other learning-based methods have been proposed based on bi-level optimization (e.g., [10]), adversarially learned priors (e.g., [12]), and autoregressive priors (e.g., [13]). Consequently, the integration of learning-based methods into physical inverse problems remains a fertile area of research. There are many directions for improvement, including recovery fidelity, computational and memory efficiency, robustness, interpretability, and ease-of-use.

This article focuses on “plug-and-play” (PnP) algorithms [15], which alternate image denoising with forward-model based signal recovery. PnP algorithms facilitate the use of state-of-the-art image models through their manifestations as image denoisers, whether patch-based (e.g., [16]) or deep neural network (DNN) based (e.g., [17] [18]). The fact that PnP algorithms decouple image modeling from forward modeling has advantages in compressive MRI, where the forward model can change significantly among different scans due to variations in the coil sensitivity maps, sampling patterns, and image resolution. Furthermore, fully sampled k-space MRI data is not needed for PnP; the image denoiser can be learned from MRI image patches, or possibly even magnitude-only patches. The objective of this article is two-fold: i) to review recent advances in plug-and-play methods, and ii) to discuss their application to compressive MRI image reconstruction.

The remainder of the paper is organized as follows. We first detail the inverse problem encountered in MRI reconstruction. We then review several PnP methods, where the unifying commonality is to iteratively call a denoising subroutine as one step of a larger optimization-inspired algorithm. Next, we describe how the result of the PnP method can be interpreted as a solution to an equilibrium equation, allowing convergence analysis from the equilibrium perspective. Finally, we present illustrative examples of PnP methods applied to MRI image recovery.

II Image recovery in compressive MRI

In this section, we describe the standard linear inverse problem formulation in MRI. We acknowledge that more sophisticated formulations exist (see, e.g., [19] for a more careful modeling of physics effects). Briefly, the measurements are samples of the Fourier transform of the image, where the Fourier domain is often referred to as “k-space.” The transform can be taken across two or three spatial dimensions and includes an additional temporal dimension in dynamic applications. Furthermore, measurements are often collected in parallel from $C\geq 1$ receiver coils. In dynamic parallel MRI with Cartesian sampling, the time-$t$ k-space measurements from the $i^{\text{th}}$ coil take the form

[TABLE]

where $\boldsymbol{x}^{(t)}\in{\mathbb{C}}^{N}$ is the vectorized 2D or 3D image at discrete time $t$, $\boldsymbol{S}_{i}\in{\mathbb{C}}^{N\times N}$ is a diagonal matrix containing the sensitivity map for the $i^{\text{th}}$ coil, $\boldsymbol{F}\in{\mathbb{C}}^{N\times N}$ is the 2D or 3D discrete Fourier transform (DFT), the sampling matrix $\boldsymbol{P}^{(t)}\in{\mathbb{R}}^{M\times N}$ contains $M$ rows of the $N\times N$ identity matrix, and $\boldsymbol{w}_{i}^{(t)}\in{\mathbb{C}}^{M}$ is additive white Gaussian noise (AWGN). Often the sampling pattern changes across frames $t$. The MRI literature often refers to $R\triangleq N/M$ as the “acceleration rate.” The AWGN assumption, which does not hold for the measured parallel MRI data, is commonly enforced using noise pre-whitening filters, which yields the model (1) but with diagonal “virtual” coil maps $\boldsymbol{S}_{i}$ [20] [21]. Additional justification of the AWGN model can be found in [22].

MRI measurements are acquired using a sequence of measurement trajectories through k-space. These trajectories can be Cartesian or non-Cartesian in nature. Cartesian trajectories are essentially lines through k-space. In the Cartesian case, one k-space dimension (i.e., the frequency encoding) is fully sampled, while the other one or two dimensions (i.e., the phase encodings) are undersampled to reduce acquisition time. Typically, one line, or “readout,” is collected after every RF pulse, and the process is repeated several times to collect adequate samples of k-space. Non-Cartesian trajectories include radial or spiral curves, which have the effect of distributing the samples among all dimensions of k-space. Compared to Cartesian sampling, non-Cartesian sampling provides more efficient coverage of k-space and yields an “incoherent” forward operator that is more conducive to compressed-sensing reconstruction  [23]. But Cartesian sampling remains the method of choice in clinical practice, due to its higher tolerance to system imperfections and an extensive record of success.

Since the sensitivity map, $\boldsymbol{S}_{i}$, is patient-specific and varies with the location of the coil with respect to the imaging plane, both $\boldsymbol{S}_{i}$ and $\boldsymbol{x}^{(t)}$ are unknown in practice. Although calibration-free methods have been proposed to estimate $\boldsymbol{S}_{i}\boldsymbol{x}^{(t)}$ (e.g., [24]) or to jointly estimate $\boldsymbol{S}_{i}$ and $\boldsymbol{x}^{(t)}$ (e.g., [25]), it is more common to first estimate $\boldsymbol{S}_{i}$ through a calibration procedure and then treat $\boldsymbol{S}_{i}$ as known in (1). Stacking $\{\boldsymbol{y}_{i}^{(t)}\}$, $\{\boldsymbol{x}^{(t)}\}$, and $\{\boldsymbol{w}_{i}^{(t)}\}$ into vectors $\boldsymbol{y}$, $\boldsymbol{x}$, and $\boldsymbol{w}$, and packing $\{\boldsymbol{P}^{(t)}\boldsymbol{F}\boldsymbol{S}_{i}\}$ into a known block-diagonal matrix $\boldsymbol{A}$, we obtain the linear inverse problem of recovering $\boldsymbol{x}$ from

[TABLE]

where $\mathcal{N}(\boldsymbol{0},\sigma^{2}\boldsymbol{I})$ denotes a circularly symmetric complex-Gaussian random vector.

III Signal Recovery and Denoising

The maximum likelihood (ML) estimate of $\boldsymbol{x}$ from $\boldsymbol{y}$ in (2) is $\widehat{\boldsymbol{x}}_{\text{\sf ml}}\triangleq\arg\max_{\boldsymbol{x}}p(\boldsymbol{y}|\boldsymbol{x})$, where $p(\boldsymbol{y}|\boldsymbol{x})$, the probability density of $\boldsymbol{y}$ conditioned on $\boldsymbol{x}$, is known as the “likelihood function.” The ML estimate is often written in the equivalent form $\widehat{\boldsymbol{x}}_{\text{\sf ml}}=\arg\min_{\boldsymbol{x}}\{-\ln p(\boldsymbol{y}|\boldsymbol{x})\}$. In the case of $\sigma^{2}$-variance AWGN $\boldsymbol{w}$, we have that $-\ln p(\boldsymbol{y}|\boldsymbol{x})=\frac{1}{2\sigma^{2}}\|\boldsymbol{y}-\boldsymbol{Ax}\|_{2}^{2}+\mathrm{const}$, and so $\widehat{\boldsymbol{x}}_{\text{\sf ml}}=\arg\min_{\boldsymbol{x}}\|\boldsymbol{y}-\boldsymbol{Ax}\|_{2}^{2}$, which can be recognized as least-squares estimation. Although least-squares estimation can give reasonable performance when $\boldsymbol{A}$ is tall and well conditioned, this is rarely the case under moderate to high acceleration (i.e., $R>2$). With acceleration, it is critically important to exploit prior knowledge of signal structure.

The traditional approach to exploiting such prior knowledge is to formulate and solve an optimization problem of the form

[TABLE]

where the regularization term $\phi(\boldsymbol{x})$ encodes prior knowledge of $\boldsymbol{x}$. In fact, $\widehat{\boldsymbol{x}}$ in (3) can be recognized as the maximum a posteriori (MAP) estimate of $\boldsymbol{x}$ under the prior density model $p(\boldsymbol{x})\propto\exp(-\phi(\boldsymbol{x}))$. To see why, recall that the MAP estimate maximizes the posterior distribution $p(\boldsymbol{x}|\boldsymbol{y})$. That is, $\widehat{\boldsymbol{x}}_{\text{\sf map}}\triangleq\arg\max_{\boldsymbol{x}}p(\boldsymbol{x}|\boldsymbol{y})=\arg\min_{\boldsymbol{x}}\{-\ln p(\boldsymbol{x}|\boldsymbol{y})\}$. Since Bayes’ rule implies that $\ln p(\boldsymbol{x}|\boldsymbol{y})=\ln p(\boldsymbol{y}|\boldsymbol{x})+\ln p(\boldsymbol{x})-\ln p(\boldsymbol{y})$, we have

[TABLE]

Recalling that the first term in (3) (i.e., the “loss” term) was observed to be $-\ln p(\boldsymbol{y}|\boldsymbol{x})$ (plus a constant) under AWGN noise, the second term in (3) must obey $\phi(\boldsymbol{x})=-\ln p(\boldsymbol{x})+\mathrm{const}$. We will find this MAP interpretation useful in the sequel.

It is not easy to design good regularizers $\phi$ for use in (3). They must not only mimic the negative log signal-prior, but also enable tractable optimization. One common approach is to use $\phi(\boldsymbol{x})=\lambda\|\boldsymbol{\Psi x}\|_{1}$ with $\boldsymbol{\Psi}^{\text{{H}}}$ a tight frame (e.g., a wavelet transform) and $\lambda>0$ a tuning parameter [26]. Such regularizers are convex, and the $\ell_{1}$ norm rewards sparsity in the transform outputs $\boldsymbol{\Psi x}$ when used with the quadratic loss. One could go further and use the composite penalty $\phi(\boldsymbol{x})=\sum_{l=1}^{D}\lambda_{l}\|\boldsymbol{\Psi}_{l}\boldsymbol{x}\|_{1}$. Due to the richness of data structure in MRI, especially for dynamic applications, utilizing multiple ($D>1$) linear sparsifying transforms has been shown to improve recovery quality [27], but tuning multiple regularization weights $\{\lambda_{l}\}$ is a non-trivial problem [28].

Particular insight comes from considering the special case of $\boldsymbol{A}=\boldsymbol{I}$, where the measurement vector in (2) reduces to an AWGN-corrupted version of the image $\boldsymbol{x}$, i.e.,

[TABLE]

The problem of recovering $\boldsymbol{x}$ from noisy $\boldsymbol{z}$, known as “denoising,” has been intensely researched for decades. While it is possible to perform denoising by solving a regularized optimization problem of the form (3) with $\boldsymbol{A}=\boldsymbol{I}$, today’s state-of-the-art approaches are either algorithmic (e.g., [16]) or DNN-based (e.g., [17] [18]). This begs an important question: can these state-of-the-art denoisers be leveraged for MRI signal reconstruction, by exploiting the connections between the denoising problem and (3)? As we shall see, this is precisely what the PnP methods do.

IV Plug-and-Play Methods

In this section, we review several approaches to PnP signal reconstruction. What these approaches have in common is that they recover $\boldsymbol{x}$ from measurements $\boldsymbol{y}$ of the form (2) by iteratively calling a sophisticated denoiser within a larger optimization or inference algorithm.

IV-A Prox-based PnP

To start, let us imagine how the optimization in (3) might be solved. Through what is known as “variable splitting,” we could introduce a new variable, $\boldsymbol{v}$, to decouple the regularizer $\phi(\boldsymbol{x})$ from the data fidelity term $\frac{1}{2\sigma^{2}}\|\boldsymbol{y}-\boldsymbol{A}\boldsymbol{x}\|_{2}^{2}$. The variables $\boldsymbol{x}$ and $\boldsymbol{v}$ could then be equated using an external constraint, leading to the constrained minimization problem

[TABLE]

Equation (6) suggests an algorithmic solution that alternates between separately estimating $\boldsymbol{x}$ and estimating $\boldsymbol{v}$, with an additional mechanism to asymptotically enforce the constraint $\boldsymbol{x}=\boldsymbol{v}$.

The original PnP method [15] is based on the alternating direction method of multipliers (ADMM) [29]. For ADMM, (6) is first reformulated as the “augmented Lagrangian”

[TABLE]

where $\boldsymbol{\lambda}$ are Lagrange multipliers and $\eta>0$ is a penalty parameter that affects the convergence speed of the algorithm, but not the final solution. With $\boldsymbol{u}\triangleq\eta\boldsymbol{\lambda}$, (7) can be rewritten as

[TABLE]

ADMM solves (8) by alternating the optimization of $\boldsymbol{x}$ and $\boldsymbol{v}$ with gradient ascent of $\boldsymbol{u}$, i.e.,

[TABLE]

where $\boldsymbol{h}(\boldsymbol{z};\eta)$ and $\operatorname{prox}_{\phi}(\boldsymbol{z};\eta)$, known as “proximal maps” (see the tutorial [30]) are defined as

[TABLE]

Under some weak technical constraints, it can be proven [29] that when $\phi$ is convex, the ADMM iteration (9) converges to $\widehat{\boldsymbol{x}}$, the global minimum of (3) and (6). If $\boldsymbol{x}_{0}$ is a fixed point of the recursion (9), then the initialization $\boldsymbol{v}_{0}=\boldsymbol{x}_{0}$ and $\boldsymbol{u}_{0}=\frac{\eta}{\sigma^{2}}\boldsymbol{A}^{\text{{H}}}(\boldsymbol{y}-\boldsymbol{Ax}_{0})$ will guarantee that $\boldsymbol{x}_{k}=\boldsymbol{x}_{0}$ for all $k>0$.

From the discussion in Section III, we immediately recognize $\operatorname{prox}_{\phi}(\boldsymbol{z};\eta)$ in (10) as the MAP denoiser of $\boldsymbol{z}$ under AWGN variance $\eta$ and signal prior $p(\boldsymbol{x})\propto\exp(-\phi(\boldsymbol{x}))$. The key idea behind the original PnP work [15] was, in the ADMM recursion (9), to “plug in” a powerful image denoising algorithm like “block-matching and 3D filtering” (BM3D) [16] in place of the proximal denoiser $\operatorname{prox}_{\phi}(\boldsymbol{z};\eta)$ from (10). If the plug-in denoiser is denoted by “$\boldsymbol{f}$,” then the PnP ADMM algorithm becomes

[TABLE]

A wide variety of empirical results (see, e.g., [15, 31, 32] [33]) have demonstrated that, when $\boldsymbol{f}$ is a powerful denoising algorithm like BM3D, the PnP algorithm (13) produces far better recoveries than the regularization-based approach (9). For parallel MRI, the advantages of PnP ADMM were demonstrated in [34].Although the value of $\eta$ does not change the fixed point of the standard ADMM algorithm (9), it affects the fixed point of the PnP ADMM algorithm (13) through the ratio $\sigma^{2}/\eta$ in (12).

The success of PnP methods raises important theoretical questions. Since $\boldsymbol{f}$ is not in general the proximal map of any regularizer $\phi$, the iterations (13) may not minimize a cost function of the form in (3), and (13) may not be an implementation of ADMM. It is then unclear if the iterations (13) will converge. And if they do converge, it is unclear what they converge to. The consensus equilibrium framework, which we discuss in Section V, aims to provide answers to these questions.

The use of a generic denoiser in place of a proximal denoiser can be translated to non-ADMM algorithms, such as FISTA  [35], primal-dual splitting (PDS)  [36], and others, as in [37, 38, 39]. Instead of optimizing $\boldsymbol{x}$ as in (13), PnP FISTA [37] uses the iterative update

[TABLE]

where (14a) is a gradient descent (GD) step on the negative log-likelihood $\frac{1}{2\sigma^{2}}\|\boldsymbol{y}-\boldsymbol{Ax}\|^{2}$ at $\boldsymbol{x}=\boldsymbol{s}_{k-1}$ with step-size \eta{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in(0,\sigma^{2}\|\boldsymbol{A}\|_{2}^{-2})}, (14b) is the plug-in replacement of the usual proximal denoising step in FISTA, and (14c) is an acceleration step, where it is typical to use $q_{k}=(1+\sqrt{1+4q_{k-1}^{2}})/2$ and $q_{0}=1$. If $\boldsymbol{x}_{0}$ is a fixed point of the recursion (14), then the initialization $\boldsymbol{s}_{0}=\boldsymbol{x}_{0}$ will guarantee that $\boldsymbol{x}_{k}=\boldsymbol{x}_{0}$ for all $k>0$.

PnP PDS [38] uses the iterative update

[TABLE]

where $\eta>0$ is a stepsize, $\gamma\in(0,1)$ is a relaxation parameter chosen such that $\gamma\leq\eta/(\eta+\sigma^{2}\|\boldsymbol{A}\|_{2}^{-2})$, and $\boldsymbol{f}(\boldsymbol{z})$ in (15a) is the plug-in replacement of the usual proximal denoiser $\operatorname{prox}_{\phi}(\boldsymbol{z};\eta)$. If $\boldsymbol{x}_{0}$ is a fixed point of the recursion (15), then the initialization $\boldsymbol{v}_{0}=\boldsymbol{Ax}_{0}-\boldsymbol{y}$ will guarantee that $\boldsymbol{x}_{k}=\boldsymbol{x}_{0}$ for all $k>0$.

Comparing PnP ADMM (13) to PnP FISTA (14) and PnP PDS (15), one can see that they differ in how the data fidelity term $\frac{1}{2\sigma^{2}}\|\boldsymbol{y}-\boldsymbol{Ax}\|^{2}$ is handled: PnP ADMM uses the proximal update (12), while PnP FISTA and PnP PDS use the GD step (14a). In most cases, solving the proximal update (12) is much more computationally costly than taking a GD step (14a). Thus, with ADMM, it is common to approximate the proximal update (12) using, e.g., several iterations of the conjugate gradient (CG) algorithm or GD, which should reduce the per-iteration complexity of (13) but may increase the number of iterations. But even with these approximations of (12), PnP ADMM is usually close to “convergence” after 10-50 iterations (e.g., see Figure 4).

An important difference between the aforementioned flavors of PnP is that the stepsize $\eta$ is constrained in FISTA but not in ADMM or PDS. Thus, PnP FISTA restricts the range of reachable fixed points relative to PnP ADMM and PnP PDS.

IV-B The balanced FISTA approach

In Section III, when discussing the optimization problem (3), the regularizer $\phi(\boldsymbol{x})=\lambda\|\boldsymbol{\Psi x}\|_{1}$ was mentioned as a popular option, where $\boldsymbol{\Psi}$ is often a wavelet transform. The resulting optimization problem,

[TABLE]

is said to be stated in “analysis” form [2]. The proximal denoiser associated with (16) has the form

[TABLE]

When $\boldsymbol{\Psi}$ is orthogonal, it is well known that $\operatorname{prox}_{\phi}(\boldsymbol{z};\eta)=\boldsymbol{f}_{\text{\sf tdt}}(\boldsymbol{z};\lambda\eta)$, where

[TABLE]

is the “transform-domain thresholding” denoiser with [\textrm{soft-thresh}(\boldsymbol{u},\tau)]_{n}\triangleq\max\big{\{}0,\frac{|u_{n}|-\tau}{|u_{n}|}\big{\}}u_{n}. The denoiser (18) is very efficient to implement, since it amounts to little more than computing forward and reverse transforms.

In practice, (16) yields much better results with non-orthogonal $\boldsymbol{\Psi}$, such as when $\boldsymbol{\Psi}^{\text{{H}}}$ is a tight frame (see, e.g., the references in [40]). In the latter case, $\boldsymbol{\Psi}^{\text{{H}}}\boldsymbol{\Psi}=\boldsymbol{I}$ with tall $\boldsymbol{\Psi}$. But, for general tight frames $\boldsymbol{\Psi}^{\text{{H}}}$, the proximal denoiser (17) has no closed-form solution. What if we simply plugged the transform-domain thresholding denoiser (18) into an algorithm like ADMM or FISTA? How can we interpret the resulting approach? Interestingly, as we describe below, if (18) is used in PnP FISTA, then it does solve a convex optimization problem, although one with a different form than (3). This approach was independently proposed in [26] and [40], where in the latter it was referred to as “balanced FISTA” (bFISTA) and applied to parallel cardiac MRI. Notably, bFISTA was proposed before the advent of PnP FISTA. More details are provided below.

The optimization problem (16) can be stated in constrained “synthesis” form as

[TABLE]

where $\boldsymbol{\alpha}$ are transform coefficients. Introducing a new parameter $\beta$, the problem (19) can be extended to a 1-parameter family of unconstrained problems

[TABLE]

with projection matrix $\boldsymbol{P}_{\boldsymbol{\Psi}}^{\perp}\triangleq\boldsymbol{I}-\boldsymbol{\Psi\Psi}^{\text{{H}}}$, and with (19) corresponding to the limiting case of (20) with $\beta\rightarrow\infty$. In practice, it is not possible to take $\beta\rightarrow\infty$ and, for finite values of $\beta$, the problems (19) and (20) are not equivalent. However, problem (20) under finite $\beta$ is interesting to consider in its own right, and it is sometimes referred to as the “balanced” approach  [41]. If we solve (20) using FISTA with step-size $\eta>0$ (recall (14a)) and choose the particular value $\beta=1/\eta$ then, remarkably, the resulting algorithm takes the form of PnP FISTA (14) with $\boldsymbol{f}(\boldsymbol{z})=\boldsymbol{f}_{\text{\sf tdt}}(\boldsymbol{z};\lambda)$. This particular choice of $\beta$ is motivated by computational efficiency (since it leads to the use of $\boldsymbol{f}_{\text{\sf tdt}}$) rather than recovery performance. Still, as we demonstrate in Section VI, it performs relatively well.

IV-C Regularization by denoising

Another PnP approach, proposed by Romano, Elad, and Milanfar in [42], recovers $\boldsymbol{x}$ from measurements $\boldsymbol{y}$ in (2) by finding the $\widehat{\boldsymbol{x}}$ that solves the optimality condition111We begin our discussion of RED by focusing on the real-valued case, as in [42] and [43], but later we extend the RED algorithms to the complex-valued case of interest in MRI.

[TABLE]

where $\boldsymbol{f}$ is an arbitrary (i.e., “plug in”) image denoiser and $\eta>0$ is a tuning parameter. In [42], several algorithms were proposed to solve (21). Numerical experiments in [42] suggest that, when $\boldsymbol{f}$ is a sophisticated denoiser (like BM3D) and $\eta$ is well tuned, the solutions $\widehat{\boldsymbol{x}}$ to (21) are state-of-the-art, similar to those of PnP ADMM.

The approach (21) was coined “regularization by denoising” (RED) in [42] because, under certain conditions, the $\widehat{\boldsymbol{x}}$ that solve (21) are the solutions to the regularized least-squares problem

[TABLE]

where the regularizer $\phi_{\text{\sf red}}$ is explicitly constructed from the plug-in denoiser $\boldsymbol{f}$. But what are these conditions? Assuming that $\boldsymbol{f}$ is differentiable almost everywhere, it was shown in [43] that the solutions of (21) correspond to those of (22) when i) $\boldsymbol{f}$ is locally homogeneous222Locally homogeneous means that (1+\epsilon)\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{f}\big{(}(1+\epsilon)\boldsymbol{x}\big{)} for all $\boldsymbol{x}$ and sufficiently small nonzero $\epsilon$. and ii) $\boldsymbol{f}$ has a symmetric Jacobian matrix (i.e., $[J\boldsymbol{f}(\boldsymbol{x})]^{\text{{T}}}=J\boldsymbol{f}(\boldsymbol{x})~{}\forall\boldsymbol{x}$). But it was demonstrated in [43] that these properties are not satisfied by popular image denoisers, such as the median filter, transform-domain thresholding, NLM, BM3D, TNRD, and DnCNN. Furthermore, it was proven in [43] that if the Jacobian of $\boldsymbol{f}$ is non-symmetric, then there does not exist any regularizer $\phi$ under which the solutions of (21) minimize a regularized loss of the form in (3).

One may then wonder how to justify (21). In [43], Reehorst and Schniter proposed an explanation for (21) based on “score matching”  [44], which we now summarize. Suppose we are given a large corpus of training images $\{\boldsymbol{x}_{t}\}_{t=1}^{T}$, from which we could build the empirical prior model

[TABLE]

where $\delta$ denotes the Dirac delta. Since images are known to exist outside $\{\boldsymbol{x}_{t}\}_{t=1}^{T}$, it is possible to build an improved prior model $\widetilde{p_{\text{\sf x}}}$ using kernel density estimation (KDE), i.e.,

[TABLE]

where $\eta>0$ is a tuning parameter. If we adopt $\widetilde{p_{\text{\sf x}}}$ as the prior model for $\boldsymbol{x}$, then the MAP estimate of $\boldsymbol{x}$ (recall (4)) becomes

[TABLE]

Because $\ln\widetilde{p_{\text{\sf x}}}$ is differentiable, the solutions to (24) must obey

[TABLE]

A classical result known as “Tweedie’s formula”  [45, 46] says that

[TABLE]

where $\boldsymbol{f}_{\text{\sf mmse}}(\cdot;\eta)$ is the minimum mean-squared error (MMSE) denoiser under the prior $\boldsymbol{x}\sim\widehat{p_{\text{\sf x}}}$ and $\eta$-variance AWGN. That is, $\boldsymbol{f}_{\text{\sf mmse}}(\boldsymbol{z})=\operatorname{E}\{\boldsymbol{x}\,|\,\boldsymbol{z}\}$, where $\boldsymbol{z}=\boldsymbol{x}+\mathcal{N}(\boldsymbol{0},\eta\boldsymbol{I})$ and $\boldsymbol{x}\sim\widehat{p_{\text{\sf x}}}$. Applying (26) to (25), the MAP estimate $\widehat{\boldsymbol{x}}$ under the KDE prior $\widetilde{p_{\text{\sf x}}}$ obeys

[TABLE]

which matches the RED condition (21) when $\boldsymbol{f}=\boldsymbol{f}_{\text{\sf mmse}}(\cdot;\eta)$. Thus, if we could implement the MMSE denoiser $\boldsymbol{f}_{\text{\sf mmse}}$ for a given training corpus $\{\boldsymbol{x}_{t}\}_{t=1}^{T}$, then RED provides a way to compute the MAP estimate of $\boldsymbol{x}$ under the KDE prior $\widetilde{p_{\text{\sf x}}}$.

Although the MMSE denoiser $\boldsymbol{f}_{\text{\sf mmse}}$ can be expressed in closed form (see [43, eqn. (67)]), it is not practical to implement for large $T$. Thus the question remains: Can the RED approach (21) also be justified for non-MMSE denoisers $\boldsymbol{f}$, especially those that are not locally homogeneous or Jacobian-symmetric? As shown in [43], the answer is yes. Consider a practical denoiser $\boldsymbol{f}_{\boldsymbol{\theta}}$ parameterized by tunable weights $\boldsymbol{\theta}$ (e.g., a DNN). A typical strategy is to choose $\boldsymbol{\theta}$ to minimize the mean-squared error on $\{\boldsymbol{x}_{t}\}_{t=1}^{T}$, i.e., set $\widehat{\boldsymbol{\theta}}=\operatorname*{arg\,min}_{\boldsymbol{\theta}}\operatorname{E}\{\|\boldsymbol{x}-\boldsymbol{f}_{\boldsymbol{\theta}}(\boldsymbol{z})\|^{2}\}$, where the expectation is taken over $\boldsymbol{x}\sim\widehat{p_{\text{\sf x}}}$ and $\boldsymbol{z}=\boldsymbol{x}+\mathcal{N}(\boldsymbol{0},\eta\boldsymbol{I})$. By the MMSE orthogonality principle, we have

[TABLE]

and so we can write

[TABLE]

where (30) follows from (26). Equation (30) says that choosing $\boldsymbol{\theta}$ to minimize the MSE is equivalent to choosing $\boldsymbol{\theta}$ so that $\frac{1}{\eta}(\boldsymbol{f}_{\boldsymbol{\theta}}(\boldsymbol{z})-\boldsymbol{z})$ best matches the “score” $\nabla\ln\widetilde{p_{\text{\sf x}}}(\boldsymbol{z};\eta)$. This is an instance of the “score matching” framework, as described in [44].

In summary, the RED approach (21) approximates the KDE-MAP approach (25)-(27) by using a plug-in denoiser $\boldsymbol{f}$ to approximate the MMSE denoiser $\boldsymbol{f}_{\text{\sf mmse}}$. When $\boldsymbol{f}=\boldsymbol{f}_{\text{\sf mmse}}$, RED exactly implements MAP-KDE, but with a practical $\boldsymbol{f}$, RED implements a score-matching approximation of MAP-KDE. Thus, a more appropriate title for RED might be “score matching by denoising.”

Comparing the RED approach from this section to the prox-based PnP approach from Section IV-A, we see that RED starts with the KDE-based MAP estimation problem (24) and replaces the $\widehat{p_{\text{\sf x}}}$-based MMSE denoiser $\boldsymbol{f}_{\text{\sf mmse}}$ with a plug-in denoiser $\boldsymbol{f}$, while PnP ADMM starts with the $\phi$-based MAP estimation problem (3) and replaces the $\phi$-based MAP denoiser $\operatorname{prox}_{\phi}$ from (10) with a plug-in denoiser $\boldsymbol{f}$. It has recently been demonstrated that, when the prior is constructed from image examples, MAP recovery often leads to sharper, more natural looking image recoveries than MMSE recovery [47] [48, 49, 50]. Thus it is interesting that RED offers an approach to MAP-based recovery using MMSE denoising, which is much easier to implement than MAP denoising [47] [48, 49, 50].

Further insight into the difference between RED and prox-based PnP can be obtained by considering the case of symmetric linear denoisers, i.e., $\boldsymbol{f}(\boldsymbol{z})=\boldsymbol{Wz}$ with $\boldsymbol{W}=\boldsymbol{W}^{\text{{T}}}$, where we will also assume that $\boldsymbol{W}$ is invertible. Although such denoisers are far from state-of-the-art, they are useful for interpretation. It is easy to show [51] that $\boldsymbol{f}(\boldsymbol{z})=\boldsymbol{Wz}$ is the proximal map of $\phi(\boldsymbol{x})=\frac{1}{2\eta}\boldsymbol{x}^{\text{{T}}}(\boldsymbol{W}^{-1}-\boldsymbol{I})\boldsymbol{x}$, i.e., that $\operatorname{prox}_{\phi}(\boldsymbol{z};\eta)=\boldsymbol{Wz}$, recalling (10). With this proximal denoiser, we know that the prox-based PnP algorithm solves the optimization problem

[TABLE]

Meanwhile, since $\boldsymbol{f}(\boldsymbol{z})=\boldsymbol{Wz}$ is both locally homogeneous and Jacobian-symmetric, we know from (22) that the RED under this $\boldsymbol{f}$ solves the optimization problem

[TABLE]

By comparing (31) and (32), we see a clear difference between RED and prox-based PnP. Using $\boldsymbol{u}\triangleq\boldsymbol{W}^{1/2}\boldsymbol{x}$, it can be seen that we can rewrite (32) as

[TABLE]

which has the same regularizer as (31) but a different loss and additional post-processing. Section V-B compares RED to prox-based PnP from yet another perspective: consensus equilibrium.

So far, we have described RED as solving for $\widehat{\boldsymbol{x}}$ in (21). But how exactly is this accomplished? In the original RED paper [42], three algorithms were proposed to solve (21): GD, inexact ADMM, and a “fixed point” heuristic that was later recognized [43] as a special case of the proximal gradient (PG) algorithm  [52, 53]. Generalizations of PG RED were proposed in [43]. The fastest among them is the accelerated-PG RED algorithm, which uses the iterative update333Although the development of RED up to this point has assumed real-valued quantities for simplicity, the RED equations (34) and (35) may be used with complex-valued quantities.

[TABLE]

where $\boldsymbol{h}$ was defined in (12), line (34b) uses the same acceleration as PnP FISTA (14b), and $L>0$ is a design parameter that can be related to the Lipschitz constant of $\phi_{\text{\sf red}}(\cdot)$ from (22) (see [43, Sec.V-C]). When $L=1$ and $q_{k}=1~{}\forall k$, (34) reduces to the “fixed point” heuristic from [42]. To reduce the implementation complexity of $\boldsymbol{h}$, one could replace (34a) with the GD step

[TABLE]

which achieves a similar complexity reduction as when going from PnP ADMM to PnP FISTA (as discussed in Section IV-A). The result would be an “accelerated GD”  [54] form of RED. Convergence of the RED algorithms will be discussed in Section V-B.

IV-D Approximate Message Passing

The approximate message passing (AMP) algorithm  [55] estimates the signal $\boldsymbol{x}$ using the iteration

[TABLE]

starting with $\boldsymbol{v}_{0}=\boldsymbol{0}$ and some $\boldsymbol{x}_{0}$. In (36), $M$ and $N$ refer to the dimensions of $\boldsymbol{A}\in{\mathbb{R}}^{M\times N}$, $J$ denotes the Jacobian operator, and the denoiser $\boldsymbol{f}_{k}$ may change with iteration $k$. In practice, the trace of the Jacobian in (36a) can be approximated using  [56]

[TABLE]

using small $\epsilon>0$ and random $\boldsymbol{p}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$. The original AMP paper  [55] considered only real-valued quantities and proximal denoising (i.e., $\boldsymbol{f}_{k}(\boldsymbol{z})=\operatorname{prox}_{\phi}(\boldsymbol{z};\eta_{k})$ for some $\eta_{k}>0$); the complex-valued case was first described in [57] and the PnP version of AMP shown in (36) was proposed in  [58][59] under the name “denoising AMP.” This approach was applied to MRI in [60], but the results are difficult to reproduce.

The rightmost term in (36a) is known as the “Onsager correction term,” with roots in statistical physics  [61]. If we remove this term from (36a), the resulting recursion is no more than unaccelerated FISTA (also known as ISTA  [62]), i.e., (14) with $q_{k}=1~{}\forall k$ and a particular stepsize $\eta$. The Onsager term is what gives AMP the special properties discussed below.

When $\boldsymbol{A}$ is i.i.d. Gaussian and the denoiser $\boldsymbol{f}_{k}$ is Lipschitz, the PnP AMP iteration (36) has remarkable properties in the large-system limit (i.e., $M,N\rightarrow\infty$ with fixed $M/N$), as proven in [63] for the real-valued case. First,

[TABLE]

That is, the denoiser input $\boldsymbol{z}_{k}$ behaves like an AWGN corrupted version of the true image $\boldsymbol{x}^{0}$. Moreover, as $N\rightarrow\infty$, the mean-squared of AMP’s estimate $\boldsymbol{x}_{k}$ is exactly predicted by the “state evolution”

[TABLE]

Furthermore, if $\boldsymbol{f}_{k}$ is the MMSE denoiser (i.e., $\boldsymbol{f}_{k}(\boldsymbol{z})=\operatorname{E}\{\boldsymbol{x}\,|\,\boldsymbol{z}=\boldsymbol{x}^{0}+\mathcal{N}(\boldsymbol{0},\eta_{k}\boldsymbol{I})\}$) and the state-evolution (39) has a unique fixed-point, then $\mathcal{E}(\eta_{k})$ converges to the minimum MSE. Or, if $\boldsymbol{f}_{k}$ is the proximal denoiser $\boldsymbol{f}_{k}(\boldsymbol{z})=\operatorname{prox}_{\phi}(\boldsymbol{z};\eta_{k})$ and the state-evolution (39) has a unique fixed-point, then AMP’s $\boldsymbol{x}_{k}$ converges to the MAP solution, even if $\phi$ is non-convex.

The remarkable properties of AMP are only guaranteed to manifest with large i.i.d. Gaussian $\boldsymbol{A}$. For this reason, a “vector AMP” (VAMP) algorithm was proposed  [64] with a rigorous state-evolution and denoiser-input property (38) that hold under the larger class of right-rotationally invariant (RRI)444Matrix $\boldsymbol{A}$ is said to be RRI if its singular-value decomposition $\boldsymbol{A}=\boldsymbol{USV}^{\text{{H}}}$ has Haar $\boldsymbol{V}$, i.e., $\boldsymbol{V}$ distributed uniformly over the group of orthogonal matrices for any $\boldsymbol{U}$ and $\boldsymbol{S}$. Note that i.i.d. Gaussian is a special case of RRI where $\boldsymbol{U}$ is also Haar and the singular values in $\boldsymbol{S}$ have a particular distribution. matrices $\boldsymbol{A}$. Later, a PnP version of VAMP was proposed  [65] and rigorously analyzed under Lipschitz $\boldsymbol{f}_{k}$  [66]. VAMP can be understood as the symmetric variant  [67] of ADMM with adaptive penalty $\eta$  [68].

Both AMP and VAMP assume that the matrix $\boldsymbol{A}$ is a typical realization of a random matrix drawn independently of the signal $\boldsymbol{x}$, so that—with high probability—multiplication-by-$\boldsymbol{A}$ has the effect of randomizing $\boldsymbol{x}$. For example, if $\boldsymbol{A}$ is i.i.d. Gaussian or RRI, then the right singular-vector matrix $\boldsymbol{V}$ is Haar distributed, and consequently $\boldsymbol{V}^{\text{{H}}}\boldsymbol{x}$ will be uniformly distributed on the sphere of radius $\|\boldsymbol{x}\|$ for any $\boldsymbol{x}$, which means that $\boldsymbol{V}^{\text{{H}}}\boldsymbol{x}$ will behave like an i.i.d. Gaussian vector at high dimensions. In MRI, however, $\boldsymbol{A}$ is Fourier structured, and MRI images $\boldsymbol{x}$ are also Fourier structured in that they have much more energy near the origin of k-space, with the result that $\boldsymbol{V}^{\text{{H}}}\boldsymbol{x}$ is far from i.i.d. Gaussian. To address this issue, a plug-and-play VAMP approach was proposed in [69] that recovers the wavelet coefficients of the image rather than the image itself. Importantly, it exploits the fact that—for natural images—the energy of the wavelet coefficients is known to decrease exponentially across the wavelet bands. By exploiting this additional structure, VAMP can properly handle the joint structure of $\boldsymbol{A}$ and $\boldsymbol{x}$ that manifests in MRI. In the denoising step, the estimated wavelet coefficients are first transformed to the spatial domain, then denoised, and then transformed back to the wavelet domain. These transformations are made efficient by the $O(N)$ complexity of the wavelet transform and its inverse. Additional details are provided in [69].

V Understanding PnP through Consensus Equilibrium

The success of the PnP methods in Section IV raises important theoretical questions. For example, in the case of PnP ADMM, if the plug-in denoiser $\boldsymbol{f}$ is not the proximal map of any regularizer $\phi$, then it is not clear what cost function is being minimized (if any) or whether the algorithm will even converge. Similarly, in the case of RED, if the plug-in denoiser $\boldsymbol{f}$ is not the MMSE denoiser $\boldsymbol{f}_{\text{\sf mmse}}$, then RED no longer solves the MAP-KDE problem, and it is not clear what RED does solve, or whether a given RED algorithm will even converge. In this section, we show that many of these questions can be answered through the consensus equilibrium (CE) framework [70, 51, 39, 43] [71]. We start by discussing CE for the PnP approaches from Section IV-A and follow with a discussion of CE for the RED approaches from Section IV-C.

V-A Consensus equilibrium for prox-based PnP

Let us start by considering the PnP ADMM algorithm (13). Rather than viewing (13) as minimizing some cost function, we can view it as seeking a solution $(\widehat{\boldsymbol{x}}_{\text{\sf pnp}},\widehat{\boldsymbol{u}}_{\text{\sf pnp}})$ to

[TABLE]

which, by inspection, must hold when (13) is at a fixed point. Not surprisingly, by setting $\boldsymbol{x}_{k}=\boldsymbol{x}_{k-1}$ in the PnP FISTA algorithm (14), it is straightforward to show that it too seeks a solution to (40). It is easy to show that the PnP PDS algorithm [38] seeks the same solution. With (40), the goal of the prox-based PnP algorithms becomes well defined! The pair (40) reaches a consensus in that the denoiser $\boldsymbol{f}$ and the data fitting operator $\boldsymbol{h}$ agree on the output $\widehat{\boldsymbol{x}}_{\text{\sf pnp}}$. The equilibrium comes from the opposing signs on the correction term $\widehat{\boldsymbol{u}}_{\text{\sf pnp}}$: the data-fitting subtracts it while the denoiser adds it.

Applying the $\boldsymbol{h}$ expression from (12) to (40a), we find that

[TABLE]

where $\boldsymbol{y}-\boldsymbol{A}\widehat{\boldsymbol{x}}_{\text{\sf pnp}}$ is the k-space measurement error and $\boldsymbol{A}^{\text{{H}}}(\boldsymbol{y}-\boldsymbol{A}\widehat{\boldsymbol{x}}_{\text{\sf pnp}})$ is its projection back into the image domain. We now see that $\widehat{\boldsymbol{u}}_{\text{\sf pnp}}$ provides a correction on $\widehat{\boldsymbol{x}}_{\text{\sf pnp}}$ for any components that are inconsistent with the measurements, but it provides no correction for errors in $\widehat{\boldsymbol{x}}_{\text{\sf pnp}}$ that lie outside the row-space of $\boldsymbol{A}$. Plugging (41) back into (40b), we obtain the image-domain fixed-point equation

[TABLE]

If $\boldsymbol{y}=\boldsymbol{Ax}+\boldsymbol{w}$, as in (2), and we define the PnP error as $\widehat{\boldsymbol{e}}_{\text{\sf pnp}}\triangleq\widehat{\boldsymbol{x}}_{\text{\sf pnp}}-\boldsymbol{x}$, then (42) implies

[TABLE]

which says that the error $\widehat{\boldsymbol{e}}_{\text{\sf pnp}}$ combines with the k-space measurement noise $\boldsymbol{w}$ in such a way that the corresponding image-space correction $\widehat{\boldsymbol{u}}_{\text{\sf pnp}}=\frac{\eta}{\sigma^{2}}\boldsymbol{A}^{\text{{H}}}(\boldsymbol{w}-\boldsymbol{A}\widehat{\boldsymbol{e}}_{\text{\sf pnp}})$ is canceled by the denoiser $\boldsymbol{f}(\cdot)$.

By viewing the goal of prox-based PnP as solving the equilibrium problem (40), it becomes clear that other solvers beyond ADMM, FISTA, and PDS can be used. For example, it was shown in [70] that the PnP CE condition (40) can be achieved by finding a fixed-point of the system555The paper [70] actually considers the consensus equilibrium among $N>1$ agents, whereas here we consider the simple case of $N=2$ agents.

[TABLE]

There exist many algorithms to solve (44). For example, one could use the Mann iteration  [30]

[TABLE]

when $\mathcal{F}$ is nonexpansive. The paper [70] also shows that this fixed point is equivalent to the solution of $\mathcal{F}(\underline{\boldsymbol{z}})=\mathcal{G}(\underline{\boldsymbol{z}})$, in which case Newton’s method or other root-finding methods could be applied.

The CE viewpoint also provides a path to proving the convergence of the PnP ADMM algorithm. Sreehari et al. [31] used a classical result from convex analysis to show that a sufficient condition for convergence is that i) $\boldsymbol{f}$ is non-expansive, i.e., $\|\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{y})\|\leq\|\boldsymbol{x}-\boldsymbol{y}\|$ for any $\boldsymbol{x}$ and $\boldsymbol{y}$, and ii) $\boldsymbol{f}(\boldsymbol{x})$ is a sub-gradient of some convex function, i.e., there exists $\varphi$ such that $\boldsymbol{f}(\boldsymbol{x})\in\partial\varphi(\boldsymbol{x})$. If these two conditions are met, then PnP ADMM (13) will converge to a global solution. Similar observations were made in other recent studies, e.g., [51] [72]. That said, Chan et al. [32] showed that many practical denoisers do not satisfy these conditions, and so they designed a variant of PnP ADMM in which $\eta$ is decreased at every iteration. Under appropriate conditions on $\boldsymbol{f}$ and the rate of decrease, this latter method also guarantees convergence, although not exactly to a fixed point of (40) since $\eta$ is no longer fixed.

Similar techniques can be used to prove the convergence of other prox-based PnP algorithms. For example, under certain technical conditions, including non-expansiveness of $\boldsymbol{f}$, it was established [39] that PnP FISTA converges to the same fixed-point as PnP ADMM.

V-B Consensus equilibrium for RED

Just as the prox-based PnP algorithms in Section IV-A can be viewed as seeking the consensus equilibrium of (40), it was shown in [43] that the proximal-gradient and ADMM-based RED algorithms seek the consensus equilibrium $(\widehat{\boldsymbol{x}}_{\text{\sf red}},\widehat{\boldsymbol{u}}_{\text{\sf red}})$ of

[TABLE]

where $\boldsymbol{h}$ was defined in (12) and $L$ is the algorithmic parameter that appears in (34).666The parameter $L$ also manifests in ADMM RED, as discussed in [43]. Since (47) takes the same form as (40), we can directly compare the CE conditions of RED and prox-based PnP.

Perhaps a more intuitive way to compare the CE conditions of RED and prox-based PnP follows from rewriting (47b) as $\widehat{\boldsymbol{x}}_{\text{\sf red}}=\boldsymbol{f}(\widehat{\boldsymbol{x}}_{\text{\sf red}})+L\widehat{\boldsymbol{u}}_{\text{\sf red}}$, after which the RED CE condition becomes

[TABLE]

which involves no inverse operations. In the typical case of $L=1$, we see that (48) matches (40), except that the correction $\widehat{\boldsymbol{u}}_{\text{\sf red}}$ is added after denoising in (48b) and before denoising in (40b).

Yet another way to compare the CE conditions of RED and prox-based PnP is to eliminate the $\widehat{\boldsymbol{u}}_{\text{\sf red}}$ variable. Solving (48a) for $\widehat{\boldsymbol{u}}_{\text{\sf red}}$ gives

[TABLE]

which mirrors the expression for $\widehat{\boldsymbol{u}}_{\text{\sf pnp}}$in (41). Then plugging $\widehat{\boldsymbol{u}}_{\text{\sf red}}$ back into (48b) and rearranging, we obtain the fixed-point equation

[TABLE]

or equivalently

[TABLE]

which says that the data-fitting correction (i.e., the left side of (51)) must balance the denoiser correction (i.e., the right side of (51)).

The CE framework also facilitates the convergence analysis of RED algorithms. For example, using the Mann iteration  [30], it was proven in [43] that when $\boldsymbol{f}$ is nonexpansive and $L>1$, the PG RED algorithm converges to a fixed point.

VI Demonstration of PnP in MRI

VI-A Parallel cardiac MRI

We now demonstrate the application of PnP methods to parallel cardiac MRI. Because the signal $\boldsymbol{x}$ is a cine (i.e., a video) rather than a still image, there are relatively few options available for sophisticated denoisers. Although algorithmic denoisers like BM4D  [73] have been proposed, they tend to run very slowly, especially relative to the linear operators $\boldsymbol{A}$ and $\boldsymbol{A}^{\text{{H}}}$. For this reason, we first trained an application specific CNN denoiser for use in the PnP framework. The architecture of the CNN denoiser, implemented and trained in PyTorch, is shown in Figure 1.

For training, we acquired 50 fully sampled, high-SNR cine datasets from eight healthy volunteers. Thirty three of those were collected on a 3 T scanner777The 3 T scanner was a Magnetom Prisma Fit from Siemens Healthineers in Erlangen, Germany and the 1.5 T scanner was a Magnetom Avanto from Siemens Healthineers in Erlangen, Germany. and the remaining 17 were collected on a 1.5 T scanner. Out of the 50 datasets, 28, 7, 7, and 8 were collected in the short-axis, two-chamber, three-chamber, and four-chamber view, respectively. The spatial and temporal resolutions of the images ranged from 1.8 mm to 2.5 mm and from 34 ms to 52 ms, respectively. The images size ranged from $160\times 130$ to $256\times 208$ pixels and the number of frames ranged from $15$ to $27$. For each of the 50 datasets, the reference image series was estimated as the least-squares solution to (1), with the sensitivity maps $\boldsymbol{S}_{i}$ estimated from the time-averaged data using ESPIRiT  [74]. We added zero-mean, complex-valued i.i.d. Gaussian noise to these “noise-free” reference images to simulate noisy images with SNR of 24 dB. Using a fixed stride of $30\times 30\times 10$, we decomposed the images into patches of size $55\times 55\times 15$. The noise-free and corresponding noisy patches were assigned as output and input to the CNN denoiser, with the real and imaginary parts processed as two separate channels. All 3D convolutions were performed using $3\times 3\times 3$ kernels. There were 64 filters of size $3\times 3\times 3\times 2$ in the first layer, 64 filters of size $3\times 3\times 3\times 64$ in the second through fourth layers, and 2 filters of size $3\times 3\times 3\times 64$ in the last layer. We set the minibatch size to four and used the Adam optimizer  [75] with a learning rate of $1\times 10^{-4}$ over 400 epochs. The training process was completed in 12 hours on a workstation equipped with a single NVIDIA GPU (GeForce RTX 2080 Ti).

For testing, we acquired four fully sampled cine datasets from two different healthy volunteers, with two image series in the short-axis view, one image series in the two-chamber view, and one image series in the four-chamber view. The spatial and temporal resolutions of the images ranged from 1.9 mm to 2 mm and from 37 ms to 45 ms, respectively. For the four datasets, the space-time signal vector, $\boldsymbol{x}$, in (2) had dimensions of $192\times 144\times 25$, $192\times 144\times 25$, $192\times 166\times 16$, and $192\times 166\times 16$, respectively, with the last dimension representing the number of frames. The datasets were retrospectively downsampled at acceleration rates, $R$, of 6, 8, and 10 using pseudo-random sampling [76]. A representative sampling pattern used to undersample one of the datasets is shown in Figure 2. The data were compressed to $C=12$ virtual coils for faster computation  [21]. The measurements were modeled as described in (1), with the sensitivity maps, $\boldsymbol{S}_{i}$, estimated from the time-averaged data using ESPIRiT.

For compressive MRI recovery, we used PnP ADMM from (13) with $\boldsymbol{f}(\cdot)$ as the CNN-based denoiser described above; we will refer to the combination as PnP-CNN. We employed a total of 100 ADMM iterations, and in each ADMM iteration, we performed four steps of CG to approximate (12), for which we used $\sigma^{2}=1=\eta$. (See Figure 5 for the effect of $\sigma^{2}/\eta$ on the final NMSE and the convergence rate.)We compared this PnP method to three CS-based methods: CS with an undecimated wavelet transform (CS-UWT), CS with total variation (TV),888Note that sometimes UWT and TV are combined [1]. and a low-rank plus sparse (L+S) method (see, e.g., [77]). We also compared to PnP-UWT and the transform-learning [14]method LASSI [78].

For PnP-UWT, we used PnP FISTA from (14) with $\boldsymbol{f}(\cdot)$ implemented as $\boldsymbol{f}_{\text{\sf tdt}}$ given in (18), i.e., bFISTA. A three-dimensional single-level Haar UWT was used as $\boldsymbol{\Psi}$ in (18). For CS-TV, we used a 3D finite-difference operator for $\boldsymbol{\Psi}$ in the regularizer $\phi(\boldsymbol{x})=\|\boldsymbol{\Psi x}\|_{1}$, and for CS-UWT, we used the aforementioned UWT instead. For both CS-TV and CS-UWT, we used monotone FISTA [79] to solve the resulting convex optimization problem (3). For L+S, the method by Otazo et al. [77] was used. The regularization weights for CS-UWT, PnP-UWT, CS-TV, and L+S were manually tuned to maximize the reconstruction SNR (rSNR)999rSNR is defined as $\|\boldsymbol{x}\|^{2}/\|\widehat{\boldsymbol{x}}-\boldsymbol{x}\|^{2}$, where $\boldsymbol{x}$ is the true image and $\widehat{\boldsymbol{x}}$ is the estimate. for Dataset #3 at $R=10$. For LASSI we used the authors’ implementation at https://gitlab.com/ravsa19/lassi, and we did our best to manually tune all available parameters.

The rSNR values are summarized in Table I. For all four datasets and three acceleration rates, PnP-CNN exhibited the highest rSNR with respect to the fully sampled reference. Also, compared to the CS methods and PnP-UWT, which uses a more traditional denoiser based on soft-thresholding of UWT coefficients, PnP-CNN was better at preserving anatomical details of the heart; see Figure 3. The performance of PnP-UWT was similar to that of CS-UWT. Figure 4 plots NMSE as a function of the number of iterations for the CS and PnP methods. Since the CS methods were implemented using CPU computation and the PnP methods were implemented using GPU computation, a direct runtime comparison was not possible. We did, however, compare the per-iteration runtime of PnP ADMM for two different denoisers: the CNN and UWT-based $\boldsymbol{f}_{\text{\sf tdt}}$ described earlier in this section. When the CNN denoiser was replaced with the UWT-based $\boldsymbol{f}_{\text{\sf tdt}}$, the per-iteration runtime changed from 2.05 s to 2.1 s, implying that the two approaches have very similar computational costs.

For PnP-CNN, Figure 5 shows the dependence of the final NMSE ($=\text{rSNR}^{-1}$) and of the convergence rate on $\sigma^{2}/\eta$ for one of the testing datasets included in this study. Overall, final NMSE varies less than 0.5 dB for $\sigma^{2}/\eta\in[0.5,2]$ for all four datasets and all three acceleration rates. Figure 6 compares CG and GD when solving (12). To this end, NMSE vs. runtime is plotted for different numbers of CG or GD inner-iterations for Dataset #3 at $R=10$. For GD, the step-size was manually optimized. Figure 6 suggests that it is best to use a $1$ to $4$ inner iterations of either GD or CG; using more inner iterations slows the convergence rate without improving the final performance.

The results in this section, although preliminary, highlight the potential of PnP methods for MRI recovery of cardiac cines. By optimizing the denoiser architecture, the performance of PnP-CNN may be further improved.

VI-B Single-coil fastMRI knee data

In this section, we investigate recovery of 2D knee images from the single-coil fastMRI dataset [80]. This dataset contains fully-sampled k-space data that are partitioned into 34 742 training slices and 7 135 testing slices. The Cartesian sampling patterns from [80] were used to achieve acceleration rate $R=4$.

We evaluated PnP using the ADMM algorithm with a learned DnCNN [17] denoiser. To accommodate complex-valued images, DnCNN was configured with two input and two output channels. The denoiser was then trained using only the central slices of the 3 T scans without fat-suppression from the training set, comprising a total of 267 slices (i.e., $<1\%$ of the total training data). The training-noise variance and the PnP ADMM tuning parameter $\sigma^{2}/\eta$ were manually adjusted in an attempt to maximize rSNR.

PnP was then compared to the TV and U-Net baseline methods described and configured in [80]. For example, 128 channels were used for the U-Net’s first layer, as recommended in [80]. We then trained three versions of the U-Net. The first version was trained on the full fastMRI training set101010The full fastMRI training set includes 1.5 T and 3 T scans, with and without fat suppression, and an average of 36 slices per volume. with random sampling masks. The second U-Net was trained on the full fastMRI training set, but with a fixed sampling mask. The third U-Net was trained using only the central slices of the 3 T scans without fat-suppression (i.e., the same data used to train the DnCNN denoiser) and with a fixed sampling mask.

To evaluate performance, we used the central slices of the non-fat-suppressed 3 T scans from the validation set, comprising a total of 49 slices. The evaluation considered both random sampling masks and the same fixed mask used for training. The resulting average rSNR and SSIM scores are summarized in Table II. The table shows that PnP-CNN performed similarly to the U-Nets and significantly better than TV. In particular, PnP-CNN achieved the highest rSNR score with both random and fixed testing masks, and the U-Net gave slightly higher SSIM scores in both tests. Among the U-Nets, the version trained with a fixed sampling mask and full data gave the best rSNR and SSIM performance when testing with the same mask, but its performance dropped considerable when testing with random masks. Meanwhile, the U-Net trained with the smaller data performed significantly worse than the other U-Nets, with either fixed or random testing masks. And although this latter U-Net used exactly the same training data as the PnP-CNN method, it was not competitive with PnP-CNN. Although preliminary, these results suggest that i) PnP methods are much less sensitive to deviations in the forward model between training and testing, and that ii) PnP methods are effective with relatively small training datasets.

VII Conclusion

PnP methods present an attractive avenue for compressive MRI recovery. In contrast to traditional CS methods, PnP methods can exploit richer image structure by using state-of-the-art denoisers. To demonstrate the potential of such methods for MRI reconstruction, we used PnP to recover cardiac cines and knee images from highly undersampled datasets. With application-specific CNN-based denoisers, PnP was able to significantly outperform traditional CS methods and to perform on par with modern deep-learning methods, but with considerably less training data. The time is ripe to investigate the potential of PnP methods for a variety of MRI applications.