Iterative Linear System Solver Research Articles

AI‐enhanced iterative solvers for accelerating the solution of large‐scale parametrized systems

SummaryRecent advances in the field of machine learning open a new era in high performance computing for challenging computational science and engineering applications. In this framework, the use of advanced machine learning algorithms for the development of accurate and cost‐efficient surrogate models of complex physical processes has already attracted major attention from scientists. However, despite their powerful approximation capabilities, surrogate model predictions are still far from being near to the ‘exact’ solution of the problem. To address this issue, the present work proposes the use of up‐to‐date machine learning tools in order to equip a new generation of iterative solvers of linear equation systems, capable of very efficiently solving large‐scale parametrized problems at any desired level of accuracy. The proposed approach consists of the following two steps. At first, a reduced set of model evaluations is performed using a standard finite element methodology and the corresponding solutions are used to establish an approximate mapping from the problem's parametric space to its solution space using a combination of deep feedforward neural networks and convolutional autoencoders. This mapping serves as a means of obtaining very accurate initial predictions of the system's response to new query points at negligible computational cost. Subsequently, an iterative solver inspired by the Algebraic Multigrid method in combination with Proper Orthogonal Decomposition, termed POD‐2G, is developed that successively refines the initial predictions of the surrogate model towards the exact solution. The application of POD‐2G as a standalone solver or as preconditioner in the context of preconditioned conjugate gradient methods is demonstrated on several numerical examples of large scale systems, with the results indicating its strong superiority over conventional iterative solution schemes.

On Krylov methods for large-scale CBCT reconstruction.

Krylov subspace methods are a powerful family of iterative solvers for linear systems of equations, which are commonly used for inverse problems due to their intrinsic regularization properties. Moreover, these methods are naturally suited to solve large-scale problems, as they only require matrix-vector products with the system matrix (and its adjoint) to compute approximate solutions, and they display a very fast convergence. Even if this class of methods has been widely researched and studied in the numerical linear algebra community, its use in applied medical physics and applied engineering is still very limited. e.g. in realistic large-scale computed tomography (CT) problems, and more specifically in cone beam CT (CBCT). This work attempts to breach this gap by providing a general framework for the most relevant Krylov subspace methods applied to 3D CT problems, including the most well-known Krylov solvers for non-square systems (CGLS, LSQR, LSMR), possibly in combination with Tikhonov regularization, and methods that incorporate total variation regularization. This is provided within an open source framework: the tomographic iterative GPU-based reconstruction toolbox, with the idea of promoting accessibility and reproducibility of the results for the algorithms presented. Finally, numerical results in synthetic and real-world 3D CT applications (medical CBCT andμ-CT datasets) are provided to showcase and compare the different Krylov subspace methods presented in the paper, as well as their suitability for different kinds of problems.

Cross-Interactive Residual Smoothing for Global and Block Lanczos-Type Solvers for Linear Systems with Multiple Right-Hand Sides

Global and block Krylov subspace methods are efficient iterative solvers for large sparse linear systems with multiple right-hand sides. However, global or block Lanczos-type solvers often exhibit large oscillations in the residual norms and may have a large residual gap relating to the loss of attainable accuracy of the approximations. Conventional residual smoothing schemes suppress these oscillations but cannot improve the attainable accuracy, whereas a recent residual smoothing scheme enables the improvement of the attainable accuracy for single right-hand side Lanczos-type solvers. The underlying concept of this scheme is that the primary and smoothed sequences of the approximations and residuals influence one another, thereby avoiding the severe propagation of rounding errors. In the present study, we extend this cross-interactive residual smoothing to the case of solving linear systems with multiple right-hand sides. The resulting smoothed methods can reduce the residual gap with a low additional cost compared to their original counterparts. We demonstrate the effectiveness of the proposed approach through rounding error analysis and numerical experiments.

On adapting Nesterov's scheme to accelerate iterative methods for linear problems

AbstractNesterov's well‐known scheme for accelerating gradient descent in convex optimization problems is adapted to accelerating stationary iterative solvers for linear systems. Compared with classical Krylov subspace acceleration methods, the proposed scheme requires more iterations, but it is trivial to implement and retains essentially the same computational cost as the unaccelerated method. An explicit formula for a fixed optimal parameter is derived in the case where the stationary iteration matrix has only real eigenvalues, based only on the smallest and largest eigenvalues. The fixed parameter and corresponding convergence factor are shown to maintain their optimality when the iteration matrix also has complex eigenvalues that are contained within an explicitly defined disk in the complex plane. A comparison to Chebyshev acceleration based on the same information of the smallest and largest real eigenvalues (dubbed restricted information Chebyshev acceleration) demonstrates that Nesterov's scheme is more robust in the sense that it remains optimal over a larger domain when the iteration matrix does have some complex eigenvalues. Numerical tests validate the efficiency of the proposed scheme. This work generalizes and extends the results of [14, lemmas 3.1 and 3.2 and theorem 3.3].

Recent advances in multi-frequency analysis with the boundary element method

The boundary element discretization of the acoustic Helmholtz equation yields frequency-dependent linear systems, which require efficient solution strategies for multi-frequency analysis. We briefly review recent advances in this field covering frequency approximation and extraction of coefficient matrices, frequency sampling strategies, iterative solvers for frequency-dependent linear systems and model order reduction techniques.

Resilience and fault tolerance in high-performance computing for numerical weather and climate prediction

Progress in numerical weather and climate prediction accuracy greatly depends on the growth of the available computing power. As the number of cores in top computing facilities pushes into the millions, increased average frequency of hardware and software failures forces users to review their algorithms and systems in order to protect simulations from breakdown. This report surveys hardware, application-level and algorithm-level resilience approaches of particular relevance to time-critical numerical weather and climate prediction systems. A selection of applicable existing strategies is analysed, featuring interpolation-restart and compressed checkpointing for the numerical schemes, in-memory checkpointing, user-level failure mitigation and backup-based methods for the systems. Numerical examples showcase the performance of the techniques in addressing faults, with particular emphasis on iterative solvers for linear systems, a staple of atmospheric fluid flow solvers. The potential impact of these strategies is discussed in relation to current development of numerical weather prediction algorithms and systems towards the exascale. Trade-offs between performance, efficiency and effectiveness of resiliency strategies are analysed and some recommendations outlined for future developments.

A Posteriori Error Estimates for Multilevel Methods for Graph Laplacians

In this paper, we study a posteriori error estimators which aid multilevel iterative solvers for linear systems of graph Laplacians. In earlier works such estimates were computed by solving a pertu...

GPBi‐CGstab(L): A Lanczos‐type product method unifying Bi‐CGstab(L) and GPBi‐CG

SummaryLanczos‐type product methods (LTPMs), in which the residuals are defined by the product of stabilizing polynomials and the Bi‐CG residuals, are effective iterative solvers for large sparse nonsymmetric linear systems. Bi‐CGstab(L) and GPBi‐CG are popular LTPMs and can be viewed as two different generalizations of other typical methods, such as CGS, Bi‐CGSTAB, and Bi‐CGStab2. Bi‐CGstab(L) uses stabilizing polynomials of degree L, while GPBi‐CG uses polynomials given by a three‐term recurrence (or equivalently, a coupled two‐term recurrence) modeled after the Lanczos residual polynomials. Therefore, Bi‐CGstab(L) and GPBi‐CG have different aspects of generalization as a framework of LTPMs. In the present paper, we propose novel stabilizing polynomials, which combine the above two types of polynomials. The resulting method is referred to as GPBi‐CGstab(L). Numerical experiments demonstrate that our presented method is more effective than conventional LTPMs.

Probabilistic linear solvers: a unifyingview

Several recent works have developed a new, probabilistic interpretationfor numerical algorithms solving linear systems in which the solution is inferred ina Bayesian framework, either directly or by inferring the unknown action of thematrix inverse. These approaches have typically focused on replicating the behaviourof the conjugate gradient method as a prototypical iterative method. In this work,surprisingly general conditions for equivalence of these disparate methods arepresented. We also describe connections between probabilistic linear solvers andprojection methods for linear systems, providing a probabilistic interpretation of afar more general class of iterative methods. In particular, this provides such aninterpretation of the generalised minimum residual method. A probabilistic view ofpreconditioning is also introduced. These developments unify the literature onprobabilistic linear solvers and provide foundational connections to the literatureon iterative solvers for linear systems.

Implementation and performance evaluation of a communication-avoiding GMRES method for stencil-based code on GPU cluster

In this study, a communication-avoiding generalized minimum residual method (CA-GMRES) is implemented on a hybrid CPU–GPU cluster, targeted for the performance acceleration of iterative linear system solver in the gyrokinetic toroidal five-dimensional Eulerian code (GT5D). In the GT5D, its sparse matrix-vector multiplication operation (SpMV) is performed as a 17-point stencil-based computation. The specialized part for the GT5D is only in the SpMV, and the other parts are usable also for other application program codes. In addition to the CA-GMRES, we implement and evaluate a modified variant of CA-GMRES (M-CA-GMRES) proposed in the previous study Idomura et al. (in: Proceedings of the 8th workshop on latest advances in scalable algorithms for large-scale systems (ScalA ’17), 2017. https://doi.org/10.1145/3148226.3148234) to reduce the amount of floating-point calculations. This study demonstrates that beneficial features of the CA-GMRES are in its minimum number of collective communications and its highly efficient calculations based on dense matrix–matrix operations. The performance evaluation is conducted on the Reedbush-L GPU cluster, which contains four NVIDIA Tesla P100 (Pascal GP100) GPUs per compute node. The evaluation results show that the M-CA-GMRES or CA-GMRES for the GT5D is advantageous over the GMRES or the generalized conjugate residual method (GCR) on GPU clusters, especially when the problem size (vector length) is large so that the cost of the SpMV is less dominant. The M-CA-GMRES is 1.09 ×, 1.22 × and 1.50 × faster than the CA-GMRES, GCR and GMRES, respectively, when 64 GPUs are used.

Survey and Review

The Survey and Review article in this issue is “A Class of Iterative Solvers for the Helmholtz Equation: Factorizations, Sweeping Preconditioners, Source Transfer, Single Layer Potentials, Polarized Traces, and Optimized Schwarz Methods,” by Martin J. Gander and Hui Zhang. The Helmholtz equation appears in the study of time-harmonic waves \(u(x) \exp(-iømega t)\) arising in acoustics, electromagnetism, and other fields. In its simplest form the equation reads \((\Delta+k^2) u = f\), where (\Delta\) denotes the Laplacian operator and the constant \(k\) is the wavenumber, and therefore is reminiscent of the Poisson equation (\Delta u =f\). However, this similarity is only superficial, and the numerical solution of boundary value problems, often posed in unbounded domains, involving the Helmholtz equation is challenging. Classical iterative solvers are not suitable, and the recent literature has devoted much effort to the development of novel solvers. New iterative methods have been suggested by authors with a wide range of backgrounds and motivations and invoke different discretization strategies; some of them have been formulated at the linear algebra level and others at the level of the differential equation. Gander and Zhang survey the new developments using a unifying point of view. The main insight is given in section 2: the forward and backward substitution sweeps of the block LU factorization of the discrete system correspond at the continuous level to sweeps from left to right and right to left of an optimal Schwarz domain decomposition method. The article will be useful not only to those working directly with the Helmholtz and related partial differential equations, but also to anybody interested in iterative solvers for large linear systems.

Algebraic multigrid methods for dual mortar finite element formulations in contact mechanics

SummaryThis paper proposes novel strategies to enable multigrid preconditioners within iterative solvers for linear systems arising from contact problems based on mortar finite element formulations. The so‐called dual mortar approach that is exclusively employed here allows for an easy condensation of the discrete Lagrange multipliers. Therefore, it has the advantage over standard mortar methods that linear systems with a saddle‐point structure are avoided, which generally require special preconditioning techniques. However, even with the dual mortar approach, the resulting linear systems turn out to be hard to solve using iterative linear solvers. A basic analysis of the mathematical properties of the linear operators reveals why the naive application of standard iterative solvers shows instabilities and provides new insights of how contact modeling affects the corresponding linear systems. This information is used to develop new strategies that make multigrid methods efficient preconditioners for the class of contact problems based on dual mortar methods. It is worth mentioning that these strategies primarily adapt the input of the multigrid preconditioners in a way that no contact‐specific enhancements are necessary in the multigrid algorithms. This makes the implementation comparably easy. With the proposed method, we are able to solve large contact problems, which is an important step toward the application of dual mortar–based contact formulations in the industry. Numerical results are presented illustrating the performance of the presented algebraic multigrid method.

Numerical Simulation of Shielding Current Density in HTS Film by Using Acceleration Technique With $\mathcal {H}$-Matrix Method

Numerical techniques are proposed for fast solving a linear system appearing in the shielding current analysis in a cracked high-temperature superconducting (HTS) film. After discretized with respect to time and space, an initial-boundary-value problem of the shielding current density reduces to a problem in which a linear system is solved at each iteration of the Newton method. Although the coefficient matrix of the system is dense, a fast iterative linear-system solver can be realized by implementing the H-matrix method to matrix-vector multiplications. On the other hand, the CPU time required for the iterative solver increases with the number of elements adjacent to cracks. To remove such unfavorable characteristics, the QR factorization is applied to a part of the linear system so that all unknowns other than correction of the current vector potential may be eliminated from the linear system. The performance of the techniques is evaluated and, as an application of the techniques, detectability of multiple cracks in an HTS film is investigated numerically.

A comparison of two preconditioner algorithms within the representer-based four-dimensional variational data assimilation system for the Navy coastal ocean model

ABSTRACTFour-dimensional variational (4D-Var) data assimilation for operational systems requires the solution of large linear systems that are poorly conditioned in general. In addition to efficient iterative solvers for linear systems, using a good preconditioner is required to guarantee an acceptable solution with a small number of iterations. We consider the assimilation of ocean observations using a weak constraints 4D-Var for the Navy Coastal Ocean Model (NCOM), based on the representer method. Two methods of preconditioning the linear system are implemented, namely the scaling of the linear system by the square root of the observations error variances, and the approximating stabilised representer matrix based on the computation of some representer functions. We evaluate their convergence using criteria such as the norm of the residuals and of the gradient of the cost function, the analysis error evaluated at the observation locations, and finally the convergence of the sequence of analyses. Results from all criteria show that the fastest convergence is achieved with the rescaling preconditioner when the conjugate gradient used to solve the linear system is equipped with a suitable inner product instead of the Euclidean.

Performance-Based Numerical Solver Selection in the Lighthouse Framework

Scientific and engineering computing rely heavily on linear algebra for large-scale data analysis, modeling and simulation, machine learning, and other applied problems. Sparse linear system solution often dominates the execution time of such applications, prompting the ongoing development of highly optimized iterative algorithms and high-performance parallel implementations. In the Lighthouse project, we enable application developers with varied backgrounds to readily discover and effectively apply the best available numerical software for their problems, aiming to maximize both developer productivity and application performance. Lighthouse is a search-based expert system built on a software taxonomy that combines expert knowledge, machine learning--based classification of existing numerical software collections, and automated code generation and optimization. In this paper we present the integration of PETSc and Trilinos iterative solvers for sparse linear systems into the Lighthouse framework. In addition to functional information in the taxonomy, we have created a comprehensive machine learning--based workflow for the automated classification of sparse solvers, which can be generalized to other types of rapidly evolving numerical methods. We present a comparative analysis of the solver classification results for a varied set of input problems and machine learning methods, achieving up to 93% accuracy in identifying the best-performing linear solution methods in PETSc and Trilinos.

Exploiting task and data parallelism in ILUPACK’s preconditioned CG solver on NUMA architectures and many-core accelerators

We present specialized implementations of the preconditioned iterative linear system solver in ILUPACK for Non-Uniform Memory Access (NUMA) platforms and many-core hardware co-processors based on the Intel Xeon Phi and graphics accelerators. For the conventional x86 architectures, our approach exploits task parallelism via the OmpSs runtime as well as a message-passing implementation based on MPI, respectively yielding a dynamic and static schedule of the work to the cores, with different numeric semantics to those of the sequential ILUPACK. For the graphics processor we exploit data parallelism by off-loading the computationally expensive kernels to the accelerator while keeping the numeric semantics of the sequential case.

Fast Iterative Solvers for Large Compressed-Sparse Row Linear Systems on Graphics Processing Unit

Engineering problems involve the solution of large sparse linear systems, and require therefore fast and high performance algorithms for algebra operations such as dot product, and matrix-vector multiplication. During the last decade, graphics processing units have been widely used. In this paper, linear algebra operations on graphics processing unit for single and double precision (with real and complex arithmetic) are analyzed in order to make iterative Krylov algorithms efficient compared to central processing units implementation. The performance of the proposed method is evaluated for the Laplace and the Helmholtz equations. Numerical experiments clearly show the robustness and effectiveness of the graphics processing unit tuned algorithms for compressed-sparse row data storage.

Fast Iterative Solvers for Large Compressed-Sparse Row Linear Systems on Graphics Processing Unit

Engineering problems involve the solution of large sparse linear systems, and require therefore fast and high performance algorithms for algebra operations such as dot product, and matrix-vector multiplication. During the last decade, graphics processing units have been widely used. In this paper, linear algebra operations on graphics processing unit for single and double precision (with real and complex arithmetic) are analyzed in order to make iterative Krylov algorithms efficient compared to central processing units implementation. The performance of the proposed method is evaluated for the Laplace and the Helmholtz equations. Numerical experiments clearly show the robustness and effectiveness of the graphics processing unit tuned algorithms for compressed-sparse row data storage.

3D Laplace-domain waveform inversion using a low-frequency time-domain modeling algorithm

We have developed a Laplace-domain full-waveform inversion technique based on a time-domain finite-difference modeling algorithm for efficient 3D inversions. Theoretically, the Laplace-domain Green’s function multiplied by a constant can be obtained regardless of the frequency content in the time-domain source wavelet. Therefore, we can use low-frequency sources and large grids for efficient modeling in the time domain. We Laplace-transform time-domain seismograms to the Laplace domain and calculate the residuals in the Laplace domain. Then, we back-propagate the Laplace-domain residuals in the time domain using a predefined time-domain source wavelet with the amplitude of the residuals. The back-propagated wavefields are transformed to the Laplace domain again to update the velocity model. The inversion results are long-wavelength velocity models on large grids similar to those obtained by the original approach based on Laplace-domain modeling. Inversion examples with 2D Gulf of Mexico field data revealed that the method yielded long-wavelength velocity models comparable with the results of the original Laplace-domain inversion methods. A 3D SEG/EAGE salt model example revealed that the 3D Laplace-domain inversion based on time-domain modeling method can be more efficient than the inversion based on Laplace-domain modeling using an iterative linear system solver.

Comparative study of conjugate gradient algorithms performance on the example of steady-state axisymmetric heat transfer problem

Abstract The finite element method (FEM) is one of the most frequently used numerical methods for finding the approximate discrete point solution of partial differential equations (PDE). In this method, linear or nonlinear systems of equations, comprised after numerical discretization, are solved to obtain the numerical solution of PDE. The conjugate gradient algorithms are efficient iterative solvers for the large sparse linear systems. In this paper the performance of different conjugate gradient algorithms: conjugate gradient algorithm (CG), biconjugate gradient algorithm (BICG), biconjugate gradient stabilized algorithm (BICGSTAB), conjugate gradient squared algorithm (CGS) and biconjugate gradient stabilized algorithm with l GMRES restarts (BICGSTAB(l)) is compared when solving the steady-state axisymmetric heat conduction problem. Different values of l parameter are studied. The engineering problem for which this comparison is made is the two-dimensional, axisymmetric heat conduction in a finned circular tube.