Extra Proximal-Gradient Network with Learned Regularization for Image Compressive Sensing Reconstruction

Abstract

Learned optimization algorithms are promising approaches to inverse problems by leveraging advanced numerical optimization schemes and deep neural network techniques in machine learning. In this paper, we propose a novel deep neural network architecture imitating an extra proximal gradient algorithm to solve a general class of inverse problems with a focus on applications in image reconstruction. The proposed network features learned regularization that incorporates adaptive sparsification mappings, robust shrinkage selections, and nonlocal operators to improve solution quality. Numerical results demonstrate the improved efficiency and accuracy of the proposed network over several state-of-the-art methods on a variety of test problems.

1. Introduction

Recent years have witnessed the substantial success of deep neural networks (DNN) in a large variety of real-world applications [1,2,3,4,5,6,7,8,9]. Equipped with proven expressive power, DNNs can be used to approximate highly complicated functions provided a sufficient amount of data [10]. However, training DNNs as end-to-end black-boxes can be extremely data demanding, rendering DNNs difficult to interpret, generalize, and sensitive to noise and outliers. To overcome these issues, learned optimization algorithms (LOAs) have started to gain attention as they are designed to combine the interpretable mechanism of optimization algorithms and the expressive power of DNNs. One of the most important applications of LOAs is solving the inverse problem of general form

$$\underset{\mathit{x}}{min}f(\mathit{x};\mathit{y})+g\left(\mathit{x}\right),$$

(1)

where f is the data fidelity term determined by the data formation and noise distribution that relate the target solution $\mathit{x}$ and the given measurement data $\mathit{y}$, and g is the critical (possibly nonconvex) regularization term that promotes the desired solution $\mathit{x}$, as f is often underdetermined and the data $\mathit{y}$ can be incomplete and noisy. In classical approaches to inverse problems, the regularization g is often handcrafted based on human heuristics and limited experience, which can be overly simplified and not capable to capture the intrinsic complex features of the solution. LOAs, on the other hand, allow the regularization g to be learned from training data and hence can result in significant improvement over the handcrafted regularizations.

Our goal in this paper is to propose an efficient extra proximal gradient algorithm that employs the Nesterov’s acceleration technique and the extra gradient scheme, and unroll this algorithm into a deep neural network called the extra proximal gradient network (EPGN) to solve a class of inverse problems (1). Motivated by the least absolute shrinkage and selection operator (LASSO) [11,12,13], our EPGN implicitly adopts an $l_1$-type regularization in (1) with a nonlinear sparsification mapping learned from data. The proximal operator of this regularization is elaborated by several linear convolutions, nonlinear activation functions, and shrinkage operations for robust sparse feature selection in EPGN. As our focus application is in image reconstruction, we also incorporate a nonlocal feature selection component into the learned regularization to leverage similar patterns within images and improve reconstruction quality. The proposed EPGN combines the advantages of the accelerated extra gradient scheme, the sparsity promoting nonlinear transforms, and the nonlocal feature selections. As a consequence, our EPGN is efficient, robust, and accurate in a variety of image reconstruction problems as demonstrated by the numerical experiments.

2. Related Work

One of the early LOAs is the learned iterative shrinkage thresholding algorithm (LISTA) for solving $l_1$ regularized linear inversion [14]. LISTA maps the standard ISTA optimization algorithm to a recurrent neural network (RNN) with certain layer weights learned from training data to improve the performance. The asymptotic linear convergence rate for LISTA is established in [15,16]. Several variations of LISTA are proposed for image reconstruction with regularizations based on low rank or group sparsity [17], $l_0$ minimization [18], and learned approximate message passing [19]. These LOA methods employ handcrafted regularizations and require a closed-form solution of the proximal operator of the regularization term. The idea of LISTA is also extended to solve composite problems with linear constraints, called differentiable linearized alternating direction method of multipliers (D-LADMM) [20], which exhibits an asymptotic linear convergence rate.

To learn more general and adaptive regularization function in (1), the other group of LOAs is proposed to solve inverse problem (1) with learnable regularization. A straightforward approach in this group uses deep convolutional neural network (CNN), denoted by ${\mathit{h}}_k(·)$, to replace the proximal operator ${prox}_{{\alpha }_{k}g}$ [21] of the unknown regularization term g in the proximal gradient update:

$$\begin{array}{c}\hfill {\mathit{x}}_{k+1}={prox}_{{\alpha }_{k}g}\left({\mathit{b}}_k\right),\end{array}$$

(2)

where ${\alpha }_k>0$ is the step size in the kth iteration, ${\mathit{b}}_k:={\mathit{x}}_k-{\alpha }_k\nabla f({\mathit{x}}_k;\mathit{y})$, and ${prox}_g(·)$ is defined by

$${prox}_g\left(\mathit{b}\right):=\underset{\mathit{x}}{argmin}\frac{1}{2}{\parallel \mathit{x}-\mathit{b}\parallel }^2+g\left(\mathit{x}\right).$$

(3)

Therefore, one avoids explicit formation of the regularization g, but creates a neural network with prescribed K phases, where each phase mimics one iteration of the proximal gradient method such as (2) to compute ${\mathit{b}}_k$ as above and ${\mathit{x}}_k={\mathit{h}}_k\left({\mathit{b}}_k\right)$. The CNN ${\mathit{h}}_k$ can also be cast as a residual network (ResNet) [22] to represent the discrepancy between ${\mathit{b}}_k$ and the improved ${\mathit{x}}_k$ [23]. Such a paradigm is also embedded into half quadratic splitting [23], ADMM [24], and primal dual methods [25] to replace the proximal operator in the subproblems. To improve over the generic black-box CNNs above, several LOA methods are proposed to unroll numerical optimization algorithms such as deep neural networks so as to preserve their efficient structures with proven efficiency, such as the ADMM-Net [26] and ISTA-Net [27]. These methods also prescribe the phase number K and map each iteration of the corresponding numerical algorithm to one phase of the network, and learn specific components of the network using training data.

3. Extra Proximal Gradient Network

In this section, we propose a novel parameter-efficient deep neural network architecture to solve the inverse problem (1) with regularization learned from data. To this end, we first introduce the accelerated extra proximal gradient algorithm that combines Nesterov’s acceleration technique and the extra proximal gradient update in Section 3.1. In Section 3.2, we mimic this algorithm to construct the proposed extra proximal gradient network (EPGN), where Nesterov’s acceleration step corresponds to a simple linear combination layer in EPGN to boost convergence, and the extra proximal gradient structure induces a predictor–corrector update scheme with efficient utilization of network parameters in EPGN. For the image reconstruction applications considered in our experiment part, we integrate mixing layers into EPGN to combine the local and nonlocal image features for enhanced reconstruction quality in Section 3.3. Additional details of the EPGN training process are provided in Section 3.4.

3.1. Extra Proximal Gradient Algorithm

The extra gradient method proposed in the seminal work [28] has attracted significant interest in optimization in recent years. It has been extended to solve variational inequality problems [29] and convex/nonconvex composite optimization problems [30] with theoretical performance guaranteed. Extra gradient algorithms use an additional gradient step in a first-order optimization algorithm to improve the convergence results. This can also be interpreted as a predictor–corrector scheme to speed up convergence. The following two variants of the original extra gradient algorithm are closely related to our proposed extra proximal gradient algorithm. The first one is the extended extra gradient method in [30], which uses extra proximal gradient steps at each iteration to solve nonconvex composite minimization problem (1) by

$$\begin{array}{cc}\hfill {\mathit{x}}_{k+\frac{1}{2}}& ={prox}_{{\alpha }_{k}g}({\mathit{x}}_k-{\alpha }_k\nabla f({\mathit{x}}_k;\mathit{y})),\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill {\mathit{x}}_{k+1}& ={prox}_{{\beta }_{k}g}({\mathit{x}}_k-{\beta }_k\nabla f({\mathit{x}}_{k+\frac{1}{2}};\mathit{y})).\hfill \end{array}$$

(4b)

The second one is the convex accelerated extra gradient algorithm developed in [31], which integrates Nesterov’s accelerated gradient method for smooth convex optimization [32] into the extra gradient scheme. Different from the classical extra gradient method, this algorithm evaluates gradients in both steps at an interpolation of the previous two iterates rather than the previous iterate only. Recall that Nesterov’s acceleration technique [32] for minimizing smooth convex function f is given by

$$\begin{array}{cc}\hfill {\tilde{\mathit{x}}}_k& ={\mathit{x}}_k+{\gamma }_k({\mathit{x}}_k-{\mathit{x}}_{k-1}),\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill {\mathit{x}}_{k+1}& ={\tilde{\mathit{x}}}_k-\alpha \nabla f({\tilde{\mathit{x}}}_k;\mathit{y}),\hfill \end{array}$$

(5b)

which performs a momentum structure (5a) to improve the convergence rate of standard gradient methods. For nonconvex problems, a monitor mechanism that tunes ${\gamma }_k$ adaptively can be introduced to remedy convergence issue [33]. Motivated by this acceleration technique, we propose to combine (4) and (5) and introduce the accelerated extra proximal gradient updating scheme summarized in Algorithm 1 to solve inverse problems of form (1). In Algorithm 1, ${\alpha }_k$ and ${\beta }_k$ are step sizes, and ${\gamma }_k$ is the momentum coefficient in the kth iteration.

Algorithm 1: Accelerated Extra Proximal Gradient Algorithm.

Input: Data $\mathit{y}$ and initialization ${\mathit{x}}_0={\mathit{x}}_{-\frac{1}{2}}$. Output:$\mathit{x}={\mathit{x}}_K$. For $k=0,1,2,\cdots ,K-1$, do $$\begin{array}{}\hfill \mathrm{(6a)}& {\tilde{\mathit{x}}}_k={\mathit{x}}_k+{\gamma }_k({\mathit{x}}_k-{\mathit{x}}_{k-\frac{1}{2}}),\hfill \mathrm{(6b)}& {\mathit{b}}_{k+\frac{1}{2}}={\tilde{\mathit{x}}}_k-{\alpha }_k\nabla f({\tilde{\mathit{x}}}_k;\mathit{y}),\hfill \mathrm{(6c)}& {\mathit{x}}_{k+\frac{1}{2}}={prox}_{{\alpha }_{k}g}\left({\mathit{b}}_{k+\frac{1}{2}}\right),\hfill \mathrm{(6d)}& {\widehat{\mathit{x}}}_k={\mathit{x}}_{k+\frac{1}{2}}+{\gamma }_k({\mathit{x}}_{k+\frac{1}{2}}-{\mathit{x}}_k),\hfill \mathrm{(6e)}& {\mathit{b}}_{k+1}={\widehat{\mathit{x}}}_k-{\beta }_k\nabla f({\widehat{\mathit{x}}}_k;\mathit{y}),\hfill \mathrm{(6f)}& {\mathit{x}}_{k+1}={prox}_{{\beta }_{k}g}\left({\mathit{b}}_{k+1}\right).\end{array}$$

3.2. Extra Proximal Gradient Network (EPGN)

We now cast Algorithm 1 as an LOA by mapping its iterations to the phases of a deep neural network. To this end, we select a phase number K (value to be specified in our experiment), and construct a deep neural network with K phases where each phase performs the updates described in (6). More specifically, we retain the same updates (6a), (6b), (6d) and (6e) in the kth phase of the network. As a result, (6a) and (6d) are simple linear combination layers to integrate the momentum term for acceleration, and (6b) and (6e) are gradient updates for improved fitting to the data. The parameters ${\alpha }_k,{\beta }_k,{\gamma }_k$ are all to be learned for every phase k (we set ${\alpha }_k={\beta }_k$ for simplicity). The remaining updates (6c) and (6f) are replaced by a robust implicit ResNet-type update (we will make it explicitly computable later):

$$\begin{array}{cc}\hfill {\mathit{x}}_{k+l}& ={\mathit{b}}_{k+l}+{\mathit{r}}_k\left({\mathit{x}}_{k+l}\right),\hfill \end{array}$$

(7)

where $l=1/2,1$, and the residual mapping ${\mathit{r}}_k$ plays a critical role of regularization that improves the quality of output ${\mathit{x}}_{k+l}$ in each phase k. In this paper, we parameterize ${\mathit{r}}_k$ as a composition of two nonlinear mappings, denoted by ${\mathcal{G}}_k$ and ${\tilde{\mathcal{G}}}_k$, such that:

$$\begin{array}{cc}\hfill {\mathit{r}}_k\left({\mathit{x}}_{k+l}\right)& ={\tilde{\mathcal{G}}}_k\circ {\mathcal{G}}_k\left({\mathit{x}}_{k+l}\right).\hfill \end{array}$$

(8)

In the remainder of this subsection, we show the details of the CNN structures of these two nonlinear mappings ${\mathcal{G}}_k$ and ${\tilde{\mathcal{G}}}_k$ and how to make the implicit residual update (7) explicit by leveraging the robust shrinkage selection operator.

3.2.1. Nonlinear Feature Extraction Operator ${\mathcal{G}}_k$

We parametrize the nonlinear operator ${\mathcal{G}}_k$ as a multilayer convolutional network of the following structure:

$${\mathcal{G}}_k\left(\mathit{x}\right)={\mathit{B}}_k\sigma \left({\mathit{A}}_k{\mathit{D}}_k\mathit{x}\right),$$

(9)

where ${\mathit{D}}_k$ and ${\mathit{A}}_k$ are two linear convolutional operations that generate and convolve the local features of the input $\mathit{x}$, $\sigma$ is a nonlinear activation function set to the rectified linear unit (ReLU) (i.e., $\sigma \left(\mathit{x}\right)=max(\mathit{x},0)$), and ${\mathit{B}}_k$ is another linear convolution that fuses the activated local features. All the linear mappings ${\mathit{A}}_k$, ${\mathit{B}}_k$, and ${\mathit{D}}_k$ are realized as $3\times 3$ convolutions. Hence, the size of the receptive field (RF) [34] of ${\mathcal{G}}_k$ is $7\times 7$.

The purpose of ${\mathcal{G}}_k$ is to extract the main features of its input, such that these features can be easily refined by a robust feature selector. To this end, we employ the soft shrinkage selection operator in LASSO, which is the proximal operator of the $l_1$ norm and proved to be effective in selecting sparse outstanding features and suppressing noises of its input. More specifically, we consider ${\mathcal{G}}_k\left(\mathit{b}\right)$ as the features prepared to be further pruned by shrinkage (as $\mathit{b}$ is obtained by direct gradient update (6b) and (6e) which may contain undesired artifacts, and hence the features ${\mathcal{G}}_k\left(\mathit{b}\right)$ need further refinement), and ${\mathcal{G}}_k\left(\mathit{x}\right)$ as the refined feature obtained by pruning ${\mathcal{G}}_k\left(\mathit{b}\right)$ using shrinkage. In other words, we expect ${\mathcal{G}}_k\left({\mathit{x}}_{k+l}\right)$ to be

$${\mathcal{S}}_k\left({\mathcal{G}}_k\left({\mathit{b}}_{k+l}\right)\right)={prox}_{{\theta }_k{\parallel ·\parallel }_1}\left({\mathcal{G}}_k\left({\mathit{b}}_{k+l}\right)\right),$$

(10)

where the shrinkage threshold ${\theta }_k>0$ is also to be learned with ${\mathit{A}}_k$, ${\mathit{B}}_k$, and ${\mathit{D}}_k$. Note that the component-wise shrinkage operator ${\mathcal{S}}_k$ in (10) has a closed form solution as ${\left[{\mathcal{S}}_k\left(\mathit{z}\right)\right]}_i=max\left(\right|z_i|-{\theta }_k,0)·z_i/\left|z_i\right|$ for each component $z_i$ of $\mathit{z}$. To further increase our network capacity, we set the convolution ${\mathit{B}}_k$ in ${\mathcal{G}}_k$ to contain $N_f$ kernels ($N_f$ is set to 32 by default), each of size $3\times 3$ in our implementation, hence ${\mathcal{G}}_k\left(\mathit{x}\right)$ has $N_f$ channels at each pixel of $\mathit{x}$. The shrinkage operator ${\mathcal{S}}_k$ in (10) is applied channel-wise with varying ${\theta }_{k,j}$ where $j=1,\cdots ,N_f$. Hence, the learnable parameters of ${\mathcal{G}}_k$ include one convolution ${\mathit{D}}_k$ with $N_f$ kernels of size $3\times 3$ and convolutions ${\mathit{A}}_k$ and ${\mathit{B}}_k$ with $N_f$ kernels of size $3\times 3\times N_f$, and those of ${\mathcal{S}}_k$ are the shrinkage thresholds ${\mathit{\theta }}_k=\{{\theta }_{k,j}:j\in \left[N_f\right]\}$.

3.2.2. Nonlinear Residual Resembling Operator ${\tilde{\mathcal{G}}}_k$

Based on (8), the purpose of the nonlinear operator $\tilde{\mathcal{G}}$ is to resemble the residual term using the refined feature ${\mathcal{G}}_k\left(\mathit{x}\right)$, we can interpret ${\mathcal{G}}_k$ and ${\tilde{\mathcal{G}}}_k$ respectively as encoder and decoder in a symmetric form. More specifically, we parametrize ${\tilde{\mathcal{G}}}_k\left(\mathit{x}\right)$ as ${\tilde{\mathit{D}}}_k{\tilde{\mathit{A}}}_k\sigma \left({\tilde{\mathit{B}}}_k\mathit{x}\right)$, where ${\tilde{\mathit{A}}}_k$, ${\tilde{\mathit{B}}}_k$, and ${\tilde{\mathit{D}}}_k$ are all $3\times 3$ convolutional operators, and ${\tilde{\mathit{D}}}_k$ compresses the $N_f$ channels back to 1 channel according to ${\mathit{D}}_k$ in implementation.

Combining the parametrized nonlinear operators ${\mathcal{G}}_k$ and ${\tilde{\mathcal{G}}}_k$ and the shrinkage operator ${\mathcal{S}}_k$ into (8), we obtain an explicit update rule of (7) given by

$$\begin{array}{c}\hfill {\mathit{x}}_{k+l}={\mathit{b}}_{k+l}+{\tilde{\mathcal{G}}}_k\circ {\mathcal{S}}_k\circ {\mathcal{G}}_k\left({\mathit{b}}_{k+l}\right).\end{array}$$

(11)

This update rule is employed in (6c) and (6f) for $l=1/2,1$ respectively in the kth phase.

To summarize, our proposed extra proximal gradient network (EPGN) of a prescribed K phases is constructed by unrolling Algorithm 1, where each phase executes (6) but with (6c) and (6f) substituted by (11). The flowchart of the computation in EPGN is shown in Figure 1. The proposed EPGN not only inherits the advantages of Algorithm 1 but also employs learnable feature selection operations. Hence, EPGN combines the following properties: (i) the simple linear momentum layers (6a) and (6d) for improved convergence; (ii) the extra proximal gradient updates (6b), (6e), and (11) mimic the predictor–corrector scheme for parameter-efficient network structure; (iii) learnable feature extraction, selection, and residual resembling operators $({\mathcal{G}}_k,{\tilde{\mathcal{G}}}_k,{\mathcal{S}}_k)$ in (11) to effectively improve solution quality as the input data flows through EPGN.

3.3. EPGN with Nonlocal Operator (NL-EPGN)

Nonlocal methods have proven effective for image reconstruction problems, such as in variational methods [35] and nonvariational approaches such as the notable BM3D algorithm [36]. Nonlocal operators can significantly improve image quality as they use image patches located in different regions to exploit the inherent self-similarity of images. Recently, the success of the nonlocal methods has motivated the investigation of the architecture of DNNs that have the ability to capture long-distance dependencies of the image. The deep network architecture for gray-scale and color image denoising in [37] is inspired by the projected gradient algorithm for solving a common variational image restoration model with a learnable nonlocal regularization. The nonlocal neural network proposed in [38] can be viewed as a generalization of the classical nonlocal mean in [35] that computes the response at a position as a weighted average of the image intensities at all positions. The weights implicitly depend on the feature maps in the patches with the size determined by the receptive fields.

To exploit repeated features in images for enhanced reconstruction quality, we adopt the idea in [38]. However, unlike [38] which only relies on nonlocal features, our NL-EPGN fuses local and nonlocal features of images using a combination operator learned through training data. More specifically, we propose NL-EPGN to integrate a nonlocal operator ${\mathcal{N}}_k$ into the residual operator in (11), so that the features refined by the shrinkage operation ${\mathcal{S}}_k$ can be passed to ${\mathcal{N}}_k$ to leverage nonlocal features in images:

$$\begin{array}{cc}\hfill {\mathit{x}}_{k+l}& ={\mathit{b}}_{k+l}+{\tilde{\mathcal{G}}}_k\circ {\mathcal{N}}_k\circ {\mathcal{S}}_k\circ {\mathcal{G}}_k\left({\mathit{b}}_{k+l}\right).\hfill \end{array}$$

(12)

The operator ${\mathcal{N}}_k$ contains two main components: a nonlocal block ${\mathcal{M}}_k$ that extracts nonlocal features of the input, and a nonlinear layer that combines the local and nonlocal features. The details of these two components are given as follows.

3.3.1. Nonlocal Feature Extraction Block ${\mathcal{M}}_k$

Our design of the nonlocal feature extraction block ${\mathcal{M}}_k\left({\mathit{x}}_k\right)$ follows the work [38] which computes a weighted average of features at all locations in an image. More precisely, let ${\left[\mathit{z}\right]}_j$ denote the input feature vector at position j and ${\left[\mathit{v}\right]}_i$ the response vector at position i of an image, then the nonlocal block ${\mathcal{M}}_k$ computes ${\left[\mathit{v}\right]}_i$ by:

$${\left[\mathit{v}\right]}_i={\sum }_j{w}_{ij}{\left[\phi \left(\mathit{z}\right)\right]}_j,$$

(13)

where the function $\phi$ computes a representation of the input signal at position j, and $w_{ij}$ is the normalized weight depending on the similarity between ${\left[\mathit{z}\right]}_i$ and ${\left[\mathit{z}\right]}_j$. The mapping $\phi$ corresponds to a learnable matrix ${\mathit{W}}^{\phi }$ (implemented as $1\times 1$ convolution). The weights are computed by embedded Gaussian:

$$w_{ij}=\frac{exp\left({\left[{\mathit{W}}^{\alpha }\mathit{z}\right]}_i^{\top }{\left[{\mathit{W}}^{\beta }\mathit{z}\right]}_j\right)}{{\sum }_{j}exp\left({\left[{\mathit{W}}^{\alpha }\mathit{z}\right]}_i^{\top }{\left[{\mathit{W}}^{\beta }\mathit{z}\right]}_j\right)},$$

(14)

where both ${\mathit{W}}^{\alpha }$ and ${\mathit{W}}^{\beta }$ are implemented as $N_f/2$ convolutional filters of kernel size $1\times 1$. We employ the bottleneck structure to reduce computation [38]. Hence, the nonlocal block ${\mathcal{M}}_k$ in phase k is implemented as $\mathit{v}={\mathcal{M}}_k\left(\mathit{z}\right)$ where

$${\mathcal{M}}_k\left(\mathit{z}\right)=\mathrm{softmax}\left({\left[{\mathit{W}}_k^{\alpha }\mathit{z}\right]}^{\top }\left[{\mathit{W}}_k^{\beta }\mathit{z}\right]\right){\mathit{W}}_k^{\phi }\left(\mathit{z}\right).$$

(15)

3.3.2. Local and Nonlocal Combination Layer

We propose to use a learnable combination layer of form $\sigma \left({\mathcal{C}}_k[\mathit{z},\mathit{v}]\right)$ to merge the input local feature $\mathit{z}$ and nonlocal feature $\mathit{v}$ obtained by nonlocal block ${\mathcal{M}}_k$ in (15). That is, the nonlocal operator ${\mathcal{N}}_k$ is defined by

$${\mathcal{N}}_k\left(\mathit{z}\right)=\sigma \left({\mathcal{C}}_k[\mathit{z},{\mathcal{M}}_k\left(\mathit{z}\right)]\right).$$

(16)

In the kth phase, the inputs of ${\mathcal{N}}_k$ are $\mathit{z}={\mathcal{S}}_k\circ {\mathcal{G}}_k\left({\mathit{b}}_{k+l}\right)$ for $l=1/2,1$ as shown in (12), $[·,·]$ stands for the concatenation operator at each pixel, and ${\mathcal{C}}_k$ corresponds to a set of learnable weight vectors which project the concatenated vector to a scalar at each pixel (implemented as $1\times 1$ convolution). The flowchart of the nonlocal operator is shown in Figure 2.

3.4. Network Training

As discussed above, the proposed EPGN (or NL-EPGN) consists of K phases, where each phase imitates one iteration in the accelerated extra gradient Algorithm 1 with proximal steps (6b) and (6e) replaced by (11) (or (12) for NL-EPGN). The flowchart of variables in the kth phase is shown in Figure 3. The parameters in the kth phase are collectively denoted by ${\Theta }_k$, which includes the feature extraction operator ${\mathcal{G}}_k=[{\mathit{A}}_k,{\mathit{B}}_k,{\mathit{D}}_k]$, the residual resembling operator ${\tilde{\mathcal{G}}}_k=[{\tilde{\mathit{A}}}_k,{\tilde{\mathit{B}}}_k,{\tilde{\mathit{D}}}_k]$, the nonlocal operator $[{\mathit{W}}_k^{\alpha },{\mathit{W}}_k^{\beta },{\mathit{W}}_k^{\phi },{\mathcal{C}}_k]$, the momentum coefficient ${\gamma }_k$, and the shrinkage thresholds ${\mathit{\theta }}_k$. Let $\Theta =\{{\Theta }_k:0\le k\le K-1\}$ be the set of all network parameters. Then, given N training data pairs of form $\{({\mathit{x}}^{\left(i\right)},{\mathit{y}}^{\left(i\right)}):1\le i\le N\}$, where ${\mathit{y}}^{\left(i\right)}$ is the input measurement data and ${\mathit{x}}^{\left(i\right)}$ is the corresponding ground truth of the ith pair, we define the loss function of $\Theta$ as

$$L(\Theta )=\frac{1}{N}\sum _{i=1}^N{\parallel {\mathit{x}}_K({\mathit{y}}^{\left(i\right)};\Theta )-{\mathit{x}}^{\left(i\right)}\parallel }_2^2,$$

(17)

where ${\mathit{x}}_K(\mathit{y};\Theta )$ denotes the output of the EPGN (i.e., the output of the last, Kth phase) parametrized by $\Theta$ given input data $\mathit{y}$. The optimal network parameter ${\Theta }^{*}$ is obtained by minimizing the loss function (17) in the training process. After training, the EPGN with parameter ${\Theta }^{*}$ serves as a feed-forward neural network that can reconstruct high-quality image $\mathit{x}$ given new measurement data $\mathit{y}$ on the fly.

4. Numerical Experiments

In this section, we evaluate the performance of the proposed EPGN and NL-EPGN on several inverse problems in imaging reconstruction applications. We focus on the reconstruction problem in compressive sensing in our experiments, however, the proposed method can be easily adapted to other image reconstruction problems by changing the data-fidelity term accordingly. All the experiments are implemented, trained, and tested in the TensorFlow framework [39] on a desktop with an Nvidia GTX-1080Ti GPU and 11 GB of graphics card memory (NVIDIA Corporation, Santa Clara, CA, USA).In all tests, the network parameters $\Theta$ of EPGN/NL-EPGN are initialized using the Xavier method [40] and trained with the Adam optimizer [41] with learning rate $1\times {10}^{-4}$ for 200 epochs. To evaluate the reconstruction quality, we use the average peak signal-to-noise ratio (PSNR).

4.1. Nature Images Compressive Sensing

We first test EPGN on the compressive sensing (CS) image reconstruction problem. In our experiment we use the 91 Images dataset for training and Set11 for testing [42]. For a fair comparison, we follow the same data preparation and result evaluation procedures in [27]. The ground truth data $\{{\mathit{x}}^{\left(i\right)}:1\le i\le N\}$ contains $N=88,912$ image patches with luminance components that are all randomly cropped into size $33\times 33$ from 91 Images dataset. We then generate a matrix with random Gaussian entries of size $10\%n$ and $25\%n$, where $n={33}^2$, and orthogonalize the rows. Then the measurement data for training is $\{{\mathit{y}}^{\left(i\right)}=\mathrm{\Psi }{\mathit{x}}^{\left(i\right)}:1\le i\le N\}$. The testing data Set11 preparation follows the same procedure as training data.

4.1.1. Comparison with Existing Methods

We set the phase number $K=9$ for EPGN and $K=7$ for NL-EPGN in this test (as shown in Figure 4 where the PSNRs of the networks become saturated).

Table 1. Natural image CS reconstruction results by existing methods and the proposed EPGN (with 9 phases) and NL-EPGE (with 7 phases) on dataset Set11 with CS ratios of 10% and 25%. Table shows the average PSNR (dB) of the comparison methods.

Method	10%		25%
	PSNR	SSIM	PSNR	SSIM
TVAL3 [43]	22.99	0.3758	27.92	0.6238
D-AMP [44]	22.64	-	28.46	-
IRCNN [23]	24.02	-	30.07	-
ReconNet [42]	24.28	0.6406	25.60	0.7589
DR${}^2$-Net [45]	24.32	0.7175	28.66	0.8432
ISTA-Net+ [27]	26.64	0.8036	32.57	0.9237
DPA-Net [46]	26.99	0.8354	31.74	0.9238
EPGN (9-phase)	27.12	0.8893	32.87	0.9611
NL-EPGN (7-phase)	27.33	0.8956	33.02	0.9623

shows the comparison of the average PSNRs of the images reconstructed by EPGN/NL-EPGN versus several state-of-the-art image reconstruction methods, namely TVAL3 [43], D-AMP [44], IRCNN [23], ReconNet [42], DR²-Net [45], ISTA-Net⁺ [27], and DPA-Net [46], where the first two are classical optimization-based methods, and the last five are deep learning-based methods. The PSNR results of the first four methods and ISTA-Net⁺ in Table 1 are quoted from [27]. We observe that EPGN and NL-EPGN outperform all aforementioned algorithms, whereas NL-EPGN obtains the highest accuracy.

4.1.2. Reconstruction Quality Assessment

Compared to the state-of-the-art ISTA-Net⁺ [27], both EPGN and NL-EPGN obtain better reconstruction results with a similar number of parameters as shown in

Table 2. Compressive sensing reconstruction performance comparison of the 9-phase ISTA-Net${}^{+}$, 15-phase ISTA-Net${}^{+}$, 9-phase EPGN, and 7-phase NL-EPGN on Set11 with a CS ratio of $25\%$ on the number of network parameters (# PARM), average PSNR in dB with standard deviation over the reconstructed images, and average reconstruction time (Time) of one image in second.

Network (# Phase)	# PARM	PSNR (dB)	Time (s)
ISTA-Net${}^{+}$ (9)	336,978	32.57 ± 2.20	0.084
ISTA-Net${}^{+}$ (15)	561,630	32.60 ± 2.19	0.103
EPGN (9)	337,275	32.87 ± 2.24	0.110
NL-EPGN (7)	290,997	33.02 ± 2.05	0.802

. In particular, Figure 5 shows the reconstructed butterfly image with a CS ratio of 10%, from which we can see that the 9-phase EPGN and 7-phase NL-EPGN can both capture the inconspicuous detail of the butterfly wings at the lower left part in the zoomed-in images. Similarly, Figure 6 shows the reconstructed cameraman image with a CS ratio of 25%, where the 9-phase EPGN and 7-phase NL-EPGN have fewer artifacts in the background compared to ISTA-Net⁺, as observed in the lower right area of the zoomed-in images. Figure 7 presents the reconstruction results of the Barbara image in Set11 from ISTA-Net⁺,the 9-phase EPGN, and the 7-phase NL-EPGN with a CS ratio of 10%. We can observe that the texture pattern of the scarf is better preserved by NL-EPGN due to the nonlocal operator.

4.1.3. Parameter Efficiency

The number of network parameters in each phase of ISTA-Net${}^{+}$ is $37,442$ [27]. The number of trainable parameters of each phase in EPGN is $\{{\mathcal{G}}_k+{\tilde{\mathcal{G}}}_k+{\gamma }_k+{\alpha }_k+{\mathit{\theta }}_k=32\times 3\times 3\times (1+32\times 2)+32\times 3\times 3\times (32\times 2+1)+1+2+32=37,475\}$. Similarly, the number of learnable parameters of each phase in EPGN is $41,571$. Therefore, the number of network parameters in each phase of ISTA-Net${}^{+}$, EPGN, and NL-EPGN are very similar (NL-EPGN is about $10.9\%$ more than ISTA-Net${}^{+}$ and EPGN). Figure 4 shows the reconstruction PSNR of these three methods versus phase number, from which we observe that NL-EPGN becomes saturated with phase number $K\ge 7$, whereas EPGN with $K\ge 9$ and ISTA-Net+ with $K\ge 11$. Nevertheless, as shown in Table 2, a 7-phase NL-EPGN has fewer network parameters than a 9-phase EPGN but achieves even higher PSNR.

We compare the reconstruction results of EPGN and ISTA-Net${}^{+}$ in a range of different phase numbers with a CS ratio of $25\%$, as shown in Figure 4. We observe that the PSNR values improve as the phase number increases and become saturated after $K\ge 9$. EPGN achieves a 0.3 dB higher PSNR on average than ISTA-Net${}^{+}$. To further demonstrate the superiority of the extra proximal-gradient method over extending network depth, we compare the 9-phase EPGN with the 15-phase ISTA-Net${}^{+}$, as shown in Table 2. Compared to the 15-phase ISTA-Net${}^{+}$ which extends the depth of the network by simply adding more phases, the 9-phase EPGN achieves better accuracy (0.27 dB higher) using much fewer parameters and similar reconstruction time. We compare the reconstruction performance of EPGN and NL-EPGN with CS ratios of 10% and 25%, the 7-phase NL-EPGN outperforms the 9-phase EPGN by 0.21 dB and 0.15 dB respectively, as shown in Table 1. We also compare the reconstruction results of NL-EPGN and EPGN in a range of different phase numbers with a CS ratio of 25%. The results are shown in Figure 4. We observe that NL-EPGN achieves an average of 0.2 dB PSNR better than EPGN. It is interesting that the PSNR of NL-EPGN shows no significant improvement after $K=7$. As shown in Table 2, the reconstruction time of the 7-phase NL-EPGN is approximate 7 to 8 times that of the 9-phase EPGN due to the time complexity of the nonlocal operator. However, the effect of the nonlocal operator is remarkable, NL-EPGN with 7 phases has a 0.15 dB PSNR improvement with 13.7% fewer parameters compared to EPGN with 9 phases. In Figure 8, we show the PSNR versus epoch using the proposed NL-EPGN and the state-of-the-art method ISTA-Net${}^{+}$ for image reconstruction with a CS ratio of 10% and phase number 3. While both networks gradually improve PSNR with more epochs, NL-EPGN appears to be significantly more effective than ISTA-Net${}^{+}$ during training as the former produces reconstructions with much higher PSNR.

4.2. MR Images Compressive Sensing

We also test the performance of EPGN on compressive sensing reconstruction of brain MR images [47] (CS-MRI). We randomly selected 100 and 50 images for training and testing, respectively, and cropped every image to the size of $190\times 190$. In the CS-MRI problem, the data fidelity is $f(\mathit{x};\mathit{y})={\parallel \mathrm{\Phi }\mathit{x}-\mathit{y}\parallel }_2^2$, where $\mathrm{\Phi }=P\mathcal{F}$, P is a binary selection matrix representating the sampling trajectory, and $\mathcal{F}$ is the discrete Fourier transform. We compare EPGN with ISTA-Net${}^{+}$ [27] on the same MRI data set. The experimental results on various undersampling ratios of radial masks are summarized in

Table 3. PSNR (dB) of reconstructions obtained by ISTA-Net${}^{+}$ and EPGN on MR images using radial masks with sampling ratios of 10%, 20%, and 30%.

Method	10%	20%	30%
ISTA-Net${}^{+}$	33.49	40.66	44.70
EPGN	33.70	40.94	45.45

. Here, we set the phase number of ISTA-Net${}^{+}$ and EPGN to 15 and 11 respectively. It is obvious that EPGN outperforms ISTA-Net${}^{+}$ for each undersampling ratio.

5. Concluding Remarks

We presented a novel deep neural network architecture, called the extra proximal gradient network (EPGN), to solve a general class of inverse problems with a focus on image reconstruction applications. EPGN imitates the accelerated extra proximal gradient algorithm and features a learned regularization that incorporates adaptive sparsification mappings, robust shrinkage selections, and the combination of local and nonlocal operators for improved solution quality and network parameter efficiency. Extensive numerical experiments show that EPGN outperforms several existing state-of-the-art methods on a variety of image reconstruction problems.