Indefinite Systems Of Linear Equations Research Articles

An efficient algorithm for the minimal least squares solution of linear systems with indefinite symmetric matrices

In this work, a new algorithm for solving symmetric indefinite systems of linear equations is presented. It factorizes the matrix into the form LDLt using Jacobi rotations in order to increase the pivot´s absolute value. Furthermore, Rook´s pivoting strategy is also adapted and implemented. In determinate compatible systems, the computational cost of the algorithm was similar to the cost of the Bunch-Kaufman method, but the error was approximately 50 % smaller for intermediate and large matrices, regardless of the condition number of the coefficient matrix. Furthermore, unlike Bunch-Kaufman, the new algorithm calculates with little additional cost the fundamental basis of the null space, and obtains the minimal least squares and minimum norm solutions. In minimal least squares with minimum norm problems, the new algorithm was compared with the LAPACK Complete Orthogonal Decomposition algorithm, among others. The obtained error with both algorithms was similar but the computational cost was at least 20 % smaller with the new algorithm, even though the Complete Orthogonal Decomposition is implemented in a blocked form.

Uzawa methods for a class of block three‐by‐three saddle‐point problems

SummaryIn this work, we consider numerical methods for solving a class of block three‐by‐three saddle‐point problems, which arise from finite element methods for solving time‐dependent Maxwell equations and some other applications. The direct extension of the Uzawa method for solving this block three‐by‐three saddle‐point problem requires the exact solution of a symmetric indefinite system of linear equations at each step. To avoid heavy computations at each step, we propose an inexact Uzawa method, which solves the symmetric indefinite linear system in some inexact way. Under suitable assumptions, we show that the inexact Uzawa method converges to the unique solution of the saddle‐point problem within the approximation level. Two special algorithms are customized for the inexact Uzawa method combining the splitting iteration method and a preconditioning technique, respectively. Numerical experiments are presented, which demonstrated the usefulness of the inexact Uzawa method and the two customized algorithms.

Preordering saddle‐point systems for sparse LDLT factorization without pivoting

SummaryThis paper focuses on efficiently solving large sparse symmetric indefinite systems of linear equations in saddle‐point form using a fill‐reducing ordering technique with a direct solver. Row and column permutations partition the saddle‐point matrix into a block structure constituting a priori pivots of order 1 and 2. The partitioned matrix is compressed by treating each nonzero block as a single entry, and a fill‐reducing ordering is applied to the corresponding compressed graph. It is shown that, provided the saddle‐point matrix satisfies certain criteria, a block LDLT factorization can be computed using the resulting pivot sequence without modification. Numerical results for a range of problems from practical applications using a modern sparse direct solver are presented to illustrate the effectiveness of the approach.

Improving the stability and robustness of incomplete symmetric indefinite factorization preconditioners

SummarySparse symmetric indefinite linear systems of equations arise in numerous practical applications. In many situations, an iterative method is the method of choice but a preconditioner is normally required for it to be effective. In this paper, the focus is on a class of incomplete factorization algorithms that can be used to compute preconditioners for symmetric indefinite systems. A limited memory approach is employed that incorporates a number of new ideas with the goal of improving the stability, robustness, and efficiency of the preconditioner. These include the monitoring of stability as the factorization proceeds and the incorporation of pivot modifications when potential instability is observed. Numerical experiments involving test problems arising from a range of real‐world applications demonstrate the effectiveness of our approach.

Solving dense symmetric indefinite systems using GPUs

SummaryThis paper studies the performance of different algorithms for solving a dense symmetric indefinite linear system of equations on multicore CPUs with a Graphics Processing Unit (GPU). To ensure the numerical stability of the factorization, pivoting is required. Obtaining high performance of such algorithms on the GPU is difficult because all the existing pivoting strategies lead to frequent synchronizations and irregular data accesses. Until recently, there has not been any implementation of these algorithms on a hybrid CPU/GPU architecture. To improve their performance on the hybrid architecture, we explore different techniques to reduce the expensive data transfer and synchronization between the CPU and GPU, or on the GPU (e.g., factorizing the matrix entirely on the GPU or in a communication‐avoiding fashion). We also study the performance of the solver using iterative refinements along with the factorization without pivoting combined with the preprocessing technique based on random butterfly transformations, or with the mixed‐precision algorithm where the matrix is factorized in single precision. This randomization algorithm only has a probabilistic proof on the numerical stability, and for this paper, we only focused on the mixed‐precision algorithm without pivoting. However, they demonstrate that we can obtain good performance on the GPU by avoiding the pivoting and using the lower precision arithmetics, respectively. As illustrated with the application in acoustics studied in this paper, in many practical cases, the matrices can be factorized without pivoting. Because the componentwise backward error computed in the iterative refinement signals when the algorithm failed to obtain the desired accuracy, the user can use these potentially unstable but efficient algorithms in most of the cases and fall back to a more stable algorithm with pivoting only in the case of the failure. Copyright © 2017 John Wiley & Sons, Ltd.

More on Generalizations and Modifications of Iterative Methods for Solving Large Sparse Indefinite Linear Systems

Continuing from the works of Li et al. (2014), Li (2007), and Kincaid et al. (2000), we present more generalizations and modifications of iterative methods for solving large sparse symmetric and nonsymmetric indefinite systems of linear equations. We discuss a variety of iterative methods such as GMRES, MGMRES, MINRES, LQ-MINRES, QR MINRES, MMINRES, MGRES, and others.

New updates of incomplete LU factorizations and applications to large nonlinear systems

In this paper, we address the problem of preconditioning sequences of large sparse indefinite systems of linear equations and present two new strategies to construct approximate updates of factorized preconditioners. Both updates are based on the availability of an incomplete factorization for one matrix of the sequence and differ in the approximation of the so-called ideal update. For a general treatment, an incomplete LU (ILU) factorization is considered, but the proposed approaches apply to incomplete factorizations of symmetric matrices as well. The first strategy is an approximate diagonal update of the ILU factorization; the second strategy relies on banded approximations of the factors in the ideal update. The efficiency and reliability of the proposed preconditioners are shown in the solution of nonlinear systems of equations by preconditioned Newton–Krylov methods. Nearly matrix-free implementations of the updating strategy are provided, and numerical experiments are carried out on application problems.

An improved two‐grid preconditioner for the solution of three‐dimensional Helmholtz problems in heterogeneous media

SUMMARYIn this paper, we address the solution of three‐dimensional heterogeneous Helmholtz problems discretized with second‐order finite difference methods with application to acoustic waveform inversion in geophysics. In this setting, the numerical simulation of wave propagation phenomena requires the approximate solution of possibly very large indefinite linear systems of equations. For that purpose, we propose and analyse an iterative two‐grid method acting on the Helmholtz operator where the coarse grid problem is solved inaccurately. A cycle of a multigrid method applied to a complex shifted Laplacian operator is used as a preconditioner for the approximate solution of this coarse problem. A single cycle of the new method is then used as a variable preconditioner of a flexible Krylov subspace method. We analyse the properties of the resulting preconditioned operator by Fourier analysis. Numerical results demonstrate the effectiveness of the algorithm on three‐dimensional applications. The proposed numerical method allows us to solve three‐dimensional wave propagation problems even at high frequencies on a reasonable number of cores of a distributed memory computer. Copyright © 2012 John Wiley & Sons, Ltd.

Performance of Zero-Level Fill-In Preconditioning Techniques for Iterative Solutions with Geotechnical Applications

AbstractBiot's symmetric indefinite linear systems of equations are commonly encountered in finite-element computations of geotechnical problems. The development of efficient solution methods for Biot's linear systems of equations is of practical importance to geotechnical software packages. In conjunction with the Krylov-subspace iterative method symmetric quasi-minimal residual (SQMR), some zero-level fill-in incomplete factorization preconditioning techniques including a symmetric successive overrelaxation (SSOR) type method and several zero-level incomplete LU [ILU(0)] methods are investigated and compared for Biot's symmetric indefinite linear systems of equations. Numerical experiments are carried out based on three practical geotechnical problems. Numerical results indicate that ILU(0) preconditioners are classical and generally efficient when adequately stabilized. However, the tunnel problem provides a counterexample demonstrating that ILU(0) preconditioners cannot be fully stabilized by prelimin...

On a new 3D primal-mixed finite element approach for thermal stress analysis of multi-layered geometrically multiscale structures

This paper presents a new reliable fully three-dimensional time efficient primal-mixed finite element approach with continuous primal and dual variables in geometrically multiscale thermoelasticity. The semi-coupling between thermal and mechanical physical fields is achieved straightforwardly via essential boundary condition per stress, and without consistency error. The direct sparse solver and matrix scaling routine are used for the solution of resulting large scale indefinite systems of linear equations. It will be shown that present solid finite element HC8/27 passes the first and the second stability condition (inf-sup test) for highly distorted finite elements with aspect ratio up to 7 orders of magnitude, for both, compressible and nearly incompressible materials. A number of pathological benchmark model problems, with material interfaces or coatings, with geometrical scale resolutions up to 8 orders of magnitude and aspect ratio of finite elements up to 7 orders of magnitude, are examined to test the robustness and execution times. It is shown that by rapid varying of spatial scale over local heterogeneities, the singularity of stress is captured without oscillation. It is shown that, if needed, present approach can simulate the simplified theories, as beam and plate theories, if the same restrictions on the stress tensor components are imposed. The new definition of multiscale reliability is given.

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.

Mixed finite element methods for incompressible flow: Stationary Stokes equations

AbstractIn this article, we develop and analyze a mixed finite element method for the Stokes equations. Our mixed method is based on the pseudostress‐velocity formulation. The pseudostress is approximated by the Raviart‐Thomas (RT) element of index k ≥ 0 and the velocity by piecewise discontinuous polynomials of degree k. It is shown that this pair of finite elements is stable and yields quasi‐optimal accuracy. The indefinite system of linear equations resulting from the discretization is decoupled by the penalty method. The penalized pseudostress system is solved by the H(div) type of multigrid method and the velocity is then calculated explicitly. Alternative preconditioning approaches that do not involve penalizing the system are also discussed. Finally, numerical experiments are presented. © 2009 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2010

Convergence analysis of a balancing domain decomposition method for solving a class of indefinite linear systems

AbstractA variant of balancing domain decomposition method by constraints (BDDC) is proposed for solving a class of indefinite systems of linear equations of the form (K−σ2M)u=f, which arise from solving eigenvalue problems when an inverse shifted method is used and also from the finite element discretization of Helmholtz equations. Here, both K and M are symmetric positive definite. The proposed BDDC method is closely related to the previous dual–primal finite element tearing and interconnecting method (FETI‐DP) for solving this type of problems (Appl. Numer. Math. 2005; 54:150–166), where a coarse level problem containing certain free‐space solutions of the inherent homogeneous partial differential equation is used in the algorithm to accelerate the convergence. Under the condition that the diameters of the subdomains are small enough, the convergence rate of the proposed algorithm is established, which depends polylogarithmically on the dimension of the individual subdomain problems and which improves with a decrease of the subdomain diameters. These results are supported by numerical experiments of solving a two‐dimensional problem. Copyright © 2009 John Wiley & Sons, Ltd.

Fast solvers with block‐diagonal preconditioners for linear FEM–BEM coupling

AbstractThe purpose of this paper is to present optimal preconditioned iterative methods to solve indefinite linear systems of equations arising from symmetric coupling of finite elements and boundary elements. This is a block‐diagonal preconditioner together with a conjugate residual method and a preconditioned inner–outer iteration. We prove the efficiency of these methods by showing that the number of iterations to preserve a given accuracy is bounded independent of the number of unknowns. Numerical examples underline the efficiency of these methods. Copyright © 2008 John Wiley & Sons, Ltd.

Sampling from Gaussian Markov random fields conditioned on linear constraints

Gaussian Markov random fields (GMRFs) are important modeling tools in statistics. They are often utilised to model spatially structured uncertainty, seasonal variation and other trends in the data. These last two examples of GMRFs are part of a larger class of GMRFs conditioned on linear constraints. Performing Monte Carlo Markov Chain inference on these models requires a large number of samples from GMRFs conditioned on linear constraints. Therefore it is vital to have fast and efficient methods for performing these samples. This article presents three Krylov subspace methods for sampling from a GMRF conditioned on linear constraints based on solving a Karush--Kuhn--Tucker, or saddle point, system. References H. Rue and L. Held. Gaussian Markov Random Fields: Theory and Applications . Chapman and Hall/CRC, 2005. USA. H. Rue. {Fast Sampling of {G}aussian {M}arkov Random Fields}. J. R. Statist. Soc. B , 63:325--338, 2001. doi:10.1111/1467-9868.00288 N. Cressie. Reply to {Wahba}. The American Statistician , 44(3):256--258, 1990. M. Benzi, G. Golub, and J. Liesen. Numerical solutions of saddle point problems. Acta Numerica , 14:1--137, 2005. doi:10.1017/S0962492904000212 D. P. Simpson, I. W. Turner, and A. N. Pettitt. Fast sampling from Gaussian Markov random field using Krylov subspace approaches, Scandinavian Journal of Statistics , Submitted, 2008. Y. Saad. Iterative Methods for Sparse Linear Systems . PWS Publishing Company, Boston, 1993. C. C. Paige and M. A. Saunders. Solution of sparse indefinite systems of linear equations. SIAM J. Numerical Analysis , 12:617--629, 1975. doi:10.1137/0712047 K. Burrage and J. Erhel. On the performance of various adaptive preconditioned {GMRES} strategies. Numerical Linear Algebra with Applications , 5(2):101--121, 1998. 3.0.CO;2-1>doi:0.1002/(SICI)1099-1506(199803/04)5:2 3.0.CO;2-1 A. N. Pettitt, I. S. Weir, and A. G. Hart. A conditional autoregressive {Gaussian} process for irregularly spaced multivariate data with application to modelling large sets of binary data. Statistics and Computing , 12, 2002. doi:10.1023/A:1020792130229

SIGEST

In discussions of the tremendous advances in scientific computation over the past few decades, it has often been written and demonstrated that as large as the increases in computer speed have been, they have been matched by improvements from more efficient algorithms. For many problem classes, over the last two to four decades, each has led to speedups of many, many orders of magnitude. The SIGEST paper in this issue, from the SIAM Journal on Scientific Computing, is an example of this tremendous type of algorithmic improvement. It also is an excellent example of how increasingly significant algorithmic advances are obtained by drawing upon a broad range of sophisticated techniques. In “On Large-Scale Diagonalization Techniques for the Anderson Model of Localization” by O. Schenk, M. Bollhöfer, and R. Römer, the authors demonstrate a greatly improved computational approach for an important, common problem: computing a few interior eigenvalues and eigenvectors of large, sparse, real symmetric, and highly indefinite matrices. By “large” they mean matrices with millions to hundreds of millions of rows and columns. The computational improvements that are demonstrated extend to a factor of a thousand or more over previous approaches, while keeping memory requirements modest. The application class is the Anderson model of localization, an important problem in quantum physics with a wide range of applications, and it appears likely that the approach can be applied successfully to large symmetric indefinite eigenvalue problems more generally. An interesting aspect of the paper is that is draws heavily upon recent advances in solving large sparse symmetric indefinite systems of linear equations. Solution of such systems is central to each of the eigenvalue solvers that the paper considers, two Lanczos-based methods and the Jacobi–Davidson method. In considering how to solve the linear systems in a manner that exploits sparsity as well as symmetry, the paper combines a diverse range of topics that includes preprocessing the matrix with symmetric weighted matchings, solving the linear system with Krylov subspace methods, and accelerating the linear system solution with multilevel preconditioners based upon incomplete factorizations. This makes for a very fitting paper for SIAM Review, highlighting a variety of important topics in numerical computation. We expect that a wide range of SIAM readers will find this to be a very readable paper that provides an interesting glimpse into a range of modern research approaches.

Preconditioners for solving the indefinite linear system for higher order H (curl) discretisations

AbstractThis paper presents numerical results for a preconditioning technique for the iterative solution of the indefinite linear system of equations obtained when a higher order H (curl) conforming finite element approximation is adopted for the vector wave equation. Numerical results are presented for a model problem that includes a singularity at the origin and dicretisations that include graded meshes and uniform p. (© 2008 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Construction of a discrete divergence-free basis through orthogonal factorization in $${\mathcal{H}}$$ -arithmetic

In computational fluid dynamics, linear constraints on the fluid velocity lead to challenging indefinite linear systems of equations. In this paper, we propose to compute an approximation to the constrained linear space of divergence-free functions using hierarchical matrix techniques. This approach will yield a data-sparse, well-conditioned basis of the desired subspace in almost optimal computational complexity which is confirmed by numerical tests. The novelty of this paper lies in the application of hierarchical matrix techniques to orthogonal factorization as well as the construction of an explicit approximation to the subspace basis.

Lanczos Conjugate-Gradient Method and Pseudoinverse Computation on Indefinite and Singular Systems

This paper extends some theoretical properties of the conjugate gradient-type method FLR (Ref. 1) for iteratively solving indefinite linear systems of equations. The latter algorithm is a generalization of the conjugate gradient method by Hestenes and Stiefel (CG, Ref. 2). We develop a complete relationship between the FLR algorithm and the Lanczos process, in the case of indefinite and possibly singular matrices. Then, we develop simple theoretical results for the FLR algorithm in order to construct an approximation of the Moore-Penrose pseudoinverse of an indefinite matrix. Our approach supplies the theoretical framework for applications within unconstrained optimization.

Regularized conjugate gradient method for skew-symmetric indefinite system of linear equations and applications

A class of regularized conjugate gradient methods is presented for solving large-scale sparse system of linear equations of which the coefficient matrix is an ill-conditioned skew-symmetric indefinite matrix. The convergence is proved and the possible choices of the parameters involved in the new methods are discussed in detail. Preliminary numerical computations show that the numerical behaviors of the new methods are superior to those of some standard Krylov subspace methods, such as CGNE, CGS, GMRES etc.