Preconditioners For Matrices Research Articles

Preconditioned discontinuous Galerkin method and convection‐diffusion‐reaction problems with guaranteed bounds to resulting spectra

AbstractThis paper focuses on the design, analysis and implementation of a new preconditioning concept for linear second order partial differential equations, including the convection‐diffusion‐reaction problems discretized by Galerkin or discontinuous Galerkin methods. We expand on the approach introduced by Gergelits et al. and adapt it to the more general settings, assuming that both the original and preconditioning matrices are composed of sparse matrices of very low ranks, representing local contributions to the global matrices. When applied to a symmetric problem, the method provides bounds to all individual eigenvalues of the preconditioned matrix. We show that this preconditioning strategy works not only for Galerkin discretization, but also for the discontinuous Galerkin discretization, where local contributions are associated with individual edges of the triangulation. In the case of nonsymmetric problems, the method yields guaranteed bounds to real and imaginary parts of the resulting eigenvalues. We include some numerical experiments illustrating the method and its implementation, showcasing its effectiveness for the two variants of discretized (convection‐)diffusion‐reaction problems.

Preconditioning Sparse Matrices with Alternating and Multiplicative Operator Splittings

We present an algebraic framework for operator splitting preconditioners for general sparse matrices. The framework leads to four different approaches: two with alternating splittings and two with a multiplicative ansatz. The ansatz generalizes ADI and ILU methods to multiple factors and to a more general factor form. The factors may be computed directly from the matrix coefficients or adaptively by incomplete sparse inversions. The special case of tridiagonal splittings is examined in more detail. We decompose the adjacency graph of the sparse matrix into multiple (almost) disjoint linear forests, and each linear forest (a union of disjoint paths) leads to a tridiagonal splitting. We obtain specialized variants of the four general approaches. Parallel implementations for all steps are provided on a GPU. We demonstrate the effectiveness and efficiency of these preconditioners combined with GMRES on various matrices.

A Novel Embedded Domain Decomposition Method for Electromagnetic Simulation of Structures in Inhomogeneous Medium

A novel embedded domain decomposition method (EDDM) is proposed to simulate electromagnetic scattering from structures in an inhomogeneous medium. In this method, the inhomogeneous medium (e.g., layered substrate) is set as the background subdomain, and the integrated structures (e.g., metallic/dielectric components) with properly defined buffer regions (BRs) are set as the embedded subdomains, where the meshes of two can be completely independent and arbitrarily overlapping. The major contribution of this work lies in the following three aspects. First, to allow for arbitrary inclusion of the embedded subdomains in an inhomogeneous background subdomain, an adaptive BR method is developed. The self-coupling equation is rederived to ensure the consistency of the BR and the inhomogeneous background, and a new mutual-coupling equation is introduced to account for the material and/or conductor differences. Second, the boundary element method (BEM) accelerated by the multilevel fast multipole algorithm (MLFMA) is utilized to truncate the finite-element region so that the computation capability can be enhanced and the extra air box setting of the background subdomain in the original EDDM can be avoided. Third, for structure adjustments such as translations and rotations, we further propose an inherited calculation to avoid mesh regeneration and repeated computation of the self-coupling and preconditioning matrices and also reduce the iteration counts and solution time. Numerical examples have shown the accuracy, robustness, and efficiency of the proposed method.

Varying Infimum Gradient Descent Algorithm for Agent-Sever Systems Using Different Order Iterative Preconditioning Methods

In the traditional gradient descent (T-GD) algorithm, the convergence rate is strongly depend on the condition number of the information matrix: a larger condition number leads to a poor optimal convergence factor infimum <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$\mu _{\text{op}}$</tex-math></inline-formula> , which sets a convergence rate ceiling. That is, once the information matrix is fixed, the convergence factor of the T-GD algorithm reaches at most the infimum <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$\mu _{\text{op}}$</tex-math></inline-formula> . This article studies a varying infimum gradient descent algorithm, which can move down the infimum by using different order iterative preconditioning methods, as follows: first, for infinite iterative algorithm, the infimum becomes smaller and smaller with the increased iteration numbers; second, for finite iterative algorithm, the infimum is equal to zero, and the parameter estimates can be obtained in only one iteration; third, construct an adaptive interval between zero and <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$\mu _{\text{op}}$</tex-math></inline-formula> , which can establish a link between the least squares and T-GD algorithms. Based on the varying infimum gradient descent algorithm, researchers can adaptively choose preconditioning matrices for different kinds of models on a case by case basis. The convergence analysis and simulation examples show effectiveness of the proposed algorithms.

A shooting approach to the scaled boundary finite element equations of elastodynamics in the frequency domain

Dealing with numerical analysis of problems, especially ones with semi-infinite boundaries, scaled boundary finite element method has emerged as one of the efficient tools for the task. Combining the exactness of strong forms with the flexibility of weak formulations makes the method an improvement to its predecessors. Problem with the method arises when the analytical solution of the semi-discretized system is not available, which is the case for numerous problems. In the most recent attempt to solve the issue, a shooting method was proposed for elastostatic problems. Generality of the method removes any concerns regarding the type of governing equations since it no longer needs any analytical solutions. In this study, the method is further developed to deal with elastodynamic problems of viscoelastic media in the frequency domain, for which the closed-form solution is not available. Adopting the shooting method, the two-point boundary value problem of each subdomain is treated as a second-order initial value problem. Boundary conditions of the bounded subdomains are formulated by enforcing equilibrium at both ends. On the other hand, in contrast to linear elasticity, neither forces nor displacements vanish at infinity. As a result, a special treatment is required to deal with the boundary conditions of the unbounded subdomains. A Runge–Kutta–Nyström method, which consists of an embedded pair of adjacent order, is utilized to integrate the strong form of oscillatory equations enabling local error estimation and adaptive step-size control. This guarantees efficiency since the integration proceeds not by equal step-sizes, but by equal user-defined local errors. In the first step, analytical solution of a viscoelastic half-space subjected to a harmonic strip load is used for the verification process. First, the half-space is discretized by Voronoi cells and each truncated edge is defined using an unbounded subdomain. Next, one bounded subdomain is used to define the whole bounded area, and one unbounded subdomain resembles the radiation condition of the far-field boundaries. For the cases involving only bounded subdomains, three other numerical tests, including a perforated plate, an L-shaped panel, and a bi-material specimen are employed, and the effect of discretization method on the results is analyzed. The method is a huge step forward since the size of subdomains are not restricted, far-field boundaries can be arbitrary-shaped, hysteretic damping can be naturally included, no additional unknowns are introduced, no compound matrices are defined, no preconditioning matrices are needed, and a simple unified approach is used for both bounded and unbounded subdomains.

Optimal rotated block-diagonal preconditioning for discretized optimal control problems constrained with fractional time-dependent diffusive equations

For a class of optimal control problems constrained with certain time- and space-fractional diffusive equations, by making use of mixed discretizations of temporal finite-difference and spatial finite-element schemes along with Lagrange multiplier approach, we obtain specially structured block two-by-two linear systems. We demonstrate positive definiteness of the coefficient matrices of these discrete linear systems, construct rotated block-diagonal preconditioning matrices, and analyze spectral properties of the corresponding preconditioned matrices. Both theoretical analysis and numerical experiments show that the preconditioned Krylov subspace iteration methods, when incorporated with these rotated block-diagonal preconditioners, can exhibit optimal convergence property in the sense that their convergence rates are independent of both discretization stepsizes and problem parameters, and their computational workloads are linearly proportional with the number of discrete unknowns.

An efficient preconditioned variant of the PSS preconditioner for generalized saddle point problems

For the generalized saddle point problems from the Stokes equations, we propose an efficient preconditioned variant of the positive semidefinite and skew-Hermitian (PVPSS) preconditioner. The new preconditioner is established by adopting matrix preconditioning strategy and relaxation technique for the PSS one derived by Pan et al. (2006). Compared with the PSS one, the PVPSS preconditioner is much closer to the coefficient matrix and easier to be implemented if proper preconditioning matrices are adopted. We prove the convergence of the PVPSS iteration method under some restrictions and discuss the spectral properties of the PVPSS preconditioned matrix. Meanwhile, the implementation and a practical way to choose the parameter of the PVPSS preconditioner are discussed. Comparisons between the PVPSS preconditioner and some existing ones are also given. Numerical experiments are carried out to illustrate that the proposed preconditioner is effective for the generalized saddle point problems and outperforms several other commonly used ones.

A note on parallel preconditioning for all-at-once evolutionary PDEs

McDonald, Pestana, and Wathen [SIAM J. Sci. Comput., 40 (2018), pp. A1012–A1033] present a method for preconditioning time-dependent PDEs via an approximation by a nearby time-periodic problem, that is, they employ circulant-related matrices as preconditioners for the non-symmetric block Toeplitz matrices which arise from an all-at-once formulation. They suggest that such an approach might be efficiently implemented in parallel. In this short article, we present parallel numerical results for their preconditioner which exhibit strong scaling. We also extend their preconditioner via a Neumann series approach which also allows for efficient parallel execution. Results are shown for both parabolic and hyperbolic PDEs. Our simple implementation (in C++ and MPI) is available at the Git repository PARALAAOMPI.

Condition number analysis of flow fields arising from CFD simulations

Condition numbers of pressure and velocity matrices obtained from an unstructured finite volume discretization of the three dimensional Navier–Stokes equations are studied for two distinct flow configurations. The exact condition number is calculated based on singular values. The continuous estimation of condition number of the velocity matrices with respect to the time varying velocity field is achieved through the Gershgorin theoretical bounds for the eigenvalues. The focus of the present study is on the analysis of flow dynamics with respect to the degree of ill-conditioning of the linear systems. For the two benchmark problems, condition numbers of pressure and velocity matrices are presented. The pressure matrix is shown to be more ill-conditioned than velocity and requires strong preconditioning. For a 3D lid-driven cavity, the appearance and disappearance of corner vortices at higher Reynolds number is clearly reflected in the condition number variation with time. In a second application of flow past a circular cylinder, condition number variation with time during the initial phase clearly reveals the onset of vortex shedding. The condition number based on Gershgorin bounds is also used to switch the SGS and ILU preconditioners for velocity matrices when they are well conditioned. Consequently, an overall reduction in the simulation time is demonstrated.

On Partial Cholesky Factorization and a Variant of Quasi-Newton Preconditioners for Symmetric Positive Definite Matrices

This work studies limited memory preconditioners for linear symmetric positive definite systems of equations. Connections are established between a partial Cholesky factorization from the literature and a variant of Quasi-Newton type preconditioners. Then, a strategy for enhancing the Quasi-Newton preconditioner via available information is proposed. Numerical experiments show the behaviour of the resulting preconditioner.

BFGS‐like updates of constraint preconditioners for sequences of KKT linear systems in quadratic programming

SummaryWe focus on efficient preconditioning techniques for sequences of Karush‐Kuhn‐Tucker (KKT) linear systems arising from the interior point (IP) solution of large convex quadratic programming problems. Constraint preconditioners (CPs), although very effective in accelerating Krylov methods in the solution of KKT systems, have a very high computational cost in some instances, because their factorization may be the most time‐consuming task at each IP iteration. We overcome this problem by computing the CP from scratch only at selected IP iterations and by updating the last computed CP at the remaining iterations, via suitable low‐rank modifications based on a BFGS‐like formula. This work extends the limited‐memory preconditioners (LMPs) for symmetric positive definite matrices proposed by Gratton, Sartenaer and Tshimanga in 2011, by exploiting specific features of KKT systems and CPs. We prove that the updated preconditioners still belong to the class of exact CPs, thus allowing the use of the conjugate gradient method. Furthermore, they have the property of increasing the number of unit eigenvalues of the preconditioned matrix as compared with the generally used CPs. Numerical experiments are reported, which show the effectiveness of our updating technique when the cost for the factorization of the CP is high.

Threshold incomplete factorization constraint preconditioners for saddle-point matrices

This paper presents a drop-threshold incomplete LD-1LT (δ) factorization constraint preconditioner for saddle-point systems using a threshold parameter δ. A transformed saddle-point matrix is partitioned into a block structure with blocks of order 1 and 2 constituting ‘a priori pivots’. Based on these pivots an incomplete LD-1LT (δ) factorization constraint preconditioner is computed that approaches an exact form as δ approaches zero. We prove that both the exact and incomplete factorizations exist such that the entries of the constraint block remain unaltered in the triangular factors. Numerical results are presented for validation.

Fast multipole preconditioners for sparse matrices arising from elliptic equations

Among optimal hierarchical algorithms for the computational solution of elliptic problems, the fast multipole method (FMM) stands out for its adaptability to emerging architectures, having high arithmetic intensity, tunable accuracy, and relaxable global synchronization requirements. We demonstrate that, beyond its traditional use as a solver in problems for which explicit free-space kernel representations are available, the FMM has applicability as a preconditioner in finite domain elliptic boundary value problems, by equipping it with boundary integral capability for satisfying conditions at finite boundaries and by wrapping it in a Krylov method for extensibility to more general operators. Here, we do not discuss the well developed applications of FMM to implement matrix-vector multiplications within Krylov solvers of boundary element methods. Instead, we propose using FMM for the volume-to-volume contribution of inhomogeneous Poisson-like problems, where the boundary integral is a small part of the overall computation. Our method may be used to precondition sparse matrices arising from finite difference/element discretizations, and can handle a broader range of scientific applications. It is capable of algebraic convergence rates down to the truncation error of the discretized PDE comparable to those of multigrid methods, and it offers potentially superior multicore and distributed memory scalability properties on commodity architecture supercomputers. Compared with other methods exploiting the low-rank character of off-diagonal blocks of the dense resolvent operator, FMM-preconditioned Krylov iteration may reduce the amount of communication because it is matrix-free and exploits the tree structure of FMM. We describe our tests in reproducible detail with freely available codes and outline directions for further extensibility.

Matrix-free monolithic homotopy continuation with application to computational aerodynamics

A matrix-free monolithic homotopy continuation algorithm is developed which allows for approximate numerical solutions to nonlinear systems of equations without the need to solve a linear system, thereby avoiding the formation of any Jacobian or preconditioner matrices. The algorithm can converge from an arbitrary starting guess, under suitable conditions, and can give a sufficiently accurate approximation to the converged solution such that a rapid locally convergent method such as Newton’s method will converge successfully. Several forms of the algorithm are presented, as are augmentations to the algorithms which can lead to improved efficiency or stability. The method is validated and the stability and efficiency are investigated numerically based on a computational aerodynamics flow solver.

Preconditioning Techniques Based on the Birkhoff–von Neumann Decomposition

Abstract We introduce a class of preconditioners for general sparse matrices based on the Birkhoff–von Neumann decomposition of doubly stochastic matrices. These preconditioners are aimed primarily at solving challenging linear systems with highly unstructured and indefinite coefficient matrices. We present some theoretical results and numerical experiments on linear systems from a variety of applications.

Preconditioners for hierarchical matrices based on their extended sparse form

AbstractIn this paper we consider linear systems with dense-matrices which arise from numerical solution of boundary integral equations. Such matrices can be well-approximated with ℋ

MPBEC, a Matlab Program for Biomolecular Electrostatic Calculations

One of the most used and efficient approaches to compute electrostatic properties of biological systems is to numerically solve the Poisson–Boltzmann (PB) equation. There are several software packages available that solve the PB equation for molecules in aqueous electrolyte solutions. Most of these software packages are useful for scientists with specialized training and expertise in computational biophysics. However, the user is usually required to manually take several important choices, depending on the complexity of the biological system, to successfully obtain the numerical solution of the PB equation. This may become an obstacle for researchers, experimentalists, even students with no special training in computational methodologies. Aiming to overcome this limitation, in this article we present MPBEC, a free, cross-platform, open-source software that provides non-experts in the field an easy and efficient way to perform biomolecular electrostatic calculations on single processor computers. MPBEC is a Matlab script based on the Adaptative Poisson–Boltzmann Solver, one of the most popular approaches used to solve the PB equation. MPBEC does not require any user programming, text editing or extensive statistical skills, and comes with detailed user-guide documentation. As a unique feature, MPBEC includes a useful graphical user interface (GUI) application which helps and guides users to configure and setup the optimal parameters and approximations to successfully perform the required biomolecular electrostatic calculations. The GUI also incorporates visualization tools to facilitate users pre- and post-analysis of structural and electrical properties of biomolecules. Program summaryProgram title: MPBECCatalogue identifier: AEXW_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEXW_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: UTSA licenseNo. of lines in distributed program, including test data, etc.: 1505214No. of bytes in distributed program, including test data, etc.: 115585507Distribution format: tar.gzProgramming language: Matlab R2010a and higher versions.Computer: Any computer at least with 2.5 GHz speed.Operating system: Windows, Linux Fedora, and Mac OSX.RAM: At least 2 GB of free memory.Classification: 3, 10.External routines: Java, Jmol, pdb2pqr, Propka.Nature of problem: Numerical solution of the linearized PB equation.Solution method: MPBEC uses sparse matrix calculations and iterative linear solvers (gmres, minres, and bicgstab) to solve the discretized (box method) linear PB equation. The stability and the convergence of these iterative solvers are improved by using the incomplete LU factorization which generates pre-conditioner matrices L (lower triangular matrix) and U (upper triangular matrix) at low computational cost.Unusual features: MPBEC incorporates a GUI to provide non-experts in computational biophysics a user-friendly and intuitive tool to obtain biomolecular electrostatic calculations. The GUI provides helpful information about how to fill out the input data by moving the mouse pointer over the corresponding text or blank box. The GUI tests all the input data before running MPBEC to make sure that the solver is properly configured and facilitates the user with the pre- and post-analysis of the biomolecular calculations.Additional comments: MPBEC checks available and required RAM memory before allocation to alert users when the available RAM memory is insufficient to perform the required calculations depending basically on the number of grid points and the size and number of atoms of the biomolecule(s). Therefore, MPBEC might be used with lower RAM memory and processor speed conditions than those mentioned above. The distribution file for this program is over 115 Mbytes and therefore is not delivered directly when download or Email is requested. Instead a html file giving details of how the program can be obtained is sent.Running time: The software runs on single processor computers at low-to-moderate computational cost depending on the computer performance, the biomolecule size, the grid resolution, and required calculations.

On Block Preconditioners for Pde-Constrained Optimization Problems

Recently, Bai proposed a block-counter-diagonal and a block-counter-triangular preconditioning matrices to precondition the GMRES method for solving the structured system of linear equations arising from the Galerkin finite-element discretizations of the distributed control problems in (Computing 91 (2011) 379-395). He analyzed the spectral properties and derived explicit expressions of the eigenvalues and eigenvectors of the preconditioned matrices. By applying the special structures and properties of the eigenvector matrices of the preconditioned matrices, we derive upper bounds for the 2-norm condition numbers of the eigenvector matrices and give asymptotic convergence factors of the preconditioned GMRES methods with the block-counter-diagonal and the block-counter-triangular preconditioners. Experimental results show that the convergence analyses match well with the numerical results.

Preconditioning Techniques for Diagonal-times-Toeplitz Matrices in Fractional Diffusion Equations

The fractional diffusion equation is discretized by an implicit finite difference scheme with the shifted Grunwald formula, which is unconditionally stable. The coefficient matrix of the discretized linear system is equal to the sum of a scaled identity matrix and two diagonal-times-Toeplitz matrices. Standard circulant preconditioners may not work for such Toeplitz-like linear systems. The main aim of this paper is to propose and develop approximate inverse preconditioners for such Toeplitz-like matrices. An approximate inverse preconditioner is constructed to approximate the inverses of weighted Toeplitz matrices by circulant matrices, and then combine them together row-by-row. Because of Toeplitz structure, both the discretized coefficient matrix and the preconditioner can be implemented very efficiently by using fast Fourier transforms. Theoretically, we show that the spectra of the resulting preconditioned matrices are clustered around one. Thus Krylov subspace methods with the proposed preconditione...

Block-triangular preconditioning methods for linear third-order ordinary differential equations based on reduced-order sinc discretizations

By applying the reduced-order sinc discretization to the two-point boundary value problem of a linear third-order ordinary differential equation, we can obtain a block two-by-two system of linear equations, with each block of its coefficient matrix being a combination of Toeplitz and diagonal matrices. This class of linear systems can be effectively solved by Krylov subspace iteration methods such as GMRES and BiCGSTAB. We construct block-triangular preconditioning matrices to accelerate the convergence rates of the Krylov subspace iteration methods, and demonstrate that the eigenvalues of certain approximations to the preconditioned matrices are uniformly bounded within a rectangle, being independent of the size of the discrete linear system, on the complex plane. In addition, we use numerical examples to show the effectiveness of the proposed preconditioning methods.