Preconditioned Conjugate Gradient Algorithm Research Articles

Regularized Least Squares Estimating Sensitivity for Self-calibrating Parallel Imaging

Calibration of the spatial sensitivity functions of coil arrays is a crucial element in parallel magnetic resonance imaging (pMRI). The self-calibrating technique for sensitivity extraction has complemented the common calibration technique that uses a separate pre-scan. In order to improve the accuracy of sensitivity estimate from small number of self-calibrating data, which is extracted from a fully sampled central region of a variable-density k-space acquisition in self-calibrating parallel images, a novel scheme for estimating the sensitivity profiles is proposed in the paper. On consideration of truncation error and measurement errors in self-calibrating data, the issue of calculating sensitivity would be formulated as a regularized least squares estimation problem, which is solved by the preconditioned conjugate gradients algorithm. When applying the estimated coil sensitivity to reconstruct full field-of-view(FOV) image from the under-sampling simulated and in vivo data, the normalized signal-to-noise ratio (NSNR) of reconstruction image is evidently improved, and meanwhile the normalized mean squared error (NMSE) is remarkably reduced, especially when a rather large accelerate factor is used.

Preconditioned non-linear conjugate gradient method for frequency domain full-waveform seismic inversion

We present preconditioned non-linear conjugate gradient algorithms as alternatives to the Gauss-Newton method for frequency domain full-waveform seismic inversion. We designed two preconditioning operators. For the first preconditioner, we introduce the inverse of an approximate sparse Hessian matrix. The approximate Hessian matrix, which is highly sparse, is constructed by judiciously truncating the Gauss-Newton Hessian matrix based on examining the auto-correlation and cross-correlation of the Jacobian matrix. As the second preconditioner, we employ the approximation of the inverse of the Gauss-Newton Hessian matrix. This preconditioner is constructed by terminating the iteration process of the conjugate gradient least-squares method, which is used for inverting the Hessian matrix before it converges. In our preconditioned non-linear conjugate gradient algorithms, the step-length along the search direction, which is a crucial factor for the convergence, is carefully chosen to maximize the reduction of the cost function after each iteration. The numerical simulation results show that by including a very limited number of non-zero elements in the approximate Hessian, the first preconditioned non-linear conjugate gradient algorithm is able to yield comparable inversion results to the Gauss-Newton method while maintaining the efficiency of the un-preconditioned non-linear conjugate gradient method. The only extra cost is the computation of the inverse of the approximate sparse Hessian matrix, which is less expensive than the computation of a forward simulation of one source at one frequency of operation. The second preconditioned non-linear conjugate gradient algorithm also significantly saves the computational expense in comparison with the Gauss-Newton method while maintaining the Gauss-Newton reconstruction quality. However, this second preconditioned non-linear conjugate gradient algorithm is more expensive than the first one.

FSAI preconditioned CG algorithm combined with GPU technique for the finite element analysis of electromagnetic scattering problems

In order to efficiently solve the large sparse complex linear system arising from the vector finite element method (vector FEM) in electromagnetic scattering problems, the factorized sparse approximate inverse (FSAI) algorithm and the programmable graphics processing unit (GPU) are employed in the context of the conjugate gradient (CG) iterative method. The combination of the FSAI with the GPU technique has two advantages. Firstly, the convergence rate of the CG algorithm is significantly accelerated. Secondly, the calculation of the sparse matrix vector product (SMVP) in the FSAI preconditioned CG algorithm is accelerated by harnessing the tremendous parallel processing capacity of the GPU. Numerical experiments indicate that the FSAI preconditioned CG algorithm enhanced by the GPU technique is very effective and can reduce both the number of iterations and the computational time significantly.

A Novel Preconditioned Conjugate Gradient Image Reconstruction Algorithm for Electrical Capacitance Tomography System

To solve the‘soft-field’nature and the ill-posed problem in electrical capacitance tomography (ECT)technology, a novel preconditioned conjugate gradient algorithm for electrical capacitance tomography is presented. The analysis of the basic principles of electrical capacitance tomography, based on the given preconditioned conjugate gradient method and calculation formula of the iteration steps of the ECT and to explore the feasibility of application of the algorithm, the algorithm to meet the convergence conditions and the reconstruction image error. Experimental results and simulation data indicate that the algorithm can provide high quality images and favorable stabilization compared with LBP，conjugate gradient algorithms in simple flow pattern and this new algorithm presents a feasible and effective way to research on image reconstruction algorithm for Electrical Capacitance Tomography System.

Fast Bayesian Estimation for VARFIMA Processes with Stable Errors

We present an approach via a multivariate preconditioned conjugate gradient (MPCG) algorithm for Bayesian inference for vector ARFIMA models with sub-Gaussian stable errors. This approach involves solution of a block-Toeplitz system, and treating the unobserved process history and the underlying positive stable process as unknown parameters in the joint posterior. We use Gibbs sampling with the Metropolis-Hastings algorithm. We illustrate our approach on time series of daily average temperatures measured over several years at different U.S. cities.

Updating incomplete factorization preconditioners for shifted linear systems arising in a wind model

The efficiency of a finite element mass consistent model for wind field adjustment depends on the stability parameter α which allows adjustment from a strictly horizontal wind to a pure vertical one. Each simulation with the wind model leads to the resolution of a linear system of equations, the matrix of which depends on a function ε ( α ) , i.e., ( M + ε N ) x ε = b ε , where M and N are constant, symmetric and positive definite matrices with the same sparsity pattern for a given level of discretization. The estimation of this parameter may be carried out by using genetic algorithms. This procedure requires the evaluation of a fitness function for each individual of the population defined in the searching space of α , that is, the resolution of one linear system of equations for each value of α . Preconditioned Conjugate Gradient algorithm (PCG) is usually applied for the resolution of these types of linear systems due to its good convergence results. In order to solve this set of linear systems, we could either construct a different preconditioner for each of them or use a single preconditioner constructed from the first value of ε to solve all the systems. In this paper, an intermediate approach is proposed. An incomplete Cholesky factorization of matrix A ε is constructed for the first linear system and it is updated for each ε at a low computational cost. Numerical experiments related to realistic wind field are presented in order to show the performance of the proposed preconditioning strategy.

Accelerated scaled memoryless BFGS preconditioned conjugate gradient algorithm for unconstrained optimization

An accelerated scaled memoryless BFGS preconditioned conjugate gradient algorithm for solving unconstrained optimization problems is presented. The basic idea is to combine the scaled memoryless BFGS method and the preconditioning technique in the frame of the conjugate gradient method. The preconditioner, which is also a scaled memoryless BFGS matrix, is reset when the Beale–Powell restart criterion holds. The parameter scaling the gradient is selected as a spectral gradient. For the steplength computation the method has the advantage that in conjugate gradient algorithms the step lengths may differ from 1 by two order of magnitude and tend to vary unpredictably. Thus, we suggest an acceleration scheme able to improve the efficiency of the algorithm. Under common assumptions, the method is proved to be globally convergent. It is shown that for uniformly convex functions the convergence of the accelerated algorithm is still linear, but the reduction in the function values is significantly improved. In mild conditions the algorithm is globally convergent for strongly convex functions. Computational results for a set consisting of 750 unconstrained optimization test problems show that this new accelerated scaled conjugate gradient algorithm substantially outperforms known conjugate gradient methods: SCALCG [3–6], CONMIN by Shanno and Phua (1976, 1978) [42,43], Hestenes and Stiefel (1952) [25], Polak–Ribiére–Polyak (1969) [32,33], Dai and Yuan (2001) [17], Dai and Liao (2001) ( t = 1 ) [14], conjugate gradient with sufficient descent condition [7], hybrid Dai and Yuan (2001) [17], hybrid Dai and Yuan zero (2001) [17], CG_DESCENT by Hager and Zhang (2005, 2006) [22,23], as well as quasi-Newton LBFGS method [26] and truncated Newton method by Nash (1985) [27].

AN ADDITIVE POISSON DATA MODELING OF ATTENUATION MAP RECONSTRUCTION FROM POST-INJECTION POSITRON EMISSION TOMOGRAPHY TRANSMISSION SCAN

We present an additive Poisson data modeling of attenuation map reconstruction from post-injection positron emission tomography transmission scan. The transmission scan data are modeled as a collection of independent Poisson observations which consist of true transmission, cross-contamination emission and random coincidences. The mean values of cross-contamination emission and random data are estimated by using Luk's approximation and Zaers' method, respectively. The unknown attenuation map is also modeled using a Gibbs image prior. Under the paradigm of Bayesian statistics, a pre-conditioned conjugate gradient algorithm is then used to compute a maximum a posteriori estimate of the attenuation map by maximizing over the posterior density. Evaluations are conducted on clinical data collected from a CTI EXACT HR+ scanner. Experimental results have indicated that the reconstructed results of attenuation map by the proposed reconstruction algorithm match well with theoretical values of soft tissue, bone and lung under 511 keV photons.

Computing options for multiple‐trait test‐day random regression models while accounting for heat tolerance

Data included 90,242,799 test day records from first, second and third parities of 5,402,484 Holstein cows and 9,326,754 animals in the pedigree. Additionally, daily temperature humidity indexes (THI) from 202 weather stations were available. The fixed effects included herd test day, age at calving, milking frequency and days in milk classes (DIM). Random effects were additive genetic, permanent environment and herd-year and were fit as random regressions. Covariates included linear splines with four knots at 5, 50, 200 and 305 DIM and a function of THI. Mixed model equations were solved using an iteration on data program with a preconditioned conjugate gradient algorithm. Preconditioners used were diagonal (D), block diagonal due to traits (BT) and block diagonal due to traits and correlated effects (BTCORR). One run included BT with a 'diagonalized' model in which the random effects were reparameterized for diagonal (co)variance matrices among traits (BTDIAG). Memory requirements were 8.7 Gb for D, 10.4 Gb for BT and BTDIAG, and 24.3 Gb for BTCORR. Computing times (rounds) were 14 days (952) for D, 10.7 days (706) for BT, 7.7 days (494) for BTDIAG and 4.6 days (289) for BTCORR. The convergence pattern was strongly influenced by the choice of fixed effects. When sufficient memory is available, the option BTCORR is the fastest and simplest to implement; the next efficient method, BTDIAG, requires additional steps for diagonalization and back-diagonalization.

Short communication: Genetic trends of milk yield under heat stress for US Holsteins

Data included 90,242,799 test-day milk records from 5,402,484 Holstein cows in the first 3 parities and 9,326,754 animals in the pedigree. Additionally, daily temperature-humidity indexes from 202 weather stations were available. Analyses were done by a random regression model in which each parity was treated as a separate trait and that accounted for heat stress. The fixed effects included herd test-day, age at calving, milking frequency, and days in milk classes. Random effects included additive genetic, permanent environment, and herd-year effects, all fit as random regressions. Five covariates in the random regressions included linear splines with 4 knots at 5, 50, 200, and 305 DIM and a function of a temperature-humidity index (THI). Mixed model equations were solved by using an iteration on data approach with a preconditioned conjugate gradient algorithm. Genetic trends for daily milk yield in absence of heat stress (intercept) were 0.140 kg/yr, 0.172 kg/yr, and 0.168 kg/yr for the first, second, and third parity, respectively. Genetic trends for decline of milk yield at temperature of 5°C THI over the threshold of sensitivity to heat stress were −0.002 kg/yr, −0.035 kg/yr, and −0.038 kg/yr, for first, second, and third parity, respectively. Genetic profiles were created by contrasting the 100 most and 100 least heat-tolerant bulls for the official proofs. The most heat-tolerant bulls transmitted lower production and dairy form but higher fertility, productive life, and type, especially udder and locomotion traits. In later parities, the type advantages were smaller. Test-day records capture only a fraction of information due to heat stress, and the real trends for heat stress may be stronger. Studies on heat stress for production should include records on later parities.

Preconditioned conjugate gradient algorithms for nonconvex problems with box constraints

The paper describes new conjugate gradient algorithms for large scale nonconvex problems with box constraints. In order to speed up convergence the algorithms employ scaling matrices which transform the space of original variables into the space in which Hessian matrices of the problem’s functionals have more clustered eigenvalues. This is done by applying limited memory BFGS updating matrices. Once the scaling matrix is calculated, the next few conjugate gradient iterations are performed in the transformed space. The box constraints are treated efficiently by the projection. We also present a limited memory quasi-Newton method which is a special version of our general algorithm. The presented algorithms have strong global convergence properties, in particular they identify constraints active at a solution in a finite number of iterations. We believe that they are competitive to the L-BFGS-B method and present some numerical results which support our claim.

A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations

We discuss the scalable parallel solution of the Poisson equation within a Particle-In-Cell (PIC) code for the simulation of electron beams in particle accelerators of irregular shape. The problem is discretized by Finite Differences. Depending on the treatment of the Dirichlet boundary the resulting system of equations is symmetric or ‘mildly’ nonsymmetric positive definite. In all cases, the system is solved by the preconditioned conjugate gradient algorithm with smoothed aggregation (SA) based algebraic multigrid (AMG) preconditioning. We investigate variants of the implementation of SA-AMG that lead to considerable improvements in the execution times. We demonstrate good scalability of the solver on distributed memory parallel processor with up to 2048 processors. We also compare our iterative solver with an FFT-based solver that is more commonly used for applications in beam dynamics.

A parallel multigrid preconditioner for the simulation of large fracture networks

Computational modeling of a fracture in disordered materials using discrete lattice modelsrequires the solution of a linear system of equations every time a new lattice bond isbroken. Solving these linear systems of equations successively is the most expensive part offracture simulations using large three-dimensional networks. In this paper, wepresent a parallel multigrid preconditioned conjugate gradient algorithm to solvethese linear systems. Numerical experiments demonstrate that this algorithmperforms significantly better than the algorithms previously used to solve thisproblem.

Acceleration of the direct reconstruction of linear parametric images using nested algorithms

Parametric imaging using dynamic positron emission tomography (PET) provides important information for biological research and clinical diagnosis. Indirect and direct methods have been developed for reconstructing linear parametric images from dynamic PET data. Indirect methods are relatively simple and easy to implement because the image reconstruction and kinetic modeling are performed in two separate steps. Direct methods estimate parametric images directly from raw PET data and are statistically more efficient. However, the convergence rate of direct algorithms can be slow due to the coupling between the reconstruction and kinetic modeling. Here we present two fast gradient-type algorithms for direct reconstruction of linear parametric images. The new algorithms decouple the reconstruction and linear parametric modeling at each iteration by employing the principle of optimization transfer. Convergence speed is accelerated by running more sub-iterations of linear parametric estimation because the computation cost of the linear parametric modeling is much less than that of the image reconstruction. Computer simulation studies demonstrated that the new algorithms converge much faster than the traditional expectation maximization (EM) and the preconditioned conjugate gradient algorithms for dynamic PET.

A separable nonlinear model for the multicommodity flow problem

This article studies the convex nonlinear multicommodity network flow problem with separable increasing cost. A numerical implementation of the primal-dual interior-point method is designed to solve the problem. At each iteration of the interior-point method, the corresponding linear problem is solved, if it is expressed by the normal equations, by using the approximate inverse algorithm (AINV) alone or combined with the preconditioned conjugate gradient algorithm. If the linear system is expressed as an augmented indefinite system, it is solved by using the AINV algorithm combined with an indefinite preconditioned conjugate gradient algorithm. Numerical experiments are conducted for networks of different dimensions and numbers of commodities and for some nonlinear costs. The computational results show the effectiveness of the interior-point method for this class of network problems.

Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion

Quantitative imaging of the elastic properties of the subsurface at depth is essential for civil engineering applications and oil- and gas-reservoir characterization. A realistic synthetic example provides for an assessment of the potential and limits of 2D elastic full-waveform inversion (FWI) of wide-aperture seismic data for recovering high-resolution P- and S-wave velocity models of complex onshore structures. FWI of land data is challenging because of the increased nonlinearity introduced by free-surface effects such as the propagation of surface waves in the heterogeneous near-surface. Moreover, the short wavelengths of the shear wavefield require an accurate S-wave velocity starting model if low frequencies are unavailable in the data. We evaluated different multiscale strategies with the aim of mitigating the nonlinearities. Massively parallel full-waveform inversion was implemented in the frequency domain. The numerical optimization relies on a limited-memory quasi-Newton algorithm thatoutperforms the more classic preconditioned conjugate-gradient algorithm. The forward problem is based upon a discontinuous Galerkin (DG) method on triangular mesh, which allows accurate modeling of free-surface effects. Sequential inversions of increasing frequencies define the most natural level of hierarchy in multiscale imaging. In the case of land data involving surface waves, the regularization introduced by hierarchical frequency inversions is not enough for adequate convergence of the inversion. A second level of hierarchy implemented with complex-valued frequencies is necessary and provides convergence of the inversion toward acceptable P- and S-wave velocity models. Among the possible strategies for sampling frequencies in the inversion, successive inversions of slightly overlapping frequency groups is the most reliable when compared to the more standard sequential inversion of single frequencies. This suggests that simultaneous inversion of multiple frequencies is critical when considering complex wave phenomena.

Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods

SUMMARY We present a full seismic waveform tomography for upper-mantle structure in the Australasian region. Our method is based on spectral-element simulations of seismic wave propagation in 3-D heterogeneous earth models. The accurate solution of the forward problem ensures that waveform misfits are solely due to as yet undiscovered Earth structure and imprecise source descriptions, thus leading to more realistic tomographic images and source parameter estimates. To reduce the computational costs, we implement a long-wavelength equivalent crustal model. We quantify differences between the observed and the synthetic waveforms using time–frequency (TF) misfits. Their principal advantages are the separation of phase and amplitude misfits, the exploitation of complete waveform information and a quasi-linear relation to 3-D Earth structure. Fr´ echet kernels for the TF misfits are computed via the adjoint method. We propose a simple data compression scheme and an accuracy-adaptive time integration of the wavefields that allows us to reduce the storage requirements of the adjoint method by almost two orders of magnitude. To minimize the waveform phase misfit, we implement a pre-conditioned conjugate gradient algorithm. Amplitude information is incorporated indirectly by a restricted line search. This ensures that the cumulative envelope misfit does not increase during the inversion. An efficient pre-conditioner is found empirically through numerical experiments. It prevents the concentration of structural heterogeneity near the sources and receivers. We apply our waveform tomographic method to ≈1000 high-quality vertical-component seismograms, recorded in the Australasian region between 1993 and 2008. The waveforms comprise fundamental- and higher-mode surface and long-period S body waves in the period range from 50 to 200 s. To improve the convergence of the algorithm, we implement a 3-D initial model that contains the long-wavelength features of the Australasian region. Resolution tests indicate that our algorithm converges after around 10 iterations and that both long- and short-wavelength features in the uppermost mantle are well resolved. There is evidence for effects related to the non-linearity in the inversion procedure. After 11 iterations we fit the data waveforms acceptably well; with no significant further improvements to be expected. During the inversion the total fitted seismogram length increases by 46 per cent, providing a clear indication of the efficiency and consistency of the iterative optimization algorithm. The resulting SV -wave velocity model reveals structural features of the Australasian upper mantle with great detail. We confirm the existence of a pronounced low-velocity band along the eastern margin of the continent that can be clearly distinguished against Precambrian Australia and the microcontinental Lord Howe Rise. The transition from Precambrian to Phanerozoic Australia (the Tasman Line) appears to be sharp down to at least 200 km depth. It mostly occurs further east of where it is inferred from gravity and magnetic anomalies. Also clearly visible are the Archean and Proterozoic cratons, the northward continuation of the continent and anomalously low S-wave velocities in the upper mantle in central Australia.

Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information

Currently, genomic evaluations use multiple-step procedures, which are prone to biases and errors. A single-step procedure may be applicable when genomic predictions can be obtained by modifying the numerator relationship matrix A to H = A + AΔ, where AΔ includes deviations from expected relationships. However, the traditional mixed model equations require H−1, which is usually difficult to obtain for large pedigrees. The computations with H are feasible when the mixed model equations are expressed in an alternate form that also applies for singular H and when those equations are solved by the conjugate gradient techniques. Then the only computations involving H are in the form of Aq or AΔq, where q is a vector. The alternative equations have a nonsymmetric left-hand side. Computing AΔq is inexpensive when the number of nonzeros in AΔ is small, and the product Aq can be calculated efficiently in linear time using an indirect algorithm. Generalizations to more complicated models are proposed. The data included 10.2million final scores on 6.2million Holsteins and were analyzed by a repeatability model. Comparisons involved the regular and the alternative equations. The model for the second case included simulated AΔ. Solutions were obtained by the preconditioned conjugate gradient algorithm, which works only with symmetric matrices, and by the bi-conjugate gradient stabilized algorithm, which also works with nonsymmetric matrices. The convergence rate associated with the nonsymmetric solvers was slightly better than that with the symmetric solver for the original equations, although the time per round was twice as much for the nonsymmetric solvers. The convergence rate associated with the alternative equations ranged from 2 times lower without AΔ to 3 times lower for the largest simulated AΔ. When the information attributable to genomics can be expressed as modifications to the numerator relationship matrix, the proposed methodology may allow the upgrading of an existing evaluation to incorporate the genomic information.

Multigrid preconditioned conjugate-gradient solver for mixed finite-element method

The mixed finite-element approximation to a second-order elliptic PDE results in a saddle-point problem and leads to an indefinite linear system of equations. The mixed system of equations can be transformed into coupled symmetric positive-definite matrix equations, or a Schur complement problem, using block Gauss elimination. A preconditioned conjugate-gradient algorithm is used for solving the Schur complement problem. The mixed finite-element method is closely related to the cell-centered finite difference scheme for solving second-order elliptic problems with variable coefficients. For the cell-centered finite difference scheme, a simple multigrid algorithm can be defined and used as a preconditioner. For distorted grids, an additional iteration is needed. Nested iteration with a multigrid preconditioned conjugate gradient inner iteration results in an effective numerical solution technique for the mixed system of linear equations arising from a discretization on distorted grids. Numerical results show that the preconditioned conjugate-gradient inner iteration is robust with respect to grid size and variability in the hydraulic conductivity tensor.

On the preconditioned conjugate gradient solution of a Stokes problem with Robin-type boundary conditions

We discuss the solution by a preconditioned conjugate gradient algorithm of a Stokes problem with Robin-type boundary conditions. This Stokes problem is encountered when applying an appropriate operator-splitting scheme to the time-discretization of a system modeling the interaction of an incompressible viscous fluid with a deformable thin elastic structure. The main contribution of this Note is the identification of a preconditioner operating in the pressure space that reduces substantially the number of iterations when compared to a conjugate gradient algorithm equipped with the canonical scalar product of L 2 . The results of numerical experiments show the validity of our approach. To cite this article: R. Glowinski, G. Guidoboni, C. R. Acad. Sci. Paris, Ser. I 347 (2009).