An Empirical Comparison of Graph Laplacian Solvers

Solving Laplacian linear systems is an important task in a variety of practical and theoretical applications. This problem is known to have solutions that perform in linear times polylogarithmic work in theory, but these algorithms are difficult to implement in practice. We examine existing solution techniques in order to determine the best methods currently available and for which types of problems are they useful. We perform timing experiments using a variety of solvers on a variety of problems and present our results. We discover differing solver behavior between web graphs and a class of synthetic graphs designed to model them. 1 Graph Laplacians and SDD Systems The Laplacian matrix of a weighted, undirected graph G is defined as LG = DG−AG, where DG is the diagonal matrix containing the sum of incident edge weights and AG is the weighted adjacency matrix. LG is a symmetric, positive semidefinite, diagonally dominant M-matrix, with a nullspace containing the constant vector. Solving linear systems on the Laplacians of structured graphs, such as two and three dimensional meshes, has long been important in finite element analysis (with applications in electrical and thermal conductivity, and fluid flow modeling [5]) and image processing (with applications in image segmentation, inpainting, regression, and classification [22]). More recently, solving linear systems on the Laplacians of large graphs without mesh-like structure has ∗Sandia is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000. †Partially supported by Intel Corp. and Contract #DEAC02-05CH11231 from DOE Office of Science. emerged as an important computational task in network analysis (with applications to maximum flow problems [8], graph sparsification [26], and spectral clustering [19].) Interest in solving Laplacians is closely associated with interest in solving slightly more general symmetric diagonally dominant (SDD) matrices. A matrix A is SDD if A = A and ai,i ≥ ∑ i6=j |ai,j |. One way of solving these problems is to first reduce them to a Laplacian matrix [3, 13] and then solve them with a Laplacian solver. Some of the Laplacian linear systems in the applications above derive from SDD linear systems. Since the reduction from an SDD matrix to a graph Laplacian is asymptotically cheap, the theoretically best SDD solvers rely on Laplacian solvers. In the last few decades the theoretical computer science community has done significant work on asymptotically efficient Laplacian solvers. Spielman and Teng [27] showed how to solve these problems in linear times polylogarithmic work, later improved upon by Koutis, Miller, and Peng [21], but these algorithms have complicated descriptions and suspected large constants, preventing their practical implementation. An algorithm proposed by Kelner et al. [18] has the advantage of a simple algorithm description, but initial experimental results [4, 16] suggest it will have trouble competing in practice. To our knowledge there is not a comprehensive experimental study of available solvers. Our goal is to provide an initial understanding of available solvers, and to determine where ideas from more theoretical solvers could be leveraged to make improvements in practice. We believe it is important to procure a set of test problems from relevant applications, and a set of performance metrics for evaluating current and future Laplacian solver performance on these test problems. We hope to understand how Laplacian solver performance varies depending on the category of test problem. This is ongoing work and we envision this paper as a living document that will iteratively improve to incorporate new solvers, better test problems, and more useful performance analysis. To that end, the authors welcome feedback on how to improve our understanding of Laplacian solvers. 2 Test Graphs/Matrices For this first round of experiments we only use graphs/matrices that come to us as graph Laplacians, and adjacency matrices that can be easily converted to graph Laplacians. We use graphs with between 10 thousand edges and 10 million edges. Larger graphs are left for future parallel experiments. We use only the largest connected component of each graphs to ensure nullspace dimension of one. We put together four sets of test graphs/matrices in these experiments from three different sources. The first source is the University of Florida (UF) sparse matrix collection [10]. We divide these graphs into two categories, 2D/3D mesh-like structural problems, and graphs with irregular degree distributions that come from web or citation networks. For the most part these details are specified in the UF collection. There is a set of graphs from the DIMACS10 challenge subcollection which appear to have 2D/3D structure even though they are not classified as such that we include under 2D/3D graphs. Only unweighted adjacency matrices from the UF collection were used and converted to Laplacians. In some cases directed graphs were symmetrized by adding the transpose. The second set of graphs are Block Two-Level Erdös Rényi (BTER) graphs generated with the Feastpack graph generator [20, 25]. This generator was designed to produce graphs with degree distributions and clustering properties similar to the real world networks included in the irregular UF set. A graph is generated based on the following input parameters: approximate number of vertices, average degree, maximum degree, maximum clustering coefficient, and global clustering coefficient. It gives a choice between a generalized log normal distribution and a discrete power law distribution; we always choose the former. We attempted to generate graphs with realistic parameters as described by Kolda et al. [20], modeling them after the graphs we use from the UF irregular test set. We do this to compare performance between the original graphs and the synthetic graphs designed to model them. The third set of graphs arises from image segmentation problems. An image can be considered as a graph by treating pixels as vertices, with positive weight edges between them measuring dissimilarity. We use Felzenszwalb and Huttenlocher’s image segmentation code [12] to produce graphs from images. Graphs are named after the subject of the original image: cats, cities, food, and space images. These graphs are all 9 point meshes. Their interest for our experiments is not their structure, but rather their edge weights as this is our only test set of weighted graphs.