Accelerated MR Fingerprinting with Low-Rank and Generative Subspace Modeling

MR Fingerprinting experiments involve acquiring a sequence of contrast-weighted MR images $\{I_{m}(\mathbf{x})\}_{m=1}^{M}$, in which each image $I_{m}(\mathbf{x})$ can be parameterized as

for $m=1,\cdots,M$, where $\rho(\mathbf{x})$ denotes the proton density, $T_{1}(\mathbf{x})$ denotes the spin-lattice relaxation time, $T_{2}(\mathbf{x})$ denotes the spin-spin relaxation time, and $\phi_{m}(\cdot)$ denotes the contrast-weighting function at the $m$th repetition time (TR), which is determined by the magnetization dynamics governed by the Bloch equation [40].

In this work, we consider a discrete image model which assumes that each contrast-weighted image $I_{m}(\mathbf{x})$ can be completely represented by its values at $N$ spatial locations $\left\{\mathbf{x}_{n}\right\}_{n=1}^{N}$. Accordingly, we can represent the contrast-weighted time-series images $\left\{I_{m}(\mathbf{x})\right\}_{m=1}^{M}$ by the following Casorati matrix [41]:

In MR Fingerprinting, the imaging equation can be written as

for $m=1,\cdots,M$ and $c=1,\cdots,N_{c}$. Here $\mathbf{I}_{m}\in\mathbb{C}^{N}$ denotes the contrast-weighted MR image at the $m$th TR, $\mathbf{F}_{m}\in\mathbb{C}^{P_{m}\times N}$ denotes the undersampled Fourier-encoding matrix at the $m$th TR, $\mathbf{S}_{c}\in\mathbb{C}^{N\times N}$ is a diagonal matrix whose diagonal entries contain the coil sensitivities from the $c$th coil, $\mathbf{d}_{m,c}\in\mathbb{C}^{P_{m}}$ denotes the measured $\mathbf{k}$-space data from the $c$th receiver coil at the $m$th TR, and $\mathbf{n}_{m,c}\in\mathbb{C}^{P_{m}}$ denotes the measurement noise that is assumed to be complex Gaussian noise.

For notational simplicity, we can rewrite (2) in a more compact form as follows:

where $\mathbf{d}_{c}\in\mathbb{C}^{P}$ contains all the measured $\mathbf{k}$-space data from the $c$th coil, $\mathbf{F}\in\mathbb{C}^{\tilde{N}\times N}$ denotes the fully-sampled Fourier encoding matrix, $\Omega(\cdot):\mathbb{C}^{\tilde{N}\times M}\rightarrow\mathbb{C}^{P}$ denotes the undersampling operator that acquires the $\mathbf{k}$-space data for each contrast-weighted image and then concatenates them into the data vector $\mathbf{d}_{c}$, $\mathbf{n}_{c}\in\mathbb{C}^{P}$ contains the measurement noise from the $c$th coil, and $P=\sum_{m=1}^{M}P_{m}$.

In MR Fingerprinting experiments, $\mathbf{k}$-space data are highly-undersampled, i.e., $P\ll NM$. A significant technical challenge arises from reconstructing contrast-weighted time-series images from highly-undersampled data. To ensure high-quality reconstruction, we impose a low-rank constraint on $\mathbf{C}$, i.e., $\mathrm{rank}(\mathbf{C})\leq L\ll\mathrm{min}\left\{N,M\right\}$, by exploiting the strong spatiotemporal correlation of contrast-weighted images in MR Fingerprinting. The low-rank constraint reduces the number of degrees of freedom in $\mathbf{C}$ to $L(M+N-L)$ [42], which is often much less than the total number of entries of $\mathbf{C}$. This enables image reconstruction from highly-undersampled $\mathbf{k}$-space data.

While there are many ways of imposing a low-rank constraint [42], here we use a low-rank matrix factorization [43, 44, 30], i.e., $\mathbf{C}=\mathbf{U}\mathbf{V}$, where $\mathbf{U}\in\mathbb{C}^{N\times L}$ and $\mathbf{V}\in\mathbb{C}^{L\times M}$. Note that the columns of $\mathbf{U}$ and the rows of $\mathbf{V}$ respectively span the spatial subspace and the temporal subspace of $\mathbf{C}$. Further, we can incorporate a subspace constraint by pre-estimating $\mathbf{V}$ from an ensemble of magnetization evolutions simulated based on spin physics [15] or some physically-acquired navigator data [45].

While the low-rank and subspace reconstruction significantly outperforms the dictionary-matching based conventional reconstruction [1], it often suffers from an ill-conditioned model fitting issue, due to the undersampling of the temporal subspace of the model. This degrades the performance of the low-rank and subspace reconstruction [15], as the data acquisition length becomes short and/or SNR becomes low.

To address the issue, here we incorporate a deep generative prior into the low-rank and subspace reconstruction. While the effectiveness of a pre-trained deep generative prior [46] has been demonstrated in solving inverse problems, it often requires large training data sets, which are often very expensive to acquire in MR Fingerprinting. Recently, it has been demonstrated that with the inductive bias of a deep convolutional architecture, an untrained generative neural networks (e.g., deep image prior [47] and deep decoder [48, 49]) can serve as an effective spatial prior. Here we use an untrained generative neural network to represent the spatial subspace of a low-rank model, which serves as a regularizer for the low-rank and subspace reconstruction. Mathematically, we have $\mathbf{U}=\mathbf{G}_{\boldsymbol{\theta}}(\mathbf{z})$, where $\mathbf{G}_{\boldsymbol{\theta}}(\mathbf{z}):\mathbb{R}^{N_{0}}\rightarrow\mathbb{R}^{N\times L}$ is a generative neural network, $\mathbf{z}\in\mathbb{R}^{N_{0}}$ is a low-dimensional latent random vector, and $\boldsymbol{\theta}$ contains the trainable parameters of the neural network. Here we illustrate an example generative convolutional neural network architecture in Fig. 1, which consists of convolutional, up-sampling, ReLu, and batch normalization layers.

Incorporating the low-rank and subspace model with the deep generative prior, the image reconstruction problem can be formulated as follows:

where $\hat{\mathbf{V}}$ denotes the pre-estimated temporal subspace of the low-rank model by singular value decomposition. After solving (4), we can reconstruct the contrast-weighted time-series images by $\hat{\mathbf{I}}=\hat{\mathbf{U}}\hat{\mathbf{V}}$, from which we then estimate MR tissue parameters via a voxel-wise pattern matching. Next, we describe an optimization algorithm to solve (4).

The proposed formulation in (4) results in a non-convex optimization problem. Here we describe an algorithm based on variable splitting and the ADMM method [34, 33] to solve this problem. First, we introduce a set of auxiliary variables $\left\{\mathbf{H}_{c}\right\}_{c=1}^{N_{c}}$ to split the coil sensitivities from the objective function. This converts (4) into the following equivalent constrained optimization problem:

for $c=1,\cdots,N_{c}$.

Second, we form the augmented Lagrangian associated with (5) as follows:

where $\left\{\mathbf{\Lambda}_{c}\right\}_{c=1}^{N_{c}}$ and $\mathbf{\Gamma}$ are the Lagrangian multipliers, $\mu_{1},\mu_{2}>0$ are the penalty parameters, and the operator $\textbf{Re}(\cdot)$ takes the real part of a complex number.

Next, we iteratively minimize (6) by solving the following three subproblems:

and updating the Lagrangian multipliers as follows:

until the change of the solution is less than a pre-specified tolerance or the maximum number of iterations is reached. Next, we describe the detailed solutions to the subproblems in (7) - (9).

1) Solution to (7): Note that the subproblem in (7) can be written as

where $\mathbf{Q}_{c}^{(k)}=\mathbf{S}_{c}\mathbf{U}^{(k)}\hat{\mathbf{V}}-\mathbf{\Lambda}_{c}^{(k)}/{\mu_{1}}$. This is a linear least-squares problem separable for each $\mathbf{H}_{c}$, i.e.,

for $c=1,\cdots,N_{c}$. Further, it can be shown that (13) is separable for each column of $\mathbf{H}_{c}$, i.e.,

for $m=1,\cdots,M$, where $\mathbf{h}_{m,c}$ and $\mathbf{q}^{(k)}_{m,c}$ are the $m$th columns of $\mathbf{H}_{c}$ and $\mathbf{Q}^{(k)}_{c}$, respectively. The solution to (14) can be written as

where $\mathds{1}_{N}\in\mathbb{C}^{N\times N}$ is the identity matrix. While (15) can be solved iteratively (e.g., with the conjugate gradient algorithm [50]), we can solve it in a more efficient way with a direct matrix inversion, given that the same coefficient matrix $(\mathbf{F}_{m}^{H}\mathbf{F}_{m}+\mu_{1}\mathds{1}_{N}/2)$ is shared by a number of linear systems of equations at different iterations of the algorithm¹¹1In standard MR Fingerprinting acquisitions (e.g., [1, 51]), the same $\mathbf{k}$-space trajectory repeats at different TRs. Thus, the same coefficient matrix $(\mathbf{F}_{m}^{H}\mathbf{F}_{m}+\mu_{1}\mathds{1}_{N}/2)$ also occurs at different $m$ in each iteration of the algorithm..

Specifically, consider the inversion of the coefficient matrix of (15), and, by applying the matrix inversion lemma [52], we have

Note that $\mathbf{F}_{m}\mathbf{F}^{H}_{m}\in\mathbb{C}^{P_{m}\times P_{m}}$ is a small-scale matrix that can be pre-computed and saved, given that $\mathbf{k}$-space data are highly-undersampled for each $\mathbf{I}_{m}$ (i.e., $P_{m}\ll\tilde{N}$). Here we only need to compute the matrix $\left[\left(\mu_{1}/2\right)^{2}\mathds{1}_{N}+\mu_{1}/2\mathbf{F}_{m}\mathbf{F}_{m}^{H}\right]^{-1}$ once and save it for later use at different iterations of the algorithm. Consequently, the matrix inversion in (LABEL:eq:subproblem1-4) only encompasses two non-uniform Fourier transforms and some small-scale matrix-vector multiplications that can be calculated very efficiently.

2) Solution to (8): The subproblem in (8) can be written as

which can be solved by

where

Noting that the coefficient matrix $\left(\sum_{c=1}^{N_{c}}{\mathbf{S}}_{c}^{H}{\mathbf{S}}_{c}+\mu_{2}/\mu_{1}\mathds{1}_{N}\right)$ is diagonal, we can solve (18) very efficiently via an entry-by-entry inversion.

3) Solution to (9): The subproblem in (9) can be written as

This is a large-scale nonlinear optimization problem similar to a neural network training problem. Here we solve it by the Adam algorithm [53].

For clarity, we summarize the overall algorithm in Algorithm 1. Given that (5) is a non-convex optimization problem, the performance of the proposed algorithm depends on initialization. Here we initialize $\mathbf{U}$ using the low-rank and subspace reconstruction, which is often computationally efficient. Additionally, we initialize the Lagrangian multipliers $\mathbf{\Lambda}_{c}$ and $\mathbf{\Gamma}$ with zero matrices and vector, and initialize $\boldsymbol{\theta}$ with a random vector. In the next section, we will illustrate the performance of the above initialization scheme.

We conducted simulations on a numerical brain phantom to evaluate the proposed method. The numerical phantom simulates a single-channel MR Fingerprinting experiment, which has the same setup as the one in [9]. Specifically, we took the $T_{1}$, $T_{2}$, and proton density maps from the Brainweb database [54] as the ground truth. We performed Bloch simulations to generate contrast-weighted time-series images by an inversion recovery fast imaging with steady state precession (IR-FISP) sequence [51]. Here we used the same set of flip angles and repetition times as in [51] for Bloch simulations. We set the field-of-view (FOV) as $300\times 300$ mm² and the corresponding matrix size as $256\times 256$.

We acquired highly-undersampled $\mathbf{k}$-space data using one spiral interleaf per TR and a set of fully-sampled $\mathbf{k}$-space data consists of $48$ spiral interleaves. We added complex white Gaussian noise to measured data. Here we define the signal-to-noise ratio (SNR) as $\mathrm{SNR}=20\log_{10}(s/\sigma)$, where $s$ denotes the average signal intensity within a region of the white matter in the first contrast-weighted image, and $\sigma$ denotes the standard deviation of the noise.

We performed image reconstruction using the low-rank and subspace method [15] and the proposed method. Here we manually selected the rank $L=6$ for both methods for optimized performance. For the proposed method, we employed a generative convolutional neural network architecture illustrated in Fig. 1. We solved the subproblem (19) using the Adam algorithm with the learning rate $=0.01$ and the number of iterations $=300$. In addition, we manually selected the penalty parameters of the ADMM algorithm for optimized performance. We implemented both reconstruction methods on a Linux server with dual Intel Xeon Platinum $8362$R processors (each with $2.80$ GHz and $32$ cores), $2.00$ TB RAM, and dual NVIDIA A$100$ GPUs (each with $80$ GB RAM) running MATLAB^® R2022b and Deep Learning Toolbox^TM.

After image reconstruction, we estimated MR tissue parameter maps by dictionary matching. Here the same dictionary was used for both methods based on the following discretization scheme: 1) the range of possible $T_{1}$ values was in $[100,3000]$ ms, with an incremental step of $10$ ms in the range of $[100,1500]$ ms and an incremental step of $20$ ms in the range of $[1501,3000]$ ms; and 2) the range of possible $T_{2}$ values was in $[20,350]$ ms, with an incremental step of $1$ ms in the range of $[20,200]$ ms and an incremental step of $2$ ms in the range of $[201,350]$ ms.

To assess the accuracy of reconstructed contrast-weighted MR images or parameter maps, we used the normalized root-mean-square error (NRMSE) defined as $\|\mathbf{g}-\hat{\mathbf{g}}\|_{2}/\|\mathbf{g}\|_{2}$,²²2The regions associated with the background, skull, scalp, and cerebrospinal fluid were not included into the NRMSE calculation, given that they are often not regions of interest for neuroimaging studies. where $\mathbf{g}$ and $\hat{\mathbf{g}}$ respectively denote the ground truth and the reconstruction.

We evaluated the performance of the proposed method for in vivo MR Fingerprinting experiments. Specifically, we conducted in vivo experiments on a 3T Siemens Magnetom Prisma scanner (Siemens Medical Solutions, Erlangen, Germany) equipped with a 20-channel receiver coil. A healthy volunteer was scanned using the IR-FISP sequence [51] with the approval from the Institutional Review Board. For in vivo experiments, we used the same set of acquisition parameters and spiral trajectories as in [51], and other relevant imaging parameters include FOV$\ =300\times 300$ mm², matrix size$\ =256\times 256$, and slice thickness $\ =5$ mm.

Apart from the undersampled MR Fingerprinting experiments, we also acquired a set of fully-sampled data with the acquisition length $M=1000$ following the same procedure as in [9], and then estimated $T_{1}$, $T_{2}$, and proton density maps as our references. Moreover, we conducted a calibration scan using a vendor-provided gradient echo sequence to estimate the coil sensitivities maps.

We performed the low-rank and subspace reconstruction and the proposed method for in vivo data. For the above two methods, we used the same rank value (i.e., $L=8$). In addition, we used the same dictionary for subspace estimation and dictionary matching as described in the simulations for both methods. For the proposed method, we used the same neural network architecture as described in the simulations, and manually selected the penalty parameters for optimized performance.

Fig. 9 shows the results from the low-rank and subspace reconstruction and the proposed method at the acquisition length $M=500$. As can be seen, the proposed method outperforms the low-rank and subspace reconstruction for the $T_{1}$, $T_{2}$, and proton density maps, which is consistent with the simulation results. This further confirms the role of a deep generative neural network as a regularizer for improving the low-rank and subspace reconstruction.