Preconditioned Conjugate Gradient Iteration Research Articles

Logarithmic total variation regularization via preconditioned conjugate gradient method for sparse reconstruction of bioluminescence tomography

Background and objectiveBioluminescence Tomography (BLT) is a powerful optical molecular imaging technique that enables the noninvasive investigation of dynamic biological phenomena. It aims to reconstruct the three-dimensional spatial distribution of bioluminescent sources from optical measurements collected on the surface of the imaged object. However, BLT reconstruction is a challenging ill-posed problem due to the scattering effect of light and the limitations in detecting surface photons, which makes it difficult for existing methods to achieve satisfactory reconstruction results. In this study, we propose a novel method for sparse reconstruction of BLT based on a preconditioned conjugate gradient with logarithmic total variation regularization (PCG-logTV). MethodThis PCG-logTV method incorporates the sparsity of overlapping groups and enhances the sparse structure of these groups using logarithmic functions, which can preserve edge features and achieve more stable reconstruction results in BLT. To accelerate the convergence of the algorithm solution, we use the preconditioned conjugate gradient iteration method on the objective function and obtain the reconstruction results. We demonstrate the performance of our proposed method through numerical simulations and in vivo experiment. Results and conclusionsThe results show that the PCG-logTV method obtains the most accurate reconstruction results, and the minimum position error (LE) is 0.254mm, which is 26%, 31% and 34% of the FISTA (0.961), IVTCG (0.81) and L1-TV (0.739) methods, and the root mean square error (RMSE) and relative intensity error (RIE) are the smallest, indicating that it is closest to the real light source. In addition, compared with the other three methods, the PCG-logTV method also has the highest DICE similarity coefficient, which is 0.928, which means that this method can effectively reconstruct the three-dimensional spatial distribution of bioluminescent light sources, has higher resolution and robustness, and is beneficial to the preclinical and clinical studies of BLT.

Unleashing the Power of Graph Spectral Sparsification for Power Grid Analysis via Incomplete Cholesky Factorization

Graph spectral sparsification based preconditioning technique has shown promising results for power grid analysis. However, the conventional methods converge slowly for high accuracy requirement. In this work, we propose an efficient approach to address this issue. Instead of using Cholesky factorization, we employ incomplete Cholesky factorization to factorize the spectral sparsifier. We also propose a concept of graph spectral pattern, which can further reduce the preconditioned conjugate gradient (PCG) iterations using less number of nonzeros. Experiments show that under 10-6 relative tolerance, our proposed preconditioning technique achieves 1.17× speedup compared to AMGPCG in average; compared to the conventional spectral sparsification based preconditioning techniques, our proposed approach achieves up to 8.53× speedup of the factorization, 8.74× speedup of the PCG iteration, and 5.6× speedup of the total time. Moreover, the speedup of the total time continues to enlarge for higher accuracy requirement, e.g., 10-12. Last, but not the least, our method is compatible with existing graph spectral sparsification algorithms for power grid analysis.

Semiexplicit time integration of a reduced magnetic vector potential magneto‐quasistatic field formulation

AbstractIn this article, a reduced vector potential formulation for transient magneto‐quasistatic field problems is presented. The finite element method for spatial discretization yields a system of differential algebraic equations (DAE). Applying a generalized Schur complement the DAE system decomposes into a system of ordinary differential equations (ODE) and an algebraic equation. In order to avoid computationally expensive linearization schemes to solve nonlinear eddy current problems the ODE system is integrated in time using a semiexplicit Euler scheme. The vector potential is split into a source vector potential and a reduced vector potential which describes the eddy current and magnetization effects. The reduced formulation is compared with the underlying nonreduced formulation regarding the overall performance. In addition, the proper orthogonal decomposition method is applied to the iterative preconditioned conjugate gradient method to minimize the number of iterations per time step.

Fast solution of three‐dimensional elliptic equations with randomly generated jumping coefficients by using tensor‐structured preconditioners

AbstractIn this paper, we propose and analyze the numerical algorithms for fast solution of periodic elliptic problems in random media in , . Both the two‐dimensional (2D) and three‐dimensional (3D) elliptic problems are considered for the jumping equation coefficients built as a checkerboard type configuration of bumps randomly distributed on a large , or lattice, respectively. The finite element method discretization procedure on a 3D uniform tensor grid is described in detail, and the Kronecker tensor product approach is proposed for fast generation of the stiffness matrix. We introduce tensor techniques for the construction of the low Kronecker rank spectrally equivalent preconditioner in a periodic setting to be used in the framework of the preconditioned conjugate gradient iteration. The discrete 3D periodic Laplacian pseudo‐inverse is first diagonalized in the Fourier basis, and then the diagonal matrix is reshaped into a fully populated third‐order tensor of size . The latter is approximated by a low‐rank canonical tensor by using the multigrid Tucker‐to‐canonical tensor transform. As an example, we apply the presented solver in numerical analysis of stochastic homogenization method where the 3D elliptic equation should be solved many hundred times, and where for every random sampling of the equation coefficient one has to construct the new stiffness matrix and the right‐hand side. The computational characteristics of the presented solver in terms of a lattice parameter and the grid‐size, , in both 2D and 3D cases are illustrated in numerical tests. Our solver can be used in various applications where the elliptic problem should be solved for a number of different coefficients for example, in many‐particle dynamics, protein docking problems or stochastic modeling.

Learning a preconditioner to accelerate compressed sensing reconstructions in MRI.

PurposeTo learn a preconditioner that accelerates parallel imaging (PI) and compressed sensing (CS) reconstructions.MethodsA convolutional neural network (CNN) with residual connections was used to train a preconditioning operator. Training and validation data were simulated using 50% brain images and 50% white Gaussian noise images. Each multichannel training example contains a simulated sampling mask, complex coil sensitivity maps, and two regularization parameter maps. The trained model was integrated in the preconditioned conjugate gradient (PCG) method as part of the split Bregman CS method. The acceleration performance was compared with that of a circulant PI‐CS preconditioner for varying undersampling factors, number of coil elements and anatomies.ResultsThe learned preconditioner reduces the number of PCG iterations by a factor of 4, yielding a similar acceleration as an efficient circulant preconditioner. The method generalizes well to different sampling schemes, coil configurations and anatomies.ConclusionIt is possible to learn adaptable preconditioners for PI and CS reconstructions that meet the performance of state‐of‐the‐art preconditioners. Further acceleration could be achieved by optimizing the network architecture and the training set. Such a preconditioner could also be integrated in fully learned reconstruction methods to accelerate the training process of unrolled networks.

Tensor method for optimal control problems constrained by fractional three‐dimensional elliptic operator with variable coefficients

AbstractWe introduce the tensor numerical method for solving optimal control problems that are constrained by fractional two‐ (2D) and three‐dimensional (3D) elliptic operators with variable coefficients. We solve the governing equation for the control function which includes a sum of the fractional operator and its inverse, both discretized over large 3D spacial grids. Using the diagonalization of the arising matrix‐valued functions in the eigenbasis of the one‐dimensional Sturm–Liouville operators, we construct the rank‐structured tensor approximation with controllable precision for the discretized fractional elliptic operators and the respective preconditioner. The right‐hand side in the constraining equation (the optimal design function) is supposed to be represented in a form of a low‐rank canonical tensor. Then the equation for the control function is solved in a tensor structured format by using preconditioned CG iteration with the adaptive rank truncation procedure that also ensures the accuracy of calculations, given an ‐threshold. This method reduces the numerical cost for solving the control problem to (plus the quadratic term with a small weight), which outperforms traditional approaches with complexity in the 3D case. The storage for the representation of all 3D nonlocal operators and functions involved is also estimated by . This essentially outperforms the traditional methods operating with fully populated matrices and vectors in . Numerical tests for 2D/3D control problems indicate the almost linear complexity scaling of the rank truncated preconditioned conjugate gradient iteration in the univariate grid size n.

Practical implementation of genetic groups in single-step genomic evaluations with Woodbury matrix identity–based genomic relationship inverse

The growing amount of genomic information in dairy cattle has increased computational and modeling challenges in the single-step evaluations. The computational challenges are due to the dense inverses of genomic (G) and pedigree (A22) relationship matrices of genotyped animals in the single-step mixed model equations. An equivalent mixed model equation is given by single-step genomic BLUP that are based on the T matrix (ssGTBLUP), where these inverses are avoided by expressing G-1 through a product of 2 rectangular matrices, and (A22)-1 through sparse matrix blocks of the inverse of full relationship matrix A-1. A proper way to account genetic groups through unknown parent groups (UPG) after the Quaas-Pollak transformation (QP) is one key factor in a single-step model. When the UPG effects are incompletely accounted, the iterative solving method may have convergence problems. In this study, we investigated computational and predictive performance of ssGTBLUP with residual polygenic (RPG) effect and UPG. The QP transformation used A-1 and, in the complete form, T and (A22)-1 matrices as well. The models were tested with official Nordic Holstein milk production test-day data and model. The results show that UPG can be easily implemented in ssGTBLUP having RPG. The complete QP transformation was computationally feasible when preconditioned conjugate gradient iteration and iteration on data without explicitly setting up G or A22 matrices were used. Furthermore, for good convergence of the preconditioned conjugate gradient method, a complete QP transformation was necessary.

A new parameter-free method for Toeplitz systems of weakly nonlinear equations

Based on the fact that the linear term is strongly dominant over the nonlinear term, by using the approximate inverse-free preconditioned conjugate gradient (AIPCG) iteration technique, we establish the Picard-AIPCG iteration method to solve Toeplitz systems of weakly nonlinear equations. Since the storage and the accurate computation of Jacobian matrix are not necessary in this method and the only need is to solve the constant Toeplitz sub-systems, there may be a great saving in computer storage and calculation workload when facing the practical problems. The global convergence has been proved under some suitable conditions. Some numerical examples demonstrate the feasibility and effectiveness of this new iteration.

FeGRASS: Fast and Effective Graph Spectral Sparsification for Scalable Power Grid Analysis

Graph spectral sparsification aims to find a ultrasparse subgraph which can preserve the spectral properties of the original graph. The subgraph can be leveraged to construct a preconditioner to speed up the solution of the original graph’s Laplacian matrix. In this work, we propose feGRASS, a fast and effective graph spectral sparsification approach for the problem of large-scale power grid analysis and other problems with similar graphs. The proposed approach is based on two novel concepts: 1) effective edge weight and 2) spectral edge similarity. The former takes advantage of node degrees and breadth-first-search (BFS) distances, which leads to a scalable algorithm for generating low-stretch spanning trees (LSSTs). Then, the latter concept is leveraged during the recovery of spectrally critical off-tree edges to produce spectrally similar subgraphs. Compared with the most recent competitor <xref ref-type="bibr" rid="ref1" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">[1]</xref> , the proposed approach is much faster for producing high-quality spectral sparsifiers. Extensive experimental results have been demonstrated to illustrate the superior efficiency of a preconditioned conjugate gradient (PCG) algorithm based on the proposed approach, for solving large power grid problems and many other real-world graph Laplacians. For instance, a power grid matrix with 60 million unknowns and 260 million nonzeros can be solved (at a 1E-3 accuracy level) within 196 s and 12 PCG iterations, on a single CPU core.

An augmented wavelet reconstructor for atmospheric tomography

Atmospheric tomography, i.e., the reconstruction of the turbulence profile in the atmosphere, is a challenging task for adaptive optics (AO) systems for the next generation of extremely large telescopes. Within the AO community, the solver of first choice is the so-called Matrix Vector Multiplication (MVM) method, which directly applies the (regularized) generalized inverse of the system operator to the data. For small telescopes this approach is feasible, however, for larger systems such as the European Extremely Large Telescope (ELT), the atmospheric tomography problem is considerably more complex, and the computational efficiency becomes an issue. Iterative methods such as the Finite Element Wavelet Hybrid Algorithm (FEWHA) are a promising alternative. FEWHA is a wavelet-based reconstructor that uses the well-known iterative preconditioned conjugate gradient (PCG) method as a solver. The number of floating point operations and the memory usage are decreased significantly by using a matrix-free representation of the forward operator. A crucial indicator for the real-time performance are the number of PCG iterations. In this paper, we propose an augmented version of FEWHA, where the number of iterations is decreased by 50% using an augmented Krylov subspace method. We demonstrate that a parallel implementation of augmented FEWHA allows the fulfilment of the real-time requirements of the ELT.

Hierarchical Interpolative Factorization Preconditioner for Parabolic Equations

This note proposes an efficient preconditioner for solving linear and semi-linear parabolic equations. With the Crank-Nicholson time stepping method, the algebraic system of equations at each time step is solved with the conjugate gradient method, preconditioned with hierarchical interpolative factorization. Stiffness matrices arising in the discretization of parabolic equations typically have large condition numbers, and therefore preconditioning becomes essential, especially for large time steps. We propose to use the hierarchical interpolative factorization as the preconditioning for the conjugate gradient iteration. Computed only once, the hierarchical interpolative factorization offers an efficient and accurate approximate inverse of the linear system. As a result, the preconditioned conjugate gradient iteration converges in a small number of iterations. Compared to other classical exact and approximate factorizations such as Cholesky or incomplete Cholesky, the hierarchical interpolative factorization can be computed in linear time and the application of its inverse has linear complexity. Numerical experiments demonstrate the performance of the method and the reduction of conjugate gradient iterations.

Q-Compensated Denoising of Seismic Data

It is widely known that strong noise can decrease the quality of seismic data. However, the anelastic attenuation could be more important to account for the weak amplitude and low quality of seismic data. Here, we develop an inversion framework to simultaneously compensate for the attenuation of seismic data and remove noise, thereby enhancing the quality of seismic data. Instead of directly applying a compensation operator to the input seismic data, we formulate an inverse problem that connects the sparse reflectivity model and the raw seismic data via the convolution and attenuation functions. The random noise is assumed to be the unpredicted part of the forward modeling process. We use the L2-norm regularization for the data misfit and impose a sparsity constraint onto the reflectivity series, e.g., using the L1-norm constraint. We use an iterative preconditioned conjugate gradient method to solve the L1-norm constrained least-squares optimization problem and obtain the reflectivity series. The denoised and compensated data are obtained by applying the convolution operator to the reflectivity. We use several synthetic and field seismic data to illustrate the effectiveness of the presented method.

A GPU-based algorithm for efficient LES of high Reynolds number flows in heterogeneous CPU/GPU supercomputers

Αn optimized MPI+OpenACC implementation model that performs efficiently in CPU/GPU systems using large-eddy simulation is presented. The code was validated for the simulation of wave boundary-layer flows against numerical and experimental data in the literature. A direct Fast-Fourier-Transform-based solver was developed for the solution of the Poisson equation for pressure taking advantage of the periodic boundary conditions. This solver was optimized for parallel execution in CPUs and outperforms by 10 times in computational time a typical iterative preconditioned conjugate gradient solver in GPUs. In terms of parallel performance, an overlapping strategy was developed to reduce the overhead of performing MPI communications using GPUs. As a result, the weak scaling of the algorithm was improved up to 30%. Finally, a large-scale simulation (Re = 2 × 105) using a grid of 4 × 108 cells was executed, and the performance of the code was analyzed. The simulation was launched using up to 512 nodes (512 GPUs + 6144 CPU-cores) on one of the current top 10 supercomputers of the world (Piz Daint). A comparison of the overall computational time showed that the GPU version was 4.2 times faster than the CPU one. The parallel efficiency of this strategy (47%) is competitive compared with the state-of-the-art CPU implementations, and it has the potential to take advantage of modern supercomputing capabilities.

Numerical study in stochastic homogenization for elliptic partial differential equations: Convergence rate in the size of representative volume elements

SummaryWe describe the numerical scheme for the discretization and solution of 2D elliptic equations with strongly varying piecewise constant coefficients arising in the stochastic homogenization of multiscale composite materials. An efficient stiffness matrix generation scheme based on assembling the local Kronecker product matrices is introduced. The resulting large linear systems of equations are solved by the preconditioned conjugate gradient iteration with a convergence rate that is independent of the grid size and the variation in jumping coefficients (contrast). Using this solver, we numerically investigate the convergence of the representative volume element (RVE) method in stochastic homogenization that extracts the effective behavior of the random coefficient field. Our numerical experiments confirm the asymptotic convergence rate of systematic error and standard deviation in the size of RVE rigorously established in Gloria et al. The asymptotic behavior of covariances of the homogenized matrix in the form of a quartic tensor is also studied numerically. Our approach allows laptop computation of sufficiently large number of stochastic realizations even for large sizes of the RVE.

Multigrid preconditioners for optimal control problems with stochastic elliptic PDE constraints

ABSTRACT In this work, we construct multigrid preconditioners to be used in the solution process of pathwise optimal control problems constrained by elliptic partial differential equations with random coefficients. We use a sparse-grid collocation method to discretize in the stochastic space and multigrid techniques in the physical space. Numerical results show that the proposed preconditioners lead to significant computational savings, with the number of preconditioned conjugate gradient iterations decreasing as the resolution in the physical space increases.

Azimuthal Anisotropy of the North American Upper Mantle Based on Full Waveform Inversion

AbstractA new azimuthal anisotropy model for the North American and Caribbean Plates, namely, , is constructed based on full waveform inversion and records from the USArray and other temporary/permanent networks deployed in the study region. A total of 180 earthquakes and 4,516 seismographic stations are employed in the inversion to simultaneously constrain radially and azimuthally anisotropic model parameters: , , , and , within the crust and mantle. Thirty‐two preconditioned conjugate gradient iterations have been utilized to minimize frequency‐dependent phase discrepancies between observed and predicted seismograms for three‐component short‐period (15–40 s) body waves and long‐period (25–100 s) surface waves. Model exhibits complicated variations in anisotropic fabrics underneath the western and eastern United States, especially at depths shallower than 100 km. For instance, the fast axis orientations in model suggest the presence of trench‐perpendicular mantle flows underneath the Cascadia Subduction Zone and also follow the strikes of the Snake River Plain, the Ouachita Orogenic Front, and the Grenville and Appalachian Orogenic Belts. The amplitudes of azimuthal anisotropy reduce to around 1% at depths greater than 200 km, and the orientations are subparallel to the global plate motion directions to the east of the Rocky Mountain, except for large discrepancies in central and eastern Canada. At a depth of 700 km, the fast axes change along the trajectory of the Farallon slab underneath the Great Lakes region and Gulf of Mexico, which might indicate the development of 2‐D poloidal‐mode mantle flows perpendicular to the strike of the sinking slab within the uppermost lower mantle. Comparisons between model with a western U.S. model from ambient noise tomography and SKS splitting measurements demonstrate a relatively good agreement for the fast axis orientations, considering the usage of different data sets and imaging techniques. However, the absolute magnitude of azimuthal anisotropy in model might be underestimated, especially at greater depths, given the poor agreement on the amplitudes of predicted and observed SKS splitting times. At the current stage, the agreement among different azimuthal anisotropy models at global and continental scales is still poor even for the United States with a dense station coverage.

Hydraulic fracturing monitoring: New concept of electromagnetics linked to elastic modeling

We have developed a simulation of the electromagnetic (EM) response and its joint interpretation with fracking monitoring simulation during hydraulic fracturing in an unconventional reservoir. A multiphysics workflow is presented, using a criterion based on a breakdown pressure to generate and propagate the hydraulic fracturing, where the pressure response and EM response were jointly constructed. The approximate solution of Maxwell equations was obtained using a mixed finite-element method (FEM) combined with a leapfrog time-stepping procedure. The spatial discretization of the Maxwell’s equations is achieved using the mixed finite-element spaces of Nedelec, in which the electrical fields have a continuous tangential component across the edges of the computational mesh. Leapfrog time stepping was used for the time discretization, so that the fully discrete sequence is first solving for the electric field using preconditioned conjugate gradient iteration and then solving for the magnetic field. A stability analysis of this FEM method is provided. The results indicate that the EM response might be sensitive enough to be monitored and the magnetic field correlates better with the saturation distribution than the electric field. Furthermore, due to its relation to water saturation and effective connectivity, it is shown that EM monitoring yields additional information to determine the stimulated reservoir volume. An error analysis of estimations from EM and the actual saturation size indicates that the petrophysical assumption is perhaps the most sensitive part of this study.

An Operator Smoothing with ILU(0) for Aggregation-based Algebraic Multigrid

Due to the potential optimal convergence, algebraic multigrid is widely used in solving large sparse linear systems. The aggregation-based version is one of the widely used methods for its cheap cost and easy implementation. But its convergence is often slow compared to other versions. In this paper, a smoothing technique based on incomplete LU factorization without fill-in is presented to the operators on each level. Each operator is approximately factorized, and the derived lower and upper triangular factors are approximately inverted, which are applied to the operator from both sides to improve the diagonal dominance, and then the effectiveness of the smoothing process and the accuracy of the grid transfer operators. The numerical results show that when incorporated into the preconditioned conjugate gradient iteration, the convergence rate is greatly improved, and though the time used for setup is larger, the time used for iteration and the overall time can also be reduced for large-scale systems.

An Improved Domain Decomposition Method for Drop Impact Reliability Analysis of 3-D ICs

Drop test is usually adopted in the integrated circuit (IC) testing to estimate the shock resistance capability of IC packaging. Generally, it is very time-consuming for the drop test simulation, and therefore the fast numerical approach is needed to reduce the computational cost. In this paper, we propose an efficient drop test simulator for through-silicon-via (TSV)-based 3-D IC to simulate its mechanical behaviors under drop impact. The proposed simulator is based on the idea of domain decomposition (DD) to improve the condition number of the coefficient matrix of the solver. We further develop two efficient numerical techniques in the simulator to improve its efficiency. First, a second-order formulation is proposed for the initial solution selection in preconditioned conjugate gradient (PCG) solver, which can efficiently reduce the number of PCG iterations at each time point during simulation. Second, an equivalent structure-embedded model is proposed and applied in the first Schwarz iteration to efficiently reduce the number of Schwarz iterations in DD method. Numerical experiments show that the proposed drop test simulator can achieve $7.03\times $ speedup in comparison with the conventional finite-element method-based solver. It is demonstrated in this paper through several examples with multichip layers, of which each chip layer consists of an $8\times 8$ TSV array, that the proposed simulator can be widely applied to reliability analysis of 3-D ICs under drop impact.

A new aggregation algorithm based on coordinates partitioning recursively for algebraic multigrid method

Aggregation based algebraic multigrid is widely used to solve sparse linear systems, due to its potential to achieve asymptotic optimal convergence and cheap cost to set up. In this kind of method, it is vital to construct coarser grids based on aggregation. In this paper, we provide a new aggregation method based on coordinates partitioning recursively. The adjacent graph of the original coefficient matrix is partitioned into sub-graphs and each sub-graph is recursively partitioned until the minimal number of nodes over the sub-graphs on some level is small enough. In this way, a hierarchy of grids can be constructed from top to bottom, which is completely different from the classical schemes. The results from the solution of model partial differential equations with the preconditioned conjugate gradient iteration show that the new algorithm has better performance and is more robust than the widely used classical algorithms in most cases.