Fast Sequential Importance Sampling to Estimate the Graph Reliability Polynomial

Abstract

The reliability polynomial of a graph counts its connected subgraphs of various sizes. Algorithms based on sequential importance sampling (SIS) have been proposed to estimate a graph’s reliability polynomial. We develop an improved SIS algorithm for estimating the reliability polynomial. The new algorithm runs in expected time O(mlogn α(m,n)) at worst and ≈m in practice, compared to Θ(m ²) for the previous algorithm. We analyze the error bounds of this algorithm, including comparison to alternative estimation algorithms. In addition to the theoretical analysis, we discuss methods for estimating the variance and describe experimental results on a variety of random graphs.

Appendices

Appendix A: Radius-Two Search for Bridges

In addition to detecting degree-one and degree-two vertices, Algorithm 2 conducts, in some circumstances, a breath-first search around a modified vertex v to find any easy bridges. The key idea of the search is that we assume that, if the vertex v is connected to a radius-two vertex, then it is connected to the entire graph. We search for edges which disconnect v from the radius-two vertices.

To begin, we let V ₁ be all vertices connected to v, and let V ₂=V−V ₁−v. A bridge can arise if there is a unique w∈V ₁ that is connected to V ₂. If furthermore this w is not connected to any vertex in V ₁ and if w is connected by a single edge to v, then the edge 〈v,w〉 is a bridge. If, on the other hand, w has a unique edge e connecting it to V ₂, then e is a bridge.

The following algorithm finds such bridges:

• 1. Initially set w=∅

• 2. For each x∈V ₁:

  • 3. For each distinct vertex y connected to x:

    • 4. If y∈V ₂ and w=∅, set w=x and continue with the next x.

    • 5. If y∈V ₂ and w≠∅, we find no bridges. Halt.

• 6. If w=∅, then the graph is disconnected. Return.

• 7. If there is a single edge connecting v to w, it is a bridge. Contract it and return.

• 8. If there is a single edge 〈w,u〉 with u∈V ₂, then it is a bridge; contract it and return.

• 9. No bridges have been found. Terminate.

This code requires a data structure which supports efficient enumeration of distinct vertices. A mere list of edges (which may be multi-edges) is not enough. With this data structure, this code runs in time $O(V_{1}^{2})$. Because Algorithm 2 only conducts this search if the degree of v is bounded by a constant $d_{\operatorname{max}}$, the total cost of these searches can be regarded as O(1). In practice, we did not implement this data structure; this abandons the guarantee on maximum running time but in practice is faster.

There is one further improvement we can make here. If in either step (8) or (9) there is a double edge e ₁,e ₂ such that simultaneously removing e ₁ and e ₂ disconnects the graph, then we can simplify the graph by contracting the two edges e ₁ and e ₂ and inserting a loop of weight w(e ₁)+w(e ₂).

Appendix B: Estimating Variance

In Sect. 6, we examined worst-case variance across all possible graphs of a given size. In this section, we consider variance from a different angle. Suppose we are given a graph G and use the algorithms we have described to estimate its reliability generator. How reliable are these estimates? In this section, we describe a method for estimating relative variance.

To begin, let us consider top-down algorithms such as Algorithm 0 or 1. In such algorithms, we may remove an edge e with probability P(e); if the edge is removed, then the corresponding estimate is multiplied by 1/P(e). P(e) is essentially the importance function attached to each edge. For Algorithm 0, P(e)=1/m for all edges; for Algorithm 1, P(e)=1/|D(G)| if e is a non-bridge, and P(e)=0 otherwise.

In such cases, the straight-forward way to estimate E[F(G)] and E[F(G)²] would be to take T samples F(G) and extract the sample mean and variance. Let $\hat{\mu}$ and $\hat{\sigma}^{2}$ be the sample estimates obtained in this way. We have that

so these estimates may be very inaccurate when σ is large. For the high-order coefficients of the reverse reliability generator, this makes it exponentially expensive to estimate μ,σ accurately. This is a very dangerous situation: when the estimates F(G) are very inaccurate, it is also very difficult to determine how inaccurate they are.

A better method of estimating accuracy is to estimate E[F(G)] and E[F(G)²] using importance sampling. Instead of choosing the edge e to remove with the indicated probability P(e), we use an alternative importance-sampling probability distribution

$$P'(e) = \frac{\mbox{\# of spanning trees of $G-e$}}{E-V+1}. $$

We emphasize that there are two distinct probability distributions in play here. P(e) is the probability distribution used by algorithm F itself; P′(e) is the probability distribution we are imposing in order to measure algorithm F.

This gives us the following SIS algorithm F ^(d) for estimating E[F(G)^d]:

Algorithm

F ^(d)

• 1. Set m ₀=1

• 2. For k=1,…,K

  • 3. Choose an edge e∈G with probability P′(e).

  • 4. Set $m_{k} \leftarrow m_{k-1} \frac{1}{P(e)^{d-1} P'(e)}$. Set G←G−e.

• 5. Return the estimate

  $$\hat{r}^{(k)} = \sum_{k=0}^{K} \frac{m_k x^k}{(k!)^d} $$

It is easy to verify that E[F ^(d)(G)]=E[F(G)^d].

The importance function P′(e) ensures that algorithm F ⁽¹⁾ estimates the high-order coefficient of E[F(G)] (the number of spanning trees of G) exactly (variance zero), although the estimates for the lower-order coefficients have higher variance. The estimate of the high-order coefficients of E[F(G)]² also has greatly improved variance, although the low-order coefficients are worse. We believe it is much safer to have an error estimate which is accurate when the error is large.

The estimate of σ ² provided by F ⁽²⁾ is typically the right order of magnitude, but for low-order coefficients it is not very precise. The sample variance is a better estimator for these low-order coefficients. To use the sample variance safely, we first estimate σ ² using F ⁽²⁾ for all coefficients. For those coefficients for which $\hat{\sigma}^{2}$ is sufficiently small, we replace our estimate of σ ² with the sample variance.

Note that to use this method, we must compute the importance function P′(e) for all edges e∈G at every stage of the recursion. To describe how this is done, we begin by recalling the Kirchoff matrix-tree theorem used to count the number of spanning trees of a graph G. Let A _G be the adjacency matrix of G, and let D be a diagonal matrix whose ith entry is the degree of vertex v _i . The Kirchoff formula states that κ(G)=detL, where L=(D−A _G )₁₁ denotes the minor of D−A _G obtained by removing the first row and column.

The SIS algorithm F ^(d) would appear to be very expensive if it computed κ naively from the Kirchoff formula. There are m iterations, in each of which we must κ(G−e) for each edge. Evaluating a determinant costs O(n ³) work and O(n ²) memory. In all, naive computations would cost O(m ² n ³) work and O(n ²) memory.

Reference [4] describes a variety of algorithms to compute κ and to update this structure as edges of G are removed. We will discuss one variant of those algorithms. The main idea is to keep track of λ(e)=κ(G−e)/κ(G) for each edge e∈G. These are precomputed for the original graph G. Next, when we update G by removing an edge e, we must update λ to the new λ′. We will find the following notation convenient. If e=〈i,j〉 is any edge in G, we define the column vector δ _e to be the vector which is +1 in coordinate i, is −1 in coordinate j, and is zero elsewhere. Observe that when edge e is removed from G, the matrix L changes by $\delta_{e} \delta_{e}^{T}$:

$$L_{G - e} = L_{G} + \delta_e \delta_e^T $$

Now let us examine how to update λ′:

As described in [4], this leads to the following technique for updating λ. Whenever we remove edge e from G, we use a sparse matrix inversion algorithm such as conjugate gradient to compute $u = L^{-1}_{G} \delta_{e}$. We then update λ(d) for each edge d∈G by

$$\lambda(d) \leftarrow\lambda(d) - (\delta_d^T u)^2 / \lambda(e) $$

In total, this technique allows to reduce the memory requirement to O(m) and the complexity to O(mn). Furthermore, in practice the conjugate gradient method needs far less than m  n iterations to compute an adequate approximation, so in practice the work requirement per iteration is more like ≈mlogn. This is still much more expensive than computing F itself, however it allows us to achieve accurate estimates for E[F(G)] and E[F(G)²] on graphs of interest. In many cases, straightforward sampling of F itself would have required infeasible computation to obtain accurate results.