To accelerate iterative, non-Cartesian image reconstruction with a circulant majorizer designed based on the k-space trajectory.

MRI scans can be accelerated by using non-Cartesian trajectories. The scans can be further accelerated by subsampling k-space and reconstructing images by solving an optimization problem of the form,¹

$$\hat{x}=\underset{x}{\text{argmin}} \frac{1}{2} \left \| y-Ax\right \|_2^2+\beta \left \|Rx \right \|_1,$$

where $$$A$$$ is the system matrix incorporating the non-Cartesian trajectory, $$$x$$$ is the object being reconstructed, $$$R$$$ is a sparsity-promoting transform (e.g., undecimated Haar wavelets), and $$$\beta$$$ is a regularization parameter. Majorize-minimize algorithms require upper bounding $$$A'A$$$ of the data fit term with a majorizing matrix, $$$M_f$$$ (i.e., $$$M_f \succeq A'A$$$ where $$$N \succeq 0$$$ implies that $$$N$$$ is positive semidefinite). The standard approach is to use $$$M_f=LI$$$, where $$$L$$$ is the maximum eigenvalue of $$$A'A$$$.² This is a loose bound for $$$A'A$$$ due to the high density of samples near the center of k-space. Previous methods accelerated convergence by running power iterations to get step sizes for each frequency subband of the wavelet transform,³ but these methods used orthogonal wavelets and a moving target cost function with weaker convergence guarantees.

We derive a generalization of the previous work (BARISTA/FISTA^2,4) for non-Cartesian sampling by calculating a circulant majorizing matrix, $$$M_f$$$. We assume that we are working in the SENSE setting (i.e., $$$A=FS$$$ where $$$F$$$ is the non-Cartesian Fourier encoding operator and $$$S$$$ is a matrix with sensitivity maps) where $$$R$$$ is an undecimated wavelet transform. We first majorize $$$F'F$$$ with a circulant matrix, $$$M_f$$$, by computing $$$D_f=\text{diag}(QFe_0)$$$, where $$$Q$$$ is a DFT matrix and $$$e_0$$$ is an impulse. We then increase the minimum entry in $$$D_f$$$ until $$$M_f=Q^{-1}D_fQ \succeq F'F$$$ as determined via power iteration.⁵ We found that this was also a majorizing matrix for $$$A'A$$$ in numerical experiments when the sensitivity maps are calculated with a square root sum of squares normalization. Then, we design a diagonal matrix, $$$D_R$$$, such that $$$D_R \succeq RM_f^{-1}R^T$$$. $$$D_R$$$ governs the algorithm step sizes in the regularizer domain. Note that $$$R^T=[R_1^T ... R_B^T]$$$ for an undecimated wavelet transform with $$$B$$$ subbands. For each subband, we have $$$R_b=Q^{-1}\Lambda_b Q$$$, where $$$\Lambda_b$$$ is the frequency response for the filter in the $$$b$$$th subband. Inserting this into $$$RM_f^{-1}R^T$$$ reveals that we can apply Gershgorin's Theorem to build $$$D_R=\text{diag}(d_{R,b}I_b)$$$, where $$$I_b$$$ is an identity matrix of the size corresponding to the $$$b$$$th subband and $$$d_{R,b}=\text{max}_n(\sum_{p=1}^B |\lambda_{n,b} \lambda_{n,p}^*/d_{n,f}|)$$$, where $$$\lambda_{n,b}$$$ is the $$$n$$$th entry in $$$\Lambda_b$$$ and $$$d_{n,f}$$$ is the $$$n$$$th entry in $$$D_f$$$. Within the BARISTA⁴ analysis algorithm (a FISTA variant), this gives the analysis denoising steps of

$$q^{(j+1)}=P_{P^M}(v^{(j)} - \beta^{-1} D_R^{-1}Rx^{(k,j+1)}),$$

$$x^{(k,j+1)}=c^{(k)}-\beta M_f^{-1}R^Tv^{(j)},$$

$$c^{(k)}=x^{(k)} - M_f^{-1}A'(Ax^{(k)}-y).$$

We compared the speed of this new algorithm to that of FISTA with the Lipschitz constant in numerical experiments. The experiments simulated an interleaved spiral acquisition with 48 leaves for a 128 by 128 image and 8 receive coils. 12 of the 48 leaves were used for reconstruction. The regularization parameter was chosen manually.

We calculated a "converged" solution, $$$x^{(\infty)}$$$, by running many thousands of iterations of FISTA. Figure 1 shows the converged image ($$$x^{(\infty)}$$$), while Figures 2 and 3 show $$$x^{(60)}$$$ of FISTA and the proposed method, respectively. The proposed method has more substantially reduced subsampling artifacts compared to FISTA after 60 iterations. We plot $$$\xi(k) = \frac{\left \| x^{(k)} - x^{(\infty)} \right \|}{\left \| x^{(\infty)} \right \|}$$$, the norm residual to convergence vs. time in Figure 4 to compare the convergence speed of the algorithms. The proposed method has about a 2-3 factor increase in convergence speed over FISTA.

We have shown accelerated convergence by upper bounding $$$A'A$$$ in the non-Cartesian, SENSE parallel MRI setting a circulant majorizer. The circulant majorizer has a frequency response that is related to the density compensation function of the k-space trajectory. We further showed that this property can be exploited in the setting where the regularizing matrix, $$$R$$$, is an undecimated Haar wavelet transform. This gives different step sizes for each frequency subband of the Haar wavelet transform, and the algorithm can be extended to any case where $$$R$$$ can be decomposed into $$$B$$$ filters. In the future we plan to examine strategies for creating smaller $$$M_f$$$ matrices that more tightly bound $$$A'A$$$.