Sparse Matrix Computations Research Articles

Bandwidth of WK-recursive networks and its sparse matrix computation

Bandwidth of WK-recursive networks and its sparse matrix computation

SpChar: Characterizing the sparse puzzle via decision trees

Sparse matrix computation is crucial in various modern applications, including large-scale graph analytics, deep learning, and recommender systems. The performance of sparse kernels varies greatly depending on the structure of the input matrix, making it difficult to gain a comprehensive understanding of sparse computation and its relationship to inputs, algorithms, and target machine architecture. Despite extensive research on certain sparse kernels, such as Sparse Matrix-Vector Multiplication (SpMV), the overall family of sparse algorithms has yet to be investigated as a whole. This paper introduces SpChar, a workload characterization methodology for general sparse computation. SpChar employs tree-based models to identify the most relevant hardware and input characteristics, starting from hardware and input-related metrics gathered from Performance Monitoring Counters (PMCs) and matrices. Our analysis enables the creation of a characterization loop that facilitates the optimization of sparse computation by mapping the impact of architectural features to inputs and algorithmic choices. We apply SpChar to more than 600 matrices from the SuiteSparse Matrix collection and three state-of-the-art Arm Central Processing Units (CPUs) to determine the critical hardware and software characteristics that affect sparse computation. In our analysis, we determine that the biggest limiting factors for high-performance sparse computation are (1) the latency of the memory system, (2) the pipeline flush overhead resulting from branch misprediction, and (3) the poor reuse of cached elements. Additionally, we propose software and hardware optimizations that designers can implement to create a platform suitable for sparse computation. We then investigate these optimizations using the gem5 simulator to achieve a significant speedup of up to 2.63× compared to a CPU where the optimizations are not applied.

Optimization of Sparse Matrix Computation for Algebraic Multigrid on GPUs

AMG is one of the most efficient and widely used methods for solving sparse linear systems. The computational process of AMG mainly consists of a series of iterative calculations of generalized sparse matrix-matrix multiplication (SpGEMM) and sparse matrix-vector multiplication (SpMV). Optimizing these sparse matrix calculations is crucial for accelerating solving linear systems. In this paper, we first focus on optimizing the SpGEMM algorithm in AmgX, a popular AMG library for GPUs. We propose a new algorithm called SpGEMM-upper, which achieves an average speedup of 2.02x on Tesla V100 and 1.96x on RTX 3090 against the original algorithm. Next, through experimental investigation, we conclude that no single SpGEMM library or algorithm performs optimally for most sparse matrices, and the same holds true for SpMV. Therefore, we build machine learning-based models to predict the optimal SpGEMM and SpMV used in the AMG calculation process. Finally, we integrate the prediction models, SpGEMM-upper, and other selected algorithms into a framework for adaptive sparse matrix computation in AMG. Our experimental results prove that the framework achieves promising performance improvements on the test set.

Numerical Solution of Two-Dimensional Nonlinear Unsteady Advection-Diffusion-Reaction Equations with Variable Coefficients

The advection-diffusion-reaction (ADR) equation is a fundamental mathematical model used to describe various processes in many different areas of science and engineering. Due to wide applicability of the ADR equation, finding accurate solution is very important to better understand a physical phenomenon represented by the equation. In this study, a numerical scheme for solving two-dimensional unsteady ADR equations with spatially varying velocity and diffusion coefficients is presented. The equations include nonlinear reaction terms. To discretize the ADR equations, the Crank–Nicolson finite difference method is employed with a uniform grid. The resulting nonlinear system of equations is solved using Newton’s method. At each iteration of Newton’s method, the Gauss–Seidel iterative method with sparse matrix computation is utilized to solve the block tridiagonal system and obtain the error correction vector. The consistency and stability of the numerical scheme are investigated. MATLAB codes are developed to implement this combined numerical approach. The validation of the scheme is verified by solving a two-dimensional advection-diffusion equation without reaction term. Numerical tests are provided to show the good performances of the proposed numerical scheme in simulation of ADR problems. The numerical scheme gives accurate results. The obtained numerical solutions are presented graphically. The result of this study may provide insights to apply numerical methods in solving comprehensive models of physical phenomena that capture the underlying situations.

Advancing Finite Difference Solutions for Two‐Dimensional Incompressible Navier–Stokes Equations Using Artificial Compressibility Method and Sparse Matrix Computation

In this article, a numerical scheme for solving two‐dimensional (2D) time‐dependent incompressible Navier–Stokes equations is presented. The artificial compressibility technique is used to incorporate a time derivative of pressure term to the continuity equation. It is employed for pressure–velocity coupling. The scheme consists of backward difference approximation for time derivatives and central difference approximation for spatial derivatives, implemented on a collocated grid. The discretization of the differential equations yields a system of algebraic equations with a block coefficient matrix. To solve this system efficiently, matrix inversion with sparse matrix computation is employed. The proposed numerical scheme is applied to solve three flow problems (lid‐driven cavity flow, rectangular channel flow, and Taylor–Green vortex problem) to validate the accuracy and applicability of the scheme. The results affirm the scheme’s capability to provide precise approximations for solutions to the Navier–Stokes equations. With slight modifications, the scheme can be applied to solve various flow problems with high accuracy, less memory usage, and reduced computational time.

Modeling the Interplay between Loop Tiling and Fusion in Optimizing Compilers Using Affine Relations

Loop tiling and fusion are two essential transformations in optimizing compilers to enhance the data locality of programs. Existing heuristics either perform loop tiling and fusion in a particular order, missing some of their profitable compositions, or execute ad-hoc implementations for domain-specific applications, calling for a generalized and systematic solution in optimizing compilers. In this article, we present a so-called basteln (an abbreviation for backward slicing of tiled loop nests) strategy in polyhedral compilation to better model the interplay between loop tiling and fusion. The basteln strategy first groups loop nests by preserving their parallelism/tilability and next performs rectangular/parallelogram tiling to the output groups that produce data consumed outside the considered program fragment. The memory footprints required by each tile are then computed, from which the upward exposed data are extracted to determine the tile shapes of the remaining fusion groups. Such a tiling mechanism can construct complex tile shapes imposed by the dependences between these groups, which are further merged by a post-tiling fusion algorithm for enhancing data locality without losing the parallelism/tilability of the output groups. The basteln strategy also takes into account the amount of redundant computations and the fusion of independent groups, exhibiting a general applicability. We integrate the basteln strategy into two optimizing compilers, with one a general-purpose optimizer and the other a domain-specific compiler for deploying deep learning models. The experiments are conducted on CPU, GPU, and a deep learning accelerator to demonstrate the effectiveness of the approach for a wide class of application domains, including deep learning, image processing, sparse matrix computation, and linear algebra. In particular, the basteln strategy achieves a mean speedup of 1.8× over cuBLAS/cuDNN and 1.1× over TVM on GPU when used to optimize deep learning models; it also outperforms PPCG and TVM by 11% and 20%, respectively, when generating code for the deep learning accelerator.

A Tiny Accelerator for Mixed-Bit Sparse CNN Based on Efficient Fetch Method of SIMO SPad

Convolution neural networks (CNNs) have been implemented with custom hardware on edge devices since its algorithms were successful in many artificial intelligence applications. Although lots of unstructured pruning and mix-bit quantization algorithms have been proposed to successfully compress CNNs, there are few hardware accelerators which can support both sparse and mix-bit CNNs. Besides, sparse matrix computation consumes lots of hardware resources such as registers or BRAM to fetch the needed input activations into processing element (PE). This paper presents a tiny accelerator for mixed-bit sparse CNNs featuring a novel scheme of single vector-based compressed sparse filter (CSF) method and single input multiple output scratch pad (SIMO SPad) to effectively compress weight and fetch the needed input activation. SIMO SPad is shared by multiple PEs, which saves 13.34% CLB LUTs, 46.24% CLB Registers. Furthermore, the accelerator supports mixed-bit sparse computation to obtain better accuracy and performance. When tested on VGG16, compared with 8-bit non-sparsity baselines, the performance of mixed-bit sparsity on Cifar10 and ImageNet improved by 4.85× and 3.33×, respectively, with small accuracy decrease degradation. Compared to state-of-the-art accelerators, the accelerator achieves 1.40× to 2.98× greater DSP efficiency, and offers 1.91× greater energy efficiency.

Mapping tree‐shaped workflows on systems with different memory sizes and processor speeds

AbstractDirected acyclic graphs are commonly used to model scientific workflows, by expressing dependencies between tasks, as well as the resource requirements of the workflow. As a special case, rooted directed trees occur in several applications, for instance in sparse matrix computations. Since typical workflows are modeled by large trees, it is crucial to schedule them efficiently, so that their execution time (or makespan) is minimized. Furthermore, it is usually beneficial to distribute the execution on several compute nodes, hence increasing the available memory, and allowing us to parallelize parts of the execution. To exploit the heterogeneity of modern clusters in this context, we investigate the partitioning and mapping of tree‐shaped workflows on two types of target architecture models: in AM1, each processor can have a different memory size, and in AM2, each processor can also have a different speed (in addition to a different memory size). We design a three‐step heuristic for AM1, which adapts and extends previous work for homogeneous clusters [Gou C, Benoit A, Marchal L. Partitioning tree‐shaped task graphs for distributed platforms with limited memory. IEEE Trans Parallel Dist Syst 2020; 31(7): 1533–1544]. The changes we propose concern the assignment to processors (accounting for the different memory sizes) and the availability of suitable processors when splitting or merging subtrees. For AM2, we extend the heuristic for AM1 with a two‐phase local search approach. Phase A is a swap‐based hill climber, while (the optional) Phase B is inspired by iterated local search. We evaluate our heuristics for AM1 and AM2 with extensive simulations, and we demonstrate that exploiting the heterogeneity in the cluster significantly reduces the makespan, compared to the state of the art for homogeneous processors.

Editorial for the special issue on architecture, algorithms and applications of high performance sparse matrix computations

Editorial for the special issue on architecture, algorithms and applications of high performance sparse matrix computations

Longitudinal Modeling of Age-Dependent Latent Traits with Generalized Additive Latent and Mixed Models

We present generalized additive latent and mixed models (GALAMMs) for analysis of clustered data with responses and latent variables depending smoothly on observed variables. A scalable maximum likelihood estimation algorithm is proposed, utilizing the Laplace approximation, sparse matrix computation, and automatic differentiation. Mixed response types, heteroscedasticity, and crossed random effects are naturally incorporated into the framework. The models developed were motivated by applications in cognitive neuroscience, and two case studies are presented. First, we show how GALAMMs can jointly model the complex lifespan trajectories of episodic memory, working memory, and speed/executive function, measured by the California Verbal Learning Test (CVLT), digit span tests, and Stroop tests, respectively. Next, we study the effect of socioeconomic status on brain structure, using data on education and income together with hippocampal volumes estimated by magnetic resonance imaging. By combining semiparametric estimation with latent variable modeling, GALAMMs allow a more realistic representation of how brain and cognition vary across the lifespan, while simultaneously estimating latent traits from measured items. Simulation experiments suggest that model estimates are accurate even with moderate sample sizes.

Parallel Sparse Computation Toolkit

This paper presents a new software framework for solving large and sparse linear systems on current hybrid architectures, from small servers to high-end supercomputers, embedding multi-core CPUs and Nvidia GPUs at the node level. The framework has a modular structure and is composed of three main components, which separate basic functionalities for managing distributed sparse matrices and executing some sparse matrix computations involved in iterative Krylov projection methods, eventually exploiting multi-threading and CUDA-based programming models, from the functionalities for setup and application of different types of one-level and multi-level algebraic preconditioners.

Inverting the discrete curl operator: A novel graph algorithm to find a vector potential of a given vector field

We provide a novel framework to compute a discrete vector potential of a given discrete vector field on arbitrary polyhedral meshes. The framework exploits the concept of acyclic matching, a combinatorial tool at the core of discrete Morse theory. We introduce the new concept of complete acyclic matchings and we show that they give the same end result of Gaussian elimination. Basically, instead of doing costly row and column operations on a sparse matrix, we compute equivalent cheap combinatorial operations that preserve the underlying sparsity structure. Currently, the most efficient algorithms proposed in literature to find discrete vector potentials make use of tree-cotree techniques. We show that they compute a special type of complete acyclic matchings. Moreover, we show that the problem of computing them is equivalent to the problem of deciding whether a given mesh has a topological property called collapsibility. This fact gives a topological characterization of well-known termination problems of tree-cotree techniques. We propose a new recursive algorithm to compute discrete vector potentials. It works directly on basis elements of 1- and 2-chains by performing elementary Gaussian operations on them associated with acyclic matchings. However, the main novelty is that it can be applied recursively. Indeed, the recursion process allows us to sidetrack termination problems of the standard tree-cotree techniques. We tested the algorithm on pathological triangulations with known topological obstructions. In all tested problems we observe linear computational complexity as a function of mesh size. Moreover, the algorithm is purely graph-based so it is straightforward to implement and does not require specialized external procedures. We believe that our framework could offer new perspectives to sparse matrix computations.

Validation of elastic wave arrival detection method based on use of sparse matrix computation

Abstract In the field of non-destructive investigations, the arrival time of the elastic waves is the important parameters to visualize the locations of the failures and/or elastic wave velocity distributions that are generally used for the soundness evaluation. The accuracy of the visualization depends on the accuracy of the arrival times, and it is demanded that the arrival times are detected automatically with required accuracy from huge amount of the elastic wave measurement data for the practical applications. In this paper, a method of the arrival time detection is proposed on the basis of the sparse matrix computation. The root mean square voltages in the noise of the wave form, is used as an index to separate the elastic waves between the noise and the signal, and the boundary of the index between the noise and the signal in the wave forms is determined in the sparse matrix. Therefore, the arrival times located in the boundary of the noise and it is detected by sparse matrix computations. The proposed method is verified on the model tests and compared with the conventional method, and it was confirmed that the proposed method detects more arrival time accurately in comparison with the conventional method.

Scalable Predictions for Spatial Probit Linear Mixed Models Using Nearest Neighbor Gaussian Processes.

Spatial probit generalized linear mixed models (spGLMM) with a linear fixed effect and a spatial random effect, endowed with a Gaussian Process prior, are widely used for analysis of binary spatial data. However, the canonical Bayesian implementation of this hierarchical mixed model can involve protracted Markov Chain Monte Carlo sampling. Alternate approaches have been proposed that circumvent this by directly representing the marginal likelihood from spGLMM in terms of multivariate normal cummulative distribution functions (cdf). We present a direct and fast rendition of this latter approach for predictions from a spatial probit linear mixed model. We show that the covariance matrix of the cdf characterizing the marginal cdf of binary spatial data from spGLMM is amenable to approximation using Nearest Neighbor Gaussian Processes (NNGP). This facilitates a scalable prediction algorithm for spGLMM using NNGP that only involves sparse or small matrix computations and can be deployed in an embarrassingly parallel manner. We demonstrate the accuracy and scalability of the algorithm via numerous simulation experiments and an analysis of species presence-absence data.

Coarse-Grained Pruning of Neural Network Models Based on Blocky Sparse Structure.

Deep neural networks may achieve excellent performance in many research fields. However, many deep neural network models are over-parameterized. The computation of weight matrices often consumes a lot of time, which requires plenty of computing resources. In order to solve these problems, a novel block-based division method and a special coarse-grained block pruning strategy are proposed in this paper to simplify and compress the fully connected structure, and the pruned weight matrices with a blocky structure are then stored in the format of Block Sparse Row (BSR) to accelerate the calculation of the weight matrices. First, the weight matrices are divided into square sub-blocks based on spatial aggregation. Second, a coarse-grained block pruning procedure is utilized to scale down the model parameters. Finally, the BSR storage format, which is much more friendly to block sparse matrix storage and computation, is employed to store these pruned dense weight blocks to speed up the calculation. In the following experiments on MNIST and Fashion-MNIST datasets, the trend of accuracies with different pruning granularities and different sparsity is explored in order to analyze our method. The experimental results show that our coarse-grained block pruning method can compress the network and can reduce the computational cost without greatly degrading the classification accuracy. The experiment on the CIFAR-10 dataset shows that our block pruning strategy can combine well with the convolutional networks.

Dimension Reduction and Clustering Models for Single-Cell RNA Sequencing Data: A Comparative Study.

With recent advances in single-cell RNA sequencing, enormous transcriptome datasets have been generated. These datasets have furthered our understanding of cellular heterogeneity and its underlying mechanisms in homogeneous populations. Single-cell RNA sequencing (scRNA-seq) data clustering can group cells belonging to the same cell type based on patterns embedded in gene expression. However, scRNA-seq data are high-dimensional, noisy, and sparse, owing to the limitation of existing scRNA-seq technologies. Traditional clustering methods are not effective and efficient for high-dimensional and sparse matrix computations. Therefore, several dimension reduction methods have been introduced. To validate a reliable and standard research routine, we conducted a comprehensive review and evaluation of four classical dimension reduction methods and five clustering models. Four experiments were progressively performed on two large scRNA-seq datasets using 20 models. Results showed that the feature selection method contributed positively to high-dimensional and sparse scRNA-seq data. Moreover, feature-extraction methods were able to promote clustering performance, although this was not eternally immutable. Independent component analysis (ICA) performed well in those small compressed feature spaces, whereas principal component analysis was steadier than all the other feature-extraction methods. In addition, ICA was not ideal for fuzzy C-means clustering in scRNA-seq data analysis. K-means clustering was combined with feature-extraction methods to achieve good results.

AMG based on compatible weighted matching for GPUs

We describe main issues and design principles of an efficient implementation, tailored to recent generations of Nvidia Graphics Processing Units (GPUs), of an Algebraic MultiGrid (AMG) preconditioner previously proposed by one of the authors and already available in the open-source package BootCMatch: Bootstrap algebraic multigrid based on Compatible weighted Matching for standard CPUs. The AMG method relies on a new approach for coarsening sparse symmetric positive definite (s.p.d.) matrices, named coarsening based on compatible weighted matching. It exploits maximum weight matching in the adjacency graph of the sparse matrix, driven by the principle of compatible relaxation, providing a suitable aggregation of unknowns which goes beyond the limits of the usual heuristics applied in the current methods. We adopt an approximate solution of the maximum weight matching problem, based on a recently proposed parallel algorithm, referred to as the Suitor algorithm, and show that it allows us to obtain good quality coarse matrices for our AMG on GPUs. We exploit inherent parallelism of modern GPUs in all the kernels involving sparse matrix computations both for the setup of the preconditioner and for its application in a Krylov solver, outperforming preconditioners available in the original sequential CPU code as well as the single node Nvidia AmgX library. Results for a large set of linear systems arising from discretization of scalar and vector partial differential equations (PDEs) are discussed.

The Sparse Polyhedral Framework: Composing Compiler-Generated Inspector-Executor Code

Irregular applications such as big graph analysis, material simulations, molecular dynamics simulations, and finite element analysis have performance problems due to their use of sparse data structures. Inspector–executor strategies improve sparse computation performance through parallelization and data locality optimizations. An inspector reschedules and reorders data at runtime, and an executor is a transformed version of the original computation that uses the newly reorganized schedules and data structures. Inspector–executor transformations are commonly written in a domain-specific or even application-specific fashion. Significant progress has been made in incorporating such inspector–executor transformations into existing compiler transformation frameworks, thus enabling their use with compile-time transformations. However, composing inspector–executor transformations in a general way has only been done in the context of the Sparse Polyhedral Framework (SPF). Though SPF enables the general composition of such transformations, the resulting inspector and executor performance suffers due to missed specialization opportunities. This paper reviews the history and current state of the art for inspector–executor strategies and reviews how the SPF enables the composition of inspector–executor transformations. Further, it describes a research vision to combine this generality in SPF with specialization to achieve composable and high performance inspectors and executors, producing a powerful compiler framework for sparse matrix computations.

On Improving Computational Efficiency of Simplified Fluid Flow Models

Single-core computational efficiency of several iterative methods for numerical solution of nonsymmetric linear systems originating from simplified fluid flow models is evaluated together with different preconditioning techniques in terms of solution behaviour and the overall rate of convergence. A previously developed simplified 3D CFD model, which involves the solution of linear systems of orders usually from units to tens of thousands, is used to generate test cases. The respective Java software application employs the Parallel Colt linear algebra library. This is to ensure that the corresponding sparse matrix computations needed in different solution and preconditioning methods are carried out in as efficient a manner as possible to minimize the influence of the computer implementation itself. Measures (warm-up phase etc.) are taken to eliminate the effects of JVM specific behaviour such as JIT compilation.

Sparse matrix computation for air quality forecast data assimilation

In this paper, we study the ensemble Kalman filter (EnKF) method for chemical species simulation in air quality forecast data assimilation. The main contribution of this paper is that we study the sparse observation data and make use of the matrix structure of the EnKF update equations to design an algorithm for the purpose of computing the analysis of chemical species in an air quality forecast system efficiently. The proposed method can also handle the combined observations from multiple chemical species together. We applied the proposed method and tested its performance in real air quality data assimilation. Numerical examples are presented to demonstrate the efficiency of the proposed computation method for EnKF updating and the effectiveness of the proposed method for NO2, NO, CO, SO2, O3, PM2.5, and PM10 prediction in air quality forecast data assimilation.