Dense Recurrent Neural Networks for Inverse Problems: History-Cognizant Unrolling of Optimization Algorithms

Inverse problems in medical imaging applications incorporate domain-specific knowledge about the forward encoding operator in a regularized reconstruction framework. Recently physics-driven deep learning (DL) methods have been proposed to use neural networks for data-driven regularization. These methods unroll iterative optimization algorithms to solve the inverse problem objective function, by alternating between domain-specific data consistency and data-driven regularization via neural networks. The whole unrolled network is then trained end-to-end to learn the parameters of the network. Due to simplicity of data consistency updates with gradient descent steps, proximal gradient descent (PGD) is a common approach to unroll physics-driven DL reconstruction methods. However, PGD methods have slow convergence rates, necessitating a higher number of unrolled iterations, leading to memory issues in training and slower reconstruction times in testing. Inspired by efficient variants of PGD methods that use a history of the previous iterates, we propose a history-cognizant unrolling of the optimization algorithm with dense connections across iterations for improved performance. In our approach, the gradient descent steps are calculated at a trainable combination of the outputs of all the previous regularization units. We also apply this idea to unrolling variable splitting methods with quadratic relaxation. Our results in reconstruction of the fastMRI knee dataset show that the proposed history-cognizant approach reduces residual aliasing artifacts compared to its conventional unrolled counterpart without requiring extra computational power or increasing reconstruction time.

Subsequently, iterative algorithms are employed to solve the regularized inverse problem by alternating between enforcing physics-driven data consistency and regularization using preselected priors.
Recently, deep learning (DL) has emerged as an alternative state-of-the art technique for solving such inverse problems in imaging applications . Several of the earlier works concentrated on training convolutional neural networks (CNN) to solve the problem without incorporating explicit knowledge of the acquisition physics, typically requiring re-training when acquisition parameters are modified [18][19][20][21][22][23][24][25][26][27][28][29][30]. Another line of work focuses on a physics-driven formulation, where the known forward model is incorporated into training, utilizing domain knowledge to solve the corresponding inverse problem [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45]51]. In physics-driven approaches, a conventional iterative optimization problem, which lands itself to a regularized least-squares formulation for modalities such as MRI and CT, is unrolled for a fixed number of iterations. This unrolled network alternates between data consistency and regularization steps similar to the conventional algorithm, where the regularizers are implemented via neural networks. The unrolled network is trained end-to-end to learn the network parameters, which predominantly characterize the regularization units.
In this study, we seek to develop a new methodology for unrolling inverse problem optimization algorithms into dense recurrent neural networks (RNN) for improved reconstruction performance. Motivated by accelerated PGD approaches, Fig. 1: (a) The conventional unrolling for iterative algorithms for inverse problems and (b) the proposed history-cognizant unrolling leading to a dense recurrent neural network architecture. In both cases, the algorithms are unrolled for T iterations, where each unrolled iteration consists of proximal (or regularization) unit (P), and a data consistency unit (DC). The P unit is implemented implicitly via a neural network, while the DC unit corresponds to a gradient descent step for the proximal gradient descent algorithms and typically to conjugate gradient for variable splitting methods. For the proposed history-cognizant unrolling, an additional 1 × 1 convolutional layer combines all previous proximal operator outputs prior to the DC unit, where C denotes concatenation.
which calculate the gradient descent at a combined history of previous iterations rather than the immediate previous iteration only [52][53][54][55][56][57][58], we propose a history-cognizant unrolling of the optimization algorithms with skip connections across unrolled iterations, leading to a dense RNN architecture. Skip connections have been previously shown to facilitate information flow through neural networks in architectures such as ResNet [59] and DenseNet [60]. Different than these previous neural network design problems, we propose to utilize such connections in an unrolled RNN architecture for solving inverse problems, which readily comes with its domainspecific design considerations. We compare the conventional unrolling and the proposed history-cognizant unrolling of the PGD algorithm in accelerated multi-coil MRI reconstruction using the fastMRI knee dataset [61], showing both visual and quantitative improvement in reconstruction quality with the proposed methodology. We also extend the same unrolling strategy to other optimization algorithms and investigate the performance improvement.

A. Forward Model and Inverse Problem
Let y denote a vector of noisy measurements from an imaging system and x be the corresponding underlying image. The forward model for this system is generally given as y = A(x) + n, where A is called a forward operator and n is measurement noise. In a number of medical imaging scenarios, such as MRI and CT reconstruction, the forward model is linear, leading to a simplified model where A is a linear forward operator and n is measurement noise, which will be the focus of the rest of the study. The corresponding inverse problem is typically ill-conditioned. Thus, regularization is frequently used, leading to a regularized least-squares method: arg min where the first term enforces consistency with measurement data and R(·) denotes a regularizer. In classical reconstruction approaches, regularizers have included Tikhonov [1,2], edgepreserving functions [62] or sparsity-promoting terms in predefined domains [3,4,6,7]. A variety of algorithms can iteratively solve the optimization problem in (2) such as PGD [63], VS [64][65][66] and PD [67,68] methods.

B. Physics-Based DL Reconstruction by Unrolling PGD
Several physics-based DL methods have unrolled a PGDbased method [32][33][34][35][36] for a fixed number of iterations. Such PGD-based unrolling alternates between a proximal operation that is implicitly defined by a learned neural network, and a gradient descent step to enforce data consistency, as follows: where µ i is the gradient descent step size at unrolled iteration i, and the proximal operation in (3a) is performed by a neural network. This leads to the schematic in Figure 1a, where the unrolled RNN consists of multiple blocks of proximal (P) and data consistency (DC) units. The former is a neural network and the latter implements a gradient descent step based on the known forward encoding operator. The unrolled RNN is then trained end-to-end to learn the mapping from the acquired sensor-domain data to reference images for medical imaging reconstruction.

C. Proposed History-Cognizant PGD Unrolling with Dense Recurrent Neural Networks
One of the major appeals of the PGD unrolling scheme is its computational simplicity. However, PGD algorithm itself has a slow convergence rate [62], potentially requiring a large number of iterations (i.e., a large number of blocks in Figure  1a) when unrolled into a neural network for sufficient reconstruction quality. In turn, this increases memory requirements and impedes training complexity. On the other hand, several methods have been proposed to accelerate convergence rate of PGD [52][53][54][55][56][57][58]. The underlying theme for these methods is to calculate the gradient step not only at the most recent proximal operation output, but at a linear combination of all past such outputs.
We propose to adapt this general accelerated PGD approach for network unrolling. In this setting, the accelerated PGD methods involve the following steps: where F(·) forms a linear combination of the previous proximal operator outputs. As a special case, Nesterov's method uses the previous two outputs [53], while more general versions have also been studied, in which all the past iterations are used [54,58]. For our unrolled network, we propose to learn F(·) during training instead of specifying a pre-determined form for the linear combination of previous iterates F(·) as in Nesterov's method [53]. This history-cognizant unrolling of the accelerated PGD algorithms leads to a dense recurrent neural network architecture (Dense-RNN) with skip connections across all unrolled iterations, as depicted in Figure 1b. In addition to the P and DC units to perform proximal operation (or regularization) and data consistency updates, each unrolled iteration also consists of a 1 × 1 convolutional layer implementing F(·) on the concatenation (denoted by C) of all past proximal operator outputs prior to the DC unit. We refer to this method as history-cognizant PGD (HC-PGD) in the rest of this study.

III. MULTI-COIL MRI RECONSTRUCTION EXPERIMENTS
AND RESULTS The proposed history-cognizant algorithm unrolling approach was compared to conventional algorithm unrolling in multi-coil MRI reconstruction. In a multi-coil MRI system with n c coils, y is the sub-sampled noisy data from all coils and x is the underlying image to be recovered. The forward linear operator A : C M ×N → C P in this system includes sensitives of the receiver coil array and a partial Fourier matrix for sub-sampling in k-space, i.e. spatial Fourier domain [1].

A. Knee MRI Datasets
Knee MRI datasets were obtained from the New York University (NYU) fastMRI initiative database [61] to investigate accelerated MRI reconstruction using conventional unrolled RNN and Dense-RNN architectures. Data were acquired without any acceleration on a clinical 3T system (Magnetom Skyra, Siemens, Germany) with a 15-channel knee coil. Data acquisition was performed using 2D turbo spin-echo sequences in coronal orientation with proton-density (Coronal PD) and proton-density weighted with fat suppression (Coronal PD-FS) weightings. Relevant imaging parameters were: resolution = 0.49 × 0.44mm 2 , slice thickness = 3mm , matrix size = 320 × 368 for both datasets.
All data were uniformly sub-sampled at an acceleration rate of 4 by keeping the central 24 lines using the mask provided in the fastMRI database [61] retrospectively. Uniform undersampling for 2D acquisitions is typically considered a more challenging setup for DL reconstruction than random undersampling [31]. A 24 × 24 central window was used to generate coil sensitivity maps using ESPIRiT [69]. The training sets consisted of 300 slices from 15 subjects for each sequence. Testing was performed on 10 different subjects for both sequences.

B. Implementation Details
The optimization algorithms were unrolled for T = 10 iterations. The proposed HC-PGD algorithm was implemented using the Dense-RNN architecture in Figure 1b. Proximal operator unit outputs are first concatenated and then combined using convolutional layers with a 1 × 1 kernel size and 2 output channels (real and imaginary components) before being inputted to the subsequent data consistency unit. On the other hand, the PGD algorithm was implemented using a conventional architecture that only includes P and DC units, without skip connections, concatenation of previous iterates or 1 × 1 convolutions.
The same CNN was used to implement the proximal operation for both PGD and HC-PGD. A residual network (ResNet) that has been established in a different regression problem was utilized [70]. This CNN consisted of 15 sub-residual blocks each having two convolutional layers of kernel size = 3×3 and output channels = 64. The first and second convolutional layers in each sub-residual block are followed by a ReLU activation function and a scaling factor of 0.1, respectively. In addition to sub-residual blocks, two input and output convolutional layers without any nonlinear activation function are to match number of desired channels. For both algorithm unrolling schemes, the proximal units shared parameters across unrolled iterations [39], leading to a total of 592,138 and 592,358 parameters for PGD and HC-PGD, respectively.
All networks were trained end-to-end, as summarized in Figure 2, by using Adam optimizer to minimize where x i ref , y i and A i are the reference SENSE-1 image, undersampled data and forward encoding operator for slice i, respectively and N is the number of training slices. f (y i , A i ; θ) denotes the network output for slice i with the corresponding measurements and forward encoding operator, while θ contains the trainable parameters of the networks. A normalized 1 − 2 loss [49] was used for L(·, ·): with u and v representing reference fully-sampled and network output images, respectively. Training was performed over 100 epochs. All algorithms were implemented using TensorFlow, and processed on a workstation with an Intel E5-2640V3 CPU (2.6GHz and 256GB memory), and two NVIDIA Tesla V100 GPUs with 32GB memory.

C. Performance of HC-PGD versus PGD Unrolling
Experiments were conducted on both Coronal PD and PD-FS datasets to compare the reconstruction quality of the conventional PGD unrolling and the proposed HC-PGD unrolling with the Dense-RNN architecture. A total of 380 slices from 10 subjects (excluded from training) were utilized for testing. Reconstructed testing data were evaluated quantitatively using normalized mean square error (NMSE) and structural similarity index (SSIM) metrics with respect to the fully-sampled reference image. Figure 3 depicts a representative slice from the Coronal PD dataset reconstructed using DL-MRI with PGD and HC-PGD unrolling methods. The first row shows the reference fully-sampled image, as well as the zero-filled acquired image, which serves as the input x (0) to the network. Difference images with respect to the reference image, scaled by a factor of 10, are also provided to facilitate comparison by showing errors. DL-MRI based on the proposed HC-PGD unrolling was able to suppress the strong aliasing artifacts successfully, while residual aliasing is visible on the DL-MRI reconstructions based on conventional PGD unrolling. For both methods, there is a slight degree of blurring in the final images, consistent with previous reports on DL-MRI methods with similar training loss functions [31]. Quantitative evaluations, shown in the lower left corner of the images, also confirm the improvement using the proposed history-cognizant unrolling with Dense-RNN architecture compared to the conventional algorithm unrolling.
The reconstructed images for a representative slice from the coronal PD-FS dataset in addition to the difference images scaled by a factor of 10, as well as the reference and zero-filled images are shown in Figure 4. There are slight residual aliasing artifacts in the DL-MRI reconstruction with PGD unrolling. The proposed HC-PGD unrolling successfully removes these artifacts further. The quantitative metrics also show improvement for the reconstruction results in moving from PGD to HC-PGD unrolling. Figure 5 summarizes the quantitative analysis of the DL-MRI reconstructions using NMSE and SSIM metrics for the PGD and HC-PGD unrollings for both the coronal PD and PD-FS datasets. The box plots depict the interquartile range with the center at the median of the metric for the corresponding reconstruction methods. Both metrics are visibly improved over the whole test sets when using the proposed history-cognizant algorithm unrolling over the conventional unrolling approach. Further statistical tests were performed using the Wilcoxon signed rank test with a significance level of P < 0.05, indicating that all differences are statistically significant based on both metrics and in both datasets.

D. Insensitivity of Dense-RNN Improvement to the Regularization CNN Choice
In order to verify that the improvements between PGD and HC-PGD unrolling are not limited to the choice of CNN used Fig. 3: A representative slice from the Coronal PD dataset with 4-fold acceleration rate using a uniform undersampling pattern, reconstructed with DL-MRI reconstruction based on the conventional PGD and the proposed HC-PGD unrolling, as well as error images with respect to the reference fully-sampled image scaled by a factor of 10. Zerofilled image corresponds to the input to the network, x (0) . The green rectangle in the reference image mark the zoom-up area shown in the second and last rows for each method. Residual aliasing artifacts, which are visible in the PGD unrolling results, have been suppressed using the HC-PGD scheme. Improvements in the quantitative metrics, shown on the lower-left corner of the images, also align with these visual observations. for the proximal operation, the two approaches in Figure 1 were also implemented using a U-Net architecture for the regularization CNN. This was based on the U-net used in Fig. 4: A representative slice from the Coronal PD-FS dataset with a 4-fold acceleration rate using a uniform undersampling pattern, reconstructed with conventional PGD and the proposed HC-PGD methods, along with the error images scaled by a factor of 10. The green rectangle in the reference image mark the zoom-up area. Proposed HC-PGD unrolling successfully removes some residual artifacts that are visible in the PGD. The quantitative metrics align with this observation. [61], which includes contracting and expanding paths up to 256 channels without batch normalization. Figure 6 depicts the results of this experiment on the coronal PD datasets, where the proposed history-cognizant unrolling outperforms the conventional unrolling methodology both visually and quantitatively. This highlights that the improvement from the

E. Performance Evaluation for Random Undersampling
In order to evaluate performance of the proposed historycognizant unrolling approach in random undersampling, we utilized an observation from a previous study [31,71], showing that physics-based networks trained on uniformly undersampled data can be transferred to random undersampling patterns at the same rate directly, even though the opposite is not necessarily true. Thus, we employed the networks trained on uniformly undersampled data, as detailed in Section III-B, to reconstruct test slices which were randomly undersampled at a similar 4-fold rate. The corresponding results, shown in Figure 7 indicate that the improvement trend remains similar to the uniform sampling scenario with the proposed HC-PGD method outperforming the conventional PGD in suppressing residual artifacts as observed in both Coronal PD and PD-FS datasets, which also align with the quantitative evaluation metrics in the images.

F. Extensions to Other Optimization Algorithms
While PGD unrolling is commonly used for DL-MRI reconstruction [32][33][34][35][36], alternative algorithms have also been utilized [37][38][39][40][41][42]. One such method for solving the objective function in Equation (2) is VS with quadratic penalty, which decouples data consistency and regularization terms into two blocks of variables [62] by introducing auxiliary variables to formulate a constrained objective function. This is then relaxed to an unconstrained problem by enforcing similarity between these variables via a quadratic penalty: arg min where z is an auxiliary variable and β is the quadratic penalty parameter. When this algorithm is unrolled as a neural network [38][39][40][41][42], reconstruction alternates between a proximal operation that is again learned by a neural network and data consistency, similar to PGD unrolling, as follows: where regularization update in (8a) is performed by a neural network and β was absorbed into the regularizer for notational convenience. The data consistency update in (8b) is computationally challenging in many practical scenarios, such as multicoil MRI reconstruction. Thus, this data consistency update can be solved iteratively via conjugate gradient (CG) to avoid matrix inversion [39].
As an extension to the HC-PGD approach, VS with quadratic relaxation can also be unrolled in a history-cognizant manner using the Dense-RNN in Figure 1, building on the optimization approaches in [56,57]. This will be referred to as the history-cognizant VS method (HC-VS). We note that the main difference to HC-PGD is that the DC units in HC-VS now implement an unrolled CG algorithm for a different least squares problem instead of a gradient step. To test the applicability of the history-cognizant unrolling scheme to different baseline algorithms, both VS and HC-VS were also implemented on the coronal PD and PD-FS datasets, using the setup described in Section III-B. In each data consistency unit, CG was implemented with 10 internal iterations.

IV. DISCUSSION
In this study, we proposed a history-cognizant method for unrolling optimization problems to solve regularized inverse problems in medical imaging using physics-based DL reconstruction. Our unrolling method builds on optimization approaches that utilize a history of the previous iterates via their linear combinations in updating the current estimate [52][53][54][55][56][57]. While these optimization procedures have been traditionally proposed to improve convergence rates or provide stronger theoretical guarantees, in the context of algorithm unrolling for physics-based DL reconstruction, they are used to enhance reconstruction quality without incurring additional computational complexity. Our method was pre-dominantly motivated by accelerated PGD approaches [52][53][54][55] and was shown to improve upon conventional PGD unrolling for multi-coil MRI reconstruction for multiple scenarios. We also explored how the unrolling can be applied to VS with quadratic relaxation for solving a regularized least squares problem, and showed that it can improve over the conventional unrolling approach.
The proposed history-cognizant unrolling methodology leads to several skip connections across iterations, leading to a densely-connected unrolled neural network. Dense neural networks have been shown to improve training by facilitating information flow through the network in forward propagation and by relieving issues concerned with gradient vanishing in backward propagation [60]. One of the advantages of these dense connections in the proposed unrolling is that the only computational overhead compared to conventional unrolling is the calculation of the 1 × 1 convolutions from the previous iterations. Thus, once trained, the computational time is increased <1.5% compared to the conventional counterpart, as depicted in Table I. One of the most commonly used accelerated PGD approaches is based on Nesterov acceleration [53]. In this case, the gradient descent is evaluated at a linear combination of the previous two iterates weighted by a step size µ and 1 − µ respectively. Theoretically, this approach achieves the optimal convergence rate for first-order methods [53]. As such, Nesterov-type unrolling for DL reconstruction was readily explored in [36] for coded-illumination microscopy, where the unrolled accelerated PGD network used the history of the two previous outputs of the proximal units. However, the main focus of this work was developing an illumination pattern design algorithm for a fixed reconstruction algorithm. On the other hand, when focusing on the reconstruction aspect for inverse problems, we found that the Nesterov-type unrolling was outperformed by the full history-cognizant Dense-RNN architecture in our initial experiments (not shown). Furthermore, the full history-cognizant unrolling subsumes Nesterov-type as a special case. Thus, we focused on this full history-cognizant case throughout our study to maximize the performance benefits.
In our experiments, the history-cognizant unrolling visibly improved the performance of PGD-based methods compared to its conventional counterpart for all datasets. In the VS setting, especially for the coronal PD datasets, the improvement is less pronounced. Coronal PD datasets have higher SNR compared to PD-FS datasets [31]. In this higher SNR case, the conventional unrolling of the VS algorithm readily provides a high-quality reconstruction. Thus, additional improvement with HC-VS is less apparent than in the PGD setting.
Our experiments in multi-coil MRI reconstruction pre-dominantly utilized uniform undersampling patterns. While random undersampling patterns are commonly used in literature for training and testing DL-based reconstruction, uniform undersampling remains the most clinically used approach, and was considered to be a harder problem compared to random undersampling in [31,71], as detailed in Section III-E.
To the best of our knowledge, the Dense-RNN architecture that arises from the proposed history-cognizant unrolling have not been used for regression problems in general, and inverse problems in particular. However, Dense-RNNs have been independently proposed for different problems, such as scene parsing or labeling [72] where semantic information of an image segment is to be extracted. In such applications, RNNs are used for each image unit to receive the dependencies from other units through recurrent information forwarding between adjacent units, which decays for further units. In this setting, Dense-RNNs improve labeling by capturing longrange semantic dependencies among image units, highlighting a similar memory dependence as in this study.

V. CONCLUSION
We proposed a history-cognizant approach to unroll a given optimization algorithm for solving inverse problems in medical imaging using physics-driven deep learning. The proposed approach leads to an unrolled dense recurrent neural network architecture. From a network-design perspective, this consequent architecture enhances training and performance by facilitating information flow through the network. Accelerated MRI reconstruction experiments show that the proposed unrolling approach may considerably enhance reconstruction quality without substantially modifying computational complexity.
ACKNOWLEDGMENT MRI data were obtained from the NYU fastMRI initiative database [61] which was acquired with the relevant institutional review board approvals as detailed in [61]. NYU fastMRI investigators provided data but did not participate in analysis or writing of this report. A listing of NYU fastMRI investigators, subject to updates, can be found at fastmri.med.nyu.edu.