Computing Effective Resistances on Large Graphs Based on Approximate Inverse of Cholesky Factor ^††thanks: This work is supported by NSFC under grant No. 62090025, and Beijing Science and Technology Plan (No. Z221100007722025).

Suppose $G=(V,E,w)$ denotes a weighted undirected graph, where $V$ and $E$ are the sets of vertices (nodes) and edges, $w$ is a positive weight function. Let $n=|V|$ and $m=|E|$. The incidence matrix $B$ is defined to be an $m\times n$ matrix such that each row of $B$ corresponds to an edge in $E$ and each column of $B$ corresponds to a node in $V$. The entries of $B$ satisfy:

Let $W$ denote an $m\times m$ diagnoal matrix with $W(e,e)=w(e)$. Then the Laplacian matrix $L_{G}\in\mathbb{R}^{n\times n}$ is defined as:

Given two nodes $p$ and $q$, the effective resistance $R_{G}(p,q)$ across $p$ and $q$ refers to the voltage difference between $p$ and $q$ when a unit current flows into $G$ through $p$ and leaves through $q$. Formally, it can be defined as:

Here $L_{G}^{\dagger}$ denotes pseudo-inverse of $L_{G}$ and $e_{p,q}=e_{p}-e_{q}$, where $e_{p}$ is the $p$-th column of the identy matrix. $L_{G}$ is singular as the smallest eigenvalue is $0$. To handle the singularity of $L_{G}$, we introduce a ground node and connect it to a randomly selected node in each connected component of $G$. It corresponds to adding some small positive values to the diagonal elements of $L_{G}$. To simplify the notations, we still use $L_{G}$ to denote the resulted symmetric diagonally dominant (SDD) matrix.

Below we briefly review the power grid reduction method proposed in [8], which employs effective resistances based graph sparsification.

Power grid reduction [17] aims to reduce the original large power grid to a small one which contains all the port nodes. A port node is defined to be a node which is connected to a voltage or current source. The voltage drops of port nodes can be obtained by analyzing the reduced model, enabling more efficient analysis for large power grids. Schur complement based method [14] is adopted to eliminate the non-port nodes without loss of accuracy but tends to generate much denser models. To address this issue, a graph sparsification based power grid reduction approach was introduced [8], which consists of circuits partitioning, Schur complement based reduction, effective resistances based port merging and effective resistances based graph sparsification. We remark that all the port nodes have important physical information and should be preserved. Although at least half of the port nodes were eliminated in [8], the algorithm can be easily extended to the case where all the port nodes should be kept. We summarize the modified version as Alg. 1.

Among all these steps, computing effective resistances is the most time-consuming one. So, to obtain a fast power grid reduction algorithm, more efficient methods for computing effective resistances are demanded.

Computing effective resistances accurately is computational challenging. For each query $(p,q)$, computing $R_{G}(p,q)$ requires solving linear equations whose coefficient matrix is $L_{G}$. It should be noted that many applications require computing effective resistances for every edge $(p,q)\in E$. Solving $|E|$ linear equations takes at least $\Omega(|E|^{2})$ time, which can be prohibitive for large-scale problems.

To compute effective resistances more efficiently, a random projection based method was proposed in [4]. It is based on the fact that the effective resistance $R_{G}(p,q)$ can be written as the distance between two vectors:

For each node $p$, $W^{1/2}BL_{G}^{\dagger}e_{p}$ is an $m$-dimensional vector. With the Johnson-Lindenstraus Lemma, these vectors can be projected into a $k$-dimensional space, such that:

where $k=O(logm)$ and $Q\in R^{k\times m}$ is a random matrix whose elements are uniformly sampled from $\pm\frac{1}{\sqrt{k}}$. Note that using (5), only $k$ linear equations need to be solved to construct the matrix $QW^{1/2}BL_{G}^{\dagger}$ and then each effective resistance query can be answered in $O(k)=O(logm)$ time. A more practical version of this algorithm was presented in [1], which incorporates the random projection based method and fast linear equation solver [18]. However, due to the big hidden constants, it still takes a huge amount of time to compute effective resistances on large graphs with reasonably high accuracy.

Suppose $L_{G}$ is factorized with Cholesky factorization:

where $L$ is a lower triangular matrix. Then the effective resistance of query $(p,q)$ can be written as:

Eq. 7 shows that the effective resistance $R_{G}(p,q)$ can be written as the distance between the $p$-th column and the $q$-th column of $L^{-1}$. If $L^{-1}$ is available, the effective resistances can be computed efficiently, but computing and storing the dense $L^{-1}$ explicitly can be prohibitive for large-scale problems. So, our idea is to derive a sparse and approximate inverse of $L$.

Based on the observation that most elements in $L^{-1}$ are very small, we develop a method for computing a sparse approximation of $L^{-1}$. Let $Z=L^{-1}$ and $z_{j}$ be the $j$-th column of $Z$. Recall that $L$ is the Cholesky factor of $L_{G}$, it can be shown that all the diagonal elements in $L$ are positive and all the off-diagonal elements in $L$ are nonpositive [19].

The matrix $Z$ has some useful structural properties which are summarized as the following lemma. Here a matrix is nonnegative means that all the elements in it are nonnegative.

Eq. 8 suggests that we can first compute the $n$-th column of $L^{-1}$ and then compute the ($n$-1)-th, ($n$-2)-th, $\cdots$, and the 1st columns one by one. Suppose we have some sparse approximations to $z_{i}$, denoted by $\tilde{z}_{i}$. Then $z_{j}$ can be computed approximately by:

The $z_{j}^{*}$ can be computed efficiently because the $\tilde{z}_{i}$s are sparse. Note that many elements in $z_{j}^{*}$ are very small, which can be set to $0$. Here we adopt a prunning strategy to control the error in the sense of $1$-norm. Note that for a vector $x$, the $1$-norm of $x$, denoted by $\|x\|_{1}$, is defined as the sum of absolute value of each element. Let $trunc_{k}(x)$ denote the vector obtained by setting the $k$ smallest elements (in the sense of absolute values) in $x$ to $0$. Our strategy is to find the largest $k$ and to set $\tilde{z}_{j}=trunc_{k}(z_{j}^{*})$, such that:

where $\epsilon$ is a user-defined threshold. It can be implemented by first sorting the nonzero elements in $z_{j}^{*}$ and then finding the largest $k$ which satisfies the above constraint. We summarize the algorithm for computing sparse approximate inverse of Cholesky factor as Alg. 2.

Now we analyze the errors caused by approximating $L^{-1}$ with $\tilde{Z}$. Formally, we give a theorem below, which relates the approximation error $\|z_{p}-\tilde{z}_{p}\|_{1}$ to the depth of node $p$ in the filled graph. The filled graph of $G$ refers to the undirected graph corresponding to the matrix $L$ [19]. Let $G_{L}=(V,F)$ denote the filled graph, where $F=\{(i,j)|i\neq j~{}and~{}L_{i,j}\neq 0\}$. The depth of node $p$, denoted by $depth(p)$, is defined as:

Note that for most graphs that stem from real world, the maximum node depth in the filled graph is not very large, as reported in the next section. By setting a sufficiently small $\epsilon$, $\tilde{Z}$ approximates $L^{-1}$ very well.

Based on the approximate inverse technique, the effective resistance $R_{G}(p,q)$ can be computed as:

Now we analyze the errors of effective resistances. Let $z_{p,q}=z_{p}-z_{q}$, $\tilde{z}_{p,q}=\tilde{z}_{p}-\tilde{z}_{q}$, $\Delta z_{p}=\tilde{z}_{p}-z_{p}$, $\Delta z_{q}=\tilde{z}_{q}-z_{q}$ and $\Delta z_{p,q}=\tilde{z}_{p,q}-{z}_{p,q}$. Eq. 12 indicates that:

We can assume that $\|\Delta z_{p,q}\|_{1}$ is very small and $\Delta z_{p,q}$ is close to $\mathbf{0}$. Then by ignoring the second-order terms, we have:

Let $\tilde{R}_{p,q}=\|\tilde{z}_{p,q}\|_{2}^{2}$ denote the approximate effective resistance, then we obtain:

If we denote $\alpha_{p,q}=\frac{2\|{z}_{p,q}\|_{1}(\|z_{p}\|_{1}depth(p)+\|z_{q}\|_{1}depth(q))}{\|z_{p,q}\|_{2}^{2}}$, then

This implies that the relative error of effective resistance scales linearly with the parameter $\epsilon$. The smaller $\epsilon$ is, the closer $\tilde{R}_{p,q}$ is to $R_{p,q}$.

Computing effective resistances using (22) requires Cholesky factorization on $L_{G}$. However, for large-scale graphs, especially for graphs that stem from social networks, Cholesky factorization on $L_{G}$ can be very expensive. In this work, we propose to use incomplete Cholesky factorization instead. In incomplete Cholesky factorization, some fill-ins with very small absolute values are dropped, which corresponds to set some branches with large resistances to open and does not introduce large errors to effective resistances. We summarize the overall algorithm for computing effective resistances as follows.

The time complexity of the proposed algorithm is closely related to the number of nonzeros in $\tilde{Z}$. In our experiments, the number of nonzeros in $\tilde{Z}$ is about $Cnlogn$ where $C$ is a small constant (e.g. $C<20$). The reason behind this phenomenon may be the decay property of $L^{-1}$ [20]. Here we just assume the average number of nonzeros in one column of $\tilde{Z}$ is $O(logn)$. Incomplete Cholesky factorization takes $O(n)$ time and the number of nonzeros in $L$ is $O(n)$. Each iteration of Alg. 2 takes $O(\frac{nnz(L)}{n}logn+logn\cdot loglogn)=O(logn\cdot loglogn)$ time and the time complexity of Alg. 2 is $O(nlogn\cdot loglogn)$. So the overall time complexity of Alg. 3 is $O(nlogn\cdot loglogn)+O(|Q_{r}|logn)$.

In this subsection, we compare the proposed algorithm (Alg. 3) with the random projection based method [1], whose results are obtained by running the codes shared on [21]. We do not compare the algorithms in [2, 3] because these algorithms can only handle unweighted graphs. The results are listed in Table I. The test cases cover a great variety of graphs obtained from social networks, finite element analysis and circuit simulation problems [22, 23, 24]. For each case, effective resistances of all edges are computed, i.e. the set of effective resistance queries $Q_{r}=E$.

$T$ denotes the runtime for computing effective resistances. $E_{a}$ and $E_{m}$ denote the average and the maximum relative errors, which are estimated by randomly selecting 1000 edges, computing accurate effective resistances for these edges and then computing the errors. For Alg. 3, $T$ includes the time for incomplete Cholesky factorization, computing approximate inverse and computing effective resistances. The drop tolerance in incomplete Cholesky factorization is set to 1E-3 and the parameter $\epsilon$ in Alg. 2 is also set to 1E-3. We also record the maximum depth of the filled graph, which is defined in the last section and denoted by $dpt$, and the number of nonzeros in the random projection matrix $Q$ and in the appproximate inverse matrix $\tilde{Z}$, both of which are divided by $nlogn$. “-” means that it takes more than 10 hours.

From the results we see that, the proposed Alg. 3 achieves an average 168X speedup over the random projection based method [1]. Although the number of nonzeros in the random projection matrix $Q$ is much larger than that in the approximate inverse matrix $\tilde{Z}$, the proposed algorithm shows one to two orders of magnitude improvement in the avarage relative error and also significant reduction in the maximum relative error. For the largest case named “thupg10”, with 6.0E7 nodes and 1.0E8 edges, effective resistances for up to 1.0E8 node pairs can be computed within just 481 seconds, while the average relative error is about 0.17%, which demonstrates the extremely high scalability of the proposed algorithm.

Combining the proposed algorithm for computing effective resistances (Alg. 3) with the graph sparsification based power grid reduction framework (Alg. 1), we obtain a fast power grid reduction algorithm. Test cases are from the well-known IBM power grid benchmarks [23]. For graph partitioning, we use the widely adopted graph partitioner METIS [25] and the number of blocks are set to $\frac{\#ports}{50}$. For transient analysis, each case is simulated for $1000$ fixed-size time steps and both original models and reduced models are analyzed with the direct solver CHOLMOD [26] (performing just once matrix factorization). For DC incremental analysis, $10\%$ blocks of each benchmark are modified to mimick the design process where the initial power grid is modified to fix violations. We compare three power grid reduction methods, which only differ in computing effective resistances. The first method computes effective resistances accurately, while the other two compute effective resistances approximately, using the random projection based approach [1] and the proposed approach (Alg. 3) respectively.

The results are listed in Table II. $|V|$ and $|E|$ denote the number of nodes and resistors. $T_{\mathit{red}}$ denotes the time for power grid reduction. Note that in the application of DC incremental analysis, only the modified blocks need to be reduced, so the $T_{\mathit{red}}$ in incremental analysis is just about 10% of that in transient analysis. $T_{\mathit{tr}}$ and $T_{\mathit{inc}}$ are the time for transient analysis and DC incremental analysis, respectively. Err is the average absolute error and Rel is the relative error, which is computed by dividing the Err by the maximum voltage drop.

From the results we see that, with the proposed algorithm for computing effective resistances (Alg. 3), the time for power grid reduction can be reduced largely, showing an average 6.4X speedup over the power grid reduction approach based on accurate effective resistances. In terms of the reduction accuracy, there is almost no increase in reduction errors. This validates the efficiency and the effectiveness of the proposed power grid reduction approach. On the other hand, with the random projection based method, it still takes a large amount of time to generate reduced models. Meanwhile, the reduction accuracy is also deteriorated, which is due to the large errors of effective resistances.

As for the total time of transient analysis, which includes the time for power grid reduction and the time for analyzing the reduced power grid, the proposed method achieves 1.7X speedup averagely over the power grid reduction approach based on accurate effective resistances. Compared with analyzing the original power grid, the proposed method shows an average 2.9X speedup, with the relative error below 1.51% for all cases. The transient waveforms of node n0_20706300_8937900 and node n1_29561400_9521100 in case “ibmpg3t” are plotted in Fig. 1. They validate the accuracy of transient simulation using the proposed fast power grid reduction method.

In the application of DC incremental analysis, the proposed method shows significant improvement in the incremental power grid reduction stage, thus achieves an average 2.5X speedup for the total time over the power grid reduction method based on accurate effective resistances. Compared with analyzing the modified power grid directly using CHOLMOD (”Original” in Table II), the proposed method shows 4.2X speedups on average, while the relative error is below 1.34% for all cases.