* [10] H. V. Nguyen, V.-D. Nguyen, O. A. Dobre, S. K. Sharma, S. Chatzinotas, B. Ottersten, and O.-S. Shin, “On the spectral and energy efficiencies of full-duplex cell-free massive MIMO,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1698–1718, 2020.
  * [11] P. Anokye, R. K. Ahiadormey, and K.-J. Lee, “Full-duplex cell-free massive MIMO with low-resolution ADCs,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 11, pp. 12 179–12 184, 2021.
  * [12] A. Chowdhury and C. R. Murthy, “Comparative study of dynamic TDD with full-duplex in cell-free massive MIMO systems,” in _Proc. NCC, Guwahati, India_ , 2023.
  * [13] Q. Ding, Y. Lian, and Y. Jing, “Performance analysis of full-duplex massive MIMO systems with low-resolution ADCs/DACs over rician fading channels,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 7, pp. 7389–7403, 2020.
  * [14] X. Xia, P. Zhu, J. Li, H. Wu, D. Wang, Y. Xin, and X. You, “Joint user selection and transceiver design for cell-free with network-assisted full duplexing,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 12, pp. 7856–7870, 2021.
  * [15] K. Shen and W. Yu, “Fractional programming for communication systems—part I: Power control and beamforming,” _IEEE Trans. on Signal Process._ , vol. 66, no. 10, pp. 2616–2630, 2018.
  * [16] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, _Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers_. now, 2011.
  * [17] T. Riihonen, S. Werner, and R. Wichman, “Mitigation of loopback self-interference in full-duplex MIMO relays,” _IEEE Trans. Signal Process._ , vol. 59, no. 12, pp. 5983–5993, 2011.
  * [18] M. Garey and D. Johnson, _Computers and Intractability: A Guide to the Theory of NP-completeness_ , ser. Mathematical Sciences Series. Freeman, 1979.
  * [19] D. Brélaz, “New methods to color the vertices of a graph,” _Commun. ACM_ , vol. 22, no. 4, p. 251–256, apr 1979.
  * [20] R. Sabbagh, C. Pan, and J. Wang, “Pilot allocation and sum-rate analysis in cell-free massive MIMO systems,” in _Proc. ICC_ , 2018, pp. 1–6.
  * [21] F. Riera-Palou, G. Femenias, A. G. Armada, and A. Pérez-Neira, “Clustered cell-free massive MIMO,” in _IEEE Globecom Workshops_ , 2018, pp. 1–6.
  * [22] H. Liu, J. Zhang, X. Zhang, A. Kurniawan, T. Juhana, and B. Ai, “Tabu-search-based pilot assignment for cell-free massive MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 2, pp. 2286–2290, 2020.
  * [23] Y. Zhang, H. Cao, P. Zhong, C. Qi, and L. Yang, “Location-based greedy pilot assignment for cell-free massive MIMO systems,” in _Proc. IEEE ICCC_ , 2018, pp. 392–396.
  * [24] T. C. Mai, H. Q. Ngo, M. Egan, and T. Q. Duong, “Pilot power control for cell-free massive MIMO,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 11, pp. 11 264–11 268, 2018.
  * [25] A. Lozano, A. Tulino, and S. Verdu, “Multiple-antenna capacity in the low-power regime,” _IEEE Trans. Inf. Theory_ , vol. 49, no. 10, pp. 2527–2544, 2003.
  * [26] Y. Gao, H. Vinck, and T. Kaiser, “Massive MIMO antenna selection: Switching architectures, capacity bounds, and optimal antenna selection algorithms,” _IEEE Trans. Signal Process._ , vol. 66, no. 5, pp. 1346–1360, 2018.
  * [27] R. Nikbakht, R. Mosayebi, and A. Lozano, “Uplink fractional power control and downlink power allocation for cell-free networks,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 6, pp. 774–777, 2020.
  * [28] R. Nikbakht and A. Lozano, “Uplink fractional power control for cell-free wireless networks,” in _Proc. IEEE ICC_ , 2019, pp. 1–5.
  * [29] G. Interdonato, P. Frenger, and E. G. Larsson, “Scalability aspects of cell-free massive MIMO,” in _Proc. IEEE ICC_ , 2019, pp. 1–6.


[1]

[cor1]Corresponding author.

# An Optimal Control Framework for Joint-channel Parallel MRI Reconstruction
without Coil Sensitivities

Wanyu Bian<EMAIL_ADDRESS>Yunmei Chen Department of Mathematics,
University of Florida, Gainesville, FL 32601, USA  Xiaojing Ye Department of
Mathematics and Statistics, Georgia State University, Atlanta, GA 30303, USA

###### Abstract

_Goal_ : This work aims at developing a novel calibration-free fast parallel
MRI (pMRI) reconstruction method incorporate with discrete-time optimal
control framework. The reconstruction model is designed to learn a
regularization that combines channels and extracts features by leveraging the
information sharing among channels of multi-coil images. We propose to recover
both magnitude and phase information by taking advantage of structured
convolutional networks in image and Fourier spaces. _Methods_ : We develop a
novel variational model with a learnable objective function that integrates an
adaptive multi-coil image combination operator and effective image
regularization in the image and Fourier spaces. We cast the reconstruction
network as a structured discrete-time optimal control system, resulting in an
optimal control formulation of parameter training where the parameters of the
objective function play the role of control variables. We demonstrate that the
Lagrangian method for solving the control problem is equivalent to back-
propagation, ensuring the local convergence of the training algorithm.
_Results_ : We conduct a large number of numerical experiments of the proposed
method with comparisons to several state-of-the-art pMRI reconstruction
networks on real pMRI datasets. The numerical results demonstrate the
promising performance of the proposed method evidently. _Conclusion_ : The
proposed method provides a general deep network design and training framework
for efficient joint-channel pMRI reconstruction. _Significance_ : By learning
multi-coil image combination operator and performing regularizations in both
image domain and k-space domain, the proposed method achieves a highly
efficient image reconstruction network for pMRI.

###### keywords:

Parallel MRIReconstruction Discrete-time optimal controlResidual learning.

##  1 Introduction

Magnetic resonance imaging (MRI) is one of the most prominent medical imaging
technologies with extensive clinical applications. In clinical applications,
an advanced medical MRI technique known as parallel MRI (pMRI) is widely used.
PMRI surrounds the scanned objects with multiple receiver coils and collects
k-space (Fourier) data in parallel. PMRI can reduce the data acquisition time
and has become the state-of-the-art technology in modern MRI applications. To
accelerate the scan process, partial data acquisitions that increase the
spacing between read-out lines in k-space are implemented in pMRI. However,
this results in aliasing artifacts, and a proper image reconstruction process
is necessary to recover the high-quality artifact-free images from the partial
data.

Two major approaches are commonly addressed to image reconstruction in pMRI:
the first approach is k-space method which interpolates the missing k-space
data using the sampled ones across multiple receiver coils [15], such as the
generalized auto-calibrating partially parallel acquisition (GRAPPA) [18] and
simultaneous acquisition of spatial harmonics (SMASH) [49]. The other approach
is the class of image space method that eliminate the aliasing artifacts in
the image domain by solving a system of equations that relate the image to be
reconstructed and partial k-space data through coil sensitivities, such as
SENSitivity Encoding (SENSE) [45]. Coil sensitivity maps are indispensable and
required to be accurately pre-estimated in traditional SENSE-based methods.
Traditional pMRI reconstruction methods in image space follows SENSE-based
framework, which is formulated as an optimization problem that minimizes a
summation of a data fidelity term and a weighted regularization term. The
detailed explanation about this formulation can be refer to [30].

In recent years, we have witnessed fast developments of SENSE-based pMRI
reconstruction incorporate with deep-learning based methods [30, 36, 53, 20,
38]. There are two critical issues that need to be carefully addressed. The
first issue is on the choice of regularization including the regularization
weight. The regularization term is of paramount importance to the severely
ill-posed inverse problem of pMRI reconstruction due to the significant
undersampling of k-space data. In the past decades, most traditional image
reconstruction methods employ handcrafted regularization terms, such as the
total variation (TV). In recent years, a class of unrolling methods that mimic
classical optimization schemes are developed, where the regularization is
realized by deep networks whose parameters are learned from data. However, the
learned regularization is often cast as a black-box that is difficult to
interpret, and the training can be very data demanding and prone to
overfitting especially when the networks are over-parameterized [31, 4, 23].

The second issue is due to the unavailability of accurate coil sensitivities
$\\{\mathbf{s}_{i}\\}$ in practice. Inaccurate coil sensitivity maps lead to
severe biases that degrade the quality of reconstructed $\mathbf{v}$. One way
to eliminate this issue is to (regularize and) reconstruct multi-coil images,
and combine channels into a full-body image in the final step by taking some
hand-crafted methods such as the root of sum of squares (RSS). Different from
RSS, our method proposed a learnable multi-coil combination operator to
combine channel-wise multi-coil images.

In this paper, we tackle the aforementioned issues in a discrete-time optimal
control framework to optimize the variational pMRI reconstruction model. We
highlight several main features of our framework as follows.

  1. 1.

Unlike most existing methods which regularize and reconstruct multi-coil
images, we employ regularization in both image and Fourier spaces to improve
reconstruction quality.

  2. 2.

Our method advocates a learned adaptive combination operator that first merges
multi-coil images into a full-body image with a complete field of view (FOV),
followed by an effective regularization on this image. This is in sharp
contrast to existing methods which only combine reconstructed multi-coil
images in the final step, whereas our regularizer leverages the combination
operator in each iteration which improves parameter efficiency.

  3. 3.

We employ a complex-valued neural network as the coil combination operator to
recover both magnitude and phase information of pMRI images when coil
sensitivity is unavailable. This combination method benefits from the coil
information shared among multiple channels, which is distinct from most hand-
crafted coil combination methods.

  4. 4.

We propose a novel deep reconstruction network whose structure is determined
by the discrete-time optimal control system for minimizing the objective
function, which yields an optimal control formulation where the parameters of
the combination and regularization operators play the role of control
variables of the discrete dynamic system. The optimal value of these
parameters is obtained by a Lagrangian method which can be implemented using
back-propagation.

We consider two clinical pMRI sequences of knee images and verified the
effective performance of the proposed combination operator, different initial
reconstructions, complex convolutions, and domain-hybrid network in the
Ablation Studies. The proposed network recovers both magnitude and phase
information of pMRI images. The effect of the aforementioned techniques
demonstrate evident improvement of reconstruction quality and parameter
efficiency using our method. For reproducing the experiment, our code is
available at https://github.com/1lol/pMRI_optimal_control.

This remainder of the paper is organized as follows: In section 2, we provide
an overview of recent developments in pMRI, cross-domain reconstructions,
complex-valued CNNs, and optimal control inspired deep training models that
related to our work. We present our proposed problem settings and
reconstruction network architecture in detail in Section 3. Extensive
numerical experiments and analyses on a variety of clinical pMRI data are
presented in Section 4. Section 5 concludes this paper.

##  2 Related work

In recent years, we have witnessed fast developments of medical imaging
incorporate with deep-learning based methods [24, 26, 25]. Most existing deep-
learning based methods rendering end-to-end neural networks mapping from the
partial k-space data to the reconstructed images [58, 66, 46, 39, 67]. These
approaches require an excessive amount of training data, and the designed
networks are cast as black boxes whose underlying mechanism can be very
difficult to interpret. To mitigate this issue, a number of unrolling methods
were proposed to map existing optimization algorithms to structured networks
where each phase of the networks correspond to one iteration of an
optimization algorithm [21, 2, 3, 63, 10, 62, 65, 48]. In what follows, we
focus on the recent developments in deep-learning based image reconstruction
methods for pMRI.

Variational Network (VN) [21] was introduced to unroll the gradient descent
algorithm as a reconstruction network which requires precalculated
sensitivities $\\{\mathbf{s}_{i}\\}$ as input. MoDL [3] proposed a weight
sharing strategy in a recursive network to learn the regularization parameters
by unrolling the conjugate gradient method. Several methods explored different
strategies to avoid using pre-calculated coil sensitivity maps for pMRI
reconstruction. Blind-PMRI-Net [40] proposed pMRI model by regularizing
sensitivity maps and MR image, where their network alternatively estimates
coil images, sensitivities and single-body image by three subnets. De-
Aliasing-Net [9] proposed a de-aliasing reconstruction model with that applied
split Bregman iteration algorithm without explicit coil sensitivity
calculation. The de-aliasing network explored cross-correlation among channels
and spatial redundancy which provoked a desirable performance. LINDBERG [57]
explored calibration-free pMRI technique which uses adaptive sparse coding to
obtain joint-sparse representation precisely by equipping a joint sparsity
regularization to extract desirable cross-channel relationship. This work
proposed to alternatively update sparse representation, sensitivity encoded
images, and K-space data. Adaptive-CS-Net [44] is a leading method in 2019
fastMRI challenge [64] that unrolled modified ISTA-Net+ [65]. The proposed
calibration-free pMRI method distincts from above related works in terms of
the learnable multi-coil combination operator to adaptively combine channels
of the updated multi-coil images through iterations.

Recently, cross-domain methods exhibits its significance in medical imaging
[16, 60, 52, 41, 51, 28, 68, 50]. KIKI-net [16] iteratively applied k-space
CNN, image domain CNN and interleaved data consistency operation for single-
coil image reconstruction. CDF-Net [41] further shows adding communication
between spatial and frequency domain gives a boost in performance. Their
results indicated that domain-specific network has individual strong points
and disadvantages in restoring tissue-structure. Our reconstruction model is
inspired of cross-domain reconstruction, the difference is that we solve for
the reconstruction model with cross-domain regularization functions through a
learnable optimization algorithm instead of an end-to-end network.

Our network applies complex-valued convolutions and activation functions. MRI
data are complex-valued, and the phase signals also carry important
pathological information such as in quantitative susceptibility mapping [13,
47]. Cole $et$ $al.$ [13] investigated the performance of complex-valued
convolution and activation functions has better reconstruction over the model
with real-valued convolution in various network architectures. DeepcomplexMRI
[58] was developed to recover multi-coil images by unrolling an end-to-end
complex-valued network.

In supervised learning, deep residual neural networks can be approximated as
discretizations of a classical optimal control problem of a dynamical system,
where training parameters can be viewed as control variables [61, 32, 35].
Control inspired learning algorithms introduced a new family of network
training models which connect with dynamical systems. Pontryagin’s maximum
principle (PMP) [19] was explored as necessary optimality conditions for the
optimal control [33], [34], these works devise the discrete method of
successive approximations (MSA) [11] and its variance for solving PMP. Neural
ODE [8] models the continuous dynamics of hidden states by some certain types
of neural networks such as ResNet, the forward propagation is equivalent to
one step of discretatized ordinary differential equations (ODE). Inspired by
Neural ODE [8], Chen $et$ $al.$ [7] modeled ODE-based deep network for MRI
reconstruction. In this paper, we model the optimization trajectory as a
discrete dynamic process from the view of the method of Lagrangian
Multipliers.

The present work is a substantial extension of the preliminary work in [5]
using domain-hybrid network with a trained initialization to solve for an
optimal control pMRI joint-channel reconstruction problem when coil-
sensitivity is unavailable. More comprehensive empirical study is conducted in
this work.

##  3 Proposed Method

###  3.1 Background

PMRI as well as general MRI reconstruction can be formulated as an inverse
problem. Consider a pMRI system with $c$ receiver coils acquiring 2D MR images
at resolution $m\times n$ (we treat a 2D image
$\mathbf{v}\in\mathbb{C}^{m\times n}$ and its column vector form
$\mathbf{v}\in\mathbb{C}^{mn}$ interchangeably hereafter). Let
$\mathbf{P}\in\mathbb{R}^{p\times mn}(p\leq mn)$ be the binary matrix
representing the undersampling mask with $p$ sampled locations in k-space,
$\mathbf{s}_{i}\in\mathbb{C}^{mn}$ the coil sensitivity, and
$\mathbf{f}_{i}\in\mathbb{C}^{p}$ the _partial_ k-space data at the $i$-th
receiver coil for $i=1,\dots,c$. The partial data $\mathbf{f}_{i}$ and the
image $\mathbf{v}$ are related by
$\mathbf{f}_{i}=\mathbf{P}\mathbf{F}(\mathbf{s}_{i}\odot\mathbf{v})+\mathbf{n}_{i}$,
where $\mathbf{s}_{i}$ is the (unknown) sensitivity map at the $i$-th coil and
$\odot$ denotes entrywise product of two matrices,
$\mathbf{F}\in\mathbb{C}^{mn\times mn}$ stands for the (normalized) discrete
Fourier transform that maps an image to its Fourier coefficients, and
$\mathbf{n}_{i}$ is the unknown acquisition noise in k-space at the $i$-th
receiver coil. Then the variational model for image reconstruction can be cast
as an optimization problem as follows:

$\min_{\mathbf{v}}\
\sum\nolimits^{c}_{i=1}\frac{1}{2}\|\mathbf{P}\mathbf{F}(\mathbf{s}_{i}\odot\mathbf{v})-\mathbf{f}_{i}\|^{2}+R(\mathbf{v}),$
(1)

where $\mathbf{v}\in\mathbb{C}^{mn}$ is the single full-body MR image to be
reconstructed, $R(\mathbf{v})$ is a regularization on the image $\mathbf{v}$,
and $\|\mathbf{w}\|^{2}:=\|\mathbf{w}\|_{2}^{2}=\sum_{j=1}^{n}|w_{j}|^{2}$ for
any complex vector $\mathbf{w}=(w_{1},\dots,w_{n})^{\top}\in\mathbb{C}^{n}$.
Our approach is based on uniform Cartesian k-space sampling. Table 1 displays
the notations and their descriptions that we used in the paper.

Table 1: Some notations and meanings that used throughout this paper.

Expression | Description
---|---
$c$ | total number of receiver coils
$\mathbf{u}=(\mathbf{u}_{1},\dots,\mathbf{u}_{c})$ | multi-coil MRI data
$\mathbf{f}=(\mathbf{f}_{1},\dots,\mathbf{f}_{c})$ | partial k-space measurement
$\mathbf{s}=(\mathbf{s}_{1},\dots,\mathbf{s}_{c})$ | coil sensitivity map
$\mathbf{v}$ | full fov image that need to be reconstructed
$\mathbf{u}^{*}$ | ground truth multi-coil MRI data
$\mathbf{v}^{*}$ | ground truth single body MRI data
$\mathbf{F}$ | discrete Fourier transform
$\mathbf{F}^{H}$ | inverse discrete Fourier transform
$\mathbf{P}$ | undersampling trajectory
$\mathbf{n}$ | measurement noise
$\mathbf{g}$ | algorithm unrolling network
$\mathbf{g}_{0}$ | initial network
$\mathcal{J},\mathcal{G},\tilde{\mathcal{G}},\tilde{\mathcal{J}}$ | image space convolutional operators in $\mathbf{g}$
$\mathcal{K}$ | k-space convolutional operators in $\mathbf{g}$
$\mathcal{K}_{0}$ | k-space convolutional operator in $\mathbf{g}_{0}$
$h$ | data fidelity term
$R$ | regularization
RSS | square root of the sum of squares
$t=1,\cdots,T$ | phase number
$k=1,\cdots,K$ | number of iterations for Alg (1)
$\mathbf{U}=(\mathbf{u}{(0)},\cdots,\mathbf{u}{(T)})^{\intercal}$ | collection of predicted multi-coil images
| at each phase
$\Theta=(\theta(0),\cdots,\theta(T))^{\intercal}$ | collection of control variable (parameters)
| at each phase
$\Lambda=(\lambda(0),\cdots,\lambda(T))^{\intercal}$ | collection of Lagrangian multipliers of (4)

###  3.2 Problem Settings

We propose a unified deep neural network for calibration-free pMRI
reconstruction by recovering images from individual receiver coils jointly
that does not require any knowledge of coil-wise sensitivity profile. Denote
that $\mathbf{u}_{i}$ is the MR image at the $i$-th receiver coil and hence is
related to the full body image $\mathbf{v}$ by
$\mathbf{u}_{i}=\mathbf{s}_{i}\odot\mathbf{v}$. On the other hand, the image
$\mathbf{u}_{i}$ corresponds to the partial k-space data $\mathbf{f}_{i}$ by
$\mathbf{f}_{i}=\mathbf{P}\mathbf{F}\mathbf{u}_{i}+\mathbf{n}_{i}$, and hence
the data fidelity term is formulated as least squares
$\frac{1}{2}\sum^{c}_{i=1}\|\mathbf{P}\mathbf{F}\mathbf{u}_{i}-\mathbf{f}_{i}\|^{2}$.
We also need a suitable regularization $R$ on the images
$\\{\mathbf{u}_{i}\\}$. However, these images have severely geometrically
inhomogeneous contrasts due to the physical variations in the sensitivities
across the image domain at different receiver coils. Therefore, it is more
appropriate and effective to apply regularization to the single full-body
image $\mathbf{v}$ than to individual coil images in
$(\mathbf{u}_{1},\dots,\mathbf{u}_{c})$.

To address the issue of proper regularization, we propose to learn a nonlinear
operator $\mathcal{J}$ that combines $\\{\mathbf{u}_{i}\\}$ into the image
$\mathbf{v}=\mathcal{J}(\mathbf{u})\in\mathbb{C}^{m\times n}$, where
$\mathbf{u}=(\mathbf{u}_{1},\dots,\mathbf{u}_{c})\in\mathbb{C}^{m\times
n\times c}$ represents the channel-wise multi-coil MRI data that consists of
$\mathbf{u}_{i}$ for $i=1,\cdots,c$. Then we apply a suitable $R$ to the image
$\mathbf{v}$. We also introduce a k-space $R_{f}$ on
$\mathbf{F}\mathbf{u}_{i}$ to take advantage of Fourier information and
enhance the model performance.

We denote
$\mathbf{f}=(\mathbf{f}_{1},\dots,\mathbf{f}_{c})\in\mathbb{C}^{p\times c}$ as
the partial k-space measurements at $c$ sensor coils. Suppose that we are
given $N$ data pairs $\\{(\mathbf{f}^{(j)},\mathbf{u}^{*(j)})\\}_{j=1}^{N}$
for training the network, where $\mathbf{u}^{*(j)}$ is the ground truth multi-
coil MR data with index $j\in\\{1,\cdots,N\\}$. Let $\Theta$ represents the
parameters that need to be learned from network by minimizing the loss
function $\ell$. We formulate the network training as a bilevel optimization
problem, where the lower level is to update $\mathbf{u}$ with fixed trainable
parameters $\Theta$ and upper level is to update $\Theta$ that learned from
the training data by minimizing loss function $\ell$.

$\displaystyle\min_{\Theta}\ \
\frac{1}{N}\sum^{N}_{j=1}\ell(\mathbf{u}^{(j)}_{\Theta}),\ \ \ $ (2a)
$\displaystyle\mathrm{s.t.}\ \
\mathbf{u}_{\Theta}^{(j)}=\operatorname*{arg\,min}_{\mathbf{u}^{(j)}}\phi_{\Theta}(\mathbf{u}^{(j)}),\
\ $ (2b)

with $\phi_{\Theta}$ defined as below:

$\phi(\mathbf{u}):=\frac{1}{2}\sum\limits^{c}_{i=1}\|\textbf{PF}\mathbf{u}_{i}-\mathbf{f}_{i}\|^{2}+R(\mathcal{J}(\mathbf{u}))+R_{f}(\mathbf{F}\mathbf{u}_{i}).$
(3)

The objective function $\phi$ is the variational model for pMRI
reconstruction, in which $\Theta$ is the set collects all the parameters that
learned from the regularizers $R\circ\mathcal{J}$ and $R_{f}\circ\mathbf{F}$,
so $\phi$ is depending on $\Theta$. Problem 2 is formulated in the scenario of
reconstructing the multi-coil MRI data. The final reconstruction result
$\mathbf{u}_{\Theta}$ is the forward network output that dependent on network
parameters $\Theta$.

The deep learning approach for pMRI reconstruction in lower level problem (2b)
can be cast and formulated as a discrete-time optimal control system. We use
one data sample $(\mathbf{f},\mathbf{u}^{*})$ and omit the average and data
indexes for notation simplicity. The architecture of deep unrolling method
follows the iterations of optimization algorithms and solve for the minimizer
of the following problem:

$\displaystyle\min_{\Theta}\ \ \ $ $\displaystyle\ell(\mathbf{u}_{\Theta}),$
(4a) $\displaystyle\mathrm{s.t.}\ \ \ $
$\displaystyle\mathbf{u}{(t)}=\mathbf{g}(\mathbf{u}{(t-1)},\theta(t)),\
t=1,\cdots,T,$ (4b)
$\displaystyle\mathbf{u}{(0)}=\mathbf{g}_{0}(\mathbf{f},\theta(0)),$ (4c)

where $\Theta=(\theta(0),\cdots,\theta(T))^{\intercal}$ is the collection of
control variables $\theta(t)$ at all time steps (phases) respectively. In
(4a), $\mathbf{u}_{\Theta}=\mathbf{u}(T)$, which is the output image from the
last $T$-th phase of the entire network so that we want to find optimal
$\mathbf{u}(T)$ that close to
$\operatorname*{arg\,min}_{\mathbf{u}}\phi_{\Theta}(\mathbf{u})$. Equations
(4b), (4c) are inspired by deep unrolling algorithm for solving the lower
level constraint (2b). Given $\theta(t)$, $\mathbf{u}(t)$ is the state of
updated reconstruction multi-coil data from $t$-th phase for each
$t=0,\cdots,T$. $\mathbf{g}$ is a multi-phase unrolling network inspired by
the proximal gradient algorithm, and the output of
$\mathbf{g}(\cdot)\in\mathbb{C}^{m\times n\times c}$ is the updated multi-coil
MRI data from each phase. The network $\mathbf{g}$ is the intermediate mapping
from $\mathbf{u}{(t)}$ to $\mathbf{u}{(t+1)}$ for $t=0,\cdots,T-1$, whose
structure is explained in Section 3.4 for minimizing the variational model
(3). The network $\mathbf{g}_{0}$ with initial control parameter $\theta(0)$
maps the partial k-space measurement $\mathbf{f}$ to an initial reconstruction
$\mathbf{u}{(0)}$ as the input of this optimal control system.

To summarize in brief, we solve for the minimizer (reconstruction result) of
the lower level problem (2b) in a discrete-time optimal control framework (4).
The dynamic system (4b), (4c) as the constrain of (4a) is modeled as the
optimal control system of the variational model (3).

The loss function $\ell(\mathbf{u}{(T)})$ measures the discrepancy between the
final state $\mathbf{u}{(T)}$ and the reference image $\mathbf{u}^{*}$
obtained using full k-space data in the training data set. We set the loss
function in (2a) and (4a) for the proposed method as:

$\displaystyle\ell(\mathbf{u}_{\Theta})=\ell(\mathbf{u}{(T)})=\sum\nolimits^{c}_{i=1}\gamma\|\mathbf{u}_{i}{(T)}-\mathbf{u}^{*}_{i}\|+$
(5)
$\displaystyle\||\mathcal{J}(\bar{\mathbf{u}}{(T)})|-\text{RSS}(\mathbf{u}^{*})\|+\eta\|\text{RSS}(\bar{\mathbf{u}}{(T)})-\text{RSS}(\mathbf{u}^{*})\|,$

where
$\text{RSS}(\mathbf{u}^{*})=(\sum_{i=1}^{c}|\mathbf{u}^{*}_{i}|^{2})^{1/2}\in\mathbb{R}^{m\times
n}$ is the pointwise root of sum of squares across the $c$ channels of
$\mathbf{u}^{*}$, $|\cdot|$ is the pointwise modulus, and $\gamma,\eta>0$ are
prescribed weight parameters.

The motivation of applying learnable image space regularization
$R\circ\mathcal{J}$ and k-space regularization $R_{f}\circ\mathbf{F}$ in model
(3) is explained in the following: (i) Image domain network recovers the high
spatial resolution, but may not suppress some artifacts. Frequency domain
network is more suitable to remove high-frequency artifacts. (ii) Image domain
and k-space information are equivalent due to the global linear
transformation, but adding nonlinear activations with CNNs can feasibly
improve the efficacy of network learning and boost the reconstruction
performance. We parametrize the combination operator $\mathcal{J}$ by CNNs
since the partial k-space data were scanned by multiple coil arrays, and
introduce cross-correlation among channels of coil-images which could be
compatible for CNN structure.

###  3.3 Design of $\mathbf{g}$ and $\mathbf{g}_{0}$

Denote function $h(\mathbf{u})$ as the data fidelity term, one of the famous
traditional method for solving problem
$\min_{\mathbf{u}}h(\mathbf{u})+R(\mathbf{u})$ is the _proximal gradient
algorithm_ [42]:

$\displaystyle\mathbf{b}^{l}$ $\displaystyle=\mathbf{u}^{l}-\alpha^{l}\nabla
h(\mathbf{u}^{l}),$ (6a) $\displaystyle\mathbf{u}^{l+1}$
$\displaystyle=\operatorname{\mathrm{prox}}_{\alpha^{l}R}(\mathbf{b}^{l}),$
(6b)

where $\alpha^{l}>0$ is the step size, the proximity operator
$\operatorname{\mathrm{prox}}_{\alpha R}$ defined below:

$\operatorname{\mathrm{prox}}_{\alpha
R}(\mathbf{b}):=\operatorname*{arg\,min}_{\mathbf{x}}\left\\{\frac{1}{2\alpha}\|\mathbf{x}-\mathbf{b}\|^{2}+R(\mathbf{x})\right\\},$
(7)

where $[\mathbf{x}]_{i}=\mathbf{x}_{i}\in\mathbb{C}^{m\times n}$ for any
$\mathbf{x}\in\mathbb{C}^{m\times n\times c}$. The proposed network structure
is inspired by (6) for solving (2b) which can be cast as an iterative
procedure with $T$ phases in the discrete dynamic system (4b) and (4c). We
parametrize the $t$-th phase consists of three steps:

$\displaystyle\mathbf{b}_{i}{(t)}=\mathbf{u}_{i}{(t)}-\rho_{t}\mathbf{F}^{H}\mathbf{P}^{\top}(\mathbf{P}\mathbf{F}\mathbf{u}_{i}{(t)}-\mathbf{f}_{i}),$
$\displaystyle\quad i=1,\cdots,c,$ (8a)
$\displaystyle\bar{\mathbf{u}}_{i}{(t)}=[\operatorname{\mathrm{prox}}_{\rho_{t}R(\mathcal{J}(\cdot))}(\mathbf{b}{(t)})]_{i},$
$\displaystyle\quad i=1,\cdots,c,$ (8b)
$\displaystyle{\mathbf{u}}_{i}{(t+1)}=[\operatorname{\mathrm{prox}}_{\rho_{t}R_{f}(\mathbf{F}(\cdot))}(\bar{\mathbf{u}}{(t)})]_{i},$
$\displaystyle\quad i=1,\cdots,c,$ (8c)

for $t=0,\cdots,T-1$ and
$\mathbf{b}{(t)}=(\mathbf{b}_{1}{(t)},\dots,\mathbf{b}_{c}{(t)})\in\mathbb{C}^{m\times
n\times c}$, $\rho_{t}>0$ is the step size. $F$ denotes the normalized
discrete Fourier transform and $F^{H}$ the complex conjugate transpose (i.e.,
Hermitian transpose) of $F$, here $F^{H}$ is the inverse discrete Fourier
transform.

The step (8a) computes $\mathbf{b}{(t)}$ by applying the gradient decent
algorithm to minimize the data fidelity term in (3) which is straightforward
to compute. The first proximal update step (8b) equipped with the joint
regularizer $R(\mathcal{J}(\cdot))$ and upgrades its input
$\mathbf{b}_{i}{(t)}$ to a multi-coil image $\bar{\mathbf{u}}_{i}{(t)}$.
Ideally, the regularization $R\circ\mathcal{J}$ can be parameterized as a deep
neural network whose parameters can be adaptively learned from data, however,
in such case $\operatorname{\mathrm{prox}}_{\rho_{t}R(\mathcal{J}(\cdot))}$
does not have closed form and can be difficult to compute. As an alternative,
the proximity operator
$\operatorname{\mathrm{prox}}_{\rho_{t}R(\mathcal{J}(\cdot))}$ can be directly
parametrized as a learnable denoiser and solve (8b) in each iteration. For the
similar reason, the proximity operator
$\operatorname{\mathrm{prox}}_{R_{f}(\mathbf{F}(\cdot))}$ in (8c) can also be
parametrized as CNNs, which further improves the accuracy of the k-space
measurement. In both (8b) and (8c), the regularizers $R\circ\mathcal{J}$ and
$R_{f}\circ\mathbf{F}$ extract complex features through neural networks and
the proximity points can be learned in denoising network via ResNet [22]
structure.

We frame step (8b) incorporates joint reconstruction to update coil-images via
ResNet [22]:
$\bar{\mathbf{u}}{(t)}=\mathbf{b}{(t)}+\mathcal{M}(\mathbf{b}{(t)})$, where
$\mathcal{M}$ represents a multi-layer CNN by executing approximation to the
proximal mapping in the image space. We propose to first learn a nonlinear
operator $\mathcal{J}$ that combines $\\{\mathbf{b}_{i}\\}$ into the image
$\mathbf{z}=\mathcal{J}(\mathbf{b}_{1},\dots,\mathbf{b}_{c})\in\mathbb{C}^{m\times
n}$ with homogeneous contrast. Then apply a nonlinear operator $\mathcal{G}$
on $\mathbf{z}$ with $\mathcal{G}(\mathbf{z})\in\mathbb{C}^{m\times n\times
N_{f}}$ to extract a $N_{f}$-dimensional features. The nonlinear operator
$\mathcal{J}$ contains four convolutions with an activation function in
between, each convolution obtains kernel size $3\times 3$. The nonlinear
operator $\mathcal{G}$ consists of three convolutions with filter size
$9\times 9$. For the sake of improving the capacity of the network,
$\tilde{\mathcal{J}}$ and $\tilde{\mathcal{G}}$ are employed as adjoint
operators of $\mathcal{J}$ and $\mathcal{G}$ respectively with symmetric
structure and parameters are trained separately. $\mathcal{G}\circ\mathcal{J}$
was designed in the sense of playing a role as an encoder and
$\tilde{\mathcal{J}}\circ\tilde{\mathcal{G}}$ as a decoder. Therefore, image
domain network can be parametrized as compositions of four CNN operators:
$\mathcal{M}=\tilde{\mathcal{J}}\circ\tilde{\mathcal{G}}\circ\mathcal{G}\circ\mathcal{J}$
with output $\mathcal{M}(\mathbf{b}{(t)})\in\mathbb{C}^{m\times n\times c}$.
This step (8b) outputs the multi-coil data
$\bar{\mathbf{u}}=(\bar{\mathbf{u}}_{1},\dots,\bar{\mathbf{u}}_{c})\in\mathbb{C}^{m\times
n\times c}$ from image domain network and we apply the combination operator
$\mathcal{J}$ on $\bar{\mathbf{u}}$ to obtain a full FOV MR single-channel
image $\mathbf{v}=\mathcal{J}(\bar{\mathbf{u}})\in\mathbb{C}^{m\times n}$ that
we desired for reconstruction.

Step (8c) leverages k-space information and further suppresses the high-
frequency artifacts. The output $\bar{\mathbf{u}}{(t)}$ from (8b) is passed to
a k-space domain network by a ResNet structure $\mathbf{u}{(t+1)}=\
\bar{\mathbf{u}}{(t)}+\mathbf{F}^{H}\mathcal{K}\big{(}\mathbf{F}(\bar{\mathbf{u}}{(t)})\big{)}$,
where $\mathbf{F}^{H}\mathcal{K}\mathbf{F}$ is a CNN operator to refine and
further improve the accuracy of the k-space data from each coil. The CNN
$\mathcal{K}$ consists four convolutions, the last convolution kernel numbers
meets the channel number of coil images.

Therefore, (8) is proceed in the following scheme:

$\displaystyle\mathbf{b}_{i}{(t)}=\
\mathbf{u}_{i}{(t)}-\rho_{t}\mathbf{F}^{H}\mathbf{P}^{\top}(\mathbf{P}\mathbf{F}\mathbf{u}_{i}{(t)}-\mathbf{f}_{i}),$
$\displaystyle\quad i=1,\cdots,c,$ (9a)
$\displaystyle\bar{\mathbf{u}}_{i}{(t)}=\
\mathbf{b}_{i}{(t)}+\mathcal{M}(\mathbf{b}_{i}{(t)}),$ $\displaystyle\quad
i=1,\cdots,c,$ (9b) $\displaystyle\mathbf{u}_{i}{(t+1)}=\
\bar{\mathbf{u}}_{i}{(t)}+\mathbf{F}^{H}\mathcal{K}\big{(}\mathbf{F}(\bar{\mathbf{u}}_{i}{(t)})\big{)},$
$\displaystyle\quad i=1,\cdots,c,$ (9c)

for $t=0,\cdots,T-1$. Now we can derive the function $\mathbf{g}$ described in
(4) by combining (9):

$\displaystyle\mathbf{u}_{i}{(t+1)}$
$\displaystyle=(\mathcal{N}-\rho_{t}\mathbf{F}^{H}\mathbf{P}^{\top}\mathbf{P}\mathbf{F}-\mathcal{N}(\rho_{t}\mathbf{F}^{H}\mathbf{P}^{\top}\mathbf{P}\mathbf{F}))\mathbf{u}_{i}{(t)}$
$\displaystyle+\mathbf{u}_{i}{(t)}+\rho_{t}\mathbf{F}^{H}\mathbf{P}^{\top}\mathbf{f}_{i}-\mathcal{N}(\rho_{t}\mathbf{F}^{H}\mathbf{P}^{\top}\mathbf{f}_{i})$
(10a) $\displaystyle=\mathbf{g}(\mathbf{u}_{i}{(t)}),$ (10b)

where
$\mathcal{N}=\mathcal{M}+\mathbf{F}^{H}\mathcal{K}\mathbf{F}+\mathbf{F}^{H}\mathcal{K}\mathbf{F}\mathcal{M}$,
and the function $\mathbf{g}$ maps $\mathbf{u}_{i}{(t)}$ to
$\mathbf{u}_{i}{(t+1)}$ for $t=0,\cdots,T-1$ defined in (10a).

Our proposed reconstruction network is composed of a prescribed $T$ phases,
the initial input of the network is designed as
$\mathbf{u}{(0)}=\mathbf{g}_{0}(\mathbf{f},\theta(0)):=\mathbf{F}^{H}(\mathbf{f}+\mathcal{K}_{0}(\mathbf{f}))$
where $\mathcal{K}_{0}$ is a CNN operator applied to $\mathbf{f}$ in residual
learning to interpolate the missing data and generate a pseudo full k-space.
With the chosen initial $\\{\mathbf{u}_{i}{(0)}\\}$ and partial k-space data
$\\{\mathbf{f}_{i}\\}$ as input, the network performs the update (9) in the
$t$-th phase for $t=0,\cdots,T-1$ and finally the entire network reconstructs
coil-images $\mathbf{u}{(T)}$ and $\mathcal{J}(\bar{\mathbf{u}}{(T)})$, which
is the final single-body image reconstructed as a by-product (complex-valued).

Fig. 1 displays the flowchart of our entire network procedure. The initial
reconstruction suppresses artifacts caused by undersampling. We display the
flowchart of these CNN structures of each phase in Fig. 2. We apply complex-
valued convolutions where multiplications between complex numbers are
performed and use componentwise complex-valued activation function
$\mathbb{C}$ReLU($a+ib$) = ReLU($a$) + $i$ReLU($b$) as suggested in [13]. The
Landweber update step (9a) plays a role in increasing communication between
image space and k-space. The learnable step size $\rho_{t}$ controls the speed
and stability of the convergence. In the image space denoising step (9b), the
operator $\mathcal{J}$ extracts feature across all the channels so that
spatial resolution is improved and tissue details are recovered in the
channel-combined image. The CNN $\mathcal{M}$ carrying channel-integration
$\mathcal{J}$ proceed $T$ times with shared weights in every two phases
therefore the spatial features across channels are learned in an efficient
way. However, the oscillatory artifacts could be misinterpreted as real
features, which might be sharpened. In the k-space denoising update step (9c),
the k-space network $\mathcal{K}$ is compatible with low-frequency
information, so it releases the high-frequency artifacts and recovers the
structure of the image [16, 60]. Therefore iterating the networks
$\mathcal{M}$ and $\mathcal{K}$ in different domains with their individual
effects, the performance compensates and the shortcomings of both networks
offset each other. We evaluate the effect of hybrid domain reconstruction in
ablation studies. Furthermore, the prescribed denoising networks in (9b) and
(9c) refine and update coil-images in each iteration. This iterative procedure
triggers the reconstruction quality of $\mathcal{J}(\bar{\mathbf{u}}{(t)})$
get successively enhancement.

Figure 1: The proposed framework paradigm for all phases. Figure 2: The
proposed framework paradigm at $t+1$-th phase. $W=Re(W)+iIm(W)$ representing
complex convolution filter kernels.

###  3.4 Network training from the view of the Method of Lagrangian
Multipliers (MLM)

The network parameters to be solved from (4) are
$\Theta=\\{\theta(t):t=0,\cdots,T\\}$, where
$\theta(t)=\\{\rho_{t},\mathcal{J}_{t},\mathcal{G}_{t},\tilde{\mathcal{G}}_{t},\tilde{\mathcal{J}}_{t},\mathcal{K}_{t}\\}$
for $t=1,\cdots,T$ and $\theta(0)=\mathcal{K}_{0}$.

The control problem (4) can be solved by using MLM. The corresponding
Lagrangian function is

$\displaystyle L(\textbf{U},\Theta;\Lambda)=$
$\displaystyle\ell(\mathbf{u}{(T)})+\langle\lambda(0),\mathbf{u}{(0)}-\mathbf{g}_{0}\big{(}\mathbf{f},\theta(0)\big{)}\rangle$
(11)
$\displaystyle+\sum\nolimits^{T}_{t=1}\langle\lambda(t),\mathbf{u}{(t)}-\mathbf{g}\big{(}\mathbf{u}{(t-1)},\theta(t)\big{)}\rangle,$

where $\textbf{U}=(\mathbf{u}{(0)},\cdots,\mathbf{u}{(T)})^{\intercal}$ is the
collection of all the states $\mathbf{u}{(t)}$,
$\Lambda=(\lambda(0),\cdots,\lambda(T))^{\intercal}$ are Lagrangian
multipliers of (4). The algorithm proceed in the iterative scheme to update
$\Theta^{k}$, for each training epoch $k=0,\cdots,K$,
$\Theta^{k}=(\theta^{k}(0),\cdots,\theta^{k}(T))^{\intercal}$,
$\textbf{U}^{k}=(\mathbf{u}^{k}{(0)},\cdots,\mathbf{u}^{k}{(T)})^{\intercal}$,
and $\Lambda^{k}=(\lambda^{k}(0),\cdots,\lambda^{k}(T))^{\intercal}$ .

If $(\textbf{U}^{*},\Theta^{*};\Lambda^{*})$ minimizes Lagrange function (11),
by the first order optimality condition, we have

$\displaystyle\partial_{\textbf{U}}L(\textbf{U}^{*},\Theta^{*};\Lambda^{*})=0$
(12a)
$\displaystyle\partial_{\Theta}L(\textbf{U}^{*},\Theta^{*};\Lambda^{*})=0$
(12b)
$\displaystyle\partial_{\Lambda}L(\textbf{U}^{*},\Theta^{*};\Lambda^{*})=0$
(12c)

•Fix $\Theta=\Theta^{k}$, define
$(\textbf{U}^{k},\Lambda^{k}):=\operatorname*{arg\,min}_{\textbf{U},\Lambda}L(\textbf{U},\Theta^{k};\Lambda)$,
then by the first order optimality condition for minimizing $L$ w.r.t
$\lambda(t)$ for $t=0,\cdots,T$, $(\textbf{U}^{k},\Lambda^{k})$ should satisfy

$\displaystyle\partial_{\lambda(0)}\langle\lambda(0),\mathbf{u}{(0)}-\mathbf{g}_{0}(\mathbf{f},\theta^{k}(0))\rangle\Big{|}_{(\mathbf{u}^{k}{(0)},\lambda^{k}(0))}=0$
(13)
$\displaystyle\implies\mathbf{u}^{k}{(0)}=\mathbf{g}_{0}(\mathbf{f},\theta^{k}(0))$
$\displaystyle\partial_{\lambda(t)}\langle\lambda(t),\mathbf{u}{(t)}-\mathbf{g}(\mathbf{u}{(t-1)},\theta^{k}(t))\rangle\Big{|}_{(\mathbf{u}^{k}{(t)},\lambda^{k}(t))}=0$
(14)
$\displaystyle\implies\mathbf{u}^{k}{(t)}=\mathbf{g}(\mathbf{u}^{k}{(t-1)},\theta^{k}(t)),\
t=1,\cdots,T.$

Then by the first order optimality condition for minimizing $L$ w.r.t
$\mathbf{u}{(t)}$, for $t=T$, we get

$\displaystyle\partial_{\mathbf{u}{(T)}}[\ell(\mathbf{u}{(T)})+\langle\lambda(t),\mathbf{u}{(T)}\rangle]\Big{|}_{(\mathbf{u}^{k}{(t)},\lambda^{k}(t))}=0$
(15)
$\displaystyle\implies\lambda^{k}(T)=-\partial_{\mathbf{u}{(T)}}\ell(\mathbf{u}^{k}{(T)}),$

for $t=0,\cdots,T-1$:

$\displaystyle\partial_{\mathbf{u}{(t)}}[\langle\lambda(t),\mathbf{u}{(t)}\rangle-\langle\lambda(t+1),\mathbf{g}(\mathbf{u}{(t)},\theta^{k}{(t+1)})\rangle]\Big{|}_{(\mathbf{u}^{k}{(t)},\lambda^{k}(t))}=0$
$\displaystyle\implies\lambda^{k}(t)=\langle\lambda^{k}(t+1),\partial_{\mathbf{u}{(t)}}\mathbf{g}(\mathbf{u}^{k}{(t)},\theta^{k}{(t+1)})\rangle,$
(16)

•Fix $(\textbf{U}^{k},\Lambda^{k})$ for updating $\Theta$, we compute the
gradient $\partial_{\Theta}L(\textbf{U}^{k},\Theta;\Lambda^{k})$:

$\displaystyle\partial_{\theta(t)}L(\textbf{U}^{k},\theta(t);\Lambda^{k})=\partial_{\theta(t)}[-\langle\lambda^{k}(t),\mathbf{g}(\mathbf{u}^{k}{(t-1)},\theta(t))\rangle]$
$\displaystyle=-\langle\lambda^{k}(t),\partial_{\theta(t)}\
\mathbf{g}(\mathbf{u}^{k}{(t-1)},\theta(t))\rangle,t=1,\cdots,T,$ (17a)
$\displaystyle\partial_{\theta(0)}L(\textbf{U}^{k},\theta(0);\Lambda^{k})=-\langle\lambda^{k}(0),\partial_{\theta(0)}\mathbf{g}_{0}(\mathbf{f},\theta(0))\rangle.$
(17b)

######  Theorem 1.

$\partial_{\Theta}L(\textbf{U}^{k},\Theta;\Lambda^{k})=\partial_{\Theta}\ell(\mathbf{u}^{k}{(T)}(\Theta))$.

###### Proof.

First we show the following holds for $t=0,\cdots,T$:

$\lambda^{k}(t)=-\partial_{\mathbf{u}{(t)}}\ell(\mathbf{u}^{k}{(T)})$ (18)

From (15), (18) is true when $t=T$. Suppose (18) true for
$t=\tau\in\\{1,\cdots,T\\}$. From (16), we have

$\displaystyle\lambda^{k}(\tau-1)$
$\displaystyle=\langle\lambda^{k}(\tau),\partial_{\mathbf{u}{(\tau-1)}}\
\mathbf{g}(\mathbf{u}^{k}{(\tau-1)},\theta^{k}(\tau))\rangle$
$\displaystyle=\langle-\partial_{\mathbf{u}{(\tau)}}\ell(\mathbf{u}^{k}{(T)}),\partial_{\mathbf{u}{(\tau-1)}}\mathbf{u}^{k}{(\tau)}\rangle$
(19a)
$\displaystyle=-\partial_{\mathbf{u}{(\tau-1)}}\ell(\mathbf{u}^{k}{(T)})$
(19b)

Thus, (18) holds for $t=\tau-1$. By the principle of induction, (18) is ture
for all $t=0,\cdots,T$.

Hence, (17a) reduces to

$\displaystyle\partial_{\theta(t)}$ $\displaystyle
L(\textbf{U}^{k},\theta(t);\Lambda^{k})=-\langle\lambda^{k}(t),\partial_{\theta(t)}\
\mathbf{g}(\mathbf{u}^{k}{(t-1)},\theta(t))\rangle$
$\displaystyle=\langle\partial_{\mathbf{u}{(t)}}\ell(\mathbf{u}^{k}{(T)}),\partial_{\theta(t)}\mathbf{u}^{k}{(t)}\rangle$
(20a) $\displaystyle=\partial_{\theta(t)}\ell(\mathbf{u}^{k}{(T)}),\
t=1,\cdots,T.$ (20b)

Also (17b) gives
$\partial_{\theta(0)}L(\textbf{U}^{k},\theta(0);\Lambda^{k})=\partial_{\theta(0)}\ell(\mathbf{u}^{k}{(T)})$.
Therefore, we derive
$\partial_{\Theta}L(\textbf{U}^{k},\Theta;\Lambda^{k})=\partial_{\Theta}\ell(\mathbf{u}^{k}{(T)}(\Theta))$.
∎

This theorem further implies that: Since the gradients are the same, applying
SGD Algorithms or its variance such as Adam [29] to minimize loss function
$\ell$ is equivalent to perform the same algorithm on $L$. Network training
algorithm using MLM can be summarized in Algorithm 1.

1 Hyperparameter: $K$ (#Iterations)

2 Initialize: Initial guess $\theta^{0}(t)\in\Theta^{0},t=0,\cdots,T$

3 for  _$k=0$ to $K-1$_ do

4 Set $\mathbf{u}^{k}{(0)}=\mathbf{g}_{0}(\mathbf{f},\theta(0))$

5 for  _$t=1$ to $T$_ do

6 $\mathbf{u}^{k}{(t)}=\mathbf{g}(\mathbf{u}^{k}{(t-1)},\theta^{k}(t))$

7 end for

8 Set $\lambda^{k}(T)=-\partial_{\mathbf{u}}\ell(\mathbf{u}^{k}{(T)})$

9 for  _$t=T-1$ to $0$ _ do

10
$\lambda^{k}(t)=\langle\lambda^{k}{(t+1)},\partial_{\mathbf{u}}\mathbf{g}(\mathbf{u}^{k}{(t)},\theta^{k}(t+1))\rangle$

11 end for

12 for  _$t=0$ to $T$ _ do

13 Set
$\theta^{k+1}(t)=\operatorname*{arg\,min}_{\theta(t)}L(\textbf{U}^{k},\theta(t);\Lambda^{k})$.

14 end for

15

16 end for

Output: $\theta^{K}(t),t=0,\cdots,T$.

Algorithm 1 Network Training by MLM

##  4 Experimental Results

###  4.1 Data Set

The data in our experiments was acquired by a 15-channel knee coil array with
two pulse sequences: a proton density weighting with (FSPD) and without (PD)
fat suppression in the coronal direction from https://github.com/VLOGroup/mri-
variationalnetwork. We used a regular Cartesian sampling mask with 31.56%
sampling ratio as shown in the lower-right of Fig. 4. Each of the two
sequences data includes images of 20 patients, we select 27-28 central image
slices from 19 patients, which amount to 526 images each of size $320\times
320$ as the training dataset, and 15 central image slices are picked from one
patient that is not included in the training data set as testing dataset.

###  4.2 Implementation

We evaluate classical methods GRAPPA [18], SPIRiT [37], and deep learning
methods VN [21], De-AliasingNet [9], DeepcomplexMRI [58] and Adaptive-CS-Net
[44] over the 15 testing Coronal FSPD and Coronal PD knee images with regular
Cartesian sampling in terms of PSNR, structural similarity (SSIM) [59] and
relative error RMSE. The following equations are computations of SSIM, PSNR
and RMSE between reconstruction $\mathbf{v}=|\mathcal{J}(\bar{\mathbf{u}})|$
and ground truth single-body image $\mathbf{v}^{*}$:

$SSIM=\frac{(2\mu_{\mathbf{v}}\mu_{\mathbf{v}^{*}}+C_{1})(2\sigma_{\mathbf{v}\mathbf{v}^{*}}+C_{2})}{(\mu_{\mathbf{v}}^{2}+\mu_{\mathbf{v}^{*}}^{2}+C_{1})(\sigma_{\mathbf{v}}^{2}+\sigma_{\mathbf{v}^{*}}^{2}+C_{2})},$
(21)

where $\mu_{\mathbf{v}},\mu_{\mathbf{v}^{*}}$ are local means of pixel
intensity, $\sigma_{\mathbf{v}},\sigma_{\mathbf{v}^{*}}$ denote the standard
deviation and $\sigma_{\mathbf{v}\mathbf{v}^{*}}$ is covariance between
$\mathbf{v}$ and $\mathbf{v}^{*}$, $C_{1}=(k_{1}L)^{2},C_{2}=(k_{2}L)^{2}$ are
two constants that avoid denominator to be zero, and $k_{1}=0.01,k_{2}=0.03$.
$L$ is the largest pixel value of the magnitude of coil-images.

$PSNR=20\log_{10}\big{(}\max(\lvert\mathbf{v}^{*}\rvert)\big{/}\frac{1}{N}\|\mathbf{v}^{*}-|\mathcal{J}(\bar{\mathbf{u}})|\|^{2}\big{)},$
(22)

where $N$ is the total number of pixels in the magnitude of ground truth.

$RMSE=\|\mathbf{v}^{*}-|\mathcal{J}(\bar{\mathbf{u}})|\|/\|\mathbf{v}^{*}\|.$
(23)

The relative error between the multi-coil reconstruction $\mathbf{u}$ and the
ground truth $\mathbf{u}^{*}$ is defined as

$RMSE=\sqrt{\sum\nolimits^{c}_{i=1}\|\mathbf{u}^{*}_{i}-\mathbf{u}_{i}\|^{2}/\sum\nolimits^{c}_{i=1}\|\mathbf{u}^{*}_{i}\|^{2}}.$
(24)

All the experiments are implemented and tested in TensorFlow [1] on a Windows
workstation with Intel Core i9 CPU at 3.3GHz and an Nvidia GTX-1080Ti GPU with
11GB of graphics card memory. The parameters in proposed networks are using
Xavier initialization [17] to initialize $\theta^{0}(t),t=0,\cdots,T$. Indeed,
solving line 13 in the Algorithm 1 using any stochastic gradient descent
algorithm such as Adam [29] has the same performance as minimizing loss
function w.r.t $\theta(t)$ using the same algorithm to update
$\theta^{k+1}(t)$, because of the equivalence of the gradients:
$\partial_{\Theta}L(\textbf{U}^{k},\Theta;\Lambda^{k})=\partial_{\Theta}\ell(\mathbf{u}^{k}{(T)}(\Theta))$,
which is proved by Theorem 1. TensorFlow provides optimized APIs for automatic
differentiation which helps to build highly performant input pipelines and
keeps high GPU utilization. In order to implement a more stable gradient
calculation, we replace line 13 by minimizing $\ell$ with the Adam algorithm.
The network was trained with total epochs $K=700$ to update $\Theta^{k}$. We
apply exponentially decay learning rate 0.0001,
$\beta_{1}=0.9,\beta_{2}=0.999,\epsilon=10^{-8}$ and mini-batch size of 2 is
used in Adam optimizer. The initial step size is set to $\rho_{0}=1$ for both
real and imaginary parts and we choose $\gamma=10^{-3},\eta=10^{-4}$ in (5).
The proposed network was implemented with $T=4$ and parameters trained in the
network $\mathcal{M}$ are shared for every two phases.

###  4.3 Comparison with existing methods

The average numerical performance with standard deviations of the proposed
method and several state-of-the-art methods are summarized in Table 3. The
proposed method achieves the best reconstruction quality in terms of
PSNR/SSIM/RMSE. Our method is parameter efficient due to the network
$\mathcal{M}$ share parameters in every two phases so that the entire network
reduces more than 1/4 of learnable parameters. We listed the trained parameter
numbers and inference time in table 3. Fig. 3 and Table 2 indicate that
reconstruction performance get improved progressively as $t$ increases.

Figure 3:  Reconstructed channel-combined images
$\mathcal{J}(\bar{\mathbf{u}}{(t)})$ for $t=1,2,3,4$ of proposed method. Table
2: Tested Average PSNR for FSPD data and PD data for each phase.

Phase number | PSNR of FSPD data | PSNR of PD data
---|---|---
1 | 16.2007 | 15.5577
2 | 33.7201 | 40.5228
3 | 36.9523 | 42.6623
4 | 40.7101 | 44.8120

GRAPPA and SPIRiT adopted calibration kernel size $5\times 5$ and Tikhonov
regularization in the calibration was set to be 0.01. Tikhonov regularization
in the reconstruction was set as $10^{-3}$ for implementing SPIRiT, which took
30 iterations. In training and testing of VN and DeepcomplexMRI, the network
and parameter settings are used as stated in their paper and code. We modified
De-AliasingNet by increasing filter numbers to be 64 in 20 iterations to
improve the performance for fair competition. In the CDF-Net, we set 128
channels as an initial layer with a depth of 4 in all the three Frequency-
Informed U-Nets. We choose $\sigma$ to be Softplus activation function. We
implemented Adaptive-CS-Net with a total of 15 reconstruction blocks including
2 of 16 filters with size $3\times 3$, 3 of 32 filters with size $5\times 5$,
and 2 of 64 filters with size $5\times 5$, the blocks are mixed with 2 and 3
scales, which was controlled as the largest tolerance of our GPU capacity. We
input the data consistency as prior knowledge of training data to the network,
neighboring slices with the center slice are used as input. For fair
competitions, CDF-Net and Adaptive-CS-Net have employed complex convolutions
and activation $\mathbb{C}$ReLU.

VN applied precalculated coil sensitivities, and output full-body image. De-
AliasingNet, DeepcomplexMRI, CDF-Net, and Adaptive-CS-Net both output multi-
coil images and do not require coil sensitivity maps. DeepcomplexMRI uses the
adaptive coil combination method [54], the other networks all applied RSS.
CDF-Net and proposed network perform cross-domain reconstruction, other
learning-based methods perform on image domain. Comparison between referenced
pMRI methods and proposed methods are shown in Fig. 4 for PD images and FSPD
images are in Fig. 5. We observe the evidence that deep learning-based methods
significantly outperform classical methods GRAPPA and SPIRiT in reconstruction
accuracy. In learning-based methods, the tested images from VN and De-
AliasingNet are more blurry than other ones and lost sharp details in some
complicated tissue, other methods display only slight differences in the
detail.

Table 3: Quantitative evaluations of the reconstructions on the Coronal FSPD &
PD data and reconstruction time for each of the referenced methods. Time (in
seconds) refers to the testing time for each method.

|  | FSPD data |  |  | PD data |  |  |
---|---|---|---|---|---|---|---|---
Method | PSNR | SSIM | RMSE | PSNR | SSIM | RMSE | Time | Parameters
GRAPPA [18] | 24.9251$\pm$0.9341 | 0.4827$\pm$0.0344 | 0.2384$\pm$0.0175 | 30.4154$\pm$0.5924 | 0.7489$\pm$0.0207 | 0.0984$\pm$0.0030 | 280s | N/A
SPIRiT [37] | 28.3525$\pm$1.3314 | 0.6509$\pm$0.0300 | 0.1614$\pm$0.0203 | 32.0011$\pm$0.7920 | 0.7979$\pm$0.0306 | 0.0824$\pm$0.0082 | 43s | N/A
VN [21] | 30.2588$\pm$1.1790 | 0.7141$\pm$0.0483 | 0.1358$\pm$0.0152 | 37.8265$\pm$0.4000 | 0.9281$\pm$0.0114 | 0.0422$\pm$0.0036 | 0.16s | 0.13 M
De-AliasingNet [9] | 36.1017$\pm$1.2981 | 0.8941$\pm$0.0269 | 0.0697$\pm$0.0119 | 41.2151$\pm$0.7872 | 0.9711$\pm$0.0033 | 0.0285$\pm$0.0015 | 0.92s | 0.34 M
DeepcomplexMRI [58] | 36.5706$\pm$1.0215 | 0.9008$\pm$0.0190 | 0.0654$\pm$0.0049 | 41.5756$\pm$0.6271 | 0.9679$\pm$0.0031 | 0.0274$\pm$0.0018 | 1.04s | 1.45 M
CDF-Net [41] | 40.0405$\pm$1.5518 | 0.9520$\pm$0.0182 | 0.0443$\pm$0.0070 | 43.7943$\pm$1.9888 | 0.9879$\pm$0.0026 | 0.0215$\pm$0.0042 | 0.47s | 3.44 M
Adaptive-CS-Net [44] | 40.1846$\pm$1.4780 | 0.9534$\pm$0.0175 | 0.0435$\pm$0.0065 | 44.1131$\pm$1.5596 | 0.9878$\pm$0.0026 | 0.0206$\pm$0.0027 | 1.57s | 5.59 M
Proposed | 40.7101$\pm$1.5357 | 0.9619$\pm$0.0144 | 0.0408$\pm$0.0051 | 44.8120$\pm$1.3185 | 0.9886$\pm$0.0023 | 0.0189$\pm$0.0018 | 0.52s | 2.92 M

Figure 4: Qualitative comparison results of reconstruction methods on the
Coronal PD knee image. The top row shows reconstructed images and the
referenced image, the second and third-row are corresponding zoomed-in ROIs of
the red box area (draw on the rightmost ground truth image), the fourth row
shows corresponding pointwise absolute error maps and the last row shows
corresponding values and regular Cartesian sampling (31.56% rate) mask.

Figure 5: Qualitative comparison results of reconstruction methods on the Coronal FSPD knee image. From top to bottom: reconstructed images of different methods and the reference image, corresponding zoomed-in ROIs, corresponding pointwise absolute error maps and color bar, corresponding evaluation metric values, and regular Cartesian sampling (31.56% rate) mask.  Table 4: Mean square error between $\mathbf{u}$ and $\mathbf{u}^{*}$ for PD data comparing to the methods that reconstruct multi-coil images.  Method | RMSE between $\mathbf{u}$ and $\mathbf{u}^{*}$
---|---
De-AliasingNet [9] | 0.1573
DeepcomplexMRI [58] | 0.1253
CDF-Net [41] | 0.1196
Adaptive-CS-Net [44] | 0.1096
pMRI-$\mathbb{C}$Net-K (26) | 0.0538
Proposed (5) | 0.0505

###  4.4 Ablation Studies

The experiments introduced in this section implemented only in the image
domain, in which (9) is replaced by:

$\displaystyle\mathbf{b}_{i}{(t)}$
$\displaystyle=\mathbf{u}_{i}{(t)}-\rho_{t}\mathbf{F}^{H}\mathbf{P}^{\top}(\mathbf{P}\mathbf{F}\mathbf{u}_{i}{(t)}-\mathbf{f}_{i}),$
(25a) $\displaystyle\mathbf{u}_{i}{(t+1)}$
$\displaystyle=\mathbf{b}_{i}{(t)}+\tilde{\mathcal{J}}\circ\tilde{\mathcal{G}}\circ\mathcal{S}_{\alpha_{t}}(\mathcal{G}\circ\mathcal{J}(\mathbf{b}_{i}{(t)})).$
(25b)

$\mathcal{S}_{\alpha_{t}}$ represents soft shrinkage operator, we set
$\alpha_{0}=0$ for both real and imaginary parts. Details for this model was
explained in the conference paper [5]. The proximal operator is learned in
residual update (25b) and only performs in the image domain.

We conduct a series of experiments to test the effects of several important
components in the network proposed. Table 6 displays the comparison of
proposed ablation study. All the experiments in Table 6 were implemented by
(25) except for the last one “proposed” method which unrolled (9).

Table 5: Labels of the variations of the proposed pMRI reconstruction network. Label | Meaning
---|---
-Net | Real-valued convolution/activation
-$\mathbb{C}$Net | Complex-valued convolution/activation
-RSS | Using root of sum of squares in place of $\mathcal{J}$
-ZF | Zero-filling as initial $\mathbf{u}{(0)}$
-SP | SPIRiT as initial $\mathbf{u}{(0)}$
-K |  $\mathbf{F}^{H}(\mathbf{f}+\mathcal{K}(\mathbf{f}))$ as initial $\mathbf{u}{(0)}$

To specify which of the components are employed, we append the labels shown in
Table 5 with different variations of the ablated pMRI networks. For instance,
the network with real-valued convolution and activation function with zero-
filled initialization is denoted by pMRI-Net-ZF, etc.

Table 6: Quantitative evaluations of the reconstructions on the Coronal FSPD &
PD data using different variations of the proposed methods (Labels of
variations are explained in Tabel 5. The experiments without being labeled
loss functions are trained with (27)).

|  | FSPD data |  |  | PD data |  |  |
---|---|---|---|---|---|---|---|---
Method | PSNR | SSIM | RMSE | PSNR | SSIM | RMSE | Phases T | Parameters
pMRI-$\mathbb{C}$Net-RSS | 36.4887$\pm$0.9787 | 0.9002$\pm$0.0197 | 0.0661$\pm$0.0051 | 41.2897$\pm$0.8430 | 0.9281$\pm$0.0357 | 0.0285$\pm$0.0037 | 5 | 5.03 M
pMRI-Net-ZF | 37.8475$\pm$1.2086 | 0.9212$\pm$0.0236 | 0.0568$\pm$0.0069 | 42.4333$\pm$0.8785 | 0.9793$\pm$0.0023 | 0.0249$\pm$0.0024 | 5 | 5.03 M
pMRI-Net-SP | 38.0205$\pm$0.8125 | 0.9291$\pm$0.0183 | 0.0555$\pm$0.0057 | 42.7435$\pm$0.4856 | 0.9754$\pm$0.0047 | 0.0239$\pm$0.0019 | 5 | 5.03 M
pMRI-$\mathbb{C}$Net-ZF | 38.1157$\pm$1.3776 | 0.9277$\pm$0.0257 | 0.0552$\pm$0.0085 | 42.7859$\pm$1.1285 | 0.9818$\pm$0.0026 | 0.0241$\pm$0.0045 | 5 | 5.03 M
pMRI-$\mathbb{C}$Net-SP | 38.3239$\pm$1.1305 | 0.9282$\pm$0.0269 | 0.0539$\pm$0.0075 | 42.8924$\pm$0.9336 | 0.9760$\pm$0.0054 | 0.0237$\pm$0.0034 | 5 | 5.03 M
pMRI-Net-K | 38.8717$\pm$1.1330 | 0.9389$\pm$0.0209 | 0.0504$\pm$0.0057 | 42.9060$\pm$0.8765 | 0.9802$\pm$0.0028 | 0.0236$\pm$0.0028 | 4 | 4.11 M
pMRI-$\mathbb{C}$Net-K (26) | 38.9661$\pm$1.4382 | 0.9421$\pm$0.0177 | 0.0498$\pm$0.0056 | 43.2604$\pm$0.7610 | 0.9833$\pm$0.0022 | 0.0226$\pm$0.0022 | 4 | 4.11 M
pMRI-$\mathbb{C}$Net-K | 39.3360$\pm$1.1854 | 0.9497$\pm$0.0208 | 0.0477$\pm$0.0048 | 43.5653$\pm$0.8265 | 0.9844$\pm$0.0022 | 0.0217$\pm$0.0010 | 4 | 4.11 M
Proposed (5) | 40.7101$\pm$1.5357 | 0.9619$\pm$0.0144 | 0.0408$\pm$0.0051 | 44.8120$\pm$1.3185 | 0.9886$\pm$0.0023 | 0.0189$\pm$0.0018 | 4 | 2.92 M

The complete structure of the network with the three types of initialization
is shown in Fig. 6.

Figure 6: The framework paradigm for all phases with three different initial
reconstruction methods, including zero-filled initial, SPIRiT reconstruction
initial, and the learned initial. For $t=1,\cdots,T$, each phase follows the
algorithm introduced in our previous work [5].

We consider two types of training datasets in the ablation experiments and
design a proper loss function for each of them in Section 4.4. In the first
case, the ground truth is $\mathbf{u}^{*}$, we set loss function:

$\displaystyle\ell(\mathbf{u})=$
$\displaystyle\sum\nolimits^{c}_{i=1}\gamma\|\mathbf{u}_{i}-\mathbf{u}^{*}_{i}\|+\||\mathcal{J}(\mathbf{u})|-\text{RSS}(\mathbf{u}^{*})\|$
(26)
$\displaystyle+\beta\|\text{RSS}(\mathbf{u}_{i}{(0)})-\text{RSS}(\mathbf{u}^{*}_{i})\|,$

$\mathbf{u}$ means $\mathbf{u}{(T)}$ for simplicity. If the training dataset
consists ground truth single-body image $\mathbf{v}^{*}\in\mathbb{R}^{m\times
n}=\text{RSS}(\mathbf{u}^{*})$, then we test our network performance using the
loss function:

$\ell(\mathbf{u})=\gamma\|\text{RSS}(\mathbf{u})-\mathbf{v}^{*}\|+\||\mathcal{J}(\mathbf{u})|-\mathbf{v}^{*}\|+\beta\|\text{RSS}(\mathbf{u}{(0)})-\mathbf{v}^{*}\|.$
(27)

The loss functions are indicated in Tables 6 and 4 follow the experiments. We
set $\gamma=1,\beta=10^{-3}$ in (26) and (27) in the implementations.

####  4.4.1 Combination operator $\mathcal{J}$ vs. RSS

In order to justify the effectiveness of the learned nonlinear combination
operator $\mathcal{J}$, we modified pMRI-$\mathbb{C}$Net-ZF by substituting
$\mathcal{J}$ with RSS which is widely used to combine multi-coil images into
a single-body image. The other operators $\mathcal{G},\tilde{\mathcal{G}}$ and
$\tilde{\mathcal{J}}$ remain to perform complex convolutions. Specifically,
the output of RSS:
$(\sum_{i=1}^{c}|\mathbf{u}_{i}{(t)}|^{2})^{1/2}\in\mathbb{R}^{mn}$ is a
single-channel real-valued image, which is input into both real and imaginary
parts of the nonlinear operator $\mathcal{G}$, so the output split into
complex value. We refer pMRI-Net-ZF/pMRI-$\mathbb{C}$Net-ZF as the networks
with combination operator $\mathcal{J}$ and pMRI-$\mathbb{C}$Net-RSS as the
network with RSS, all other settings remain the same as before.

The reconstructed images and average evaluation results of these two types of
networks are shown in Fig. 5 and Table 6. pMRI-Net-ZF and
pMRI-$\mathbb{C}$Net-ZF outperform pMRI-$\mathbb{C}$Net-RSS with mean
improvements of 1.36 dB and 1.63 dB in PSNR respectively. It indicates that
the learned combination operator $\mathcal{J}$ gives more favorable
performance compared with applying RSS.

####  4.4.2 Effect of initialization

The choice of the input of the reconstruction network $\mathbf{u}{(0)}$, also
has impacts on the final reconstruction quality. Instead of directly using the
partial k-space data $\mathbf{f}$ as the input of our network, we use three
different choices of input $\mathbf{u}{(0)}$: (i) Zero-filling reconstruction
$\mathbf{F}^{H}\mathbf{f}$; (ii) SPIRiT [37] reconstruction; (iii)
$\mathbf{F}^{H}(\mathbf{f}+\mathcal{K}(\mathbf{f}))$, one can treat
$\mathbf{f}+\mathcal{K}(\mathbf{f})$ as an interpolated pseudo full k-space.

We observe that from Table 6, SPIRiT initial makes a slight improvement over
the zero-filled initial, whereas the learned initial achieves the highest
reconstruction quality compared to the other two initializations. Fig. 7
displays the three types of initials. We observe that both SPIRiT and the
learned initial obtain higher spatial resolution over zero-filling. SPIRiT is
a classical k-space method, this initial did a better job on reducing the
aliasing artifacts and keeping edges compare to zero-filling and the learned
initial, but SPIRiT introduces more noise.

Comparing to zero-filling, the learned initial preserves structure features of
major tissue thanks to the k-space network $\mathcal{K}$. Comparing to the
SPIRiT initial, the learned initial reduces resolution noises in the image
space. Learning-based initial obtain a balanced performance between zero-
filling and SPIRiT in the sense of avoiding the weakness of these two
initials.

Figure 7: The first row (from left to right) shows the RSS of initial
$\bar{\mathbf{u}}{(0)}$ for pMRI-Net-ZF/pMRI-$\mathbb{C}$Net-ZF, pMRI-Net-
SP/pMRI-$\mathbb{C}$Net-SP, pMRI-Net-K/pMRI-$\mathbb{C}$Net-K, and reference
image on the Coronal PD knee image. The second row shows their pointwise error
maps and color bar.

####  4.4.3 Complex convolutions

In this experiment, we compared -Net and -$\mathbb{C}$Net. The quantitative
evaluation from Table 6 indicates that the proposed complex-valued networks
are outstanding in terms of PSNR/SSIM/RMSE over proposed real-valued networks.
Table 4 informs the complex convolutions extract the features in each
$\mathbf{u}_{i}$ and obtain the lowest RMSE. The phase image for one channel
of the reconstructed coil-image is displayed in Fig. 3. These results
demonstrate complex-valued networks are playing important roles in updating
multi-coil images and preserving phase information of each channel.

Figure 8: The first row shows phase information of one coil in the
reconstructed coil-image for GRAPPA, SPIRiT, De-AliasingNet, DeepcomplexMRI,
CDF-Net, Adaptive-CS-Net, pMRI-$\mathbb{C}$Net-K with loss function (26), the
proposed method and referenced image. The second row shows the corresponding
pointwise error maps and color bar, the maximum error is $30^{\circ}$.

####  4.4.4 Comparison with proposed image domain reconstruction and domain
hybrid reconstruction

The major difference between pMRI-$\mathbb{C}$Net-K and the proposed network
is that pMRI-$\mathbb{C}$Net-K only iterates (25b) in the image domain, while
the proposed method iterates (9) which conducts a domain hybrid
reconstruction. The denoising network in (25b) is
$\tilde{\mathcal{J}}\circ\tilde{\mathcal{G}}\circ\mathcal{S}_{\alpha_{t}}(\mathcal{G}\circ\mathcal{J})$,
and in the proposed network we use
$\mathcal{M}+\mathbf{F}^{H}\mathcal{K}\mathbf{F}+\mathbf{F}^{H}\mathcal{K}\mathbf{F}\mathcal{M}$,
where
$\mathcal{M}=\tilde{\mathcal{J}}\circ\tilde{\mathcal{G}}\circ\mathcal{G}\circ\mathcal{J}$
is the image domain network. The soft-thresholding operator
$\mathcal{S}_{\alpha_{t}}$ was eliminated in the proposed method since we
found the results of adding the soft-thresholding does not make an obvious
difference.

The average reconstruction outcomes in Table 6 and Table 4 suggest that domain
hybrid approach achieves better performance. Comparing to the
pMRI-$\mathbb{C}$Net-K, our proposed method improved 0.598 dB in PSNR, 0.003
in SSIM, and reduced 0.014 in RMSE, which shows in Table 6 for PD dataset.

##  5 Conclusion

This paper introduces a discrete-time optimal control framework for the
calibration-free pMRI reconstruction model. We apply a convolutional
combination operator to combine channels of the multi-coil images and apply a
parametrized regularization function to the channel-combined image to
reconstruct channel-wise multi-coil images. The proposed method is inspired by
the proximal gradient algorithm. The proximal point is learned by two
denoising networks, which conducts in the image domain and k-space domain. We
cast the reconstruction network as a structured discrete-time optimal control
system, resulting in an optimal control formulation of parameter training,
which provides an interpretable and high-performance deep architecture for
pMRI reconstruction. We design network training from the view of the Method of
Lagrangian Multipliers. We showed that the method of Lagrangian multipliers is
equivalent to back-propagation, and we can employ SGD based algorithms to
obtain a solution satisfying the necessary condition of the optimal control
problem. The reconstruction results are of high perceived quality demonstrate
the superior performance of the proposed pMRI-Net.

## References

  * Abadi et al. [2016] Abadi, M., et al., 2016.  Tensorflow: A system for large-scale machine learning, in: 12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation ($\\{$OSDI$\\}$ 16), pp. 265–283.
  * Adler and Öktem [2018] Adler, J., Öktem, O., 2018\.  Learned primal-dual reconstruction.  IEEE transactions on medical imaging 37, 1322–1332.
  * Aggarwal et al. [2019] Aggarwal, H.K., Mani, M.P., Jacob, M., 2019.  Modl: Model-based deep learning architecture for inverse problems.  IEEE Transactions on Medical Imaging 38, 394–405.
  * Bengio et al. [1994] Bengio, Y., Simard, P., Frasconi, P., 1994.  Learning long-term dependencies with gradient descent is difficult.  IEEE Transactions on Neural Networks 5, 157–166.  doi:10.1109/72.279181.
  * Bian et al. [2020] Bian, W., Chen, Y., Ye, X., 2020.  Deep parallel mri reconstruction network without coil sensitivities, in: Deeba, F., Johnson, P., Würfl, T., Ye, J.C. (Eds.), Machine Learning for Medical Image Reconstruction, Springer International Publishing, Cham. pp. 17–26.
  * Buonocore and Gao [1997] Buonocore, M.H., Gao, L., 1997\.  Ghost artifact reduction for echo planar imaging using image phase correction.  Magnetic resonance in medicine 38, 89–100.
  * Chen et al. [2020] Chen, E.Z., Chen, T., Sun, S., 2020.  Mri image reconstruction via learning optimization using neural odes, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 83–93.
  * Chen et al. [2018] Chen, R.T., Rubanova, Y., Bettencourt, J., Duvenaud, D., 2018\.  Neural ordinary differential equations.  arXiv preprint arXiv:1806.07366 .
  * Chen et al. [2019] Chen, Y., et al., 2019.  Model-based convolutional de-aliasing network learning for parallel mr imaging, in: Medical Image Computing and Computer Assisted Intervention – MICCAI 2019, Springer International Publishing, Cham. pp. 30–38.
  * Cheng et al. [2019] Cheng, J., et al., 2019.  Model learning: Primal dual networks for fast mr imaging, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 21–29.
  * Chernousko and Lyubushin [1982] Chernousko, F.L., Lyubushin, A.A., 1982\.  Method of successive approximations for solution of optimal control problems.  Optimal Control Applications and Methods 3, 101–114.
  * Choi et al. [2013] Choi, J., Kim, D., Oh, C., Han, Y., Park, H., 2013\.  An iterative reconstruction method of complex images using expectation maximization for radial parallel mri.  Physics in Medicine & Biology 58, 2969.
  * Cole et al. [2020] Cole, E.K., et al., 2020.  Analysis of deep complex-valued convolutional neural networks for mri reconstruction.  arXiv:2004.01738 URL: https://arxiv.org/abs/2004.01738.
  * Constantinides et al. [1997] Constantinides, C.D., Atalar, E., McVeigh, E.R., 1997.  Signal-to-noise measurements in magnitude images from nmr phased arrays.  Magnetic Resonance in Medicine 38, 852–857.  doi:10.1002/mrm.1910380524.
  * Deshmane et al. [2012] Deshmane, A., et al., 2012.  Parallel mr imaging.  Journal of Magnetic Resonance Imaging 36, 55–72.
  * Eo et al. [2018] Eo, T., et al., 2018.  Kiki-net: cross-domain convolutional neural networks for reconstructing undersampled magnetic resonance images.  Magnetic Resonance in Medicine 80, 2188–2201.  doi:10.1002/mrm.27201.
  * Glorot and Bengio [2010] Glorot, X., Bengio, Y., 2010\.  Understanding the difficulty of training deep feedforward neural networks, in: Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 249–256.
  * Griswold et al. [2002] Griswold, M.A., et al., 2002.  Generalized autocalibrating partially parallel acquisitions (grappa).  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 47, 1202–1210.
  * Halkin [1966] Halkin, H., 1966.  A maximum principle of the pontryagin type for systems described by nonlinear difference equations.  SIAM Journal on Control 4, 90–111.  doi:10.1137/0304009.
  * Hammernik et al. [2021] Hammernik, K., Schlemper, J., Qin, C., Duan, J., Summers, R.M., Rueckert, D., 2021\.  Systematic evaluation of iterative deep neural networks for fast parallel mri reconstruction with sensitivity-weighted coil combination.  Magnetic Resonance in Medicine .
  * Hammernik et al. [2018] Hammernik, K., et al., 2018.  Learning a variational network for reconstruction of accelerated mri data.  Magnetic Resonance in Medicine 79, 3055–3071.
  * He et al. [2016] He, K., Zhang, X., Ren, S., Sun, J., 2016\.  Deep residual learning for image recognition, in: 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 770–778.  doi:10.1109/CVPR.2016.90.
  * Hochreiter [1998] Hochreiter, S., 1998.  The vanishing gradient problem during learning recurrent neural nets and problem solutions.  International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems 06, 107–116.  doi:10.1142/S0218488598000094.
  * Huang [2020] Huang, F., 2020.  Medical imaging using neural networks.  US Patent App. 16/904,981.
  * Huang and Chen [2020] Huang, F., Chen, M., 2020.  Magnetic resonance imaging method and device.  US Patent 10,852,376.
  * Huang et al. [2020] Huang, F., Han, D., Mei, L., 2020.  Magnetic resonance imaging with deep neutral networks.  US Patent App. 16/735,874.
  * Islam et al. [2021] Islam, S.R., Maity, S.P., Ray, A.K., 2021.  Compressed sensing regularized calibrationless parallel magnetic resonance imaging via deep learning.  Biomedical Signal Processing and Control 66, 102399.
  * Jiang et al. [2020] Jiang, J., et al., 2020.  Unified cross-modality feature disentangler for unsupervised multi-domain mri abdomen organs segmentation, in: Medical Image Computing and Computer Assisted Intervention – MICCAI 2020, Springer International Publishing, Cham. pp. 347–358.
  * Kingma and Ba [2015] Kingma, D.P., Ba, J., 2015.  Adam: A method for stochastic optimization URL: http://arxiv.org/abs/1412.6980.
  * Knoll et al. [2020] Knoll, F., Hammernik, K., Zhang, C., Moeller, S., Pock, T., Sodickson, D.K., Akcakaya, M., 2020.  Deep-learning methods for parallel magnetic resonance imaging reconstruction: A survey of the current approaches, trends, and issues.  IEEE signal processing magazine 37, 128–140.
  * LeCun et al. [2015] LeCun, Y., Bengio, Y., Hinton, G., 2015.  Deep learning.  nature 521, 436–444.
  * [32] Li, Q., .  Dynamical systems and machine learning .
  * Li et al. [2017] Li, Q., Chen, L., Tai, C., et al., 2017.  Maximum principle based algorithms for deep learning.  arXiv preprint arXiv:1710.09513 .
  * Li and Hao [2018] Li, Q., Hao, S., 2018.  An optimal control approach to deep learning and applications to discrete-weight neural networks, in: Dy, J., Krause, A. (Eds.), Proceedings of the 35th International Conference on Machine Learning, PMLR. pp. 2985–2994.
  * Li et al. [2019] Li, Q., Lin, T., Shen, Z., 2019.  Deep learning via dynamical systems: An approximation perspective.  arXiv preprint arXiv:1912.10382 .
  * Lu et al. [2020] Lu, T., Zhang, X., Huang, Y., Guo, D., Huang, F., Xu, Q., Hu, Y., Ou-Yang, L., Lin, J., Yan, Z., et al., 2020\.  pfista-sense-resnet for parallel mri reconstruction.  Journal of Magnetic Resonance 318, 106790.
  * Lustig and Pauly [2010] Lustig, M., Pauly, J.M., 2010\.  Spirit: Iterative self-consistent parallel imaging reconstruction from arbitrary k-space.  Magnetic Resonance in Medicine 64, 457–471.  doi:10.1002/mrm.22428.
  * Lv et al. [2021] Lv, J., Wang, C., Yang, G., 2021.  Pic-gan: A parallel imaging coupled generative adversarial network for accelerated multi-channel mri reconstruction.  Diagnostics 11, 61\.
  * Mardani et al. [2019] Mardani, M., et al., 2019.  Deep generative adversarial neural networks for compressive sensing mri.  IEEE Transactions on Medical Imaging 38, 167–179.
  * Meng et al. [2019] Meng, N., et al., 2019.  A prior learning network for joint image and sensitivity estimation in parallel mr imaging, in: Medical Image Computing and Computer Assisted Intervention – MICCAI 2019, Springer International Publishing, Cham. pp. 732–740.
  * Nitski et al. [2020] Nitski, O., et al., 2020.  Cdf-net: Cross-domain fusion network for accelerated mri reconstruction, in: Medical Image Computing and Computer Assisted Intervention – MICCAI 2020, Springer International Publishing, Cham. pp. 421–430.
  * Parikh and Boyd [2014] Parikh, N., Boyd, S., 2014.  Proximal algorithms.  Foundations and Trends in optimization 1, 127–239.
  * Pelc et al. [1991] Pelc, N.J., Herfkens, R.J., Shimakawa, A., Enzmann, D.R., et al., 1991\.  Phase contrast cine magnetic resonance imaging.  Magnetic resonance quarterly 7, 229–254.
  * Pezzotti et al. [2020] Pezzotti, N., et al., 2020.  An adaptive intelligence algorithm for undersampled knee mri reconstruction: Application to the 2019 fastmri challenge.  arXiv:2004.07339 .
  * Pruessmann et al. [1999] Pruessmann, K.P., et al., 1999.  Sense: sensitivity encoding for fast mri.  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 42, 952–962.
  * Quan et al. [2018] Quan, T.M., Nguyen-Duc, T., Jeong, W.K., 2018.  Compressed sensing mri reconstruction using a generative adversarial network with a cyclic loss.  IEEE Transactions on Medical Imaging 37, 1488–1497.
  * Sandino et al. [2020] Sandino, C.M., et al., 2020.  Compressed sensing: From research to clinical practice with deep neural networks: Shortening scan times for magnetic resonance imaging.  IEEE Signal Processing Magazine 37, 117–127.
  * Schlemper et al. [2018] Schlemper, J., et al., 2018.  A deep cascade of convolutional neural networks for dynamic mr image reconstruction.  IEEE Transactions on Medical Imaging 37, 491–503.
  * Sodickson and Manning [1997] Sodickson, D.K., Manning, W.J., 1997\.  Simultaneous acquisition of spatial harmonics (smash): fast imaging with radiofrequency coil arrays.  Magnetic resonance in medicine 38, 591–603.
  * Souza and Frayne [2019] Souza, R., Frayne, R., 2019\.  A hybrid frequency-domain/image-domain deep network for magnetic resonance image reconstruction, in: 2019 32nd SIBGRAPI Conference on Graphics, Patterns and Images (SIBGRAPI), IEEE. pp. 257–264.
  * Souza et al. [2019] Souza, R., Lebel, R.M., Frayne, R., 2019.  A hybrid, dual domain, cascade of convolutional neural networks for magnetic resonance image reconstruction, PMLR, London, United Kingdom. pp. 437–446.
  * Sriram et al. [2020] Sriram, A., et al., 2020.  Grappanet: Combining parallel imaging with deep learning for multi-coil mri reconstruction, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 14315–14322.
  * Tavaf et al. [2021] Tavaf, N., Torfi, A., Ugurbil, K., Van de Moortele, P.F., 2021\.  Grappa-gans for parallel mri reconstruction.  arXiv preprint arXiv:2101.03135 .
  * Walsh et al. [2000a] Walsh, D., Gmitro, A., Marcellin, M., 2000a.  Adaptive reconstruction of phased array mr imagery.  Magnetic Resonance in Medicine 43, 682–690.
  * Walsh et al. [2000b] Walsh, D.O., Gmitro, A.F., Marcellin, M.W., 2000b.  Adaptive reconstruction of phased array mr imagery.  Magnetic Resonance in Medicine 43, 682–690.  doi:10.1002/(SICI)1522-2594(200005)43:5<682::AID-MRM10>3.0.CO;2-G.
  * Wang et al. [2009] Wang, K., Du, J., O’Halloran, R., Fain, S., Kecskemeti, S., Wieben, O., Johnson, K.M., Mistretta, C., 2009\.  Ultrashort te spectroscopic imaging (utesi) using complex highly-constrained backprojection with local reconstruction (hypr lr).  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 62, 127–134.
  * Wang et al. [2017] Wang, S., Tan, S., Gao, Y., Liu, Q., Ying, L., Xiao, T., Liu, Y., Liu, X., Zheng, H., Liang, D., 2017\.  Learning joint-sparse codes for calibration-free parallel mr imaging.  IEEE transactions on medical imaging 37, 251–261.
  * Wang et al. [2020a] Wang, S., et al., 2020a.  Deepcomplexmri: Exploiting deep residual network for fast parallel mr imaging with complex convolution.  Magnetic Resonance Imaging 68, 136 – 147.
  * Wang et al. [2004] Wang, Z., et al., 2004.  Image quality assessment: from error visibility to structural similarity.  IEEE transactions on image processing 13, 600–612.
  * Wang et al. [2020b] Wang, Z., et al., 2020b.  Ikwi-net: A cross-domain convolutional neural network for undersampled magnetic resonance image reconstruction.  Magnetic Resonance Imaging 73, 1–10.
  * Weinan [2017] Weinan, E., 2017.  A proposal on machine learning via dynamical systems.  Communications in Mathematics and Statistics 5, 1–11.
  * Yang et al. [2016] Yang, Y., et al., 2016.  Deep admm-net for compressive sensing mri, in: Lee, D.D., Sugiyama, M., Luxburg, U.V., Guyon, I., Garnett, R. (Eds.), Advances in Neural Information Processing Systems 29. Curran Associates, Inc., pp. 10–18.
  * Yang et al. [2020] Yang, Y., et al., 2020.  Admm-csnet: A deep learning approach for image compressive sensing.  IEEE Transactions on Pattern Analysis and Machine Intelligence 42, 521–538.
  * Zbontar et al. [2018] Zbontar, J., et al., 2018.  fastMRI: An open dataset and benchmarks for accelerated MRI.  arXiv:1811.08839 URL: https://arxiv.org/abs/1811.08839.
  * Zhang and Ghanem [2018] Zhang, J., Ghanem, B., 2018\.  Ista-net: Interpretable optimization-inspired deep network for image compressive sensing, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1828–1837.
  * Zhou et al. [2019] Zhou, Z., et al., 2019.  Parallel imaging and convolutional neural network combined fast mr image reconstruction: Applications in low-latency accelerated real-time imaging.  Medical Physics 46, 3399–3413.
  * Zhu et al. [2018] Zhu, B., et al., 2018.  Image reconstruction by domain-transform manifold learning.  Nature 555, 487–492.  doi:10.1038/nature25988.
  * Zhu et al. [2020] Zhu, Y., et al., 2020.  Cross-domain medical image translation by shared latent gaussian mixture model, in: Medical Image Computing and Computer Assisted Intervention – MICCAI 2020, Springer International Publishing, Cham. pp. 379–389.