Accelerating Non-Cartesian MRI Reconstruction Convergence using k-space Preconditioning

TL;DR

This paper introduces a k-space preconditioning method for MRI reconstruction that accelerates convergence without sacrificing accuracy or increasing per-iteration computation, demonstrated to converge in about ten iterations.

Contribution

The paper presents a novel dual formulation-based k-space preconditioning approach that improves MRI reconstruction speed without added computational complexity.

Findings

Converges in about ten iterations in practice.

Outperforms existing methods in convergence speed.

Does not require inner loops or extra computation per iteration.

Abstract

We propose a k-space preconditioning formulation for accelerating the convergence of iterative Magnetic Resonance Imaging (MRI) reconstructions from non-uniformly sampled k-space data. Existing methods either use sampling density compensations which sacrifice reconstruction accuracy, or circulant preconditioners which increase per-iteration computation. Our approach overcomes both shortcomings. Concretely, we show that viewing the reconstruction problem in the dual formulation allows us to precondition in k-space using density-compensation-like operations. Using the primal-dual hybrid gradient method, the proposed preconditioning method does not have inner loops and are competitive in accelerating convergence compared to existing algorithms. We derive l2-optimized preconditioners, and demonstrate through experiments that the proposed method converges in about ten iterations in practice.

Tables

Table 1. TABLE I: Computation time for constructing preconditioners

	Liver	Cardiac	Brain
Circulant precond.	0.0974 s	0.0502 s	0.0147 s
MC k-space precond.	0.231 s	0.117 s	0.0334 s

Equations

$${\mathbf{y}}_c[i]=\frac{1}{\sqrt{N}}\sum _{n=0}^{N-1}{\mathbf{s}}_c[n]\mathbf{x}[n]e^{-ı2\pi f_{i}n/N}+{\mathbf{w}}_c[i]$$

$$\mathbf{y}=\mathbf{Ax}+\mathbf{w}$$

$$\underset{\mathbf{x}}{\min }\frac{1}{2}\Vert \mathbf{Ax}-\mathbf{y}{\Vert }_2^2+g(\mathbf{x})$$

$${\mathbf{x}}^{(k+1)}={\text{prox}}_{\alpha g}({\mathbf{x}}^{(k)}-\alpha {\mathbf{A}}^H(\mathbf{A}{\mathbf{x}}^{(k)}-\mathbf{y}))$$

$${\mathbf{x}}^{(k+1)}={\text{prox}}_{\alpha g}({\mathbf{x}}^{(k)}-\alpha {\mathbf{A}}^H\mathbf{D}(\mathbf{A}{\mathbf{x}}^{(k)}-\mathbf{y}))$$

$$\underset{\mathbf{x}}{\min }\frac{1}{2}\Vert {\mathbf{D}}^{1/2}(\mathbf{Ax}-\mathbf{y}){\Vert }_2^2+g(\mathbf{x})$$

$${\mathbf{x}}^{(k+1)}={\text{prox}}_{\alpha g,\mathbf{P}}({\mathbf{x}}^{(k)}-\alpha \mathbf{P}{\mathbf{A}}^H(\mathbf{A}{\mathbf{x}}^{(k)}-\mathbf{y}))$$

$$\mathbf{P}=\mathbf{F}\text{diag}(\mathbf{p}){\mathbf{F}}^H$$

$${\text{prox}}_{\alpha g,\mathbf{P}}(\mathbf{z})=\underset{\mathbf{x}}{\mathrm{argmin}}\frac{1}{2\alpha }\Vert {\mathbf{P}}^{-1/2}\mathbf{x}-\mathbf{z}{\Vert }_2^2+g(\mathbf{x}).$$

$$\underset{\mathbf{u}}{\max }-\left(\frac{1}{2}\Vert \mathbf{u}{\Vert }_2^2-\Re ⟨\mathbf{u},\mathbf{y}⟩+g^{\ast }(-{\mathbf{A}}^H\mathbf{u})\right)$$

$${\mathbf{u}}^{(k+1)}$$

$${\mathbf{x}}^{(k+1)}$$

$${\bar{\mathbf{x}}}^{(k+1)}$$

$${\sigma }^{(k)}{\tau }^{(k)}{\lambda }_{\text{max}}(\mathbf{PA}{\mathbf{A}}^H)<1,$$

$${\theta }^{(k)}$$

$${\sigma }^{(k+1)}$$

$${\tau }^{(k+1)}$$

$$g(\mathbf{x})=h(\mathbf{Gx}).$$

$${\mathbf{u}}^{(k+1)}$$

$${\mathbf{v}}^{(k+1)}$$

$${\mathbf{x}}^{(k+1)}$$

$${\bar{\mathbf{x}}}^{k+1}$$

$${\sigma }^{(k)}{\tau }^{(k)}({\lambda }_{\text{max}}(\mathbf{PA}{\mathbf{A}}^H)+{\lambda }_{\text{max}}(\mathbf{G}{\mathbf{G}}^H))<1.$$

$$\mathbf{p}=\underset{\mathbf{p}}{\mathrm{argmin}}{\Vert \text{diag}(\mathbf{p})\mathbf{A}{\mathbf{A}}^H-\mathbf{I}\Vert }_F^2.$$

$${\mathbf{p}[i]}^{-1}=\frac{{\sum }_{j=1}^M|{\mathbf{a}}_i^H{\mathbf{a}}_j{|}^2}{\Vert {\mathbf{a}}_i{\Vert }_2^2}.$$

$${\mathbf{p}[i]}^{-1}$$

$$=\frac{1}{N}\sum _{j=1}^M{\left|\frac{\sin (\pi (f_i-f_j))}{\sin (\pi (f_i-f_j)/N)}\right|}^2$$

$${{\mathbf{p}}_c[i]}^{-1}=\frac{1}{\Vert {\mathbf{s}}_c{\Vert }_2^{2}N}\sum _{j=1}^M\sum _{d=1}^C{\left|\sum _{n=0}^{N-1}{\mathbf{s}}_c[n]{\mathbf{s}}_d^{\ast }[n]e^{-ı2\pi (f_i-f_j)n/N}\right|}^2$$

$${\mathbf{r}}_{cd}[k]$$

$${\left|\sum _{n=0}^{N-1}{\mathbf{s}}_c[n]{\mathbf{s}}_d^{\ast }[n]e^{-ı2\pi (f_i-f_j)n/N}\right|}^2$$

$$=\sum _{k=-N+1}^{N-1}{\mathbf{r}}_{cd}[k]e^{-ı2\pi (f_i-f_j)k/N}$$

$$\mathbf{h}[k]$$

$${{\mathbf{p}}_c[i]}^{-1}$$

$$\underset{\mathbf{p}}{\min }\frac{1}{2}{\Vert \text{diag}(\mathbf{p})\mathbf{A}{\mathbf{A}}^H-\mathbf{I}\Vert }_F^2$$

$$\underset{\mathbf{p}}{\min }\frac{1}{2}\sum _{i=1}^M\sum _{j=1}^M{\left|\mathbf{p}[i]{\mathbf{a}}_i^H{\mathbf{a}}_j-{\mathbf{\delta }}_{ij}\right|}^2,$$

Full text

I Abstract

We propose a k-space preconditioning formulation for accelerating the convergence of iterative Magnetic Resonance Imaging (MRI) reconstructions from non-uniformly sampled k-space data. Existing methods either use sampling density compensations which sacrifice reconstruction accuracy, or circulant preconditioners which increase per-iteration computation. Our approach overcomes both shortcomings. Concretely, we show that viewing the reconstruction problem in the dual formulation allows us to precondition in k-space using density-compensation-like operations. Using the primal-dual hybrid gradient method, the proposed preconditioning method does not have inner loops and are competitive in accelerating convergence compared to existing algorithms. We derive $\ell 2$-optimized preconditioners, and demonstrate through experiments that the proposed method converges in about ten iterations in practice.

II Introduction

Non-Cartesian trajectories provides many advantages over Cartesian sampling based on their unique properties. Spiral [1, 20] and cones trajectories [10], for example, can be designed to traverse k-space efficiently, which make them suitable for fast imaging applications, including coronary imaging [20], and arterial spin labeled perfusion imaging [2]. Many non-Cartesian trajectories, such as radial [17] and projection reconstruction [9] naturally sample low-frequency regions densely, which can provide auto-calibration regions for parallel imaging (PI), and robustness to motion for dynamic applications. Such variable density sampling [30] property is also more adapted to signal energy than uniform sampling, which results in less coherent undersampling artifacts in the wavelet transform domain. Hence, variable density non-Cartesian trajectories are often used with compressed sensing (CS) [18].

On the other hand, reconstructions from non-Cartesian trajectories, especially with PI, are more complex and time-consuming than from Cartesian trajectories. The long reconstruction time is one reason that has limited the clinical adoption of non-Cartesian trajectories. Iterative reconstructions, such as CG-SENSE [25], have to be used for PI, which can often take many iterations to converge. In comparison, the Cartesian SENSE method [26] has an analytic solution that can be efficiently solved in a single step.

One way to make non-Cartesian PI/CS reconstructions more efficient is reducing the number of iterations. In general, the slow convergence of iterative methods is due to the ill-conditioning of the reconstruction problem. For non-Cartesian imaging, such ill-conditioning comes from the variable density sampling distribution in k-space. This often shows up in images as blurring artifacts when the reconstruction has not yet converged. Slow convergence is even more significant for 3D acquisitions and CS reconstructions. For instance, Figure 1 shows the iteration progression for a $\ell 1$-wavelet regularized reconstruction of a 3D ultra-short echo-time (UTE) radial acquisition using the Fast Iterative Soft-Thresholding Algorithm (FISTA) [4], and primal dual hybrid gradient method (PDHG) [6] (also known as the Chambolle-Pock method). Even after 100 iterations, the reconstructed image still displays significant blurring due to slow convergence.

Density compensation [12, 20, 11, 23] is often used as a heuristic to compensate for slow convergence in non-Cartesian iterative reconstruction. It was originally developed for gridding reconstruction, and was mostly designed for Nyquist-sampled trajectories. The use of density compensation in iterative PI reconstruction was first introduced by Pruessmann et al. [25]. While their work showed that in practice density compensation can speed up convergence, reconstruction error was also increased. This is because the data consistency for densely sampled regions is weighted down in the objective function (more detail in Section IV-A).

An alternative to density compensation is preconditioning. Preconditioning has the advantage of preserving the original objective function and hence does not affect the reconstruction accuracy. Many techniques [28, 27, 33, 21, 16, 29] have been proposed for MRI iterative reconstruction as described in detail in Section IV. However, a drawback of existing methods is that they increase the per-iteration computation. In particular, most existing preconditioners have circulant structures, and require at least two additional FFTs per iteration. Moreover, all prior methods require inner loops in their algorithms for non-Cartesian reconstructions, which result in more parameters to tune and can often incur additional computational overhead from initializing inner loop variables.

In this article, we present a framework for speeding up convergence that combines the computational efficiency of density compensation, and the objective preserving property of preconditioning. Similar to the work of Trzasko et al. [29], we consider using efficient operations in k-space for preconditioning. Our contribution is to recognize that a diagonal preconditioner can be applied in k-space more generally by viewing the objective function in the dual formulation. Using PDHG [6], the resulting method with preconditioning does not have inner loops, so it has a similar computational complexity as the vanilla proximal gradient method. Moreover, instead of using off-the-shelf density compensation factors, we derive $\ell 2$-optimized diagonal preconditioners for the MRI forward model. We demonstrate through experiments that the proposed diagonal preconditioner speeds up iterative reconstruction for non-Cartesian imaging, with $\ell 2$, $\ell 1$-wavelet, and total variation regularizations.

III Problem Setup

Throughout this article, we consider the following discrete multi-channel MRI forward model, in which we are given an $N$-size image $\mathbf{x}\in\mathbb{C}^{N}$, $C$-channel sensitivity maps $\{\mathbf{s}_{c}\in\mathbb{C}^{N}\}_{c=1}^{C}$, white Gaussian noise vectors $\{\mathbf{w}_{c}\in\mathbb{C}^{M}\}_{c=1}^{C}$, and k-space measurements $\{\mathbf{y}_{c}\in\mathbb{C}^{M}\}_{c=1}^{C}$ with sampling points $\{f_{i}\}_{i=1}^{M}$ such that

[TABLE]

for $i\in\{1,\ldots,M\}$, and $c\in\{1,\ldots,C\}$. For notation simplicity, we focus on one-dimensional signals. The above model can be succinctly represented as a linear model:

[TABLE]

where $\mathbf{y}\in\mathbb{C}^{MC}$ and $\mathbf{w}\in\mathbb{C}^{MC}$ are stacked versions of $\{\mathbf{y}_{c}\in\mathbb{C}^{M}\}_{c=1}^{C}$ and $\{\mathbf{w}_{c}\in\mathbb{C}^{M}\}_{c=1}^{C}$.

Given the acquired k-space measurements $\mathbf{y}$, we consider the following regularized least squares problem to reconstruct the image:

[TABLE]

where $g(\mathbf{x})$ is the regularization function. Specifically, we will consider three regularization functions as concrete examples: $\ell 2$-norm $(\lambda\|\mathbf{x}\|_{2}^{2}/2)$, $\ell 1$-wavelet $(\lambda\|\mathbf{W}\mathbf{x}\|_{1})$, where $\mathbf{W}\in\mathbb{C}^{N\times N}$ is a unitary wavelet transform operator, and total variation $(\lambda\|\mathbf{G}\mathbf{x}\|_{1})$, where $\mathbf{G}\in\mathbb{C}^{N\times N}$ is a first-order finite difference operator with periodic boundary extension.

Since the image size $N$ is on the order of tens of thousands or more, the above reconstruction problem is in practice only solved approximately using first-order gradient methods. In the following section, we will focus on the proximal gradient method as an example to illustrate the advantages and disadvantages of using density compensation and preconditioners to accelerate convergence.

The proximal gradient method when applied to objective function (3) gives the following update for the $k$th iteration:

[TABLE]

where $\text{prox}_{\alpha g}(\mathbf{z})=\text{argmin}_{\mathbf{x}}\|\mathbf{x}-\mathbf{z}\|_{2}^{2}/(2\alpha)+g(\mathbf{x})$.

The convergence rate depends only on $\mathbf{A}^{H}\mathbf{A}$. More concretely, when $\mathbf{A}$ is not singular, then the step-size $\alpha$ can be chosen so that the convergence rate is inversely proportional to the condition number of $\mathbf{A}^{H}\mathbf{A}$. When $\mathbf{A}$ is singular, then the step-size can be chosen so that the convergence rate is inversely proportional to the maximum eigenvalue of $\mathbf{A}^{H}\mathbf{A}$. For variable density sampling, the condition number or the maximum eigenvalue of $\mathbf{A}^{H}\mathbf{A}$ is much higher than for uniform density sampling for a given undersampling factor and hence results in slow convergence.

IV Prior Arts on Density Compensation and Image Domain Preconditioning

IV-A Density Compensation

One effective heuristic to accelerate convergence for non-Cartesian imaging is incorporating density compensation factors during iterations. Given a diagonal matrix $\mathbf{D}\in\mathbb{C}^{MC\times MC}$ with density compensation factor as diagonals, the heuristic modifies the proximal gradient method as follows:

[TABLE]

The use of density compensation in iterative PI reconstruction was first introduced by Pruessmann et al. [25], and was shown to speed up convergence in practice. Computationally, incorporating density compensation in each iteration costs an additional $O(MC)$ multiplications, adding very little overhead to the overall iteration. However, the main drawback is that such k-space weighting is known to increase reconstruction errors, as implicitly it is solving for a weighted objective function:

[TABLE]

Note that data consistency is weighed down in densely sampled regions. Measurements are essentially thrown away for convergence, resulting in increased reconstruction error, and noise coloring.

IV-B Image-domain Preconditioning

An alternative is to use preconditioning, which only affects the convergence, but not the objective function. Since the objective function is not changed, there is no error penalty for using preconditioners.

In particular, given a preconditioner $\mathbf{P}\in\mathbb{C}^{N\times N}$, the preconditioned proximal gradient method applies:

[TABLE]

The preconditioner $\mathbf{P}$ should be designed to approximate the (pseudo) inverse of $\mathbf{A}^{H}\mathbf{A}$ such that the condition number or maximum eigenvalue of $\mathbf{P}\mathbf{A}^{H}\mathbf{A}$ is much lower than that of $\mathbf{A}^{H}\mathbf{A}$.

Many works have proposed the use of preconditioning to accelerate MRI reconstructions. In particular, it was first described by Sutton et al. [28] for single-channel non-Cartesian imaging in the presence of field inhomogeneities. It was further explored by Ramani et al. [27] for PI-CS reconstructions. Their method leveraged a circulant preconditioner developed by Yagle [34] for Toeplitz systems. Weller et al. [33] considered the non-Cartesian $\ell 1$-SPIRiT [19] method and used an $\ell 2$-optimal circulant preconditioner developed by Chan [8]. Muckley et al. [21] considered FISTA [4] and designed a circulant preconditioner that majorizes the sensing matrix motivated by the convergence criterion. Finally, Koolstra et al. [16] considered the split-Bregman method for Cartesian PI-CS reconstructions and presented a circulant preconditioner that incorporates multi-channel sensitivity maps in the construction of their proposed preconditioner.

A main drawback of the above mentioned preconditioning methods is that they all increase per-iteration computational complexity. This is because to compensate ill conditioning from variable density in k-space, existing preconditioners have to use circulant operators, which cost two FFTs per iteration. That is, existing preconditioners are of the form,

[TABLE]

where $\mathbf{p}\in\mathbb{C}^{N}$ is a Fourier weighting vector, and $\mathbf{F}\in\mathbb{C}^{N\times N}$ is the unitary discrete Fourier transform operator.

A more subtle issue is that all of them require inner loops in their algorithms when incorporating CS with non-Cartesian MRI. In particular, the proximal operator has to be modified to incorporate the preconditioner, which requires inner iterations to solve even when the proximal operator is simple:

[TABLE]

In summary, although existing preconditioners have shown that they can accelerated convergence, their shortcoming lies in the per-iteration increase in complexity.

V Diagonal k-space Preconditioning

Ideally, we would like to develop a preconditioning method that can achieve the computational efficiency of density compensation without changing the objective function.

Recently, Trzasko et al. [29] showed that through an algebraic manipulation, a diagonal preconditioner can be applied in k-space for the least squares sub-problem of ADMM [5]. This enables a different mechanism for preconditioning. In particular, they show that it is possible to use efficient operations in k-space for preconditioning. However, their formulation still required inner loops to solve for the sub-problem.

Here we show that k-space preconditioning without inner loops is achievable by looking at the convex dual problem. Since the reconstruction problem (3) is unconstrained, it must satisfy strong duality. Its corresponding dual problem (see the Supporting Materials for a derivation using the Lagrangian) is given by:

[TABLE]

where $\mathbf{u}\in\mathbb{C}^{m}$ is the dual variable and $g^{*}$ denotes the convex conjugate function of $g$.

Our key observation is that because the dual variable resides in k-space, it is now possible to perform preconditioning in k-space on the dual problem.

In general, one would require primal-dual methods to solve for the primal and dual problems at the same time. In this article, we opt for the PDHG [6] method for an algorithm without inner loops. Other primal-dual reconstruction methods, such as those described in the work of Komodakis et al. [15], can also be used.

Before describing the algorithm, we note that the $\ell 2$-regularized reconstruction is a special case that can efficiently recover the primal variable from the dual problem. While this property is not used in our experiments, we show that this is connected to Trzasko et al. [29] in the Supporting Materials.

V-A Preconditioned PDHG for simple regularizers

Let us first consider the case when the regularizer $g(\mathbf{x})$ is simple, i.e. its proximal operator is easy to compute. An example is the $\ell 1$ wavelet regularization function $g(\mathbf{x})=\lambda\|\mathbf{W}\mathbf{x}\|_{1}$, where $\mathbf{W}$ is the wavelet transform operator.

Following [6] and [24], for each iteration $k$, the diagonal preconditioned version of PDHG for simple proximal operators with a diagonal preconditioner $\mathbf{P}\in\mathbb{C}^{MC\times MC}$ is given by,

[TABLE]

where $\bar{\mathbf{x}}^{k}$ and $\theta^{k}$ are the extrapolated primal variable and extrapolation parameter to provide acceleration, and $\tau^{k}$ and $\sigma^{k}$ are the primal and dual step-size respectively such that

[TABLE]

where $\lambda_{\text{max}}$ denotes the maximum eigenvalue.

One drawback of primal-dual methods is that there are many ways of choosing the step-sizes, which can affect the convergence rate. In this article, we adopt a simple scheme, which reflects that in practice it is impossible to fine-tune parameters to each dataset. We fix $\sigma^{(0)}=1$ and choose $\tau^{(0)}=1/\lambda_{\text{max}}(\mathbf{P}\mathbf{A}\mathbf{A}^{H})$. Since the problem is dual strongly convex [6] (the function involving the dual variable $\mathbf{u}$ is the $\ell 2$ norm), acceleration can be obtained by choosing the subsequent step-sizes appropriately as:

[TABLE]

which can be derived from Theorem 5.1 in [7] with $\mu_{i}=\min_{j}\mathbf{p}_{j}/\mathbf{p}_{i}$ and $\tilde{\sigma}_{(0)}=\min_{j}\mathbf{p}_{j}$.

V-B Preconditioned PDHG for simple regularizers composed with linear operators

Next, we consider regularization functions which consist of a simple function $h$ composed with a linear operator $\mathbf{G}$:

[TABLE]

One example is the total variation regularization $g(\mathbf{x})=\lambda\|\mathbf{G}\mathbf{x}\|_{1}$ with $\mathbf{G}$ being the finite-difference operator.

In this case, PDHG can be modified to perform:

[TABLE]

where $\tau^{(k)}$ and $\sigma^{(k)}$ are the primal and dual step-sizes respectively such that

[TABLE]

Note that unlike the previous case, the problem is no longer dual strongly convex (the function involving the dual variable $\mathbf{v}$ is $h$, which in general is not strongly convex), so the same acceleration scheme cannot be used. We choose $\theta^{(k)}=1$, $\sigma^{(k)}=1$ and $\tau^{(k)}=1/(\lambda_{\text{max}}(\mathbf{P}\mathbf{A}\mathbf{A}^{H})+\lambda_{\text{max}}(\mathbf{G}\mathbf{G}^{H}))$ for all $k$.

VI $\ell 2$ optimized diagonal k-space preconditioners

Now that we know how to precondition in k-space, it becomes clear from the dual problem that the preconditioner should be designed to precondition the matrix $\mathbf{A}\mathbf{A}^{H}$. In this article, we consider diagonal preconditioners to approximate the inverse of the normal operator $\mathbf{A}\mathbf{A}^{H}$ in the least squares sense. The diagonal structure is desired because we want to apply the preconditioner efficiently in k-space, similarly to density compensation. The least squares design, on the other hand, is used here so that we can efficiently compute the preconditioner.

Concretely, we consider a diagonal preconditioner $\mathbf{p}\in\mathbb{C}^{MC}$ such that,

[TABLE]

Let $\mathbf{a}_{i}\in\mathbb{C}^{N}$ denote the $i$th row vector of $\mathbf{A}$. As shown in Appendix A, the general expression for the inverse of the diagonal preconditioner is given by:

[TABLE]

To further look into the preconditioner, we first consider the single-channel case. In this case, $\mathbf{a}_{i}[n]=e^{-\imath 2\pi f_{i}n/N}/\sqrt{N}$, and $\|\mathbf{a}_{i}\|_{2}^{2}=1$. Then the diagonal preconditioner at k-space position $i$, is given by:

[TABLE]

For Cartesian trajectories, the frequency spacing $f_{i}-f_{j}$ are all integers, and hence $\mathbf{p}[i]=1$ for all $i$, which matches our expectation that the condition number for single channel Cartesian imaging systems cannot be improved. For non-Cartesian trajectories, the diagonal preconditioner can be interpreted as calculating density from the sinc squared kernel $|\sin(\pi f)/(N\sin(\pi f/N))|^{2}$.

Moving on to multi-channel, for k-space position $i$ and coil $c$, the row vector is given by $\mathbf{a}_{ci}[n]=\mathbf{s}_{c}[n]e^{-\imath 2\pi f_{i}n/N}/\sqrt{N}$. Hence, we obtain,

[TABLE]

Here we pause to note that the multi-channel preconditioner design is different from density compensation calculations in that we incorporate coil sensitivity maps. One downside is that the multi-channel preconditioner has to be calculated whenever the coil sensitivity maps and the sampling trajectory change. For many clinical applications, the coil sensitivity maps are computed from a pre-scan or estimated from the first scan and used multiple times for a sequence of scans. In this case, the overhead of computing the preconditioner becomes negligible. For applications in which this overhead matters, the single-channel preconditioner (10) can be used as an approximation. As shown in our experiments in Section VII-B, we found that the single-channel preconditioner is competitive in accelerating convergence.

Since the multi-channel preconditioner has to be computed whenever the coil sensitivity maps change, its computation time should not be impractically long. A direct summation implementation takes $O(M^{2}NC^{2})$ computation. In the following, we show that using Fourier transform properties, we can reduce the computational complexity to $O(C^{2}N\log N+CM)$, which makes it comparable to common calibration methods, such as ESPIRiT [31]. Figure 2 provides a high-level diagram of the overall process.

VI-A Efficient computation of the multi-channel preconditioner

First, we note that we can express the squared terms for the multi-channel preconditioner with cross-correlations, which can be computed in $O(C^{2}N\log N)$ using FFTs. Let us define,

[TABLE]

Then

[TABLE]

Next, we note that the preconditioner can be expressed in terms of convolution with the point spread function, which can be approximated using NUFFT with $O(N\log N+M)$ computational complexity. Let us define

[TABLE]

then,

[TABLE]

The final step involves $C$ NUFFTs on the point-wise multiplication of $\mathbf{r}$, and $\mathbf{h}$. Hence we obtain the overall computational complexity to be $O(C^{2}N\log N+CM)$.

VII Experiments

In the spirit of reproducible research, we provide a software package in Python to reproduce the results described in this article. The software package can be downloaded from:

https://github.com/mikgroup/kspace_precond.git

For each regularization, we evaluated the proposed method on three 2D non-Cartesian datasets: 1) a liver dataset acquired with a stack-of-stars trajectory, 2) a brain dataset acquired with a ramp-sampled UTE radial trajectory, and 3) a cardiac dataset with a variable density spiral trajectory. We also applied on one 3D UTE dataset to illustrate the additional benefit of using preconditioners on 3D datasets. These datasets are described in more detail in Section VII-A. The regularization parameters were manually selected to provide reconstructions with good subjective image quality.

We consider three different k-space preconditioners to use with PDHG for comparisons: 1) the Pipe-Menon density compensation factor [23] which was originally designed for gridding reconstruction, 2) single-channel diagonal k-space preconditioner, Eq. (10), which will be abbreivated as SC k-space preconditioner, and 3) multi-channel diagonal k-space preconditioner Eq. (11), which will be abbreivated as MC k-space preconditioner.

For $\ell 2$-regularized reconstruction, conjugate gradient (CG) with and without circulant preconditioning, and PDHG with various k-space preconditionings were applied and compared with $\lambda=0.01$. For $\ell 1$-wavelet regularized reconstruction, FISTA, ADMM with and without circulant preconditioning, and PDHG with various k-space preconditionings were applied and compared with $\lambda=0.001$. The Daubechies-4 wavelet transform was used. For total variation regularized reconstruction, ADMM with and without circulant preconditioning, and PDHG with various k-space preconditionings were applied and compared with $\lambda=0.001$. Anisotropic total variation along horizontal and vertical directions were used. For ADMM, adaptive parameter selection as detailed in [5] was used for the convergence parameter $\rho$ with the initial value set to 1. CG with a tolerance of $1$ were used to solve the sub-problems in ADMM. The parameters for ADMM were selected by sweeping various values over the liver datasets without preconditioning as shown in Supplementary Figure S1 and S2. Similar behaviours were seen for ADMM with circulant preconditioning in Supplementary Figure S3 and S4.

The circulant preconditioner used for CG and ADMM is designed with a least squares metric that takes into account of the sensitivity maps and the regularization function. For single-channel non-Cartesian imaging, the preconditioner coincides with Chan’s preconditioner [8] used in Weller et al. [33]. And for multi-channel Cartesian imaging with total variation regularization, the preconditioner has similar structure as Koolstra et al’s precondtioner [16]. Appendix B contains the detailed derivation.

All methods were implemented in Python using the software packages NumPy [32] and CuPy [22] on a workstation with four Nvidia Titan Xp GPUs. All operations, except the wavelet transform, were run on a single GPU. Normalized coil sensitivity maps were obtained using ESPIRiT [31] on the gridded low frequency k-space. The NUFFT operations were implemented following Beatty et al. [3] with an oversampling ratio of 1.25 and an interpolation kernel width of 4. The maximum eigenvalue $\lambda_{\text{max}}(\mathbf{P}\mathbf{A}\mathbf{A}^{H})$ used for PDHG was estimated using the power method with 30 iterations.

For the 2D datasets, all methods were run for 30 iterations (outer iterations for ADMM), and the objective values were computed for each iteration. For the 3D dataset, all methods were run for 1000 iterations. For $\ell 2$-regularized reconstructions, the computation time for constructing the circulant preconditioner and the proposed preconditioner was also recorded.

Finally, we note that algorithms compared in this work have subtle differences than those used in prior works. This is because each prior work considers a different reconstruction setting and it becomes impossible to evaluate each method in the same way. In particular, a monotonic version of FISTA was used in Ramani and Fessler [27], which has higher computational complexity than the vanilla FISTA compared in this work. In addition, both Ramani and Fessler [27] and Weller et al. [33] consider more complex regularization functions than the ones in this work. Koolstra et al. also considered Cartesian imaging with both $\ell 1$-wavelet and total-variation regularizations, which can exhibit very different convergence behaviors.

VII-A Dataset Details

The liver dataset was acquired with a stack-of-stars trajectory using a 3D T1-FFE sequence (TR/TE 4.35 ms/1.20ms, resolution $1\times 1\times 1.5$ mm3, field-of-view $40\times 40\times 12.5$ cm3, number of spokes $378$, reconstruction matrix size of $401\times 401$ with an effective undersampling factor of about 1.66). The sequence was implemented on a 3T MR system (Philips Healthcare) equipped with a 16-channel torso coil. The center slice was extracted after taking an inverse FFT along the slice direction for the experiments.

The cardiac dataset was acquired with a variable density spiral trajectory on a 1.5 T GE scanner (GE Healthcare, Waukesha, WI) with an 8-channel cardiac coil and the HeartVista RTHawk platform (HeartVista, Los Altos, CA). The trajectory consists of 3 interleaves with 3996 readout points and an effective undersampling factor of about 8. It has an FOV of $32\times 32$ cm2, a reconstruction matrix size of $320\times 320$ and TR of 25.8 ms.

The brain dataset was acquired with a centered-out radial trajectory on a 7.0 T GE clinical scanner (GE Healthcare, Waukesha, WI) with 8-channel head coil. The trajectory has 256 readout points and 754 half-spokes. The following prescribed parameters were used: flip angle of 5 degree, field-of-view $20\times 20$ cm2, in-plane resolution $1\times 1$ mm2, reconstruction matrix size of $255\times 200$, and TE/TR = 3.4 ms/2 s.

The 3D UTE dataset was acquired on a 3 T GE clinical scanner (GE Healthcare, Waukesha, WI) with 8-channel body coil with an optimized bit-reversed ordered radial trajectory [14] using the sequence described in [13]. The following prescribed parameters were used: FOV of $32\times 32\times 32$ cm3, reconstruction matrix size of $356\times 206\times 211$ with an effective undersampling factor of about 3, flip angle of 4 degrees, 1.25 mm isotropic resolution, sampling bandwidth of 62.5 kHz, and readout duration of 1 ms. 75,800 spokes were acquired.

VII-B Results

Figure 3 shows the iteration progression for the $\ell 2$-regularized reconstruction of the liver dataset. Both visually and quantitatively in terms of computation time, the proposed preconditioning methods (SC and MC) along with the vanilla CG converge faster than other methods. The per-iteration computation time for CG with circulant preconditioning is noticeably longer than for other methods, which results in slower convergence rate in terms of computation time. This is also supported by the Supplementary Figures S5 and S6.

Table I shows the computation time for constructing the preconditioners. The construction of the proposed preconditioner is about twice as slow as constructing the circulant preconditioner, but is still in a reasonable range. For applications in which the coil sensitivity maps are calculated from the pre-scan or estimated from the first scan and used multiple times for a sequence of scans, this overhead can be negligible. For applications in which the recalculation time matters, the SC kspace precondioner is preferred as it can be precomputed.

Figure 4 shows the iteration progression for $\ell 1$-wavelet regularized reconstruction of the cardiac dataset. Both visually and quantitatively in terms of computation time, the proposed preconditioning methods (SC and MC) converge faster than other methods. The figure also shows that ADMM with circulant preconditioning is slower than ADMM without preconditioning, which can be due to the additional FFTs. While PDHG with density compensation accelerates convergence with respect to the vanilla PDHG, it also shows excessive oscillations. Supplementary Figures S7 and S8 also support these observations.

Figure 5 shows the iteration progression for total variation regularized reconstruction of the liver dataset. Similar to the $\ell 1$-wavelet regularized reconstruction experiment, the proposed preconditioning methods converge faster than other methods. While ADMM with and without circulant preconditionings are competitive in terms of iteration number, their increased per-iteration time make them slightly slower in reaching an approximate optimal point. Other datasets shown in Supplementary Figures S10 and S9 support the above observations as well.

Finally, the iteration progression for the 3D UTE dataset was shown earlier in Figure 1. Both FISTA and PDHG exhibit extreme blurring even after 100 iterations. In contrast, PDHG with the proposed preconditioner converges in about ten iterations, both visually and quantitatively in terms of minimizing the objective value. This shows that the proposed method can offer an order magnitude speedup in 3D than in 2D.

VIII Discussion

In this article, we presented a preconditioning method using the convex dual formulation. This enables the use of efficient k-space operations as preconditioners and does not modify the objective function. Through experiments, we have demonstrated that the proposed technique indeed accelerates the convergence of non-Cartesian reconstructions.

In particular, we compared the performance of the proposed preconditioning to circulant preconditioning with CG for $\ell 2$-regularized reconstructions and with ADMM for $\ell 1$ wavelet and total variation regularized reconstructions. The main advantage of the proposed preconditioning lies in the per-iteration computation time. k-space preconditioning is faster than circulant preconditioning, and adds very little computational overhead. This is expected as the circulant preconditioning requires two additional FFTs per iteration, whereas the proposed k-space diagonal preconditioning requires only element-wise multiplications. Moreover, in our experiments, ADMM always incurs additional computational overhead compared to other methods without inner loops. We conjecture that this may be due to additional memory copying and variable initialization in the inner loops. For example, residual terms in CG need to be initialized before each inner loop. On the other hand, we note that the computational disadvantage of inner loops can be reduced by pre-calculating the residual outside CG based on the solution found in the previous outer iteration.

Although the MC k-space preconditioner achieves the fastest convergence most of the time, the SC k-space preconditioner provides similar or sometimes even faster convergence acceleration. Since the MC k-space preconditioner cannot be precomputed, we recommend using the SC k-space preconditioner for general settings. In applications in which the same sensitivity maps and trajectory are used repeatedly, then the MC k-space preconditioner can provide a slightly faster convergence behavior with less oscillation. We have also investigated using the Pipe-Menon density compensation factor as a k-space preconditioner. While in some cases it accelerates convergence, Figures 4 and 5, and Supplementary Figures S6, S9, and S10 show that it can introduce excessive oscillations, which slow down the overall convergence. Based on this, we recommend the SC k-space preconditioner over density compensation factors as k-space preconditioners.

Finally, the experiment with the 3D UTE dataset in Figure 1 shows that the method offers orders of magnitude speedup for 3D datasets. This is expected because 3D trajectories have a higher variation in k-space density than 2D trajectories. In particular, the proposed method converged in about ten iterations, whereas other methods did not, even after hundreds of iterations. We note that most experiments in this article are done on 2D datasets because of their fast evaluation time. In terms of practical application, we expect the proposed preconditioning methods to be vastly more effective and useful for 3D non-Cartesian imaging.

IX Conclusion

We have shown a method to speed up non-Cartesian iterative reconstruction that retains the per-iteration computational efficiency of density compensation and reconstruction accuracy of preconditioning methods. In contrast to most existing preconditioning methods, the proposed technique does not increase the per-iteration computation time compared to vanilla iterative methods, such as the conjugate gradient method. With the proposed preconditioning, iterative reconstruction for non-Cartesian imaging can be accelerated without compromises.

Appendix A Derivation for the $\ell 2$ optimized diagonal preconditioner

We are interested in solving the following minimization problem:

[TABLE]

Expanding the objective function element-by-element, we obtain

[TABLE]

where $\bm{\delta}$ is the Dirac delta function.

Taking the gradient with respect to $\mathbf{p}_{i}$, setting it to zero and re-arranging, we have,

[TABLE]

Appendix B Derivation for the $\ell 2$ optimized Circulant Preconditioner

The circulant preconditioner we consider in this article minimizes the following problem:

[TABLE]

Each element of $\mathbf{A}^{H}\mathbf{A}$ is given by,

[TABLE]

and the circulant matrix $\mathbf{P}$ has the form,

[TABLE]

where $\mathbf{p}\in\mathbb{C}^{N}$ is the underlying convolution kernel, and $((~{}))_{N}$ denotes the modulo operation by $N$.

Hence, minimizing the least squares criterion results in,

[TABLE]

Let us define the autocorrelation function $\mathbf{r}\in\mathbb{C}^{2N-1}$ and the point-spread function $\mathbf{h}\in\mathbb{C}^{2N-1}$ as follows,

[TABLE]

then after some algebra, the convolution kernel can be expressed as,

[TABLE]

Derivation for the Convex Dual Problem

Here we will derive the dual problem through the Lagrangian function. Let us first introduce a variable $\mathbf{z}$ to make the objective function (3) a constrained optimization problem:

[TABLE]

Introducing a Lagrangian variable $\mathbf{u}$ gives us,

[TABLE]

Switching the min and the max, and minimizing over $\mathbf{z}$ gives us $\mathbf{z}=\mathbf{y}-\mathbf{u}$. Substituting and re-arranging, we obtain,

[TABLE]

Using the definition of the conjugate function $g^{*}(\mathbf{x}^{*})=\max_{\mathbf{x}^{*}}\langle\mathbf{x}^{*},\mathbf{x}\rangle-g(\mathbf{x})$, we have,

[TABLE]

Connection to Trzasko et al. [29]

Here we re-derive the result in Trzasko et al. through convex duality. Let us consider $g(\mathbf{x})=\frac{\lambda}{2}\|\mathbf{x}\|_{2}^{2}$, then the dual problem is given by,

[TABLE]

which has the optimality condition

[TABLE]

Hence, we can precondition in k-space by preconditioning the dual variable by solving:

[TABLE]

In general, the primal and dual variables $\mathbf{x}$ and $\mathbf{u}$ are connected with the following relationship:

[TABLE]

where $\partial g(\mathbf{x})$ denotes the sub-differential of $g$ at $\mathbf{x}$.

Since $\partial g(\mathbf{x})=\{\lambda\mathbf{x}\}$, from the primal dual relationship (12) we can recover the primal variable by performing,

[TABLE]

The above method is precisely what Trzasko et al. [29] proposed for the $\ell 2$-regularized sub-problem within ADMM.

Supporting Figures