Preconditioned Conjugate Gradient Algorithm Research Articles

High-efficiency computation for electromagnetic forming process: An explicit-implicit GPU approach

Low computational efficiency is one of the crucial challenges in electromagnetic forming process simulation. However, the huge gap in switching from the electromagnetic field to mechanics hinders the development of high-efficiency parallel coupling. A graphics processing unit (GPU)-based explicit-implicit coupling method is proposed to address the high-computational cost in simulations of the electromagnetic forming process, in which the large deformation (explicit) and electromagnetic field (implicit) are intimately coupled. The electromagnetic-mechanical coupling is achieved by constructing five parallel algorithms, which include the parallel operations of vector-vector, matrix-matrix, matrix-vector, vector assembly, and matrix-free assembly parallel preconditioned conjugate gradient algorithm. The high-efficiency GPU parallel computing is realized through the special data storage structure and parallel strategies. The explicit-implicit GPU method is validated against the experimental results, and several examples are presented for simulations of the electromagnetic forming process to illustrate both the accuracy and high-efficiency parallel acceleration performance of the present program.

Multi-dimensional visual data completion via weighted hybrid graph-Laplacian

Graph Laplacian has been extensively studied as a priori to assist in visual completion of partially observed entries. However, a simple treatment using first-order operation may be suboptimal since the inherent smooth property within the visual data is ill-considered. To relieve this issue, in this paper, we propose an improved priori, namely weighted hybrid graph Laplacian (WHGL), to better characterize the spatially hierarchical knowledge hidden in the multi-dimensional data. On the whole, two individual layers are presented to encode the structured properties of the given data. In the first layer, the underlying adjacency matrix of the local and nonlocal relationships within the data is composed to capture the pixel-level structure information. Global smoothness is considered as the second layer to reveal the piecewise correlations of different visual elements. Mathematically, we generalize the graph Laplacian in the first- and second-order to construct the respective layers. Finally, a re-weighting strategy is used to mitigate the discontinuity of the reconstructed data under extreme cases, e.g. 95% voxels missing. With an initial step borrowed from some off-the-shelf methods, our model can be solved under limited iterations by the Preconditioned Conjugate Gradients algorithm. The experimental results demonstrate the superiority of the proposed method.

Three-Dimensional Direct Current Resistivity Finite Element Method Forward Modelling Based on Jacobi Preconditioned Conjugate Gradient Algorithm

The calculation speed and accuracy of resistivity forward model are important for direct current exploration. This study proposes a calculation method for the forward model. First, we structured the unstructured grid to enable dissection of the three-dimensional finite element space. We then solved the large sparse linear equations, using the preconditioned conjugate gradient method. The diagonal matrix in the Jacobi iteration was used as the pre-processed matrix, which had the characteristics of a convenient inversion solution and consumed less memory; therefore, it substantially increased the speed to solve large linear equations. The forward modelling performed using synthetic data shows that the preconditioned conjugate gradient algorithm exhibits fast calculation speed and high accuracy. These results are in accordance with the potential distribution rule of the underground electric field. The results have value for three-dimensional direct current exploration in the field of engineering. Furthermore, it provides a reference for solving large linear equations of the finite element method in other fields.

Two New Preconditioned Conjugate Gradient Methods for Minimization Problems

In application to general function, each of the conjugate gradient and Quasi-Newton methods has particular advantages and disadvantages. Conjugate gradient (CG) techniques are a class of unconstrained optimization algorithms with strong local and global convergence qualities and minimal memory needs. Quasi-Newton methods are reliable and efficient on a wide range of problems and they converge faster than the conjugate gradient method and require fewer function evaluations but they have the disadvantage of requiring substantially more storage and if the problem is ill-conditioned, they may take several iterations. A new class has been developed, termed preconditioned conjugate gradient (PCG) method. It is a method that combines two methods, conjugate gradient and Quasi-Newton. In this work, two new preconditioned conjugate gradient algorithms are proposed namely New PCG1 and New PCG2 to solve nonlinear unconstrained optimization problems. A new PCG1 combines conjugate gradient method Hestenes-Stiefel (HS) with new self-scaling symmetric Rank one (SR1), and a new PCG2 combines conjugate gradient method Hestenes-Stiefel (HS) with new self-scaling Davidon, Flecher and Powell (DFP). The algorithm uses the strong Wolfe line search condition. Numerical comparisons with standard preconditioned conjugate gradient algorithms show that for these new algorithms, computational scheme outperforms the preconditioned conjugate gradient.

Algorithm Research Based on the PDE Sensitivity Filter

The PDE sensitivity filter can eliminate the checkerboard patterns and numerical instability existing in the topology optimization results of continuum structures, and the essence of the PDE sensitivity filter is the Helmholtz partial differential equation with Neumann boundary conditions. To solve the large-scale PDE sensitivity filter problem, the conjugate gradient algorithm, the multigrid algorithm and the multigrid preconditioned conjugate gradient algorithm were used to solve the algebraic equations obtained by finite element analysis, and the effects of accuracy, filter radius and grid numbers on the efficiency of topology optimization were studied. The results show that, compared with the conjugate gradient algorithm and the multi-grid algorithm, the multi-grid preconditioned conjugate gradient algorithm has the least number of iterations and the shortest running time, which greatly improves the efficiency of topology optimization.

A cascade of preconditioned conjugate gradient networks for accelerated magnetic resonance imaging

Background and objective: Recent unfolding based compressed sensing magnetic resonance imaging (CS-MRI) methods only reinterpret conventional CS-MRI optimization algorithms and, consequently, inherit the weaknesses of the alternating optimization strategy. In order to avoid the structural complexity of the alternating optimization strategy and achieve better reconstruction performance, we propose to directly optimize the ℓ1 regularized convex optimization problem using a deep learning approach. Method: In order to achieve direct optimization, a system of equations solving the ℓ1 regularized optimization problem is constructed from the optimality conditions of a novel primal-dual form proposed for the effective training of the sparsifying transform. The optimal solution is obtained by a cascade of unfolding networks of the preconditioned conjugate gradient (PCG) algorithm trained to minimize the mean element-wise absolute difference (ℓ1 loss) between the terminal output and ground truth image in an end-to-end manner. The performance of the proposed method was compared with that of U-Net, PD-Net, ISTA-Net+, and the recently proposed projection-based cascaded U-Net, using single-coil knee MR images of the fastMRI dataset. Results: In our experiment, the proposed network outperformed existing unfolding-based networks and the complex version of U-Net in several subsampling scenarios. In particular, when using the random Cartesian subsampling mask with 25 % sampling rate, the proposed model outperformed PD-Net by 0.76 dB, ISTA-Net+ by 0.43 dB, and U-Net by 1.21 dB on the positron density without suppression (PD) dataset in term of peak signal to noise ratio. In comparison with the projection-based cascade U-Net, the proposed algorithm achieved approximately the same performance when the sampling rate was 25% with only 1.62% number of network parameters at the cost of a longer reconstruction time (approximately twice). Conclusion: A cascade of unfolding networks of the PCG algorithm was proposed to directly optimize the ℓ1 regularized CS-MRI optimization problem. The proposed network achieved improved reconstruction performance compared with U-Net, PD-Net, and ISTA-Net+, and achieved approximately the same performance as the projection-based cascaded U-Net while using significantly fewer network parameters.

Pipelined Preconditioned Conjugate Gradient Methods for real and complex linear systems for distributed memory architectures

Preconditioned Conjugate Gradient (PCG) is a popular method for solving large and sparse linear systems of equations. The performance of PCG at scale is affected due to the costly global synchronization steps that arise in dot-products on distributed memory systems. Pipelined PCG (PIPECG) removes the costly global synchronization steps from PCG by only executing a single non-blocking allreduce per iteration and overlapping it with independent computations.In our previous work, we have developed a novel pipelined PCG algorithm called PIPECG-OATI (One Allreduce per Two Iterations) for real linear systems which executes a single non-blocking allreduce per two iterations and provides a large overlap of global communication with independent computations at higher number of cores. Our method achieves this overlap by using iteration combination and by introducing new recurrence and non-recurrence computations. We implement optimizations in the PIPECG-OATI method to use cache memory efficiently.In this work, we present PIPECG-OATI-c method for linear systems with complex Hermitian positive definite and complex symmetric matrices. We compare our method with various pipelined CG methods on a variety of problems and demonstrate that our method always gives the least run times. We performed experiments with our method using 20M and 30M unknowns on up to 16K cores and obtained up to 2.48X performance improvement over PCG and 2.14X performance improvement over PIPECG methods. We also experimented with up to 1-billion unknowns on 16K cores, the largest problem size explored for the CG problem, to our knowledge, and obtained about 25% improvement over PCG.

Least-Squares Reverse Time Migration for Reflection-Angle-Dependent Reflectivity

Least-squares reverse time migration (LSRTM) can estimate high-quality reflectivity of subsurface medium from seismic data. However, the subsurface reflectivity depends on reflection angles, and its variations over reflection angles are extremely important because they can be used to estimate sub-surface physical properties for seismic interpretation. We present a new formulation of the LSRTM method that can estimate reflection-angle-dependent reflectivity from seismic data. We derive a forward modeling operator which predicts the reflection data without calculating the reflection angles, and verify that it approximately equals to the reflection-angle-dependent wave-equation-based Kirchhoff modeling operator under the assumption that the velocity perturbation is small and the reflection angle is smaller than the critical angle. Based on the proposed modeling operator associated with the adjoint of the angle-dependent wave-equation-based Kirchhoff modeling operator, we reformulate LSRTM as an inverse problem to invert for reflection-angle-dependent reflectivity using a preconditioned conjugate gradient algorithm. The algorithm uses a low-rank filter as the preconditioner to attenuate migration artifacts. Imaging tests on synthetic and field seismic data are used to verify validity and superiority of the proposed method. The tests illustrate that the proposed method can produce the reflection-angle-dependent reflectivity with much higher signal-to-noise ratio, resolution and amplitude fidelity than reverse time migration. Compared with conventional LSRTM, it can produce more focused stacked image when the migration velocity contains errors. Moreover, conventional LSRTM only produces the angle-independent reflectivity, whereas the proposed method has the feasibility to produce the reflection-angle-dependent reflectivity.

On numerical solution of nonlinear parabolic multicomponent diffusion-reaction problems

Abstract This work continues our previous analysis concerning the numerical solution of the multi-component mass transfer equations. The present test problems are two-dimensional, parabolic, non-linear, diffusion- reaction equations. An implicit finite difference method was used to discretize the mathematical model equations. The algorithm used to solve the non-linear system resulted for each time step is the modified Picard iteration. The numerical performances of the preconditioned conjugate gradient algorithms (BICGSTAB and GMRES) in solving the linear systems of the modified Picard iteration were analysed in detail. The numerical results obtained show good numerical performances.

Gravity Gradient Inversion of Gravity Field and Steady‐State Ocean Circulation Explorer Satellite Data for the Lithospheric Density Structure in the Qinghai‐Tibet Plateau Region and the Surrounding Regions

AbstractThe continuous uplift and evolution of the Qinghai‐Tibet Plateau region (TPR) represents one of the most important Cenozoic geological events on Earth. The deformation mechanism of the lithosphere in the TPR is still under debate. Relevant viewpoints can be informed by investigations of the lithospheric density structure in the Qinghai‐Tibet Plateau and surrounding regions (TP‐SRs). Considering investigating the fine‐scale lithospheric density structure of the large‐scale areas, the GOCE satellite data covering entire TP‐SRs are used in this study. GOCE gravity gradient data are first subjected to several corrections. Then, based on gravity gradient curvature (SI) calculation complemented with the preconditioned conjugate gradient algorithm, the distribution of SI and the lithospheric density structure at depths of 0–80 km in TP‐SRs are obtained. The effects of temperature differences and material composition on the lithospheric density distribution are further distinguished. The results show that in the crust, the density anomalies in the TPR are blocked by the surrounding hard blocks. In the upper mantle, the distribution of the density anomalies in the TPR is predominantly affected by the northward subduction of the Indian Plate and the high heat flux environment. The Longmen Shan Thrust Belt is found to have penetrated the lower crust but has not reached the depth of the upper mantle. The density anomalies in the basins on the northern side of the TPR exhibit features indicative of independent stable blocks.

Elastic Least-Squares Reverse-Time Migration Based on a Modified Acoustic-Elastic Coupled Equation for OBS Four-Component Data

Elastic least-squares reverse-time migration (ELSRTM) yields subsurface high-resolution elastic images from multicomponent seismic data. For ocean bottom seismic (OBS) acquisitions, both pressure and displacement components are measured using four-component (4C) modern receivers. However, these 4C data cannot be effectively simulated and migrated using conventional ELSRTM approaches due to the limitation of the standard elastic wave equation. For this reason, we propose a new ELSRTM method based on a modified acoustic-elastic coupled (AEC) equation. Owing to this equation, either the prediction or the migration of OBS 4C data can be implemented using a one-time wavefield simulation. Incorporation with a perturbation imaging condition, the proposed method can invert for P- and S-wave velocity and density and generate amplitude-preserving images. Besides, it can provide subsurface impedance reflectivity models by stacking the obtained multi-parameter images. In the least-squares migration, a preconditioned conjugate gradient algorithm is implemented based on a multi-parameter diagonal Hessian. Numerical experiments on an uncorrelated layer model and a portion of the Southern Yellow China sea model can demonstrate the effectiveness of this method on improving imaging resolution and accelerating convergence rate.

Three-dimensional forward modelling of gravity field vector and its gradient tensor using the compact difference schemes

SUMMARY The traditional gravity forward modelling methods for solving partial differential equations (PDEs) only can yield second-order accuracy. When computing the gravity field vector and gradient tensor from the obtained potential, those numerical differentiation approaches will inevitably lose accuracy. To mitigate this issue, we propose an efficient and accurate 3-D forward modelling algorithm based on a fourth-order compact difference scheme. First, a 19-point fourth-order compact difference scheme with general meshsizes in x-, y- and z-directions is adopted to discretize the governing 3-D Poisson’s equation. The resulting symmetric positive-definite linear systems are solved by the pre-conditioned conjugate gradient algorithm. To obtain the first-order (i.e. the gravity field vector) and second-order derivatives (i.e. the gravity gradient tensor) with fourth-order accuracy, we seek to solve a sequence of tridiagonal linear systems resulting from the above mentioned finite difference discretization by using fast Thomas algorithm. Finally, two synthetic models and a real topography relief are used to verify the accuracy of our method. Numerical results show that our method can yield a nearly fourth-order accurate approximation not only to the gravitational potential, but also to the gravity field vector and its gradient tensor, which clearly demonstrates its superiority over the traditional PDE-based methods.

Inversion of high-amplitude magnetic total field anomaly: an application to the Mengku iron-ore deposit, northwest China

In magnetic prospecting, the total field anomaly formula that represents the projection of the magnetic anomaly vector on the geomagnetic field is widely used because it simplifies the calculation of forward modelling and inversion of magnetic data. However, the projection anomaly yields errors relative to the true observed magnetic anomaly, especially for high-amplitude magnetic anomalies such as in iron orebody and unexploded ordnance prospecting. In this study, we analyse the difference between the projection anomaly and observed modulus difference anomaly with physical parameters, and propose to directly invert for the modulus difference anomaly by constructing a nonlinear matrix equation between the model corrections and data corrections. The inversion is then implemented using a preconditioned conjugate gradient algorithm. Synthetic and field magnetic data were used to test the inversion method. Comparison of the two types of total field anomalies shows that the error of the projection anomaly increased with increasing total-field magnetic anomaly. When the total-field magnetic anomaly was < 5,000 nT, the difference between the projection anomaly and modulus difference anomaly results can be ignored. For high-amplitude magnetic anomalies, the modulus difference anomaly inversion produced more accurate representations of both the shape and location of the magnetic sources.

Susceptibility and remanent magnetization inversion of magnetic data with a priori information of the Köenigsberger ratio

SUMMARY Magnetic susceptibility and natural remanent magnetization of rocks are useful parameters to study geological structures and geodynamic processes. Traditional widely used algorithms for the inversion of magnetic data can recover the distribution of the apparent susceptibility or total magnetization intensity, but do not provide information on the remanent magnetization. In this paper, we propose a framework to directly invert for the magnetic susceptibility and the natural remanent magnetization vector using surface or airborne magnetic data, assuming that the Köenigsberger ratio of the rock is known or approximately deducible. The susceptibility and remanence are computed using two different approaches: (1) the susceptibility, intensity, and direction of the remanent magnetization are continuously recovered for each discretized cell and (2) the remanence direction is assumed to be uniform in each subzone and is iteratively computed as discrete values. Both processes are implemented using the preconditioned conjugate gradient algorithm. The method is tested on three synthetic models and one field data set from the Zaohuohexi iron-ore deposit, Qinghai Province, northwest China. The results of the continuous inversion show the trend of the remanent magnetization directions, while the discrete inversion yields more specific values. This inversion framework can determine the source bodies’ geometry and position, and also provide superposed and comprehensive information on the natural remanent magnetization, which may be useful to investigate geological bodies bearing stable primary remanent magnetization.

Convolution operators for visual tracking based on spatial–temporal regularization

In recent years, the method based on discriminative correlation filter has been shown excellent performance in short-term visual tracking. However, discriminative correlation filter-based method heavily suffers from the problem of the multiple peaks and model drift in responds maps incurred by occlusion and rotation. To solve the above problem, we proposed convolution operators for visual tracking based on spatial–temporal regularization. Firstly, we add spatial–temporal regularization in loss function, which will guarantee continuity of the model in time. And we use preconditioned conjugate gradient algorithm to obtain filter coefficients. Secondly, we proposed channel reliability to estimate quality of the learned filter and fuse the different reliability coefficients to weight response map in location. We set a threshold to reduce the number of iteration in location and accelerate the compute speed of algorithm. Finally, we use two different correlation filters to estimate location and scale of target, respectively. Extensively experiment in five video sequences show that our tracker has been significantly improved performance in case of occlusion and rotation. The AUC in success plot improves 33.2% than ECO-HC and 41.5% than STRCF, respectively.

Application of Tikhonov Regularization for 1D Geothermal Heat Flux Ill-Posed Inverse Problem: A Case Study on Chad Sedimentary Basin

Subsurface heat flux information is important in geothermal exploration. With the information, geophysicists can map exactly the thermal potential in a particular area. Based on the surface heat flux, inverse modeling produces the 1D subsurface heat flux distribution. However, inverse problems in the geothermal system are generally ill-posed. Small changes in the data can cause large changes in the solution and the solution may not be unique. To solve the mentioned non-linear and ill-posed equation above, Tikhonov regularization is a choice for stabilizing the inverse calculation. This paper demonstrates how Tikhonov regularization is useful to solve subsurface heat flux distribution both in the synthetic model and real model. Based on surface heat flux distribution from the direct problem, the preconditioned conjugate gradient algorithm calculates the subsurface heat flux. With the correct choice of the regularization parameter, the inverse model fits the initial model. For the testing purposes in real-world conditions, Chad sedimentary basin located in Chad and Nigeria is used as a model. A high geothermal gradient is found in this area. Therefore, geothermal explorations are on the rise recently. Its thermal conductivity, heat production, and stratigraphy data from previous researches provide information about the initial model. The heat flux curve generated from inversion matches the initial noisy model with the error of around 10−9 mW/m2. Therefore, to answer the increasing energy demand, this method can be highly applicable to future geothermal prospecting.

Shot-record extended model domain preconditioners for least-squares migration

We have studied the preconditioned conjugate gradient (CG) algorithm in the context of shot-record extended model domain least-squares migration. The CG algorithm is a powerful iterative technique that can solve the least-squares migration problem efficiently; however, to see the merits of least-squares migration, one needs to apply the algorithm for several iterations. Generally speaking, the convergence rate of the CG algorithm depends on the condition number of the operator. Preconditioners are a family of operators that are easy to build and invert. Proper preconditioners can cluster the eigenvalues of the original operator; hence, they reduce the condition number of the operator that one wishes to invert. Accordingly, preconditioning the operator can, in theory, improve the convergence rate of the algorithm. In least-squares migration, the diagonal scaling of the Hessian and the approximated inverse of the Hessian are proven to work well as a preconditioner. We develop and apply two types of preconditioners for the shot-record extended model domain least-squares migration problem. The first preconditioner belongs to the diagonal scaling category, and a second preconditioner is a filter-based approach, which approximates the partial Hessian operators by local convolutional filters. The goal is to increase the convergence rate of the shot-record extended model domain least-squares migration using the reformulated cost function with a preconditioned operator. Experiments with a synthetic Sigsbee model and a real data example from the Gulf of Mexico, Mississippi Canyon data set, indicate that preconditioning the linear system of the equations improves the convergence rate of the algorithm.

Inversion of the density structure of the lithosphere in the North China Craton from GOCE satellite gravity gradient data

Lithospheric thinning in the North China Craton (NCC), one of the oldest cratons on Earth, has been a focus of geoscientists across the globe for a long time. In this study, seismic tomographic P-wave velocity variations are converted to density perturbations as the initial constraint values, and the four components of the GOCE Level-2 gravity gradient product are subjected to topographic effect correction, relief correction on the Moho and sedimentary layer boundary, and long-wavelength correction. In addition, the anomalous gravity gradient effect resulting from the uncertainty of the Moho and sedimentary layer depth is considered. Considering the effects of temperature variations in the lithosphere, which leads to the non-uniform density distribution, the anomalous gravity gradient effect resulting from the temperature variations is corrected at the first time in this area. Inverse computation is performed based on the preconditioned conjugate gradient algorithm. Because the Lagrange empirical parameter in the algorithm exhibits uncertainty in data inversion, the regularization parameter, which is typically an empirical parameter, is replaced with the value of the inflection point of the L-curve. The results show the followings. The non-uniform density distribution in the lithosphere is related not only to the composition but also to the internal temperature variations in the lithosphere. The density of the lithosphere in the NCC is significantly non-uniform in both the horizontal and vertical directions and exhibits a notable segment-wise spatial distribution pattern. Two tectonic units, namely the Taihangshan tectonic zone and the Linfen–Weihe graben, constitute a central gravity gradient transition zone. There is a significant difference in the density distribution on the two sides of this central transition zone.

Non-overlapping domain decomposition solution schemes for structural mechanics isogeometric analysis

Isogeometric analysis (IGA) is a novel computer aided engineering technique that addresses diverse problems in computational mechanics [1–4], all under the exact geometric representation. Apart from the exact geometric representation, the high continuity of IGA shape functions enhances the accuracy and robustness of the method. However, the price to pay is that the resulting matrices are denser, with increased bandwidth and overlapping which make the solution of large-scale problems more computationally intensive. As a result, effective solution techniques are still considered an open issue for further research. In this paper, an innovative family of solution schemes based on domain decomposition methods (DDM) is proposed, that significantly reduces the computational cost. Specifically, the solution of the global system is performed with the preconditioned conjugate gradient algorithm (PCG) whose preconditioning step is evaluated with a dual DDM where special care is taken to avoid the overlapped subdomains which are inherent in decomposed IGA formulation due to the increased continuity of the shape functions.

Massively Parallel Implementation of Divide-and-Conquer Jacobi Iterations Using Particle-Mesh Ewald for Force Field Polarization.

To accelerate the evaluation of the self-consistent polarization in large condensed-phase systems with polarizable force fields, the new divide-and-conquer Jacobi iterations (DC-JI) solver is adapted for periodic boundary conditions with particle-mesh Ewald and implemented in a massively parallel fashion within the Tinker-HP software package. DC-JI captures the mutual polarization of close-range interactions within subclusters of atoms using Cholesky decomposition and couples in the polarization effects between these clusters iteratively. Iterative convergence is accelerated with direct inversion of the iterative subspace (DIIS) extrapolation. Compared to widely used preconditioned conjugate gradient (PCG) or conventional Jacobi iterations (JI/DIIS) algorithms, DC-JI/DIIS solves the polarization equations ∼20-30% faster in protein systems ranging from ∼10,000-175,000 atoms run on hundreds of processor cores. This translates to ∼10-15% speed-ups in the number of nanoseconds of simulation time that can be achieved per day. Not only is DC-JI/DIIS faster than PCG, but it also gives more energetically robust solutions for a given convergence threshold. These improvements make numerically robust polarizable force field simulations more computationally tractable for chemical systems of interest.