Plug-and-Play Deep Energy Model for Inverse problems

TL;DR

This paper introduces an energy-based Plug-and-Play framework for image recovery that learns a CNN-based prior with convergence guarantees, enabling more complex priors and improved MRI reconstruction performance.

Contribution

It presents the first energy-based PnP formulation using a CNN score function with convergence guarantees, allowing for more expressive priors and better image reconstruction.

Findings

Improved MRI reconstruction performance over traditional PnP methods.

Energy-based formulation ensures convergence even without contraction constraints.

CNN score function models complex image priors effectively.

Abstract

We introduce a novel energy formulation for Plug- and-Play (PnP) image recovery. Traditional PnP methods that use a convolutional neural network (CNN) do not have an energy based formulation. The primary focus of this work is to introduce an energy-based PnP formulation, which relies on a CNN that learns the log of the image prior from training data. The score function is evaluated as the gradient of the energy model, which resembles a UNET with shared encoder and decoder weights. The proposed score function is thus constrained to a conservative vector field, which is the key difference with classical PnP models. The energy-based formulation offers algorithms with convergence guarantees, even when the learned score model is not a contraction. The relaxation of the contraction constraint allows the proposed model to learn more complex priors, thus offering improved performance over traditional PnP schemes. Our experiments in magnetic resonance image reconstruction demonstrates the improved performance offered by the proposed energy model over traditional PnP methods.

Tables

Table 1. Table I: Comparison of the proposed Energy based PnP method for two different acceleration in the presence of Gaussian noise of std = 0.01 absent 0.01 =0.01 .

Algorithm	Avg. PSNR $\pm$ std
	Acceleration $6$x	Acceleration $2$x
EPnP GD	37.67 $\pm$1.19	41.47 $\pm$ 1.18
RED	34.93 $\pm$ 1.71	40.95 $\pm$ 1.41
SCORE	34.95 $\pm$ 1.61	40.43 $\pm$ 1.58

Equations

$$\mathbf{A}\mathbf{x}+\mathbf{n}=\mathbf{b}$$

$${\mathbf{x}}^{\ast }=\arg \underset{\mathbf{x}}{\min }\underset{-\log p(\mathbf{b}|\mathbf{x})}{\underbrace{\frac{1}{2{\eta }^2}\Vert \mathbf{A}\mathbf{x}-\mathbf{b}{\Vert }_2^2}}+\underset{-\log p_{\mathrm{data}}(\mathbf{x})}{\underbrace{\phi (\mathbf{x})}}.$$

$$p_{\theta }(\mathbf{x})=\frac{1}{{\mathbf{Z}}_{\theta }}\exp \left(-\frac{{\Psi }_{\theta }(\mathbf{x})}{{\sigma }^2}\right)$$

$$f(\mathbf{x})=\frac{1}{2{\eta }^2}\Vert \mathbf{A}\mathbf{x}-\mathbf{b}{\Vert }_2^2+\frac{1}{{\sigma }^2}{\Psi }_{\theta }(\mathbf{x})$$

$${\nabla }_{\mathbf{x}}\left(-\log p_{\theta }(\mathbf{x})\right)=\frac{{\nabla }_{\mathbf{x}}{\Psi }_{\theta }(\mathbf{x})}{{\sigma }^2}=\frac{\mathbf{H}(\mathbf{x})}{{\sigma }^2}$$

$${\theta }^{\ast }=\arg \underset{\theta }{\min }{\mathbb{E}}_{\mathbf{x}\sim p_{\mathrm{data}}}\Vert F_{\theta }(\mathbf{x})-{\nabla }_{\mathbf{x}}\log p_{\mathrm{data}}(\mathbf{x}){\Vert }^2.$$

$${\theta }^{\ast }$$

$$\log p_{\theta }(\mathbf{x})={\int }_C{F}_{\theta }(\mathbf{z})d\mathbf{z}+\log p_{\theta }(\mathbf{a})$$

$$\begin{array}[]{ll}{\boldsymbol{x}}_{n+1}={\boldsymbol{x}}_{n}-\gamma\underbrace{\left(\dfrac{{\mathbf{A}}^{T}({\mathbf{A}}{\boldsymbol{x}}_{n}-{\boldsymbol{b}})}{\eta^{2}}+\dfrac{H({\boldsymbol{x}}_{n})}{\sigma^{2}}\right)}_{\nabla f(\mathbf{x}_{n})}.\end{array}$$

$$L_{\mathrm{eq}}=\frac{1}{{\eta }^2}+\frac{L}{{\sigma }^2}.$$

Full text

I Introduction

The recovery of an image from a few corrupted measurements is a common problem in several applications including denoising, deblurring, and MRI image reconstruction. The maximum a posteriori (MAP) framework formulates image recovery as the minimization of a cost function that is the sum of two terms. The first term is a measure of the consistency of the image with the measurements, while the second regularization term depends on the available prior information about the image. Classical approaches rely on handcrafted priors, while recent methods use deep learned plug and play (PnP) priors that offer improved performance over traditional approaches [1, 2, 3].

While the empirical performance of PnP methods is superior to compressed sensing and low-rank methods, a challenge with current PnP methods is the lack of a well-defined energy formulation, typically enjoyed by classical methods. While some approaches including the RED [4] introduced an energy formulation, subsequent analysis showed that many of the assumptions in RED are not satisfied by CNN based denoisers [5]; most PnP methods are now understood from a consensus equilibrium perspective [3]. Another limitation is the requirement that the denoiser is a contraction, which is needed for the algorithm to converge. This constraint is often enforced by spectral normalization of the individual layers [6]. In addition to restricting the possible CNN architectures, our experiments show that this constraint often restricts the CNN from learning accurate prior models, which translates to lower performance.

The main focus of this paper is to introduce a novel plug and play energy formulation to overcome the above challenges. In particular, we use a convolution neural network (CNN) to represent the negative log density of the images. We compute the gradient of the above CNN, which models the gradient of the negative log-likelihood of the data, which is often termed as the score. The score model is pre-trained using noise score matching [7]. The main distinction of this approach from traditional PnP models is that score function is constrained to be a conservative vector field, which satisfies the property that its line integral is independent of the specific path. The learned energy function can be used in arbitrary inverse problems by combining it with the data likelihood term to obtain the posterior distribution, which is minimized using gradient descent.

We also introduce a novel convergence guarantee, which shows that the cost will decrease monotonically when the step-size of the algorithm is sufficiently low. The convergence guarantees are valid even when the score function is not a contraction. Traditional PnP models that learn non-conservative score functions impose an Lipschitz bound on the score function to guarantee convergence. However, this constraint may restrict the type of prior densities that can be learned. Because our network is not required to be a contraction, it can learn more complex energy functions, which translate to improved performance. We determine the utility of the proposed scheme in the recovery of MRI data from highly undersampled measurements.

II Proposed Method

II-A Problem Formulation

Let us consider the problem of recovering an image ${\boldsymbol{x}}\in\mathbb{R}^{m}$ from its corrupted linear measurements ${\boldsymbol{b}}\in\mathbb{R}^{n}$:

[TABLE]

where ${\mathbf{A}}\in\mathbb{R}^{n\times m}(n\leq m)$ is a known linear transformation and ${\boldsymbol{n}}\in\mathbb{N}(0,\eta^{2}{\mathbf{I}})$ is the additive white Gaussian noise. Then the MAP estimate of $\mathbf{x}$ is given by :

[TABLE]

II-B PnP energy model and the resulting score function

In this paper, we use a parametric model to represent $p_{\rm data}({\boldsymbol{x}})$:

[TABLE]

where $\theta$ denote the parameters of the model and ${\boldsymbol{Z}}_{\theta}$ is a normalizing constant. We model $\Psi_{\theta}({\boldsymbol{x}}):\mathbb{R}^{n}\rightarrow\mathbb{R}$ by a neural network. Here, $\sigma^{2}>0$ is the training noise variance, described later. Using the above prior in (2), the MAP objective simplifies to :

[TABLE]

We note that the gradient of the log prior $\nabla_{{\boldsymbol{x}}}\left({-\log p_{\theta}({\boldsymbol{x}})}\right):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$:

[TABLE]

is often referred to as the score function of the prior. For example, assume that $\Psi_{\theta}({\boldsymbol{x}})=\eta(\mathbf{W}{\boldsymbol{x}})$, where $\mathbf{W}$ is a $\mathbb{R}^{m\times m}$ vector and $\eta$ is non-linear operator, applied on each entry. Then $\nabla_{{\boldsymbol{x}}}\Psi_{\theta}({\boldsymbol{x}})$ is specified by $\mathbf{W}^{T}\eta^{\prime}(\mathbf{W}{\boldsymbol{x}})$, which is a two layer neural network with shared weights. Fig. 1 shows two different example models. From Fig.1 it is observed that when the convolutional layers use larger strides, the architecture resembles a UNET with skip connections. The output of the encoder is a scalar that yields the energy $\Psi_{\theta}({\boldsymbol{x}})$, while the output of the decoder is the score. However, unlike a regular UNET, the network of the encoder and decoder are shared. In particular, the network is constrained such that $\nabla{\mathcal{E}}_{\theta}$ is a conservative vector field, which ensures that the resulting score network can be interpreted as the derivative of a well-defined energy function ${\mathcal{E}}_{\theta}({\boldsymbol{x}})$.

II-C Learning the energy using denoising score matching

Score matching [7] aims to learn a parametric model $F_{\theta}({\boldsymbol{x}})=\nabla_{x}\log p_{\theta}({\boldsymbol{x}})$ from training data $\mathbf{x}\sim p_{\rm data}({\boldsymbol{x}})$ such that it matches the score of the proposed distribution $\nabla_{{\boldsymbol{x}}}\log p({\boldsymbol{x}})$

[TABLE]

The above cost function is essentially the Fisher divergence between $p_{\rm data}({\boldsymbol{x}})$ and $p_{\theta}(x)$. Because this cost is challenging to evaluate without the explicit knowledge of $p_{\rm data}({\boldsymbol{x}})$, the DSM approach instead learns the score of $\tilde{{\boldsymbol{x}}}={\boldsymbol{x}}+\mathbf{n}$, where $\mathbf{n}\sim\mathcal{N}(0,\sigma^{2}{\mathbf{I}})$. Note that $\tilde{{\boldsymbol{x}}}\sim p_{\rm data}({\boldsymbol{x}})\ast\mathcal{N}(0,\sigma^{2}{\mathbf{I}})$, which is a smoothed version of $p_{\rm data}({\boldsymbol{x}})$ [7]:

[TABLE]

In this approach, zero-mean Gaussian white noise of variance $\sigma$ is added to each training data point ${\boldsymbol{x}}$ to yield its corresponding perturbed sample $\tilde{{\boldsymbol{x}}}={\boldsymbol{x}}+\sigma{\boldsymbol{z}}$ where ${\boldsymbol{z}}\in\mathcal{N}(0,{\mathbf{I}})$. As $\sigma\rightarrow 0$, this approach minimizes the problem in (6).

We set $F_{\theta}=-\frac{\mathbf{H}_{\theta}({\boldsymbol{x}})}{\sigma^{2}}$ to obtain

[TABLE]

While the approach in (7) is very similar to classical DSM [7], the main distinction is the representation of $F_{\theta}$ as the gradient vector field of $\frac{{\Psi}_{\theta}({\boldsymbol{x}})}{\sigma^{2}}$. Constraining the score function as a conservative vector field allows us to use the Fundamental theorem of line integrals [8] to represent the energy at any given point $x$ as the line integral:

[TABLE]

where $\mathbf{a}$ is an arbitrary point and $C$ is any curve between ${\boldsymbol{a}}$ and ${\boldsymbol{x}}$. Because the score function is conservative, the line integral in (8) is independent of the path taken. We note that this property is not valid for arbitrary residual denoisers used in classical PnP models.

II-D Optimization Algorithm

In this work, we assume the Lipschitz constant of the score network $\Psi_{\theta}$ is upper-bounded by $L$, which may be greater than one. We use the CLIP approach [9] to estimate the approximate constant. One may also use the product of the spectral norms to obtain an upperbound for $L$. We propose to minimize the cost function in (4) using the steepest descent algorithm:

[TABLE]

We note that the fixed point of the algorithm satisfies $\nabla f(x)=0$, which correspond to the minimum of (4).

Theorem II.1**.**

Let the maximum eigen-value of $\mathbf{A}^{T}\mathbf{A}$ is one. Then thegradient descent algorithm specified by (9) will monotonically converge to a minimum of (4), when the step-size $\gamma\leq 1/L_{\rm eq}$, where:

[TABLE]

The proof is omitted due to space constraints. Note that the Lipschitz constant of $\nabla f(x)$ is bounded by (10). The result is based on the well-known result that a Lipschitz constrained function can be upper-bounded by a quadratic function; the quadratic term depends on the Lipschitz constant. The pseudo code for the algorithm is shown below. We compute $L$ after training, which allows us to choose a step-size that can guarantee convergence.

We note that constraining the network to learn a conservative score function ensures that the gradient magnitude decreases as one approaches the minimum, thus guaranteeing convergence. We note that the above convergence result is only applicable to energy based models; the gradient vector field learned by classical PnP methods is not constrained to be conservative. Hence, it may learn high magnitude gradients, even close to the minima. Classical methods constrain the Lipschitz constant of the score function to guarantee convergence. However, this approach often translates to networks with lower performance.

[TABLE]

III Experiments and Results

III-A Experimental setup

In this paper, we evaluate the performance of the proposed method in the context of recovering parallel MRI data from highly undersampled measurements. Here, the matrix ${\mathbf{A}}$ in (1) will be equal to ${\mathbf{A}}={\boldsymbol{S}}{\mathbf{F}}{\boldsymbol{C}}$, where ${\boldsymbol{S}}$ is the sampling matrix, ${\mathbf{F}}$ is Fourier Transform, and ${\boldsymbol{C}}$ is the coil sensitivity map. We use the publicly available dataset [10], which consists of images with a matrix size of $256\times 232$. The training and test dataset consists of $360$ and $160$ slices, respectively. The dataset was obtained using a 12-channel head coil and CSM was pre-computed using ESPIRIT algorithm. We evaluated the proposed method using variable-density Cartesian random sampling mask with two different undersampling factors and in the presence of Gaussian noise with standard deviation $\eta=0.01$.

III-B Architecture of the networks

We implemented the function ${\boldsymbol{\Psi}}_{\theta}({\boldsymbol{x}})$ in (3) using a five layer convolutional neural network that consists of $64$ channels per layer and $3\times 3$ filters, followed by a linear layer. After each convolutional layer, ReLU activations were used. As discussed in Section II.A., the function ${\boldsymbol{H}}({\boldsymbol{x}})$ was realized using the chain rule, which guarantees the learned score to be a conservative vector field. The architecture of an example network is shown in Fig. 1. We note that the weights of the encoder and the decoder are shared. We train ${\boldsymbol{H}}({\boldsymbol{x}})$ as a noise estimator as in (7) with $\sigma=0.01$. The Lipschitz constant of the denoiser was estimated using the technique described in [9] and was used to determine the step size $\gamma$ of the steepest-descent algorithm in (9). The algorithm was run until $\dfrac{|f_{\rm{MAP}}({\boldsymbol{x}}_{n+1})-f_{\rm{MAP}}({\boldsymbol{x}}_{n})|}{|f_{\rm{MAP}}({\boldsymbol{x}}_{n})|}\leq 10^{-8}$ was satisfied. The proposed algorithm is refered to as EPnP-GD.

The proposed algorithm is compared with the steepest descent based RED algorithm (SD-RED) algorithm [4] and a score-based model. We use a ten layer CNN with 64 channels per layers with 3x3 filters. The number of network parameters is approximately the same as the energy network used above. We use spectral normalization [6] to ensure that the networks are contractions. All the networks were trained on the training datasets described above.

III-C Results

Following training, the performance of the algorithms was evaluated on the test dataset. Table. I shows the results for two different acceleration rates in the presence of Gaussian noise with standard deviation of $0.01$. For fair comparisons, all algorithms where initialized with ${\mathbf{A}}^{T}{\boldsymbol{b}}$. From Table I, we observe that EPnP GD offers improved performance compared to RED and SCORE at both accelerations. We note that RED and traditional SCORE models require the Lipschitz constant of the networks to be bounded by one, while the EPnP-GD approach does not require this constraint. We attribute the improved performance to the relaxation of the Lipschitz constraint. The improved results can also be appreciated from the reconstructions shown in Fig. 2. We note that EPnP offers reconstructions with reduced noise artifacts compared to RED and SCORE models.

We note that the proposed model is associated with a well-defined cost function unlike score models. We show the convergence of the cost function in (4) for different acceleration and slices in Fig. 3. We show the convergence plot of the EPnP-GD algorithm in Fig. 3 for ten different test slices wherein From the figure it can be observed that the EPnP-GD decreases the cost function monotonically.

IV Conclusion

In this paper, we introduced a novel energy formulation for PnP framework. The proposed framework represented the log prior of the data using a deep CNN, where the gradient of the prior modeled the score function. The resulting score network resembles a U-NET whose encoder and decoder weights are shared. The parameters of the network are learned using denoising score matching. A steepest descent algorithm was used to apply the learned model to inverse problems. The proposed algorithm is guaranteed to converge monotonically to a minimum of the MAP objective, even when when the score function is not a contraction. The preliminary results in the context of MRI shows that the relaxation of the contraction constraint translates to improved performance.