Parallel sparsity patterns for factored incomplete inverse matrices

Factored approximate inverse preconditioners with dynamic sparsity patterns

Robotics applications have hard time constraints and heavy computational burdens that can greatly benefit from domain-specific hardware accelerators. For the latency-critical problem of robot motion planning and control, there exists a performance gap of at least an order of magnitude between joint actuator response rates and state-of-the-art software solutions. Hardware acceleration can close this gap, but it is essential to define automated hardware design flows to keep the design process agile as applications and robot platforms evolve. To address this challenge, we introduce robomorphic computing: a methodology to transform robot morphology into a customized hardware accelerator morphology. We (i) present this design methodology, using robot topology and structure to exploit parallelism and matrix sparsity patterns in accelerator hardware; (ii) use the methodology to generate a parameterized accelerator design for the gradient of rigid body dynamics, a key kernel in motion planning; (iii) evaluate FPGA and synthesized ASIC implementations of this accelerator for an industrial manipulator robot; and (iv) describe how the design can be automatically customized for other robot models. Our FPGA accelerator achieves speedups of 8× and 86× over CPU and GPU when executing a single dynamics gradient computation. It maintains speedups of 1.9× to 2.9× over CPU and GPU, including computation and I/O round-trip latency, when deployed as a coprocessor to a host CPU for processing multiple dynamics gradient computations. ASIC synthesis indicates an additional 7.2× speedup for single computation latency. We describe how this principled approach generalizes to more complex robot platforms, such as quadrupeds and humanoids, as well as to other computational kernels in robotics, outlining a path forward for future robomorphic computing accelerators.

Parallel algorithms for computing sparse approximations to the inverse of a sparse matrix either use a prescribed sparsity pattern for the approximate inverse or attempt to generate a good pattern as part of the algorithm. This paper demonstrates that, for PDE problems, the patterns of powers of sparsified matrices (PSMs) can be used a priori as effective approximate inverse patterns, and that the additional effort of adaptive sparsity pattern calculations may not be required. PSM patterns are related to various other approximate inverse sparsity patterns through matrix graph theory and heuristics concerning the PDE's Green's function. A parallel implementation shows that PSM-patterned approximate inverses are significantly faster to construct than approximate inverses constructed adaptively, while often giving preconditioners of comparable quality.

Consistent Bayesian sparsity selection for high-dimensional Gaussian DAG models with multiplicative and beta-mixture priors

We consider the joint sparse estimation of regression coefficients and the covariance matrix for covariates in a high-dimensional regression model, where the predictors are both relevant to a response variable of interest and functionally related to one another via a Gaussian directed acyclic graph (DAG) model. Gaussian DAG models introduce sparsity in the Cholesky factor of the inverse covariance matrix, and the sparsity pattern in turn corresponds to specific conditional independence assumptions on the underlying predictors. A variety of methods have been developed in recent years for Bayesian inference in identifying such network-structured predictors in regression setting, yet crucial sparsity selection properties for these models have not been thoroughly investigated. In this paper, we consider a hierarchical model with spike and slab priors on the regression coefficients and a flexible and general class of DAG-Wishart distributions with multiple shape parameters on the Cholesky factors of the inverse covariance matrix. Under mild regularity assumptions, we establish the joint selection consistency for both the variable and the underlying DAG of the covariates when the dimension of predictors is allowed to grow much larger than the sample size. We demonstrate that our method outperforms existing methods in selecting network-structured predictors in several simulation settings.

A recursive algorithm for decomposition and creation of the inverse of the genomic relationship matrix

We propose a sparse algebra for samplet compressed kernel matrices to enable efficient scattered data analysis. We show that the compression of kernel matrices by means of samplets produces optimally sparse matrices in a certain S-format. The compression can be performed in cost and memory that scale essentially linearly with the number of data points for kernels of finite differentiability. The same holds true for the addition and multiplication of S-formatted matrices. We prove that the inverse of a kernel matrix, given that it exists, is compressible in the S-format as well. The use of selected inversion allows to directly compute the entries in the corresponding sparsity pattern. Moreover, S-formatted matrix operations enable the efficient, approximate computation of more complicated matrix functions such as Aα\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${\\varvec{A}}^\\alpha $$\\end{document} or exp(A)\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\\exp ({\\varvec{A}})$$\\end{document} of a matrix A\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${\\varvec{A}}$$\\end{document}. The matrix algebra is justified mathematically by pseudo differential calculus. As an application, we consider Gaussian process learning algorithms for implicit surfaces. Numerical results are presented to illustrate and quantify our findings.

During the last decades explicit preconditioning methods have gained interest among the scientific community, due to their efficiency for solving large sparse linear systems in conjunction with Krylov subspace iterative methods. The effectiveness of explicit preconditioning schemes relies on the fact that they are close approximants to the inverse of the coefficient matrix. Herewith, we propose a Generic Approximate Sparse Inverse (GenASPI) matrix algorithm based on ILU(0) factorization. The proposed scheme applies to matrices of any structure or sparsity pattern unlike the previous dedicated implementations. The new scheme is based on the Generic Approximate Banded Inverse (GenAbI), which is a banded approximate inverse used in conjunction with Conjugate Gradient type methods for the solution of large sparse linear systems. The proposed GenASPI matrix algorithm, is based on Approximate Inverse Sparsity patterns, derived from powers of sparsified matrices and is computed with a modified procedure based on the GenAbI algorithm. Finally, applicability and implementation issues are discussed and numerical results along with comparative results are presented.

Dubois, Greenbaum and Rodrigue proposed using a truncated Neumann series as an approximation to the inverse of a matrix A for the purpose of preconditioning conjugate gradient iterative approximations to $Ax = b$. If we assume that A has been symmetrically scaled to have unit diagonal and is thus of the form $(I - G)$, then the Neumann series is a power series in G with unit coefficients. The incomplete inverse was thought of as a replacement of the incomplete Cholesky decomposition suggested by Meijerink and van der Vorst in the family of methods ICCG $(n)$. The motivation for the replacement was the desire to have a preconditioned conjugate gradient method which only involved vector operations and which utilized long vectors. We here suggest parameterizing the incomplete inverse to form a preconditioning matrix whose inverse is a polynomial in G. We then show how to select the parameters to minimize the condition number of the product of the polynomial and $(I - G)$. Theoretically the resulting algorithm is the best of the class involving polynomial preconditioners. We also show that polynomial preconditioners which minimize the mean square error with respect to a large class of weight functions are positive definite. We give recurrence relations for the computation of both classes of polynomial preconditioners.

ABSTRACTThe solution of large sparse linear systems is an essential part in many scientific fields. Numerical solution of such systems is usually performed using iterative methods in conjunction with effective preconditioning schemes. A new class of factored approximate inverses is proposed, namely adaptive factored incomplete inverse matrices, which are computed using a recursive Schur complement‐based approach. This class of approximate inverses does not require the knowledge of a sparsity pattern, which is formed adaptively during computation. For this to be possible, a flexible sparse storage scheme was designed. Several numerical dropping strategies are showcased and discussed, along with the effects of reordering to the preconditioner density and the corresponding convergence behavior. Dropping based on monitoring the growth of elements in the incomplete LU factorization is also presented and discussed along with improvements during computation of the successive Schur complements. A static sparsity pattern variant is also provided and discussed. Implementation details and analysis, for computing the proposed scheme, are given along with the computational complexity and memory requirements. Numerical results depicting the effectiveness and applicability of the proposed scheme are also presented.

Covariance estimation and selection for high-dimensional multivariate datasets is a fundamental problem in modern statistics. Gaussian directed acyclic graph (DAG) models are a popular class of models used for this purpose. Gaussian DAG models introduce sparsity in the Cholesky factor of the inverse covariance matrix, and the sparsity pattern in turn corresponds to specific conditional independence assumptions on the underlying variables. A variety of priors have been developed in recent years for Bayesian inference in DAG models, yet crucial convergence and sparsity selection properties for these models have not been thoroughly investigated. Most of these priors are adaptations/generalizations of the Wishart distribution in the DAG context. In this paper, we consider a flexible and general class of these “DAG-Wishart” priors with multiple shape parameters. Under mild regularity assumptions, we establish strong graph selection consistency and establish posterior convergence rates for estimation when the number of variables $p$ is allowed to grow at an appropriate subexponential rate with the sample size $n$.

We analyze the proximal Newton method for minimizing a sum of a self-concordant function and a convex function with an inexpensive proximal operator. We present new results on the global and local convergence of the method when inexact search directions are used. The method is illustrated with an application to L1-regularized covariance selection, in which prior constraints on the sparsity pattern of the inverse covariance matrix are imposed. In the numerical experiments the proximal Newton steps are computed by an accelerated proximal gradient method, and multifrontal algorithms for positive definite matrices with chordal sparsity patterns are used to evaluate gradients and matrix-vector products with the Hessian of the smooth component of the objective.

We review progress in kinetic plasma modeling by electromagnetic particle-in-cell (EM-PIC) algorithms on unstructured grids. These algorithms are implemented in modular CONPIC and BORPIC C++ codes which integrate a matrix-free explicit finite-element (FE) Maxwell solver based on a parallel sparse-approximate inverse (SPAI) algorithm and a first-principles charge-conserving scatter algorithm to transfer charged particle information into dynamic variables on the grid. The Maxwell solver of the EM-PIC algorithm utilizes a mixed FE basis and discretizes the time-dependent coupled first-order Maxwell's system explicitly. The explicit solver approximates the inverse FE system matrix ("mass" matrix) using hierarchical sparsity patterns based on the sparsity pattern of the original matrix. The resulting algorithm effectively accounts for multiscale plasma phenomena. We discuss the application of the developed EM-PIC algorithm to the analysis of laboratory plasmas, vacuum electronic devices for generation of high power microwave signals, and RF electronics multipactor effects in harsh space environments.

We review progress in kinetic plasma modeling by electromagnetic particle-in-cell (EM-PIC) algorithms on unstructured grids. These algorithms are implemented in modular CONPIC and BORPIC C++ codes that integrate a matrix-free explicit finite-element (FE) Maxwell solver based on a parallel sparse-approximate inverse (SPAI) algorithm and a first-principles charge-conserving scatter algorithm to transfer the charged particle information into dynamic variables on the grid. The Maxwell solver of the EM-PIC algorithm utilizes a mixed FE basis and discretizes the time-dependent coupled first-order Maxwell’s system explicitly. The explicit solver approximates the inverse FE system matrix (“mass” matrix) using hierarchical sparsity patterns based on the sparsity pattern of the original matrix. The resulting algorithm effectively accounts for multiscale plasma phenomena. We discuss the application of the developed EM-PIC algorithm to the analysis of laboratory plasmas, vacuum electronic devices for generation of high-power microwave signals, and RF electronics multipactor effects and apply the algorithm to the analysis of multipactor effects and its mitigation in coaxial cables.

Using recursion to compute the inverse of the genomic relationship matrix