A Parallel Solver for Graph Laplacians

TL;DR

This paper introduces a parallel multigrid solver for graph Laplacians that efficiently handles large, irregular graphs in distributed memory environments, significantly outperforming existing methods.

Contribution

It presents a novel parallel aggregation and low degree elimination algorithms tailored for irregular graphs, enabling scalable solutions for massive graph Laplacian problems.

Findings

Scales to 576 processes and 1.7 billion edges

Outperforms existing parallel solvers in speed and efficiency

Effective for irregular degree graphs in distributed settings

Abstract

Problems from graph drawing, spectral clustering, network flow and graph partitioning can all be expressed in terms of graph Laplacian matrices. There are a variety of practical approaches to solving these problems in serial. However, as problem sizes increase and single core speeds stagnate, parallelism is essential to solve such problems quickly. We present an unsmoothed aggregation multigrid method for solving graph Laplacians in a distributed memory setting. We introduce new parallel aggregation and low degree elimination algorithms targeted specifically at irregular degree graphs. These algorithms are expressed in terms of sparse matrix-vector products using generalized sum and product operations. This formulation is amenable to linear algebra using arbitrary distributions and allows us to operate on a 2D sparse matrix distribution, which is necessary for parallel scalability. Our solver outperforms the natural parallel extension of the current state of the art in an algorithmic comparison. We demonstrate scalability to 576 processes and graphs with up to 1.7 billion edges.

Figures

Equations

$$A_{ij}=\{\begin{array}{ll}{w}_{ij}& \text{if }(i,j)\in E\\ 0& \text{otherwise}\end{array}$$

$$L$$

$$D_{ij}$$

$$=\frac{-\text{work}}{{\log }_{10}\Delta r},r=b-Lx,\Delta r=\frac{\Vert r_{\text{final}}\Vert }{\Vert r_{\text{initial}}\Vert }$$

$$=\frac{\text{total FLOPS}}{\text{FLOPS to compute residual on finest level}}$$

$$(Av{)}_i:=\sum _j{A}_{ij}{v}_j$$

$$(\hat{A}\oplus .\otimes v{)}_i:=\underset{j}{⨁}{\hat{A}}_{ij}\otimes v_j$$

$$\hat{A}\in {\mathbb{B}}^{m\times n},$$

$$\otimes :\mathbb{B}\times \mathbb{C}\to \mathbb{D},$$

$$a{\otimes }_{\text{elim}}c=\{\begin{array}{ll}c& \text{if }a\ne 0\\ \varnothing & \text{otherwise}\end{array}$$

$$x{\oplus }_{\text{elim}}y=\{\begin{array}{ll}x& \text{if }\text{ hash}(x)\le \text{hash}(y)\\ y& \text{if }\text{ hash}(x)>\text{hash}(y)\end{array}$$

$$L_l=\Pi \left(\begin{array}{cc}{L}_{\mathcal{FF}}& L_{\mathcal{F}\mathcal{C}}\\ L_{\mathcal{F}\mathcal{C}}^T& L_{\mathcal{CC}}\end{array}\right){\Pi }^T$$

$$L_l=\Pi \left(\begin{array}{cc}I& \\ L_{\mathcal{F}\mathcal{C}}^T{L}_{\mathcal{FF}}^{-1}& I\end{array}\right)\left(\begin{array}{cc}{L}_{\mathcal{FF}}& \\ & L_{l+1}\end{array}\right)\left(\begin{array}{cc}I& L_{\mathcal{FF}}^{-1}{L}_{\mathcal{F}\mathcal{C}}\\ 0& I\end{array}\right){\Pi }^T$$

$$L_{l+1}=L_{\mathcal{CC}}-L_{\mathcal{F}\mathcal{C}}^T{L}_{\mathcal{FF}}^{-1}{L}_{\mathcal{F}\mathcal{C}}.$$

$$P_l=\Pi \left(\begin{array}{c}-L_{\mathcal{FF}}^{-1}{L}_{\mathcal{F}\mathcal{C}}\\ I\end{array}\right).$$

$$\begin{array}{rl}{L}_l^{-1}& =\Pi \left(\begin{array}{cc}I& -L_{\mathcal{FF}}^{-1}{L}_{\mathcal{F}\mathcal{C}}\\ 0& I\end{array}\right)\left(\begin{array}{cc}{L}_{\mathcal{FF}}^{-1}& \\ & L_{l+1}^{-1}\end{array}\right)\left(\begin{array}{cc}I& \\ -L_{\mathcal{F}\mathcal{C}}^T{L}_{\mathcal{FF}}^{-1}& I\end{array}\right){\Pi }^T\\ & =P{L}_{l+1}^{-1}{P}^T+\Pi \left(\begin{array}{cc}{L}_{\mathcal{FF}}^{-1}& 0\\ 0& 0\end{array}\right){\Pi }^T\end{array}$$

$$\mathcal{F}\mathrm{smooth}(x,b)=\Pi \left(\begin{array}{cc}{L}_{\mathcal{FF}}^{-1}& 0\\ 0& 0\end{array}\right){\Pi }^{T}b.$$

$$L$$

$$x$$

$$C_{ij}$$

$$S_{ij}$$

$$\mathit{status}_{i}\coloneqq\begin{dcases*}(\mathit{Seed},i)&Vertex $i$ is a $\mathit{Seed}$\\ (\mathit{Undecided},i)&Vertex $i$ has not yet joined an aggregate\\ (\mathit{Decided},j)&Vertex $i$ is aggregated with $\mathit{Seed}$ $j$\\ \end{dcases*}$$

$$w{\otimes }_{\text{agg}}(state,i)$$

$$({state}_a,i_a,w_a)$$

$$\{\begin{array}{ll}({state}_a,i_a,w_a)& \text{if }{state}_a={state}_b\text{ \& }{w}_a\ge w_b\\ ({state}_b,i_b,w_b)& \text{if }{state}_a={state}_b\text{ \& }{w}_a<w_b\\ ({state}_a,i_a,w_a)& \text{if }{state}_a>{state}_b\\ ({state}_b,i_b,w_b)& \text{if }{state}_a<{state}_b\end{array}$$

$$R_{ij}$$

$$P$$

$$\frac{\text{nnz}(L)}{\text{TDA}\cdot \text{number of processes}}\text{TDA}:=\frac{-\text{time}}{{\log }_{10}\Delta r}$$

$$\frac{\text{nnz}(L)}{\frac{\text{solve time}}{\text{work unit}}\cdot \text{number of processes}}$$

Full text

Graph Laplacian; Unsmoothed Aggregation multigrid; Distributed Memory; Iterative Solver

††conference: Platform for Advanced Scientific Computing; July 02–04; Basel, SwitzerlandOur solver is available at https://github.com/ligmg/ligmg

1. Introduction

Graph Laplacians arising from large irregular degree graphs present a variety of challenges for typical multigrid solvers. Irregular graphs often present a broad degree distribution and relatively low diameter. Although coarse grid complexity is relatively predictable and can be controlled directly in multigrid methods for PDE problems (typically bounded degree and large diameter), that is not the case for irregular graphs. Consequently, methods such as smoothed aggregation lead to loss of sparsity in coarse grids (Livne and Brandt, 2012). Standard coarsening strategies for parallel aggregation either use some form of vertex partitioning or a maximal independent set (Gee et al., [n. d.]; Bell et al., 2012; Brannick et al., 2013). Vertices in irregular graphs cannot even be partitioned among processes (we use an edge distribution for load balancing reasons), much less partitioned into aggregates. A maximal independent set is not viable for irregular graphs or those with small diameter. For example, a single vertex with degree of one million will produce unacceptably rapid coarsening. Even if high-degree vertices are removed, social network graphs such as the Facebook graph have very small maximal independent sets (Backstrom et al., 2012; Boldi and Vigna, 2012), which also lead to unacceptably rapid coarsening. To avoid this, any viable coarsening scheme must permit adjacent root vertices (as MIS forms the smallest possible aggregates without adjacent roots). We present a new parallel algorithm for choosing aggregates that accurately captures low energy components while maintaining coarse grid sparsity. We also develop a parallel algorithm for finding and eliminating low degree vertices, a technique introduced for a sequential multigrid algorithm by Livne and Brandt (Livne and Brandt, 2012), that is important for irregular graphs.

While matrices can be distributed using a vertex partition (a 1D/row distribution) for PDE problems, this leads to unacceptable load imbalance for irregular graphs. We represent matrices using a 2D distribution (a partition of the edges) that maintains load balance for parallel scaling (Buluç and Gilbert, 2011). Many algorithms that are practical for 1D distributions are not feasible for more general distributions. Our new parallel algorithms are distribution agnostic.

2. Background

2.1. Graph Laplacian Matrix

Graphs arise in many areas to describe relationships between objects. They can be people and friendships in a social network, computers and network connections in a local area network, cities and the roads that connect them, etc. We start with a matrix representation of a graph.

A weighted graph $G=(V,E,w)$ (where $V$ are vertices, $E$ are edges and $w$ are edge weights) can be expressed as an adjacency matrix $A$:

[TABLE]

The Laplacian matrix $L$ can then be expressed as:

[TABLE]

In this paper we consider positively weighted ($w\geq 0$), undirected ($(i,j)\in E\iff(j,i)\in E$ and $w_{ij}=w_{ji}$) graphs. Graph Laplacians of this type have a couple of properties:

• •

Column and row sums are zero.

• •

Off diagonal entries are negative.

• •

Diagonal entries are positive.

• •

$L$ is symmetric positive semi-definite.

In this paper, we assume that the graph $G$ is connected. A connected graph is one where each vertex can reach every other vertex by traversing edges. For example, in a connected social network, every person has a series of friendships that link them to any other person. Assuming a connected graph simplifies our setup phase because we do not need to keep track of each component.

For a connected graph, the corresponding Laplacian matrix has a null space spanned by the constant vector. The dimension one nullspace is easier to keep track of and orthogonalize against. Given the eigensystem $Lu_{i}=\lambda_{i}u_{i}$, the eigenvalues $0=\lambda_{0}\leq\lambda_{1}\leq\lambda_{2}\leq...\leq\lambda_{n-1}$ are non-negative and real. The number of eigenvalues $\lambda_{i}=0$ is equal to the number of connected components in the graph. Because we are only considering connected graphs, $\lambda_{1}>0$. The eigenvector, $u_{1}$, associated with the second smallest eigenvalue approximates the sparsest cut in the graph. The sign of $u_{1}$ determines which side of the cut each vertex belongs to (Spielman, 2007).

Some applications of the graph Laplacian include:

• •

Graph partitioning As mentioned above, the second smallest eigenvalue of the graph Laplacian approximates the sparsest cut in the graph.

• •

Graph drawing Like graph partitioning, spectral graph drawing relies on finding eigenvalues and eigenvectors of the Laplacian to embed the graph in a two dimensional space for visualization. Each vertex, $i$, is placed at $(u_{1}(i),u_{2}(i))$. Other graph drawing techniques, such as the maxent-stress model, rely on solving the graph Laplacian (Gansner et al., 2013).

• •

PDEs on unstructured meshes Some discretizations of partial differential equations on unstructured meshes result in Laplacian matrices (Spielman, 2010). If the mesh has large variation in vertex degree, solution techniques for arbitrary graph Laplacians might outperform those designed for fixed degree or more structured meshes.

• •

Electrical flow A circuit of resistors with current inputs and sinks between them can be modeled as a graph Laplacian. Each vertex in the graph is a current input or sink, and edges correspond to resistors with resistance equal to reciprocal resistance. If $r$ is a vector where $r_{i}$ is the current input or draw, then solving $p=L^{-1}r$ gives the potential $p_{i}$ at vertex $i$ (Spielman, 2010).

See Spielman’s article on applications of the graph Laplacian for more details (Spielman, 2010).

All these applications (either directly or indirectly) require applying the action of the inverse of the graph Laplacian ($L^{-1}$). We refer to computing $x$ in $Lx=b$ as solving the graph Laplacian $L$. We refer to a method of solving $Lx=b$ as a solver. To compute small eigenvectors and eigenvalues for applications such as graph partitioning, an eigensolver must be used. Typically, these eigensolvers iterate on the inverse matrix (in this case $L^{-1}$) or on a shifted inverted matrix. The majority of time is often spent in the linear solve (i.e., computing $L^{-1}x$).

There are many different kinds of graphs: road networks, computer networks, social networks, power grids, document relations, etc. We focus our performance evaluation on social network graphs because they often provide the most difficulty for existing multigrid solvers. Structured graphs, such as PDEs discretized on a grid, provide more opportunities for a solver to cut corners. Social networks have a couple of properties that make them more difficult to solve than other graphs:

• •

Social network graphs are sparse: the number of edges is roughly a constant factor of the number of vertices.

• •

Edges in these networks follow a power-law degree distribution (these graphs are called scale-free). A small number of vertices have a large number of neighbors, whereas the rest have a relatively small degree. This presents a challenge when determining how to distribute work in parallel.

• •

Social network graphs have high connectivity. Most vertices can be reached in a small number of hops from any other vertex. This causes large fill in for techniques such as LU factorization and smoothed aggregation multigrid.

2.2. Related Work

A variety of different solvers that have been proposed for solving large sparse SPD systems. Some are general purpose solvers, whereas others are tailored specifically for graph Laplacians.

2.2.1. Direct Solvers

A variety of direct solvers capable of solving sparse systems. Some of these solvers, such as SuperLU_DIST (Li and Demmel, 2003) and MUMPS (Amestoy et al., 2001), function in distributed memory. However, none of these solvers performs well on large graph Laplacian systems because “small” vertex separators (George and Liu, 1981) do not exist. Because we are interested in very large systems, we require that our solver scales linearly in the number of nonzeros in the matrix, which these direct solvers do not.

2.2.2. Simpler Preconditioners

Iterative solvers such as conjugate gradients, coupled with simple preconditioners like Jacobi or incomplete Cholesky, sometimes perform well on highly irregular graphs. For example, Jacobi (diagonal) preconditioning is often sufficient for social network graphs with small diameter and Incomplete Cholesky may provide more robustness, but tends to exhibit poor parallel scalability.

2.2.3. Theoretical Solvers

A variety of theoretical Laplacian solvers have been proposed in literature starting with Spielman and Teng’s 2003 paper (Spielman and Teng, 2004). To our knowledge, no working implementation of this algorithm exists. Since then, many more theoretical linear solvers have been proposed. Most use either a support graph or low stretch spanning tree sparsifier as a preconditioner for a Krylov iterative solver. Several serial implementations of these ideas exist in “Laplacians.jl” a package written by Daniel Spielman (Spielman et al., 2017).

Kelner et al. later proposed a simple and novel technique with a complexity bound of $O(m\log^{2}{n}\log{\log{n}}\log{(\epsilon^{-1})})$ (Kelner et al., 2013)($n$ is the number of vertices, $m$ is the number of edges, $\epsilon$ is the solution tolerance). An implementation of this algorithm exists but appears to not be practical (as of yet) (Boman et al., 2015).

2.2.4. Practical Serial Solvers

Three practical graph Laplacian solvers have been proposed: Koutis and Miller’s Combinatorial Multigrid (CMG) (Koutis et al., 2011), Livne and Brandt’s Lean Algebraic Multigrid (LAMG) (Livne and Brandt, 2012), and Napov and Notay’s Degree-aware Rooted Aggregation (DRA) (Napov and Notay, 2016). All use multigrid techniques to solve the Laplacian problem.

Combinatorial Multigrid, like much of the theoretical literature, takes a graph theoretic approach. It constructs a multilevel preconditioner using clustering on a modified spanning tree (Koutis et al., 2011). Its main focus is on problems arising in imaging applications. The spanning tree construction and clustering presented in CMG do not lend themselves to a simple parallel implementation.

Lean Algebraic Multigrid uses a more standard AMG approach with modifications suited for Laplacian matrices. Notably, it employs unsmoothed aggregation tailored to scale-free graphs, a specialized distance function, and a Krylov method to accelerate solutions on each of the multigrid levels. These changes are not rooted in graph theory but produce good empirical results. Empirically, LAMG is slightly slower than CMG but more robust (Livne and Brandt, 2012). LAMG’s partial elimination procedure and clustering process are both inherently serial.

Degree-aware Rooted Aggregation applies a similar partial elimination technique to LAMG, except it is limited to degree 1 vertices. Its performance relies on a combination of unsmoothed aggregation based on vertex degree and a multilevel Krylov method called K-cycles(Napov and Notay, 2016). A K-cycle is a multigrid W-cycle (one can imagine a K-cycle with different cycle index, but we only consider K-cycles with a cycle index of 2) with Krylov acceleration applied at each level. DRA’s aggregation (like LAMG’s) is inherently serial and would require modifications for parallelism.

2.3. Algebraic Multigrid

The three solvers mentioned previously have one thing in common: they use some sort of algebraic multigrid method. Algebraic multigrid (AMG) is a family of techniques for solving linear systems of the form $Lx=b$. Specifically, algebraic multigrid constructs a hierarchy of approximations to $L$ using only the values in the matrix $L$. These approximations, $\{L=L_{0},L_{1},L_{2},L_{3},...\}$, are successively coarser: $\text{size}(L_{l})>\text{size}(L_{l+1})$. For each level $l$ in the hierarchy, we have a restriction operator $R_{l}$ that transfers a residual from level $l$ to level $l+1$ and a prolongation $P_{l}$ that transfers a solution from level $l+1$ to level $l$.

Algorithm 1 depicts a multigrid cycle with cycle index $\gamma$ ($\gamma=1$ is called a V-cycle, and $\gamma=2$ is a W-cycle). Repeated application of this algorithm $x^{k+1}\leftarrow\textsc{mgcycle}(0,x^{k},b)$ often converges to a solution $Lx^{*}=b$, reaching a given tolerance in a number of iterations that is bounded independent of problem size. Fast convergence depends on sufficiently accurate restriction and prolongation operators complemented by pre-smoothing and post-smoothing that provide local relaxation. Typical smoother choices are Gauss-Seidel, Jacobi, and Chebyshev iteration. Note that pre- and post-smoothing may use different smoothers or different numbers of smoothing iterations. Given a restriction and prolongation, we identify “low frequencies” as those functions that can be accurately transferred to a coarse space and back. A smoother need only be stable on such functions but must reduce the error uniformly for all “high frequencies”—those which cannot be accurately transferred.

There are two main ways of constructing $\{L_{1},L_{2},L_{3},...\}$: classical AMG and aggregation-based AMG. Classical (or Ruge-Stüben) AMG constructs $R$ and $P$ using a coarse-fine splitting: the coarse grid is a subset of the degrees of freedom (or “points”) of the fine grid (Ruge and Stüben, 1987). $P$ keeps values at coarse points and extends them via a partition of unity to neighboring fine points. Then, $R=P^{T}$. In aggregation-based AMG, degrees of freedom are clustered into aggregates. An aggregate on the fine level becomes a point on the coarse level. $R$ and $P$ usually take some weighted average of points in each aggregate to the coarse level. Aggregation-based multigrid methods typically smooth $R$ and $P$ for better performance (Vaněk et al., 1996). In both classical and aggregation-based AMG, the coarse level matrix is usually constructed via a Galerkin product: $L_{l+1}=R_{l}L_{l}P_{l}$.

2.3.1. Measuring Performance

In order to evaluate the performance of our solver relative to the state of the art, we need some performance metric. We could use runtime, but it is dependent on the implementation of the algorithms (for example, a MATLAB implementation would probably be slower than a C++ implementation). Instead, we will use work per digit of accuracy (WDA) to measure performance.

[TABLE]

WDA measures how much work is required to reduce the residual by an order of magnitude. $r$ is the residual norm, and $\Delta r$ is the change in residual norm from an initial solution to a final solution. Work is expressed in terms of the number of FLOPS required for a solve divided by that required to compute a residual on the finest level. This measure is also proportional to required memory transfers and is typically approximated in terms of the number of nonzeros in the sparse operators $L_{l},R_{l},P_{l}$. WDA only measures the efficiency of a single solve; it does not take the setup phase into account. WDA also does not account for parallel scalability.

2.3.2. AMG as a Preconditioner

AMG can be used as a solver by itself, but it often provides faster and more robust convergence when used as a preconditioner for a Krylov method. In order to solve $Lx=b$, the Krylov method is applied to $MLx=Mb$, where $M$ is a single AMG cycle, such as a V-cycle or W-cycle. The convergence of the Krylov method will depend on the spectrum of $ML$, converging in a number of iterations bounded by square root of the condition number. We use our multigrid solver as a preconditioner to conjugate gradients, or in the case of K-cycles, flexible CG (see §3.7).

3. Main Contribution

Our solver uses an unsmoothed aggregation based multigrid technique with low degree elimination. Its notable features are: 1. a matrix distribution that improves parallel performance but increases complexity of aggregation and elimination algorithms 2. a parallel elimination algorithm using this matrix distribution 3. a parallel aggregation algorithm also using this matrix distribution. We started building our solver by analyzing the performance of LAMG and its potential for parallelism.

3.1. Issues With a Parallel Implementation of LAMG

To understand what parts of LAMG we could adapt to a distributed memory setting, we ran LAMG on 110 graphs from the University of Florida Sparse Matrix Collection (Davis and Hu, 2011) while varying its parameters for cycle index, smoother and iterate recombination (we would like to also vary the aggregation routine, but LAMG’s energy unaware aggregation is not working in the MATLAB implementation). LAMG uses a cycle index of 1.5, which is half way between a V-cycle and a W-cycle. Iterate recombination is a multilevel Krylov method that chooses the optimal search direction from the previous solution guesses at each level. It is similar in nature to K-cycles (see section 3.7). Our goal was to find which parts of LAMG had the largest effect on solver performance.

Figure 1 shows the result of these tests. Visually, Gauss Seidel smoothing and iterate recombination together give the smallest WDA. Cycle index has a very small effect relative to smoother and iterate recombination.

We applied a statistical analysis to LAMG’s performance with its features on and off. The largest change in solver performance came from using Gauss-Seidel smoothing over Jacobi smoothing. On social network graphs, Gauss-Seidel smoothing often performs much better than Jacobi smoothing. Iterate recombination is the second most important factor. It helps improve performance on high diameter graphs such as road networks. The least important factor is the 1.5 cycle index. Although it does slightly improve WDA, it is much less important than Gauss-Seidel smoothing and iterate recombination.

Given the importance of Gauss-Seidel smoothing and iterate recombination, we would like to use them in a parallel implementation. However, there are a couple of challenges with using these features in a parallel implementation of LAMG:

1. (1)

The power-law vertex degree distribution can cause large work and communication imbalances. 2. (2)

LAMG’s low degree elimination is a sequential algorithm. 3. (3)

LAMG’s energy-based aggregation is a sequential algorithm. 4. (4)

Multilevel Krylov acceleration can be a parallel bottleneck. 5. (5)

W-cycles exhibit poor parallel performance. 6. (6)

Gauss-Seidel smoothing is infeasible in parallel.

3.2. 2D Matrix Distribution

Our initial implementation used a vertex distribution of the graph: each processor owns some number of rows in the Laplacian matrix. This initial implementation scaled very poorly as we increased the number of processes. One process had many more edges than all the other processes, causing a work imbalance. Using a partitioner can temporarily alleviate this problem, but in the limiting case (where each process owns a single vertex), processes with hubs will have significantly more work than those without. Research on scaling sparse matrix operations on adjacency matrices suggests that a more sophisticated distribution of matrix entries is important. We use CombBLAS, which has demonstrated the scalability and load balancing benefits of a 2D matrix distribution (Buluç and Gilbert, 2011). A 2D matrix distribution can be thought of as a partition of graph edges instead of vertices. Computational nodes (or processes) form a 2D grid over the matrix (called a processor grid). Each process is given a block of the matrix corresponding to its position in the grid (see Figure 2). Vectors can either be distributed across all processes or just processes on the diagonal. In our implementation, we found performance did not change with vector distribution, so we distribute them across diagonal entries in the processor grid. The disadvantage of a 2D distribution is that it has higher constant factors and poorer data locality.

CombBLAS expresses graph algorithms in the language of linear algebra. Instead of the usual multiplication and addition used in matrix products, CombBLAS allows the user to use custom multiplication and addition operations. We will follow the $\oplus.\otimes$ notation used by GraphBLAS (the standardization effort of CombBLAS)(Kepner et al., 2016). Here, $\oplus$ specifies the custom addition operator and $\otimes$ the multiplication. While an $m\times n$ matrix $A$ maps from $\mathbb{R}^{n}$ to $\mathbb{R}^{m}$, a generalized matrix $\hat{A}$ maps from $\mathbb{C}^{n}$ to $\mathbb{D}^{m}$, where $\mathbb{C}$ and $\mathbb{D}$ may be different.

[TABLE]

[TABLE]

Our operators are similar in structure to a semiring, but they permit different types of elements in the input matrix and the input and output vectors. If we structure our algorithms in terms of generalized matrix-vector products, we can piggyback on the proven performance of CombBLAS. However, if we cannot express all parts of an algorithm using this linear algebraic approach, we must keep in mind that no computational node has a complete view of all the edges to and from any vertex. Implementing an arbitrary vertex neighborhood operation, such as choosing the median of neighbors, would require potentially non-scalable custom communication.

3.3. Random Vertex Ordering

A 2D matrix distribution alleviates communication bottlenecks and some load balancing difficulties, but the processes responsible for diagonal blocks are often found to have many more nonzeros than typical off-diagonal blocks. For example, social network often have a couple of large degree “hubs” that are connected to many other vertices. These vertices correspond to an almost dense column and row in the graph Laplacian. Often a single process will end up with a few hubs and have 10x (or more) edges than other processes. A simple technique to better balance the workload is to randomly order vertices. This trades data locality for better load distribution. More sophisticated techniques exist for 2D matrix partitioning (Boman et al., 2013), but we found that a random distribution is sufficient for acceptable load balance. We found that random vertex ordering increased not only asymptotic parallel scalability but also performance for relatively small process counts. We apply this randomization only to the input matrix; we do not re-randomize at coarser levels.

3.4. Parallel Low-Degree Elimination

Low-degree elimination greatly reduces problem complexity in graph Laplacians, especially those arising from social networks. Like LAMG, we eliminate vertices of degree 4 or less.

The main difficulty in adapting low-degree elimination to a distributed memory system is deciding which vertices to eliminate. If we had a vertex centric distribution, each process could locally decide which of its local vertices to eliminate. However, we have a 2D edge distribution, so we will instead structure our elimination in terms of linear algebra and allow CombBLAS to do the heavy lifting.

Our algorithm for low-degree elimination is detailed in Algorithm 2. It essentially boils down to two steps. First, mark all vertices of degree 4 or less as candidates for elimination. Then, for each candidate, check whether it has the lowest hash value among all neighboring candidates. If it has the lowest hash value, it will be eliminated. The hash value is a hash of the vertex’s id. We use a hash of the id instead of the id itself in order to prevent biases that might occur when using a non-random matrix ordering.

In linear algebraic terms, choosing vertices to eliminate can be accomplished by creating a vector that marks each vertex as a candidate or not (line 2 of Algorithm 2) and then multiplying said vector with the Laplacian matrix using a custom $\otimes$ and $\oplus$ (line 2). The $\otimes_{\text{elim}}$ filters out matrix entries that are not candidates and neighbors:

[TABLE]

where $c=\varnothing$ is used to indicate that a vertex should not be considered.

The $\oplus_{\text{elim}}$ chooses the candidate with the smallest hashed id,

[TABLE]

where $\text{hash}(\varnothing)=\infty$.

We use the Laplacian $L$ so that the neighborhood of each vertex contains itself. Each entry $z_{i}$ corresponds to the neighbor of $v_{i}$ that is a candidate and has the lowest hashed id. If $z_{i}=i$, then we know that $v_{i}$ is the candidate with the smallest hash among its neighbors (line 4) and can be eliminated. Let $\mathcal{F}$ denote eliminated vertices and $\mathcal{C}$ be vertices that have not been eliminated. Following Livne and Brandt (Livne and Brandt, 2012), we express elimination in the language of multigrid. We first introduce a permutation $\Pi$ of the degrees of freedom such that

[TABLE]

which admits the block factorization

[TABLE]

in terms of the Schur complement

[TABLE]

Note that $L_{\mathcal{FF}}$ is diagonal so its inverse is also diagonal and that we can alternately express $L_{l+1}=P_{l}^{T}L_{l}P_{l}$ in terms of the prolongation

[TABLE]

Inverting the block factorization yields

[TABLE]

which is an additive 2-level method with $\mathcal{F}$-point smoother

[TABLE]

Rather than applying $L_{l+1}^{-1}$ exactly, we approximate it by continuing the multigrid cycle.

The ids of eliminated vertices are broadcast down processor rows and columns. Each process constructs entries of $L_{\mathcal{F}\mathcal{C}}$ and $L_{\mathcal{F}\mathcal{F}}^{-1}$ that depend on its local entries in $L_{l}$. These constructed entries are then scattered to the processes that own them in $P_{l}$. Alternatively, we could construct $\Pi$, $L_{\mathcal{F}\mathcal{C}}$, and $L_{\mathcal{F}\mathcal{F}}^{-1}$ explicitly and use them to build $P_{l}$.

Our candidate selection scheme is not as powerful as the serial LAMG scheme. The serial scheme will eliminate every other vertex of a chain. In the best case we do the same, but in the worst case we eliminate only one vertex if the hash values of vertices in the chain are in sequential order. To address this issue, we can run low-degree elimination multiple times in a row to eliminate more of the graph. In practice, we find one iteration is sufficient to remove most of the low degree structure.

We apply low-degree elimination before every aggregation level and only if more than 5% of the vertices will be eliminated.

3.5. Parallel Aggregation

Our parallel aggregation algorithm (Algorithm 3) uses a strength of connection metric, $S$, to determine how to form aggregates. It indirectly determines how likely any two vertices will be clustered together. We use the affinity strength of connection metric proposed by Livne and Brandt in the LAMG paper (Livne and Brandt, 2012). To construct the strength of connection matrix, $S$, we smooth four vectors three times each. The relevant parts of the smoothed vectors are broadcast down communication grid rows and columns. Each process then constructs its local part of the matrix using the affinity metric

[TABLE]

[TABLE]

where $A$ is the adjacency matrix. The smoothing we use is three iterations of Jacobi smoothing. The total cost of creating $S$ is 4 vectors $*$ 3 smoothing iterations, for 12 matrix-vectors multiplies total. Currently, we smooth each vector separately and do not exploit any of the parallelism available in smoothing multiple vectors together.

The construction of $C_{ij}$ is entirely local because $C$ has the same distribution as $A$. Constructing $S$ requires communication along processor rows and columns to find the largest nonzero of $C$ for each column and row in the matrix. Note that $S$ has [math] diagonal.

Our aggregation algorithm uses a voting scheme in which each vertex votes for which of its neighbors it would like to aggregate with. A vector $status$ contains a state $\in\{\mathit{Seed},\mathit{Undecided},\mathit{Decided}\}$ and an index for each vertex,

[TABLE]

Initially, $\mathit{status}_{i}=(\mathit{Undecided},i)$. In each voting iteration, each $\mathit{Undecided}$ vertex either aggregates with a neighboring $\mathit{Seed}$ (and becomes $\mathit{Decided}$) or votes for a neighboring $\mathit{Undecided}$ vertex to become a $\mathit{Seed}$. If a vertex is voted for enough times, it will turn into a $\mathit{Seed}$. The strength of connection matrix determines which neighboring vertex is aggregated to or voted for. For the first round of aggregation, we only consider very strong connections ($>0.5$). We gradually reduce this bound in subsequent iterations (Livne and Brandt apply a similar technique in LAMG (Livne and Brandt, 2012)). The vertices’ choice of neighbor is expressed as a matrix-vector product $S\>{\oplus_{\text{agg}}.\otimes_{\text{agg}}}\>\mathit{status}$ with

[TABLE]

[TABLE]

where $\mathit{Seed}>\mathit{Undecided}>\mathit{Decided}$. Note that the input vector contains pairs, whereas the output vector contains 3-tuples.

The votes for each vertex are tallied using a sparse reduction that has the same communication structure as a matrix-vector product. Our implementation uses an allreduce because it has lower constants and is not a bottleneck. Using the tallied votes, we update the status vector with new roots and aggregates. The vote counts are persisted across voting iterations so that vertices with low degree may eventually become seeds. We choose to do 10 voting iterations, and we convert $\mathit{Undecided}$ vertices to $\mathit{Seed}$s if they receive 8 or more votes. Both these numbers are arbitrary. In practice, we find that performance is not sensitive to increasing or decreasing these constants by moderate amounts.

The result of aggregation is a distributed vector $v$, where vertex $i$ is part of aggregate $v_{i}$. We perform a global reordering so that aggregates are numbered starting at [math]. We construct $R$ by inserting $R_{v_{j}j}=1$, where $j$ is in the locally owned portion of $v$, then scattering to a balanced 2D distribution. This 2D distribution is similar to $L$, except each process has a rectangular local block instead of a square one.

[TABLE]

3.6. Smoothing

In general, Gauss-Seidel smoothing is the best performing smoother on graph Laplacians (section 3.1 provides more details on Gauss-Seidel vs Jacobi performance). However, its parallel performance on graph Laplacians is very poor. Most processes have an overwhelming amount of connections that reference values outside of the local block of the matrix, and the graph cannot be colored with a reasonable number of colors. Our resulting choice of smoother is Chebyshev/Jacobi smoothing (Adams et al., 2003) because it is stronger than (weighted) Jacobi with equivalent parallel performance. Instead of applying $k$ iterations of a smoother, we use one application of degree $k$ Chebyshev smoothing. We choose our lower and upper bounds because $.3$ and $1.1$ times the largest eigenvalue as estimated by $10$ Arnoldi iterations (we include these iterations in our setup cost).

3.7. K-cycles

We would like to include some form of multilevel Krylov acceleration because it improves solver robustness (see section 3.1 for more details). We implemented K-cycles as described by Notay and Vassilevski in “Recursive Krylov-based multigrid cycles” (Notay and Vassilevski, 2008) and used in Napov and Notay’s DRA (Napov and Notay, 2016). At each aggregation level in our multigrid hierarchy, we perform a number of Flexible Conjugate Gradient (FCG) (Notay, 2000) iterations using the rest of the hierarchy as a preconditioner. We do not apply Krylov acceleration to elimination levels because they exactly interpolate the solution from the coarse grid.

4. Numerical Results

Our solver uses a V-cycle or K-cycle with one iteration of degree 2 Chebyshev smoothing before restriction and one iteration of degree 2 Chebyshev smoothing after prolongation. The V-cycle is used as a preconditioner for Conjugate Gradient (FCG when using K-cycles). Our solver is implemented in C++ and uses CombBLAS (Buluç and Gilbert, 2011) for sparse linear algebra and PETSc (Balay et al., 2016; Balay et al., 1997) for Chebyshev smoothing, eigenvalue estimation, and Krylov methods.

Our numerical tests were run on NERSC’s Edison and Cori clusters. Edison is a Cray XC30 supercomputer with 24 “Ivy Bridge” Intel Xeon E5-2695 v2 cores per node and a Cray Aries interconnect. Cori is Cray supercomputer with 36 “Haswell” Intel Xeon E5-2698 v3 cores per node and a Cray Aries interconnect. For each test, we run four MPI processes per physical node in order to obtain close to peak bandwidth.

We solve to a relative tolerance of $10^{-8}$. We use a random right hand side with the constant vector projected out. We also tested with a right hand side composed of low eigenmodes but did not notice any difference in performance compared to a random right hand side. The coarsest level size is set to have no more than 1000 nonzeros in $L$. Our solver uses somewhere from 20 to 40 levels depending on the problem. In order to ensure enough work per process, we use a smaller number of processors on coarser levels if the amount of work drops below a threshold (if $\text{nnz}(L)/10000<$ number of processors). The coarsest level always ends up on a single process.

4.1. Comparison to Serial

We compared our parallel solver to Livne and Brandt’s serial LAMG implementation in MATLAB (Livne and Brandt, 2012, 2011). It is hard to fairly compare single threaded performance between our solver and LAMG. One or the other could be better optimized, or choice of programming language might make a difference. To provide a fair comparison, we measure the work per digit of accuracy of each solver.

In general, our method exhibited poorer convergence factors than the serial LAMG implementation. We expect that our solver performs worse because we have made multiple concessions for parallel scalability. These changes are:

1. (1)

No energy-based aggregation 2. (2)

Chebyshev smoothing (versus Gauss-Seidel) 3. (3)

No 1.5 cycle index 4. (4)

No multilevel Krylov acceleration (for V-cycles)

Figure 3 shows the performance (measured in terms of WDA) of serial LAMG and our solver. The fourth line shows LAMG with all of its parallel-unfriendly features enabled. The third line shows LAMG without these features (but using LAMG’s standard aggregation and elimination). The first line is our solver without K-cycles and the second line with K-cycles. Our solver has a higher median WDA and variance than LAMG with all features enabled. Our solver is not as robust and has more outliers. Most of these outliers are road networks. A couple of graphs have a fairly high WDA with our solver but are solved quickly by LAMG. This is expected because we have made concessions in order to achieve parallel performance. However, our solver performs much better than LAMG with parallel friendly features (Jacobi smoothing and no recombination).

Also interesting to note is the small difference in WDA of our solver with and without K-cycles. K-cycles have low variance, and hence are more robust, but median performance is not much improved. However, this small gain in WDA is not worth the high parallel cost of K-cycles (as seen in Figure 4). On the hollywood graph, K-cycles are clearly slower and scale worse than V-cycles. K-cycles need inner products on every level of the cycle (as with Krylov smoothers), which require more communication. The W-cycle structure of K-cycles also causes a parallel bottleneck. The coarsest level is visited $2^{10}-2^{11}$ times, which results in lots of sequential solves and data redistributions. Because K-cycles appear worse than V-cycles (especially when compared to the marginal robustness they provide), we use V-cycles by default in our solver. All following performance results for our solver use V-cycles.

4.2. Strong Scaling

To measure the scalability of our approach, we measured strong scaling on four real world social network graphs using up to 576 processors. We choose these four graphs because they are some of the largest real world irregular graphs that we can find and are infeasible to solve on a single process. Many of the graphs in the 110 we use for serial experiments are small enough that solving them on a single process is fast enough.

Our largest graph, com-friendster, (from the Stanford Large Network Dataset Collection (Leskovec and Krevl, 2014; Yang and Leskovec, 2015)) has 3.6 billion nonzero entries in its Laplacian matrix. Solving a Laplacian of this size on a single node is infeasible (even if it could fit in memory, the solve time would be much too long). Figure 5 shows the efficiency, measured as

[TABLE]

versus solve time (where TDA is time per digit of accuracy). We do not use a percentage for efficiency because the choice of base efficiency makes it difficult to compare different solves on different graphs. A horizontal line indicates that the solver is scaling optimally (increasing machine size moves to the left). The largest graph (com-friendster) appears to have better scaling than the smaller graphs. com-orkut and hollywood have a much higher efficiency than com-friendster and com-lj because com-orkut and hollywood have a lower WDA. If we normalize by WDA on both axes, so the y axis becomes

[TABLE]

(this is similar to a “per iteration” metric), then we get Figure 6. Normalizing by WDA compares efficiency independent of how difficult it is to solve each graph. The smaller graphs are solved faster but scale poorly with increasing number of processes. com-friendster takes longer to solve because it is larger and efficiency is somewhat lower (likely due to poor cache behavior with the random ordering) but scales better than the smaller graphs. The poor scaling for smaller graphs is explained by having less work per process, resulting in a solve phase dominated by communication.

Our high setup cost (relative to solve time) also increases the number of repeated solves necessary. The majority of the setup phase is spent in the finest level triple product ($P^{T}LP$), which is handled by a generic CombBLAS matrix-matrix product using $+.*$. A matrix-matrix product that exploits the special structure of $R$ and $P$ would lower setup times. Still, these setup times are reasonable and can be amortized over multiple solve phases (when possible).

For a complete view of the scalability of our solver, we would like to measure weak scaling. However, it is difficult to find a fair way to measure weak scaling on graphs. We could generate a series of increasingly large random graphs, but in practice, most solvers perform much better on random graphs than real world graphs.

5. Conclusions

We have presented a distributed memory graph Laplacian solver for social network graphs. Our solver uses a 2D matrix distribution combined with parallel elimination and unsmoothed aggregation to demonstrate parallel performance on up to 576 processes. Our novel aggregation algorithm can handle arbitrary matrix distributions and forms accurate aggregates while controlling coarse grid complexity on a variety of irregular graphs. The parallel elimination algorithm presented uses generalized matrix products to find elimination candidates independent of matrix distribution. To our knowledge, this is the first distributed memory multigrid solver for graph Laplacians. It enables solving graph Laplacians that would be infeasible to solve on a single computer.

Our solver’s robustness is behind that of serial LAMG (as outlined in section 4.1) but outperforms the natural parallel extensions of LAMG (No iterate recombination, Jacobi smoothing, and cycle index 1). Further improvement is a topic for future research. Our solver performs well on social network graphs, but WDA is sometimes high for more regular graphs such as road networks. Many such graphs also admit a vertex partition with low edge cut, in which case our 2D distribution may be unnecessary. Extending our solver to handle scaled graph Laplacians (used by some applications) and unsymmetric graphs are other areas for future research.

Acknowledgements.

This material is based upon work supported by the Sponsor U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under Award Number Grant #DE-SC0016140. This research used resources of the Sponsor National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. Grant #DE-AC02-05CH11231.