Network Reliability, in depth: two-terminal reliability, exact factoring, path and cut bounds, and Monte Carlo with importance sampling

The model: graphs whose edges fail

Take an undirected graph G with vertices V and edges E. Each edge e is up with probability p_e and down with probability q_e = 1 - p_e, independently of every other edge. Nodes are perfectly reliable for now. A state is the set of edges that are up, and there are 2^|E| of them.

Three standard measures differ only in which nodes must stay connected. Two-terminal reliability R_st is the probability that s and t are connected by up edges. K-terminal reliability requires a set K of nodes to be mutually connected, and all-terminal reliability requires the whole graph to stay connected. In practice you usually want the complement, unreliability U = 1 - R, because well-designed networks have R very close to 1 and the interesting number is how close.

Two structural facts do most of the work. A path of up edges from s to t is sufficient for success, and a cut of down edges separating s from t is sufficient for failure. Every method below either enumerates states, decomposes on edges, or exploits paths and cuts.

Why exact computation is hard

Summing the probabilities of all 2^|E| states is exact and hopeless beyond about 30 edges. Better exact algorithms exist, but no polynomial one is expected: two-terminal reliability is #P-complete (Valiant, 1979), and Provan and Ball (1983) showed that K-terminal and all-terminal reliability are #P-complete too, even on restricted graph classes. #P-complete problems count solutions, and an efficient exact algorithm would imply efficient algorithms for counting problems believed to be intractable.

Approximation behaves differently for the different measures. For all-terminal unreliability, Karger gave a fully polynomial randomised approximation scheme by showing that when the graph is likely to disconnect, it almost always does so along one of a polynomial number of near-minimum cuts. No such scheme is known for two-terminal reliability in general, which is why the practical toolkit is exact methods on small graphs, bounds, and carefully designed simulation.

Exact: factoring by contraction and deletion

The factoring theorem, also called pivotal decomposition, conditions on one edge e: either it is up, in which case its endpoints can be merged (contraction), or it is down, in which case it can be deleted.

R(G) = p_e R(G / e) + q_e R(G - e)

Recursion ends when s and t have been merged (R = 1) or no edges remain (R = 0). Edges whose endpoints are already merged become self-loops and are dropped. The implementation below represents contraction as a relabelling of vertices, and memoises on the (labels, remaining edges) pair so that identical subproblems are solved once.

from functools import lru_cache

def exact_factoring(n, edges, s, t):
    """Two-terminal reliability by R = p*R(G/e) + (1-p)*R(G-e).
    edges: list of (u, v, p_up) on vertices 0..n-1."""
    @lru_cache(maxsize=None)
    def rec(labels, es):
        if labels[s] == labels[t]:
            return 1.0
        es = tuple(e for e in es if labels[e[0]] != labels[e[1]])  # drop self-loops
        if not es:
            return 0.0
        (u, v, p), rest = es[0], es[1:]
        a, b = labels[u], labels[v]
        merged = tuple(a if x == b else x for x in labels)        # contract e
        return p * rec(merged, rest) + (1 - p) * rec(labels, rest)
    return rec(tuple(range(n)), tuple(edges))

Worst-case cost is still exponential, but production solvers make factoring practical for graphs of a few dozen edges by applying reductions before each split: a chain of series edges collapses to one edge with probability p₁p₂, a pair of parallel edges to one with probability 1 - q₁q₂, and degree-two vertices other than s and t disappear. Choosing the pivot edge well, for example one adjacent to s, also shrinks the tree.

Worked example: the bridge network

The bridge network has four nodes and five edges: s connects to a and b, a and b connect to t, and a cross edge joins a and b. It is the smallest graph that is not series-parallel, so it is the classic test of factoring. Pivot on the cross edge e3. If e3 is up, a and b merge, leaving two parallel pairs in series: (1 - q²)². If it is down, two parallel two-edge paths remain: 1 - (1 - p²)². Expanding gives the polynomial

R = 2p² + 2p³ - 5p⁴ + 2p⁵

With p = 0.9 on every edge, both the polynomial and the code above give 0.97848, so the unreliability is 0.02152. The two dominant failure modes are the two size-two cuts, {e1, e2} around s and {e4, e5} around t, each failing with probability 0.01.

Bounds from disjoint paths and cuts

Bounds cost polynomial time and catch mistakes in every other method. Both follow from independence. If you can find k edge-disjoint s-t paths, the network fails only if every one of them fails, and they fail independently because they share no edge:

R ≥ 1 - ∏_paths (1 - ∏_{e in path} p_e)

Likewise, if you can find edge-disjoint s-t cuts, the network works only if every cut keeps at least one edge up:

R ≤ ∏_cuts (1 - ∏_{e in cut} q_e)

For the bridge network with p = 0.9, the paths s-a-t and s-b-t give a lower bound of 1 - (1 - 0.81)² = 0.9639, and the cuts {e1, e2} and {e4, e5} give an upper bound of (1 - 0.01)² = 0.9801. The true value 0.97848 sits between them. Disjoint paths can be found with unit-capacity max flow and disjoint cuts by peeling layers of a BFS from s. The Esary-Proschan bounds use all minimal paths and cuts and are tighter, at the cost of enumerating them.

Crude Monte Carlo with union-find

Crude Monte Carlo samples a state by flipping a biased coin per edge, checks s-t connectivity, and averages. Union-find makes each check nearly linear in the number of edges, and samples are independent, so the estimator parallelises perfectly.

import math, random

def find(parent, x):
    while parent[x] != x:
        parent[x] = parent[parent[x]]   # path halving
        x = parent[x]
    return x

def connected(n, edges, up, s, t):
    parent = list(range(n))
    for (u, v, _), ok in zip(edges, up):
        if ok:
            a, b = find(parent, u), find(parent, v)
            if a != b:
                parent[a] = b
    return find(parent, s) == find(parent, t)

def crude_mc(n, edges, s, t, m, rng):
    """Estimate unreliability and its standard error from m samples."""
    fails = 0
    for _ in range(m):
        up = [rng.random() < p for _, _, p in edges]
        if not connected(n, edges, up, s, t):
            fails += 1
    u = fails / m
    return u, math.sqrt(u * (1 - u) / m)

On the bridge network with p = 0.9, 20,000 samples with seed 7 estimate U = 0.02175 with standard error 0.00103, within one standard error of the exact 0.02152.

The weakness is rare failure. The relative error of the estimate is roughly 1 / sqrt(m U), so estimating U = 10^-6 to 10% needs about 10⁸ samples. Raise every edge to p = 0.999 and the exact unreliability falls to 2.002 x 10^-6; 20,000 crude samples with seed 11 see zero failures and report an estimate of 0 with a standard error of 0, which looks like certainty and is simply ignorance.

Rare failures: importance sampling

Importance sampling draws states from a distribution in which failures are common, then corrects with the likelihood ratio. Sample each edge down with an inflated probability q', and weight each failing sample by the product over edges of q_e/q' for down edges and p_e/(1 - q') for up edges. The weighted average is an unbiased estimate of U.

def is_mc(n, edges, s, t, m, rng, q_tilt):
    """Unbiased unreliability estimate sampling failures at rate q_tilt."""
    acc = acc2 = 0.0
    for _ in range(m):
        w, up = 1.0, []
        for _, _, p in edges:
            if rng.random() < q_tilt:
                up.append(False); w *= (1 - p) / q_tilt
            else:
                up.append(True); w *= p / (1 - q_tilt)
        x = 0.0 if connected(n, edges, up, s, t) else w
        acc += x; acc2 += x * x
    mean = acc / m
    return mean, math.sqrt((acc2 / m - mean * mean) / m)

With q' = 0.3 and the same 20,000 samples and seed, the estimate is 1.99 x 10^-6 with standard error 5.5 x 10^-8, about 3% relative error. The tilt must be chosen with care: too small and failures stay rare, too large and the weights become wildly variable, which shows up as a standard error that jumps between runs. A good heuristic is to set q' so the expected number of down edges is close to the size of the minimum cut, since that is how the network most likely fails. Cut-based methods go further by sampling only states in which some near-minimum cut is down.

Modelling real networks

Real networks violate the tidy model in three ways, and each has a standard fix. Node failures: split each node v into v-in and v-out joined by an edge carrying the node's reliability, so a switch failure becomes an edge failure. Shared risk: two fibres in one conduit, or two links on one line card, fail together. Model each shared-risk group as a component whose failure removes all its edges, and sample components rather than edges; ignoring this is the most common way a reliability estimate comes out orders of magnitude optimistic. Repair: availability over time depends on repair rates as well as failure rates, so convert each component to a steady-state availability, MTBF divided by MTBF plus MTTR, before using it as p.

A leaf-spine fabric shows the stakes. Two leaves each connect to four spines, every link with availability 0.999. The leaves stay connected unless all four two-hop paths fail, so U = (1 - 0.999²)⁴, about 1.6 x 10^-11. The factoring code agrees to four significant figures only: computing 1 - R in floating point when R is within 10^-11 of 1 loses most of the digits, so compute unreliability directly when it is tiny. Put all four spines in one power domain with availability 0.9999 and the unreliability becomes about 10^-4, more than six orders of magnitude worse. The shared domain, not the link count, sets the number. The same reasoning applies to rails in a GPU fabric and to the cut structure studied in the Stoer-Wagner minimum cut.

Failure modes

• Independence assumed silently. Correlated failures dominate real outages; model shared-risk groups explicitly.
• Reporting R instead of U. 0.99999 and 0.9999 look alike; 10^-5 and 10^-4 do not. Report unreliability and compute it directly.
• Zero-failure Monte Carlo. No observed failures is not an estimate. Report the rule-of-three bound, 3/m at 95% confidence, or switch to importance sampling.
• Unstable weights. An importance sampler whose standard error varies by a factor of ten across seeds is not converged; reduce the tilt or use cut-based sampling.
• Memo blow-up. Exact factoring without reductions on a 60-edge graph will exhaust memory. Reduce first and cap the recursion.
• Disconnected after a cut vertex. If s and t sit on opposite sides of an articulation point, its reliability multiplies everything; find these with Tarjan bridges and articulation points first.

Trade-offs: choosing a method

In design work the usual loop is: bound first to see whether the design is in the right range, simulate to get a number with an error bar, and run exact factoring on reduced subnetworks to validate the simulator.

What to do next

1. Implement the factoring function and reproduce 0.97848 on the bridge network.
2. Check it against brute-force state enumeration on random graphs of up to 15 edges.
3. Add series and parallel reductions and measure how far they extend the exact range.
4. Implement the disjoint path and cut bounds and assert every estimate lies within them.
5. Run crude Monte Carlo and importance sampling at p = 0.999 and compare their errors.
6. Rebuild your real topology with node splitting and shared-risk groups before trusting any number.
7. Report unreliability, with a confidence interval, rather than reliability.