A power sparse approximate inverse preconditioning procedure for large sparse linear systems

Motivated by the Cayley–Hamilton theorem, a novel adaptive procedure, called a Power Sparse Approximate Inverse (PSAI) procedure, is proposed that uses a different adaptive sparsity pattern selection approach to constructing a right preconditioner M for the large sparse linear system Ax=b. It determines the sparsity pattern of M dynamically and computes the n independent columns of M that is optimal in the Frobenius norm minimization, subject to the sparsity pattern of M. The PSAI procedure needs a matrix–vector product at each step and updates the solution of a small least squares problem cheaply. To control the sparsity of M and develop a practical PSAI algorithm, two dropping strategies are proposed. The PSAI algorithm can capture an effective approximate sparsity pattern of A^?1 and compute a good sparse approximate inverse M efficiently. Numerical experiments are reported to verify the effectiveness of the PSAI algorithm. Numerical comparisons are made for the PSAI algorithm and the adaptive SPAI algorithm proposed by Grote and Huckle as well as for the PSAI algorithm and three static Sparse Approximate Inverse (SAI) algorithms. The results indicate that the PSAI algorithm is at least comparable to and can be much more effective than the adaptive SPAI algorithm and it often outperforms the static SAI algorithms very considerably and is more robust and practical than the static ones for general problems. Copyright © 2008 John Wiley & Sons, Ltd.     

稀疏近似逆与多层块ILU预条件技术

设计了一种求解一般稀疏线性方程组的健壮且有效的可并行化预条件子,这种预条件子涉及在多层块ILU预条件子(BILUM)中使用稀疏近似逆(AINV)技术．所得的预条件子保持了BILUM的健壮性,它比标准的BILUM预条件子有两点优势:控制稀疏性的能力和增强了并行性．数值例子显示了新预条件子的有效性和效率．     

SPAI Preconditioners for HPC Applications

Iterative methods to solve systems of linear equations Ax = b usually require preconditioners M to speed convergence. But the calculation of many preconditioners can be notoriously sequential. The sparse approximate inverse preconditioner (SPAI) has particular potential for high performance computing [1]. We have ported the SPAI algorithm to graphical processing units (GPUs) within NVIDIA's CUSP library [2] for sparse linear algebra. GPUs perform well on dense linear algebra problems where data resides for long periods on the device. Since the underlying minimization problems are independent, they are mapped to separate thread-blocks, and an optimized QR algorithm implemented using NVIDIA's CUBLAS library is employed on each. Traditionally the challenge has been to determine a sparsity pattern Sp( M ) of the preconditioner dynamically [3], which reduces the condition number of MA to a degree where a preconditioned iterative solver such as GMRES becomes computationally viable. Due to the extremely high performance of the GPU, it is possible to consider initial sparsity patterns much denser than have been previously considered. We therefore consider a fixed sparsity patterns to simplify the GPU implementation. We evaluate the performance of the resulting preconditioner on a standard set of sparse matrices, and compare SPAI to other preconditioners. (© 2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Sparse bivariate polynomial factorization

Motivated by Sasaki’s work on the extended Hensel construction for solving multivariate algebraic equations, we present a generalized Hensel lifting, which takes advantage of sparsity, for factoring bivariate polynomial over the rational number field. Another feature of the factorization algorithm presented in this article is a new recombination method, which can solve the extraneous factor problem before lifting based on numerical linear algebra. Both theoretical analysis and experimental data show that the algorithm is efficient, especially for sparse bivariate polynomials.     

A multilevel Crout ILU preconditioner with pivoting and row permutation

In this paper, we present a new incomplete LU factorization using pivoting by columns and row permutation. Pivoting by columns helps to avoid small pivots and row permutation is used to promote sparsity. This factorization is used in a multilevel framework as a preconditioner for iterative methods for solving sparse linear systems. In most multilevel incomplete ILU factorization preconditioners, preprocessing (scaling and permutation of rows and columns of the coefficient matrix) results in further improvements. Numerical results illustrate that these preconditioners are suitable for a wide variety of applications. Copyright © 2007 John Wiley & Sons, Ltd.     

Sparse approximate solutions to stochastic Galerkin equations

In this Note, we formulate a sparse Krylov-based algorithm for solving large-scale linear systems of algebraic equations arising from the discretization of randomly parametrized (or stochastic) elliptic partial differential equations (SPDEs). We analyze the proposed sparse conjugate gradient (CG) algorithm within the framework of inexact Krylov subspace methods, prove its convergence and study its abstract computational cost. Numerical studies conducted on stochastic diffusion models show that the proposed sparse CG algorithm outperforms the classical CG method when the sought solutions admit a sparse representation in a polynomial chaos basis. In such cases, the sparse CG algorithm recovers almost exactly the sparsity pattern of the exact solutions, which enables accelerated convergence. In the case when the SPDE solution does not admit a sparse representation, the convergence of the proposed algorithm is very similar to the classical CG method.     

Efficient computation of sparse approximate inverses

We investigate different methods for computing a sparse approximate inverse M for a given sparse matrix A by minimizing ∥AM − E∥ in the Frobenius norm. Such methods are very useful for deriving preconditioners in iterative solvers, especially in a parallel environment. We compare different strategies for choosing the sparsity structure of M and different ways for solving the small least squares problem that are related to the computation of each column of M. Especially we show how we can take full advantage of the sparsity of A. © 1998 John Wiley & Sons, Ltd.     

A PROCESS FOR SOLVING A FEW EXTREME EIGENPAIRS OF LARGE SPARSE POSITIVE DEFINITE GENERALIZED EIGENVALUE PROBLEM

1. IntroductionWe are concerned in this work with finding a few extreme eigenvalues and theircorresponding eigenvectors of a generalized large scale eigenvalue problem in which thematrices are sparse and symmetric positive definite.Although finding a few extreme eigenpairs is of interest both in theory and practice,there are only few usable and efficient methods up to now. Reinsch and Baner ([12]),suggested a oR algorithm with Newton shift for the standard eigenproblem which included an ingen…     

Optimization of unconstrained functions with sparse Hessian matrices—Quasi-Newton methods

Newton-type methods and quasi-Newton methods have proven to be very successful in solving dense unconstrained optimization problems. Recently there has been considerable interest in extending these methods to solving large problems when the Hessian matrix has a known a priori sparsity pattern. This paper treats sparse quasi-Newton methods in a uniform fashion and shows the effect of loss of positive-definiteness in generating updates. These sparse quasi-Newton methods coupled with a modified Cholesky factorization to take into account the loss of positive-definiteness when solving the linear systems associated with these methods were tested on a large set of problems. The overall conclusions are that these methods perform poorly in general—the Hessian matrix becomes indefinite even close to the solution and superlinear convergence is not observed in practice. Research for this paper was performed at the Department of Operations Research, Stanford, CA 94305. The research was partially supported by the Department of Energy Contract AM03-76SF00326. PA# DE-AT03-76ER72018: Office of Naval Research Contract N00014-75-C-0267; National Science Foundation Grants MCS-7681259, MCS-7926009 and ECS-8012974.     

Computational experience with methods for estimating sparse hessians for nonlinear optimization

The paper compares a factorized sparse quasi-Newton update of Goldfarb with a nonfactorized BFGS sparse update of Shanno on a series of test problems, with numerical results strongly favoring the unfactorized update. Analysis of Goldfarb's method is done to explain the poor numerical performance. Two specific conjugate gradient methods for solving the required systems of linear equations with the unfactorized update are described and tested.This research was supported by the National Science Foundation under Grant No. MCS-77-07327     

On iterative techniques for computing flow in large two-dimensional discrete fracture networks

Computation of flow in discrete fracture networks often involves solving for hydraulic head values at all intersection points of a large number of stochastically generated fractures inside a bounded domain. For large systems, this approach leads to the generation of problems involving highly sparse matrices which must be solved iteratively. Distributions of fracture lengths spanning over several orders of magnitude, and the randomness of fracture orientations and locations, lead to coefficient matrices that are devoid of any regular structure in the sparsity pattern. In addition to the rapid increase in computational effort with increase in the size of the fracture network, the spread in the distribution of fracture parameters, such as length and transmissivity, dramatically influences the convergence behavior of the system of linear equations. An overview of the discrete fracture network (DFN) methodology for computation of flow is presented along with a comparative study of various Krylov subspace iterative methods for the resulting class of sparse matrices. The rate of convergence of the iterative techniques is found to exhibit a systematic pattern with respect to changes in statistical parameters of the stochastically generated fracture networks. Salient features of the observed trends in the convergence pattern are discussed and guidelines for design of DFN algorithms are provided.     

Solving sparse linear systems by an abs-method that corresponds toLU-decomposition

Abaffy, Broyden, and Spedicato (ABS) have recently proposed a class of direct methods for solving nonsparse linear systems. It is the purpose of this paper to demonstrate that with proper choice of parameters, ABS methods exploit sparsity in a natural way. In particular, we study the choice of parameters which corresponds to an LU-decomposition of the coefficient matrix. In many cases, the fill-in, represented by the nonzero elements of the deflection matrix, uses less storage than the corresponding fill-in of an explicit LU factorization. Hence the above can be a viable and simple method for solving sparse linear systems. A simple reordering scheme for the coefficient matrix is also proposed for the purpose of reducing fill-in of the deflection matrices.     

On the evaluation of general sparse hybrid linear solvers

General sparse hybrid solvers are commonly used kernels for solving wide range of scientific and engineering problems. This work addresses the current problems of efficiently solving general sparse linear equations with direct/iterative hybrid solvers on many core distributed clusters. We briefly discuss the solution stages of Maphys, HIPS, and PDSLin hybrid solvers for large sparse linear systems with their major algorithmic differences. In this category of solvers, different methods with sophisticated preconditioning algorithms are suggested to solve the trade off between memory and convergence. Such solutions require a certain hierarchical level of parallelism more suitable for modern supercomputers that allow to scale for thousand numbers of processors using Schur complement framework. We study the effect of reordering and analyze the performance, scalability as well as memory for each solve phase of PDSLin, Maphys, and HIPS hybrid solvers using large set of challenging matrices arising from different actual applications and compare the results with SuperLU_DIST direct solver. We specifically focus on the level of parallel mechanisms used by the hybrid solvers and the effect on scalability. Tuning and Analysis Utilities (TAU) is employed to assess the efficient usage of heap memory profile and measuring communication volume. The tests are run on high performance large memory clusters using up to 512 processors.     

预处理CG算法解油藏模拟问题的有效性比较

1引言 在大型科学和工程计算问题的实际应用中，经常会遇到求解除椭圆型或抛物线型偏微分方程问题。经差分法或有限元方法离散化后得到一个大型稀疏线性方程组。本文比较了几     

Regularizers for structured sparsity

We study the problem of learning a sparse linear regression vector under additional conditions on the structure of its sparsity pattern. This problem is relevant in machine learning, statistics and signal processing. It is well known that a linear regression can benefit from knowledge that the underlying regression vector is sparse. The combinatorial problem of selecting the nonzero components of this vector can be “relaxed” by regularizing the squared error with a convex penalty function like the ?₁ norm. However, in many applications, additional conditions on the structure of the regression vector and its sparsity pattern are available. Incorporating this information into the learning method may lead to a significant decrease of the estimation error. In this paper, we present a family of convex penalty functions, which encode prior knowledge on the structure of the vector formed by the absolute values of the regression coefficients. This family subsumes the ?₁ norm and is flexible enough to include different models of sparsity patterns, which are of practical and theoretical importance. We establish the basic properties of these penalty functions and discuss some examples where they can be computed explicitly. Moreover, we present a convergent optimization algorithm for solving regularized least squares with these penalty functions. Numerical simulations highlight the benefit of structured sparsity and the advantage offered by our approach over the Lasso method and other related methods.     

A High-Quality Preconditioning Technique for Multi-Length-Scale Symmetric Positive Definite Linear Systems

We study preconditioning techniques used in conjunction with the conjugate gradient method for solving multi-length-scale symmetric positive definite linear systems originating from the quantum Monte Carlo simulation of electron interaction of correlated materials. Existing preconditioning techniques are not designed to be adaptive to varying numerical properties of the multi-length-scale systems. In this paper, we propose a hybrid incomplete Cholesky (HIC) preconditioner and demonstrate its adaptivity to the multi-length-scale systems. In addition, we propose an extension of the compressed sparse column with row access (CSCR) sparse matrix storage format to efficiently accommodate the data access pattern to compute the HIC preconditioner. We show that for moderately correlated materials, the HIC preconditioner achieves the optimal linear scaling of the simulation. The development of a linear-scaling preconditioner for strongly correlated materials remains an open topic.     

The generalized residual cutting method and its convergence characteristics

Iterative methods and especially Krylov subspace methods (KSM) are a very useful numerical tool in solving for large and sparse linear systems problems arising in science and engineering modeling. More recently, the nested loop KSM have been proposed that improve the convergence of the traditional KSM. In this article, we review the residual cutting (RC) and the generalized residual cutting (GRC) that are nested loop methods for large and sparse linear systems problems. We also show that GRC is a KSM that is equivalent to Orthomin with a variable preconditioning. We use the modified Gram–Schmidt method to derive a stable GRC algorithm. We show that GRC presents a general framework for constructing a class of “hybrid” (nested) KSM based on inner loop method selection. We conduct numerical experiments using nonsymmetric indefinite matrices from a widely used library of sparse matrices that validate the efficiency and the robustness of the proposed methods.     

An efficient analyse phase for element problems

The analyse phase of a sparse direct solver for symmetrically structured linear systems of equations is used to determine the sparsity pattern of the matrix factor. This allows the subsequent numerical factorisation and solve phases to be executed efficiently. Many direct solvers require the system matrix to be in assembled form. For problems arising from finite element applications, assembling and then using the system matrix can be costly in terms of both time and memory. This paper describes and implements a variant of the work of Gilbert, Ng and Peyton for matrices in elemental form. The proposed variant works with an equivalent matrix that avoids explicitly assembling the system matrix and exploits supervariables. Numerical experiments using problems from practical applications are used to demonstrate the significant advantages of working directly with the elemental form. Copyright © 2012 John Wiley & Sons, Ltd.     

Exploiting sparsity in linear and nonlinear matrix inequalities via positive semidefinite matrix completion

A basic framework for exploiting sparsity via positive semidefinite matrix completion is presented for an optimization problem with linear and nonlinear matrix inequalities. The sparsity, characterized with a chordal graph structure, can be detected in the variable matrix or in a linear or nonlinear matrix-inequality constraint of the problem. We classify the sparsity in two types, the domain-space sparsity (d-space sparsity) for the symmetric matrix variable in the objective and/or constraint functions of the problem, which is required to be positive semidefinite, and the range-space sparsity (r-space sparsity) for a linear or nonlinear matrix-inequality constraint of the problem. Four conversion methods are proposed in this framework: two for exploiting the d-space sparsity and the other two for exploiting the r-space sparsity. When applied to a polynomial semidefinite program (SDP), these conversion methods enhance the structured sparsity of the problem called the correlative sparsity. As a result, the resulting polynomial SDP can be solved more effectively by applying the sparse SDP relaxation. Preliminary numerical results on the conversion methods indicate their potential for improving the efficiency of solving various problems.