text:
\section{Introduction}
\label{sec:Introduction}

\vspace{-0.5cm}

Magnetic Resonance Fingerprinting (MRF)~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low,liao20173d,cao2017robust} has emerged as a promising Quantitative Magnetic Resonance Imaging (QMRI) approach, with the capability of providing multiple tissue's intrinsic spin parameters simultaneously, such as the spin-lattice magnetic relaxation time (T1) and the spin-spin magnetic relaxation time (T2). Based on the fact that the response from each tissue with respect to a given pseudo-random pulse sequence is unique, MRF exploits pseudo-randomized acquisition parameters to create unique temporal signal signatures, analogous to a "fingerprint", for different tissues. A dictionary matching operation is then performed to map an inquiry temporal signature to the best matching entry in a precomputed dictionary, leading to multiple tissue parameters directly.

The temporal signatures are generated by varying the acquisition parameters of a pseudo-random excitation pulse sequence, such as repetition time (TR), time of echo (TE), and radio frequency flip angle (FA) over time. The dictionary is composed of a large number of entries that are usually simulated by the Bloch equations given pseudo-random pulse sequences. Each entry represents a unique temporal signature associated with a specific tissue and its quantitative parameters, such as the T1 and T2 relaxation times. Thus, once the best matching (i.e. most correlated) entry is found, it directly leads to multiple quantitative parameters simultaneously via a lookup-table operation.

MRI physics and physiological constraints make the MR scanning procedure time-consuming. To shorten acquisition time, subsampling is commonly performed in k-space (a.k.a conjugate Fourier transform domain) in order to reduce the number of samples and accelerate imaging speed. However, such k-space subsampling results in temporal signatures that are corrupted by aliasing, blurring and noise. This hampers the accuracy associated with estimation of the tissue parameters using a dictionary matching procedure. In order to alleviate the impact of such distortion and corruption, de-aliasing operations are often exploited to restore cleaner signatures before performing signature-to-parameter mapping. Therefore, MRF reconstruction usually involves two operations: signature restoration and parameter restoration.

Inspired by the successful application of sparsity-driven image processing approaches in MRI reconstruction~\cite{lustig2008compressed,weizman2015compressed,weizman2016reference,eldar2015sampling}, several works~\cite{davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low} suggest to incorporate prior knowledge such as sparsity and low-rank to attenuate distortion and corruption, improving the signature restoration performance, during the initial MRF reconstruction stage. This is then followed by a dictionary matching operation, performing mapping from purified temporal signatures to tissue's quantitative parameters. However, such dictionary matching based signature-to-parameter mapping exhibits several drawbacks~\cite{cohen2018mr,hoppe2017deep}. Since the simulated dictionary and lookup-table contain a finite number of elements, they can only cover a limited number of discrete values for each type of tissue parameter. We refer to the difference between a continuous-valued tissue parameter and its closest available discrete value on a lattice as the discretization error. For example, a pair of dictionary and lookup-table that contain 101 elements will lead to a discretization error of maximum 25 ms if they cover the range of 0 ms - 5000 ms with a fixed interval of 50 ms for a specific tissue parameter, e.g. T1. To reduce the discretization error, a huge dictionary that is composed of a large number of entries is needed to represent tissues with fine granularity over the entire value range of target tissue parameters. However, storing a large dictionary becomes prohibitively memory-consuming, as the dictionary size and density often increase exponentially with the number of tissue parameters. Specifically, the number of entries in a dictionary will be $P^s$ for $s$ parameters each containing $P$ values, since every combination of these $s$ parameters determines a specific tissue which is characterized by a specific signature. For example, given T1, T2 relaxation times, i.e. $s=2$, if each of them contains $1000$ values, the dictionary will have $1000^2$ entries. In addition, finding the best matching entry becomes computationally intense for a large dictionary, considerably limiting the inference speed.

In this paper, we propose an alternative approach to dictionary matching based on deep neural networks (a.k.a. deep learning)~\cite{lecun2015deep,goodfellow2016deep}, which we refer to as HYDRA: HYbrid Deep magnetic ResonAnce fingerprinting. The motivation derives from the fact that a well designed and tuned deep neural network is capable of approximating complex functions, leading to state-of-the-art results in a number of tasks such as image classification, super-resolution, speech recognition, and more~\cite{krizhevsky2012imagenet,he2016deep,he2016identity,dong2016image,kim2016deeply,gehring2017convolutional,hinton2012deep}. Recent work~\cite{cohen2018mr,hoppe2017deep} proposed to exploit neural networks to replace the dictionary and lookup-table used in conventional MRF reconstruction approaches. These proposed neural networks suffer from two limitations: First, these approaches are based on neural network models containing only 3-layers, thus suffer from limited capacity of capturing complex mapping functions. Second, these methods focused exclusively on parameter restoration stage (the second stage in MRF reconstruction), but not on signature restoration (the first stage in MRF reconstruction). Therefore, these techniques rely on fully-sampled data instead of typically available sub-sampled k-space data.

Different from Cohen et al.'s fully-connected feed-forward neural network~\cite{cohen2018mr}, and Hoppe et al.'s vanilla convolutional neural network (CNN)~\cite{hoppe2017deep}, the proposed HYDRA involves both a signature restoration and a parameter restoration phase. Signature restoration is implemented using a low-rank based de-aliasing method adapted from Mazor et al.~\cite{mazor2018low} while parameter restoration is implemented using a deep nonlocal residual convolutional neural network developed for this purpose. Our key contributions with respect to prior work are:

\vspace{-0.3cm}

\begin{itemize}
	\setlength\itemsep{0.0em}
	\item
	HYDRA is, to the best of our knowledge, the first deep network approach to combine model-based de-aliasing and learning-based parameter mapping. HYDRA eliminates the requirement for the memory and time-consuming dictionary matching operation, thus significantly improving inference speed without compromising on reconstruction performance.
	\item
	A 1D nonlocal residual convolutional neural network is designed to capture the mappings from temporal MRF signals to tissue parameters. Owing to residual learning and a self-attention mechanism, our network is deeper and more sophisticated than competing network models. This allows to capture complex parameter mappings more effectively, and output continuous parameters to alleviate discretization issues.
	\item
	The designed network is trained on synthesized MRF data simulated with the Bloch equations, but is still applicable to anatomical data. This contributes to eliminating the requirement for a large amount of real MRF data.
	\item
	A low-rank based de-aliasing technique is developed in order to take advantage of temporal similarity for signature restoration.
	\item
	The low-rank based signature restoration is organically combined with the learning-based parameter restoration to achieve fast and accurate MRF reconstruction. Such strategy enables HYDRA to handle both fully-sampled k-space data and more importantly sub-sampled k-space data.
	\item
	A series of numerical experiments are conducted to evaluate the proposed approach on both synthetic and anatomical data. The results demonstrate improved inference speed, accuracy and discretization errors over competing methods~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low,hoppe2017deep,cohen2018mr}.
\end{itemize}

\vspace{-0.3cm}

The rest of the paper is organized as follows. In Section~\ref{sec:Methods}, we formulate the MRF reconstruction problem, introduce related methods, and present our approach, involving the use of a low-rank based signature restoration procedure together with a deep network for parameter restoration. Section~\ref{sec:Experiments} is devoted to experimental results, followed by a discussion in Section~\ref{sec:Discussion} and a conclusion in Section~\ref{sec:Conclusion}.

\section{Materials and Methods}
\label{sec:Methods}

\subsection{The MRF Problem Formulation}
\label{ssec:ProblemForm}

\vspace{-0.2cm}

MRF data is composed of multiple frames sampled in k-space over time. A series of such frames are vectorized and then stacked together along the temporal dimension to construct a measurement matrix $\mathbf{Y} \in \mathcal{C}^{Q \times L}$, where $Q$ is the number of k-space samples in each frame, and $L$ is the number of frames. Due to k-space subsampling, every column vector $\mathbf{Y}_{:,i}$ represents a subsampled Fourier transform of a vectorized image frame $\mathbf{X}_{:,i}$:
\begin{equation}
\label{Eq:Y_X}
\mathbf{Y} = [\mathbf{Y}_{:,1}, \cdots, \mathbf{Y}_{:,L} ]
= [ F_u\{\mathbf{X}_{:,1} \}, \cdots, F_u\{\mathbf{X}_{:,L} \} ] \,,
\end{equation}
where $F_u\{\cdot \}$ denotes a subsampled 2D Fourier transform.

Each column $\mathbf{X}_{:,i}$ represents a MR contrast acquired with RF sequence parameters:
\begin{equation}
\label{Eq:TRE}
\boldsymbol{\Theta}^{TRE}_{:,i} = [TR^i, TE^i, FA^i]^T, \; i \in [1,L]
\end{equation}
where $TR^i$ and $TE^i$ denote the repetition time and echo time, respectively, and $FA^i$ denotes the flip angle of the RF pulse during sampling the $i$-th contrast. Every row $\mathbf{X}_{j,:}$ represents a temporal signature, i.e. temporal signal evolution of a specific tissue at the $j$-th image pixel. The signature depends on the tissue's relaxation times, such as T1 and T2, grouped as a row vector:
\begin{equation}
\label{Eq:T12}
\boldsymbol{\Theta}^{T12}_{j,:} = [T1^j, T2^j], \; j \in [1,N]
\end{equation}
where, $N$ denotes the number of pixels in each image frame.
Note that, $j$ is the spatial index while $i$ is the temporal index throughout. Given RF sequence parameters $\boldsymbol{\Theta}^{TRE}$, and parameters $\boldsymbol{\Theta}^{T12}_{j,:}$ of a specific tissue, its temporal signature $\mathbf{X}_{j,:}$ can be derived as:
\begin{equation}
\label{Eq:X_Bloch}
\mathbf{X}_{j,:} = f(\boldsymbol{\Theta}^{T12}_{j,:}, \boldsymbol{\Theta}^{TRE})
\end{equation}
where $f(\cdot)$ denotes the Bloch equations. This MR contrast matrix $\mathbf{X}$ is associated with the k-space measurements $\mathbf{Y}$ per column by the subsampled Fourier transform, and it is related to tissue parameters $\boldsymbol{\Theta}^{T12}$ per row by the Bloch equations, as illustrated in Fig.~\ref{Fig:XY}.

Given RF sequence parameters $\boldsymbol{\Theta}^{TRE}$ and k-space measurements $\mathbf{Y}$, the goal of MRF reconstruction is to estimate the tissue parameters $\boldsymbol{\Theta}^{T12}$. Typically, the image stack $\mathbf{X}$ is first reconstructed from $\mathbf{Y}$, referred to as signature restoration, and then mapped to tissue parameters $\boldsymbol{\Theta}^{T12}$ via dictionary matching, referred to as parameter restoration~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low}. This process is illustrated in Fig.~\ref{Fig:DictMatching}.

The dictionary is a collection of temporal signatures that are usually simulated by the Bloch equations for various typical tissues, given the pseudo-random RF pulse sequences and tissue parameters. Given an inquiry temporal signature, dictionary matching computes the inner product between the temporal signature with each dictionary entry, selecting the entry in the dictionary exhibiting the highest correlation with the inquiry one as the best matching signature. Once the best entry is found, it directly leads to multiple tissue parameters, such as T1, T2, simultaneously, via searching a lookup-table.

Let $\mathbf{LUT} \in \mathcal{R}^{K \times 2} $ denote a lookup-table composed of $K$ tissues, each containing 2 parameters, i.e., T1 and T2 relaxation times%
\footnote{Note that the off resonance parameter, which appeared in the original MRF paper~\cite{ma2013magnetic}, has been omitted here, since the sequence used in our experiments is derived from the FISP sequence, which is insensitive to off resonance effects~\cite{jiang2015mr,mazor2018low}.}. Let $\mathbf{D} \in \mathcal{C}^{K \times L} $ denote the corresponding dictionary simulated using Bloch equations given the RF sequence parameters $\boldsymbol{\Theta}^{TRE}$, formulated as $\mathbf{D}_{k,:} = f(\mathbf{LUT}_{k,:}, \boldsymbol{\Theta}^{TRE})$. Since each temporal signature $\mathbf{D}_{k,:}$ is linked with the $k$-th tissue's parameters $\mathbf{LUT}_{k,:}$, the choice of a large dictionary size $K$ can in principle provide enough granularity to capture a range of possible tissue values.

In conclusion, existing MRF reconstruction approaches involve two stages: signature restoration and parameter restoration, that can be succinctly written as
\begin{equation}
\label{Eq:OverallFormulation}
\boldsymbol{\Theta}^{T12}_{j,:} = g(h(\mathbf{Y})_{j,:} | \boldsymbol{\Theta}^{TRE} ), \; j \in [1, N],
\end{equation}
where the function $\mathbf{X} = h(\mathbf{Y})$ represents the signature restoration operation such as sparsity or low-rank based de-aliasing and denoising methods, whereas $ g(\mathbf{X}_{j,:} | \boldsymbol{\Theta}^{TRE} )$ denotes the parameter restoration operation, such as dictionary matching based methods~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low}.

Our approach aims to perform signature restoration via low-rank based de-aliasing and parameter restoration via a neural network in order to achieve improved MRF reconstruction performance. We highlight that our method only requires a simulated dictionary during network training. Once the network is trained, the dictionary is not needed anymore. In addition, our approach also eliminates a simulated dictionary for signature restoration, which is a key difference from FLOR~\cite{mazor2016low,mazor2018low} during signature restoration.

\begin{algorithm}[t]
	\small
	\caption{Original MRF method~\cite{ma2013magnetic}}
	\label{Alg:MRF_DictMatching}
	\textbf{Input:}

	A set of subsampled k-space images: $\mathbf{Y}$

	A pre-simulated dictionary: $\mathbf{D}$

	An appropriate lookup-table: $\mathbf{LUT}$

	\textbf{Output:}

	Magnetic parameter maps: $\widehat{T}_1$, $\widehat{T}_2$ %

	\textbf{Step 1. Restore signatures:}

	$$\widehat{\mathbf{X}}_{:,i} = F_u^H\{ \mathbf{Y}_{:,i} \}, \, \forall i$$

	\textbf{Step 2. Restore parameters for every $j$ via dictionary matching:}

	$$
	\widehat{k_j} = \arg \underset{k}{\max}
	\frac{\textbf{Re} \left< \mathbf{D}_{k,:}, \widehat{\mathbf{X}}_{j,:} \right> }{\| \mathbf{D}_{k,:} \|_2^2}
	\; , \; \quad
	\widehat{T}_1^j,  \widehat{T}_2^j = \mathbf{LUT}(\widehat{k_j})
	$$
\end{algorithm}

\subsection{Previous Methods}
\label{ssec:RelatedWork}

\vspace{-0.5cm}

\mypar{Dictionary Matching based MRF approaches}
The original MRF reconstruction algorithm~\cite{ma2013magnetic} is based on dictionary matching, as presented in Algorithm~\ref{Alg:MRF_DictMatching}. It finds the best matching dictionary entry for the acquired temporal signature according to their inner product and then searches the lookup-table to obtain corresponding tissue parameters. Here, $F_u^H\{\cdot\}$ denotes the inverse Fourier transform operating on the zero filled k-space data where zeros are filled at the unknown frequencies and symbol $\textbf{Re} \left< \mathbf{a}, \mathbf{b} \right>$ represents the real part of the inner product of two vectors $\mathbf{a}$ and $\mathbf{b}$.

Exploiting the nature of signals, by using appropriate prior knowledge, can often contribute to improved signal processing performance. In this spirit, later works suggested to incorporate sparsity in MRF reconstruction to further improve performance, inspired by successful applications of sparsity in MRI reconstruction~\cite{lustig2008compressed,weizman2015compressed,weizman2016reference}. Davies et al.~\cite{davies2014compressed} proposed BLoch response recovery via Iterative Projection (BLIP) which exploits sparsity in the dictionary domain. BLIP consists of iterating between two main steps: (a) a gradient step which enforces consistency with the measurements, based on the Projected Landweber Algorithm (PLA) generalized from the iterative hard thresholding method; (b) a projection which matches each row of $\mathbf{X}$ to a single dictionary atom. Instead of exploiting sparsity in the dictionary domain, Wang et al.~\cite{wang2016magnetic} suggested to leverage sparsity in the wavelet domain of each imaging frame, $\mathbf{X}_{:,i}$. They further replaced the Euclidean norm with the Mahalanobis distance for dictionary matching. Considering that adjacent MR image frames along the temporal dimension should exhibit high resemblance, Mazor et al.~\cite{mazor2016low,mazor2018low} proposed a magnetic resonance Fingerprint with LOw-Rank prior for reconstructing the image stack and quantitative parameters, referred to as FLOR, which achieved state-of-the-art performance. The algorithm, described in Algorithm~\ref{Alg:FLOR}, relies on two priors: a low rank prior on the the matrix $\mathbf{X}$, and the fact that the rows of $\mathbf{X}$ lie in the column space of the dictionary $\mathbf{D}$.

\begin{algorithm}[thb]
	\small
	\caption{FLOR~\cite{mazor2018low}}
	\label{Alg:FLOR}
	\textbf{Input:}

	A set of subsampled k-space images: $\mathbf{Y}$;
	A pre-simulated dictionary: $\mathbf{D}$;
	An appropriate lookup-table: $\mathbf{LUT}$;
	Parameters $\mu$ for gradient step and $\lambda$ for regularization

	\textbf{Output:}

	Magnetic parameter maps: $\widehat{T}_1$, $\widehat{T}_2$

	\textbf{Initialization:}

	$\widehat{\mathbf{X}}^0=0$, $\mathbf{P}=\mathbf{D}^\dagger \mathbf{D}$, where $\mathbf{D}^\dagger$ is the pseudo-inverse of $\mathbf{D}$.

	\textbf{Step 1. Restore signatures via iterating until convergence:}

	\begin{itemize}
		\item
		Gradient step for every $i$:
		\begin{equation} \label{Eq:GradientStep}
		\widehat{\mathbf{Z}}_{:,i}^{t+1} = \widehat{\mathbf{X}}_{:,i}^t
		- \mu F_u^H\{F_u\{\widehat{\mathbf{X}}_{:,i}^t\} - \mathbf{Y}_{:,i}\}
		\end{equation}
		where the superscript $^t$ represents the index of iterations.
		\item
		Project onto the dictionary subspace:
		\begin{equation} \label{Eq:ProjectOnDict}
		[\mathbf{U}, \mathbf{S}, \mathbf{V}] = \mathrm{svd}(\widehat{\mathbf{Z}}^{t+1} \mathbf{P})
		\end{equation}
		where $\mathrm{svd}$ denotes the singular-value decomposition operation, and $\mathbf{S} = diag(\{\sigma_j\})$ is a rectangular diagonal matrix with singular values $\{\sigma_j\}$ on its diagonal.
		%
		\item
		Soft-threshold the non-zero singular values with  $\lambda \mu$ and reconstruct signatures $\widehat{\mathbf{X}}^{t+1}$:
		%
		\begin{equation} \label{Eq:Soft-threshold}
		\begin{split}
		\sigma_j' = \max \{ \sigma_j - \lambda \mu, 0 \}
		\; , \; \quad
		\widehat{\mathbf{X}}^{t+1} = \mathbf{U} \mathbf{S}' \mathbf{V}^H
		\end{split}
		\end{equation}
		where $\mathbf{S}' = diag(\{\sigma_j'\})$.
	\end{itemize}

	\textbf{Step 2. Restore parameters for every $j$ via dictionary matching:}

	$$
	\widehat{k_j} = \arg \underset{k}{\max}
	\frac{\textbf{Re} \left< \mathbf{D}_{k,:}, \widehat{\mathbf{X}}_{j,:} \right> }{\| \mathbf{D}_{k,:} \|_2^2}
	\; , \; \quad
	 \widehat{T}_1^j,  \widehat{T}_2^j = \mathbf{LUT}(\widehat{k_j})
	$$
\end{algorithm}

\mypar{Learning-based MRF approaches}
The above techniques all use dictionary matching to perform mapping from temporal signatures to tissue parameters. Therefore, these methods suffer from drawbacks such as discretization error, slow inference speed and memory-consuming storage. In order to alleviate these issues, recent works~\cite{cohen2018mr,hoppe2017deep} propose to exploit neural networks to replace dictionaries and lookup-tables used in conventional MRF reconstruction approaches. Cohen et al. suggest a fully-connected feed-forward neural network (FNN)~\cite{cohen2018mr}. Since the input layer of the FNN is fully connected with the input temporal signature, the number of neurons in the input layer corresponds to the length of the input temporal signature. This makes the network structure less flexible, as a FNN network trained on temporal signatures with a certain length is not applicable to temporal signatures with a different length. In addition, the fully-connected structure results in rapid increase in the number of parameters along with the growth of depth, making the network more susceptible to overfitting. Hoppe et al.~\cite{hoppe2017deep} propose a 3-layer vanilla CNN for parameter restoration. Both \cite{cohen2018mr} and \cite{hoppe2017deep} focus exclusively on learning the signature-to-parameter mapping from a pair of dictionary and lookup-table simulated using the Bloch equations. During the validation, they assume that clean temporal signatures are available as input into the trained networks. However, since temporal signatures obtained from k-space subsampled MRF data are always contaminated by aliasing and noise, their approaches, when applied directly in such k-space subsampling situations, suffer from heavy artifacts introduced during the signature restoration phase, leading to poor performance.

\vspace{-0.2cm}

\subsection{Proposed Methods}
\label{ssec:ProposedMethods}

\vspace{-0.2cm}

The proposed hybrid deep magnetic resonance fingerprinting (HYDRA) approach, summarized in Algorithm \ref{Alg:LowRank_ResCNN}, consists of two stages: signature restoration and parameter restoration, (see also \eqref{Eq:OverallFormulation}). As illustrated in Fig.~\ref{Fig:ResNet-CNN}, a low-rank based de-aliasing method is used to restore signatures, and then a 1D nonlocal residual convolutional neural network is used to map each restored signature to corresponding tissue parameters.

In particular, given $\boldsymbol{\Theta}^{TRE}$ and k-space samples $\mathbf{Y}$, in our proposed approach, the function $\mathbf{X} = h(\mathbf{Y})$ in~\eqref{Eq:OverallFormulation} represents a signature restoration operation using low-rank based de-aliasing techniques without requiring a dictionary. The function $\boldsymbol{\Theta}^{T12}_{j,:} = g(\mathbf{X}_{j,:} | \boldsymbol{\Theta}^{TRE} )$ in~\eqref{Eq:OverallFormulation} represents a parameter restoration operation that exploits a trained neural network to map each restored signature $\mathbf{X}_{j,:}$ to corresponding tissue parameters $\boldsymbol{\Theta}^{T12}_{j,:}$ directly. In the subsequent sections, we provide a detailed description of both stages of our technique.

\vspace{-0.2cm}

\subsubsection{Low-rank based signature restoration}
\label{sssec:Low-rankSignature}

\vspace{-0.2cm}

Since MRF data consists of multiple frames exhibiting temporal similarity, the imaging contrasts matrix $\mathbf{X}$ is typically a low-rank matrix~\cite{mazor2018low}. Therefore, $h(\cdot)$ leverages a low-rank prior for denoising and de-aliasing, formulated as
\begin{equation}
\label{Eq:LowRank_Constraint}
\begin{array}{rl}
h(\mathbf{Y}) =
\arg \underset{\mathbf{X}}{\min}
&
\frac{1}{2} \sum_i \|\mathbf{Y}_{:,i} - F_u\{ \mathbf{X}_{:,i} \} \|_2^2
\\
\text{s.t.}
&
\textrm{rank}(\mathbf{X}) < r
\end{array}
\end{equation}
where the parameter $r$ is the rank of the matrix, a fixed pre-chosen parameter. Since typically $r$ is not known in advance, we consider a relaxed regularized version:

\vspace{-0.5cm}

\begin{equation}
\label{Eq:LowRank_NuclearRegularization}
h(\mathbf{Y}) = \arg \underset{\mathbf{X}}{\min}
\frac{1}{2} \sum\nolimits_i \|\mathbf{Y}_{:,i} - F_u\{ \mathbf{X}_{:,i} \} \|_2^2
+ \lambda \, \|\mathbf{X} \|_*
\end{equation}

\vspace{-0.2cm}

\noindent where $\|\mathbf{X} \|_*$ denotes the nuclear norm~\cite{cai2010singular} of $\mathbf{X}$, defined as the sum of the singular values of $\mathbf{X}$, and $\lambda$ is the Lagrangian multiplier manually selected for balancing data fidelity and the rank. Problem~\eqref{Eq:LowRank_NuclearRegularization} can be solved using the incremental subgradient proximal method~\cite{sra2012optimization}, similar as to FLOR~\cite{mazor2018low}. The procedure for solving \eqref{Eq:OverallFormulation} is shown in Algorithm \ref{Alg:LowRank_ResCNN}.

One of differences from FLOR~\cite{mazor2018low} is the fact that we removed the operation of projecting the temporal signal onto a dictionary. This allows to eliminate the requirement for a simulated dictionary in the signature restoration stage, which also alleviates the memory consumption issue. In addition, the computational complexity is reduced by $ N \cdot L^2$ floating-point operations in each iteration, where $L$ is the dimension of a dictionary element, and $N$ is the number of pixels in each image frame. On the other hand, the gained benefits are at the price of requiring more iterations to converge. Another difference from FLOR~\cite{mazor2018low} is that we exploit a network, instead of dictionary matching, for signature-to-parameter mapping.

\begin{algorithm}[t]
	\small
	\caption{Proposed MRF reconstruction approach: HYDRA}
	\label{Alg:LowRank_ResCNN}

	\textbf{Input:}

	A set of subsampled k-space images: $\mathbf{Y}$

	The trained network: $g$

	Parameters $\mu$ for gradient step and $\lambda$ for regularization

	\textbf{Output:}
	Magnetic parameter maps $\widehat{T}_1$, $\widehat{T}_2$ %

	\textbf{Initialization:}  $\widehat{\mathbf{X}}^0 = 0$

	\textbf{Step 1. Restore signatures via iterating until convergence:}

	\begin{itemize}
		\item
		Gradient step for every $i$, the same as \eqref{Eq:GradientStep}.
		\item
		Perform SVD:
		\begin{equation*}
		[\mathbf{U}, \mathbf{S}, \mathbf{V}] = \mathrm{svd}(\widehat{\mathbf{Z}}^{t+1})
		\end{equation*}
		\item
		Soft-threshold the non-zero singular values $\{\sigma_j\}$ of $\mathbf{S}$ with parameter $\lambda \mu$ and reconstruct signatures $\widehat{\mathbf{X}}^{t+1}$, the same as \eqref{Eq:Soft-threshold}.
	\end{itemize}

	\textbf{Step 2. Restore parameters for every $j$ via the trained network:}

	$$ \widehat{T}_1^j, \widehat{T}_2^j = g(\widehat{\mathbf{X}}_{j,:})$$

\end{algorithm}

\subsubsection{Learning-based parameter restoration}
\label{sssec:ParamRestore}

\vspace{-0.2cm}

Once the imaging contrasts matrix $\mathbf{X}$ is recovered from the k-space samples $\mathbf{Y}$, each temporal signature $\mathbf{X}_{j,:}$ is input into the trained network for parameter restoration, formulated as:
\begin{equation}
\label{Eq:mapping}
\boldsymbol{\Theta}^{T12}_{j,:} = g(\mathbf{X}_{j,:}  | \boldsymbol{\Theta}^{TRE}) , \; j \in [1, N]
\end{equation}
where $g(\cdot)$ denotes the trained network, $\boldsymbol{\Theta}^{TRE}$ denotes the fixed RF sequence parameters. We next describe the network structure, training and testing procedures.

\mypar{Network structure}

The proposed network has a 1D nonlocal residual CNN architecture with short-cuts for residual learning and nonlocal operations for achieving a self-attention mechanism. As illustrated in Fig.~\ref{Fig:ResNet-CNN}, it starts with two 1D convolutional layers before connecting with 4 residual / non-local operation blocks, and finally ends with a global-average-pooling layer followed by a fully-connected layer. Every residual block is followed by a non-local operation block. Four such blocks are interspersed with each other.

Each residual block contains a max-pooling layer with stride 2, two convolution layers and a shortcut that enforces the network to learn the residual content. The filter size is set to be equal to 21 throughout convolutional layers. The number of channels, a.k.a feature maps, in the first two convolutional layers is set to 16 and then is doubled in subsequent four residual blocks until 128 in the final residual block. The size of feature maps in the next block halves in contrast with the previous one due to max-pooling. In this way, we gradually reduce temporal resolution while extract more features to increase content information.

Inspired by the self-attention scheme and nonlocal neural networks~\cite{vaswani2017attention,wang2018non}, non-local operations are incorporated into the designed network to achieve the attention mechanism, in order to capture long-range dependencies with fewer layers. In contrast to the progressive behavior of convolutional operations that process one local neighborhood at a time, the non-local operations compute the response at a position as a weighted sum of the features at all positions in the feature maps. Formally, the nonlocal operation is formulated as~\cite{wang2018non}:
\begin{equation}
	\mathbf{y}_i = \frac{1}{C(\mathbf{x})} \sum_{\forall j} f(\mathbf{x}_i, \mathbf{x}_j) g(\mathbf{x}_j) \,.
	\label{Eq:NonlocalOperation}
\end{equation}
Here, $i$ is the index of an output position and $j$ is the index that enumerates all possible temporal positions, $\mathbf{x}$ is the input temporal signal or its features and $\mathbf{y}$ is the output signal of the same size as $\mathbf{x}$. A pairwise function $f$ computes a scalar between $i$ and all $j$ to represent the affinity relationship of these two positions. The unary function $g$ computes a representation of the input signal at position $j$. The response is normalized by a factor $C(\mathbf{x})$. There exists a few instantiations for function $f$ and $g$. For simplicity, the unary function $g$ is chosen as a linear embedding: $g(\mathbf{x}_j) = \mathbf{W}_g (\mathbf{x}_j)$, where $\mathbf{W}_g$ is a weight matrix to be learned. Regarding the affinity matrix $f$, we adopt the embedded Gaussian to compute similarity in an embedding space, which is formulated as:
$f(\mathbf{x}_i, \mathbf{x}_j) = e^{{\theta(\mathbf{x}_i)}^\top \phi(\mathbf{x}_j)}.$ Here, $\theta(\mathbf{x}_i) = \mathbf{W}_\theta \mathbf{x}_i$ and $\phi(\mathbf{x}_j) = \mathbf{W}_\phi \mathbf{x}_j$ are two learned embeddings. The normalization factor is set as $C(\mathbf{x}) = \sum_{\forall j} f(\mathbf{x}_i, \mathbf{x}_j)$. For a given $i$, $\frac{1}{C(\mathbf{x})} f(\mathbf{x}_i, \mathbf{x}_j)$ becomes the softmax computation along the dimension $j$, which leads to $\mathbf{y} = \textit{softmax}(\mathbf{x}^\top \mathbf{W}_\theta^\top \mathbf{W}_\phi \mathbf{x}) g(\mathbf{x})$, which is the self-attention form~\cite{vaswani2017attention}.

The non-local behavior in \eqref{Eq:NonlocalOperation} is due to the fact that all positions $(j)$ are considered in the operation. As a comparison, a convolutional operation sums up the weighted input in a local neighborhood~\cite{wang2018non}. It implies that the non-local operation directly captures long-range dependencies in the temporal dimension via computing interactions between any two points, regardless of their positional distance. In this way, the network is able to extract global features and take advantage of the full receptive field in each layer.

The global-average-pooling layer is used to average each feature map in order to integrate information in every channel for improved robustness to corrupted input data. This global-average-pooling layer also reduces the number of parameters significantly, thus lessening the computation cost as well as preventing over-fitting. The last fully-connected layer outputs estimated parameters -- T1 and T2 relaxation times. The designed network contains around 0.27 million parameters. The weights are initialized using He-normal-distribution~\cite{he2015delving}. The max-norm kernel constraint~\cite{srivastava2014dropout} is exploited to regularize the weight matrix directly in order to prevent over-fitting. The designed network can also be adapted for various MRF sequences, such as the original MRF sequence -- inversion-recovery balanced steady state free-precession (IR-bSSFP) sequence, that depends also on the intrinsic df parameter. It is possible to adjust the number of outputs to adapt to more parameters, such as proton density, B0.

To summarize, our network is motivated and inspired by recent successful applications of convolutional neural networks and variants. Convolutional neural networks have been proved to be a powerful model to capture useful features from signals and images. By introducing convolution, local receptive field and weight sharing design, a CNN is capable of taking advantage of local spatial coherence and translation invariance characteristics in the input signal, thus become especially well suited to extract relevant information at a low computational cost~\cite{krizhevsky2012imagenet,he2016deep,he2016identity,dong2016image,kim2016deeply,gehring2017convolutional}. On the other hand, the residual network architecture~\cite{he2016deep,he2016identity} provides an effective way to design and train a deeper model, since it alleviates the gradient vanishing or exploding problems by propagating gradients throughout the model via short-cuts, a.k.a skip connections. By leveraging non-local operation based attention mechanism, neural networks are endowed with capability of extracting global features and capturing long-range dependencies.

\vspace{-0.2cm}

\mypar{Network training}

The designed network is trained on a synthesized dictionary $\mathbf{D}$ and corresponding lookup-table $\mathbf{LUT}$ to learn the signature-to-parameter mappings $\mathbf{LUT}_{k,:} = g(\mathbf{D}_{k,:} | \boldsymbol{\Theta}^{TRE})$.

The training dataset is synthesized as follows. First, we determine the range of tissue parameters. For example, one may set T1 relaxation times to cover a range of [1, 5000] ms and T2 relaxation times to cover a range of [1, 2000] ms with an increment of 10 ms for both. Thus, the T1 and T2 values constitute a grid with dimension 500 $\times$ 200, in which each point represents a specific combination of T1 and T2 values, and hence characterizes a specific tissue. Points corresponding to T1 < T2 have been excluded as such combinations have no physical meaning. All the valid points are stacked together to generate a lookup-table. For instance, the above setting for T1 and T2 leads to a lookup-table of dimension $80100 \times 2$. The RF pulse sequences used in our work are fast imaging with steady state precession (FISP) pulse sequences with parameters that have been used in previous publications in the field of MRF~\cite{mazor2018low,jiang2015mr,cao2017robust}. Given the lookup-table and RF pulse sequences, dictionary entries can be synthesized by solving the Bloch equations using the extended phase graph formalism~\cite{hennig1991echoes,weigel2015extended}.

When the training dataset is ready, the dictionary entries are used as input signals and corresponding lookup-table entries serve as the groundtruth. All the dictionary entries are input into the designed network batch by batch which outputs estimated parameters. The root mean square errors (RMSE) of the outputs are calculated with respect to corresponding groundtruth. The resulting RMSE loss is then backpropagated from the output layer to the first layer to update the weights and bias by using Adam~\cite{kingma2014adam} as the optimization algorithm. More training details are provided in the subsequent experiment section. Once the training procedure is completed, given an inquiry signal evolution $\mathbf{X}_{j,:}$, it is able to map such a time sequence directly to corresponding tissue parameters, as formulated in \eqref{Eq:mapping}, implying that no dictionary or lookup-table are required during the inference. Since we only need to store the trained network which is a compact model, it consumes less memory than storing the dictionary and lookup-table.

We emphasize that even though the network is trained on a grid of tissue values, it is expected to capture the mapping function from temporal signatures to tissue parameters. Thus the trained network is capable of outputting tissue values not existing in the grid of training values. Detailed results can be found in Fig.~\ref{Fig:ResNet-CNN_1Dtest} and Table~\ref{Tab:TestingSynthetic}. This feature is favorable, as it implies that well designed and trained networks have an ability to overcome discretization issues. The overall procedures for solving \eqref{Eq:LowRank_NuclearRegularization} and \eqref{Eq:mapping} are shown in Algorithm \ref{Alg:LowRank_ResCNN}.

\section{Experimental Results}
\label{sec:Experiments}

\vspace{-0.3cm}

In this section, we conduct a series of experiments to evaluate our approach, comparing it with other state-of-the-art MRF methods~\cite{ma2013magnetic,davies2014compressed,mazor2018low,hoppe2017deep,cohen2018mr}.

The experiments are categorized into a few types: training, testing on synthetic data, testing on anatomical data using variable density Gaussian sampling patterns and spiral sampling patterns at different sampling ratios and number of time frames, as described in Table~\ref{Tab:ExperimentsTypes}. For the network training, synthesized temporal signatures, i.e. simulated dictionary entries of $\mathbf{D}$ shown as Fig.~\ref{Fig:SynSignatures}, are used as input signals and corresponding parameter values in the lookup-table $\mathbf{LUT}$ serve as the groundtruth. The proposed network is trained to capture the signature-to-parameter mappings. For testing on synthetic data, synthesized temporal signatures in $\mathbf{X}$ are used as input signals and corresponding parameter values in $\boldsymbol{\Theta}^{T12}$ serve as the groundtruth. The aim is to test the parameter restoration performance only. For testing on anatomical data, the k-space measurements $\mathbf{Y}$ which are derived from the Fourier transform (for Gaussian patterns) or non-uniform FFT (for spiral trajectories)~\cite{fessler2003nonuniform} of $\mathbf{X}$, are used as input and corresponding parameter values in $\boldsymbol{\Theta}^{T12}$ serve as reference. When there is no k-space subsampling, the aim is to test the parameter restoration performance only. When there exists k-space subsampling, the aim is to test the overall performance, including both signature restoration and parameter restoration. More detailed descriptions are provided in each subsection.

\vspace{-0.2cm}

\subsection{Training}
\label{ssec:Exp_Training}

\vspace{-0.2cm}

As mentioned in Section~\ref{sssec:ParamRestore}, the designed network is trained on a pair of synthesized dictionary $\mathbf{D}$ and lookup-table $\mathbf{LUT}$, simulated using Bloch equations and FISP pulse sequences~\cite{jiang2015mr,mazor2018low}.

The FISP pulse sequence used in our experiments was designed with parameters $\boldsymbol{\Theta}^{TRE}_{:,i} = [TR^i, TE^i, FA^i]^T, i \in [1,L]$ that have been used in previous publications in the field of MRF~\cite{mazor2018low,jiang2015mr,cao2017robust}. The echo time $TE^i$ was constant of 2ms. The repetition time $TR^i$ was randomly varied in the range of 11.5 - 14.5 ms with a Perlin noise pattern. All the flip angles $FA^i, i \in [1,L]$ constituted a sinusoidal variation in the range of 0 - 70 degrees to ensure smoothly varying transient state of the magnetization, as shown in Figure~\ref{Fig:FA_TR}.

For the range of tissue parameters, T1 relaxation times are set to cover a range of [1, 5000] ms and T2 relaxation times to cover a range of [1, 2000] ms with an increment of 10 ms for both. Such parameter ranges cover the relaxation time values that can be commonly found in a brain scan~\cite{vymazal1999t1}. All the valid combinations of T1 and T2 values are stacked together, generating a lookup-table $\mathbf{LUT}$ of dimension $K \times 2$ where $K=80100$. Given the lookup-table and RF pulse sequences, dictionary entries are synthesized by solving the Bloch equations using the extended phase graph formalism, leading to a dictionary of dimension $K \times L$ where $L = 200$ or $1000$ is the number of time frames.

When the training dataset is ready, the dictionary entries are used as input signals and corresponding lookup-table entries serve as the groundtruth to train the designed network, as mentioned in Section~\ref{sssec:ParamRestore}. The model was trained for 50 epochs. It takes around 30 seconds for running one epoch on average, thus around 25 minutes for completing 50 epochs, on a NVIDIA GeForce GTX 1080 Ti GPU. In each training epoch, 20\% of the training samples are separated aside for validation dataset. The learning rate decays from 1e-2 to 1e-6 every 10 epochs. Each batch was experimentally set to contain 256 time-sequences in order to balance the convergence rate and weights updating rate well. For comparison purposes, we also implemented Hoppe et al.'s CNN referring to~\cite{hoppe2017deep}, and Cohen et al.'s FNN referring to~\cite{cohen2018mr} with the same structure and parameters as specified in their papers. Then we use the same GPU and training dataset to train their networks with specified learning rate and number of epochs until convergence.

We adopt a few widely used metrics, such as root mean square error (RMSE), signal-to-noise ratio (SNR) and peak signal-to-noise ratio (PSNR) to evaluate the image quality quantitatively. The definitions of RMSE, SNR and PSNR are given as follows:
\begin{align}
\textrm{RMSE} &= \sqrt{
	\frac{
		\| \mathbf{X} - \widehat{\mathbf{X}} \|_F^2
	}
	{ N } } \,,
\label{Eq:Def_RMSE}
\\
\textrm{SNR} &= 20 \log_{10} \frac{\| \mathbf{X} \|_F^2}{\textrm{RMSE}} \,,
\label{Eq:Def_SNR}
\\
\textrm{PSNR} &= 20 \log_{10} \frac{\textrm{PeakVal}}{\textrm{RMSE}} \,,
\label{Eq:Def_PSNR}
\end{align}
where matrices $\mathbf{X}$ and $\widehat{\mathbf{X}}$ denote the ground truth signal and its reconstructed version, respectively, $N$ denotes the total number of elements in the signal and $\| \cdot \|_F$ denotes the Frobenius norm. PeakVal stands for the pixel peak value in an image, e.g., 1 for a normalized signal.

\subsection{Testing on synthetic dataset}
\label{ssec:Testingonsynthetic}

\vspace{-0.2cm}

In this subsection, we evaluate the performance of HYDRA on a synthetic testing dataset. The procedures of constructing a synthetic testing dataset is similar to the construction of the training dataset: 500 different T1 values are randomly selected from 1 - 5000 ms, while 200 different T2 values are randomly selected from 1 - 2000 ms, using random permutation based on uniformly distributed pseudorandom numbers. All the valid combinations from the selected T1 and T2 values are stacked together, generating a parameter matrix $\boldsymbol{\Theta}^{T12}$ of dimension $80000 \times 2$ with $N=80000$. The RF pulse sequences are the same as in the training stage. Given the parameter matrix and RF pulse sequences, input signal signatures are synthesized by solving the Bloch equations using the extended phase graph formalism, leading to a signature matrix $\mathbf{X}$ of dimension $N \times L = 80000 \times 200$, with each row representing a temporal signature corresponding to a specific combination of T1 and T2 values. The signature matrix $\mathbf{X}$ and parameter matrix $\boldsymbol{\Theta}^{T12}$ constitute the synthetic testing dataset, with $\mathbf{X}$ as input and $\boldsymbol{\Theta}^{T12}$ as the groundtruth.

We input the synthetic testing signatures $\mathbf{X}$ into Hoppe et al.'s CNN~\cite{hoppe2017deep}, Cohen et al.'s FNN~\cite{cohen2018mr}, and the network of HYDRA to compare the outputs with groundtruth T1 and T2 values in $\boldsymbol{\Theta}^{T12}$. We also compare with dictionary matching methods~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low} which exploit the same dictionary $\mathbf{D}$ and lookup-table $\mathbf{LUT}$ to find the best matching entry for each signature in $\mathbf{X}$ and then estimate parameter values by searching the lookup-table. As shown in Table~\ref{Tab:TestingSynthetic}, Table~\ref{Tab:TestingSynthetic_Contin}, Fig.~\ref{Fig:1DtestReal_Continuous} and Fig.~\ref{Fig:ResNet-CNN_1Dtest}, the estimated parameter values using the proposed network obtained outstanding agreement with the groundtruth, yielding higher PSNR, SNR and smaller RMSE than the dictionary matching method~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low}, as well as competing networks~\cite{hoppe2017deep,cohen2018mr}.

In particular, to illustrate in detail how well neural networks tackle the discretization issue inherent to dictionary matching, we show the testing performance on continuous-valued T1, T2 parameters which have small intervals, e.g. 0.5ms, that is 20 times smaller than the training grid intervals 10ms, between neighboring values in Table~\ref{Tab:TestingSynthetic_Contin}. Since these values and their corresponding MRF signatures do not exist in the training dictionary and lookup-table%
\footnote{
	As mentioned in the experiment setting in Section~\ref{ssec:Exp_Training}, in the training dataset, T1 relaxation times are set to cover a range of [1, 5000] ms and T2 relaxation times to cover a range of [1, 2000] ms with an increment of 10 ms for both, that is, T1 values = \{1, 11, 21, $\cdots$, 4991\}, and T2 values = \{1, 11, 21, $\cdots$, 1991\}.
},
the dictionary matching methods report a T1 and T2 value -- the closest discretized value present in the dictionary -- that can be quite distinct from the groundtruth. In contrast, the various neural network approaches can potentially learn an underlying mapping from the temporal signatures to the respective T1 and T2 values, leading to estimates that are much closer to the groundtruth. Interestingly, our approach outperforms previous networks~\cite{hoppe2017deep,cohen2018mr} as shown in Table~\ref{Tab:TestingSynthetic_Contin} and Fig.~\ref{Fig:1DtestReal_Continuous}. Evidently, neural networks demonstrate much better robustness to discretization issues, leading to improved parameter restoration in comparison to dictionary based methods.

Another impressive advantage of HYDRA is the fast inference speed. HYDRA takes only 8.2 s to complete the mapping operation for eighty thousand temporal signatures, that is, 53$\times$ faster than dictionary matching. Furthermore, the inference speed of HYDRA is subject to the network topology. That is, once the network structure is fixed, the complexity is fixed. In contrast, the complexity of dictionary matching is limited by the dictionary density. This implies that our advantage will be more prominent in comparison with competing techniques using a dictionary with higher density.

\subsection{Testing on anatomical dataset}

\vspace{-0.3cm}

In this subsection, we evaluate our approach on an anatomical testing dataset. We construct the dataset from brain scans that were acquired with GE Signa 3T HDXT scanner from a healthy subject.\footnote{The experiment procedures involving human subjects described in this paper were approved by the Institutional Review Board of Tel-Aviv Sourasky Medical Center, Israel.} Since there are no groundtruth parameter values for the T1 and T2 parameter maps, we obtain gold standard data by acquiring Fast Imaging Employing Steady-state Acquisition (FIESTA) and Spoiled Gradient Recalled Acquisition in Steady State (SPGR) images, at 4 different flip angles (3$^\circ$,5$^\circ$,12$^\circ$ and 20$^\circ$), and implementing corrections~\cite{liberman2014t1} followed by DESPOT1 and DESPOT2~\cite{deoni2005high} algorithms. The constructed gold standard T1, T2 parameter maps have a dimension of $128 \times 128$ for each map, accordingly leading to a parameter matrix $\boldsymbol{\Theta}^{T12}$ of size $16384 \times 2$ by stacking vectorized T1, T2 maps together. Based on the parameter matrix $\boldsymbol{\Theta}^{T12}$ and pre-defined RF pulse sequences, we generate temporal signatures using Bloch equations, the same mechanism as generating the synthetic testing dataset, leading to a signature matrix $\mathbf{X}$ of dimension $N \times L = 16384 \times 200$. The signature matrix $\mathbf{X}$ and parameter matrix $\boldsymbol{\Theta}^{T12}$ constitute the anatomical testing dataset, with $\mathbf{X}$ as input and $\boldsymbol{\Theta}^{T12}$ as the gold standard reference.

Note that, since the gold standard T1, T2 maps exhibit spatial structures in the image domain, the resulting signature matrix $\mathbf{X}$ can be regarded as a stack of $L=200$ vectorized image frames, where each frame exhibits specific spatial structures. Therefore, it makes sense to perform Fourier transform and k-space subsampling for each column of $\mathbf{X}$ to get k-space measurements $\mathbf{Y}$. This is the key difference between the anatomical dataset and the synthetic dataset.

We first explore the case with full k-space sampling in order to evaluate the parameter restoration performance of HYDRA. Then, we consider situations with k-space subsampling in order to evaluate both the signature restoration and the parameter restoration performance of HYDRA.

\vspace{-0.2cm}

\subsubsection{Full k-space sampling}

\vspace{-0.2cm}

In the first case, the fully-sampled k-space measurements $\mathbf{Y}$, derived from the Fourier transform of $\mathbf{X}$, are used as input to obtain the estimated $\boldsymbol{\Theta}^{T12}$. This is equivalent to inputting $\mathbf{X}$ into the network of HYDRA, or performing dictionary matching based on $\mathbf{X}$ directly, since the inverse Fourier transform of the fully-sampled measurements $\mathbf{Y}$ is exactly the same as $\mathbf{X}$. The aim is to test the parameter restoration performance only. In the experiment, corresponding parameter values in $\boldsymbol{\Theta}^{T12}$ serve as the gold standard reference. For comparison, dictionary matching methods~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low} exploit the same dictionary $\mathbf{D}$ and lookup-table $\mathbf{LUT}$ as in our training stage to find the best matching entry and estimate parameter values for each signature in $\mathbf{X}$.

Visual and quantitative results are shown in Fig.~\ref{Fig:ResNet-CNN_1DtestReal}, Fig.~\ref{Fig:ResNet-CNN_test} and Table~\ref{Tab:TestingPhantom}. It can be seen that our basic version of HYDRA outperforms dictionary matching~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low}, yielding better visual and quantitative performance, e.g., 7.9 dB SNR gains for T2 mapping. The RMSE of T2 mapping is also reduced to 2.498 from 6.252, accordingly. Our nonlocal version of HYDRA achieves even better performance, leading to 10 dB SNR gains with RMSE as small as 1.86. This is owing to the advantage that the trained network is a powerful function approximator, which is able to output well-estimated parameter values based on learnt mapping functions, even though these values do not exist in the training dictionary and lookup-table. In contrast, dictionary matching only matches signatures to discrete parameters existing in the training dataset. In other words, if there are no exact matching dictionary element and parameter values for an inquiry MRF signature, it will find adjacent values as approximations, thus introducing discretization error. On the other hand, the advantage of HYDRA over dictionary matching on T1 mapping is not as significant as on T2 mapping quantitatively. But the visual improvements are evident. A similar trend is observed when comparing our network with competing networks such as Hoppe et al.'s CNN~\cite{hoppe2017deep} and Cohen et al.'s FNN~\cite{cohen2018mr}. In addition, HYDRA takes around 2 s to accomplish the mapping for 16384 signatures, 40$\times$ faster than dictionary matching~\cite{ma2013magnetic,jiang2015mr,davies2014compressed,wang2016magnetic,mazor2016low,mazor2018low}.

\vspace{-0.2cm}

\subsubsection{k-space subsampling using Gaussian patterns}

\vspace{-0.2cm}

In k-space subsampling situations, the developed low-rank based de-aliasing method is applied to restore the signature matrix $\mathbf{X}$ from the measurements matrix $\mathbf{Y}$. Then, the reconstructed $\mathbf{X}$ is used as input into the network for parameter mapping to obtain the corresponding tissue parameter values. In the experiments, the sub-sampling factor $\beta$ is set to be 70\% and 15\%. For $\beta$ = 15\%, 15\% k-space data is acquired by a series of 2D random Gaussian sampling patterns, shown in Fig.~\ref{Fig:SamplePatterns}, leading to a k-space measurement matrix $\mathbf{Y}$ of size $Q \times L = 2458 \times 200$. Similarly, $\beta$ = 70\% gives rise to a k-space measurement matrix $\mathbf{Y}$ of size $Q \times L = 11469 \times 200$. A larger $\lambda $ enforces lower rank for the restored signature matrix $\mathbf{X}$ to strengthen the de-aliasing effect, while a smaller $\lambda $ encourages $\mathbf{X}$ to have a subsampled Fourier transform that approximates the k-space measurements matrix $\mathbf{Y}$ better. Therefore, we tried a range of values from 1 to 20 for $\lambda $ and experimentally select the best one $\lambda = 5$. Since the low-rank based signature restoration involves gradient descent steps, a larger step size $\mu$ accelerates gradient descent speed, but tends to result in oscillation or even divergence, while a smaller $\mu$ leads to a slower convergence. We experimentally find that $\mu = 1$ gives a good balance. The same k-space measurements $\mathbf{Y}$ are also used by dictionary-matching based methods~\cite{ma2013magnetic,davies2014compressed,mazor2018low} for comparison, and the same signature restoration approach is used to convert $\mathbf{Y}$ onto $\mathbf{X}$ for learning based methods~\cite{hoppe2017deep,cohen2018mr}. The aim is to evaluate the overall performance on both signature restoration and parameter restoration.

Quantitative performance is shown in Table~\ref{Tab:TestingPhantomSubsamplingNonlocal}. Note that the advantage of learning-based methods over dictionary matching degrades when the subsampling factor increases. This is due to the fact that the restored signatures from highly subsampled k-space data exhibit deviations and distortions, thus leading to poorer input for the trained networks. In spite of this, the proposed approach outperforms the dictionary matching based methods~\cite{ma2013magnetic,davies2014compressed} with significant gains, and also yields better or comparable performance as the state-of-the-art methods FLOR~\cite{mazor2018low}, CNN~\cite{hoppe2017deep} and FNN~\cite{cohen2018mr}. In addition, it takes around 23s for low-rank based signature restoration and less than 3s for network based parameter restoration. Thus, the total time cost is around 26s, almost 4.8$\times$ faster than FLOR~\cite{mazor2018low}. Furthermore, the speed of our method is 60$\times$ faster than FLOR~\cite{mazor2018low} for parameter restoration.

We compared the performance with/without nonlocal operations in our developed network. The results in Table~\ref{Tab:TestingPhantom} and \ref{Tab:TestingPhantomSubsamplingNonlocal} show that the proposed network with nonlocal operations based self-attention scheme outperforms the basic counterpart. In particular, the nonlocal version achieves 6 dB gains in terms of SNR over the basic version for T2 mapping. Such significant improvement demonstrates the benefits of capturing long-range dependencies and global features using the nonlocal operation based attention scheme.

We also investigated the performance with respect to the number of time frames. In particular, we increased $L$ from 200 to 1000 and kept other experiment settings the same as before. The quantitative results are shown in Table~\ref{Tab:GaussSubsamp_L1000}. It is noticed that given more time frames, all the methods show better performance. Moreover, the performance of learning-based methods, including~CNN~\cite{hoppe2017deep}, FNN~\cite{cohen2018mr} and HYDRA, improve more than model-based techniques~\cite{ma2013magnetic,davies2014compressed,mazor2018low}. In particular, our approach outperforms competing algorithms quantitatively in terms of PSRN, SNR, and RMSE, as well as demonstrates visual advantage, as shown in Figure~\ref{Fig:KSpaceGaussRatio_15} and Figure~\ref{Fig:GaussSubsample0_15_L1000}.

\subsubsection{k-space subsampling using spiral trajectories}

We carried out additional experiments with widely used non-Cartesian sampling patterns -- variable density spiral trajectories~\cite{lee2003fast,mazor2018low}. A set of spiral trajectories used in the experiments are shown in Figure~\ref{Fig:SpiralSamplePatterns}. They have FOV of 24 and rotation angle difference of 7.5 degrees between any two adjacent spirals to spread out the alias artifacts. Given such spiral trajectories, data were subsampled to acquire 1488 k-space samples in each time frame, leading to a subsampling ratio of 9\% which is defined by the number of acquired samples in the k-space domain divided by the number of pixels in a frame. This setting closely matches the original MRF paper~\cite{ma2013magnetic} where each single spiral trajectory samples 1450 k-space points (leading to a subsampling ratio around 9\%) and any two adjacent spiral trajectories have a rotation angle of 7.5 degrees.

In the case of spiral subsampling, during the signature restoration, SParse Uniform ReSampling (SPURS) algorithm~\cite{kiperwas2017spurs} was exploited to implement nonuniform Fourier transform between k-space domain and image domain, as SPURS has proved to achieve smaller approximation errors while maintaining low computational cost comparing with other resampling methods, such as nonuniform-FFT algorithm~\cite{fessler2003nonuniform} and regularized Block Uniform ReSampling (rBURS)~\cite{rosenfeld2002new}. In the experiments, 1000 density variable spiral trajectories were used for k-space subsampling, leading to 1000 time frames. The quantitative and qualitative reconstruction results demonstrate that HYDRA outperforms competing methods with smaller estimation errors, as shown in Table~\ref{Tab:SpiralSubsamp0_09_L1000}, Figure~\ref{Fig:SpiralSubsample0_09_L1000} and Figure~\ref{Fig:KSpaceSpiralRatio_09}.

\vspace{-0.3cm}

\section{Discussion}
\label{sec:Discussion}

\vspace{-0.3cm}

\subsection{Relation to previous works}

\vspace{-0.3cm}

Our low-rank based signature restoration method is adapted from FLOR~\cite{mazor2018low} by removing the operation of projecting the temporal signal onto a dictionary. Thus, the signature restoration does not require a simulated dictionary, and saves computational cost. Although recent works~\cite{cohen2018mr,hoppe2017deep} exploit neural networks to perform parameter mapping, replacing dictionaries and lookup-tables used in conventional MRF reconstruction approaches, our technique is different from these methods~\cite{cohen2018mr,hoppe2017deep}. We design a deep nonlocal residual CNN for capturing signature-to-parameter mapping which is organically combined with low-rank based de-aliasing techniques for signature restoration. In this way, our algorithm can bypass some of the issues associated with other techniques:
(1) The input dimension issue. The proposed approach can ingest temporal signatures with different lengths without the need to change the structure of the network. This is due to the fact that we rely on convolutional neural networks (CNNs) rather than fully-connected neural networks (FNNs) such as the model used in \cite{cohen2018mr}.
(2) The k-space subsampling issue. The proposed approach involves a hybrid of a neural network with a low-rank based de-aliasing approach. Thus it is able to deal with correlations both over time and space via exploiting low-rank regularization and convolution operation. This enables our work to handle k-space subsampling situations.
(3) The complex mappings issue. By exploiting a residual network structure, our method can be successfully extended to deeper levels and thus obtain a better capacity to learn complex signature-to-parameter mapping functions.
(4) Distortion and corruption issue. Due to the subsampling in k-space, the restored temporal signatures suffer from local distortion and corruption. Such deviation may lead to performance degeneration in the second stage. By incorporating non-local operations in the network design, our method is able to capture global features and find most relevant components for inference, thereby reducing interference of local distortion and corruption.

\vspace{-0.3cm}

\subsection{Computational complexity}

\vspace{-0.3cm}

HYDRA involves two main stages: the low-rank based signature restoration stage and the network based parameter restoration stage. Even though the time cost for parameter restoration is longer than previous methods~\cite{cohen2018mr,hoppe2017deep}, the time cost in the this stage is only a small fraction of the total time consumption, as the computational complexity is dominated by the signature restoration stage. In other words, the computational burden of HYDRA lies in the SVD calculation in the first stage. Hence, fast SVD methods can be employed to dramatically improve the efficiency of signature restoration.

\subsection{Model storage requirements}

Regarding the storage requirement (in double precision), HYDRA needs only 2.1 megabytes to store the network with 0.5 million parameters, while it requires 108 megabytes to store a simulated dictionary of size $80100 \times 200$ and 551 megabytes for size $80100 \times 1000$. Note that the dictionary volume will grow exponentially with the number of parameters, but the space required for storing a network is not strictly limited by the dictionary density once the topology of the network is fixed, thus significantly alleviating the storage burden inherent to the exponential growth of multi-dimensional dictionaries.

\subsection{Impact of providing continuous T1/T2 values}

Providing continuous T1/T2 values is an advantage of neural network based parameter mapping over dictionary matching. This property may find promising applications in some practical scenarios, for exampling, monitoring sensitive changing of pathology condition over time, such as multiple sclerosis~\cite{manfredonia2007normal,papadopoulos2010ti}, stroke~\cite{bernarding2000histogram}, and treatment responses~\cite{mcsheehy2010quantified,weidensteiner2014tumour}, where the differences in T1 and T2 values between healthy and diseased tissues or between disease stages could be very small~\cite{jiang2017repeatability}. On the other hand, to fulfil this potential of network based MRF techniques, prerequisites on the accuracy and precision of MRI measurements are needed. Taking T1/T2 quantification as an example, even for the inversion recovery spin echo (IR-SE) / multiple single-echo spin echo MRI sequences which are considered as the gold standard for T1/T2 quantification, there exist variations of 2\% - 9\% on the measured relaxation times~\cite{jiang2017repeatability}. Such anatomical measurement uncertainties and model imperfections may weaken the advantage and clinical impact of providing continuous T1/T2 values using network based MRF techniques to some extent. Therefore, improving the accuracy of gold standard approaches in the future would contribute to making the most of the potential of neural networks in the MRF domain.

\section{Conclusion}
\label{sec:Conclusion}

\vspace{-0.2cm}

We proposed a hybrid deep MRF approach which combines low-rank based signature restoration with learning-based parameter restoration. In our approach, a low-rank based de-aliasing method is used to restore clean signatures from subsampled k-space measurements. Then, a 1D deep nonlocal residual CNN is developed for efficient signature-to-parameter mapping, replacing the time-consuming dictionary matching operation in conventional MRF techniques. Our approach demonstrates advantages in terms of inference speed, accuracy and storage requirements over competing MRF methods as no dictionary is needed for recovery.


text:
\section*{Introduction}

Many-body perturbation theory for condensed-matter physics allows the
Green's function of a system of interacting electrons to be formulated in
a picture of independent quasiparticles moving in an effective potential.
The key quantity that incorporates the contributions of dynamic exchange
and correlation to this effective potential is the self-energy operator
$\Sigma$, which in general is both nonlocal and energy dependent. It is
itself a functional of the Green's function that can formally be
expressed through an infinite series of Feynman diagrams. In practice,
however, approximations for the functional form are required, of which the
most popular is the $GW$ approximation introduced by Hedin.\cite{Hed65}
This approach replaces the infinite series of Feynman diagrams for the
self-energy operator by a single Fock-like diagram that is the product of
the Green's function $G$ and the dynamically screened Coulomb interaction
$W$, hence the name of the scheme. Originally the $GW$ approximation was
devised to be applied self-consistently on an equal footing with the
Hartree and Hartree-Fock approximations, in the sense that the Green's
function used to generate the self-energy be identical to that obtained
from it, but the computational cost proved prohibitive. Thus $GW$
self-energies are traditionally calculated from a zeroth order Hartree or
LDA Green's function, although a small degree of self-consistency is
sometimes introduced by rigidly shifting the zeroth order Green's function
on the energy axis in such a way that its chemical potential agrees with
that of the $GW$ Green's function derived from it. It is in this fashion
that the $GW$ approximation has in the past been successfully applied to
a wide range of materials including semiconductors,\cite{Hyb85,God86}
simple metals,\cite{Nor87} and transition metals.\cite{Ary92} The first
fully self-consistent $GW$ calculations for model systems were reported only
recently,\cite{Gro95,Shi96} but these consistently showed an undesired loss
of structural features in the Green's function in disagreement with
experiment.\cite{Bar96}

As future $GW$ calculations for realistic materials are therefore unlikely
to adopt a fully self-consistent approach, the question of particle number
conservation, long since a subject of debate in this context, has gained
renewed significance. Since the imaginary part of the Green's function
is directly related to the density of states as well as to the real-space
charge density, it must integrate to the correct number of electrons
contained in the system, a condition frequently used to determine the
numerical accuracy in practical implementations. However, while the fully
self-consistent $GW$ scheme is known to conserve the exact particle
number,\cite{Bay61} the same, while tacitly assumed, has never actually been
proven for the approximation used in practice, in which the self-energy is
calculated from a zeroth order Green's function. Numerical simulations
appeared to corroborate this positive assumption, with deviations below one
percent fully within the range expected due to systematic numerical errors,
and this has already prompted occasional confirmative conjectures in the
literature.\cite{Ary92} However, in this paper we shall show that upon
elimination of these errors, which arise from numerical integration,
transformation to Fourier space with a finite broadening of the
quasiparticle peaks, alignment of the chemical potentials only within
second order perturbation theory, etc., there still remains a genuine
albeit numerically small violation of particle number conservation, as was
previously demonstrated by accurate molecular calculations for the
analogous scheme based on time-dependent Hartree-Fock rather than Hartree
theory that is frequently employed in quantum chemistry.\cite{Res76} We
also demonstrate that restoring a degree of self-consistency by means of an
appropriate shift in the zeroth order Green's function, which in practice
is often omitted because of its negligible influence on the band structure,
further significantly reduces the deviation. The significance of this result
is that it provides a way to minimize the fundamental limit on the
accuracy with which quantities such as the charge density or the total
energy can be calculated in the $GW$ approximation.

\section*{Description of model system}

For our proof we consider as an analytically solvable counter-example a
four-site Hubbard cluster with tetrahedral symmetry, populated by two
electrons. To the best of our knowledge this is in fact the first model
system forwarded in the literature for which the Green's function and
particle number in the $GW$ approximation can be calculated strictly
analytically, without any numerical errors. Its significance for theoretical
investigations therefore stretches beyond the objective of this paper, and
for the benefit of the reader we will thus present all derivations in
sufficient detail. The Hamiltonian is
\begin{equation}
{\cal H} = \left(\! \epsilon - \frac{U}{2} \right) \sum_{{\bf R},\sigma}
\hat{n}_{{\bf R}\sigma} -t \!\!\sum_{{\bf R},{\bf R}',\sigma}\!\!
c^\dagger_{{\bf R}\sigma} c_{{\bf R}'\sigma} + U \sum_{\bf R}
\hat{n}_{{\bf R}\uparrow} \hat{n}_{{\bf R}\downarrow} ,
\end{equation}
where $c^\dagger_{{\bf R}\sigma}$ and $c_{{\bf R}\sigma}$ are the creation
and annihilation operators for an electron at site ${\bf R}$ with spin
$\sigma$, and $\hat{n}_{{\bf R}\sigma} = c^\dagger_{{\bf R}\sigma}
c_{{\bf R}\sigma}$. We set $\epsilon = 4 t$ and $t = 1$ while leaving the
onsite interaction strength $U$ variable. To construct the $GW$ self-energy
we start from a zeroth order Green's function in the Hartree approximation,
noting that the exchange-correlation potential in a corresponding
density-functional treatment is a mere constant due to spatial symmetries
and would not affect the following arguments. Spatial symmetries and
degeneracy in $\sigma$ also require a uniform fractional site occupation of
one half and accordingly a Hartree potential of $U/2$ on all sites.
Analytic diagonalization of the Hartree $4 \times 4$ Hamiltonian matrix
${\cal H}^{\rm H}_{{\bf RR}'} = \epsilon \delta_{{\bf RR}'} - t$ for each
spin orientation then yields a nondegenerate ground state at zero energy and
a threefold degenerate excited state at energy $\epsilon$. With the electrons
in the ground state $| 0 \rangle$ the zeroth order Green's function becomes
\begin{equation}
G^{\rm H}_{{\bf RR}'}(\omega) = \frac{\langle {\bf R} | 0 \rangle \langle 0
| {\bf R}' \rangle}{\omega - i\eta} + \sum_{\nu=1}^3 \frac{\langle {\bf R} |
\nu \rangle \langle \nu | {\bf R}' \rangle}{\omega - \epsilon + i\eta} ,
\end{equation}
where we have omitted the spin index. $\eta$ denotes a positive
infinitesimal. Using the relation $\langle {\bf R} | 0 \rangle \langle 0 |
{\bf R}' \rangle = 1/4$ for the components of the ground-state vector we
proceed to calculating the polarization propagator in the random-phase
approximation (RPA), defined through
\begin{eqnarray}
P^{\rm RPA}_{{\bf RR}'}(\omega)
&=& -2 \frac{i}{2\pi} \int\!\! G^{\rm H}_{{\bf RR}'}(\omega+\omega')
G^{\rm H}_{{\bf R}'{\bf R}}(\omega) \,d\omega' \nonumber \\
&=& \frac{1}{2} \sum_{\nu=1}^3 \langle {\bf R} | \nu \rangle
\langle \nu | {\bf R}' \rangle \!\left\{ \frac{1}{\omega - \epsilon + i\eta}
- \frac{1}{\omega + \epsilon - i\eta} \right\} , \nonumber \\
&& \mbox{}
\end{eqnarray}
including a factor 2 for the spin summation. The polarization propagator is
diagonal in eigenvector space, and can thus be analytically inverted to
yield the screened interaction
\begin{eqnarray}
W^{\rm RPA}_{{\bf RR}'}(\omega)
&=& U \delta_{{\bf RR}'} + U \sum_{{\bf R}''} P^{\rm RPA}_{{\bf RR}''}(\omega)
W^{\rm RPA}_{{\bf R}''{\bf R}'}(\omega) \nonumber \\
&=& U \delta_{{\bf RR}'} + \frac{\epsilon U^2}{2 z} \sum_{\nu=1}^3
\langle {\bf R} | \nu \rangle \langle \nu | {\bf R}' \rangle \nonumber \\
&& \times \left\{ \frac{1}{\omega - z + i\eta}
- \frac{1}{\omega + z - i\eta} \right\} ,
\end{eqnarray}
with $z = [\epsilon (\epsilon + U)]^{1/2}$. A similar analytic expression
for the screening in a related model system was previously given in Ref.\
\onlinecite{Ver95}. The self-energy in the $GW$ approximation is a convolution
of the zeroth order Green's function and the RPA screened interaction,
\begin{eqnarray} \label{sigma}
\Sigma^{GW}_{{\bf RR}'}(\omega)
&=& \frac{i}{2\pi} \int\!\! G^{\rm H}_{{\bf RR}'}(\omega-\omega')
W^{\rm RPA}_{{\bf RR}'}(\omega') e^{i \delta \omega'} \,d\omega' \nonumber \\
&=& -\frac{U}{4} \delta_{{\bf RR}'} + \frac{3 \epsilon U^2}{8 z}
\frac{\langle {\bf R} | 0 \rangle \langle 0 | {\bf R}' \rangle}{\omega - z
- \epsilon + i\eta} \nonumber \\
&&+ \frac{\epsilon U^2}{8 z} \sum_{\nu=1}^3 \langle {\bf R} | \nu \rangle
\langle \nu | {\bf R}' \rangle \nonumber \\
&& \times \left\{ \frac{2}{\omega - z - \epsilon + i\eta}
+ \frac{1}{\omega + z - i\eta} \right\} ,
\end{eqnarray}
where $\delta$ denotes a positive infinitesimal. We have written
$\Sigma^{GW}$ in such a way as to emphasize that it is diagonal in the
eigenvectors $| \nu \rangle$ of the initial Hartree system. On the other
hand, diagonalization of the full Hamiltonian including the self-energy
yields the quasiparticle states of the interacting electron system, so for
our model the two are in fact identical, although the corresponding energy
eigenvalues are not. In this way we can calculate the chemical potential
$\mu$ of the interacting electron system exactly from its true quasiparticle
properties. By definition the chemical potential is identical to the energy
eigenvalue of the highest occupied quasiparticle state, which is zero for
the original Hartree system, and at the $GW$ level is given implicitly
through the self-energy correction,
\begin{equation}
\mu = 0 + \Sigma^{GW}_{\nu=0}(\mu-\tilde{\omega}) .
\end{equation}
Here we have allowed for the possibility of using a self-energy derived
from a zeroth order Green's function whose chemical potential has been
shifted by $\tilde{\omega}$ on the energy axis. To simulate the effect of
a self-consistent calculation we determine the shift by requiring the
chemical potentials of the shifted zeroth order Green's function and the
$GW$ Green's function obtained from it to be identical. While this equation
is usually solved within second order perturbation theory, the simple form
of the self-energy for our model (\ref{sigma}) allows us to derive the
exact analytic solution
\begin{equation} \label{shift}
\tilde{\omega} = -\frac{U}{4} - \frac{3 \epsilon U^2}{8 z (z+\epsilon)} ,
\end{equation}
which correctly approaches zero as $U \to 0$. Accordingly we solve Dyson's
equation
\begin{equation} \label{dyson}
G^{GW}_{\nu}(\omega) = G^{\rm H}_{\nu}(\omega) + G^{\rm H}_{\nu}(\omega)
\Sigma^{GW}_{\nu}(\omega-\tilde{\omega}) G^{GW}_{\nu}(\omega)
\end{equation}
in eigenvector space both for $\tilde{\omega} = 0$ and $\tilde{\omega}$ as
in Eq.\ (\ref{shift}), and compare the results. In this diagonal form
Dyson's equation is analytically solvable. For $\nu = 0$ the self-energy
(\ref{sigma}) contains one pole, which adds a satellite to the
quasiparticle peak of the Hartree Green's function and yields a quadratic
equation for the positions of the poles of $G^{GW}$. Similarly the
quasiparticle and satellite structure of the $\nu > 0$ matrix elements is
obtained from the zeroes of a third order polynomial, reflecting the richer
spectrum in the self-energy. The $GW$ Green's function therefore takes the
analytic form
\begin{eqnarray} \label{greensfun}
\lefteqn{ G^{GW}_{{\bf RR}'}(\omega+\tilde{\omega}) } \nonumber \\
&=& \langle {\bf R} | 0 \rangle \langle 0 | {\bf R}' \rangle
\frac{\omega-z-\epsilon}{(\omega-x_1-i\eta) (\omega-x_2+i\eta)} \nonumber \\
&&+ \sum_{\nu=1}^3 \langle {\bf R} | \nu \rangle \langle \nu | {\bf R}'
\rangle \nonumber \\
&&\times \frac{(\omega-z-\epsilon) (\omega+z)}{(\omega-y_1-i\eta)
(\omega-y_2+i\eta) (\omega-y_3+i\eta)} ,
\end{eqnarray}
where we have defined the symbols
\begin{equation}
x_{1,2} = \frac{\displaystyle z+\epsilon-\tilde{\omega}-\frac{U}{4}}{2}
\mp \sqrt{ \left( \frac{\displaystyle z+\epsilon+\tilde{\omega}
+\frac{U}{4}}{2} \right)^2 + \frac{3 \epsilon U^2}{8z} }
\end{equation}
together with $y_1 = -(b/3) - 2r \cos(\phi/3)$ and $y_{2,3} = -(b/3) + 2r
\cos[(\pi \mp \phi)/3]$ as well as the auxiliary quantities
$r = \pm \sqrt{|p|}$, $\phi = \arccos(q/r^3)$, $q = b^3/27 - bc/6 + d/2$,
$p = (3c-b^2)/9$, and the polynomial coefficients
\begin{mathletters}
\begin{eqnarray}
b &=& -2\epsilon + \tilde{\omega} + \frac{U}{4} , \\
c &=& -z (z+\epsilon) + \epsilon \left( \epsilon-\tilde{\omega}
-\frac{U}{4} \right) - \frac{3 \epsilon U^2}{8z} , \\
d &=& z (z+\epsilon) \left( \epsilon-\tilde{\omega}-\frac{U}{4} \right)
+\frac{(\epsilon - z) \epsilon U^2}{8z} .
\end{eqnarray}
\end{mathletters}

\section*{Results}

From the Green's function in the $GW$ approximation (\ref{greensfun}) the
particle number may be obtained by an analytic contour integration along a
path closed across the upper half plane, to sample all occupied states
below the chemical potential. By inspection we note that these are the
states at $x_1$ and $y_1$. For the total particle number we thus obtain
\begin{eqnarray} \label{number}
N &=& 2 \sum_{\bf R} \frac{1}{2 \pi i} \int\!\! G^{GW}_{\bf RR}(\omega)
e^{i \delta \omega} \,d\omega \nonumber \\
&=& \frac{2 (x_1-z-\epsilon)}{x_1-x_2} + \frac{6 (y_1-z-\epsilon)
(y_1+z)}{(y_1-y_2) (y_1-y_3)} ,
\end{eqnarray}
including a factor 2 for spin summation. If the particle number in the $GW$
approximation was conserved, $N$ would have to be a constant with a value
of two. In particular, it would also have to be independent of the
interaction strength $U$, which so far we have not specified, but an
analysis of the expression (\ref{number}) confirms that this is not the
case. In Fig.\ \ref{fig:particlenumber} we show the calculated particle
number as a function of $U$, both with and without applying the shift
$\tilde{\omega}$ in the zeroth order Green's function.
In either case the growing deviation from the true value as $U$ increases
is clearly visible, evidence of a genuine and fundamental violation
of particle number conservation in the $GW$ approximation as usually
applied. In absolute terms, however, we find that the discrepancy is much
reduced when the chemical potentials are aligned in the prescribed way, and
that it depends only very weakly on the interaction strength up to high
values of $U$. For a medium correlation of $U = 4$ the numerical deviation
amounts to an underestimation of merely 0.21\% of the true particle number.

\begin{figure}
\epsfxsize=3.25in \epsfbox{1997_PhysRevB1.eps}
\caption{Analytically calculated particle number in the $GW$ approximation
for our two-electron model system as a function of the interaction strength
$U$. Applying a rigid shift to the zeroth order Green's function such as to
align its chemical potential self-consistently with that of the $GW$
Green's function significantly reduces the violation of particle number
conservation, but fails to remove it completely.}
\label{fig:particlenumber}
\end{figure}

A comparison of the weight factors of $G^{GW}$ with those of the exact
Green's function, calculated by numerical diagonalization of the Hamiltonian
matrix, shows that in the absence of the self-consistency shift
$\tilde{\omega}$ the main error in $N$ stems from a serious overestimation of
the satellite spectrum, while the weight of the quasiparticle peak at $x_1$
deviates relatively little from the correct value up to high correlation
strength. When the self-consistency shift is applied the satellites are
also overestimated, although by a lesser amount. In this case, however, this
is the result of a balanced weight transfer from the quasiparticle peaks that
has little influence on the integrated spectral weight and yields a total
particle number in much better agreement with the correct value.

\section*{Summary}

In summary, we have presented a two-electron model system for which the
particle number from the Green's function in the $GW$ approximation can be
derived strictly analytically, without any additional inaccuracies that have
previously beset numerical calculations. Through an analysis of this model
system we have demonstrated that there is a genuine violation of particle
number conservation in the $GW$ approximation as it is usually applied,
i.e., with the self-energy calculated non-self-consistently from a zeroth
order Green's function. However, we also find that the numerical deviation
from the exact particle number can be kept low even for an extremely strong
correlation by introducing a small degree of self-consistency in the form
of a simple rigid shift of the zeroth order Green's function on the energy
axis in such a way as to align its chemical potential with that of the $GW$
Green's function, which is calculated from it through Dyson's equation.
While the effect of this shift on the quasiparticle band structure may be
small, our results clearly indicate its significance for the calculation of
quantities such as the particle number and charge density that are derived
through integration over the complete spectral function.

\acknowledgements

The author wishes to thank R.~W.\ Godby, M.~M.\ Rieger, and T.~J.\ Pollehn
for useful discussions, and acknowledges financial support from the
Deut\-scher Aka\-de\-mi\-scher Aus\-tausch\-dienst under its HSP III scheme,
the Stu\-di\-en\-stif\-tung des deut\-schen Vol\-kes, the Gott\-lieb
Daimler- und Karl Benz-Stiftung, Pembroke College Cambridge, and the
Engineering and Physical Sciences Research Council.

redpajama_set_name: RedPajamaArXiv