Sparse Linear Systems Research Articles

A meshless multi-symplectic local radial basis function collocation scheme for the “good” Boussinesq equation

A novel meshless multi-symplectic scheme is proposed for the “good” Boussinesq equation. The scheme consists of a local radial basis function (RBF) collocation method (LRBFCS) in space and a symplectic integrator in time. The LRBFCS is simple and efficient, since only a sparse banded linear system has to be solved. Moreover, it can avoid the ill-conditioned problem and shape-parameter-sensitivity of the global RBF method. The multi-symplectic LRBFCS (MLRBFCS) is more accurate than traditional methods. Numerical experiments with uniform knots and random knots are designed to verify the effectiveness of the method.

Algorithm 1021: SPEX Left LU, Exactly Solving Sparse Linear Systems via a Sparse Left-looking Integer-preserving LU Factorization

SPEX Left LU is a software package for exactly solving unsymmetric sparse linear systems. As a component of the sparse exact (SPEX) software package, SPEX Left LU can be applied to any input matrix, A , whose entries are integral, rational, or decimal, and provides a solution to the system \( Ax = b \) , which is either exact or accurate to user-specified precision. SPEX Left LU preorders the matrix A with a user-specified fill-reducing ordering and computes a left-looking LU factorization with the special property that each operation used to compute the L and U matrices is integral. Notable additional applications of this package include benchmarking the stability and accuracy of state-of-the-art linear solvers and determining whether singular-to-double-precision matrices are indeed singular. Computationally, this article evaluates the impact of several novel pivoting schemes in exact arithmetic, benchmarks the exact iterative solvers within Linbox, and benchmarks the accuracy of MATLAB sparse backslash. Most importantly, it is shown that SPEX Left LU outperforms the exact iterative solvers in run time on easy instances and in stability as the iterative solver fails on a sizeable subset of the tested (both easy and hard) instances. The SPEX Left LU package is written in ANSI C, comes with a MATLAB interface, and is distributed via GitHub, as a component of the SPEX software package, and as a component of SuiteSparse.

A local domain BEM for solving transient convection-diffusion-reaction problems

The convection-diffusion-reaction (CDR) equation has been extensively used to simulate a variety of physical phenomena. A robust numerical method for solving linear CDR problems is the Boundary Element Method (ΒΕΜ). However, the conventional BEM leads to dense coefficient matrices and as a result the memory requirements grow quadratically with respect to the number of degrees of freedom. In this work, a Local Domain BEM (LD-BEM) for solving the transient CDR equation with a constant velocity field is presented. The domain of interest is fragmented into small subdomains and the integral representation of the solution is considered separately for each of the subdomains. Eliminating the fluxes at all subdomain interfaces, the proposed LD-BEM leads to sparse linear system coefficient matrices and a reduced number of degrees of freedom. Eight numerical examples are solved to assess the efficiency and accuracy of the proposed method.

A novel localized collocation solver based on a radial Trefftz basis for thermal conduction analysis in FGMs with exponential variations

This paper introduces a novel localized radial Trefftz collocation method (LRTCM) for steady-state heat conduction analysis in 2D and 3D functionally graded materials (FGMs) with exponential variations. The LRTCM is constructed by combining the traditional radial Trefftz collocation method (RTCM) and the localized collocation scheme. In the proposed LRTCM, the radial Trefftz basis of the heat conduction equation in FGMs is derived by the variable transformation and domain mapping method, and the approximation solution at each node is expressed as a linear combination of a radial Trefftz basis in each local subdomain, thus the resulting sparse linear system naturally satisfies governing equation at interior nodes and boundary conditions at boundary nodes. The LRTCM not only inherits the advantages of the RTCM (meshlessness, no fictitious boundary and high accuracy), but also works well in large-scale structures or with more complicated boundary conditions when compared to the RTCM. Finally, the feasibility, accuracy, efficiency and stability of the proposed method are investigated by three examples.

Kaczmarz Anomaly in Tomography Problems

The Kaczmarz method is an important tool for solving large sparse linear systems that arise in computerized tomography. The Kaczmarz anomaly phenomenon has been observed recently when solving certain types of random systems. This raises the question of whether a similar anomaly occurs in tomography problems. The aim of the paper is to answer this question, to examine the extent of the phenomenon and to explain its reasons. Another tested issue is the ability of random row shuffles to sharpen the anomaly and to accelerate the rate of convergence. The results add important insight into the nature of the Kaczmarz method.

An efficient cascadic multigrid solver for 3-D magnetotelluric forward modelling problems using potentials

SUMMARY The fast and accurate 3-D magnetotelluric (MT) forward modelling is core engine of the interpretation and inversion of MT data. In this study, we develop an improved extrapolation cascadic multigrid method (EXCMG) to solve the large sparse complex linear system arising from the finite-element (FE) discretization on non-uniform orthogonal grids of the Maxwell’s equations using potentials. First, the vector Helmholtz equation and the scalar auxiliary equation are derived from the Maxwell’s equations using Coulomb-gauged potentials. The weighted residual method is adopted to discretize the weak formulation and assemble the FE equation. Secondly, carefully choosing the preconditioned complex stable bi-conjugate gradient method (BiCGStab) as multigrid smoother, we develop an improved EXCMG method on non-uniform grids to solve the resulting large sparse complex non-Hermitian linear systems. Finally, several examples including three standard testing models (COMMEMI3D-1, COMMEMI3D-2 and DTM1.0) and a topographic model are used to validate the accuracy and efficiency of the proposed multigrid solver. Numerical results show that the proposed EXCMG algorithm greatly improves the efficiency of 3-D MT forward modelling, is more efficient than some existing solvers, such as Pardiso, incomplete LU factorization preconditioned biconjugate gradients stabilized method (ILU-BiCGStab) and flexible generalized minimum residual method with auxiliary space Maxwell preconditioner (FGMRES-AMS), and capable to simulate large-scale problems with more than 100 million unknowns.

A localized spatiotemporal particle collocation method for long-time transient homogeneous diffusion analysis

This short communication proposes a novel localized spatiotemporal particle collocation method (LSPCM) to analyze long-time transient homogeneous diffusion problems. A time-dependent general solution of diffusion equation without singularity or near singularity at origin is developed. The LSPCM bypasses Laplace transformation or finite difference method to treat the time derivative of the diffusion equation, and can directly approximate the transient homogeneous diffusion problems. Each discrete node in computational domain is only supported by its nearby nodes covered by the local spatiotemporal subdomain. The physical variable of each discrete node is approximated by the linear combination of the physical variable of its nearby nodes and the time-dependent general solutions. Through satisfying boundary conditions and continuity conditions, the LSPCM eventually results in a large-scale sparse linear system which is easy to store and solve. Numerical experiments show that the LSPCM can accurately simulate long-time transient homogeneous diffusion problems by using a very long-time interval.

Random geometries for optimal control PDE problems based on fictitious domain FEMs and cut elements

This work investigates an elliptic optimal control problem defined on uncertain domains and discretized by a fictitious domain finite element method and cut elements. Key ingredients of the study are to manage cases considering the usually computationally “forbidden” combination of poorly conditioned equation system matrices due to challenging geometries, optimal control searches with iterative methods, slow convergence to system solutions on deterministic and non-deterministic level, and expensive remeshing due to geometrical changes. We overcome all these difficulties, utilizing the advantages of proper preconditioners adapted to unfitted mesh methods, improved types of Monte Carlo methods, and mainly employing the advantages of embedded FEMs, based on a fixed background mesh computed once even if geometrical changes are taking place. The sensitivity of the control problem is introduced in terms of random domains, employing a Quasi Monte Carlo method and emphasizing on a deterministic target state. The variational discretization concept is adopted, optimal error estimates for the state, adjoint state and control are derived that confirm the efficiency of the cut finite element method in challenging geometries. The performance of a multigrid scheme especially developed for unfitted finite element discretizations adapted to the optimal control problem is also tested. Some fundamental preconditioners are applied to the arising sparse linear systems coming from the discretization of the state and adjoint state variational forms in the spatial domain. The corresponding convergence rates along with the quality of the prescribed preconditioners are verified by numerical examples.

Hierarchical Orthogonal Factorization: Sparse Least Squares Problems

In this work, we develop a fast hierarchical solver for solving large, sparse least squares problems. We build upon the algorithm, spaQR (sparsified QR Gnanasekaran and Darve in SIAM J Matrix Anal Appl 43(1):94–123, 2022), that was developed by the authors to solve large sparse linear systems. Our algorithm is built on top of a Nested Dissection based multifrontal QR approach. We use low-rank approximations on the frontal matrices to sparsify the vertex separators at every level in the elimination tree. Using a two-step sparsification scheme, we reduce the number of columns and maintain the ratio of rows to columns in each front without introducing any additional fill-in. With this improvised scheme, we show that the runtime of the algorithm scales as \(\mathcal {O}(M \log N)\) and uses \(\mathcal {O}(M)\) memory to store the factorization. This is achieved at the expense of a small and controllable approximation error. The end result is an approximate factorization of the matrix stored as a sequence of sparse orthogonal and upper-triangular factors and hence easy to apply/solve with a vector. Numerical experiments compare the performance of the spaQR algorithm with direct multifrontal QR, an inner-outer iterative method, and CGLS iterative method with a diagonal preconditioner and Robust Incomplete Factorization (RIF) preconditioner.

GPU-accelerated N-1 static security analysis based on fine-grained parallelism HELM

With the increment of scale of power system, the accuracy and efficiency requirements of N-1 static security analysis (SSA) keep increasing. In order to satisfy these requirements, this paper proposed a GPU-accelerated N-1 SSA method to reduce computation burden and ensure the accuracy of alternating current power flow (ACPF) and connectivity test simultaneously. First, in the algorithm, ACPF is performed by holomorphic embedding load flow method (HELM). The construction and solution of sparse linear system both adopts fine-grained parallelism. According to the mathematical lemma of combining linear equations, this paper performs fine-grained parallelism Padé approximation of HELM. Second, considering the high time complexity of the original search-based algorithm and adjacency matrix method which is not suitable for sparse graph, this paper proposes a GPU-accelerated Union Set algorithm. Its main advantage is lower time complexity compared to other graph theory methods. And most importantly, it is very suitable for parallelization. Shared memory can be used in this algorithm well. Therefore, it is not complicated to perform a block-level parallelism in GPU-accelerated Union Set. Finally, the overall framework of GPU-accelerated N-1 SSA method is shown in this paper. It has been tested on large-scale power systems of up to 2869 buses. The computation results are compared with MATPOWER and the execution times are compared with other mainstream method to show the improvement in performance.

A Parallel Strategy for Solving Sparse Linear Systems Over Finite Fields

In this paper we describe a number of parallel techniques that were applied to the problem of finding the null-spaces of thousands of large sparse matrices. This collection of matrices were derived from the discrete logarithm problem attack over the finite field $\F_{3^{6\cdot 509}}$ recently carried out by Adj et al. in~\cite{Adj2017}. Our software library was mainly executed in the supercomputer ABACUS~\cite{abacus}, where in total $21,870$ large sparse linear algebra systems were processed. Solving those linear algebra problems involved a computational effort of over $138$ core-years, requiring a memory space of over 645 gigabytes to store the corresponding vector solutions.

Domain decomposition based preconditioner combined local low-rank approximation with global corrections

To solve general sparse linear systems, this paper presents a domain decomposition based parallel preconditioner. Vertex-based partitioning is utilized to reorder the original coefficient matrix, resulting in a s×s block structure. Here, s is the number of subdomains used in the partition. Variables corresponding to the interface nodes are obtained by solving a linear system with coefficient matrix being the Schur complement S of the reordered matrix. Combining local low-rank correction approximation with a global low-rank correction technique to approximate the inverse of S, the method presented in this paper is different with previous Schur complement based preconditioners. The global low-rank correction terms are obtained by using the information comes from the local low-rank correction terms. In addition, variables corresponding to the interior nodes are computed by solving s small linear systems in parallel. Some numerical tests are presented to show the efficiency and robustness of the proposed preconditioner.

Surrogate modeling for porous flow using deep neural networks

Neural networks are a potential method to solve partial differential equations (PDEs) without considering equation discretization, linearization, and solving large sparse linear systems. Some data-driven or physics-informed methods are studied for surrogate modeling and uncertainty quantification tasks for PDE systems. Compared with data-driven methods, physics-informed methods reduce the reliance on labeled data and have better generalizations and broader applications. However, the methods have low accuracy for solving non-stationary problems with strong nonlinearities. This paper proposes a physics-informed method to solve a porous flow equation with a single source or sink as the boundary condition. The present network framework consists of two network modules. The first module uses two hidden layers with frozen weights, obtained by pre-training a simple double hidden layer network using boundary conditions and initial conditions as loss functions. The second module is used as a residual module to learn measurable label data to improve the accuracy of the solution on the nonlinear boundary. The numerical results verify the validity of the PDE solution, and the results show that the method performs better in terms of both solution and prediction than existing methods. In addition, the method allows a good inversion of the permeability from the measurable data.

Partitioning and Reordering for Spike-Based Distributed-Memory Parallel Gauss--Seidel

Gauss--Seidel (GS) is a widely used iterative method for solving sparse linear systems of equations and also known to be effective as a smoother in algebraic multigrid methods. Parallelization of GS is a challenging task since solving the sparse lower triangular system in GS constitutes a sequential bottleneck at each iteration. We propose a distributed-memory parallel GS (dmpGS) by implementing a parallel sparse triangular solver (stSpike) based on the Spike algorithm. stSpike decouples the global triangular system into smaller systems that can be solved concurrently and requires the solution of a much smaller reduced sparse lower triangular system which constitutes a sequential bottleneck. In order to alleviate this bottleneck and to reduce the communication overhead of dmpGS, we propose a partitioning and reordering model consisting of two phases. The first phase is a novel hypergraph partitioning model whose partitioning objective simultaneously encodes minimizing the reduced system size and the communication volume. The second phase is an in-block row reordering method for decreasing the nonzero count of the reduced system. Extensive experiments on a dataset consisting of 359 sparse linear systems verify the effectiveness of the proposed partitioning and reordering model in terms of reducing the communication and the sequential computational overheads. Parallel experiments on 12 large systems using up to 320 cores demonstrate that the proposed model significantly improves the scalability of dmpGS.

Sparse Spectral-Galerkin Method on An Arbitrary Tetrahedron Using Generalized Koornwinder Polynomials

In this paper, we propose a sparse spectral-Galerkin approximation scheme for solving the second-order partial differential equations on an arbitrary tetrahedron. Generalized Koornwinder polynomials are introduced on the reference tetrahedron as basis functions with their various recurrence relations and differentiation properties being explored. The method leads to well-conditioned and sparse linear systems whose entries can either be calculated directly by the orthogonality of the generalized Koornwinder polynomials for differential equations with constant coefficients or be evaluated efficiently via our recurrence algorithm for problems with variable coefficients. Clenshaw algorithms for the evaluation of any polynomial in an expansion of the generalized Koornwinder basis are also designed to boost the efficiency of the method. Finally, numerical experiments are carried out to illustrate the effectiveness of the proposed Koornwinder spectral method.

Numerical Valuation of European and American Options under Fractional Black-Scholes Model

In this paper, we investigate the numerical valuation of European and American options under the time fractional Black-Scholes model. We first apply a coordinate stretching transformation to the asset price so that the spatial region can focus on the vicinity of singularities, which are usually found in the payoff function. The radial basis function finite difference method is used for the spatial discretization, and the improved L1 method is used to deal with the reduced order of convergence for the nonsmooth initial data. We use the operator splitting method for solving the linear complementary problem of American options. The proposed scheme leads to a sparse linear system which is trivial to solve. Moreover, the stability of the proposed numerical scheme is analyzed using Fourier analysis. Numerical experiments demonstrate the accuracy and efficiency of the proposed method.

Exactly Solving Sparse Rational Linear Systems via Roundoff-Error-Free Cholesky Factorizations

Exactly solving sparse symmetric positive definite (SPD) linear systems is a key problem in mathematics, engineering, and computer science. This paper derives two new sparse roundoff-error-free (REF) Cholesky factorization algorithms which exactly solve sparse SPD linear systems $A {x} = {b}$, where $A \in \mathbb{Q}^{n x n}$ and ${x}, {b} \in {Q}^{n x p}$. The key properties of these factorizations are that (1) they exclusively use integer-arithmetic and (2) in the bit-complexity model, they solve the linear system $A {x} = {b}$ in time proportional to the cost of the integer-arithmetic operations. Namely, the overhead related to data structures and ancillary operations (those not strictly required to perform the factorization) is subsumed by the cost of the integer-arithmetic operations that are essential/intrinsic to the factorization. Notably, to date our algorithms are the only exact algorithm for solving SPD linear systems with this asymptotically efficient complexity bound. Computationally, we show that the novel factorizations are faster than both sparse rational-arithmetic LDL and sparse exact LU factorization. Altogether, the derived sparse REF Cholesky factorizations present a framework to solve any rational SPD linear system exactly and efficiently.

Localized MFS for three‐dimensional acoustic inverse problems on complicated domains

AbstractThis paper proposes a semi‐analytical and local meshless collocation method, the localized method of fundamental solutions (LMFS), to address three‐dimensional (3D) acoustic inverse problems in complex domains. The proposed approach is a recently developed numerical scheme with the potential of being mathematically simple, numerically accurate, and requiring less computational time and storage. In LMFS, an overdetermined sparse linear system is constructed by using the known data at the nodes on the accessible boundary and by making the remaining nodes satisfy the governing equation. In the numerical procedure, the pseudoinverse of a matrix is solved via the truncated singular value decomposition, and thus the regularization techniques are not needed in solving the resulting linear system with a well‐conditioned matrix. Numerical experiments, involving complicated geometry and the high noise level, confirm the effectiveness and performance of the LMFS for solving 3D acoustic inverse problems.

A local domain boundary element method for solving the nonlinear fisher KPP diffusion-reaction equation

The nonlinear Fisher-KPP diffusion-reaction equation has been successfully employed to model competitive processes in biological populations. A robust numerical method for solving diffusion problems is the Boundary Element Method (ΒΕΜ). However, the conventional BEM suffers from large CPU memory requirements since it leads to dense coefficient matrices and for nonlinear problems, its efficiency is further reduced due to the appearance of volume integrals. In this work, a Local Domain BEM (LD-BEM) for solving the Fisher KPP equation is presented. The domain of interest is fragmented into small subdomains and the integral representation of the problem is considered separately for each of the considered subdomains. Eliminating the fluxes at all subdomain interfaces, the proposed LD-BEM leads to sparse linear system coefficient matrices and a reduced number of degrees of freedom. Five numerical examples are solved to assess the accuracy of the proposed method.

A scaled wave finite element method for computing scalar wave radiation and scattering in exterior domains

The calculation of high-frequency wave radiation in exterior domains by finite element methods can lead to large computations. In this paper, it is shown that the solution in the exterior domain can be decomposed as a series expansion of functions with an analytical part made from the product of harmonic waves and polynomials in a scaled variable and a numerical part made of a finite element approximation vector. The solution of the radiation or scattering problem can be found by solving a sparse linear system which is set from the dynamic stiffness matrices of several scaled layers around the radiating body. These dynamic stiffness matrices are classical finite element matrices obtained from any finite element software. Moreover, accurate results can be obtained from a small number of terms in the series expansion. Several examples are given to estimate the efficiency of the proposed method.