Ultrafast Preconditioned Conjugate Gradient MAP reconstruction for fully 3-D microPET

Iterative list-mode PET reconstruction requires the pre-calculation of a sensitivity image, in theory by backprojecting of all possible detector lines of response (LORs). For systems with a large number of LORs, or because of subject motion compensation, computing the exact sensitivity image is computationally expensive, so sub-sampling the LORs is attractive to reduce the work. The resulting inaccuracies in the sensitivity image propagate to the reconstructed emission image, which motivates looking at different sub-sampling schemes to reduce the errors. We extended the MOLAR reconstruction, originally implemented using uniform random (Uniform) subsampling of LORs, by adding a quasi-random (Halton) sequence based sampler and also a previously described Monte Carlo-based (MonteCarlo) non-uniform sampler. We also evaluated sampling using both techniques simultaneously (Halton+MonteCarlo). A simple post-smoothing of the sensitivity image can be used to improve variability in exchange for a loss in sensitivity spatial resolution. We included this method (Uniform-Smoothed) in our comparison. Applying the five sampling methods (Uniform, Uniform-Smoothed, Halton, MonteCarlo, Halton+MonteCarlo), we obtain region of interest coefficients of variation (CoV) for replicated reconstructions of a human brain. Relative to Uniform, we observe two to three-fold improvement in CoV using Halton with a 4% increase in time to compute the sensitivity image, and slightly better improvements for Halton+MonteCarlo with a 163% increased computational investment. Halton and Halton+MonteCarlo produced lower CoVs than Uniform-Smoothed, without smoothing-induced bias in the reconstructed images.

We describe a fast forward and back projector pair based on inverse Fourier rebinning for use in iterative image reconstruction for fully 3-D positron emission tomography (PET). The projector pair is used as part of a factored system matrix that takes into account detector-pair response by using shift-variant sinogram blur kernels, thereby combining the computational advantages of Fourier rebinning with iterative reconstruction using accurate system models. The forward projector consists of a 2-D projector, which maps 3-D images into 2-D direct sinograms, followed by exact inverse rebinning which maps the 2-D into fully 3-D sinograms. The back projector is implemented as the transpose of the forward projector and differs from the true exact rebinning operator in the sense that it does not require reprojection to compute missing lines of response (LORs). We compensate for two types of inaccuracies that arise in a cylindrical PET scanner when using inverse Fourier rebinning: 1) nonuniform radial sampling and 2) nonconstant oblique angles in the radial direction in a single oblique sinogram. We examine the effects of these corrections on sinogram accuracy and reconstructed image quality. We evaluate performance of the new projector pair for maximum a posteriori (MAP) reconstruction of simulated and in vivo data. The new projector results in only a small loss in resolution towards the edge of the field-of-view when compared to the fully 3-D geometric projector and requires an order of magnitude less computation.

Gauss-Seidel preconditioned conjugate gradient (GS-PCG) method based on element-by-element (EBE) is derived and the calculation process of finite element method (FEM) is presented in this paper. The corresponding FEM program are developed and implemented on the magnetic field calculations which produced by straight conductor carrying current and 3-phases induction motor respectively. The validity and accuracy of the method have been verified by comparing the result with the analytical solution. Compared with Jacobi preconditioned conjugate gradient (J-PCG) method based on EBE-FEM, less computational time and cost are needed for the proposed method with same accuracy. The proposed method can be applied to electromagnetic field analysis of transformer and motor, which may refer large scale numerical computation.

A preconditioned Chebyshev basis communication-avoiding conjugate gradient method (P-CBCG) is applied to the pressure Poisson equation in a multiphase thermal-hydraulic CFD code JUPITER, and its computational performance and convergence properties are compared against a preconditioned conjugate gradient (P-CG) method and a preconditioned communication-avoiding conjugate gradient (P-CACG) method on the Oakforest-PACS, which consists of 8,208 KNLs. The P-CBCG method reduces the number of collective communications with keeping the robustness of convergence properties. Compared with the P-CACG method, an order of magnitude larger communication-avoiding steps are enabled by the improved robustness. It is shown that the P-CBCG method is 1.38 $$\times $$ and 1.17 $$\times $$ faster than the P-CG and P-CACG methods at 2,000 processors, respectively.

In this study, we analyze the real-time solution of a series of augmented least squares problems, which are generated by adding information to an original least squares model repetitively. Instead of solving the least squares problems directly, we transform them into a batch of saddle point linear systems and subsequently solve the linear systems using restrictively preconditioned conjugate gradient (RPCG) methods. Approximation of the new Schur complement is generated effectively based on a previously approximated Schur complement. Owing to the variations of the preconditioned conjugate gradient method, the proposed methods generate convergence results similar to the conjugate gradient method and achieve a very fast convergent iterative sequence when the coefficient matrix is well preconditioned. Numerical tests show that the new methods are more effective than some standard Krylov subspace methods. Updated RPCG methods meet the requirement of real-time computing successfully for multifactor models.

A novel sine transform based preconditioned conjugate gradient (PCG) method for signal restoration is presented. If the bandwidth 2? +1 in the model matrix is a constant, then the asymptotic complexity for the new algorithm is O(N) per iterative step, where Nis the order of the model. This complexity is reduced by an order of magnitude in comparison with the known fast Fourier transform (FFT) technique proposed recently in the literature. Moreover, if ? has a size similar to N, then the new PCG method is accomplished by the fast sine transform (FST) and the fast Hartley transform (FHT). The computational and storage cost per preconditioning step for the new PCG method is reduced by 50% as compared to the FFT approach. To stabilize the computation, a special Tikhonov regularization is introduced. Numerical experimentations show that the new PCG method is of a speedy convergence rate. The new PCG methods maintain less computational complexity per iterative step and have a better convergence rate than the other known PCG methods.

In order to solve Toeplitz linear systems An(f)x=b generated by a nonnegative integrable function f, through use of the preconditioned conjugate gradient (PCG) method, several authors have proposed An(g) as preconditioner in the case where g is a trigonometric polynomial [10, 14, 27, 12, 28]. In preceding works, we studied the distribution and the extremal properties of the spectrum of the preconditioned matrix G=An−1(g) An(f). In this paper we prove that the union of the spectra of all the Gn is dense on the essential range of f/g, i.e.,ER(f/g) and we obtain asymptotic information about the rate of convergence of the smallest eigenvalue λln of Gn to r (and of λnn to R). As a consequence of this second order result, it is possible to handle the case where f has zeros of any order θ, through the PCG methods proposed in [10, 14]. This is a noteworthy extension since the techniques developed in [10, 14, 27, 12, 28] are shown to be effective only when f has zeros of even orders. The cost of this procedure is O(n1+c(θ) log n) arithmetic operations (ops) where the quantity c(θ) belongs to interval [0,2−1] and takes the maximum value 2−1 when f has a zero of odd order. Finally, for the special case of zeros of odd orders, we propose a further algorithm which makes use of the PCG techniques proposed in [10, 14, 27, 12, 28] for theeven order case, reducing the cost to O(n long n) ops.

High-resolution global ocean general circulation models (OGCMs) play a key role in accurate ocean forecasting. However, the models of the operational forecasting systems are still not in high resolution due to the subsequent high demand for large computation, as well as the low parallel efficiency barrier. Good scalability is an important index of parallel efficiency and is still a challenge for OGCMs. We found that the communication cost in a barotropic solver, namely, the preconditioned conjugate gradient (PCG) method, is the key bottleneck for scalability due to the high frequency of the global reductions. In this work, we developed a new algorithm—a communication-avoiding Krylov subspace method with a PCG (CA-PCG)—to improve scalability and then applied it to the Nucleus for European Modelling of the Ocean (NEMO) as an example. For PCG, inner product operations with global communication were needed in every iteration, while for CA-PCG, inner product operations were only needed every eight iterations. Therefore, the global communication cost decreased from more than 94.5% of the total execution time with PCG to less than 63.4% with CA-PCG. As a result, the execution time of the barotropic modes decreased from more than 17,000 s with PCG to less than 6000 s with CA-PCG, and the total execution time decreased from more than 18,000 s with PCG to less than 6200 s with CA-PCG. Besides, the ratio of the speedup can also be increased from 3.7 to 4.6. In summary, the high process count scalability when using CA-PCG was effectively improved from that using the PCG method, providing a highly effective solution for accurate ocean simulation.

For transonic flow calculations, a combination of the successive line over-relaxation (SLOR) and the preconditioned conjugate gradient (CG) method has been suggested by Wong and Hafez (1981). This paper studies the method of minimal residual (MR) which avoids a combined iteration. This method is closely related to the CG method, may be regarded as a first-order gradient method, and is applicable to symmetric and nonsymmetric matrices. The problem is formulated as a nonlinear mixed elliptic-hyperbolic partial differential equation which includes an artificial viscosity and a switching function which is zero in subsonic regions and nonzero in supersonic regions. Alternatives to the SLOR method which provide faster convergence rates are introduced. The preconditioned MR algorithm is developed, and transonic potential flows around NACA 0012 airfoil are calculated for different Mach numbers and angles of attack. Preliminary results are presented, demonstrating that the MR algorithm requires no parameter estimation and rapidly converges for subsonic flows.

Classical iterative methods when applied to the partial solution of the generalized eigenvalue problem Ax =λBx, may yield very poor convergence rates particularly when ill‐conditioned problems are considered. In this paper the preconditioned conjugate gradient (CG) method via the minimization of the Rayleigh quotient and the reverse power method is employed for the partial eigenproblem. The triangular splitting preconditioners employed are obtained from an incomplete Choleski factorization and a partial Evans preconditioner. This approach can dramatically improve the convergence rate of the basic CG method and is applicable to any symmetric eigenproblem in which one of the matrices A, B is positive definite. Because of the renewed interest in CG techniques for FE work on microprocessors and parallel computers, it is believed that this improved approach to the generalized eigenvalue problem is likely to be very promising.

In positron emission tomography (PET), 3D iterative image reconstruction methods have a huge computational burden. In this paper, we developed a list-mode image reconstruction method using graphics processing units (GPUs). Efficiency of acceleration for GPU implementation largely depends on the method chosen, where a reduced number of conditional statements and a reduced memory size are required. On the other hand, accurate system models are required to improve the quality of reconstructed images. Various accurate system models for conventional CPU implementation have been proposed, but these models basically require many conditional statements and huge memory size. Therefore, we developed a new system model which matches GPU implementation better. In this model, the detector response functions, which vary depending on each line of response (LOR), are pre-computed in CPUs and modeled by sixth-order polynomial functions in order to reduce the memory size occupied in GPUs. Each element of a system matrix is obtained on-the-fly in GPUs by calculating the distance between an LOR and a voxel. Therefore the developed system model enables efficient GPU implementation of the accurate system modeling with a reduced number of conditional statements and a reduced memory size. We applied the developed method to a small OpenPET prototype, in which 4-layered depth-of-interaction (DOI) detectors were used. For image reconstruction, we used the dynamic row-action maximum likelihood algorithm (DRAMA). Compared with a conventional model for GPU implementation, in which DRFs are given as a Gaussian function of fixed width, we saw no remarkable difference for DOI data, but for non-DOI data the proposed model outperformed the conventional at the peripheral region of the field-of-view. The proposed model had almost the same calculation time as the conventional model did. For further acceleration, we tried parallel GPU implementation, and we obtained 3.8-fold acceleration by using 4 GPUs.

We have developed a new method to transform any kind of 3D PET scanner data to a regular cylinder structure which consists of virtual crystals arranged in regular rings without any gaps. In this way, the 3D projection data can be stored without introducing geometrical mismatches with respect to the assigned projectors of the iterative reconstruction. In addition, the physical crystal sizes of the scanner can be taken into account during the data conversion procedure. The conversion of all measured Lines-of-Response (LOR) to a virtual, regular cylinder geometry allows to formalise the geometrical transcription while simultaneously minimising the loss of information, e.g. arc corrections can be omitted completely. Finally, the regularity of the converted projection data allows a very fast and accurate implementation of the reconstruction algorithm which is based on a memory-resident compressed system response matrix. The geometrical accuracy of the applied projectors depends only on the user-defined granularity of the virtual geometry which can be adjusted as a compromise between operating expense and precision. The ClearPET Neuro Scanner [2] with its rotating gantry and the crystal double-layer design has been used for a first application of the new universal 3D reconstruction technique. The presented results indicate the high potential of the new conversion method in terms of achievable contrast-noise ratio and resolution.

Groundwater flow simulation has become one of the top international issues in new generation of environmental applications. When managing large-scale groundwater flow problems, the intensive computational ability and large amounts of memory space required for modeling are the main bottlenecks for researchers. In order to solve three-dimensional large-scale groundwater flow problems more rapidly, the Open MP was adopted to parallelize the preconditioned conjugate gradient (PCG) algorithm in this paper. And this paper carried out a numerical experiment of the three-dimensional groundwater flow model on a computer with four cores. Based on the numerical experiment, it is found that the execution time of the original serial PCG program is about 1.74 to 2.86 times of the parallel PCG program executed with different number of threads. The experimental results also demonstrate that the PCG solver based on Open MP is an effective way for solving large-scale three-dimensional groundwater flow problem.

Utility of the preconditioned conjugate gradient algorithm with a diagonal preconditioner for solving mixed-model equations in animal breeding applications was evaluated with 16 test problems. The problems included single- and multiple-trait analyses, with data on beef, dairy, and swine ranging from small examples to national data sets. Multiple-trait models considered low and high genetic correlations. Convergence was based on relative differences between left- and right-hand sides. The ordering of equations was fixed effects followed by random effects, with no special ordering within random effects. The preconditioned conjugate gradient program implemented with double precision converged for all models. However, when implemented in single precision, the preconditioned conjugate gradient algorithm did not converge for seven large models. The preconditioned conjugate gradient and successive overrelaxation algorithms were subsequently compared for 13 of the test problems. The preconditioned conjugate gradient algorithm was easy to implement with the iteration on data for general models. However, successive overrelaxation requires specific programming for each set of models. On average, the preconditioned conjugate gradient algorithm converged in three times fewer rounds of iteration than successive overrelaxation. With straightforward implementations, programs using the preconditioned conjugate gradient algorithm may be two or more times faster than those using successive overrelaxation. However, programs using the preconditioned conjugate gradient algorithm would use more memory than would comparable implementations using successive overrelaxation. Extensive optimization of either algorithm can influence rankings. The preconditioned conjugate gradient implemented with iteration on data, a diagonal preconditioner, and in double precision may be the algorithm of choice for solving mixed-model equations when sufficient memory is available and ease of implementation is essential.

<abstract><p>The conjugate gradient (CG) method is an optimization method, which, in its application, has a fast convergence. Until now, many CG methods have been developed to improve computational performance and have been applied to real-world problems. In this paper, a new hybrid three-term CG method is proposed for solving unconstrained optimization problems. The search direction is a three-term hybrid form of the Hestenes-Stiefel (HS) and the Polak-Ribiére-Polyak (PRP) CG coefficients, and it satisfies the sufficient descent condition. In addition, the global convergence properties of the proposed method will also be proved under the weak Wolfe line search. By using several test functions, numerical results show that the proposed method is most efficient compared to some of the existing methods. In addition, the proposed method is used in practical application problems for image restoration and portfolio selection.</p></abstract>