## 4 Results  
To compare the performance of the three implementations, we measured the elapsed time per conjugate gradient (CG) iteration across the six datasets listed in Table 1. Each reported time represents the average of five runs. The results are shown in Figure 2, where the primary y-axis (left) indicates the elapsed time (in seconds) per CG iteration, and the secondary y-axis (right) denotes the speedup relative to the CPU reference implementation. The CUDA-based implementation significantly reduced the elapsed time, achieving speedup factors of between 9.7 and 13.9 (green dashed line) compared to the CPU-only version. The Matlab-GPU implementation (blue dashed line) also improved speedup but was less pronounced than the CUDA solution. The NRMSE of the CPU and the CUDA images are depicted in Table 4. They remained within the range of \(10^{-6}\) to \(10^{-4}\) for all datasets. Since these values are lower than the reconstruction kernel's tolerance of \(\epsilon_{r}=10^{-3}\), no detectable image quality degradation was observed. Figure 3 presents a roofline analysis of the \(E^{H}E\) operator, highlighting both the achievable performance and operational intensity for each operator in Algorithm 1. We relied on the Empirical Roofline Toolkit to measure realistic bandwidth and peak FP32/FP64 performance values rather than theoretical specifications. By computing the motion operator elements on the fly, the performance of operator \(T_{t}\) was significantly enhanced, reaching 105 GFLOP/s. The Fast Fourier Transform (FFT) steps, executed via the cuFFT library, achieved 310 GFLOP/s. Overall, the \(E^{H}E\) function is bandwidth-bound, exhibiting a performance of approximately 250 GFLOP/s at a high-bandwidth memory (HBM) operational intensity of 2.94 FLOP/Byte. For comparison, ref. [13] reported a performance of 16 GFLOP/s for the Matlab GPU implementation of Algorithm 1. In contrast, our native GPU implementation delivered a 15.6\(\times\) performance gain, narrowing the gap to the roofline limits that define the upper bounds of achievable performance. Figures 4 and 5 illustrate the reconstructed images for each implementation. In both the brain MRI dataset (\(D_{4}\)), subject to rigid-body motion, and the breast MRI dataset (\(D_{2}\)), subject to non-rigid respiratory and cardiac motion, the three methods effectively suppressed motion-induced artifacts.  
Figure 2: Comparison of the elapsed time per CG iteration for the three implementations.  
Figure 4: Comparison of GRICS reconstructions with the three implementations for dataset \(D_{2}\). (**A**) Uncorrected reconstruction. (**B**): The CPU-only solution. (**C**): The Matlab-GPU solution. (**D**): The CUDA solution. The three solutions markedly reduced the motion artifacts visible in the blue and orange squares.  
\begin{table}
\begin{tabular}{c c} \hline
**Dataset** & **NRMSE (CPU_Image, CUDA_Image)** \\ \hline \(D_{1}\) & \(10^{-5}\) \\ \hline \(D_{2}\) & \(10^{-6}\) \\ \hline \(D_{3}\) & \(10^{-5}\) \\ \hline \(D_{4}\) & \(10^{-4}\) \\ \hline \(D_{5}\) & \(10^{-4}\) \\ \hline \(D_{6}\) & \(10^{-4}\) \\ \hline \end{tabular}
\end{table}
Table 4: NRMSE of the CPU image and the CUDA image.  
Figure 3: Roofline analysis of the implemented CUDA kernels on NVIDIA A100. The work in [24] was used to plot the model. Dataset \(D_{1}\).