Sparse Systems Of Linear Equations Research Articles

Smoothness on rank-order path graphs and its use in compressive spectral imaging with side information

This paper proposes a novel reconstruction approach to compressive spectral imaging (CSI) with panchromatic side information, which is based on the notion of approximate rank-order statistics. To that end, we assume that the signal of interest is sufficiently smooth on an unknown graph. When restricted to the family of path graphs, we show that the best path is indeed the rank-order path graph induced by the signal. That is, the path graph whose edge structure is given by the permutation that sorts the entries of the signal in ascending order. Our goal is to show that smoothness on rank-order path graphs inferred from the rank-order statistics of a co-registered panchromatic signal can be used to find accurate spectral image estimates from a compressive snapshot of the scene. We derive theoretical properties of rank-order path graphs and give illustrative examples of their use in signal recovery from undersampled measurements. Our approach leads to solutions with a closed-form, found efficiently by iterative inversion of highly sparse systems of linear equations. We evaluate our method through an experimental demonstration and extensive simulations. Our method performs notably better against a bilateral-filter graph model, adapted to the task, and some traditional and state-of-the-art algorithms.

ParGeMSLR: A parallel multilevel Schur complement low-rank preconditioning and solution package for general sparse matrices

This paper discusses parGeMSLR, a C++/MPI software library for the solution of sparse systems of linear algebraic equations via preconditioned Krylov subspace methods in distributed-memory computing environments. The preconditioner implemented in parGeMSLR is based on algebraic domain decomposition and partitions the symmetrized adjacency graph recursively into several non-overlapping partitions via a p-way vertex separator, where p is an integer multiple of the total number of MPI processes. From a numerical perspective, parGeMSLR builds a Schur complement approximate inverse preconditioner as the sum between the matrix inverse of the interface coupling matrix and a low-rank correction term. To reduce the cost associated with the computation of the approximate inverse matrices, parGeMSLR exploits a multilevel partitioning of the algebraic domain. The parGeMSLR library is implemented on top of the Message Passing Interface and can solve both real and complex linear systems. Furthermore, parGeMSLR can take advantage of hybrid computing environments with in-node access to one or more Graphics Processing Units. Finally, the parallel efficiency (weak and strong scaling) of parGeMSLR is demonstrated on a few model problems arising from discretizations of 3D Partial Differential Equations.

AAQAL: A Machine Learning-Based Tool for Performance Optimization of Parallel SPMV Computations Using Block CSR

The sparse matrix–vector product (SpMV), considered one of the seven dwarfs (numerical methods of significance), is essential in high-performance real-world scientific and analytical applications requiring solution of large sparse linear equation systems, where SpMV is a key computing operation. As the sparsity patterns of sparse matrices are unknown before runtime, we used machine learning-based performance optimization of the SpMV kernel by exploiting the structure of the sparse matrices using the Block Compressed Sparse Row (BCSR) storage format. As the structure of sparse matrices varies across application domains, optimizing the block size is important for reducing the overall execution time. Manual allocation of block sizes is error prone and time consuming. Thus, we propose AAQAL, a data-driven, machine learning-based tool that automates the process of data distribution and selection of near-optimal block sizes based on the structure of the matrix. We trained and tested the tool using different machine learning methods—decision tree, random forest, gradient boosting, ridge regressor, and AdaBoost—and nearly 700 real-world matrices from 43 application domains, including computer vision, robotics, and computational fluid dynamics. AAQAL achieved 93.47% of the maximum attainable performance with a substantial difference compared to in practice manual or random selection of block sizes. This is the first attempt at exploiting matrix structure using BCSR, to select optimal block sizes for the SpMV computations using machine learning techniques.

Strategies of Preconditioner Updates for Sequences of Linear Systems Associated with the Neutron Diffusion

The time-dependent neutron diffusion equation approximates the neutronic power evolution inside a nuclear reactor core. Applying a Galerkin finite element method for the spatial discretization of these equations leads to a stiff semi-discrete system of ordinary differential equations. For time discretization, an implicit scheme is used, which implies solving a large and sparse linear system of equations for each time step. The GMRES method is used to solve these systems because of its fast convergence when a suitable preconditioner is provided. This work explores several matrix-free strategies based on different updated preconditioners, which are constructed by low-rank updates of a given initial preconditioner. They are two tuned preconditioners based on the bad and good Broyden’s methods, initially developed for nonlinear equations and optimization problems, and spectral preconditioners. The efficiency of the resulting preconditioners under study is closely related to the selection of the subspace used to construct the update. Our numerical results show the effectiveness of these methodologies in terms of CPU time and storage for different nuclear benchmark transients, even if the initial preconditioner is not good enough.

Incorporating Gradient Similarity for Robust Time Delay Estimation in Ultrasound Elastography.

Energy-based ultrasound elastography techniques minimize a regularized cost function consisting of data and continuity terms to obtain local displacement estimates based on the local time-delay estimation (TDE) between radio frequency (RF) frames. The data term associated with the existing techniques takes only the amplitude similarity into account and hence is not sufficiently robust to the outlier samples present in the RF frames under consideration. This drawback creates noticeable artifacts in the strain image. To resolve this issue, we propose to formulate the data function as a linear combination of the amplitude and gradient similarity constraints. We estimate the adaptive weight concerning each similarity term following an iterative scheme. Finally, we optimize the nonlinear cost function in an efficient manner to convert the problem to a sparse system of linear equations which are solved for millions of variables. We call our technique rGLUE: robust data term in GLobal Ultrasound Elastography. rGLUE has been validated using simulation, phantom, in vivo liver, and breast datasets. In all our experiments, rGLUE substantially outperforms the recent elastography methods both visually and quantitatively. For simulated, phantom, and in vivo datasets, respectively, rGLUE achieves 107%, 18%, and 23% improvements of signal-to-noise ratio (SNR) and 61%, 19%, and 25% improvements of contrast-to-noise ratio (CNR) over global ultrasound elastography (GLUE), a recently published elastography algorithm.

Solving the Eigenfrequencies Problem of Waveguides by Localized Method of Fundamental Solutions with External Source

The localized method of fundamental solutions (LMFS) is a domain-type, meshless numerical method. Compared with numerical methods that have a high grid dependence, it does not require grid generation and numerical integration, so it can effectively improve computational efficiency and avoid complex integration processes. Moreover, it is formed using the traditional method of fundamental solutions (MFS) and the localization approach. Previous studies have shown that the MFS may produce a dense and ill-conditioned matrix. However, the proposed LMFS can yield a sparse system of linear algebraic equations, so it is more suitable and effective in solving complicated engineering problems. In this article, LMFS was used to solve eigenfrequency problems in electromagnetic waves, which were controlled using two-dimensional Helmholtz equations. Additionally, the resonant frequencies of the eigenproblem were determined by the response amplitudes. In order to determine the eigenfrequencies, LMFS was applied for solving a sequence of inhomogeneous problems by introducing an external source. Waveguides with different shapes were analyzed to prove the stability of the present LMFS in this paper.

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.

Block-based spectral image reconstruction for compressive spectral imaging using smoothness on graphs.

A novel reconstruction method for compressive spectral imaging is designed by assuming that the spectral image of interest is sufficiently smooth on a collection of graphs. Since the graphs are not known in advance, we propose to infer them from a panchromatic image using a state-of-the-art graph learning method. Our approach leads to solutions with closed-form that can be found efficiently by solving multiple sparse systems of linear equations in parallel. Extensive simulations and an experimental demonstration show the merits of our method in comparison with traditional methods based on sparsity and total variation and more recent methods based on low-rank minimization and deep-based plug-and-play priors. Our approach may be instrumental in designing efficient methods based on deep neural networks and covariance estimation.

Pipelined Preconditioned Conjugate Gradient Methods for real and complex linear systems for distributed memory architectures

Preconditioned Conjugate Gradient (PCG) is a popular method for solving large and sparse linear systems of equations. The performance of PCG at scale is affected due to the costly global synchronization steps that arise in dot-products on distributed memory systems. Pipelined PCG (PIPECG) removes the costly global synchronization steps from PCG by only executing a single non-blocking allreduce per iteration and overlapping it with independent computations.In our previous work, we have developed a novel pipelined PCG algorithm called PIPECG-OATI (One Allreduce per Two Iterations) for real linear systems which executes a single non-blocking allreduce per two iterations and provides a large overlap of global communication with independent computations at higher number of cores. Our method achieves this overlap by using iteration combination and by introducing new recurrence and non-recurrence computations. We implement optimizations in the PIPECG-OATI method to use cache memory efficiently.In this work, we present PIPECG-OATI-c method for linear systems with complex Hermitian positive definite and complex symmetric matrices. We compare our method with various pipelined CG methods on a variety of problems and demonstrate that our method always gives the least run times. We performed experiments with our method using 20M and 30M unknowns on up to 16K cores and obtained up to 2.48X performance improvement over PCG and 2.14X performance improvement over PIPECG methods. We also experimented with up to 1-billion unknowns on 16K cores, the largest problem size explored for the CG problem, to our knowledge, and obtained about 25% improvement over PCG.

Efficient Simulations of Large Scale Convective Heat Transfer Problems

We describe an approach for efficient solution of large scale convective heat transfer problems, formulated as coupled unsteady heat conduction and incompressible fluid flow equations. The original problem is discretized in time using classical implicit methods, while stabilized finite elements are used for space discretization. The algorithm employed for the discretization of the fluid flow problem uses Picard's iterations to solve the arising nonlinear equations. Both problems, heat transfer and Navier-Stokes quations, give rise to large sparse systems of linear equations. The systems are solved using iterative GMRES solver with suitable preconditioning. For the incompressible flow equations we employ a special preconditioner based on algebraic multigrid (AMG) technique. The paper presents algorithmic and implementation details of the solution procedure, which is suitably tuned, especially for ill conditioned systems arising from discretizations of incompressible Navier-Stokes equations. We describe parallel implementation of the solver using MPI and elements of PETSC library. The scalability of the solver is favourably compared with other methods such as direct solvers and standard GMRES method with ILU preconditioning.

A distributed Gauss–Newton method for distribution system state estimation

State estimation techniques are expected to play a significant role in the management of future distribution systems. This paper proposes an efficient distributed algorithm for state estimation of distribution networks. The algorithm is designed to solve the sparse linear system of equations at each Gauss–Newton iteration with a multi-area architecture. According to this algorithm, the augmented matrix of the linear system is fully decomposed. Each area calculates and stores a part of the augmented matrix using only interior and border measurements. A distributed Gaussian elimination procedure is then performed, which relies solely on simple local calculation and limited data exchange between neighboring areas. The proposed scheme enables the distributed computation of the overall weighted least square estimation problem. Parallel processing abilities of area processors are fully utilized which makes the algorithm computationally efficient. Case studies performed on the unbalanced IEEE 33 and 123-bus systems demonstrate the good performance of the algorithm in terms of accuracy and efficiency.

Multiscale Cholesky preconditioning for ill-conditioned problems

Many computer graphics applications boil down to solving sparse systems of linear equations. While the current arsenal of numerical solvers available in various specialized libraries and for different computer architectures often allow efficient and scalable solutions to image processing, modeling and simulation applications, an increasing number of graphics problems face large-scale and ill-conditioned sparse linear systems --- a numerical challenge which typically chokes both direct factorizations (due to high memory requirements) and iterative solvers (because of slow convergence). We propose a novel approach to the efficient preconditioning of such problems which often emerge from the discretization over unstructured meshes of partial differential equations with heterogeneous and anisotropic coefficients. Our numerical approach consists in simply performing a fine-to-coarse ordering and a multiscale sparsity pattern of the degrees of freedom, using which we apply an incomplete Cholesky factorization. By further leveraging supernodes for cache coherence, graph coloring to improve parallelism and partial diagonal shifting to remedy negative pivots, we obtain a preconditioner which, combined with a conjugate gradient solver, far exceeds the performance of existing carefully-engineered libraries for graphics problems involving bad mesh elements and/or high contrast of coefficients. We also back the core concepts behind our simple solver with theoretical foundations linking the recent method of operator-adapted wavelets used in numerical homogenization to the traditional Cholesky factorization of a matrix, providing us with a clear bridge between incomplete Cholesky factorization and multiscale analysis that we leverage numerically.

Towards distortion-free imaging of the eye.

The high power of the eye and optical components used to image it result in "static" distortion, remaining constant across acquired retinal images. In addition, raster-based systems sample points or lines of the image over time, suffering from "dynamic" distortion due to the constant motion of the eye. We recently described an algorithm which corrects for the latter problem but is entirely blind to the former. Here, we describe a new procedure termed "DIOS" (Dewarp Image by Oblique Shift) to remove static distortion of arbitrary type. Much like the dynamic correction method, it relies on locating the same tissue in multiple frames acquired as the eye moves through different gaze positions. Here, the resultant maps of pixel displacement are used to form a sparse system of simultaneous linear equations whose solution gives the common warp seen by all frames. We show that the method successfully handles torsional movement of the eye. We also show that the output of the previously described dynamic correction procedure may be used as input for this new procedure, recovering an image of the tissue that is, in principle, a faithful replica free of any type of distortion. The method could be extended beyond ocular imaging, to any kind of imaging system in which the image can move or be made to move across the detector.

Accelerated Randomized Methods for Receiver Design in Extra-Large Scale MIMO Arrays

Massive multiple-input-multiple-output (M-MIMO) features a capability for spatial multiplexing of large number of users. This number becomes even more extreme in extra-large (XL-MIMO), a variant of M-MIMO where the antenna array is of very large size. Yet, the problem of signal processing complexity in M-MIMO is further exacerbated by the XL size of the array. The basic processing problem boils down to a sparse system of linear equations that can be addressed by the randomized Kaczmarz (RK) algorithm. This algorithm has recently been applied to devise low-complexity M-MIMO receivers; however, it is limited by the fact that certain configurations of the linear equations may significantly deteriorate the performance of the RK algorithm. In this paper, we embrace the interest in accelerated RK algorithms and introduce three new RK-based low-complexity receiver designs. In our experiments, our methods are not only able to overcome the previous scheme, but they are more robust against inter-user interference (IUI) and sparse channel matrices arising in the XL-MIMO regime. In addition, we show that the RK-based schemes use a mechanism similar to that used by successive interference cancellation (SIC) receivers to approximate the regularized zero-forcing (RZF) scheme.

Analytic Global Regularized Backscatter Quantitative Ultrasound.

Although a variety of techniques have been developed to reduce the appearance of B-mode speckle, quantitative ultrasound (QUS) aims at extracting the hidden properties of the tissue. Herein, we propose two novel techniques to accurately and precisely estimate two important QUS parameters, namely, the average attenuation coefficient and the backscatter coefficient. Both the techniques optimize a cost function that incorporates data and continuity constraint terms, which we call AnaLytical Global rEgularized BackscatteR quAntitative ultrasound (ALGEBRA). We propose two versions of ALGEBRA, namely, 1-D- and 2-D-ALGEBRA. In 1-D-ALGEBRA, the regularized cost function is formulated in the axial direction, and the QUS parameters are calculated for one line of radio frequency (RF) echo data. In 2-D-ALGEBRA, the regularized cost function is formulated for the entire image, and the QUS parameters throughout the image are estimated simultaneously. This simultaneous optimization allows 2-D-ALGEBRA to "see" all the data before estimating the QUS parameters. In both the methods, we efficiently optimize the cost functions by casting it as a sparse linear system of equations. As a result of this efficient optimization, 1-D-ALGEBRA and 2-D-ALGEBRA are, respectively, 600 and 300 times faster than optimization using the dynamic programming (DP) method previously proposed by our group. In addition, the proposed technique has fewer input parameters that require manual tuning. Our results demonstrate that the proposed ALGEBRA methods substantially outperform least-square and DP methods in estimating the QUS parameters in phantom experiments.

Accelerating advanced preconditioning methods on hybrid architectures

Many problems, in diverse areas of science and engineering, involve the solution of largescale sparse systems of linear equations. In most of these scenarios, they are also a computational bottleneck, and therefore their efficient solution on parallel architectureshas motivated a tremendous volume of research.This dissertation targets the use of GPUs to enhance the performance of the solution of sparse linear systems using iterative methods complemented with state-of-the-art preconditioned techniques. In particular, we study ILUPACK, a package for the solution of sparse linear systems via Krylov subspace methods that relies on a modern inverse-based multilevel ILU (incomplete LU) preconditioning technique.We present new data-parallel versions of the preconditioner and the most important solvers contained in the package that significantly improve its performance without affecting its accuracy. Additionally we enhance existing task-parallel versions of ILUPACK for shared- and distributed-memory systems with the inclusion of GPU acceleration. The results obtained show a sensible reduction in the runtime of the methods, as well as the possibility of addressing large-scale problems efficiently.

Low-latitude reduction-to-the-pole and upward continuation between arbitrary surfaces based on the partial differential equation framework

SUMMARY Magnetic surveys conducted in complex conditions, such as low magnetic latitudes, uneven observation surfaces, or above high-susceptibility sources, pose significant challenges for obtaining stable solutions for reduction-to-the-pole (RTP) and upward-continuation processing on arbitrary surfaces. To tackle these challenges, in this study, we propose constructing an equivalent-susceptibility model based on the partial differential equation (PDE) framework in the space domain. A multilayer equivalent-susceptibility method was used for RTP and upward-continuation operations, thus allowing for application on undulating observation surfaces and strong self-demagnetization effect in a non-uniform mesh. A novel positivity constraint is introduced to improve the accuracy and efficiency of the inversion. We analysed the effect of the depth-weighting function in the inversion of equivalent susceptibility for RTP and upward-continuation reproduction. Iterative and direct solvers were utilized and compared in solving the large, sparse, non-symmetric and ill-conditioned system of linear equations produced by PDE-based equivalent-source construction. Two synthetic models were used to illustrate the efficiency and accuracy of the proposed method in processing both ground and airborne magnetic data. Aeromagnetic and ground data collected in Brazil at a low magnetic latitude region were used to validate the proposed method for processing RTP and upward-continuation operations on magnetic data sets with strong self-demagnetization.

Local non-singular knot method for large-scale computation of acoustic problems in complicated geometries

This paper presents a local non-singular knot method (LNKM) to accurately solve the large-scale acoustic problems in complicated geometries. The LNKM is a domain-type meshless collocation method, which relies only on scattered nodes. Firstly, a series of subdomains corresponding to every nodes can be searched based on the Euclidean distance between nodes. To each subdomain, a small linear system can be yielded by using the non-singular general solutions of Helmholtz-type equations. Secondly, the unknown variables at every nodes can be explicitly expressed by the function values at their corresponding supporting nodes. Finally, a large sparse system of linear equations is formed and solved to obtain the numerical solutions of physical quantities at every nodes. The proposed LNKM is mathematically simple, numerically accurate, and more applicable to the large-scale computation. Four numerical examples conform its effectiveness and accuracy for the large-scale computation of Helmholtz-type equations in complicated geometries.

On the preconditioning of three-by-three block saddle point problems

We establish a new iterative method for solving a class of large and sparse linear systems of equations with three-by-three block coefficient matrices having saddle point structure. Convergence properties of the proposed method are studied in details and its induced preconditioner is examined for accelerating the convergence speed of generalized minimal residual (GMRES) method. More precisely, we analyze the eigenvalue distribution of the preconditioned matrix. Numerical experiments are reported to demonstrate the effectiveness of the proposed preconditioner.

Two-Level Nyström--Schur Preconditioner for Sparse Symmetric Positive Definite Matrices

Randomized methods are becoming increasingly popular in numerical linear algebra. However, few attempts have been made to use them in developing preconditioners. Our interest lies in solving large-scale sparse symmetric positive definite linear systems of equations where the system matrix is preordered to doubly bordered block diagonal form (for example, using a nested dissection ordering). We investigate the use of randomized methods to construct high quality preconditioners. In particular, we propose a new and efficient approach that employs Nystr\"om's method for computing low rank approximations to develop robust algebraic two-level preconditioners. Construction of the new preconditioners involves iteratively solving a smaller but denser symmetric positive definite Schur complement system with multiple right-hand sides. Numerical experiments on problems coming from a range of application areas demonstrate that this inner system can be solved cheaply using block conjugate gradients and that using a large convergence tolerance to limit the cost does not adversely affect the quality of the resulting Nystr\"om--Schur two-level preconditioner.