diGRASS: Directed Graph Spectral Sparsification via Spectrum-Preserving Symmetrization

Graph sparsification aims at finding a subgraph (sparsifier) that has the same set of vertices but much fewer edges than the original graph. There are several types of graph sparsifiers proposed for undirected graphs. Graph spanners [2, 3, 4, 34] are trying to preserve the pair-wise shortest-path distance between the original graph and the sparsifier. Cut sparsifiers [7, 21] are targeting preserving the cut values between cuts. Spectral sparsification methods preserve the graph spectral (structural) properties, such as distances between vertices, effective resistances, cuts in the graph, as well as the stationary distributions of Markov chains [12, 13, 40]. Therefore, spectral graph sparsification is a much stronger notion than cut sparsification, and more spectral related sparsification methods are proposed in recent years, such as spectral preservation of pseudoinverse for the graph Laplacian [29] and linear-sized sparsifier [28].

When dealing with directed graphs, it’s natural to convert directed graphs into undirected ones so that existing undirected graph algorithms can be subsequently leveraged. The related transforming procedures are called symmetrization methods. We will review three existing graph symmetrization methods, including the \(A+A^\top\) ² symmetrization, bibliometric symmetrization methods, and the random-walk symmetrization.

— \(\mathbf {A+A^\top }\) symmetrization simply ignores the edges’ directions, which is the simplest and most efficient way for directed graph symmetrization. However, edge directions may play an important role in directed graphs. As shown in Figure 1, edges \((8,1)\) and \((4,5)\) seem to have the equal importance in the symmetrized undirected graph \(\mathbf {A+A^\top }\) . However, in the original directed graph, edge \((8,1)\) is much more important than edge \((4,5)\) , since removing edge \((8,1)\) will lead to the loss of more connections in the directed graph. For example, removing edge \((4,5)\) will only affect the walks from node 4 to any other nodes as well as walks from any other nodes to node 5. However, if we remove edge \((8,1)\) in the directed graph, it will affect the walks from node 8 to any other nodes and the walks to node 1; there will also be no access from nodes \(5 ,6, 7\) , and 8 to nodes \(1, 2, 3\) , and 4.

— Bibliographic symmetrization [37] adopts \(\mathbf {AA^\top +A^\top A}\) as the adjacency matrix after symmetrization to take the in-going and out-going edges into consideration. However, it cannot be scaled to large-scale graphs since it will create much denser undirected graphs after symmetrization. Also, disconnected graphs can be created due to the \(\mathbf {AA^\top +A^\top A}\) symmetrization, as shown in Figure 1.

— Random-walk symmetrization [11] is based on random walks and allows normalized cut to be preserved after symmetrization. This is also a new symmetrization approach used in recent work for defining the Laplacian matrix of directed graphs. When defining the Laplacian matrix, we can apply the random walk for aperiodic graphs, or lazy random walk scheme for periodic graphs. In [11], Cheeger’s inequality has been extended to directed graphs and plays a significant role in spectral analysis of directed graphs. It connects Cheeger constant (conductance) with spectral properties (eigenvalues of the graph Laplacian) of a graph. It also provides the bound for the smallest eigenvalue of the directed graph Laplacian. However, the related theoretical results can only be applied to strongly-connected aperiodic directed graphs, which are rare to find in real-world applications.

Refs. [13] and [12] expanded the scope and work on not only strongly-connected graphs but also Eulerian graphs. However, there are obvious limitations with this approach. For example, a random graph needs to be converted into an Eulerian graph via an Eulerian scaling procedure by introducing additional edges, changing the directions of the edges, or reweighing the edges, which may jeopardize the original graphs’ spectral properties [13]. In addition, the Eulerian scaling is very timing consuming for large-scale graphs. Lastly, even though the complexity of algorithms is nearly linear-time, it is still not fast in practice for different applications, such as solving asymmetric linear systems, computing the stationary distribution of a Markov chain or computing expected commute time in a directed graph, and so on.

In [9], the authors design cut sparsifiers and sketches for directed graphs. Cut is a property that connects between undirected and directed graphs. How to construct cut sparsifiers for directed graphs really depends on cut balance, which is the ratio between incoming and outgoing edges in any given cut.

Undirected graph. Consider a weighted, undirected graph \({G}=({V},{E}, \omega)\) with \({n = |V|}\) and \({m = |E|}\) , where \({V}\) denotes a set of vertices, \({n}\) denotes the number of vertices, \({E}\) denotes a set of edges, \({m}\) denotes the number of edges, and \(\omega\) denotes a weight function that assigns a positive weight to each edge. The adjacency matrix of graph \({G}\) can be defined as follows:The Laplacian matrix can be computed bywhere \(\mathbf {D_G}\) is a diagonal matrix with elements \(D_{G}(p,p)=\sum \nolimits _{t\ne p}\omega (p,t)\) . For any real vector \(\mathbf {x}\in {\mathbb {R} ^n}\) , the Laplacian quadratic form of graph G is defined as \(\mathbf {{{x^\top }L_{G} x}} = \sum \nolimits _{\left({p,q} \right) \in E} {{\omega ({p,q})}{{({x(p) - x(q)})}^2}}\) . Recall that the undirected graph Laplacian is defined in Equation (2). Alternatively, the undirected graph Laplacian can also be written as [38]where matrix \(\mathbf {W}\) is the diagonal matrix with \(W(p, p)\) to be the node degree for node p and matrix \(\mathbf {B}\) is shown as

Directed graph. Consider a directed graph \(G=(V,E_G,w_G)\) with V denoting the set of vertices, \(E_G\) representing the set of directed edges, and \(w_G\) denoting the associated edge weights. Let \(n=|V|\) , \(m=|E_G|\) be the size of node and edge set. In the following, we denote the diagonal matrix by \(\mathbf {D_G}\) with \({D_G}(i,i)\) being equal to the (weighted) outdegree of node i, as well as the adjacency matrix of G by \(\mathbf {A_G}\) :Then the directed Laplacian matrix can be constructed as follows [13, 14]:

The directed graph Laplacian matrix can also be constructed as \(\mathbf {L_G=B^\top W C}\) , where \(\mathbf {W}\) and \(\mathbf {B}\) are defined the same as above, while matrix \(\mathbf {C}\) is a signed edge-vertex incidence (injection) matrix defined as follows:

For better illustration, we have summarized the most-frequently used symbols in our article in Table 1. It can be shown that any directed (undirected) graph Laplacian constructed using Equation (6) will satisfy the following properties: (I) Each column (and row) sum is equal to zero; (II) All off-diagonal elements are non-positive; (III) The Laplacian matrix is asymmetric (symmetric) and indefinite (positive semidefinite).

Spectral sparsifier was first introduced by Spielman and Teng [40]. Given an undirected graph with n vertices and m edges, a nearly-linear time algorithm was introduced for building \((1\pm \epsilon)\) spectral sparsifiers with \(O(n \log n/\epsilon ^2)\) edges in [38]. S is said to be a \((1\pm \epsilon)\) spectral sparsifier of G if the following inequality holds for any \(x\in {\mathbb {R} ^n}\) :where \(\mathbf {L_{G}}\) and \(\mathbf {L_S}\) denote the symmetric diagonally dominant (SDD) Laplacian matrices of graphs G and S, respectively. The key to the analysis of spectral sparsifier S, which is the improved construction of Equation (8), is to observe the following equation [6]:where \(\mathbf {x} \in \mathbb {R}^{n}\) . Relative condition number can be defined as \(\kappa (\mathbf {L_G},\mathbf {L_S})\le \sigma ^2\) , implying that a smaller relative condition number or \(\sigma ^2\) corresponds to a higher (better) spectral similarity between two graphs.

For directed graphs, the subgraph S can be considered spectrally similar to the original graph G if the condition number or the ratio between the largest and smallest singular values of \(\mathbf {L_S^{+} L_G}\) is close to 1, where \(\mathbf {L_S^+}\) denotes the Moore–Penrose pseudoinverse of \(\mathbf {L_S}\) . Spectral sparsification of directed graphs is equivalent to finding an ultra-sparse subgraph S such that the condition number of \(\mathbf {(L_S^{+}L_G)^\top (L_S^{+}L_G)}\) is small enough. Note that the singular values of \(\mathbf {L_S^{+} L_G}\) are the square roots of eigenvalues of \(\mathbf {(L_S^{+}L_G)^\top (L_S^{+}L_G)}\) , and \(\mathbf {(L_S^{+}L_G)^\top (L_S^{+}L_G)}\) can be written into \(\mathbf {L_G^\top (L_S L_S^\top)^{+}L_G}\) . Although \(\mathbf {L_G^\top (L_S L_S^\top)^{+}L_G}\) is not equal to \(\mathbf {(L_S L_S^\top)^{+}(L_GL_G^\top)}\) , they do share the same eigenvalues under special conditions according to the following theorem [20]:

Based on Therorem 1, \(\mathbf {L_G^\top (L_S L_S^\top)^{+}L_G}\) and \(\mathbf {(L_S L_S^\top)^{+}(L_GL_G^\top)}\) will share the same eigenvalues. Under this condition, spectral sparsification of directed graphs is equivalent to finding an ultra-sparse subgraph S such that the condition number of \(\mathbf {(L_S L_S^\top)^{+}(L_GL_G^\top)}\) is small enough. Theorem 2 shows both \(\mathbf {L_GL_G^\top }\) and \(\mathbf {L_SL_S^\top }\) are the Laplacian matrices of some undirected graphs.

\(\mathbf {L_GL_G^\top }\) symmetrization is a novel spectrum-preserving Laplacian symmetrization procedure for converting directed graphs into undirected ones. On the other hand, \(\mathbf {L_G^\top L_G}\) does not work for this purpose since \(\mathbf {L_G^\top L_G}\) does not correspond to the Laplacian of an undirected graph. Since the row sum of \(\mathbf {L_G}\) is not zero, the row sum of the \(\mathbf {L_G^\top L_G}\) will not be zero shown as follows:Although \(\mathbf {L_G^\top L_G}\) is a PSD matrix, the all-one vector is not its null space, and existing methods for spectral sparsification of undirected graphs [5, 43] cannot be exploited for sparsifying directed graphs.

In this section, we prove the existence of nearly-linear-sized spectral sparsifier for directed graphs under the condition that their corresponding undirected graphs (obtained through the proposed Laplacian symmetrization scheme) only contain non-negative edge weights.

Theorem 4 proves that there exists an undirected graph \({S_u}\) and the connection between graph \({S_u}\) and the original directed graph G. The next key step is to show that \(\mathbf {L_{S_u}}\) can be factorized into the product of a directed graph Laplacian \(\mathbf {L_S}\) and its transpose. Also, Theorem 4 does not immediately imply a sparse structure in \(\mathbf {L_S}\) . To prove the existence of nearly-linear-sized spectral sparsifier for a directed graph, we further assume the undirected graph \({G_u}\) obtained through the proposed Laplacian symmetrization only contains edges with non-negative weights. Next, the following Lemma 5 can be exploited to prove the existence of nearly-linear-sized \(\mathbf {L_S}\) based on Equation (22):

Reference [25] provides a nearly-linear time algorithm for constructing \(\mathcal {L}\) using the sparsified Cholesky factorization method, which is a process of constructing clique structure of the Schur complement. More importantly, [25] demonstrates that \(\mathcal {L}\) is a lower diagonal matrix which always corresponds to a feed-forward directed graph. As a result, we can conclude that given a Laplacian matrix of an undirected graph with non-negative edge weights, it can be always factorized into the product of a nearly-linear-sized directed Laplacian matrix and its transpose through the \(LDL^T\) decomposition (sparsified Cholesky factorization).

Combining Theorem 4 and Lemma 5 will allow us to prove the following main theorem.

However, the above theorem has its own limitations. For example, when comparing with the works of [11, 13, 14] which can preserve the cut in the sparsifier, ours cannot.

Motivated by the recent research on low-stretch spanning trees [1, 17] and spectral perturbation analysis [18, 43] for nearly-linear-time spectral sparsification of undirected graphs, we propose a practically-efficient algorithm for sparsifying general directed graphs by first constructing the initial subgraph sparsifiers of directed graphs through the following steps:

As aforementioned, when the condition number of \(\mathbf {L_{S_u}^+L_{G_u}}\) is small, the condition number of \(\mathbf {L_{S}^+L_{G}}\) will be small, which represents that graph S is a good spectral sparsifier for graph G. To this end, we will exploit the following spectral perturbation analysis framework for computing spectral sensitivity of each off-subgraph edges. For the generalized eigenvalue problemlet matrix \(\mathbf {V=[v_1, \ldots , v_n]}\) . Then \(\mathbf {v_i}\) and \(\mathbf {\lambda _i}\) can be constructed to satisfy the following orthogonality requirement:Consider the following first-order generalized eigenvalue perturbation problem:where a small perturbation \(\delta \mathbf {L_{S_u}}\) in \(\mathbf {L_{S_u}}\) is introduced, leading to the perturbed generalized eigenvalues and eigenvectors \(\lambda _i+\delta \lambda _i\) and \(\mathbf {v_i+\delta v_i}\) . By only keeping the first-order terms, Equation (29) becomesLet \(\mathbf {\delta v_i=\sum _j \psi _{i,j}v_j}\) , then Equation (30) can be expressed asBased on the orthogonality properties in Equation (28), multiplying \(\mathbf {v_i}\) to both sides of Equation (31) results inwhich further leads toThen the task of spectral sparsification of general (un)directed graphs will require to recover as few as possible off-subgraph edges to the initial directed subgraph S such that the largest eigenvalues, or the condition number of \(\mathbf {L^+_{S_u}} \mathbf {L_{G_u}}\) can be dramatically reduced. Expand \(\mathbf {\delta L_{S_u}}\) with only the first-order terms aswhere \(\mathbf {\delta L_S}= {{w_G(p,q)}\mathbf {e_{p,q}}}\mathbf {e_p}^\top\) for \({(p,q)\in E_G\setminus E_S}\) , \(\mathbf {e_p}\in \mathbb {R}^n\) denotes the vector with only the p-th element being 1 and others being 0, and \(\mathbf {e_{p,q}}=\mathbf {e_p}-\mathbf {e_q}\) . The spectral sensitivity for each off-subgraph edge \((p,q)\) can be expressed asIt is obvious that Equation (35) can be leveraged to rank the spectral importance of each off-subgraph edge. Consequently, spectral sparsification of general graphs can be achieved by only recovering a few dissimilar off-subgraph edges with large spectral sensitivity values. In this work, the following method based on t-step power iterations is proposed for efficient computation of dominant generalized eigenvectorswhere \(\mathbf { h_0}\) denotes a random vector. When the number of power iterations is small (e.g., \(t\le 3\) ), \(\mathbf {h_t}\) will be a linear combination of the first few dominant generalized eigenvectors corresponding to the largest few eigenvalues. Then the spectral sensitivity for the off-subgraph edge \((p,q)\) can be approximately computed byThe computation of \(\mathbf {h_t}\) through power iterations requires solving the linear system of equations \(\mathbf {L_{S_u} x=b}\) for t times. Note that only \(\mathbf {L_{S_u}}\) needs to be explicitly computed for generalized power iterations. The Lean Algebraic Multigrid (LAMG) [30] solver is leveraged for computing \(\mathbf {h_t}\) , which can handle undirected graphs with negative edge weights and has an empirical \(O(|E_{S_u}|)\) complexity for solving Laplacian matrices \(\mathbf {L_{S_u}}\) .

The setup phase of LAMG contains two main steps [30], as shown in Figure 3. First, a nodal elimination procedure is performed to eliminate disconnected and low-degree nodes. Next, a node aggregation procedure is applied for aggregating strongly connected nodes according to the following affinity metric \(c_{uv}\) for nodes u and v:where \(\mathbf {X} = (\mathbf {x^{(1)}, \dots , x{^{(K)}}})\) is computed by applying a few Gauss–Seidel (GS) relaxations using K initial random vectors to the linear system equation \(\mathbf {L_{S_u} x}=0\) . Let \(\mathbf {\tilde{x}}\) represent the approximation of the true solution \(\mathbf {x}\) after applying several GS relaxations to \(\mathbf {L_{S_u} x}=0\) . Due to the smoothing property of GS relaxation, the latest error can be expressed as \(\mathbf {e_s \,=\,x-\tilde{x}}\) , which will only contain the smooth components of the initial error, while the highly oscillating modes will be effectively damped out [8]. Nodes u and v are considered strongly connected to each other if \(\mathbf {X(u,:)}\) and \(\mathbf {X(v,:)}\) are highly correlated for all the K test vectors (or a larger \(c_{uv}\) value), which thus should be aggregated to form a coarse level node.

Once the multilevel hierarchical representations of the original graph (Laplacians) have been created, algebraic multigrid (AMG) solvers can be built and subsequently leveraged to solve large Laplacian matrices efficiently.

The proposed spectral sparsification algorithm will first sort all off-subgraph edges according to their spectral sensitivities in descending order \(({p_1,q_1}),({p_2,q_2}),...\) and then select top few off-subgraph edges to be recovered to the initial subgraph. To avoid recovering redundant edges into the subgraph, it is indispensable to check the edge similarities: only the edges that are not similar to each other will be added to the initial sparsifier. To this end, we exploit the following spectral embedding scheme for distinguishing off-subgraph edges leveraging approximate dominant generalized eigenvectors \(\mathbf {h_t}\) computed by Equation (36):where \((p,q_k)\) are the directed edges sharing the same head with \((p,q)\) but different tails. Then the proposed scheme for checking the spectral similarity of two off-subgraph edges will include the following steps:

Step 1: Perform t-step power iterations with \(r=O(\log n)\) initial random vectors \(\mathbf {h^{(1)}_0,\ldots ,h^{(r)}_0}\) to compute r approximate dominant generalized eigenvectors \(\mathbf {h^{(1)}_t,\ldots ,h^{(r)}_t}\) ;

Step 2: For each edge \((p,q)\) , compute an r-dimensional spectral embedding vector \(\mathbf {s_{p,q}} \in \mathbb {R}^r\) with \(s_{p,q}(r)=\psi _{p,q}(h_t^{(r)})\) ;

Step 3: Check the similarity of two off-subgraph edges \((p_i,q_i)\) and \((p_j,q_j)\) withIf \(\textrm {SpectralSim}(i,j) \lt \varrho\) for a given threshold \(\varrho\) , edge \((p_i,q_i)\) is considered spectrally dissimilar to \((p_j,q_j)\) . Given a list of candidate off-subgraph edges, Algorithm 2 is proposed for edge similarity checking.

Algorithm 1 shows the algorithm flow for directed graph sparsification, where \(\mathbf {L_G}\) is the Laplacian matrix for original graph, \(\mathbf {L_S}\) is the Laplacian matrix of initial spanning tree, \(d_{\textrm {out}}\) is the user-defined outgoing degree for nodes, and \(\lambda _{\textrm {limit}}\) is the desired maximum generalized eigenvalue. Algorithm 2 obtains the edges after checking edges similarities for the off-subgraph edges, where \(\tilde{E}_{\textrm {list}}\) is the set of off-subgraph edges; \(E_{\textrm {list}}\) is the set of edges that will be added into sparsifier; \(d_p\) is the outgoing degree for node p. The complexity has been summarized as follows:

Recent research has focused on developing more efficient algorithms for solving undirected Laplacians [22, 24]. In this work, we will focus on solving asymmetric Laplacian matrices that correspond to directed graphs: solving the following linear system equations \(\mathbf {L_G x=b}\) , where the right-hand-side (RHS) vector \(\mathbf {b}\) lies in the left singular vector space, will be equivalent to solving the following problem:Let \(\mathbf {y= L_G^{\top +}x}\) , then we will first solve \(\mathbf {L_{G_u} y=b }\) . Once \(\mathbf {y}\) is obtained, we can get the solution \(\mathbf {x=L_G^\top y}\) . Since \(\mathbf {L_{G_u}}\) is a much denser matrix, \(\mathbf {L_{G_u}}\) should not be explicitly formed when solving \(\mathbf {L_{G_u} y=b }\) . To this end, iterative methods such as the preconditioned conjugate gradient (PCG) method can be leveraged for solving \(\mathbf {L_{G_u} y=b }\) with \(\mathbf {L_{S_u}}\) as the preconditioner. Note that only \(\mathbf {L_{S_u}}\) will be explicitly computed during PCG iterations.

The directed graph sparsifier can also be directly leveraged as a preconditioner for solving \(\mathbf {L_Gx = b}\) using existing iterative methods, such as the generalized minimal residual (GMRES) method [36]. GMRES is a widely-adopted Krylov-subspace iterative method for solving asymmetric matrices. Given an initial solution vector \(\mathbf {x_0}\) , GMRES gradually improves the solution \(\mathbf {x_m}\) of the mth iteration by minimizing the residue as follows:where \(\mathbf {r_0 = b-L_Gx_0}\) , and \(\mathcal {K}_m(\mathbf {L_G,r_0}) = span\lbrace \mathbf {r_0, L_Gr_0, L_G^2r_0,\ldots , L_G^{m-1}r_0}\rbrace\) denotes the Krylov subspace.

The idea of PageRank is to give a measurement of the importance for each web page. For example, PageRank algorithm aims at finding the most popular web pages, while the personalized PageRank algorithm aims at finding the pages that users will most likely visit. To state it mathematically, the PageRank vector \(\mathbf {\pi }\) satisfies the following equation:where \(0 \lt c \lt 1\) is the damping constant, and vector \(\mathbf {v_0}\) with non-negative coordinates, satisfying \(\mathbf {1^\top v_0 = 1}\) , is the personalization vector. The original, non-personalized definition of the PageRank is described when \(\mathbf {v_0}=\frac{1}{n}\mathbf {1}\) . Meanwhile, \(\mathbf {D_G^{-1}}\) can not be defined if there exist nodes that have no outgoing edges. To deal with such situation, a self-loop with a small edge weight can be added for each node.

It has been shown that partitioning and clustering of directed graphs can play very significant roles in a variety of applications related to machine learning [31], data mining and circuit synthesis and optimization [32], and so on. However, the efficiency of existing methods for partitioning directed graphs strongly depends on the complexity of the underlying graphs [31]. For an undirected graph, the eigenvectors corresponding to the first few smallest eigenvalues can be utilized for the spectral partitioning purpose [39]. For a directed graph G on the other hand, the eigenvectors corresponding to the first few different smallest eigenvalues of Laplacian \(\mathbf {L_{G_u}}\) will be required for directed graph partitioning. The eigenvalues according to the symmetrization of the directed graph in Figure 5 have a few multiplicities, which are shown in Figure 4. The partitioning result of the directed graph in Figure 5 will depend on the eigenvectors that correspond to eigenvalues of \(\mu _1,\mu _2,\mu _4,\mu _8\) . As shown in Figure 5, the spectral partitioning results can be quite different between the directed and undirected graph with the same set of nodes and edges.

The datasets are from SuiteSparse Matrix Collection [16]. If a node has only incoming edges or is isolated with the rest of nodes, this bode will be removed from the graph. The statistics of datasets are summarized in Table 2. The detailed description for each graph is shown as follows:

Figure 6 shows the spectral sensitivities of all the off-subgraph edges (e1 to e19 represented with blue color) in both directed and undirected graphs calculated using MATLAB’s eigs function and the proposed method based on (37) using the LAMG solver, respectively. Meanwhile, the spectral sensitivities of all the off-subgraph edges (e1 to e19) with respect to the dominant eigenvalues ( \(\lambda _{max}\) or \(\lambda _{1}\) ) in both directed and undirected graphs are plotted. We observe that spectral sensitivities for directed and undirected graphs are drastically different from each other. The reason is that the spectral sensitivities for off-subgraph edges in the directed graph depend on the edge directions. It is also observed that the approximate spectral sensitivities calculated by the proposed t-step power iterations with the LAMG solver match the true solution very well for both directed and undirected graphs.

Table 3 shows comprehensive results on directed graph spectral sparsification for a variety of real-world directed graphs using the proposed method, where \(|V_G|(|E_G|)\) denotes the number of nodes (edges) for the original directed graph G; \(|E_{S^0}|\) and \(|E_S|\) denote the numbers of edges in the initial subgraph \(S^0\) and final spectral sparsifier S. Notice that we will directly apply the MATLAB’s eigs function if graph size is relatively small ( \(|E_{S^0}|\lt 1E4\) ); otherwise, we will apply LAMG solver for better efficiency when calculating the generalized eigenvector \(\mathbf {h_t}\) . Note that a small diagonal entry with value of \(1e-6\) is added to all symmetrized undirected graphs during the calculation. We report the total runtime for the eigsolver using either the LAMG solver or eigs function. \(\frac{\lambda _{max, S^{0}}}{\lambda _{max} }\) denotes the reduction rate of the largest generalized eigenvalue of \(\mathbf {L^+_{S_u}} \mathbf {L_{G_u}}\) from initial sparsifier to final sparsifier. Ablation study. Since the proposed method is iteratively adding edges for forming the sparsifier. We demonstrate the performance of the runtime and generalized eigenvalue reduction with respect to the number of added edges in the sparsifier. Figure 7 shows the runtime scalability regarding to the number of off-subgraph edges ( \(|E_{added}|\) ) added in the final sparsifier for graph “gre \(\_1107\) ” (left), “big” (middle) and “gre \(\_115\) ” (right). It shows that the runtime scales linearly with the added number of edges for all three graphs. Figure 8 shows how \(\lambda _{max}(L_{G_u}, L_{S_u})\) is changing when including different number of edges in the sparsifier. We can observe that \(\lambda _{max}\) can be efficiently reduced when adding more edges in the sparsifier, especially at the early-stage of sparsifier construction. It also demonstrates that the most spectrally-critical edges can be efficiently identified and included at the early stage comparing to the edges that are less critical.

Since there are no other existing directed graph sparsification methods to be compared, we compare our proposed method with the existing undirected graph sparsification tool GRASS [18, 19, 43]. To this end, we first convert directed graphs into undirected ones ( \(G^{\prime }_u\) ) using \(\mathbf {A+A^\top }\) symmetrization. Then undirected graph sparsifiers \(S^{\prime }_u\) will be computed by GRASS. In the last, the directed graph sparsifiers can be constructed by recovering edge directions to the undirected sparsifier \(S^{\prime }_u\) . Note that a larger diagonal entry with value of \(1e-4\) is added to all symmetrized undirected graphs during the calculation. The experimental results have been shown in Table 4, where \(\lambda _{max}\) represents the largest generalized eigenvalue between the original graph and its final sparsifier. By keeping similar numbers of edges in the sparsifiers, we observe that the proposed spectral sparsification method consistently produces much better spectral sparsifiers than GRASS. Note that for graphs “harvard500” and “wordnet3”, we cannot include more edges into the sparsifiers \(S^{\prime }_u\) using GRASS, implying that the final \(\lambda _{max}(\mathbf {{L}_{{G_u}}},\mathbf {{L}_{{S_u}}})\) cannot be further reduced; on the other hand, our method is able to further reduce its condition number, achieving a much better spectral approximation level.

Figure 9 shows the relative residual ( \(res=\mathbf {\Vert L_Gx-b\Vert /\Vert b\Vert }\) ) and runtime plots when spectral sparsifiers are applied as the preconditioners for solving the Laplacians of directed graphs “hor”, “gre \(\_\) 115” and “gre \(\_\) 185”, respectively. As observed, the performance of the PCG solver has been substantially improved by leveraging sparsifier-based preconditioners. Note that for graph “hor” the plain PCG solver without using any preconditioner cannot converge to the desired accuracy within the maximum number of iterations (500 iterations). Figure 10 shows the relative residual and runtime plots when the preconditioners obtained via Incomplete LU (ILU) factorization of the original directed graphs and their spectral sparsifiers are applied for “wordnet3”, “harvard500” and “big”, respectively. “ILU( \(\cdot\) )” and “LU( \(\cdot\) )” indicate that ILU and LU decompositions have been leveraged to construct the preconditioners, respectively. “nnz” denotes the number of nonzeros in the preconditioners. The MATLAB’s built-in functions gmres, ilu, and lu with default settings have been applied in our experiments. Note that the GMRES iterations with preconditioners show much faster convergence for all test cases. It is also observed for each test case the preconditioner computed using the directed sparsifier always has lowest number of nonzeros (nnz).

Figure 11 shows the application of the proposed directed graph sparsification for computing (personalized) PageRank vectors with \(c=0.85\) , where the correlation of (personalized) PageRank results using the original graphs (x-axis) and sparsifiers (y-axis) are plotted for graphs “ibm32” (left), “mathworks100” (middle) and “gre \(\_\) 1107” (right), respectively. Note that a few steps of GS smoothing have been applied to remove the high-frequency errors to obtain the smoothed (personalized) PageRank vectors when using the sparsified graphs. We observe that the (personalized) PageRank vectors obtained from sparsifiers can well approximate the results computed with the original graphs.

Table 5 shows the detailed partitioning results on different graphs. Since there is no clear clue for spectral directed graph partitioning, we choose to perform spectral partitioning on the symmetrized undirected graph \(G_u\) and \(S_u\) , where two-way spectral partitioning are applied by utilizing the Fiedler Vector of its Laplacian matrix. np is the number of nodes that share the different partitions when comparing the partitioning results on graph \(G_u\) and \(S_u\) , where a smaller np indicates a more similar partitioning results between two graphs, thus a better spectral similarity between the original graph and the sparsifier. \(np/|V_G|\) can be considered as the percentage of the mismatched node over all node set. cut is the cut value between two partitions, which is equivalent to the number of edges connecting two partitions. \(\theta\) is the ratio cut [42] value that can be computed with the following equation given the partition \(V_i\) and \(V_j\) :Figures 12, 13, 14, and 15 show the partitioning results on the symmetrized graph \({G_u}\) and its symmetrized sparsifier \({S_u}\) for “ibm”, “peta”, “gre \(\_\) 1107”, and “big” graphs. As observed, very similar partitioning results have been obtained, indicating well preserved spectral properties within the spectrally-sparsified directed graph.